1 /*******************************************************************************
2 * bbox.cpp
3 *
4 * This module implements the bounding box calculations.
5 * This file was written by Alexander Enzmann. He wrote the code for
6 * POV-Ray's bounding boxes and generously provided us these enhancements.
7 * The box intersection code was further hacked by Eric Haines to speed it up.
8 *
9 * Just so everyone knows where this came from, the code is VERY heavily
10 * based on the slab code from Mark VandeWettering's MTV raytracer.
11 * POV-Ray is just joining the crowd of admirers of Mark's contribution to
12 * the public domain. [ARE]
13 *
14 * ---------------------------------------------------------------------------
15 * Persistence of Vision Ray Tracer ('POV-Ray') version 3.7.
16 * Copyright 1991-2013 Persistence of Vision Raytracer Pty. Ltd.
17 *
18 * POV-Ray is free software: you can redistribute it and/or modify
19 * it under the terms of the GNU Affero General Public License as
20 * published by the Free Software Foundation, either version 3 of the
21 * License, or (at your option) any later version.
22 *
23 * POV-Ray is distributed in the hope that it will be useful,
24 * but WITHOUT ANY WARRANTY; without even the implied warranty of
25 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
26 * GNU Affero General Public License for more details.
27 *
28 * You should have received a copy of the GNU Affero General Public License
29 * along with this program. If not, see <http://www.gnu.org/licenses/>.
30 * ---------------------------------------------------------------------------
31 * POV-Ray is based on the popular DKB raytracer version 2.12.
32 * DKBTrace was originally written by David K. Buck.
33 * DKBTrace Ver 2.0-2.12 were written by David K. Buck & Aaron A. Collins.
34 * ---------------------------------------------------------------------------
35 * $File: //depot/public/povray/3.x/source/backend/bounding/bbox.cpp $
36 * $Revision: #1 $
37 * $Change: 6069 $
38 * $DateTime: 2013/11/06 11:59:40 $
39 * $Author: chrisc $
40 *******************************************************************************/
41
42 // frame.h must always be the first POV file included (pulls in platform config)
43 #include "backend/frame.h"
44 #include "backend/bounding/bbox.h"
45 #include "backend/scene/objects.h"
46 #include "backend/math/vector.h"
47 #include "backend/math/matrices.h"
48 #include "backend/scene/threaddata.h"
49 #include "base/pov_err.h"
50
51 // this must be the last file included
52 #include "base/povdebug.h"
53
54 namespace pov
55 {
56
57 const int BUNCHING_FACTOR = 4;
58 // Initial number of entries in a priority queue.
59 const int INITIAL_PRIORITY_QUEUE_SIZE = 256;
60
61 BBOX_TREE *create_bbox_node(int size);
62
63 int find_axis(BBOX_TREE **Finite, ptrdiff_t first, ptrdiff_t last);
64 void calc_bbox(BBOX *BBox, BBOX_TREE **Finite, ptrdiff_t first, ptrdiff_t last);
65 void build_area_table(BBOX_TREE **Finite, ptrdiff_t a, ptrdiff_t b, DBL *areas);
66 int sort_and_split(BBOX_TREE **Root, BBOX_TREE **&Finite, size_t *numOfFiniteObjects, ptrdiff_t first, ptrdiff_t last, size_t& maxfinitecount);
67
68 void priority_queue_insert(PriorityQueue& Queue, DBL Depth, BBOX_TREE *Node);
69
PriorityQueue()70 PriorityQueue::PriorityQueue()
71 {
72 QSize = 0;
73 Queue = reinterpret_cast<Qelem *>(POV_MALLOC(INITIAL_PRIORITY_QUEUE_SIZE * sizeof(Qelem), "priority queue"));
74 Max_QSize = INITIAL_PRIORITY_QUEUE_SIZE;
75 }
76
~PriorityQueue()77 PriorityQueue::~PriorityQueue()
78 {
79 POV_FREE(Queue);
80 }
81
Destroy_BBox_Tree(BBOX_TREE * Node)82 void Destroy_BBox_Tree(BBOX_TREE *Node)
83 {
84 if(Node != NULL)
85 {
86 if(Node->Entries > 0)
87 {
88 for(short i = 0; i < Node->Entries; i++)
89 Destroy_BBox_Tree(Node->Node[i]);
90
91 POV_FREE(Node->Node);
92
93 Node->Entries = 0;
94 Node->Node = NULL;
95 }
96
97 POV_FREE(Node);
98 }
99 }
100
Recompute_BBox(BBOX * bbox,const TRANSFORM * trans)101 void Recompute_BBox(BBOX *bbox, const TRANSFORM *trans)
102 {
103 int i;
104 VECTOR lower_left, lengths, corner;
105 VECTOR mins, maxs;
106
107 if(trans == NULL)
108 return;
109
110 Assign_BBox_Vect(lower_left, bbox->Lower_Left);
111 Assign_BBox_Vect(lengths, bbox->Lengths);
112
113 Make_Vector(mins, BOUND_HUGE, BOUND_HUGE, BOUND_HUGE);
114 Make_Vector(maxs, -BOUND_HUGE, -BOUND_HUGE, -BOUND_HUGE);
115
116 for(i = 1; i <= 8; i++)
117 {
118 Assign_Vector(corner, lower_left);
119
120 corner[X] += ((i & 1) ? lengths[X] : 0.0);
121 corner[Y] += ((i & 2) ? lengths[Y] : 0.0);
122 corner[Z] += ((i & 4) ? lengths[Z] : 0.0);
123
124 MTransPoint(corner, corner, trans);
125
126 if(corner[X] < mins[X]) { mins[X] = corner[X]; }
127 if(corner[X] > maxs[X]) { maxs[X] = corner[X]; }
128 if(corner[Y] < mins[Y]) { mins[Y] = corner[Y]; }
129 if(corner[Y] > maxs[Y]) { maxs[Y] = corner[Y]; }
130 if(corner[Z] < mins[Z]) { mins[Z] = corner[Z]; }
131 if(corner[Z] > maxs[Z]) { maxs[Z] = corner[Z]; }
132 }
133
134 // Clip bounding box at the largest allowed bounding box.
135 if(mins[X] < -BOUND_HUGE / 2) { mins[X] = -BOUND_HUGE / 2; }
136 if(mins[Y] < -BOUND_HUGE / 2) { mins[Y] = -BOUND_HUGE / 2; }
137 if(mins[Z] < -BOUND_HUGE / 2) { mins[Z] = -BOUND_HUGE / 2; }
138 if(maxs[X] > BOUND_HUGE / 2) { maxs[X] = BOUND_HUGE / 2; }
139 if(maxs[Y] > BOUND_HUGE / 2) { maxs[Y] = BOUND_HUGE / 2; }
140 if(maxs[Z] > BOUND_HUGE / 2) { maxs[Z] = BOUND_HUGE / 2; }
141
142 Make_BBox_from_min_max(*bbox, mins, maxs);
143 }
144
Recompute_Inverse_BBox(BBOX * bbox,const TRANSFORM * trans)145 void Recompute_Inverse_BBox(BBOX *bbox, const TRANSFORM *trans)
146 {
147 int i;
148 VECTOR lower_left, lengths, corner;
149 VECTOR mins, maxs;
150
151 if(trans == NULL)
152 return;
153
154 Assign_BBox_Vect(lower_left, bbox->Lower_Left);
155 Assign_BBox_Vect(lengths, bbox->Lengths);
156
157 Make_Vector(mins, BOUND_HUGE, BOUND_HUGE, BOUND_HUGE);
158 Make_Vector(maxs, -BOUND_HUGE, -BOUND_HUGE, -BOUND_HUGE);
159
160 for(i = 1; i <= 8; i++)
161 {
162 Assign_Vector(corner, lower_left);
163
164 corner[X] += ((i & 1) ? lengths[X] : 0.0);
165 corner[Y] += ((i & 2) ? lengths[Y] : 0.0);
166 corner[Z] += ((i & 4) ? lengths[Z] : 0.0);
167
168 MInvTransPoint(corner, corner, trans);
169
170 if(corner[X] < mins[X]) { mins[X] = corner[X]; }
171 if(corner[X] > maxs[X]) { maxs[X] = corner[X]; }
172 if(corner[Y] < mins[Y]) { mins[Y] = corner[Y]; }
173 if(corner[Y] > maxs[Y]) { maxs[Y] = corner[Y]; }
174 if(corner[Z] < mins[Z]) { mins[Z] = corner[Z]; }
175 if(corner[Z] > maxs[Z]) { maxs[Z] = corner[Z]; }
176 }
177
178 // Clip bounding box at the largest allowed bounding box.
179 if(mins[X] < -BOUND_HUGE / 2) { mins[X] = -BOUND_HUGE / 2; }
180 if(mins[Y] < -BOUND_HUGE / 2) { mins[Y] = -BOUND_HUGE / 2; }
181 if(mins[Z] < -BOUND_HUGE / 2) { mins[Z] = -BOUND_HUGE / 2; }
182 if(maxs[X] > BOUND_HUGE / 2) { maxs[X] = BOUND_HUGE / 2; }
183 if(maxs[Y] > BOUND_HUGE / 2) { maxs[Y] = BOUND_HUGE / 2; }
184 if(maxs[Z] > BOUND_HUGE / 2) { maxs[Z] = BOUND_HUGE / 2; }
185
186 Make_BBox_from_min_max(*bbox, mins, maxs);
187 }
188
189 // Create a bounding box hierarchy from a given list of finite and
190 // infinite elements. Each element consists of
191 //
192 // - an infinite flag
193 // - a bounding box enclosing the element
194 // - a pointer to the structure representing the element (e.g an object)
Build_BBox_Tree(BBOX_TREE ** Root,size_t numOfFiniteObjects,BBOX_TREE ** & Finite,size_t numOfInfiniteObjects,BBOX_TREE ** Infinite,size_t & maxfinitecount)195 void Build_BBox_Tree(BBOX_TREE **Root, size_t numOfFiniteObjects, BBOX_TREE **&Finite, size_t numOfInfiniteObjects, BBOX_TREE **Infinite, size_t& maxfinitecount)
196 {
197 ptrdiff_t low, high;
198 BBOX_TREE *cd, *root;
199
200 // This is a resonable guess at the number of finites needed.
201 // This array will be reallocated as needed if it isn't.
202 maxfinitecount = 2 * numOfFiniteObjects;
203
204 // Now do a sort on the objects, with the end result being
205 // a tree of objects sorted along the x, y, and z axes.
206 if(numOfFiniteObjects > 0)
207 {
208 low = 0;
209 high = numOfFiniteObjects;
210
211 while(sort_and_split(Root, Finite, &numOfFiniteObjects, low, high, maxfinitecount) == 0)
212 {
213 low = high;
214 high = numOfFiniteObjects;
215 }
216
217 // Move infinite objects in the first leaf of Root.
218 if(numOfInfiniteObjects > 0)
219 {
220 root = *Root;
221 root->Node = reinterpret_cast<BBOX_TREE **>(POV_REALLOC(root->Node, (root->Entries + 1) * sizeof(BBOX_TREE *), "composite"));
222 POV_MEMMOVE(&(root->Node[1]), &(root->Node[0]), root->Entries * sizeof(BBOX_TREE *));
223 root->Entries++;
224 cd = create_bbox_node(numOfInfiniteObjects);
225 for(size_t i = 0; i < numOfInfiniteObjects; i++)
226 cd->Node[i] = Infinite[i];
227
228 calc_bbox(&(cd->BBox), Infinite, 0, numOfInfiniteObjects);
229 root->Node[0] = cd;
230 calc_bbox(&(root->BBox), root->Node, 0, root->Entries);
231
232 // Root and first node are infinite.
233 root->Infinite = true;
234 root->Node[0]->Infinite = true;
235 }
236 }
237 else
238 {
239 // There are no finite objects and no Root was created.
240 // Create it now and put all infinite objects into it.
241
242 if(numOfInfiniteObjects > 0)
243 {
244 cd = create_bbox_node(numOfInfiniteObjects);
245 for(size_t i = 0; i < numOfInfiniteObjects; i++)
246 cd->Node[i] = Infinite[i];
247 calc_bbox(&(cd->BBox), Infinite, 0, numOfInfiniteObjects);
248 *Root = cd;
249 (*Root)->Infinite = true;
250 }
251 }
252 }
253
Build_Bounding_Slabs(BBOX_TREE ** Root,vector<ObjectPtr> & objects,unsigned int & numberOfFiniteObjects,unsigned int & numberOfInfiniteObjects,unsigned int & numberOfLightSources)254 void Build_Bounding_Slabs(BBOX_TREE **Root, vector<ObjectPtr>& objects, unsigned int& numberOfFiniteObjects, unsigned int& numberOfInfiniteObjects, unsigned int& numberOfLightSources)
255 {
256 ptrdiff_t iFinite, iInfinite;
257 BBOX_TREE **Finite, **Infinite;
258 ObjectPtr Temp;
259 size_t maxfinitecount = 0;
260
261 // Count frame level and infinite objects.
262 numberOfFiniteObjects = numberOfInfiniteObjects = numberOfLightSources = 0;
263
264 for(vector<ObjectPtr>::iterator i(objects.begin()); i != objects.end(); i++)
265 {
266 if((*i)->Type & LIGHT_SOURCE_OBJECT)
267 {
268 if((reinterpret_cast<LightSource *>(*i))->children.size() > 0)
269 {
270 Temp = (reinterpret_cast<LightSource *>(*i))->children[0];
271 numberOfLightSources++;
272 }
273 else
274 Temp = NULL;
275 }
276 else
277 Temp = (*i);
278
279 if(Temp != NULL)
280 {
281 if(Test_Flag(Temp, INFINITE_FLAG))
282 numberOfInfiniteObjects++;
283 else
284 numberOfFiniteObjects++;
285 }
286 }
287
288 // If bounding boxes aren't used we can return.
289 if(numberOfFiniteObjects + numberOfInfiniteObjects < 1)
290 return;
291
292 // This is a resonable guess at the number of finites needed.
293 // This array will be reallocated as needed if it isn't.
294 maxfinitecount = 2 * numberOfFiniteObjects;
295
296 // Now allocate an array to hold references to these finites and
297 // any new composite objects we may generate.
298 Finite = Infinite = NULL;
299
300 if(numberOfFiniteObjects > 0)
301 Finite = reinterpret_cast<BBOX_TREE **>(POV_MALLOC(maxfinitecount*sizeof(BBOX_TREE *), "bounding boxes"));
302
303 // Create array to hold pointers to infinite objects.
304 if(numberOfInfiniteObjects > 0)
305 Infinite = reinterpret_cast<BBOX_TREE **>(POV_MALLOC(numberOfInfiniteObjects*sizeof(BBOX_TREE *), "bounding boxes"));
306
307 // Init lists.
308 for(int i = 0; i < numberOfFiniteObjects; i++)
309 Finite[i] = create_bbox_node(0);
310
311 for(int i = 0; i < numberOfInfiniteObjects; i++)
312 Infinite[i] = create_bbox_node(0);
313
314 // Set up finite and infinite object lists.
315 iFinite = iInfinite = 0;
316
317 for(vector<ObjectPtr>::iterator i(objects.begin()); i != objects.end(); i++)
318 {
319 if((*i)->Type & LIGHT_SOURCE_OBJECT)
320 {
321 if((reinterpret_cast<LightSource *>(*i))->children.size() > 0)
322 Temp = (reinterpret_cast<LightSource *>(*i))->children[0];
323 else
324 Temp = NULL;
325 }
326 else
327 Temp = (*i);
328
329 if(Temp != NULL)
330 {
331 // Add object to the appropriate list.
332 if(Test_Flag(Temp, INFINITE_FLAG))
333 {
334 Infinite[iInfinite]->Infinite = true;
335 Infinite[iInfinite]->BBox = Temp->BBox;
336 Infinite[iInfinite]->Node = reinterpret_cast<BBOX_TREE **>(Temp);
337
338 iInfinite++;
339 }
340 else
341 {
342 Finite[iFinite]->BBox = Temp->BBox;
343 Finite[iFinite]->Node = reinterpret_cast<BBOX_TREE **>(Temp);
344
345 iFinite++;
346 }
347 }
348 }
349
350 // Now build the bounding box tree.
351 Build_BBox_Tree(Root, numberOfFiniteObjects, Finite, numberOfInfiniteObjects, Infinite, maxfinitecount);
352
353 // Get rid of the Finite and Infinite arrays and just use Root.
354 if(Finite != NULL)
355 POV_FREE(Finite);
356
357 if(Infinite != NULL)
358 POV_FREE(Infinite);
359 }
360
Intersect_BBox_Tree(PriorityQueue & pqueue,const BBOX_TREE * Root,const Ray & ray,Intersection * Best_Intersection,TraceThreadData * Thread)361 bool Intersect_BBox_Tree(PriorityQueue& pqueue, const BBOX_TREE *Root, const Ray& ray, Intersection *Best_Intersection, TraceThreadData *Thread)
362 {
363 int i, found;
364 DBL Depth;
365 const BBOX_TREE *Node;
366 Intersection New_Intersection;
367
368 // Create the direction vectors for this ray.
369 Rayinfo rayinfo(ray);
370
371 // Start with an empty priority queue.
372 pqueue.QSize = 0;
373 New_Intersection.Object = NULL;
374 found = false;
375
376 // Check top node.
377 Check_And_Enqueue(pqueue, Root, &Root->BBox, &rayinfo, Thread);
378
379 // Check elements in the priority queue.
380 while(pqueue.QSize != 0)
381 {
382 Priority_Queue_Delete(pqueue, &Depth, &Node);
383
384 // If current intersection is larger than the best intersection found
385 // so far our task is finished, because all other bounding boxes in
386 // the priority queue are further away.
387 if(Depth > Best_Intersection->Depth)
388 break;
389
390 // Check current node.
391 if(Node->Entries)
392 {
393 // This is a node containing leaves to be checked.
394 for (i = 0; i < Node->Entries; i++)
395 Check_And_Enqueue(pqueue, Node->Node[i], &Node->Node[i]->BBox, &rayinfo, Thread);
396 }
397 else
398 {
399 // This is a leaf so test contained object.
400 if(Find_Intersection(&New_Intersection, reinterpret_cast<ObjectPtr>(Node->Node), ray, Thread))
401 {
402 if(New_Intersection.Depth < Best_Intersection->Depth)
403 {
404 *Best_Intersection = New_Intersection;
405 found = true;
406 }
407 }
408 }
409 }
410
411 return (found);
412 }
413
Intersect_BBox_Tree(PriorityQueue & pqueue,const BBOX_TREE * Root,const Ray & ray,Intersection * Best_Intersection,const RayObjectCondition & precondition,const RayObjectCondition & postcondition,TraceThreadData * Thread)414 bool Intersect_BBox_Tree(PriorityQueue& pqueue, const BBOX_TREE *Root, const Ray& ray, Intersection *Best_Intersection, const RayObjectCondition& precondition, const RayObjectCondition& postcondition, TraceThreadData *Thread)
415 {
416 int i, found;
417 DBL Depth;
418 const BBOX_TREE *Node;
419 Intersection New_Intersection;
420
421 // Create the direction vectors for this ray.
422 Rayinfo rayinfo(ray);
423
424 // Start with an empty priority queue.
425 pqueue.QSize = 0;
426 New_Intersection.Object = NULL;
427 found = false;
428
429 // Check top node.
430 Check_And_Enqueue(pqueue, Root, &Root->BBox, &rayinfo, Thread);
431
432 // Check elements in the priority queue.
433 while(pqueue.QSize != 0)
434 {
435 Priority_Queue_Delete(pqueue, &Depth, &Node);
436
437 // If current intersection is larger than the best intersection found
438 // so far our task is finished, because all other bounding boxes in
439 // the priority queue are further away.
440 if(Depth > Best_Intersection->Depth)
441 break;
442
443 // Check current node.
444 if(Node->Entries)
445 {
446 // This is a node containing leaves to be checked.
447 for (i = 0; i < Node->Entries; i++)
448 Check_And_Enqueue(pqueue, Node->Node[i], &Node->Node[i]->BBox, &rayinfo, Thread);
449 }
450 else
451 {
452 if(precondition(ray, reinterpret_cast<ObjectPtr>(Node->Node), 0.0) == true)
453 {
454 // This is a leaf so test contained object.
455 if(Find_Intersection(&New_Intersection, reinterpret_cast<ObjectPtr>(Node->Node), ray, postcondition, Thread))
456 {
457 if(New_Intersection.Depth < Best_Intersection->Depth)
458 {
459 *Best_Intersection = New_Intersection;
460 found = true;
461 }
462 }
463 }
464 }
465 }
466
467 return (found);
468 }
469
priority_queue_insert(PriorityQueue & Queue,DBL Depth,const BBOX_TREE * Node)470 void priority_queue_insert(PriorityQueue& Queue, DBL Depth, const BBOX_TREE *Node)
471 {
472 unsigned size;
473 int i;
474 //PriorityQueue::Qelem tmp;
475 PriorityQueue::Qelem *List;
476
477 Queue.QSize++;
478
479 size = Queue.QSize;
480
481 /* Reallocate priority queue if it's too small. */
482
483 if (size >= Queue.Max_QSize)
484 {
485 /*
486 if (size >= INT_MAX/2)
487 {
488 // TODO FIXME Error("Priority queue overflow.");
489 }
490 */
491
492 Queue.Max_QSize *= 2;
493
494 Queue.Queue = reinterpret_cast<PriorityQueue::Qelem *>(POV_REALLOC(Queue.Queue, Queue.Max_QSize*sizeof(PriorityQueue::Qelem), "priority queue"));
495 }
496
497 List = Queue.Queue;
498
499 /*
500 ***
501 List[size].depth = Depth;
502 List[size].node = Node;
503
504 i = size;
505
506 while (i > 1 && List[i].depth < List[i / 2].depth)
507 {
508 tmp = List[i];
509
510 List[i] = List[i / 2];
511
512 List[i / 2] = tmp;
513
514 i = i / 2;
515 }
516 ***
517 */
518
519 i = size;
520 while(i > 1 && Depth < List[i/2].depth)
521 {
522 List[i] = List[i/2];
523 i /= 2;
524 }
525 List[i].depth = Depth;
526 List[i].node = Node;
527 }
528
529 // Get an element from the priority queue.
530 // NOTE: This element will always be the one closest to the ray origin.
Priority_Queue_Delete(PriorityQueue & Queue,DBL * Depth,const BBOX_TREE ** Node)531 void Priority_Queue_Delete(PriorityQueue& Queue, DBL *Depth, const BBOX_TREE **Node)
532 {
533 PriorityQueue::Qelem tmp;
534 PriorityQueue::Qelem *List;
535 int i, j;
536 unsigned size;
537
538 if (Queue.QSize == 0)
539 {
540 // TODO FIXME Error("priority queue is empty.");
541 }
542
543 List = Queue.Queue;
544
545 *Depth = List[1].depth;
546 *Node = List[1].node;
547
548 List[1] = List[Queue.QSize];
549
550 Queue.QSize--;
551
552 size = Queue.QSize;
553
554 i = 1;
555
556 while (2 * i <= (int)size)
557 {
558 if (2 * i == (int)size)
559 {
560 j = 2 * i;
561 }
562 else
563 {
564 if (List[2*i].depth < List[2*i+1].depth)
565 {
566 j = 2 * i;
567 }
568 else
569 {
570 j = 2 * i + 1;
571 }
572 }
573
574 if (List[i].depth > List[j].depth)
575 {
576 tmp = List[i];
577
578 List[i] = List[j];
579
580 List[j] = tmp;
581
582 i = j;
583 }
584 else
585 {
586 break;
587 }
588 }
589 }
590
Check_And_Enqueue(PriorityQueue & Queue,const BBOX_TREE * Node,const BBOX * BBox,const Rayinfo * rayinfo,TraceThreadData * Thread)591 void Check_And_Enqueue(PriorityQueue& Queue, const BBOX_TREE *Node, const BBOX *BBox, const Rayinfo *rayinfo, TraceThreadData *Thread)
592 {
593 DBL tmin, tmax;
594 DBL dmin, dmax;
595
596 if(Node->Infinite == false)
597 {
598 Thread->Stats()[nChecked]++;
599
600 if(rayinfo->nonzero[X])
601 {
602 if (rayinfo->positive[X])
603 {
604 dmin = (BBox->Lower_Left[X] - rayinfo->slab_num[X]) * rayinfo->slab_den[X];
605 dmax = dmin + (BBox->Lengths[X] * rayinfo->slab_den[X]);
606 if(dmax < EPSILON)
607 return;
608 }
609 else
610 {
611 dmax = (BBox->Lower_Left[X] - rayinfo->slab_num[X]) * rayinfo->slab_den[X];
612
613 if(dmax < EPSILON)
614 return;
615
616 dmin = dmax + (BBox->Lengths[X] * rayinfo->slab_den[X]);
617 }
618
619 if(dmin > dmax)
620 return;
621 }
622 else
623 {
624 if((rayinfo->slab_num[X] < BBox->Lower_Left[X]) ||
625 (rayinfo->slab_num[X] > BBox->Lengths[X] + BBox->Lower_Left[X]))
626 return;
627
628 dmin = -BOUND_HUGE;
629 dmax = BOUND_HUGE;
630 }
631
632 if(rayinfo->nonzero[Y])
633 {
634 if(rayinfo->positive[Y])
635 {
636 tmin = (BBox->Lower_Left[Y] - rayinfo->slab_num[Y]) * rayinfo->slab_den[Y];
637 tmax = tmin + (BBox->Lengths[Y] * rayinfo->slab_den[Y]);
638 }
639 else
640 {
641 tmax = (BBox->Lower_Left[Y] - rayinfo->slab_num[Y]) * rayinfo->slab_den[Y];
642 tmin = tmax + (BBox->Lengths[Y] * rayinfo->slab_den[Y]);
643 }
644
645 // Unwrap the logic - do the dmin and dmax checks only when tmin and
646 // tmax actually affect anything, also try to escape ASAP. Better
647 // yet, fold the logic below into the two branches above so as to
648 // compute only what is needed.
649
650 // You might even try tmax < EPSILON first (instead of second) for an
651 // early quick out.
652
653 if(tmax < dmax)
654 {
655 if(tmax < EPSILON)
656 return;
657
658 // check bbox only if tmax changes dmax
659
660 if(tmin > dmin)
661 {
662 if(tmin > tmax)
663 return;
664
665 // do this last in case it's not needed!
666 dmin = tmin;
667 }
668 else if(dmin > tmax)
669 return;
670
671 // do this last in case it's not needed!
672 dmax = tmax;
673 }
674 else if(tmin > dmin)
675 {
676 if(tmin > dmax)
677 return;
678
679 // do this last in case it's not needed!
680 dmin = tmin;
681 }
682 }
683 else if((rayinfo->slab_num[Y] < BBox->Lower_Left[Y]) ||
684 (rayinfo->slab_num[Y] > BBox->Lengths[Y] + BBox->Lower_Left[Y]))
685 return;
686
687 if(rayinfo->nonzero[Z])
688 {
689 if(rayinfo->positive[Z])
690 {
691 tmin = (BBox->Lower_Left[Z] - rayinfo->slab_num[Z]) * rayinfo->slab_den[Z];
692 tmax = tmin + (BBox->Lengths[Z] * rayinfo->slab_den[Z]);
693 }
694 else
695 {
696 tmax = (BBox->Lower_Left[Z] - rayinfo->slab_num[Z]) * rayinfo->slab_den[Z];
697 tmin = tmax + (BBox->Lengths[Z] * rayinfo->slab_den[Z]);
698 }
699
700 if(tmax < dmax)
701 {
702 if(tmax < EPSILON)
703 return;
704
705 // check bbox only if tmax changes dmax
706 if(tmin > dmin)
707 {
708 if(tmin > tmax)
709 return;
710
711 // do this last in case it's not needed!
712 dmin = tmin;
713 }
714 else if(dmin > tmax)
715 return;
716 }
717 else if(tmin > dmin)
718 {
719 if(tmin > dmax)
720 return;
721
722 // do this last in case it's not needed!
723 dmin = tmin;
724 }
725 }
726 else
727 if((rayinfo->slab_num[Z] < BBox->Lower_Left[Z]) || (rayinfo->slab_num[Z] > BBox->Lengths[Z] + BBox->Lower_Left[Z]))
728 return;
729
730 Thread->Stats()[nEnqueued]++;
731 }
732 else
733 // Set intersection depth to -Max_Distance.
734 dmin = -MAX_DISTANCE;
735
736 priority_queue_insert(Queue, dmin, Node);
737 }
738
create_bbox_node(int size)739 BBOX_TREE *create_bbox_node(int size)
740 {
741 BBOX_TREE *New;
742
743 New = reinterpret_cast<BBOX_TREE *>(POV_MALLOC(sizeof(BBOX_TREE), "bounding box node"));
744
745 New->Infinite = false;
746 New->Entries = size;
747
748 Make_BBox(New->BBox, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0);
749
750 if(size)
751 New->Node = reinterpret_cast<BBOX_TREE **>(POV_MALLOC(size*sizeof(BBOX_TREE *), "bounding box node"));
752 else
753 New->Node = NULL;
754
755 return (New);
756 }
757
758 template<int Axis>
compboxes(const void * in_a,const void * in_b)759 int CDECL compboxes(const void *in_a, const void *in_b)
760 {
761 const BBOX *a, *b;
762 BBOX_VAL am, bm;
763 typedef const BBOX_TREE *CONST_BBOX_TREE_PTR;
764
765 a = &((*reinterpret_cast<const CONST_BBOX_TREE_PTR *>(in_a))->BBox);
766 b = &((*reinterpret_cast<const CONST_BBOX_TREE_PTR *>(in_b))->BBox);
767
768 am = 2.0 * a->Lower_Left[Axis] + a->Lengths[Axis];
769 bm = 2.0 * b->Lower_Left[Axis] + b->Lengths[Axis];
770
771 if(am < bm)
772 return -1;
773 else
774 {
775 if(am == bm)
776 return 0;
777 else
778 return 1;
779 }
780 }
781
find_axis(BBOX_TREE ** Finite,ptrdiff_t first,ptrdiff_t last)782 int find_axis(BBOX_TREE **Finite, ptrdiff_t first, ptrdiff_t last)
783 {
784 int which = X;
785 ptrdiff_t i;
786 DBL e, d = -BOUND_HUGE;
787 VECTOR mins, maxs;
788 BBOX *bbox;
789
790 Make_Vector(mins, BOUND_HUGE, BOUND_HUGE, BOUND_HUGE);
791 Make_Vector(maxs, -BOUND_HUGE, -BOUND_HUGE, -BOUND_HUGE);
792
793 for(i = first; i < last; i++)
794 {
795 bbox = &(Finite[i]->BBox);
796
797 if(bbox->Lower_Left[X] < mins[X])
798 mins[X] = bbox->Lower_Left[X];
799
800 if(bbox->Lower_Left[X] + bbox->Lengths[X] > maxs[X])
801 maxs[X] = bbox->Lower_Left[X];
802
803 if(bbox->Lower_Left[Y] < mins[Y])
804 mins[Y] = bbox->Lower_Left[Y];
805
806 if(bbox->Lower_Left[Y] + bbox->Lengths[Y] > maxs[Y])
807 maxs[Y] = bbox->Lower_Left[Y];
808
809 if(bbox->Lower_Left[Z] < mins[Z])
810 mins[Z] = bbox->Lower_Left[Z];
811
812 if(bbox->Lower_Left[Z] + bbox->Lengths[Z] > maxs[Z])
813 maxs[Z] = bbox->Lower_Left[Z];
814 }
815
816 e = maxs[X] - mins[X];
817
818 if(e > d)
819 {
820 d = e;
821 which = X;
822 }
823
824 e = maxs[Y] - mins[Y];
825
826 if(e > d)
827 {
828 d = e;
829 which = Y;
830 }
831
832 e = maxs[Z] - mins[Z];
833
834 if(e > d)
835 which = Z;
836
837 return (which);
838 }
839
calc_bbox(BBOX * BBox,BBOX_TREE ** Finite,ptrdiff_t first,ptrdiff_t last)840 void calc_bbox(BBOX *BBox, BBOX_TREE **Finite, ptrdiff_t first, ptrdiff_t last)
841 {
842 ptrdiff_t i;
843 DBL tmin, tmax;
844 VECTOR bmin, bmax;
845 BBOX *bbox;
846
847 Make_Vector(bmin, BOUND_HUGE, BOUND_HUGE, BOUND_HUGE);
848 Make_Vector(bmax, -BOUND_HUGE, -BOUND_HUGE, -BOUND_HUGE);
849
850 for(i = first; i < last; i++)
851 {
852 bbox = &(Finite[i]->BBox);
853
854 tmin = bbox->Lower_Left[X];
855 tmax = tmin + bbox->Lengths[X];
856
857 if(tmin < bmin[X]) { bmin[X] = tmin; }
858 if(tmax > bmax[X]) { bmax[X] = tmax; }
859
860 tmin = bbox->Lower_Left[Y];
861 tmax = tmin + bbox->Lengths[Y];
862
863 if(tmin < bmin[Y]) { bmin[Y] = tmin; }
864 if(tmax > bmax[Y]) { bmax[Y] = tmax; }
865
866 tmin = bbox->Lower_Left[Z];
867 tmax = tmin + bbox->Lengths[Z];
868
869 if(tmin < bmin[Z]) { bmin[Z] = tmin; }
870 if(tmax > bmax[Z]) { bmax[Z] = tmax; }
871 }
872
873 Make_BBox_from_min_max(*BBox, bmin, bmax);
874 }
875
build_area_table(BBOX_TREE ** Finite,ptrdiff_t a,ptrdiff_t b,DBL * areas)876 void build_area_table(BBOX_TREE **Finite, ptrdiff_t a, ptrdiff_t b, DBL *areas)
877 {
878 ptrdiff_t i, imin, dir;
879 DBL tmin, tmax;
880 VECTOR bmin, bmax, len;
881 BBOX *bbox;
882
883 if (a < b)
884 {
885 imin = a; dir = 1;
886 }
887 else
888 {
889 imin = b; dir = -1;
890 }
891
892 Make_Vector(bmin, BOUND_HUGE, BOUND_HUGE, BOUND_HUGE);
893 Make_Vector(bmax, -BOUND_HUGE, -BOUND_HUGE, -BOUND_HUGE);
894
895 for(i = a; i != (b + dir); i += dir)
896 {
897 bbox = &(Finite[i]->BBox);
898
899 tmin = bbox->Lower_Left[X];
900 tmax = tmin + bbox->Lengths[X];
901
902 if (tmin < bmin[X]) { bmin[X] = tmin; }
903 if (tmax > bmax[X]) { bmax[X] = tmax; }
904
905 tmin = bbox->Lower_Left[Y];
906 tmax = tmin + bbox->Lengths[Y];
907
908 if (tmin < bmin[Y]) { bmin[Y] = tmin; }
909 if (tmax > bmax[Y]) { bmax[Y] = tmax; }
910
911 tmin = bbox->Lower_Left[Z];
912 tmax = tmin + bbox->Lengths[Z];
913
914 if (tmin < bmin[Z]) { bmin[Z] = tmin; }
915 if (tmax > bmax[Z]) { bmax[Z] = tmax; }
916
917 VSub(len, bmax, bmin);
918
919 areas[i - imin] = len[X] * (len[Y] + len[Z]) + len[Y] * len[Z];
920 }
921 }
922
sort_and_split(BBOX_TREE ** Root,BBOX_TREE ** & Finite,size_t * numOfFiniteObjects,ptrdiff_t first,ptrdiff_t last,size_t & maxfinitecount)923 int sort_and_split(BBOX_TREE **Root, BBOX_TREE **&Finite, size_t *numOfFiniteObjects, ptrdiff_t first, ptrdiff_t last, size_t& maxfinitecount)
924 {
925 BBOX_TREE *cd;
926 ptrdiff_t size, i, best_loc;
927 DBL *area_left, *area_right;
928 DBL best_index, new_index;
929
930 int Axis = find_axis(Finite, first, last);
931 size = last - first;
932 if(size <= 0)
933 return (1);
934
935 // Actually, we could do this faster in several ways. We could use a
936 // logn algorithm to find the median along the given axis, and then a
937 // linear algorithm to partition along the axis. Oh well.
938
939 switch(Axis)
940 {
941 case X:
942 QSORT(reinterpret_cast<void *>(&Finite[first]), size, sizeof(BBOX_TREE *), compboxes<X>);
943 break;
944 case Y:
945 QSORT(reinterpret_cast<void *>(&Finite[first]), size, sizeof(BBOX_TREE *), compboxes<Y>);
946 break;
947 case Z:
948 QSORT(reinterpret_cast<void *>(&Finite[first]), size, sizeof(BBOX_TREE *), compboxes<Z>);
949 break;
950 }
951
952 // area_left[] and area_right[] hold the surface areas of the bounding
953 // boxes to the left and right of any given point. E.g. area_left[i] holds
954 // the surface area of the bounding box containing Finite 0 through i and
955 // area_right[i] holds the surface area of the box containing Finite
956 // i through size-1.
957
958 area_left = reinterpret_cast<DBL *>(POV_MALLOC(size * sizeof(DBL), "bounding boxes"));
959 area_right = reinterpret_cast<DBL *>(POV_MALLOC(size * sizeof(DBL), "bounding boxes"));
960
961 // Precalculate the areas for speed.
962 build_area_table(Finite, first, last - 1, area_left);
963 build_area_table(Finite, last - 1, first, area_right);
964 best_index = area_right[0] * (size - 3.0);
965 best_loc = -1;
966
967 // Find the most effective point to split. The best location will be
968 // the one that minimizes the function N1*A1 + N2*A2 where N1 and N2
969 // are the number of objects in the two groups and A1 and A2 are the
970 // surface areas of the bounding boxes of the two groups.
971
972 for(i = 0; i < size - 1; i++)
973 {
974 new_index = (i + 1) * area_left[i] + (size - 1 - i) * area_right[i + 1];
975
976 if(new_index < best_index)
977 {
978 best_index = new_index;
979 best_loc = i + first;
980 }
981 }
982
983 POV_FREE(area_left);
984 POV_FREE(area_right);
985
986 // Stop splitting if the BUNCHING_FACTOR is reached or
987 // if splitting stops being effective.
988 if((size <= BUNCHING_FACTOR) || (best_loc < 0))
989 {
990 cd = create_bbox_node(size);
991
992 for(i = 0; i < size; i++)
993 cd->Node[i] = Finite[first+i];
994
995 calc_bbox(&(cd->BBox), Finite, first, last);
996 *Root = cd;
997 if(*numOfFiniteObjects >= maxfinitecount)
998 {
999 // Prim array overrun, increase array by 50%.
1000 maxfinitecount = 1.5 * maxfinitecount;
1001
1002 // For debugging only.
1003 // TODO MESSAGE Debug_Info("Reallocing Finite to %d\n", maxfinitecount);
1004 Finite = reinterpret_cast<BBOX_TREE **>(POV_REALLOC(Finite, maxfinitecount * sizeof(BBOX_TREE *), "bounding boxes"));
1005 }
1006
1007 Finite[*numOfFiniteObjects] = cd;
1008 (*numOfFiniteObjects)++;
1009
1010 return (1);
1011 }
1012
1013 sort_and_split(Root, Finite, numOfFiniteObjects, first, best_loc + 1, maxfinitecount);
1014 sort_and_split(Root, Finite, numOfFiniteObjects, best_loc + 1, last, maxfinitecount);
1015
1016 return (0);
1017 }
1018
1019 }
1020