Kstars

RangeConvex.cpp
1 //# Filename: RangeConvex.cpp # # The RangeConvex # classes are
2 //defined here. # # Author: Peter Z. Kunszt based on A. Szalay's code
3 //# # Date: October 23, 1998 # # SPDX-FileCopyrightText: 2000 Peter Z. Kunszt
4 //Alex S. Szalay, Aniruddha R. Thakar # The Johns Hopkins University #
5 //# Modification History: # # Oct 18, 2001 : Dennis C. Dinge --
6 //Replaced ValVec with std::vector
7 
8 #include "RangeConvex.h"
9 
10 #include "SpatialGeneral.h"
11 
12 #define N(n) index_->nodes_[(n)]
13 #define NC(n, m) index_->nodes_[(n)].childID_[(m)] // the children n->m
14 #define NV(m) index_->nodes_[id].v_[(m)] // the vertices of n
15 #define V(m) index_->vertices_[(m)] // the vertex vector m
16 
17 #define SGN(x) ((x) < 0 ? -1 : ((x) > 0 ? 1 : 0)) // signum
18 
19 // ===========================================================================
20 //
21 // Member functions for class RangeConvex
22 //
23 // ===========================================================================
24 
26 {
27 }
28 
29 /////////////CONSTRUCTOR FROM A TRIANGLE//////////////////
30 //
31 // Initialize convex from a triangle. The corners of these vectors
32 // form a triangle, so we just add three ZERO constraints to the convex
33 // where the direction is given by the cross product of the corners.
34 // Of course, the sign has to be determined (on which side of the triangle
35 // are we?) If the three points lie on one line, no constraints are added.
36 //
38 {
39  SpatialVector a1 = (*v2) ^ (*v3); // set directions of half-spheres
40  SpatialVector a2 = (*v3) ^ (*v1);
41  SpatialVector a3 = (*v1) ^ (*v2);
42  float64 s1 = a1 * (*v1); // we really need only the signs of these
43  float64 s2 = a2 * (*v2);
44  float64 s3 = a3 * (*v3);
45 
46  if (s1 * s2 * s3 != 0) // this is nonzero if not on one line
47  {
48  if (s1 < 0.0L)
49  a1 = (-1) * a1; // change sign if necessary
50  if (s2 < 0.0L)
51  a2 = (-1) * a2;
52  if (s3 < 0.0L)
53  a3 = (-1) * a3;
54  constraints_.push_back(SpatialConstraint(a1, 0.0)); // we don't care about the
55  constraints_.push_back(SpatialConstraint(a2, 0.0)); // order since all angles are
56  constraints_.push_back(SpatialConstraint(a3, 0.0)); // 90 degrees.
57  }
58  sign_ = zERO;
59 }
60 
61 /////////////CONSTRUCTOR FROM A RECTANGLE/////////////////
62 //
63 // Initialize convex from a rectangle. The vectors that form a rectangle
64 // may be in any order, the code finds the edges by itself.
65 // If one of the vectors lies within the triangle formed by the
66 // other three vectors, the previous constructor is used.
67 //
69  const SpatialVector *v4)
70 {
71  int i, j, k, l, m; // indices
72  // to simplify things, copy input into a 4-array
73  const SpatialVector *v[4] = { v1, v2, v3, v4 };
74  SpatialVector d[6];
75  float64 s[6][2];
76  for (i = 0, k = 0; i < 4; i++)
77  for (j = i + 1; j < 4; j++, k++) // set directions of half-spheres
78  {
79  d[k] = (*v[i]) ^ (*v[j]); // two of these are diagonals.
80  d[k].normalize();
81  for (l = 0, m = 0; l < 4; l++)
82  if (l != i && l != j)
83  s[k][m++] = d[k] * (*v[l]); // set the 'sign'
84  }
85 
86  // the sides are characterized by having both other corners
87  // to the same (inner) side. so it is easy to find the edges.
88  // again, the sign has to be taken care of -> direction of d
89  // the nice thing here is that if one of the corners is inside
90  // a triangles formed by the other three, only 3 constraints get
91  // added.
92  for (i = 0; i < 6; i++)
93  if (s[i][0] * s[i][1] > 0.0) // not >= because we don't want aligned corners
94  constraints_.push_back(SpatialConstraint((s[i][0] > 0.0 ? d[i] : (-1 * d[i])), 0.0));
95 
96  // Special cases: 1
97  // if three of the corners are aligned, we end up with
98  // only two constraints. Find the third and append it.
99  // Indeed, there are 3 identical constraints among the d[],
100  // so the first that qualifies gets appended.
101  if (constraints_.size() == 2)
102  {
103  for (i = 0; i < 6; i++)
104  if (s[i][0] == 0.0 || s[i][1] == 0.0)
105  {
106  constraints_.push_back(SpatialConstraint(((s[i][0] + s[i][1]) > 0.0 ? d[i] : (-1 * d[i])), 0.0));
107  break;
108  }
109  }
110  // Special cases: 2
111  // if all four corners are aligned, no constraints have been appended.
112  sign_ = zERO;
113 }
114 
115 /////////////ADD//////////////////////////////////////////
116 //
118 {
119  constraints_.push_back(c);
120  // order constraints: by ascending opening angle. Since we append
121  // always at the end, we only need one ordering sweep starting at
122  // the end
123  for (size_t i = constraints_.size() - 1; i > 0; i--)
124  {
125  if (constraints_[i].s_ < constraints_[i - 1].s_)
126  {
127  SpatialConstraint tmp(constraints_[i]);
128  constraints_[i] = constraints_[i - 1];
129  constraints_[i - 1] = tmp;
130  }
131  }
132 
133  if (constraints_.size() == 1) // first constraint
134  {
135  sign_ = c.sign_;
136  return;
137  }
138 
139  switch (sign_)
140  {
141  case nEG:
142  if (c.sign_ == pOS)
143  sign_ = mIXED;
144  break;
145  case pOS:
146  if (c.sign_ == nEG)
147  sign_ = mIXED;
148  break;
149  case zERO:
150  sign_ = c.sign_;
151  break;
152  case mIXED:
153  break;
154  }
155 }
156 
157 /////////////SIMPLIFY0////////////////////////////////////
158 // simplify0: simplify zERO convexes. calculate corners of convex
159 // and the bounding circle.
160 //
161 // zERO convexes are made up of constraints which are all great
162 // circles. It can happen that some of the constraints are redundant.
163 // For example, if 3 of the great circles define a triangle as the convex
164 // which lies fully inside the half sphere of the fourth constraint,
165 // that fourth constraint is redundant and will be removed.
166 //
167 // The algorithm is the following:
168 //
169 // zero-constraints are half-spheres, defined by a single normalized
170 // vector v, pointing in the direction of that half-sphere.
171 //
172 // Two zero-constraints intersect at
173 //
174 // i = +- v x v
175 // 1,2 1 2
176 //
177 // the vector cross product of their two defining vectors.
178 //
179 // The two vectors i1,2 are tested against every other constraint in
180 // the convex if they lie within their half-spheres. Those
181 // intersections i which lie within every other constraint, are stored
182 // into corners_.
183 //
184 // Constraints that do not have a single corner on them, are dropped.
185 //
186 
187 void RangeConvex::simplify0()
188 {
189  size_t i, j, k;
190  SpatialVector vi1, vi2;
191  std::vector<size_t> cornerConstr1, cornerConstr2, removeConstr;
192  std::vector<SpatialVector> corner;
193  if (constraints_.size() == 1) // for one constraint, it is itself the BC
194  {
195  boundingCircle_ = constraints_[0];
196  return;
197  // For 2 constraints, take the bounding circle a 0-constraint...
198  // this is by no means optimal, but the code is optimized for at least
199  // 3 zERO constraints... so this is acceptable.
200  }
201  else if (constraints_.size() == 2)
202  {
203  // test for constraints being identical - rule 1 out
204  if (constraints_[0].a_ == constraints_[1].a_)
205  {
206  constraints_.erase(constraints_.end() - 1);
207  boundingCircle_ = constraints_[0];
208  return;
209  }
210  // test for constraints being two disjoint half spheres - empty convex!
211  if (constraints_[0].a_ == (-1.0) * constraints_[1].a_)
212  {
213  constraints_.clear();
214  return;
215  }
216  boundingCircle_ = SpatialConstraint(constraints_[0].v() + constraints_[1].v(), 0);
217  return;
218  }
219 
220  // Go over all pairs of constraints
221  for (i = 0; i < constraints_.size() - 1; i++)
222  {
223  bool ruledout = true;
224  for (j = i + 1; j < constraints_.size(); j++)
225  {
226  // test for constraints being identical - rule i out
227  if (constraints_[i].a_ == constraints_[j].a_)
228  break;
229  // test for constraints being two disjoint half spheres - empty convex!
230  if (constraints_[i].a_ == (-1.0) * constraints_[j].a_)
231  {
232  constraints_.clear();
233  return;
234  }
235  // vi1 and vi2 are their intersection points
236  vi1 = constraints_[i].a_ ^ constraints_[j].a_;
237  vi1.normalize();
238  vi2 = (-1.0) * vi1;
239  bool vi1ok = true, vi2ok = true;
240  // now test whether vi1 or vi2 or both are inside every other constraint.
241  // if yes, store them in the corner array.
242  for (k = 0; k < constraints_.size(); k++)
243  {
244  if (k == i || k == j)
245  continue;
246  if (vi1ok && vi1 * constraints_[k].a_ <= 0.0)
247  vi1ok = false;
248  if (vi2ok && vi2 * constraints_[k].a_ <= 0.0)
249  vi2ok = false;
250  if (!vi1ok && !vi2ok)
251  break;
252  }
253  if (vi1ok)
254  {
255  corner.push_back(vi1);
256  cornerConstr1.push_back(i);
257  cornerConstr2.push_back(j);
258  ruledout = false;
259  }
260  if (vi2ok)
261  {
262  corner.push_back(vi2);
263  cornerConstr1.push_back(i);
264  cornerConstr2.push_back(j);
265  ruledout = false;
266  }
267  }
268  // is this constraint ruled out? i.e. none of its intersections
269  // with other constraints are corners... remove it from constraints_ list.
270  if (ruledout)
271  removeConstr.push_back(i);
272  }
273 
274  // Now set the corners into their correct order, which is an
275  // anti-clockwise walk around the polygon.
276  //
277  // start at any corner. so take the first.
278 
279  corners_.clear();
280  corners_.push_back(corner[0]);
281 
282  // The trick is now to start off into the correct direction.
283  // this corner has two edges it can walk. we have to take the
284  // one where the convex lies on its left side.
285  i = cornerConstr1[0]; // the i'th constraint and j'th constraint
286  j = cornerConstr2[0]; // intersect at 0'th corner
287  size_t c1(0), c2(0), k1(0), k2(0);
288  // Now find the other corner where the i'th and j'th constraints intersect.
289  // Store the corner in vi1 and vi2, and the other constraint indices
290  // in c1,c2.
291  for (k = 1; k < cornerConstr1.size(); k++)
292  {
293  if (cornerConstr1[k] == i)
294  {
295  vi1 = corner[k];
296  c1 = cornerConstr2[k];
297  k1 = k;
298  }
299  if (cornerConstr2[k] == i)
300  {
301  vi1 = corner[k];
302  c1 = cornerConstr1[k];
303  k1 = k;
304  }
305  if (cornerConstr1[k] == j)
306  {
307  vi2 = corner[k];
308  c2 = cornerConstr2[k];
309  k2 = k;
310  }
311  if (cornerConstr2[k] == j)
312  {
313  vi2 = corner[k];
314  c2 = cornerConstr1[k];
315  k2 = k;
316  }
317  }
318  // Now test i'th constraint-edge ( corner 0 and corner k ) whether
319  // it is on the correct side (left)
320  //
321  // ( (corner(k) - corner(0)) x constraint(i) ) * corner(0)
322  //
323  // is >0 if yes, <0 if no...
324  //
325  size_t c, currentCorner;
326  if (((vi1 - corner[0]) ^ constraints_[i].a_) * corner[0] > 0)
327  {
328  corners_.push_back(vi1);
329  c = c1;
330  currentCorner = k1;
331  }
332  else
333  {
334  corners_.push_back(vi2);
335  c = c2;
336  currentCorner = k2;
337  }
338  // now append the corners that match the index c until we got corner 0 again
339  // currentCorner holds the current corners index
340  // c holds the index of the constraint that has just been intersected with
341  // So:
342  // x We are on a constraint now (i or j from before), the second corner
343  // is the one intersecting with constraint c.
344  // x Find other corner for constraint c.
345  // x Save that corner, and set c to the constraint that intersects with c
346  // at that corner. Set currentcorner to that corners index.
347  // x Loop until 0th corner reached.
348  while (currentCorner)
349  {
350  for (k = 0; k < cornerConstr1.size(); k++)
351  {
352  if (k == currentCorner)
353  continue;
354  if (cornerConstr1[k] == c)
355  {
356  if ((currentCorner = k) == 0)
357  break;
358  corners_.push_back(corner[k]);
359  c = cornerConstr2[k];
360  break;
361  }
362  if (cornerConstr2[k] == c)
363  {
364  if ((currentCorner = k) == 0)
365  break;
366  corners_.push_back(corner[k]);
367  c = cornerConstr1[k];
368  break;
369  }
370  }
371  }
372  // Remove all redundant constraints
373  for (i = 0; i < removeConstr.size(); i++)
374  constraints_.erase(constraints_.end() - removeConstr[i] - 1);
375 
376  // Now calculate the bounding circle for the convex.
377  // We take it as the bounding circle of the triangle with
378  // the widest opening angle. All triangles made out of 3 corners
379  // are considered.
380  boundingCircle_.d_ = 1.0;
381  if (constraints_.size() >= 3)
382  {
383  for (i = 0; i < corners_.size(); i++)
384  for (j = i + 1; j < corners_.size(); j++)
385  for (k = j + 1; k < corners_.size(); k++)
386  {
387  SpatialVector v = (corners_[j] - corners_[i]) ^ (corners_[k] - corners_[j]);
388  v.normalize();
389  // Set the correct opening angle: Since the plane cutting
390  // out the triangle also correctly cuts out the bounding cap
391  // of the triangle on the sphere, we can take any corner to
392  // calculate the opening angle
393  float64 d = v * corners_[i];
394  if (boundingCircle_.d_ > d)
395  boundingCircle_ = SpatialConstraint(v, d);
396  }
397  }
398 }
399 
400 /////////////SIMPLIFY/////////////////////////////////////
401 // simplify: We have the following decision tree for the
402 // simplification of convexes:
403 //
404 // Always test two constraints against each other. We have
405 //
406 // * If both constraints are pOS
407 // # If they intersect: keep both
408 // # If one lies in the other: drop the larger one
409 // # Else: disjunct. Empty convex, stop.
410 //
411 // * If both constraints are nEG
412 // # If they intersect or are disjunct: ok
413 // # Else: one lies in the other, drop smaller 'hole'
414 //
415 // * Mixed: one pOS, one nEG
416 // # No intersection, disjunct: pOS is redundant
417 // # Intersection: keep both
418 // # pOS within nEG: empty convex, stop.
419 // # nEG within pOS: keep both.
420 //
421 
423 {
424  if (sign_ == zERO)
425  {
426  simplify0(); // treat zERO convexes separately
427  return;
428  }
429 
430  size_t i, j;
431  size_t clen;
432  bool redundancy = true;
433 
434  while (redundancy)
435  {
436  redundancy = false;
437  clen = constraints_.size();
438 
439  for (i = 0; i < clen; i++)
440  {
441  for (j = 0; j < i; j++)
442  {
443  int test;
444 
445  // don't bother with two zero constraints
446  if (constraints_[i].sign_ == zERO && constraints_[j].sign_ == zERO)
447  continue;
448 
449  // both pos or zero
450  if ((constraints_[i].sign_ == pOS || constraints_[i].sign_ == zERO) &&
451  (constraints_[j].sign_ == pOS || constraints_[j].sign_ == zERO))
452  {
453  if ((test = testConstraints(i, j)) == 0)
454  continue; // intersection
455  if (test < 0) // disjoint ! convex is empty
456  {
457  constraints_.clear();
458  return;
459  }
460  // one is redundant
461  if (test == 1)
462  constraints_.erase(constraints_.end() - i - 1);
463  else if (test == 2)
464  constraints_.erase(constraints_.end() - j - 1);
465  else
466  continue; // intersection
467  redundancy = true; // we did cut out a constraint -> do the loop again
468  break;
469  }
470 
471  // both neg or zero
472  if ((constraints_[i].sign_ == nEG) && (constraints_[j].sign_ == nEG))
473  {
474  if ((test = testConstraints(i, j)) <= 0)
475  continue; // ok
476  // one is redundant
477  if (test == 1)
478  constraints_.erase(constraints_.end() - 1 - j);
479  else if (test == 2)
480  constraints_.erase(constraints_.end() - 1 - i);
481  else
482  continue; // intersection
483  redundancy = true; // we did cut out a constraint -> do the loop again
484  break;
485  }
486 
487  // one neg, one pos/zero
488  if ((test = testConstraints(i, j)) == 0)
489  continue; // ok: intersect
490  if (test < 0) // neg is redundant
491  {
492  if (constraints_[i].sign_ == nEG)
493  constraints_.erase(constraints_.end() - 1 - i);
494  else
495  constraints_.erase(constraints_.end() - 1 - j);
496  redundancy = true; // we did cut out a constraint -> do the loop again
497  break;
498  }
499  // if the negative constraint is inside the positive: continue
500  if ((constraints_[i].sign_ == nEG && test == 2) || (constraints_[j].sign_ == nEG && test == 1))
501  continue;
502  // positive constraint in negative: convex is empty!
503  constraints_.clear();
504  return;
505  }
506  if (redundancy)
507  break;
508  }
509  }
510 
511  // reset the sign of the convex
512  sign_ = constraints_[0].sign_;
513  for (i = 1; i < constraints_.size(); i++)
514  {
515  switch (sign_)
516  {
517  case nEG:
518  if (constraints_[i].sign_ == pOS)
519  sign_ = mIXED;
520  break;
521  case pOS:
522  if (constraints_[i].sign_ == nEG)
523  sign_ = mIXED;
524  break;
525  case zERO:
526  sign_ = constraints_[i].sign_;
527  break;
528  case mIXED:
529  break;
530  }
531  }
532 
533  if (constraints_.size() == 1) // for one constraint, it is itself the BC
534  boundingCircle_ = constraints_[0];
535  else if (sign_ == pOS)
536  boundingCircle_ = constraints_[0];
537 }
538 
539 /////////////TESTCONSTRAINTS//////////////////////////////
540 // testConstraints: Test for the relative position of two constraints.
541 // Returns 0 if they intersect
542 // Returns -1 if they are disjoint
543 // Returns 1 if j is in i
544 // Returns 2 if i is in j
545 //
546 int RangeConvex::testConstraints(size_t i, size_t j)
547 {
548  float64 phi = ((constraints_[i].sign_ == nEG ? (-1 * constraints_[i].a_) : constraints_[i].a_) *
549  (constraints_[j].sign_ == nEG ? (-1 * constraints_[j].a_) : constraints_[j].a_));
550  phi = (phi <= -1.0L + gEpsilon ? gPi : acos(phi)); // correct for math lib -1.0
551  float64 a1 = (constraints_[i].sign_ == pOS ? constraints_[i].s_ : gPi - constraints_[i].s_);
552  float64 a2 = (constraints_[j].sign_ == pOS ? constraints_[j].s_ : gPi - constraints_[j].s_);
553 
554  if (phi > a1 + a2)
555  return -1;
556  if (a1 > phi + a2)
557  return 1;
558  if (a2 > phi + a1)
559  return 2;
560  return 0;
561 }
562 
563 /////////////INTERSECT////////////////////////////////////
564 // used by intersect.cpp application
565 //
566 void RangeConvex::intersect(const SpatialIndex *idx, HtmRange *htmrange)
567 {
568  hr = htmrange;
569  index_ = idx;
570  addlevel_ = idx->maxlevel_ - idx->buildlevel_;
571 
572  simplify(); // don't work too hard...
573 
574  if (constraints_.empty())
575  return; // nothing to intersect!!
576 
577  // Start with root nodes (index = 1-8) and intersect triangles
578  for (uint32 i = 1; i <= 8; i++)
579  {
580  testTrixel(i);
581  }
582 }
583 
584 ////////////SAVETRIXEL
585 //
586 inline void RangeConvex::saveTrixel(uint64 htmid)
587 {
588  // Some application want a complete htmid20 range
589  //
590  int level, i, shifts;
591  uint64 lo, hi;
592 #ifdef _WIN32
593  uint64 IDHIGHBIT = 1;
594  uint64 IDHIGHBIT2 = 1;
595  IDHIGHBIT = IDHIGHBIT << 63;
596  IDHIGHBIT2 = IDHIGHBIT2 << 62;
597 #endif
598 
599  for (i = 0; i < IDSIZE; i += 2)
600  {
601  if ((htmid << i) & IDHIGHBIT)
602  break;
603  }
604 
605  level = (IDSIZE - i) >> 1;
606  level -= 2;
607  if (level < olevel)
608  {
609  /* Size is the length of the string representing the name of the
610  trixel, the level is size - 2
611  */
612  shifts = (olevel - level) << 1;
613  lo = htmid << shifts;
614  hi = lo + ((uint64)1 << shifts) - 1;
615  }
616  else
617  {
618  lo = hi = htmid;
619  }
620  hr->mergeRange(lo, hi);
621 
622  return;
623 }
624 
625 /////////////TRIANGLETEST/////////////////////////////////
626 // testTrixel: this is the main test of a triangle vs a Convex. It
627 // will properly mark up the flags for the triangular node[index], and
628 // all its children
629 SpatialMarkup RangeConvex::testTrixel(uint64 id)
630 {
631  SpatialMarkup mark;
632  uint64 childID;
633  uint64 tid;
634 
635  // const struct SpatialIndex::QuadNode &indexNode = index_->nodes_[id];
636  const struct SpatialIndex::QuadNode *indexNode = &index_->nodes_[id];
637 
638  //
639  // do the face test on the triangle
640 
641  // was: mark = testNode(V(NV(0)),V(NV(1)),V(NV(2)));
642  // changed to by Gyorgy Fekete. Overall Speedup approx. 2%
643 
644  mark = testNode(id); // was:(indexNode or id);
645 
646  switch (mark)
647  {
648  case fULL:
649  tid = N(id).id_;
650 
651  saveTrixel(tid); //was: plist_->push_back(tid);
652 
653  return mark;
654  case rEJECT:
655  tid = N(id).id_;
656  return mark;
657  default:
658  // if pARTIAL or dONTKNOW, then continue, test children,
659  // but do not reach beyond the leaf nodes.
660  // If Convex is fully contained within one (sWALLOWED),
661  // we can stop looking further in another child
662 
663  // #define NC(n,m) index_->nodes_[(n)].childID_[(m)] // the children n->m
664  // childID = index_->nodes_[id].childID_[0];
665 
666  // Test how many of the four children are rejected?
667  //
668  // [ed:algo]
669  //
670  break;
671  }
672 
673  // NEW NEW algorithm Disabled when enablenew is 0
674  //
675  {
676  childID = indexNode->childID_[0];
677  if (childID != 0)
678  {
679  ////////////// [ed:split]
680  tid = N(id).id_;
681  childID = indexNode->childID_[0];
682  testTrixel(childID);
683  childID = indexNode->childID_[1];
684  testTrixel(childID);
685  childID = indexNode->childID_[2];
686  testTrixel(childID);
687  childID = indexNode->childID_[3];
688  testTrixel(childID);
689  }
690  else /// No children...
691  {
692  if (addlevel_)
693  {
694  testPartial(addlevel_, N(id).id_, V(NV(0)), V(NV(1)), V(NV(2)), 0);
695  }
696  else
697  {
698  saveTrixel(N(id).id_); // was: plist_->push_back(N(id).id_);
699  }
700  }
701  } ///////////////////////////// END OF NEW ALGORITHM returns mark below;
702 
703  /* NEW NEW NEW
704  If rejected, then we return [done]
705  If full, then we list the id (propagate to children) [done]
706 
707  If partial, then we look ahead to see how many children are rejected.
708  But ah, next iteration could benefit from having computed this already.
709 
710  If two children are rejected, then we stop
711  If one or 0 nodes are rejected, then we
712  */
713  return mark;
714 }
715 
716 /////////////TESTPARTIAL//////////////////////////////////
717 // testPartial: test a triangle's subtriangle whether they are partial.
718 // if level is nonzero, recurse.
719 //
720 void RangeConvex::testPartial(size_t level, uint64 id, const SpatialVector &v0, const SpatialVector &v1,
721  const SpatialVector &v2, int PPrev)
722 {
723  uint64 ids[4], id0;
724  SpatialMarkup m[4];
725  int P, F; // count number of partials and fulls
726 
727  SpatialVector w0 = v1 + v2;
728  w0.normalize();
729  SpatialVector w1 = v0 + v2;
730  w1.normalize();
731  SpatialVector w2 = v1 + v0;
732  w2.normalize();
733 
734  ids[0] = id0 = id << 2;
735  ids[1] = id0 + 1;
736  ids[2] = id0 + 2;
737  ids[3] = id0 + 3;
738 
739  m[0] = testNode(v0, w2, w1);
740  m[1] = testNode(v1, w0, w2);
741  m[2] = testNode(v2, w1, w0);
742  m[3] = testNode(w0, w1, w2);
743 
744  F = (m[0] == fULL) + (m[1] == fULL) + (m[2] == fULL) + (m[3] == fULL);
745  P = (m[0] == pARTIAL) + (m[1] == pARTIAL) + (m[2] == pARTIAL) + (m[3] == pARTIAL);
746 
747  //
748  // Several interesting cases for saving this (the parent) trixel.
749  // Case P==4, all four children are partials, so pretend parent is full, we save and return
750  // Case P==3, and F==1, most of the parent is in, so pretend that parent is full again
751  // Case P==2 or 3, but the previous testPartial had three partials, so parent was in an arc
752  // as opposed to previous partials being fewer, so parent was in a tiny corner...
753 
754  if ((level-- <= 0) || ((P == 4) || (F >= 2) || (P == 3 && F == 1) || (P > 1 && PPrev == 3)))
755  {
756  saveTrixel(id);
757  return;
758  }
759  else
760  {
761  // look at each child, see if some need saving;
762  for (int i = 0; i < 4; i++)
763  {
764  if (m[i] == fULL)
765  {
766  saveTrixel(ids[i]);
767  }
768  }
769  // look at the four kids again, for partials
770 
771  if (m[0] == pARTIAL)
772  testPartial(level, ids[0], v0, w2, w1, P);
773  if (m[1] == pARTIAL)
774  testPartial(level, ids[1], v1, w0, w2, P);
775  if (m[2] == pARTIAL)
776  testPartial(level, ids[2], v2, w1, w0, P);
777  if (m[3] == pARTIAL)
778  testPartial(level, ids[3], w0, w1, w2, P);
779  }
780 }
781 
782 /////////////TESTNODE/////////////////////////////////////
783 // testNode: tests the QuadNodes for intersections.
784 //
785 SpatialMarkup RangeConvex::testNode(uint64 id)
786 //uint64 id)
787 // const struct SpatialIndex::QuadNode *indexNode)
788 {
789  const SpatialVector *v0, *v1, *v2;
790  // const struct SpatialIndex::QuadNode &indexNode = index_->nodes_[id];
791  const struct SpatialIndex::QuadNode *indexNode = &index_->nodes_[id];
792  int m;
793  m = indexNode->v_[0];
794  v0 = &index_->vertices_[m]; // the vertex vector m
795 
796  m = indexNode->v_[1];
797  v1 = &index_->vertices_[m];
798 
799  m = indexNode->v_[2];
800  v2 = &index_->vertices_[m];
801 
802  // Start with testing the vertices for the QuadNode with this convex.
803  int vsum = testVertex(v0) + testVertex(v1) + testVertex(v2);
804 
805  SpatialMarkup mark = testTriangle(*v0, *v1, *v2, vsum);
806 
807  // since we cannot play games using the on-the-fly triangles,
808  // substitute dontknow with partial.
809  if (mark == dONTKNOW)
810  mark = pARTIAL;
811 
812  return mark;
813 }
814 SpatialMarkup RangeConvex::testNode(const SpatialVector &v0, const SpatialVector &v1, const SpatialVector &v2)
815 {
816  // Start with testing the vertices for the QuadNode with this convex.
817 
818  int vsum = testVertex(v0) + testVertex(v1) + testVertex(v2);
819 
820  SpatialMarkup mark = testTriangle(v0, v1, v2, vsum);
821 
822  // since we cannot play games using the on-the-fly triangles,
823  // substitute dontknow with partial.
824  if (mark == dONTKNOW)
825  mark = pARTIAL;
826 
827  return mark;
828 }
829 
830 /////////////TESTTRIANGLE//////////////////////////////////
831 // testTriangle: tests a triangle given by 3 vertices if
832 // it intersects the convex.
833 //
834 SpatialMarkup RangeConvex::testTriangle(const SpatialVector &v0, const SpatialVector &v1, const SpatialVector &v2,
835  int vsum)
836 {
837  if (vsum == 1 || vsum == 2)
838  return pARTIAL;
839 
840  // If vsum = 3 then we have all vertices inside the convex.
841  // Now use the following decision tree:
842  //
843  // * If the sign of the convex is pOS or zERO : mark as fULL intersection.
844  //
845  // * Else, test for holes inside the triangle. A 'hole' is a nEG constraint
846  // that has its center inside the triangle. If there is such a hole,
847  // return pARTIAL intersection.
848  //
849  // * Else (no holes, sign nEG or mIXED) test for intersection of nEG
850  // constraints with the edges of the triangle. If there are such,
851  // return pARTIAL intersection.
852  //
853  // * Else return fULL intersection.
854 
855  if (vsum == 3)
856  {
857  if (sign_ == pOS || sign_ == zERO)
858  return fULL;
859  if (testHole(v0, v1, v2))
860  return pARTIAL;
861  if (testEdge(v0, v1, v2))
862  return pARTIAL;
863  return fULL;
864  }
865 
866  // If we have reached that far, we have vsum=0. There is no definite
867  // decision making possible here with our methods, the markup may result
868  // in dONTKNOW. The decision tree is the following:
869  //
870  // * Test with bounding circle of the triangle.
871  //
872  // # If the sign of the convex zERO test with the precalculated
873  // bounding circle of the convex. If it does not intersect with the
874  // triangle's bounding circle, rEJECT.
875  //
876  // # If the sign of the convex is nonZERO: if the bounding circle
877  // lies outside of one of the constraints, rEJECT.
878  //
879  // * Else: there was an intersection with the bounding circle.
880  //
881  // # For zERO convexes, test whether the convex intersects the edges.
882  // If none of the edges of the convex intersects with the edges of
883  // the triangle, we have a rEJECT. Else, pARTIAL.
884  //
885  // # If sign of convex is pOS, or miXED and the smallest constraint does
886  // not intersect the edges and has its center inside the triangle,
887  // return sWALLOW. If no intersection of edges and center outside
888  // triangle, return rEJECT.
889  //
890  // # So the smallest constraint DOES intersect with the edges. If
891  // there is another pOS constraint which does not intersect with
892  // the edges, and has its center outside the triangle, return
893  // rEJECT. If its center is inside the triangle return sWALLOW.
894  // Else, return pARTIAL for pOS and dONTKNOW for mIXED signs.
895  //
896  // * If we are here, return dONTKNOW. There is an intersection with
897  // the bounding circle, none of the vertices is inside the convex and
898  // we have very strange possibilities left for pOS and mIXED signs. For
899  // nEG, i.e. all constraints negative, we also have some complicated
900  // things left for which we cannot test further.
901 
902  if (!testBoundingCircle(v0, v1, v2))
903  return rEJECT;
904 
905  if (sign_ == pOS || sign_ == mIXED || (sign_ == zERO && constraints_.size() <= 2))
906  {
907  // Does the smallest constraint intersect with the edges?
908  if (testEdgeConstraint(v0, v1, v2, 0))
909  {
910  // Is there another positive constraint that does NOT intersect with
911  // the edges?
912  size_t cIndex;
913  if ((cIndex = testOtherPosNone(v0, v1, v2)))
914  {
915  // Does that constraint lie inside or outside of the triangle?
916  if (testConstraintInside(v0, v1, v2, cIndex))
917  return pARTIAL;
918  // Does the triangle lie completely within that constr?
919  else if (constraints_[cIndex].contains(v0))
920  return pARTIAL;
921  else
922  return rEJECT;
923  }
924  else
925  {
926  if (sign_ == pOS || sign_ == zERO)
927  return pARTIAL;
928  else
929  return dONTKNOW;
930  }
931  }
932  else
933  {
934  if (sign_ == pOS || sign_ == zERO)
935  {
936  // Does the smallest lie inside or outside the triangle?
937  if (testConstraintInside(v0, v1, v2, 0))
938  return pARTIAL;
939  else
940  return rEJECT;
941  }
942  else
943  return dONTKNOW;
944  }
945  }
946  else if (sign_ == zERO)
947  {
948  if (corners_.size() > 0 && testEdge0(v0, v1, v2))
949  return pARTIAL;
950  else
951  return rEJECT;
952  }
953  return pARTIAL;
954 }
955 
956 /////////////TESTVERTEX/////////////////////////////////////
957 // testVertex: same as above, but for any spatialvector, no markup speedup
958 //
959 int RangeConvex::testVertex(const SpatialVector &v)
960 {
961  for (size_t i = 0; i < constraints_.size(); i++)
962  if ((constraints_[i].a_ * v) < constraints_[i].d_)
963  return 0;
964 
965  return 1;
966 }
967 int RangeConvex::testVertex(const SpatialVector *v)
968 {
969  for (size_t i = 0; i < constraints_.size(); i++)
970  if ((constraints_[i].a_ * *v) < constraints_[i].d_)
971  return 0;
972 
973  return 1;
974 }
975 
976 /////////////TESTHOLE/////////////////////////////////////
977 // testHole: test for holes. If there is a negative constraint whose center
978 // is inside the triangle, we speak of a hole. Returns true if
979 // found one.
980 //
981 bool RangeConvex::testHole(const SpatialVector &v0, const SpatialVector &v1, const SpatialVector &v2)
982 {
983  bool test = false;
984 
985  for (size_t i = 0; i < constraints_.size(); i++)
986  {
987  if (constraints_[i].sign_ == nEG) // test only 'holes'
988  {
989  // If (a ^ b * c) < 0, vectors abc point clockwise.
990  // -> center c not inside triangle, since vertices a,b are ordered
991  // counter-clockwise. The comparison here is the other way
992  // round because c points into the opposite direction as the hole
993 
994  if (((v0 ^ v1) * constraints_[i].a_) > 0.0L)
995  continue;
996  if (((v1 ^ v2) * constraints_[i].a_) > 0.0L)
997  continue;
998  if (((v2 ^ v0) * constraints_[i].a_) > 0.0L)
999  continue;
1000  test = true;
1001  break;
1002  }
1003  }
1004  return test;
1005 }
1006 
1007 /////////////TESTEDGE0////////////////////////////////////
1008 // testEdge0: test if the edges intersect with the zERO convex.
1009 // The edges are given by the vertex vectors e[0-2]
1010 // All constraints are great circles, so test if their intersect
1011 // with the edges is inside or outside the convex.
1012 // If any edge intersection is inside the convex, return true.
1013 // If all edge intersections are outside, check whether one of
1014 // the corners is inside the triangle. If yes, all of them must be
1015 // inside -> return true.
1016 //
1017 bool RangeConvex::testEdge0(const SpatialVector &v0, const SpatialVector &v1, const SpatialVector &v2)
1018 {
1019  // We have constructed the corners_ array in a certain direction.
1020  // now we can run around the convex, check each side against the 3
1021  // triangle edges. If any of the sides has its intersection INSIDE
1022  // the side, return true. At the end test if a corner lies inside
1023  // (because if we get there, none of the edges intersect, but it
1024  // can be that the convex is fully inside the triangle. so to test
1025  // one single edge is enough)
1026 
1027  struct edgeStruct
1028  {
1029  SpatialVector e; // The half-sphere this edge delimits
1030  float64 l; // length of edge
1031  const SpatialVector *e1; // first end
1032  const SpatialVector *e2; // second end
1033  } edge[3];
1034 
1035  // fill the edge structure for each side of this triangle
1036  edge[0].e = v0 ^ v1;
1037  edge[0].e1 = &v0;
1038  edge[0].e2 = &v1;
1039  edge[1].e = v1 ^ v2;
1040  edge[1].e1 = &v1;
1041  edge[1].e2 = &v2;
1042  edge[2].e = v2 ^ v0;
1043  edge[2].e1 = &v2;
1044  edge[2].e2 = &v0;
1045  edge[0].l = acos(v0 * v1);
1046  edge[1].l = acos(v1 * v2);
1047  edge[2].l = acos(v2 * v0);
1048 
1049  for (size_t i = 0; i < corners_.size(); i++)
1050  {
1051  size_t j = 0;
1052  if (i < corners_.size() - 1)
1053  j = i + 1;
1054  SpatialVector a1;
1055  float64 l1, l2; // lengths of the arcs from intersection to edge corners
1056  float64 cedgelen = acos(corners_[i] * corners_[j]); // length of edge of convex
1057 
1058  // calculate the intersection - all 3 edges
1059  for (size_t iedge = 0; iedge < 3; iedge++)
1060  {
1061  a1 = (edge[iedge].e) ^ (corners_[i] ^ corners_[j]);
1062  a1.normalize();
1063  // if the intersection a1 is inside the edge of the convex,
1064  // its distance to the corners is smaller than the edgelength.
1065  // this test has to be done for both the edge of the convex and
1066  // the edge of the triangle.
1067  for (size_t k = 0; k < 2; k++)
1068  {
1069  l1 = acos(corners_[i] * a1);
1070  l2 = acos(corners_[j] * a1);
1071  if (l1 - cedgelen <= gEpsilon && l2 - cedgelen <= gEpsilon)
1072  {
1073  l1 = acos(*(edge[iedge].e1) * a1);
1074  l2 = acos(*(edge[iedge].e2) * a1);
1075  if (l1 - edge[iedge].l <= gEpsilon && l2 - edge[iedge].l <= gEpsilon)
1076  return true;
1077  }
1078  a1 *= -1.0; // do the same for the other intersection
1079  }
1080  }
1081  }
1082  return testVectorInside(v0, v1, v2, corners_[0]);
1083 }
1084 
1085 /////////////TESTEDGE/////////////////////////////////////
1086 // testEdge: test if edges intersect with constraint. This problem
1087 // is solved by a quadratic equation. Return true if there is
1088 // an intersection.
1089 //
1090 bool RangeConvex::testEdge(const SpatialVector &v0, const SpatialVector &v1, const SpatialVector &v2)
1091 {
1092  for (size_t i = 0; i < constraints_.size(); i++)
1093  {
1094  if (constraints_[i].sign_ == nEG) // test only 'holes'
1095  {
1096  if (eSolve(v0, v1, i))
1097  return true;
1098  if (eSolve(v1, v2, i))
1099  return true;
1100  if (eSolve(v2, v0, i))
1101  return true;
1102  }
1103  }
1104  return false;
1105 }
1106 
1107 /////////////ESOLVE///////////////////////////////////////
1108 // eSolve: solve the quadratic eq. for intersection of an edge with a circle
1109 // constraint. Edge given by grand circle running through v1, v2
1110 // Constraint given by cIndex.
1111 bool RangeConvex::eSolve(const SpatialVector &v1, const SpatialVector &v2, size_t cIndex)
1112 {
1113  float64 gamma1 = v1 * constraints_[cIndex].a_;
1114  float64 gamma2 = v2 * constraints_[cIndex].a_;
1115  float64 mu = v1 * v2;
1116  float64 u2 = (1 - mu) / (1 + mu);
1117 
1118  float64 a = -u2 * (gamma1 + constraints_[cIndex].d_);
1119  float64 b = gamma1 * (u2 - 1) + gamma2 * (u2 + 1);
1120  float64 c = gamma1 - constraints_[cIndex].d_;
1121 
1122  float64 D = b * b - 4 * a * c;
1123 
1124  if (D < 0.0L)
1125  return false; // no intersection
1126 
1127  // calculate roots a'la Numerical Recipes
1128 
1129  float64 q = -0.5L * (b + (SGN(b) * sqrt(D)));
1130 
1131  float64 root1(0), root2(0);
1132  int i = 0;
1133 
1134  if (a > gEpsilon || a < -gEpsilon)
1135  {
1136  root1 = q / a;
1137  i++;
1138  }
1139  if (q > gEpsilon || q < -gEpsilon)
1140  {
1141  root2 = c / q;
1142  i++;
1143  }
1144 
1145  // Check whether the roots lie within [0,1]. If not, the intersection
1146  // is outside the edge.
1147 
1148  if (i == 0)
1149  return false; // no solution
1150  if (root1 >= 0.0L && root1 <= 1.0L)
1151  return true;
1152  if (i == 2 && ((root1 >= 0.0L && root1 <= 1.0L) || (root2 >= 0.0L && root2 <= 1.0L)))
1153  return true;
1154 
1155  return false;
1156 }
1157 
1158 /////////////TESTBOUNDINGCIRCLE///////////////////////////
1159 // testBoundingCircle: test for boundingCircles intersecting with constraint
1160 //
1161 bool RangeConvex::testBoundingCircle(const SpatialVector &v0, const SpatialVector &v1, const SpatialVector &v2)
1162 {
1163  // Set the correct direction: The normal vector to the triangle plane
1164  SpatialVector c = (v1 - v0) ^ (v2 - v1);
1165  c.normalize();
1166 
1167  // Set the correct opening angle: Since the plane cutting out the triangle
1168  // also correctly cuts out the bounding cap of the triangle on the sphere,
1169  // we can take any corner to calculate the opening angle
1170  float64 d = acos(c * v0);
1171 
1172  // for zero convexes, we have calculated a bounding circle for the convex.
1173  // only test with this single bounding circle.
1174 
1175  if (sign_ == zERO)
1176  {
1177  float64 tst;
1178  if (((tst = c * boundingCircle_.a_) < -1.0L + gEpsilon ? gPi : acos(tst)) > (d + boundingCircle_.s_))
1179  return false;
1180  return true;
1181  }
1182 
1183  // for all other convexes, test every constraint. If the bounding
1184  // circle lies completely outside of one of the constraints, reject.
1185  // else, accept.
1186 
1187  size_t i;
1188  for (i = 0; i < constraints_.size(); i++)
1189  {
1190  if (((c * constraints_[i].a_) < -1.0L + gEpsilon ? gPi : acos(c * constraints_[i].a_)) >
1191  (d + constraints_[i].s_))
1192  return false;
1193  }
1194  return true;
1195 }
1196 
1197 /////////////TESTEDGECONSTRAINT///////////////////////////
1198 // testEdgeConstraint: test if edges intersect with a given constraint.
1199 //
1200 bool RangeConvex::testEdgeConstraint(const SpatialVector &v0, const SpatialVector &v1, const SpatialVector &v2,
1201  size_t cIndex)
1202 {
1203  if (eSolve(v0, v1, cIndex))
1204  return true;
1205  if (eSolve(v1, v2, cIndex))
1206  return true;
1207  if (eSolve(v2, v0, cIndex))
1208  return true;
1209  return false;
1210 }
1211 
1212 /////////////TESTOTHERPOSNONE/////////////////////////////
1213 // testOtherPosNone: test for other positive constraints that do
1214 // not intersect with an edge. Return its index
1215 //
1216 size_t RangeConvex::testOtherPosNone(const SpatialVector &v0, const SpatialVector &v1, const SpatialVector &v2)
1217 {
1218  size_t i = 1;
1219  while (i < constraints_.size() && constraints_[i].sign_ == pOS)
1220  {
1221  if (!testEdgeConstraint(v0, v1, v2, i))
1222  return i;
1223  i++;
1224  }
1225  return 0;
1226 }
1227 
1228 /////////////TESTCONSTRAINTINSIDE/////////////////////////
1229 // testConstraintInside: look if a constraint is inside the triangle
1230 //
1231 bool RangeConvex::testConstraintInside(const SpatialVector &v0, const SpatialVector &v1, const SpatialVector &v2,
1232  size_t i)
1233 {
1234  return testVectorInside(v0, v1, v2, constraints_[i].a_);
1235 }
1236 
1237 /////////////TESTVECTORINSIDE////////////////////////////
1238 // testVectorInside: look if a vector is inside the triangle
1239 //
1240 bool RangeConvex::testVectorInside(const SpatialVector &v0, const SpatialVector &v1, const SpatialVector &v2,
1241  SpatialVector &v)
1242 {
1243  // If (a ^ b * c) < 0, vectors abc point clockwise.
1244  // -> center c not inside triangle, since vertices are ordered
1245  // counter-clockwise.
1246 
1247  if ((((v0 ^ v1) * v) < 0) || (((v1 ^ v2) * v) < 0) || (((v2 ^ v0) * v) < 0))
1248  return false;
1249  return true;
1250 }
void normalize()
Normalize vector length to 1.
void simplify()
Simplify the convex, remove redundancies.
The Constraint is really a cone on the sky-sphere.
QStringView level(QStringView ifopt)
KGuiItem test()
void intersect(const SpatialIndex *index, HtmRange *hr)
Intersect with index.
SpatialMarkup testTrixel(uint64 nodeIndex)
void add(SpatialConstraint &)
Add a constraint.
RangeConvex()
Default Constructor.
Definition: RangeConvex.cpp:25
This file is part of the KDE documentation.
Documentation copyright © 1996-2023 The KDE developers.
Generated on Thu Sep 21 2023 04:05:27 by doxygen 1.8.17 written by Dimitri van Heesch, © 1997-2006

KDE's Doxygen guidelines are available online.