File indexing completed on 2024-03-24 03:46:45
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 }