home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
OS/2 Shareware BBS: Graphics
/
Graphics.zip
/
povsrc31.zip
/
bsphere.c
< prev
next >
Wrap
C/C++ Source or Header
|
1999-10-20
|
14KB
|
652 lines
/*****************************************************************************
* bsphere.c
*
* This module implements the bounding sphere calculations.
*
* from Persistence of Vision(tm) Ray Tracer
* Copyright 1996,1999 Persistence of Vision Team
*---------------------------------------------------------------------------
* NOTICE: This source code file is provided so that users may experiment
* with enhancements to POV-Ray and to port the software to platforms other
* than those supported by the POV-Ray Team. There are strict rules under
* which you are permitted to use this file. The rules are in the file
* named POVLEGAL.DOC which should be distributed with this file.
* If POVLEGAL.DOC is not available or for more info please contact the POV-Ray
* Team Coordinator by email to team-coord@povray.org or visit us on the web at
* http://www.povray.org. The latest version of POV-Ray may be found at this site.
*
* This program is based on the popular DKB raytracer version 2.12.
* DKBTrace was originally written by David K. Buck.
* DKBTrace Ver 2.0-2.12 were written by David K. Buck & Aaron A. Collins.
*
*****************************************************************************/
#include "frame.h"
#include "vector.h"
#include "povproto.h"
#include "bsphere.h"
/*****************************************************************************
* Local preprocessor defines
******************************************************************************/
#define BRANCHING_FACTOR 4
/*****************************************************************************
* Local typedefs
******************************************************************************/
/*****************************************************************************
* Local variables
******************************************************************************/
static int maxelements, Axis;
/*****************************************************************************
* Static functions
******************************************************************************/
static void merge_spheres (VECTOR C, DBL *r, VECTOR C1, DBL r1, VECTOR C2, DBL r2);
static void recompute_bound (BSPHERE_TREE *Node);
static int find_axis (BSPHERE_TREE **Elements, int first, int last);
static void build_area_table (BSPHERE_TREE **Elements, int a, int b, DBL *areas);
static int sort_and_split (BSPHERE_TREE **Root, BSPHERE_TREE **Elements, int *nElem, int first, int last);
static int CDECL comp_elements (CONST void *in_a, CONST void *in_b);
/*****************************************************************************
*
* FUNCTION
*
* merge_spheres
*
* INPUT
*
* C, r - Center and radius^2 of new sphere
* C1, r1, C2, r2 - Centers and radii^2 of spheres to merge
*
* OUTPUT
*
* C, r
*
* RETURNS
*
* AUTHOR
*
* Dieter Bayer
*
* DESCRIPTION
*
* Calculate a sphere that encloses two given spheres.
*
* CHANGES
*
* Jul 1994 : Creation.
*
* Oct 1994 : Added test for enclosed spheres. Calculate optimal sphere. [DB]
*
******************************************************************************/
static void merge_spheres(VECTOR C, DBL *r, VECTOR C1, DBL r1, VECTOR C2, DBL r2)
{
DBL l, r1r, r2r, k1, k2;
VECTOR D;
VSub(D, C1, C2);
VLength(l, D);
/* Check if one sphere encloses the other. */
r1r = sqrt(r1);
r2r = sqrt(r2);
if (l + r1r <= r2r)
{
Assign_Vector(C, C2);
*r = r2;
return;
}
if (l + r2r <= r1r)
{
Assign_Vector(C, C1);
*r = r1;
return;
}
k1 = (1.0 + (r1r - r2r) / l) / 2.0;
k2 = (1.0 + (r2r - r1r) / l) / 2.0;
VLinComb2(C, k1, C1, k2, C2);
*r = Sqr((l + r1r + r2r) / 2.0);
}
/*****************************************************************************
*
* FUNCTION
*
* recompute_bound
*
* INPUT
*
* Node - Pointer to node
*
* OUTPUT
*
* Node
*
* RETURNS
*
* AUTHOR
*
* Dieter Bayer
*
* DESCRIPTION
*
* Recompute the bounding sphere of a given node in the bounding hierarchy,
* i. e. find the bounding sphere that encloses the bounding spheres
* of all nodes.
*
* NOTE: The sphere found is probably not the tightest sphere possible!
*
* CHANGES
*
* Jul 1994 : Creation.
*
* Oct 1994 : Improved bounding sphere calculation. [DB]
*
******************************************************************************/
static void recompute_bound(BSPHERE_TREE *Node)
{
short i;
DBL r2;
VECTOR C;
COOPERATE_1
Assign_Vector(C, Node->Node[0]->C);
r2 = Node->Node[0]->r2;
for (i = 1; i < Node->Entries; i++)
{
merge_spheres(C, &r2, C, r2, Node->Node[i]->C, Node->Node[i]->r2);
}
Assign_Vector(Node->C, C);
Node->r2 = r2;
}
/*****************************************************************************
*
* FUNCTION
*
* comp_elements
*
* INPUT
*
* in_a, in_b - elements to compare
*
* OUTPUT
*
* RETURNS
*
* int - result of comparison
*
* AUTHOR
*
* Dieter Bayer
*
* DESCRIPTION
*
* -
*
* CHANGES
*
* Oct 1994 : Creation. (Derived from the bounding slab creation code)
*
******************************************************************************/
static int CDECL comp_elements(CONST void *in_a, CONST void *in_b)
{
DBL am, bm;
am = (*(BSPHERE_TREE **)in_a)->C[Axis];
bm = (*(BSPHERE_TREE **)in_b)->C[Axis];
if (am < bm - EPSILON)
{
return (-1);
}
else
{
if (am > bm + EPSILON)
{
return (1);
}
else
{
return (0);
}
}
}
/*****************************************************************************
*
* FUNCTION
*
* find_axis
*
* INPUT
*
* OUTPUT
*
* RETURNS
*
* AUTHOR
*
* Dieter Bayer
*
* DESCRIPTION
*
* -
*
* CHANGES
*
* Oct 1994 : Creation. (Derived from the bounding slab creation code)
*
******************************************************************************/
static int find_axis(BSPHERE_TREE **Elements, int first, int last)
{
int which = X;
int i;
DBL e, d = - BOUND_HUGE;
VECTOR C, mins, maxs;
Make_Vector(mins, BOUND_HUGE, BOUND_HUGE, BOUND_HUGE);
Make_Vector(maxs, -BOUND_HUGE, -BOUND_HUGE, -BOUND_HUGE);
for (i = first; i < last; i++)
{
Assign_Vector(C, Elements[i]->C);
mins[X] = min(mins[X], C[X]);
maxs[X] = max(maxs[X], C[X]);
mins[Y] = min(mins[Y], C[Y]);
maxs[Y] = max(maxs[Y], C[Y]);
mins[Z] = min(mins[Z], C[Z]);
maxs[Z] = max(maxs[Z], C[Z]);
}
e = maxs[X] - mins[X];
if (e > d)
{
d = e; which = X;
}
e = maxs[Y] - mins[Y];
if (e > d)
{
d = e; which = Y;
}
e = maxs[Z] - mins[Z];
if (e > d)
{
which = Z;
}
return (which);
}
/*****************************************************************************
*
* FUNCTION
*
* build_area_table
*
* INPUT
*
* OUTPUT
*
* RETURNS
*
* AUTHOR
*
* Dieter Bayer
*
* DESCRIPTION
*
* Generates a table of bounding sphere surface areas.
*
* CHANGES
*
* Oct 1994 : Creation. (Derived from the bounding slab creation code)
*
******************************************************************************/
static void build_area_table(BSPHERE_TREE **Elements, int a, int b, DBL *areas)
{
int i, imin, dir;
DBL r2;
VECTOR C;
if (a < b)
{
imin = a; dir = 1;
}
else
{
imin = b; dir = -1;
}
Assign_Vector(C, Elements[a]->C);
r2 = Elements[a]->r2;
for (i = a; i != (b + dir); i += dir)
{
merge_spheres(C, &r2, C, r2, Elements[i]->C, Elements[i]->r2);
areas[i-imin] = r2;
}
}
/*****************************************************************************
*
* FUNCTION
*
* sort_and_split
*
* INPUT
*
* OUTPUT
*
* RETURNS
*
* AUTHOR
*
* Dieter Bayer
*
* DESCRIPTION
*
* -
*
* CHANGES
*
* Oct 1994 : Creation. (Derived from the bounding slab creation code)
*
******************************************************************************/
static int sort_and_split(BSPHERE_TREE **Root, BSPHERE_TREE **Elements, int *nElem, int first, int last)
{
int size, i, best_loc;
DBL *area_left, *area_right;
DBL best_index, new_index;
BSPHERE_TREE *cd;
Axis = find_axis(Elements, first, last);
size = last - first;
if (size <= 0)
{
return (1);
}
/*
* Actually, we could do this faster in several ways. We could use a
* logn algorithm to find the median along the given axis, and then a
* linear algorithm to partition along the axis. Oh well.
*/
QSORT((void *)(&Elements[first]), (size_t)size, sizeof(BSPHERE_TREE *), comp_elements);
/*
* area_left[] and area_right[] hold the surface areas of the bounding
* boxes to the left and right of any given point. E.g. area_left[i] holds
* the surface area of the bounding box containing Elements 0 through i and
* area_right[i] holds the surface area of the box containing Elements
* i through size-1.
*/
area_left = (DBL *)POV_MALLOC(size * sizeof(DBL), "blob bounding hierarchy");
area_right = (DBL *)POV_MALLOC(size * sizeof(DBL), "blob bounding hierarchy");
/* Precalculate the areas for speed. */
build_area_table(Elements, first, last - 1, area_left);
build_area_table(Elements, last - 1, first, area_right);
best_index = area_right[0] * (size - 3.0);
best_loc = - 1;
/*
* Find the most effective point to split. The best location will be
* the one that minimizes the function N1*A1 + N2*A2 where N1 and N2
* are the number of objects in the two groups and A1 and A2 are the
* surface areas of the bounding boxes of the two groups.
*/
for (i = 0; i < size - 1; i++)
{
new_index = (i + 1) * area_left[i] + (size - 1 - i) * area_right[i + 1];
if (new_index < best_index)
{
best_index = new_index;
best_loc = i + first;
}
}
POV_FREE(area_left);
POV_FREE(area_right);
/*
* Stop splitting if the BRANCHING_FACTOR is reached or
* if splitting stops being effective.
*/
if ((size <= BRANCHING_FACTOR) || (best_loc < 0))
{
cd = (BSPHERE_TREE *)POV_MALLOC(sizeof(BSPHERE_TREE), "blob bounding hierarchy");
cd->Entries = (short)size;
cd->Node = (BSPHERE_TREE **)POV_MALLOC(size*sizeof(BSPHERE_TREE *), "blob bounding hierarchy");
for (i = 0; i < size; i++)
{
cd->Node[i] = Elements[first+i];
}
recompute_bound(cd);
*Root = cd;
if (*nElem > maxelements)
{
/* Prim array overrun, increase array by 50%. */
maxelements = 1.5 * maxelements;
/* For debugging only. */
Debug_Info("Reallocing elements to %d\n", maxelements);
Elements = (BSPHERE_TREE **)POV_REALLOC(Elements, maxelements * sizeof(BSPHERE_TREE *), "bounding slabs");
}
Elements[*nElem] = cd;
(*nElem)++;
return (1);
}
else
{
sort_and_split(Root, Elements, nElem, first, best_loc + 1);
sort_and_split(Root, Elements, nElem, best_loc + 1, last);
return (0);
}
}
/*****************************************************************************
*
* FUNCTION
*
* Build_Bounding_Sphere_Hierarchy
*
* INPUT
*
* nElem - number of elements in Elements
* Elements - array containing nElem elements
*
* OUTPUT
*
* Root - root node of the hierarchy
*
* RETURNS
*
* AUTHOR
*
* Dieter Bayer
*
* DESCRIPTION
*
* Create the bounding sphere hierarchy for a given set of elements.
* One element consists of an element number (index into the array
* of elements; e.g. blob components, triangles) and a bounding
* sphere enclosing this element (center and squared radius).
*
* CHANGES
*
* Oct 1994 : Creation. (Derived from the bounding slab creation code)
*
******************************************************************************/
void Build_Bounding_Sphere_Hierarchy(BSPHERE_TREE **Root, int nElem, BSPHERE_TREE **Elements)
{
int low, high;
/*
* This is a resonable guess at the number of elements needed.
* This array will be reallocated as needed if it isn't.
*/
maxelements = 2 * nElem;
/*
* Do a sort on the elements, with the end result being
* a tree of objects sorted along the x, y, and z axes.
*/
if (nElem > 0)
{
low = 0;
high = nElem;
while (sort_and_split(Root, Elements, &nElem, low, high) == 0)
{
low = high;
high = nElem;
}
}
}
/*****************************************************************************
*
* FUNCTION
*
* Destroy_Bounding_Sphere_Hierarchy
*
* INPUT
*
* Node - Pointer to current node
*
* OUTPUT
*
* Node
*
* RETURNS
*
* AUTHOR
*
* Dieter Bayer
*
* DESCRIPTION
*
* Destroy bounding sphere hierarchy.
*
* CHANGES
*
* Aug 1994 : Creation.
*
* Dec 1994 : Fixed memory leakage. [DB]
*
******************************************************************************/
void Destroy_Bounding_Sphere_Hierarchy(BSPHERE_TREE *Node)
{
short i;
if (Node != NULL)
{
if (Node->Entries > 0)
{
/* This is a node. Free sub-nodes. */
for (i = 0; i < Node->Entries; i++)
{
Destroy_Bounding_Sphere_Hierarchy(Node->Node[i]);
}
POV_FREE(Node->Node);
}
POV_FREE(Node);
}
}