[LONG] Re: boundary-spheres (boxes) and 3d-transformation
From zbaker@venom.st.hmc.edu (Zach Baker)
Organization The New Venom
Date 19 May 1997 13:25:51 GMT
Summary bounding box wu graphics gems finding minimal optimal approximation
Newsgroups comp.sys.ibm.pc.demos,comp.graphics.algorithms
Message-ID <5lpkcv$8i8$1@cinenews.claremont.edu>
References 1
(crossposted to comp.graphics.algorithms from original article in
comp.sys.ibm.pc.demos)
In article <337AFCC4.2BB3@hamburg.netsurf.de>,
Nils Pipenbrinck <nils.pipenbrinck@hamburg.netsurf.de> wrote:
>[...]
>Since I do all clipping in 3d I want to find out if my boundary-sphere
>would intersect with the view-fustrum. The problem is, that I have very
>complicated transformation matrices (eg. multiple translations and
>rotations). The middlepoint of the object may be easy transformed into
>viewspace using the same code as I use for my polygon-vertices, but the
>radius of the sphere must be handled different (I think).
Well, what you can do is calculate the frustum (as a pyramid of 5
points, chop off the frustum at some Z) in the object space you want,
but you'll have to take the inverse of the matrix. Then just use the
3 points on each side of the frustum to define a plane, and find the
signed distance between the center of the bounding-sphere and the four
planes, then compare and you're done. Not brilliant but pretty simple.
>I have an idea to get around this problem: Using a boundary-box (8
>instead of 1 vertice).
>[...]
>Which way do you go?
I like just using an (oriented) boundary box...
>Oh, by the way, if anyone has code to calculate the minimal boundary box
>of an object please let me know!
...because I did something for this last week.
> I found an algorithm in Graphic Gems III, but I'm to stupid to
> calculate the eigenvectors of a matrix.
Well, it's not like it's all that easy, so don't feel bad. There's
not many situations for eigenanalysis in this type of stuff anyway...
> (found a lot of matrix-libraries on the net, but most of them even
> refuse to compile under watcom... )
Oh yes, I got the feeling they were all cribbed from Recipes and were
converted to C from Fortran. But fortunately I found a good symmetric
3x3 matrix eigenvector-finding routine from Dave Eberly's MAGIC
library (Hi, Dave! Thanks!) that is thankfully simple.
So anyway, I'm posting what I got out of Wu's Graphics Gems article,
but I have some problems with the article's algorithm:
1) It's not implemented! Come on, this is one of the ideas behind
Graphics Gems in the first place, right?! So, I can't be entirely
sure that this is what the author had in mind.
2) It seems like it should use the points that make up the 3D convex
hull of the point set. The approximation is bad if there are many
points inside the convex hull which define a different principle axis
than the convex hull of the points. If anyone has a good algorithm
for computing the points of a 3D convex hull I'd like to see it.
But it's otherwise very useful and gives good results. If you want
optimal results, however, it's not a substitute for the solution in
Preparata & Shamos, it's just a good fast approximation, as it says.
After my signature follow two C files in plain text, be sure to cut
them into 2 separate files. I use BCPL-style // comments (some call
them C++-style, yeah, whatever) so be sure to enable them; gcc uses
"-Wp,-lang-c-c++-comments" to do this.
The first file is all you need if you want to make a module that is
called with an array of points and a place to return the box's corner
point and three size vectors. It uses several math.h routines.
The second file is a cheesy PC test program. I ripped out some of the
goodies from 3dsrdr (Thanks, Jare!) to read vertex lists from 3ds
files, made a routine that prints out an asc file of a bounding box
for the file's vertices, and just used a global array of 100000
vertices to read points into.
BTW, I could have gotten this all wrong and be a trained chimp who's
learned to post on Usenet, so don't assume this will work. It's a
good policy in any case, so test it.
---
Zach Baker <zbaker@venom.st.hmc.edu>
"Some non-essential information regarding the Pentium processor..."
---
// -----------------------CUT HERE: makeobb.c-----------------------
// --- makeobb.c
/* Contains code adapted to C from MAGIC by Dave Eberly.
http://www2.cs.unc.edu/~eberly/
This is a MODIFIED version of a file from the original MAGIC
library, which is subject to the following license:
*/
/* Permission is granted to use, copy, and distribute at no
charge the MAGIC library code. Permission is also granted to
the user to modify the code. However, if the modified code is
distributed, then all such modifications must be documented
to indicate that they are not part of the original package.
The original code comes with absolutely no warranty and no
guarantee is made that the code is bug-free.
*/
/* All other code is subject to the following license:
Permission is granted to freely use, copy, modify and distribute.
THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR
IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
*/
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
/* from MAGIC */
void EigenSymmetric3(double mat[3][3], double eval[3],
double evec[3][3]);
#define X 0
#define Y 1
#define Z 2
void Zero(double pt[3]){
pt[X] = 0.0;
pt[Y] = 0.0;
pt[Z] = 0.0;
}
void MakeCovarOBB(const float pts[][3], int npts,
double boxpt[3], double boxvecs[3][3]){
double oon = 1.0/npts;
double mCoVars[3][3], mXForm[3][3];
double ptAvg[3], ptAvgSq[3], ptAvgXS[3];
double EigenVals[3];
double XfMin[3], XfMax[3];
double XfTemp;
int i;
/* check for pointless case */
if(!npts){
Zero(boxpt);
Zero(boxvecs[0]);
Zero(boxvecs[1]);
Zero(boxvecs[2]);
return;
}
/* Make averages for covariance matrix. */
/* Intitialize averages. */
Zero(ptAvg);
Zero(ptAvgSq);
Zero(ptAvgXS);
/* Sum all point data. */
for(i=0; i < npts; i++){
ptAvg[X] += pts[i][X];
ptAvg[Y] += pts[i][Y];
ptAvg[Z] += pts[i][Z];
ptAvgSq[X] += pts[i][X]*pts[i][X];
ptAvgSq[Y] += pts[i][Y]*pts[i][Y];
ptAvgSq[Z] += pts[i][Z]*pts[i][Z];
ptAvgXS[X] += pts[i][X]*pts[i][Y];
ptAvgXS[Y] += pts[i][X]*pts[i][Z];
ptAvgXS[Z] += pts[i][Y]*pts[i][Z];
}
/* Divide sums to make average. */
ptAvg[X] *= oon;
ptAvg[Y] *= oon;
ptAvg[Z] *= oon;
ptAvgSq[X] *= oon;
ptAvgSq[Y] *= oon;
ptAvgSq[Z] *= oon;
ptAvgXS[X] *= oon;
ptAvgXS[Y] *= oon;
ptAvgXS[Z] *= oon;
/* Compose covariance matrix. */
mCoVars[0][0] = ptAvgSq[X] - ptAvg[X]*ptAvg[X];
mCoVars[0][1] = ptAvgXS[X] - ptAvg[X]*ptAvg[Y];
mCoVars[0][2] = ptAvgXS[Y] - ptAvg[X]*ptAvg[Z];
mCoVars[1][0] = ptAvgXS[X] - ptAvg[Y]*ptAvg[X];
mCoVars[1][1] = ptAvgSq[Y] - ptAvg[Y]*ptAvg[Y];
mCoVars[1][2] = ptAvgXS[Z] - ptAvg[Y]*ptAvg[Z];
mCoVars[2][0] = ptAvgXS[Y] - ptAvg[Z]*ptAvg[X];
mCoVars[2][1] = ptAvgXS[Z] - ptAvg[Z]*ptAvg[Y];
mCoVars[2][2] = ptAvgSq[Z] - ptAvg[Z]*ptAvg[Z];
/* Find eigenvectors for the covariance matrix. */
/* put eigenvectors into mXForm */
EigenSymmetric3(mCoVars, EigenVals, mXForm);
/* Normalize the eigenvectors to make unit axes. */
for(i=0; i < 3; i++){
XfTemp = (float)(1.0/sqrt(mXForm[i][0]*mXForm[i][0] +
mXForm[i][1]*mXForm[i][1] +
mXForm[i][2]*mXForm[i][2]));
mXForm[i][0] *= XfTemp;
mXForm[i][1] *= XfTemp;
mXForm[i][2] *= XfTemp;
}
/* find minimum/maximum of points in new coord
system defined by the eigenvectors */
/* initialize min & max with first point's values */
XfMin[X] = pts[0][X]*mXForm[0][0] +
pts[0][Y]*mXForm[0][1] +
pts[0][Z]*mXForm[0][2];
XfMax[X] = XfMin[X];
XfMin[Y] = pts[0][X]*mXForm[1][0] +
pts[0][Y]*mXForm[1][1] +
pts[0][Z]*mXForm[1][2];
XfMax[Y] = XfMin[Y];
XfMin[Z] = pts[0][X]*mXForm[2][0] +
pts[0][Y]*mXForm[2][1] +
pts[0][Z]*mXForm[2][2];
XfMax[Z] = XfMin[Z];
/* check subsequent points */
for(i=1; i < npts; i++){
XfTemp = pts[i][X]*mXForm[0][0] +
pts[i][Y]*mXForm[0][1] +
pts[i][Z]*mXForm[0][2];
if(XfMin[X] > XfTemp)
XfMin[X] = XfTemp;
if(XfMax[X] < XfTemp)
XfMax[X] = XfTemp;
XfTemp = pts[i][X]*mXForm[1][0] +
pts[i][Y]*mXForm[1][1] +
pts[i][Z]*mXForm[1][2];
if(XfMin[Y] > XfTemp)
XfMin[Y] = XfTemp;
if(XfMax[Y] < XfTemp)
XfMax[Y] = XfTemp;
XfTemp = pts[i][X]*mXForm[2][0] +
pts[i][Y]*mXForm[2][1] +
pts[i][Z]*mXForm[2][2];
if(XfMin[Z] > XfTemp)
XfMin[Z] = XfTemp;
if(XfMax[Z] < XfTemp)
XfMax[Z] = XfTemp;
}
/* Now the orthogonal bounding box in transformed
coordinates is in XfMin and XfMax. */
/* Make averages for covariance matrix. */
/* Intitialize averages. */
Zero(ptAvg);
Zero(ptAvgSq);
/* Sum all point data. */
for(i=0; i < npts; i++){
ptAvg[X] += pts[i][X];
ptAvg[Y] += pts[i][Y];
ptAvg[Z] += pts[i][Z];
ptAvgSq[X] += pts[i][X]*pts[i][X];
ptAvgSq[Y] += pts[i][Y]*pts[i][Y];
ptAvgSq[Z] += pts[i][X]*pts[i][Y];
}
/* Divide sums to make average. */
ptAvg[X] *= oon;
ptAvg[Y] *= oon;
ptAvg[Z] *= oon;
ptAvgSq[X] *= oon;
ptAvgSq[Y] *= oon;
ptAvgSq[Z] *= oon;
mCoVars[0][0] = ptAvgSq[X] - ptAvg[X]*ptAvg[X];
mCoVars[0][1] = ptAvgSq[Z] - ptAvg[X]*ptAvg[Y];
mCoVars[1][0] = ptAvgSq[Z] - ptAvg[Y]*ptAvg[X];
mCoVars[1][1] = ptAvgSq[Y] - ptAvg[Y]*ptAvg[Y];
/* Transform XfMin back into original coordinates
to find the box's base. */
/* Use special case for orthogonal unit vector matrix:
Since det == 1, inverse is M^T. */
/* Find coordinates of base point. */
boxpt[X] = XfMin[X]*mXForm[0][0] +
XfMin[Y]*mXForm[1][0] +
XfMin[Z]*mXForm[2][0];
boxpt[Y] = XfMin[X]*mXForm[0][1] +
XfMin[Y]*mXForm[1][1] +
XfMin[Z]*mXForm[2][1];
boxpt[Z] = XfMin[X]*mXForm[0][2] +
XfMin[Y]*mXForm[1][2] +
XfMin[Z]*mXForm[2][2];
/* Use individual differences from XfMax to find
the width, length, and depth vectors, transform
them and place in boxvecs. */
XfTemp = XfMax[X] - XfMin[X];
boxvecs[0][X] = XfTemp*mXForm[0][0];
boxvecs[0][Y] = XfTemp*mXForm[0][1];
boxvecs[0][Z] = XfTemp*mXForm[0][2];
XfTemp = XfMax[Y] - XfMin[Y];
boxvecs[1][X] = XfTemp*mXForm[1][0];
boxvecs[1][Y] = XfTemp*mXForm[1][1];
boxvecs[1][Z] = XfTemp*mXForm[1][2];
XfTemp = XfMax[Z] - XfMin[Z];
boxvecs[2][X] = XfTemp*mXForm[2][0];
boxvecs[2][Y] = XfTemp*mXForm[2][1];
boxvecs[2][Z] = XfTemp*mXForm[2][2];
/* done */
}
/* beginning of MAGIC code */
typedef double REAL;
const REAL epsilon = 1e-06f;
// roots r0 < r1 < r2
#define distinctRoots 0
// root r0
#define singleRoot 1
// roots r0 = r1 < r2
#define doubleRoot01 2
// roots r0 < r1 = r2
#define doubleRoot12 4
// roots r0 = r1 = r2
#define tripleRoot 6
//===========================================================================
int CubicRoots (REAL c0, REAL c1, REAL c2, REAL *r0, REAL *r1, REAL *r2)
{
// polynomial is L^3-c2*L^2+c1*L-c0
int maxiter = 50;
REAL discr, temp, pval, pdval, b0, b1;
int i;
// find local extrema (if any) of p'(L) = 3*L^2-2*c2*L+c1
discr = c2*c2-3*c1;
if ( discr >= 0.0 ) {
discr = (REAL)sqrt(discr);
temp = (c2+discr)/3;
pval = temp*(temp*(temp-c2)+c1)-c0;
if ( pval >= 0.0 ) {
// real root occurs before the positive local maximum
(*r0) = (c2-discr)/3 - 1; // initial guess for Newton's methods
pval = 2*epsilon;
for (i = 0; i < maxiter && fabs(pval) > epsilon; i++) {
pval = (*r0)*((*r0)*((*r0)-c2)+c1)-c0;
pdval = (*r0)*(3*(*r0)-2*c2)+c1;
(*r0) -= pval/pdval;
}
// Other two roots are solutions to quadratic equation
// L^2 + ((*r0)-c2)*L + [(*r0)*((*r0)-c2)+c1] = 0.
b1 = (*r0)-c2;
b0 = (*r0)*((*r0)-c2)+c1;
discr = b1*b1-4*b0;
if ( discr < -epsilon )
{
// single root r0
return singleRoot;
}
else
{
int result = distinctRoots;
// roots r0 <= r1 <= r2
discr = (REAL)sqrt(fabs(discr));
(*r1) = (REAL)0.5*(-b1-discr);
(*r2) = (REAL)0.5*(-b1+discr);
if ( fabs((*r0)-(*r1)) <= epsilon )
{
(*r0) = (*r1);
result |= doubleRoot01;
}
if ( fabs((*r1)-(*r2)) <= epsilon )
{
(*r1) = (*r2);
result |= doubleRoot12;
}
return result;
}
}
else {
// real root occurs after the negative local minimum
(*r2) = temp + 1; // initial guess for Newton's method
pval = 2*epsilon;
for (i = 0; i < maxiter && fabs(pval) > epsilon; i++) {
pval = (*r2)*((*r2)*((*r2)-c2)+c1)-c0;
pdval = (*r2)*(3*(*r2)-2*c2)+c1;
(*r2) -= pval/pdval;
}
// Other two roots are solutions to quadratic equation
// L^2 + (r2-c2)*L + [r2*(r2-c2)+c1] = 0.
b1 = (*r2)-c2;
b0 = (*r2)*((*r2)-c2)+c1;
discr = b1*b1-4*b0;
if ( discr < -epsilon )
{
// single root
(*r0) = (*r2);
return singleRoot;
}
else
{
int result = distinctRoots;
// roots r0 <= r1 <= r2
discr = (REAL)sqrt(fabs(discr));
(*r0) = (REAL)0.5*(-b1-discr);
(*r1) = (REAL)0.5*(-b1+discr);
if ( fabs((*r0)-(*r1)) <= epsilon )
{
(*r0) = (*r1);
result |= doubleRoot01;
}
if ( fabs((*r1)-(*r2)) <= epsilon )
{
(*r1) = (*r2);
result |= doubleRoot12;
}
return result;
}
}
}
else {
// p(L) has one real root
(*r0) = c0;
pval = 2*epsilon;
for (i = 0; i < maxiter && fabs(pval) > epsilon; i++) {
pval = (*r0)*((*r0)*((*r0)-c2)+c1)-c0;
pdval = (*r0)*(3*(*r0)-2*c2)+c1;
(*r0) -= pval/pdval;
}
return singleRoot;
}
}
//---------------------------------------------------------------------------
void Tridiagonal3 (REAL mat[3][3], REAL diag[3], REAL subd[2])
{
REAL a, b, c, d, e, f, ell, q;
a = mat[0][0];
b = mat[0][1];
c = mat[0][2];
d = mat[1][1];
e = mat[1][2];
f = mat[2][2];
diag[0] = a;
subd[2] = (REAL)0;
if ( fabs(c) >= epsilon ) {
ell = (REAL)sqrt(b*b+c*c);
b /= ell;
c /= ell;
q = 2*b*e+c*(f-d);
diag[1] = d+c*q;
diag[2] = f-c*q;
subd[0] = ell;
subd[1] = e-b*q;
mat[0][0] = (REAL)1; mat[0][1] = (REAL)0; mat[0][2] = (REAL)0;
mat[1][0] = (REAL)0; mat[1][1] = b; mat[1][2] = c;
mat[2][0] = (REAL)0; mat[2][1] = c; mat[2][2] = -b;
}
else {
diag[1] = d;
diag[2] = f;
subd[0] = b;
subd[1] = e;
mat[0][0] = (REAL)1; mat[0][1] = (REAL)0; mat[0][2] = (REAL)0;
mat[1][0] = (REAL)0; mat[1][1] = (REAL)1; mat[1][2] = (REAL)0;
mat[2][0] = (REAL)0; mat[2][1] = (REAL)0; mat[2][2] = (REAL)1;
}
}
//---------------------------------------------------------------------------
void PreEigenvector (REAL diag[3], REAL subd[2], REAL eval, REAL evec[3][3],
int *numvec)
{
if ( fabs(subd[0]) >= epsilon )
{
if ( fabs(subd[1]) >= epsilon )
{
REAL diff = eval-diag[0];
evec[(*numvec)][0] = subd[0]*subd[1];
evec[(*numvec)][1] = subd[1]*diff;
evec[(*numvec)][2] = diff*(eval-diag[1])-subd[0]*subd[0];
(*numvec)++;
}
else
{
evec[(*numvec)][0] = subd[0];
evec[(*numvec)][1] = eval-diag[0];
evec[(*numvec)][2] = (REAL)0;
(*numvec)++;
if ( fabs(eval-diag[2]) < epsilon )
{
evec[(*numvec)][0] = (REAL)0;
evec[(*numvec)][1] = (REAL)0;
evec[(*numvec)][2] = (REAL)1;
(*numvec)++;
}
}
}
else
{
if ( fabs(subd[1]) >= epsilon )
{
evec[(*numvec)][0] = (REAL)0;
evec[(*numvec)][1] = subd[1];
evec[(*numvec)][2] = eval-diag[1];
(*numvec)++;
if ( fabs(eval-diag[0]) < epsilon )
{
evec[(*numvec)][0] = (REAL)1;
evec[(*numvec)][1] = (REAL)0;
evec[(*numvec)][2] = (REAL)0;
}
}
else
{
if ( fabs(eval-diag[0]) < epsilon )
{
evec[(*numvec)][0] = (REAL)1;
evec[(*numvec)][1] = (REAL)0;
evec[(*numvec)][2] = (REAL)0;
(*numvec)++;
}
if ( fabs(eval-diag[1]) < epsilon )
{
evec[(*numvec)][0] = (REAL)0;
evec[(*numvec)][1] = (REAL)1;
evec[(*numvec)][2] = (REAL)0;
(*numvec)++;
}
if ( fabs(eval-diag[2]) < epsilon )
{
evec[(*numvec)][0] = (REAL)0;
evec[(*numvec)][1] = (REAL)0;
evec[(*numvec)][2] = (REAL)1;
(*numvec)++;
}
}
}
}
//---------------------------------------------------------------------------
void EigenSymmetric3 (REAL mat[3][3], REAL eval[3], REAL evec[3][3])
{
REAL diag[3], subd[2];
REAL c0, c1, c2, sum;
REAL prevec[3][3]; // matrix P, vectors stored as rows, not columns
int k, row, col, flag, numvec;
// Householder reduction T = Q^t M Q
// Input:
// mat, symmetric 3x3 matrix M
// Output:
// mat, orthogonal matrix Q
// diag, diagonal entries of T
// subd, subdiagonal entries of T (T is symmetric)
Tridiagonal3(mat,diag,subd);
// eigenvalue construction
c0 = diag[0]*diag[1]*diag[2]-diag[0]*subd[1]*subd[1]
-diag[2]*subd[0]*subd[0];
c1 = diag[0]*diag[1]+diag[0]*diag[2]+diag[1]*diag[2]
-subd[0]*subd[0]-subd[1]*subd[1];
c2 = diag[0]+diag[1]+diag[2];
flag = CubicRoots(c0,c1,c2,eval + 0, eval + 1, eval + 2);
// eigenvector construction
numvec = 0;
switch ( flag )
{
case distinctRoots:
PreEigenvector(diag,subd,eval[0],prevec,&numvec);
PreEigenvector(diag,subd,eval[1],prevec,&numvec);
PreEigenvector(diag,subd,eval[2],prevec,&numvec);
break;
case doubleRoot01:
case doubleRoot12:
PreEigenvector(diag,subd,eval[0],prevec,&numvec);
PreEigenvector(diag,subd,eval[2],prevec,&numvec);
break;
case tripleRoot:
PreEigenvector(diag,subd,eval[0],prevec,&numvec);
break;
}
// transform back to original coordinates by E^t = PQ where the
// rows of E are the eigenvectors of M
for (row = 0; row < 3; row++)
{
for (col = 0; col < 3; col++)
{
sum = (REAL)0;
for (k = 0; k < 3; k++)
sum += prevec[row][k]*mat[k][col];
evec[row][col] = sum;
}
}
}
/* end of MAGIC code */
/* EOF */
// -----------------------CUT HERE: testobb.c-----------------------
// --- testobb.c
/* Permission is granted to freely use, copy, modify and distribute.
THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR
IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
*/
// Reads 3DStudio binary/3DS file and outputs 3DStudio ASCII/ASC
// file containing a bounding box of all vertices in the 3DS file
// to demonstrate usage of MakeCovarOBB().
// Test program adapted from public domain file 3dsrdr.c v1.3.
#include <stdlib.h>
#include <stdio.h>
// really really cheesy, I know
#define MAX_VERTS 100000
float Verts[MAX_VERTS][3];
int NumVerts = 0;
void MakeCovarOBB(const float pts[][3], int npts,
double boxpt[3], double boxvecs[3][3]);
void Read(FILE *f, long fomax);
void OutputASCBox(const double boxpts[8][3]);
void Add(double pt[3], const double ptadd[3]){
pt[0] += ptadd[0];
pt[1] += ptadd[1];
pt[2] += ptadd[2];
}
int main(int argc, char *argv[]) {
FILE *f;
long flen;
double obbvecs[3][3];
double obbpt[3];
double obbbox[8][3];
int i;
if (argc < 2) {
printf("Writes covariance bounding box of a .3DS object to "
"stdout in .ASC format.\n");
printf("Usage: %s file.3DS > file-obb.ASC\n", argv[0]);
return 1;
}
f = fopen(argv[1], "rb");
if (f == NULL) {
printf("Can't open %s!\n", argv[1]);
return 1;
}
fseek(f, 0, SEEK_END);
flen = ftell(f);
fseek(f, 0, SEEK_SET);
Read(f, flen);
// use vertices Read() to make OBB
MakeCovarOBB(Verts, NumVerts, obbpt, obbvecs);
// make 8 points that define this particular
// box for the ASC file
// no expansion of box, so there might be
// precision problems
for(i=0; i < 8; i++){
obbbox[i][0] = obbpt[0];
obbbox[i][1] = obbpt[1];
obbbox[i][2] = obbpt[2];
}
// XYZ
Add(obbbox[0], obbvecs[0]);
Add(obbbox[0], obbvecs[1]);
Add(obbbox[0], obbvecs[2]);
// YZ
Add(obbbox[1], obbvecs[1]);
Add(obbbox[1], obbvecs[2]);
// Z
Add(obbbox[2], obbvecs[2]);
// XZ
Add(obbbox[3], obbvecs[0]);
Add(obbbox[3], obbvecs[2]);
// XY
Add(obbbox[4], obbvecs[0]);
Add(obbbox[4], obbvecs[1]);
// Y
Add(obbbox[5], obbvecs[1]);
// 0
// X
Add(obbbox[7], obbvecs[0]);
// output ASC file with bounding box
OutputASCBox(obbbox);
fclose(f);
// OK
return 0;
}
void Read(FILE *f, long fomax){
unsigned short id;
unsigned int len;
long fo;
// look through blocks until end file offset reached
while(fo = ftell(f), fo < fomax){
// get block ID
if(fread(&id, sizeof(id), 1, f) != 1)
return;
// get block size
if(fread(&len, sizeof(len), 1, f) != 1 || !len)
return;
if(id == 0x4000){ // object
int c;
// skip through string identifier
while ((c = fgetc(f)) != EOF && c != '\0')
;
}
if(id == 0x4110){ // vertex list
unsigned short nv;
// get # of vertices
if(fread(&nv, sizeof(nv), 1, f) != 1)
return;
// stash in Verts
while(nv-- > 0){
if(fread(Verts[NumVerts], 3*sizeof(float), 1, f) != 1)
return;
if(NumVerts == MAX_VERTS-1)
break;
else
NumVerts++;
}
}
// descend hierarchy for appropriate composite groups
if(id == 0x3D3D || id == 0x4000 ||
id == 0x4100 || id == 0x4D4D)
Read(f, fo + len);
// go to end of block
fseek(f, fo + len, SEEK_SET);
if(feof(f) || ferror(f))
return;
}
}
void OutputASCBox(const double boxpts[8][3]){
int nv;
printf("Ambient light color: Red=0.3 Green=0.3 Blue=0.3\n\n");
printf("Named object: \"Bounding Box\"\n");
printf("Tri-mesh, Vertices: 8 Faces: 12\n");
printf("Vertex list:\n");
for(nv=0; nv < 8; nv++)
printf("Vertex %i: X:%f Y:%f Z:%f\n", nv,
boxpts[nv][0], boxpts[nv][1], boxpts[nv][2]);
printf("Face list:\n");
printf("Face 0: A:2 B:1 C:0 AB:1 BC:1 CA:1\n");
printf("Material:\"r220g150b150a100\"\n");
printf("Face 1: A:3 B:2 C:0 AB:1 BC:1 CA:1\n");
printf("Material:\"r220g150b150a100\"\n");
printf("Face 2: A:2 B:3 C:7 AB:1 BC:1 CA:1\n");
printf("Material:\"r220g150b150a100\"\n");
printf("Face 3: A:6 B:2 C:7 AB:1 BC:1 CA:1\n");
printf("Material:\"r220g150b150a100\"\n");
printf("Face 4: A:7 B:3 C:0 AB:1 BC:1 CA:1\n");
printf("Material:\"r220g150b150a100\"\n");
printf("Face 5: A:4 B:7 C:0 AB:1 BC:1 CA:1\n");
printf("Material:\"r220g150b150a100\"\n");
printf("Face 6: A:5 B:4 C:0 AB:1 BC:1 CA:1\n");
printf("Material:\"r220g150b150a100\"\n");
printf("Face 7: A:1 B:5 C:0 AB:1 BC:1 CA:1\n");
printf("Material:\"r220g150b150a100\"\n");
printf("Face 8: A:5 B:1 C:2 AB:1 BC:1 CA:1\n");
printf("Material:\"r220g150b150a100\"\n");
printf("Face 9: A:6 B:5 C:2 AB:1 BC:1 CA:1\n");
printf("Material:\"r220g150b150a100\"\n");
printf("Face 10: A:7 B:4 C:5 AB:1 BC:1 CA:1\n");
printf("Material:\"r220g150b150a100\"\n");
printf("Face 11: A:6 B:7 C:5 AB:1 BC:1 CA:1\n");
printf("Material:\"r220g150b150a100\"\n");
}
/* EOF */