home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Mega Top 1
/
os2_top1.zip
/
os2_top1
/
APPS
/
RAYTRACE
/
POV22F
/
SOURCE
/
BLOB.C
< prev
next >
Wrap
C/C++ Source or Header
|
1993-08-10
|
19KB
|
658 lines
/****************************************************************************
* blob.c
*
* This module contains the code for the blob shape.
*
* This file was written by Alexander Enzmann. He wrote the code for
* blobs and generously provided us these enhancements.
*
* from Persistence of Vision Raytracer
* Copyright 1993 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 leaving a message in CompuServe's Graphics Developer's
* Forum. The latest version of POV-Ray may be found there as well.
*
* 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"
#ifndef min
#define min(x,y) ((x)<(y)?(x):(y))
#endif
#ifndef max
#define max(x,y) ((x)>(y)?(x):(y))
#endif
METHODS Blob_Methods =
{
All_Blob_Intersections,
Inside_Blob, Blob_Normal,
Copy_Blob,
Translate_Blob, Rotate_Blob, Scale_Blob, Transform_Blob,
Invert_Blob, Destroy_Blob
};
extern long Ray_Blob_Tests, Ray_Blob_Tests_Succeeded;
extern int Shadow_Test_Flag;
static int determine_influences PARAMS((VECTOR *P, VECTOR *D, BLOB *blob, DBL mindist));
static DBL calculate_field_value PARAMS ((OBJECT *obj, VECTOR *Pos));
#ifndef Blob_Tolerance
#define Blob_Tolerance 1.0e-3
#endif
#define COEFF_LIMIT 1.0e-20
#define INSIDE_TOLERANCE 1.0e-6
/* Starting with the density function: (1-r^2)^2, we have a field
that varies in strength from 1 at r = 0 to 0 at r = 1. By
substituting r/rad for r, we can adjust the range of influence
of a particular component. By multiplication by coeff, we can
adjust the amount of total contribution, giving the formula:
coeff * (1 - (r/rad)^2)^2
This varies in strength from coeff at r = 0, to 0 at r = rad. */
void
MakeBlob(blob, threshold, bloblist, npoints, sflag)
BLOB *blob;
DBL threshold;
blobstackptr bloblist;
int npoints;
int sflag;
{
unsigned i;
DBL rad, coeff;
blobstackptr temp;
VECTOR mins, maxs;
if (npoints < 1)
Error("Need at least one component in a blob.");
blob->threshold = threshold;
blob->list = (Blob_Element **)malloc(npoints*sizeof(Blob_Element *));
if (blob->list == NULL)
MAError("blob data");
for (i=0;i < (unsigned)npoints;i++)
{
blob->list[i] = (Blob_Element *)malloc(sizeof(Blob_Element));
if (blob->list[i] == NULL)
MAError("blob data");
}
blob->count = npoints;
blob->Sturm_Flag = sflag;
/* Initialize the blob data */
for(i=0;i < (unsigned)npoints;i++)
{
temp = bloblist;
if (fabs(temp->elem.coeffs[2]) < EPSILON ||
temp->elem.radius2 < EPSILON)
{
perror("Degenerate blob element\n");
}
/* Store blob specific information */
rad = temp->elem.radius2;
rad *= rad;
coeff = temp->elem.coeffs[2];
blob->list[i]->radius2 = rad;
blob->list[i]->coeffs[2] = coeff;
blob->list[i]->coeffs[1] = -(2.0 * coeff) / rad;
blob->list[i]->coeffs[0] = coeff / (rad * rad);
blob->list[i]->pos.x = temp->elem.pos.x;
blob->list[i]->pos.y = temp->elem.pos.y;
blob->list[i]->pos.z = temp->elem.pos.z;
rad = temp->elem.radius2;
if (i == 0)
{
/* First component, just set the bounds */
Make_Vector(&mins,
temp->elem.pos.x-rad,
temp->elem.pos.y-rad,
temp->elem.pos.z-rad)
Make_Vector(&maxs,
temp->elem.pos.x+rad,
temp->elem.pos.y+rad,
temp->elem.pos.z+rad)
}
else
{
/* Check min/max on the bounds */
mins.x = min(mins.x, temp->elem.pos.x - rad);
mins.y = min(mins.y, temp->elem.pos.y - rad);
mins.z = min(mins.z, temp->elem.pos.z - rad);
maxs.x = max(maxs.x, temp->elem.pos.x + rad);
maxs.y = max(maxs.y, temp->elem.pos.y + rad);
maxs.z = max(maxs.z, temp->elem.pos.z + rad);
}
bloblist = bloblist->next;
free(temp);
}
blob->Bounds.Lower_Left = mins;
VSub(blob->Bounds.Lengths, maxs, mins);
/* Allocate memory for intersection intervals */
npoints *= 2;
blob->intervals = (Blob_Interval *)malloc(npoints*sizeof(Blob_Interval));
if (blob->intervals == NULL)
MAError("blob data");
}
/* Make a sorted list of points along the ray that the various blob
components start and stop adding their influence. It would take
a very complex blob (with many components along the current ray)
to warrant the overhead of using a faster sort technique. */
static int
determine_influences(P, D, blob, mindist)
VECTOR *P, *D;
BLOB *blob;
DBL mindist;
{
int i, j, k, cnt;
DBL b, t, t0, t1, disc;
VECTOR V;
Blob_Interval *intervals = blob->intervals;
cnt = 0;
for (i=0;i<blob->count;i++)
{
/* Use standard sphere intersection routine
to determine where the ray hits the volume
of influence of each component of the blob. */
VSub(V, blob->list[i]->pos, *P);
VDot(b, V, *D);
VDot(t, V, V);
disc = b * b - t + blob->list[i]->radius2;
if (disc < EPSILON)
continue;
disc = sqrt(disc);
t1 = b + disc;
if (t1 < mindist) t1 = 0.0;
t0 = b - disc;
if (t0 < mindist) t0 = 0.0;
if (t1 == t0) continue;
else if (t1 < t0)
{
disc = t0;
t0 = t1;
t1 = disc;
}
/* Store the points of intersection of this
blob with the ray. Keep track of: whether
this is the start or end point of the hit,
which component was pierced by the ray,
and the point along the ray that the
hit occured at. */
for (k=0;k<cnt && t0 > intervals[k].bound;k++);
if (k<cnt)
{
/* This hit point is smaller than one that
already exists - bump the rest and insert
it here */
for (j=cnt;j>k;j--)
memcpy(&intervals[j], &intervals[j-1],
sizeof(Blob_Interval));
intervals[k].type = 0;
intervals[k].index = i;
intervals[k].bound = t0;
cnt++;
for (k=k+1;k<cnt && t1 > intervals[k].bound;k++);
if (k<cnt)
{
for (j=cnt;j>k;j--)
memcpy(&intervals[j], &intervals[j-1],
sizeof(Blob_Interval));
intervals[k].type = 1;
intervals[k].index = i;
intervals[k].bound = t1;
}
else
{
intervals[cnt].type = 1;
intervals[cnt].index = i;
intervals[cnt].bound = t1;
}
cnt++;
}
else
{
/* Just plop the start and end points at
the end of the list */
intervals[cnt].type = 0;
intervals[cnt].index = i;
intervals[cnt].bound = t0;
cnt++;
intervals[cnt].type = 1;
intervals[cnt].index = i;
intervals[cnt].bound = t1;
cnt++;
}
}
return cnt;
}
/* Calculate the field value of a blob - the position vector
"Pos" must already have been transformed into blob space. */
static DBL
calculate_field_value(obj, Pos)
OBJECT *obj;
VECTOR *Pos;
{
int i;
DBL len, density;
VECTOR V;
Blob_Element *ptr;
BLOB *blob = (BLOB *)obj;
density = 0.0;
for (i=0;i<blob->count;i++)
{
ptr = blob->list[i];
VSub(V, ptr->pos, *Pos);
VDot(len, V, V);
if (len < ptr->radius2)
{
/* Inside the radius of influence of this
component, add it's contribution */
density += len * (len * ptr->coeffs[0] +
ptr->coeffs[1]) +
ptr->coeffs[2];
}
}
return density;
}
/* Generate intervals of influence of each component. After these
are made, determine their aggregate effect on the ray. As the
individual intervals are checked, a quartic is generated
that represents the density at a particular point on the ray.
After making the substitutions in MakeBlob, there is a formula
for each component that has the form:
c0 * r^4 + c1 * r^2 + c2.
In order to determine the influence on the ray of all of the
individual components, we start by determining the distance
from any point on the ray to the specified point. This can
be found using the pythagorean theorem, using C as the center
of this component, P as the start of the ray, and D as the
direction of travel of the ray:
r^2 = (t * D + P - C) . (t * D + P - C)
we insert this equation for each appearance of r^2 in the
components' formula, giving:
r^2 = D.D t^2 + 2 t D . (P - C) + (P - C) . (P - C)
Since the direction vector has been normalized, D.D = 1.
Using the substitutions:
t0 = (P - C) . (P - C),
t1 = D . (P - C)
We can write the formula as:
r^2 = t0 + 2 t t1 + t^2
Taking r^2 and substituting into the formula for this component
of the blob we get the formula:
density = c0 * (r^2)^2 + c1 * r^2 + c2,
or:
density = c0 * (t0 + 2 t t1 + t^2)^2 +
c1 * (t0 + 2 t t1 + t^2) +
c2
Expanding terms and collecting with respect to "t" gives:
t^4 * c0 +
t^3 * 4 c0 t1 +
t^2 * (c1 + 2 * c0 t0 + 4 c0 t1^2)
t * 2 (c1 t1 + 2 c0 t0 t1) +
c2 + c1*t0 + c0*t0^2
This formula can now be solved for "t" by any of the quartic
root solvers that are available.
*/
int All_Blob_Intersections(Object, Ray, Depth_Stack)
OBJECT *Object;
RAY *Ray;
ISTACK *Depth_Stack;
{
BLOB *blob = (BLOB *)Object;
DBL dist, len, *tcoeffs, coeffs[5], roots[4];
int i, j, cnt;
VECTOR P, D, V;
int root_count, in_flag;
Blob_Element *element;
DBL t0, t1, c0, c1, c2;
VECTOR IPoint, dv;
Blob_Interval *intervals = blob->intervals;
int Intersection_Found = FALSE;
Ray_Blob_Tests++;
/* Transform the ray into the blob space */
if (blob->Trans != NULL)
{
MInvTransPoint(&P, &Ray->Initial, blob->Trans);
MInvTransDirection(&D, &Ray->Direction, blob->Trans);
}
else
{
P = Ray->Initial;
D = Ray->Direction;
}
len = sqrt(D.x * D.x + D.y * D.y + D.z * D.z);
if (len == 0.0)
return 0;
else
{
D.x /= len;
D.y /= len;
D.z /= len;
}
/* Figure out the intervals along the ray where each
component of the blob has an effect. */
if ((cnt = determine_influences(&P, &D, blob, 0.01)) == 0)
/* Ray doesn't hit the sphere of influence of any of
its component elements */
return 0;
/* Clear out the coefficients */
for (i=0;i<4;i++) coeffs[i] = 0.0;
coeffs[4] = -blob->threshold;
/* Step through the list of influence points, adding the
influence of each blob component as it appears */
for (i=0,in_flag=0;i<cnt;i++)
{
if (intervals[i].type == 0)
{
/* Something is just starting to influence the ray,
so calculate its coefficients and add them
into the pot. */
in_flag++;
element = blob->list[intervals[i].index];
VSub(V, P, element->pos);
c0 = element->coeffs[0];
c1 = element->coeffs[1];
c2 = element->coeffs[2];
VDot(t0, V, V);
VDot(t1, V, D);
tcoeffs = &(element->tcoeffs[0]);
tcoeffs[0] = c0;
tcoeffs[1] = 4.0 * c0 * t1;
tcoeffs[2] = 2.0 * c0 * (2.0 * t1 * t1 + t0) + c1;
tcoeffs[3] = 2.0 * t1 * (2.0 * c0 * t0 + c1);
tcoeffs[4] = c0 * t0 * t0 + c1 * t0 + c2;
for (j=0;j<5;j++) coeffs[j] += tcoeffs[j];
}
else
{
/* We are losing the influence of a component, so
subtract off its coefficients */
tcoeffs = &(blob->list[intervals[i].index]->tcoeffs[0]);
for (j=0;j<5;j++) coeffs[j] -= tcoeffs[j];
if (--in_flag == 0)
/* None of the components are currently affecting
the ray - skip ahead. */
continue;
}
/* Figure out which root solver to use */
if (blob->Sturm_Flag == 0)
/* Use Ferrari's method */
root_count = solve_quartic(coeffs, &roots[0]);
else
/* Sturm sequences */
if (fabs(coeffs[0]) < COEFF_LIMIT)
if (fabs(coeffs[1]) < COEFF_LIMIT)
root_count = solve_quadratic(&coeffs[2], &roots[0]);
else
root_count = polysolve(3, &coeffs[1], &roots[0]);
else
root_count = polysolve(4, coeffs, &roots[0]);
/* See if any of the roots are valid */
for(j=0;j<root_count;j++)
{
dist = roots[j];
/* First see if the root is in the interval of influence of
the currently active components of the blob */
if ((dist >= intervals[i].bound) &&
(dist <= intervals[i+1].bound) &&
(dist > Blob_Tolerance))
{
VScale(IPoint, D, dist);
VAdd(IPoint, IPoint, P);
/* Transform the point into world space */
if (blob->Trans != NULL)
MTransPoint(&IPoint, &IPoint, blob->Trans);
VSub(dv, IPoint, Ray->Initial);
VLength(len, dv);
if (Point_In_Clip(&IPoint, Object->Clip))
{
push_entry(len,IPoint,Object,Depth_Stack);
Intersection_Found = TRUE;
}
}
}
}
if(Intersection_Found)
Ray_Blob_Tests_Succeeded++;
return Intersection_Found;
}
/* Calculate the density at this point, then compare to
the threshold to see if we are in or out of the blob */
int Inside_Blob (Test_Point, Object)
VECTOR *Test_Point;
OBJECT *Object;
{
VECTOR New_Point;
BLOB *blob = (BLOB *) Object;
/* Transform the point into blob space */
if (blob->Trans != NULL)
MInvTransPoint(&New_Point, Test_Point, blob->Trans);
else
New_Point = *Test_Point;
if (calculate_field_value(Object, &New_Point) >
blob->threshold-INSIDE_TOLERANCE)
return ((int) 1-blob->Inverted);
else
return ((int) blob->Inverted);
}
void Blob_Normal (Result, Object, IPoint)
OBJECT *Object;
VECTOR *Result, *IPoint;
{
VECTOR New_Point, V;
int i;
DBL dist, val;
BLOB *blob = (BLOB *) Object;
Blob_Element *temp;
/* Transform the point into the blobs space */
if (blob->Trans != NULL)
MInvTransPoint(&New_Point, IPoint, blob->Trans);
else
New_Point = *IPoint;
Make_Vector(Result, 0, 0, 0);
/* For each component that contributes to this point, add
its bit to the normal */
for(i=0;i<blob->count;i++)
{
temp = blob->list[i];
V.x = New_Point.x - temp->pos.x;
V.y = New_Point.y - temp->pos.y;
V.z = New_Point.z - temp->pos.z;
dist = (V.x * V.x + V.y * V.y + V.z * V.z);
if (dist <= temp->radius2)
{
val = -2.0 * (2.0 * temp->coeffs[0] * dist +
temp->coeffs[1]);
Result->x += val * V.x;
Result->y += val * V.y;
Result->z += val * V.z;
}
}
val = (Result->x * Result->x + Result->y * Result->y +
Result->z * Result->z);
if (val < EPSILON)
{
Result->x = 1.0;
Result->y = 0.0;
Result->z = 0.0;
}
else
{
val = 1.0 / sqrt(val);
Result->x *= val;
Result->y *= val;
Result->z *= val;
}
/* Transform back to world space */
if (blob->Trans != NULL)
MTransNormal(Result, Result, blob->Trans);
VNormalize(*Result, *Result);
}
void Translate_Blob (Object, Vector)
OBJECT *Object;
VECTOR *Vector;
{
TRANSFORM Trans;
Compute_Translation_Transform(&Trans, Vector);
Transform_Blob(Object, &Trans);
}
void Rotate_Blob (Object, Vector)
OBJECT *Object;
VECTOR *Vector;
{
TRANSFORM Trans;
Compute_Rotation_Transform(&Trans, Vector);
Transform_Blob(Object, &Trans);
}
void Scale_Blob (Object, Vector)
OBJECT *Object;
VECTOR *Vector;
{
TRANSFORM Trans;
Compute_Scaling_Transform(&Trans, Vector);
Transform_Blob(Object, &Trans);
}
void Transform_Blob(Object,Trans)
OBJECT *Object;
TRANSFORM *Trans;
{
if (((BLOB *)Object)->Trans == NULL)
((BLOB *)Object)->Trans = Create_Transform();
recompute_bbox(&Object->Bounds, Trans);
Compose_Transforms(((BLOB *)Object)->Trans, Trans);
}
void Invert_Blob(Object)
OBJECT *Object;
{
((BLOB *) Object)->Inverted = 1 - ((BLOB *)Object)->Inverted;
}
void *Copy_Blob (Object)
OBJECT *Object;
{
BLOB *New;
unsigned i, cnt;
New = Create_Blob();
*New = * ((BLOB *)Object);
New->Trans = Copy_Transform(New->Trans);
cnt = ((BLOB *)Object)->count;
New->list = (Blob_Element **)malloc(cnt * sizeof(Blob_Element *));
if (New->list == NULL)
MAError("blob data");
for (i=0;i<cnt;i++)
{
New->list[i] = (Blob_Element *)malloc(sizeof(Blob_Element));
if (New->list[i] == NULL)
MAError("blob data");
memcpy(New->list[i], ((BLOB *)Object)->list[i], sizeof(Blob_Element));
}
New->intervals = (Blob_Interval *)malloc(2*New->count*sizeof(Blob_Interval));
if (New->intervals == NULL)
MAError("blob data");
return (New);
}
/* Allocate a blob. */
BLOB *Create_Blob()
{
BLOB *New;
if ((New = (BLOB *) malloc (sizeof (BLOB))) == NULL)
MAError ("blob");
INIT_OBJECT_FIELDS(New, BLOB_OBJECT, &Blob_Methods)
New->Trans = NULL;
New->Inverted = FALSE;
New->count = 0;
New->threshold = 0.0;
New->list = NULL;
New->intervals = NULL;
New->Sturm_Flag = FALSE;
return (New);
}
void Destroy_Blob (Object)
OBJECT *Object;
{
unsigned i;
Destroy_Transform(((BLOB *)Object)->Trans);
for (i=0;i < (unsigned)((BLOB *)Object)->count;i++)
free (((BLOB *)Object)->list[i]);
free (((BLOB *)Object)->list);
free (((BLOB *)Object)->intervals);
free (Object);
}