# 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
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
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
187void RangeConvex::simplify0()
188{
189 size_t i, j, k;
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;
331 }
332 else
333 {
334 corners_.push_back(vi2);
335 c = c2;
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//
546int 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//
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//
586inline 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
630{
631 SpatialMarkup mark;
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];
683 childID = indexNode->childID_[1];
685 childID = indexNode->childID_[2];
687 childID = indexNode->childID_[3];
689 }
690 else /// No children...
691 {
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//
720void 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//
785SpatialMarkup RangeConvex::testNode(uint64 id)
786//uint64 id)
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
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}
814SpatialMarkup RangeConvex::testNode(const SpatialVector &v0, const SpatialVector &v1, const SpatialVector &v2)
815{
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//
834SpatialMarkup 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//
959int 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}
967int 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//
981bool 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//
1017bool 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//
1090bool 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.
1111bool 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//
1161bool 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//
1200bool 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//
1216size_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//
1231bool 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//
1240bool 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 simplify()
Simplify the convex, remove redundancies.
SpatialMarkup testTrixel(uint64 nodeIndex)
void intersect(const SpatialIndex *index, HtmRange *hr)
Intersect with index.
RangeConvex()
Default Constructor.
The Constraint is really a cone on the sky-sphere.
SpatialIndex is a quad tree of spherical triangles.
SpatialVector is a 3D vector usually living on the surface of the sphere.
void normalize()
Normalize vector length to 1.
QStringView level(QStringView ifopt)
KGuiItem test()
This file is part of the KDE documentation.