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 
22 HTMesh::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 
64 HTMesh::~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 
72 Trixel HTMesh::index(double ra, double dec) const
73 {
74  return (Trixel)htm->idByPoint(SpatialVector(ra, dec)) - magicNum;
75 }
76 
77 bool 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
104 void 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
116 void 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
137 void 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 
162 void 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
181 void 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, len, bufNum);
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 / (fabs(cx) + fabs(cy) + fabs(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  norm = sqrt(cx * cx + cy * cy + cz * cz);
240  double ra0 = atan2(cy, cx) / degree2Rad;
241  double dec0 = asin(cz / norm) / degree2Rad;
242 
243  if (htmDebug > 0)
244  printf("new ra, dec = (%f, %f)\n", ra0, dec0);
245 
246  SpatialVector p1(ra1, dec1);
247  SpatialVector p0(ra0, dec0);
248  SpatialVector p2(ra2, dec2);
249  RangeConvex convex(&p1, &p0, &p2);
250 
251  if (!performIntersection(&convex, bufNum))
252  printf("In intersect(%f, %f, %f, %f)\n", ra1, dec1, ra2, dec2);
253 }
254 
256 {
257  if (!validBufNum(bufNum))
258  return nullptr;
259  return m_meshBuffer[bufNum];
260 }
261 
262 int HTMesh::intersectSize(BufNum bufNum)
263 {
264  if (!validBufNum(bufNum))
265  return 0;
266  return m_meshBuffer[bufNum]->size();
267 }
268 
269 void HTMesh::vertices(Trixel id, double *ra1, double *dec1, double *ra2, double *dec2, double *ra3, double *dec3)
270 {
271  SpatialVector v1, v2, v3;
272  htm->nodeVertex(id + magicNum, v1, v2, v3);
273  *ra1 = v1.ra();
274  *dec1 = v1.dec();
275  *ra2 = v2.ra();
276  *dec2 = v2.dec();
277  *ra3 = v3.ra();
278  *dec3 = v3.dec();
279 }
int append(Trixel trixel)
add trixels to the buffer
Definition: MeshBuffer.cpp:26
int error() const
returns the number of trixels that would have overflowed the buffer.
Definition: MeshBuffer.h:49
void nodeVertex(const uint64 id, SpatialVector &v1, SpatialVector &v2, SpatialVector &v3) const
return the actual vertex vectors
uint64 idByPoint(const SpatialVector &vector) const
find a node by giving a vector.
The Constraint is really a cone on the sky-sphere.
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
float64 ra()
return ra - this norms the vector to 1 if not already done so
void intersect(const SpatialIndex *index, HtmRange *hr)
Intersect with index.
float64 dec()
return dec - this norms the vector to 1 if not already done so
int size() const
the number of trixels in the result set
Definition: MeshBuffer.h:44
MeshBuffer * meshBuffer(BufNum bufNum=0)
returns a pointer to the MeshBuffer specified by bufNum.
Definition: HTMesh.cpp:255
Trixel index(double ra, double dec) const
returns the index of the trixel that contains the specified point.
Definition: HTMesh.cpp:72
A spatial convex is composed of spatial constraints.
Definition: RangeConvex.h:59
HTMesh(int level, int buildLevel, int numBuffers=1)
constructor.
Definition: HTMesh.cpp:22
int intersectSize(BufNum bufNum=0)
returns the number of trixels in the result buffer bufNum.
Definition: HTMesh.cpp:262
void add(SpatialConstraint &)
Add a constraint.
void reset()
prepare the buffer for a new result set
Definition: MeshBuffer.h:32
This file is part of the KDE documentation.
Documentation copyright © 1996-2023 The KDE developers.
Generated on Sat Dec 9 2023 04:04:48 by doxygen 1.8.17 written by Dimitri van Heesch, © 1997-2006

KDE's Doxygen guidelines are available online.