6#include "MeshIterator.h"
8#include "SpatialVector.h"
9#include "SpatialIndex.h"
10#include "RangeConvex.h"
12#include "HtmRangeIterator.h"
23 : m_level(level), m_buildLevel(buildLevel), m_numBuffers(numBuffers), htmDebug(0)
28 if (m_buildLevel > m_level)
29 m_buildLevel = m_level;
39 for (
int i = m_level; i--;)
48 magicNum = numTrixels;
49 degree2Rad = 3.1415926535897932385E0 / 180.0;
53 if (m_meshBuffer ==
nullptr)
55 fprintf(stderr,
"Out of memory allocating %d MeshBuffers.\n", numBuffers);
58 for (
int i = 0; i < numBuffers; i++)
67 for (BufNum i = 0; i < m_numBuffers; i++)
68 delete m_meshBuffer[i];
77bool HTMesh::performIntersection(
RangeConvex *convex, BufNum bufNum)
79 if (!validBufNum(bufNum))
82 convex->setOlevel(m_level);
85 HtmRangeIterator iterator(&range);
89 while (iterator.hasNext())
91 buffer->
append((Trixel)iterator.next() - magicNum);
96 fprintf(stderr,
"%s: trixel overflow.\n", name);
106 double d = cos(radius * degree2Rad);
111 if (!performIntersection(&convex, bufNum))
112 printf(
"In intersect(%f, %f, %f)\n", ra, dec, radius);
116void HTMesh::intersect(
double ra1,
double dec1,
double ra2,
double dec2,
double ra3,
double dec3, BufNum bufNum)
118 if (fabs(ra1 - ra3) + fabs(dec1 - dec3) < eps)
121 else if (fabs(ra1 - ra2) + fabs(dec1 - dec2) < eps)
124 else if (fabs(ra2 - ra3) + fabs(dec2 - dec3) < eps)
132 if (!performIntersection(&convex, bufNum))
133 printf(
"In intersect(%f, %f, %f, %f, %f, %f)\n", ra1, dec1, ra2, dec2, ra3, dec3);
137void HTMesh::intersect(
double ra1,
double dec1,
double ra2,
double dec2,
double ra3,
double dec3,
double ra4,
138 double dec4, BufNum bufNum)
140 if (fabs(ra1 - ra4) + fabs(dec1 - dec4) < eps)
141 return intersect(ra2, dec2, ra3, dec3, ra4, dec4);
143 else if (fabs(ra1 - ra2) + fabs(dec1 - dec2) < eps)
144 return intersect(ra2, dec2, ra3, dec3, ra4, dec4);
146 else if (fabs(ra2 - ra3) + fabs(dec2 - dec3) < eps)
147 return intersect(ra1, dec1, ra2, dec2, ra4, dec4);
149 else if (fabs(ra3 - ra4) + fabs(dec3 - dec4) < eps)
150 return intersect(ra1, dec1, ra2, dec2, ra4, dec4);
158 if (!performIntersection(&convex, bufNum))
159 printf(
"In intersect(%f, %f, %f, %f, %f, %f, %f, %f)\n", ra1, dec1, ra2, dec2, ra3, dec3, ra4, dec4);
162void HTMesh::toXYZ(
double ra,
double dec,
double *x,
double *y,
double *z)
167 double sinRa = sin(ra);
168 double cosRa = cos(ra);
169 double sinDec = sin(dec);
170 double cosDec = cos(dec);
183 double x1, y1, z1, x2, y2, z2;
191 toXYZ(ra1, dec1, &x1, &y1, &z1);
192 toXYZ(ra2, dec2, &x2, &y2, &z2);
197 len += fabs(y1 - y2);
198 len += fabs(z1 - z2);
202 printf(
"htmDebug = %d\n", htmDebug);
203 printf(
"p1 = (%f, %f, %f)\n", x1, y1, z1);
204 printf(
"p2 = (%f, %f, %f)\n", x2, y2, z2);
205 printf(
"edge: %f (radians) %f (degrees)\n", edge, edge / degree2Rad);
206 printf(
"len : %f (radians) %f (degrees)\n", len, len / degree2Rad);
213 double cx = y1 * z2 - z1 * y2;
214 double cy = z1 * x2 - x1 * z2;
215 double cz = x1 * y2 - y1 * x2;
218 printf(
"cp = (%f, %f, %f)\n", cx, cy, cz);
220 double norm = edge10 / (fabs(cx) + fabs(cy) + fabs(cz));
228 printf(
"cpn = (%f, %f, %f)\n", cx, cy, cz);
236 printf(
"cpf = (%f, %f, %f)\n", cx, cy, cz);
239 norm = sqrt(cx * cx + cy * cy + cz * cz);
240 double ra0 = atan2(cy, cx) / degree2Rad;
241 double dec0 = asin(cz / norm) / degree2Rad;
244 printf(
"new ra, dec = (%f, %f)\n", ra0, dec0);
251 if (!performIntersection(&convex, bufNum))
252 printf(
"In intersect(%f, %f, %f, %f)\n", ra1, dec1, ra2, dec2);
257 if (!validBufNum(bufNum))
259 return m_meshBuffer[bufNum];
264 if (!validBufNum(bufNum))
266 return m_meshBuffer[bufNum]->
size();
269void HTMesh::vertices(Trixel
id,
double *ra1,
double *dec1,
double *ra2,
double *dec2,
double *ra3,
double *dec3)
HTMesh(int level, int buildLevel, int numBuffers=1)
constructor.
Trixel index(double ra, double dec) const
returns the index of the trixel that contains the specified point.
int intersectSize(BufNum bufNum=0)
returns the number of trixels in the result buffer bufNum.
void intersect(double ra, double dec, double radius, BufNum bufNum=0)
NOTE: The intersect() routines below are all used to find the trixels needed to cover a geometric obj...
MeshBuffer * meshBuffer(BufNum bufNum=0)
returns a pointer to the MeshBuffer specified by bufNum.
The sole purpose of a MeshBuffer is to hold storage space for the results of an HTM inetersection and...
int error() const
returns the number of trixels that would have overflowed the buffer.
int size() const
the number of trixels in the result set
void reset()
prepare the buffer for a new result set
int append(Trixel trixel)
add trixels to the buffer
A spatial convex is composed of spatial constraints.
void intersect(const SpatialIndex *index, HtmRange *hr)
Intersect with index.
void add(SpatialConstraint &)
Add a constraint.
The Constraint is really a cone on the sky-sphere.
SpatialIndex is a quad tree of spherical triangles.
uint64 idByPoint(const SpatialVector &vector) const
find a node by giving a vector.
void nodeVertex(const uint64 id, SpatialVector &v1, SpatialVector &v2, SpatialVector &v3) const
return the actual vertex vectors
SpatialVector is a 3D vector usually living on the surface of the sphere.
float64 dec()
return dec - this norms the vector to 1 if not already done so
float64 ra()
return ra - this norms the vector to 1 if not already done so