Kstars

HTMesh.cpp
1#include <cstdlib>
2#include <iostream>
3
4#include "HTMesh.h"
5#include "MeshBuffer.h"
6#include "MeshIterator.h"
7
8#include "SpatialVector.h"
9#include "SpatialIndex.h"
10#include "RangeConvex.h"
11#include "HtmRange.h"
12#include "HtmRangeIterator.h"
13
14/******************************************************************************
15 * Note: There is "complete" checking for duplicate points in the line and
16 * polygon intersection routines below. This may have a slight performance
17 * impact on indexing lines and polygons.
18 *
19 * -- James B. Bowlin
20 *****************************************************************************/
21
22HTMesh::HTMesh(int level, int buildLevel, int numBuffers)
23 : m_level(level), m_buildLevel(buildLevel), m_numBuffers(numBuffers), htmDebug(0)
24{
25 name = "HTMesh";
26 if (m_buildLevel > 0)
27 {
28 if (m_buildLevel > m_level)
29 m_buildLevel = m_level;
30 htm = new SpatialIndex(m_level, m_buildLevel);
31 }
32 else
33 {
34 htm = new SpatialIndex(m_level);
35 }
36
37 edge = 2. / 3.14; // inverse of roughly 1/4 circle
38 numTrixels = 8;
39 for (int i = m_level; i--;)
40 {
41 numTrixels *= 4;
42 edge *= 2.0;
43 }
44 edge = 1.0 / edge; // roughly length of one edge (radians)
45 edge10 = edge / 10.0;
46 eps = 1.0e-6; // arbitrary small number
47
48 magicNum = numTrixels;
49 degree2Rad = 3.1415926535897932385E0 / 180.0;
50
51 // Allocate MeshBuffers
52 m_meshBuffer = (MeshBuffer **)malloc(sizeof(MeshBuffer *) * numBuffers);
53 if (m_meshBuffer == nullptr)
54 {
55 fprintf(stderr, "Out of memory allocating %d MeshBuffers.\n", numBuffers);
56 exit(0);
57 }
58 for (int i = 0; i < numBuffers; i++)
59 {
60 m_meshBuffer[i] = new MeshBuffer(this);
61 }
62}
63
64HTMesh::~HTMesh()
65{
66 delete htm;
67 for (BufNum i = 0; i < m_numBuffers; i++)
68 delete m_meshBuffer[i];
69 free(m_meshBuffer);
70}
71
72Trixel HTMesh::index(double ra, double dec) const
73{
74 return (Trixel)htm->idByPoint(SpatialVector(ra, dec)) - magicNum;
75}
76
77bool HTMesh::performIntersection(RangeConvex *convex, BufNum bufNum)
78{
79 if (!validBufNum(bufNum))
80 return false;
81
82 convex->setOlevel(m_level);
83 HtmRange range;
84 convex->intersect(htm, &range);
85 HtmRangeIterator iterator(&range);
86
87 MeshBuffer *buffer = m_meshBuffer[bufNum];
88 buffer->reset();
89 while (iterator.hasNext())
90 {
91 buffer->append((Trixel)iterator.next() - magicNum);
92 }
93
94 if (buffer->error())
95 {
96 fprintf(stderr, "%s: trixel overflow.\n", name);
97 return false;
98 };
99
100 return true;
101}
102
103// CIRCLE
104void HTMesh::intersect(double ra, double dec, double radius, BufNum bufNum)
105{
106 double d = cos(radius * degree2Rad);
107 SpatialConstraint c(SpatialVector(ra, dec), d);
108 RangeConvex convex;
109 convex.add(c); // [ed:RangeConvex::add]
110
111 if (!performIntersection(&convex, bufNum))
112 printf("In intersect(%f, %f, %f)\n", ra, dec, radius);
113}
114
115// TRIANGLE
116void HTMesh::intersect(double ra1, double dec1, double ra2, double dec2, double ra3, double dec3, BufNum bufNum)
117{
118 if (fabs(ra1 - ra3) + fabs(dec1 - dec3) < eps)
119 return intersect(ra1, dec1, ra2, dec2);
120
121 else if (fabs(ra1 - ra2) + fabs(dec1 - dec2) < eps)
122 return intersect(ra1, dec1, ra3, dec3);
123
124 else if (fabs(ra2 - ra3) + fabs(dec2 - dec3) < eps)
125 return intersect(ra1, dec1, ra2, dec2);
126
127 SpatialVector p1(ra1, dec1);
128 SpatialVector p2(ra2, dec2);
129 SpatialVector p3(ra3, dec3);
130 RangeConvex convex(&p1, &p2, &p3);
131
132 if (!performIntersection(&convex, bufNum))
133 printf("In intersect(%f, %f, %f, %f, %f, %f)\n", ra1, dec1, ra2, dec2, ra3, dec3);
134}
135
136// QUADRILATERAL
137void HTMesh::intersect(double ra1, double dec1, double ra2, double dec2, double ra3, double dec3, double ra4,
138 double dec4, BufNum bufNum)
139{
140 if (fabs(ra1 - ra4) + fabs(dec1 - dec4) < eps)
141 return intersect(ra2, dec2, ra3, dec3, ra4, dec4);
142
143 else if (fabs(ra1 - ra2) + fabs(dec1 - dec2) < eps)
144 return intersect(ra2, dec2, ra3, dec3, ra4, dec4);
145
146 else if (fabs(ra2 - ra3) + fabs(dec2 - dec3) < eps)
147 return intersect(ra1, dec1, ra2, dec2, ra4, dec4);
148
149 else if (fabs(ra3 - ra4) + fabs(dec3 - dec4) < eps)
150 return intersect(ra1, dec1, ra2, dec2, ra4, dec4);
151
152 SpatialVector p1(ra1, dec1);
153 SpatialVector p2(ra2, dec2);
154 SpatialVector p3(ra3, dec3);
155 SpatialVector p4(ra4, dec4);
156 RangeConvex convex(&p1, &p2, &p3, &p4);
157
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);
160}
161
162void HTMesh::toXYZ(double ra, double dec, double *x, double *y, double *z)
163{
164 ra *= degree2Rad;
165 dec *= degree2Rad;
166
167 double sinRa = sin(ra);
168 double cosRa = cos(ra);
169 double sinDec = sin(dec);
170 double cosDec = cos(dec);
171
172 *x = cosDec * cosRa;
173 *y = cosDec * sinRa;
174 *z = sinDec;
175}
176
177// Intersect a line segment by forming a very thin triangle to use for the
178// intersection. Use cross product to ensure we have a perpendicular vector.
179
180// LINE
181void HTMesh::intersect(double ra1, double dec1, double ra2, double dec2, BufNum bufNum)
182{
183 double x1, y1, z1, x2, y2, z2;
184
185 //if (ra1 == 0.0 || ra1 == 180.00) ra1 += 0.1;
186 //if (ra2 == 0.0 || ra2 == 180.00) ra2 -= 0.1;
187 //if (dec1 == 0.0 ) dec1 += 0.1;
188 //if (dec2 == 0.0 ) dec2 -= 0.1;
189
190 // convert to Cartesian. Ugh.
191 toXYZ(ra1, dec1, &x1, &y1, &z1);
192 toXYZ(ra2, dec2, &x2, &y2, &z2);
193
194 // Check if points are too close
195 double len;
196 len = fabs(x1 - x2);
197 len += fabs(y1 - y2);
198 len += fabs(z1 - z2);
199
200 if (htmDebug > 0)
201 {
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);
207 }
208
209 if (len < edge10)
210 return intersect(ra1, dec1, len, bufNum); // Use circular aperture
211
212 // Cartesian cross product => perpendicular!. Ugh.
213 double cx = y1 * z2 - z1 * y2;
214 double cy = z1 * x2 - x1 * z2;
215 double cz = x1 * y2 - y1 * x2;
216
217 if (htmDebug > 0)
218 printf("cp = (%f, %f, %f)\n", cx, cy, cz);
219
220 double norm = edge10 / std::sqrt(cx * cx + cy * cy + cz * cz);
221
222 // give it length edge/10
223 cx *= norm;
224 cy *= norm;
225 cz *= norm;
226
227 if (htmDebug > 0)
228 printf("cpn = (%f, %f, %f)\n", cx, cy, cz);
229
230 // add it to (ra1, dec1)
231 cx += x1;
232 cy += y1;
233 cz += z1;
234
235 if (htmDebug > 0)
236 printf("cpf = (%f, %f, %f)\n", cx, cy, cz);
237
238 // back to spherical
239 double ra0 = atan2(cy, cx) / degree2Rad;
240 double dec0 = atan2(cz, sqrt(cx * cx + cy * cy)) / degree2Rad;
241
242 if (htmDebug > 0)
243 printf("new ra, dec = (%f, %f)\n", ra0, dec0);
244
245 SpatialVector p1(ra1, dec1);
246 SpatialVector p0(ra0, dec0);
247 SpatialVector p2(ra2, dec2);
248 RangeConvex convex(&p1, &p0, &p2);
249
250 if (!performIntersection(&convex, bufNum))
251 printf("In intersect(%f, %f, %f, %f)\n", ra1, dec1, ra2, dec2);
252}
253
255{
256 if (!validBufNum(bufNum))
257 return nullptr;
258 return m_meshBuffer[bufNum];
259}
260
261int HTMesh::intersectSize(BufNum bufNum)
262{
263 if (!validBufNum(bufNum))
264 return 0;
265 return m_meshBuffer[bufNum]->size();
266}
267
268void HTMesh::vertices(Trixel id, double *ra1, double *dec1, double *ra2, double *dec2, double *ra3, double *dec3)
269{
270 SpatialVector v1, v2, v3;
271 htm->nodeVertex(id + magicNum, v1, v2, v3);
272 *ra1 = v1.ra();
273 *dec1 = v1.dec();
274 *ra2 = v2.ra();
275 *dec2 = v2.dec();
276 *ra3 = v3.ra();
277 *dec3 = v3.dec();
278}
HTMesh(int level, int buildLevel, int numBuffers=1)
constructor.
Definition HTMesh.cpp:22
Trixel index(double ra, double dec) const
returns the index of the trixel that contains the specified point.
Definition HTMesh.cpp:72
int intersectSize(BufNum bufNum=0)
returns the number of trixels in the result buffer bufNum.
Definition HTMesh.cpp:261
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...
Definition HTMesh.cpp:104
MeshBuffer * meshBuffer(BufNum bufNum=0)
returns a pointer to the MeshBuffer specified by bufNum.
Definition HTMesh.cpp:254
The sole purpose of a MeshBuffer is to hold storage space for the results of an HTM inetersection and...
Definition MeshBuffer.h:24
int error() const
returns the number of trixels that would have overflowed the buffer.
Definition MeshBuffer.h:49
int size() const
the number of trixels in the result set
Definition MeshBuffer.h:44
void reset()
prepare the buffer for a new result set
Definition MeshBuffer.h:32
int append(Trixel trixel)
add trixels to the buffer
A spatial convex is composed of spatial constraints.
Definition RangeConvex.h:60
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
This file is part of the KDE documentation.
Documentation copyright © 1996-2024 The KDE developers.
Generated on Mon Nov 4 2024 16:38:43 by doxygen 1.12.0 written by Dimitri van Heesch, © 1997-2006

KDE's Doxygen guidelines are available online.