File indexing completed on 2024-03-24 15:17:06

0001 #include <cstdlib>
0002 #include <iostream>
0003 
0004 #include "HTMesh.h"
0005 #include "MeshBuffer.h"
0006 #include "MeshIterator.h"
0007 
0008 #include "SpatialVector.h"
0009 #include "SpatialIndex.h"
0010 #include "RangeConvex.h"
0011 #include "HtmRange.h"
0012 #include "HtmRangeIterator.h"
0013 
0014 /******************************************************************************
0015  * Note: There is "complete" checking for duplicate points in the line and
0016  * polygon intersection routines below.  This may have a slight performance
0017  * impact on indexing lines and polygons.
0018  *
0019  * -- James B. Bowlin
0020  *****************************************************************************/
0021 
0022 HTMesh::HTMesh(int level, int buildLevel, int numBuffers)
0023     : m_level(level), m_buildLevel(buildLevel), m_numBuffers(numBuffers), htmDebug(0)
0024 {
0025     name = "HTMesh";
0026     if (m_buildLevel > 0)
0027     {
0028         if (m_buildLevel > m_level)
0029             m_buildLevel = m_level;
0030         htm = new SpatialIndex(m_level, m_buildLevel);
0031     }
0032     else
0033     {
0034         htm = new SpatialIndex(m_level);
0035     }
0036 
0037     edge       = 2. / 3.14; //  inverse of roughly 1/4 circle
0038     numTrixels = 8;
0039     for (int i = m_level; i--;)
0040     {
0041         numTrixels *= 4;
0042         edge *= 2.0;
0043     }
0044     edge   = 1.0 / edge; // roughly length of one edge (radians)
0045     edge10 = edge / 10.0;
0046     eps    = 1.0e-6; // arbitrary small number
0047 
0048     magicNum   = numTrixels;
0049     degree2Rad = 3.1415926535897932385E0 / 180.0;
0050 
0051     // Allocate MeshBuffers
0052     m_meshBuffer = (MeshBuffer **)malloc(sizeof(MeshBuffer *) * numBuffers);
0053     if (m_meshBuffer == nullptr)
0054     {
0055         fprintf(stderr, "Out of memory allocating %d MeshBuffers.\n", numBuffers);
0056         exit(0);
0057     }
0058     for (int i = 0; i < numBuffers; i++)
0059     {
0060         m_meshBuffer[i] = new MeshBuffer(this);
0061     }
0062 }
0063 
0064 HTMesh::~HTMesh()
0065 {
0066     delete htm;
0067     for (BufNum i = 0; i < m_numBuffers; i++)
0068         delete m_meshBuffer[i];
0069     free(m_meshBuffer);
0070 }
0071 
0072 Trixel HTMesh::index(double ra, double dec) const
0073 {
0074     return (Trixel)htm->idByPoint(SpatialVector(ra, dec)) - magicNum;
0075 }
0076 
0077 bool HTMesh::performIntersection(RangeConvex *convex, BufNum bufNum)
0078 {
0079     if (!validBufNum(bufNum))
0080         return false;
0081 
0082     convex->setOlevel(m_level);
0083     HtmRange range;
0084     convex->intersect(htm, &range);
0085     HtmRangeIterator iterator(&range);
0086 
0087     MeshBuffer *buffer = m_meshBuffer[bufNum];
0088     buffer->reset();
0089     while (iterator.hasNext())
0090     {
0091         buffer->append((Trixel)iterator.next() - magicNum);
0092     }
0093 
0094     if (buffer->error())
0095     {
0096         fprintf(stderr, "%s: trixel overflow.\n", name);
0097         return false;
0098     };
0099 
0100     return true;
0101 }
0102 
0103 // CIRCLE
0104 void HTMesh::intersect(double ra, double dec, double radius, BufNum bufNum)
0105 {
0106     double d = cos(radius * degree2Rad);
0107     SpatialConstraint c(SpatialVector(ra, dec), d);
0108     RangeConvex convex;
0109     convex.add(c); // [ed:RangeConvex::add]
0110 
0111     if (!performIntersection(&convex, bufNum))
0112         printf("In intersect(%f, %f, %f)\n", ra, dec, radius);
0113 }
0114 
0115 // TRIANGLE
0116 void HTMesh::intersect(double ra1, double dec1, double ra2, double dec2, double ra3, double dec3, BufNum bufNum)
0117 {
0118     if (fabs(ra1 - ra3) + fabs(dec1 - dec3) < eps)
0119         return intersect(ra1, dec1, ra2, dec2);
0120 
0121     else if (fabs(ra1 - ra2) + fabs(dec1 - dec2) < eps)
0122         return intersect(ra1, dec1, ra3, dec3);
0123 
0124     else if (fabs(ra2 - ra3) + fabs(dec2 - dec3) < eps)
0125         return intersect(ra1, dec1, ra2, dec2);
0126 
0127     SpatialVector p1(ra1, dec1);
0128     SpatialVector p2(ra2, dec2);
0129     SpatialVector p3(ra3, dec3);
0130     RangeConvex convex(&p1, &p2, &p3);
0131 
0132     if (!performIntersection(&convex, bufNum))
0133         printf("In intersect(%f, %f, %f, %f, %f, %f)\n", ra1, dec1, ra2, dec2, ra3, dec3);
0134 }
0135 
0136 // QUADRILATERAL
0137 void HTMesh::intersect(double ra1, double dec1, double ra2, double dec2, double ra3, double dec3, double ra4,
0138                        double dec4, BufNum bufNum)
0139 {
0140     if (fabs(ra1 - ra4) + fabs(dec1 - dec4) < eps)
0141         return intersect(ra2, dec2, ra3, dec3, ra4, dec4);
0142 
0143     else if (fabs(ra1 - ra2) + fabs(dec1 - dec2) < eps)
0144         return intersect(ra2, dec2, ra3, dec3, ra4, dec4);
0145 
0146     else if (fabs(ra2 - ra3) + fabs(dec2 - dec3) < eps)
0147         return intersect(ra1, dec1, ra2, dec2, ra4, dec4);
0148 
0149     else if (fabs(ra3 - ra4) + fabs(dec3 - dec4) < eps)
0150         return intersect(ra1, dec1, ra2, dec2, ra4, dec4);
0151 
0152     SpatialVector p1(ra1, dec1);
0153     SpatialVector p2(ra2, dec2);
0154     SpatialVector p3(ra3, dec3);
0155     SpatialVector p4(ra4, dec4);
0156     RangeConvex convex(&p1, &p2, &p3, &p4);
0157 
0158     if (!performIntersection(&convex, bufNum))
0159         printf("In intersect(%f, %f, %f, %f, %f, %f, %f, %f)\n", ra1, dec1, ra2, dec2, ra3, dec3, ra4, dec4);
0160 }
0161 
0162 void HTMesh::toXYZ(double ra, double dec, double *x, double *y, double *z)
0163 {
0164     ra *= degree2Rad;
0165     dec *= degree2Rad;
0166 
0167     double sinRa  = sin(ra);
0168     double cosRa  = cos(ra);
0169     double sinDec = sin(dec);
0170     double cosDec = cos(dec);
0171 
0172     *x = cosDec * cosRa;
0173     *y = cosDec * sinRa;
0174     *z = sinDec;
0175 }
0176 
0177 // Intersect a line segment by forming a very thin triangle to use for the
0178 // intersection.  Use cross product to ensure we have a perpendicular vector.
0179 
0180 // LINE
0181 void HTMesh::intersect(double ra1, double dec1, double ra2, double dec2, BufNum bufNum)
0182 {
0183     double x1, y1, z1, x2, y2, z2;
0184 
0185     //if (ra1 == 0.0 || ra1 == 180.00) ra1 += 0.1;
0186     //if (ra2 == 0.0 || ra2 == 180.00) ra2 -= 0.1;
0187     //if (dec1 == 0.0 ) dec1 += 0.1;
0188     //if (dec2 == 0.0 ) dec2 -= 0.1;
0189 
0190     // convert to Cartesian. Ugh.
0191     toXYZ(ra1, dec1, &x1, &y1, &z1);
0192     toXYZ(ra2, dec2, &x2, &y2, &z2);
0193 
0194     // Check if points are too close
0195     double len;
0196     len = fabs(x1 - x2);
0197     len += fabs(y1 - y2);
0198     len += fabs(z1 - z2);
0199 
0200     if (htmDebug > 0)
0201     {
0202         printf("htmDebug = %d\n", htmDebug);
0203         printf("p1 = (%f, %f, %f)\n", x1, y1, z1);
0204         printf("p2 = (%f, %f, %f)\n", x2, y2, z2);
0205         printf("edge: %f (radians) %f (degrees)\n", edge, edge / degree2Rad);
0206         printf("len : %f (radians) %f (degrees)\n", len, len / degree2Rad);
0207     }
0208 
0209     if (len < edge10)
0210         return intersect(ra1, len, bufNum);
0211 
0212     // Cartesian cross product => perpendicular!.  Ugh.
0213     double cx = y1 * z2 - z1 * y2;
0214     double cy = z1 * x2 - x1 * z2;
0215     double cz = x1 * y2 - y1 * x2;
0216 
0217     if (htmDebug > 0)
0218         printf("cp  = (%f, %f, %f)\n", cx, cy, cz);
0219 
0220     double norm = edge10 / (fabs(cx) + fabs(cy) + fabs(cz));
0221 
0222     // give it length edge/10
0223     cx *= norm;
0224     cy *= norm;
0225     cz *= norm;
0226 
0227     if (htmDebug > 0)
0228         printf("cpn  = (%f, %f, %f)\n", cx, cy, cz);
0229 
0230     // add it to (ra1, dec1)
0231     cx += x1;
0232     cy += y1;
0233     cz += z1;
0234 
0235     if (htmDebug > 0)
0236         printf("cpf  = (%f, %f, %f)\n", cx, cy, cz);
0237 
0238     // back to spherical
0239     norm        = sqrt(cx * cx + cy * cy + cz * cz);
0240     double ra0  = atan2(cy, cx) / degree2Rad;
0241     double dec0 = asin(cz / norm) / degree2Rad;
0242 
0243     if (htmDebug > 0)
0244         printf("new ra, dec = (%f, %f)\n", ra0, dec0);
0245 
0246     SpatialVector p1(ra1, dec1);
0247     SpatialVector p0(ra0, dec0);
0248     SpatialVector p2(ra2, dec2);
0249     RangeConvex convex(&p1, &p0, &p2);
0250 
0251     if (!performIntersection(&convex, bufNum))
0252         printf("In intersect(%f, %f, %f, %f)\n", ra1, dec1, ra2, dec2);
0253 }
0254 
0255 MeshBuffer *HTMesh::meshBuffer(BufNum bufNum)
0256 {
0257     if (!validBufNum(bufNum))
0258         return nullptr;
0259     return m_meshBuffer[bufNum];
0260 }
0261 
0262 int HTMesh::intersectSize(BufNum bufNum)
0263 {
0264     if (!validBufNum(bufNum))
0265         return 0;
0266     return m_meshBuffer[bufNum]->size();
0267 }
0268 
0269 void HTMesh::vertices(Trixel id, double *ra1, double *dec1, double *ra2, double *dec2, double *ra3, double *dec3)
0270 {
0271     SpatialVector v1, v2, v3;
0272     htm->nodeVertex(id + magicNum, v1, v2, v3);
0273     *ra1  = v1.ra();
0274     *dec1 = v1.dec();
0275     *ra2  = v2.ra();
0276     *dec2 = v2.dec();
0277     *ra3  = v3.ra();
0278     *dec3 = v3.dec();
0279 }