1 /*****************************************************************************
2 * bsphere.cpp
3 *
4 * This module implements the bounding sphere calculations.
5 *
6 * from Persistence of Vision(tm) Ray Tracer version 3.6.
7 * Copyright 1991-2003 Persistence of Vision Team
8 * Copyright 2003-2004 Persistence of Vision Raytracer Pty. Ltd.
9 *---------------------------------------------------------------------------
10 * NOTICE: This source code file is provided so that users may experiment
11 * with enhancements to POV-Ray and to port the software to platforms other
12 * than those supported by the POV-Ray developers. There are strict rules
13 * regarding how you are permitted to use this file. These rules are contained
14 * in the distribution and derivative versions licenses which should have been
15 * provided with this file.
16 *
17 * These licences may be found online, linked from the end-user license
18 * agreement that is located at http://www.povray.org/povlegal.html
19 *---------------------------------------------------------------------------
20 * This program is based on the popular DKB raytracer version 2.12.
21 * DKBTrace was originally written by David K. Buck.
22 * DKBTrace Ver 2.0-2.12 were written by David K. Buck & Aaron A. Collins.
23 *---------------------------------------------------------------------------
24 * $File: //depot/povray/3.6-release/source/bsphere.cpp $
25 * $Revision: #3 $
26 * $Change: 3032 $
27 * $DateTime: 2004/08/02 18:43:41 $
28 * $Author: chrisc $
29 * $Log$
30 *****************************************************************************/
31
32 #include "frame.h"
33 #include "vector.h"
34 #include "bsphere.h"
35 #include "povray.h"
36
37 #include <algorithm>
38
39 BEGIN_POV_NAMESPACE
40
41 /*****************************************************************************
42 * Local preprocessor defines
43 ******************************************************************************/
44
45 const int BRANCHING_FACTOR = 4;
46
47
48
49 /*****************************************************************************
50 * Local typedefs
51 ******************************************************************************/
52
53
54
55 /*****************************************************************************
56 * Local variables
57 ******************************************************************************/
58
59 static int maxelements, Axis; // GLOBAL VARIABLE
60
61
62
63 /*****************************************************************************
64 * Static functions
65 ******************************************************************************/
66
67 static void merge_spheres (VECTOR C, DBL *r, VECTOR C1, DBL r1, VECTOR C2, DBL r2);
68 static void recompute_bound (BSPHERE_TREE *Node);
69 static int CDECL comp_elements (const void *in_a, const void *in_b);
70
71 /*****************************************************************************
72 *
73 * FUNCTION
74 *
75 * merge_spheres
76 *
77 * INPUT
78 *
79 * C, r - Center and radius^2 of new sphere
80 * C1, r1, C2, r2 - Centers and radii^2 of spheres to merge
81 *
82 * OUTPUT
83 *
84 * C, r
85 *
86 * RETURNS
87 *
88 * AUTHOR
89 *
90 * Dieter Bayer
91 *
92 * DESCRIPTION
93 *
94 * Calculate a sphere that encloses two given spheres.
95 *
96 * CHANGES
97 *
98 * Jul 1994 : Creation.
99 *
100 * Oct 1994 : Added test for enclosed spheres. Calculate optimal sphere. [DB]
101 *
102 ******************************************************************************/
103
merge_spheres(VECTOR C,DBL * r,VECTOR C1,DBL r1,VECTOR C2,DBL r2)104 static void merge_spheres(VECTOR C, DBL *r, VECTOR C1, DBL r1, VECTOR C2, DBL r2)
105 {
106 DBL l, r1r, r2r, k1, k2;
107 VECTOR D;
108
109 VSub(D, C1, C2);
110
111 VLength(l, D);
112
113 /* Check if one sphere encloses the other. */
114
115 r1r = sqrt(r1);
116 r2r = sqrt(r2);
117
118 if (l + r1r <= r2r)
119 {
120 Assign_Vector(C, C2);
121
122 *r = r2;
123
124 return;
125 }
126
127 if (l + r2r <= r1r)
128 {
129 Assign_Vector(C, C1);
130
131 *r = r1;
132
133 return;
134 }
135
136 k1 = (1.0 + (r1r - r2r) / l) / 2.0;
137 k2 = (1.0 + (r2r - r1r) / l) / 2.0;
138
139 VLinComb2(C, k1, C1, k2, C2);
140
141 *r = Sqr((l + r1r + r2r) / 2.0);
142 }
143
144
145
146 /*****************************************************************************
147 *
148 * FUNCTION
149 *
150 * recompute_bound
151 *
152 * INPUT
153 *
154 * Node - Pointer to node
155 *
156 * OUTPUT
157 *
158 * Node
159 *
160 * RETURNS
161 *
162 * AUTHOR
163 *
164 * Dieter Bayer
165 *
166 * DESCRIPTION
167 *
168 * Recompute the bounding sphere of a given node in the bounding hierarchy,
169 * i. e. find the bounding sphere that encloses the bounding spheres
170 * of all nodes.
171 *
172 * NOTE: The sphere found is probably not the tightest sphere possible!
173 *
174 * CHANGES
175 *
176 * Jul 1994 : Creation.
177 *
178 * Oct 1994 : Improved bounding sphere calculation. [DB]
179 *
180 ******************************************************************************/
181
recompute_bound(BSPHERE_TREE * Node)182 static void recompute_bound(BSPHERE_TREE *Node)
183 {
184 short i;
185 DBL r2;
186 VECTOR C;
187
188 Assign_Vector(C, Node->Node[0]->C);
189
190 r2 = Node->Node[0]->r2;
191
192 for (i = 1; i < Node->Entries; i++)
193 {
194 merge_spheres(C, &r2, C, r2, Node->Node[i]->C, Node->Node[i]->r2);
195 }
196
197 Assign_Vector(Node->C, C);
198
199 Node->r2 = r2;
200 }
201
202
203
204 /*****************************************************************************
205 *
206 * FUNCTION
207 *
208 * comp_elements
209 *
210 * INPUT
211 *
212 * in_a, in_b - elements to compare
213 *
214 * OUTPUT
215 *
216 * RETURNS
217 *
218 * int - result of comparison
219 *
220 * AUTHOR
221 *
222 * Dieter Bayer
223 *
224 * DESCRIPTION
225 *
226 * -
227 *
228 * CHANGES
229 *
230 * Oct 1994 : Creation. (Derived from the bounding slab creation code)
231 *
232 ******************************************************************************/
233
comp_elements(const void * in_a,const void * in_b)234 static int CDECL comp_elements(const void *in_a, const void *in_b)
235 {
236 DBL am, bm;
237
238 am = (*(BSPHERE_TREE **)in_a)->C[Axis];
239 bm = (*(BSPHERE_TREE **)in_b)->C[Axis];
240
241 if (am < bm - EPSILON)
242 {
243 return (-1);
244 }
245 else
246 {
247 if (am > bm + EPSILON)
248 {
249 return (1);
250 }
251 else
252 {
253 return (0);
254 }
255 }
256 }
257
258
259
260 /*****************************************************************************
261 *
262 * FUNCTION
263 *
264 * find_axis
265 *
266 * INPUT
267 *
268 * OUTPUT
269 *
270 * RETURNS
271 *
272 * AUTHOR
273 *
274 * Dieter Bayer
275 *
276 * DESCRIPTION
277 *
278 * -
279 *
280 * CHANGES
281 *
282 * Oct 1994 : Creation. (Derived from the bounding slab creation code)
283 *
284 ******************************************************************************/
285
find_axis(BSPHERE_TREE ** Elements,int first,int last)286 static int find_axis(BSPHERE_TREE **Elements, int first, int last)
287 {
288 int which = X;
289 int i;
290 DBL e, d = - BOUND_HUGE;
291 VECTOR C, mins, maxs;
292
293 Make_Vector(mins, BOUND_HUGE, BOUND_HUGE, BOUND_HUGE);
294 Make_Vector(maxs, -BOUND_HUGE, -BOUND_HUGE, -BOUND_HUGE);
295
296 for (i = first; i < last; i++)
297 {
298 Assign_Vector(C, Elements[i]->C);
299
300 mins[X] = min(mins[X], C[X]);
301 maxs[X] = max(maxs[X], C[X]);
302
303 mins[Y] = min(mins[Y], C[Y]);
304 maxs[Y] = max(maxs[Y], C[Y]);
305
306 mins[Z] = min(mins[Z], C[Z]);
307 maxs[Z] = max(maxs[Z], C[Z]);
308 }
309
310 e = maxs[X] - mins[X];
311
312 if (e > d)
313 {
314 d = e; which = X;
315 }
316
317 e = maxs[Y] - mins[Y];
318
319 if (e > d)
320 {
321 d = e; which = Y;
322 }
323
324 e = maxs[Z] - mins[Z];
325
326 if (e > d)
327 {
328 which = Z;
329 }
330
331 return (which);
332 }
333
334
335
336 /*****************************************************************************
337 *
338 * FUNCTION
339 *
340 * build_area_table
341 *
342 * INPUT
343 *
344 * OUTPUT
345 *
346 * RETURNS
347 *
348 * AUTHOR
349 *
350 * Dieter Bayer
351 *
352 * DESCRIPTION
353 *
354 * Generates a table of bounding sphere surface areas.
355 *
356 * CHANGES
357 *
358 * Oct 1994 : Creation. (Derived from the bounding slab creation code)
359 *
360 ******************************************************************************/
361
build_area_table(BSPHERE_TREE ** Elements,int a,int b,DBL * areas)362 static void build_area_table(BSPHERE_TREE **Elements, int a, int b, DBL *areas)
363 {
364 int i, imin, dir;
365 DBL r2;
366 VECTOR C;
367
368 if (a < b)
369 {
370 imin = a; dir = 1;
371 }
372 else
373 {
374 imin = b; dir = -1;
375 }
376
377 Assign_Vector(C, Elements[a]->C);
378
379 r2 = Elements[a]->r2;
380
381 for (i = a; i != (b + dir); i += dir)
382 {
383 merge_spheres(C, &r2, C, r2, Elements[i]->C, Elements[i]->r2);
384
385 areas[i-imin] = r2;
386 }
387 }
388
389
390
391 /*****************************************************************************
392 *
393 * FUNCTION
394 *
395 * sort_and_split
396 *
397 * INPUT
398 *
399 * OUTPUT
400 *
401 * RETURNS
402 *
403 * AUTHOR
404 *
405 * Dieter Bayer
406 *
407 * DESCRIPTION
408 *
409 * -
410 *
411 * CHANGES
412 *
413 * Oct 1994 : Creation. (Derived from the bounding slab creation code)
414 *
415 ******************************************************************************/
416
sort_and_split(BSPHERE_TREE ** Root,BSPHERE_TREE *** Elements,int * nElem,int first,int last)417 static int sort_and_split(BSPHERE_TREE **Root, BSPHERE_TREE ***Elements, int *nElem, int first, int last)
418 {
419 int size, i, best_loc;
420 DBL *area_left, *area_right;
421 DBL best_index, new_index;
422 BSPHERE_TREE *cd;
423
424 Axis = find_axis(*Elements, first, last);
425
426 size = last - first;
427
428 if (size <= 0)
429 {
430 return (1);
431 }
432
433 Do_Cooperate(1);
434
435 /*
436 * Actually, we could do this faster in several ways. We could use a
437 * logn algorithm to find the median along the given axis, and then a
438 * linear algorithm to partition along the axis. Oh well.
439 */
440
441 QSORT((void *)(*Elements + first), size, sizeof(BSPHERE_TREE *), comp_elements);
442
443 /*
444 * area_left[] and area_right[] hold the surface areas of the bounding
445 * boxes to the left and right of any given point. E.g. area_left[i] holds
446 * the surface area of the bounding box containing Elements 0 through i and
447 * area_right[i] holds the surface area of the box containing Elements
448 * i through size-1.
449 */
450
451 area_left = (DBL *)POV_MALLOC(size * sizeof(DBL), "blob bounding hierarchy");
452 area_right = (DBL *)POV_MALLOC(size * sizeof(DBL), "blob bounding hierarchy");
453
454 /* Precalculate the areas for speed. */
455
456 build_area_table(*Elements, first, last - 1, area_left);
457 build_area_table(*Elements, last - 1, first, area_right);
458
459 best_index = area_right[0] * (size - 3.0);
460
461 best_loc = - 1;
462
463 /*
464 * Find the most effective point to split. The best location will be
465 * the one that minimizes the function N1*A1 + N2*A2 where N1 and N2
466 * are the number of objects in the two groups and A1 and A2 are the
467 * surface areas of the bounding boxes of the two groups.
468 */
469
470 for (i = 0; i < size - 1; i++)
471 {
472 new_index = (i + 1) * area_left[i] + (size - 1 - i) * area_right[i + 1];
473
474 if (new_index < best_index)
475 {
476 best_index = new_index;
477 best_loc = i + first;
478 }
479 }
480
481 POV_FREE(area_left);
482 POV_FREE(area_right);
483
484 /*
485 * Stop splitting if the BRANCHING_FACTOR is reached or
486 * if splitting stops being effective.
487 */
488
489 if ((size <= BRANCHING_FACTOR) || (best_loc < 0))
490 {
491 cd = (BSPHERE_TREE *)POV_MALLOC(sizeof(BSPHERE_TREE), "blob bounding hierarchy");
492
493 cd->Entries = (short)size;
494
495 cd->Node = (BSPHERE_TREE **)POV_MALLOC(size*sizeof(BSPHERE_TREE *), "blob bounding hierarchy");
496
497 for (i = 0; i < size; i++)
498 {
499 cd->Node[i] = (*Elements)[first+i];
500 }
501
502 recompute_bound(cd);
503
504 *Root = cd;
505
506 if (*nElem >= maxelements)
507 {
508 /* Prim array overrun, increase array by 50%. */
509
510 maxelements = 1.5 * maxelements;
511
512 /* For debugging only. */
513
514 Debug_Info("Reallocing elements to %d\n", maxelements);
515
516 *Elements = (BSPHERE_TREE **)POV_REALLOC(*Elements, maxelements * sizeof(BSPHERE_TREE *), "bounding slabs");
517 }
518
519 (*Elements)[*nElem] = cd;
520
521 (*nElem)++;
522
523 return (1);
524 }
525 else
526 {
527 sort_and_split(Root, Elements, nElem, first, best_loc + 1);
528
529 sort_and_split(Root, Elements, nElem, best_loc + 1, last);
530
531 return (0);
532 }
533 }
534
535
536
537 /*****************************************************************************
538 *
539 * FUNCTION
540 *
541 * Build_Bounding_Sphere_Hierarchy
542 *
543 * INPUT
544 *
545 * nElem - number of elements in Elements
546 * Elements - array containing nElem elements
547 *
548 * OUTPUT
549 *
550 * Root - root node of the hierarchy
551 *
552 * RETURNS
553 *
554 * AUTHOR
555 *
556 * Dieter Bayer
557 *
558 * DESCRIPTION
559 *
560 * Create the bounding sphere hierarchy for a given set of elements.
561 * One element consists of an element number (index into the array
562 * of elements; e.g. blob components, triangles) and a bounding
563 * sphere enclosing this element (center and squared radius).
564 *
565 * CHANGES
566 *
567 * Oct 1994 : Creation. (Derived from the bounding slab creation code)
568 *
569 ******************************************************************************/
570
Build_Bounding_Sphere_Hierarchy(BSPHERE_TREE ** Root,int nElem,BSPHERE_TREE *** Elements)571 void Build_Bounding_Sphere_Hierarchy(BSPHERE_TREE **Root, int nElem, BSPHERE_TREE ***Elements)
572 {
573 int low, high;
574
575 /*
576 * This is a resonable guess at the number of elements needed.
577 * This array will be reallocated as needed if it isn't.
578 */
579
580 maxelements = 2 * nElem;
581
582 /*
583 * Do a sort on the elements, with the end result being
584 * a tree of objects sorted along the x, y, and z axes.
585 */
586
587 if (nElem > 0)
588 {
589 low = 0;
590 high = nElem;
591
592 while (sort_and_split(Root, Elements, &nElem, low, high) == 0)
593 {
594 low = high;
595 high = nElem;
596
597 Do_Cooperate(0);
598 }
599 }
600 }
601
602
603
604 /*****************************************************************************
605 *
606 * FUNCTION
607 *
608 * Destroy_Bounding_Sphere_Hierarchy
609 *
610 * INPUT
611 *
612 * Node - Pointer to current node
613 *
614 * OUTPUT
615 *
616 * Node
617 *
618 * RETURNS
619 *
620 * AUTHOR
621 *
622 * Dieter Bayer
623 *
624 * DESCRIPTION
625 *
626 * Destroy bounding sphere hierarchy.
627 *
628 * CHANGES
629 *
630 * Aug 1994 : Creation.
631 *
632 * Dec 1994 : Fixed memory leakage. [DB]
633 *
634 ******************************************************************************/
635
Destroy_Bounding_Sphere_Hierarchy(BSPHERE_TREE * Node)636 void Destroy_Bounding_Sphere_Hierarchy(BSPHERE_TREE *Node)
637 {
638 short i;
639
640 if (Node != NULL)
641 {
642 if (Node->Entries > 0)
643 {
644 /* This is a node. Free sub-nodes. */
645
646 for (i = 0; i < Node->Entries; i++)
647 {
648 Destroy_Bounding_Sphere_Hierarchy(Node->Node[i]);
649 }
650
651 POV_FREE(Node->Node);
652 }
653
654 POV_FREE(Node);
655 }
656 }
657
658 END_POV_NAMESPACE
659