g_superell.c

Go to the documentation of this file.
00001 /*                    G _ S U P E R E L L . C
00002  * BRL-CAD
00003  *
00004  * Copyright (c) 1985-2006 United States Government as represented by
00005  * the U.S. Army Research Laboratory.
00006  *
00007  * This library is free software; you can redistribute it and/or
00008  * modify it under the terms of the GNU Lesser General Public License
00009  * as published by the Free Software Foundation; either version 2 of
00010  * the License, or (at your option) any later version.
00011  *
00012  * This library is distributed in the hope that it will be useful, but
00013  * WITHOUT ANY WARRANTY; without even the implied warranty of
00014  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00015  * Library General Public License for more details.
00016  *
00017  * You should have received a copy of the GNU Lesser General Public
00018  * License along with this file; see the file named COPYING for more
00019  * information.
00020  */
00021 
00022 /** @addtogroup g_  */
00023 
00024 /*@{*/
00025 /** @file g_superell.c
00026  * Intersect a ray with a Superquadratic Ellipsoid.
00027  *
00028  *      NOTICE: this primitive is incomplete and should beconsidered
00029  *              experimental.  this primitive will exhibit several
00030  *              instabilities in the existing root solver method.
00031  *
00032  *
00033  *
00034  *  Authors -
00035  *      Christopher Sean Morrison (Programming)
00036  *      Edwin O. Davisson (Mathematics)
00037  *
00038  *  Source -
00039  *      SECAD/VLD Computing Consortium, Bldg 394
00040  *      The U. S. Army Ballistic Research Laboratory
00041  *      Aberdeen Proving Ground, Maryland  21005
00042  *
00043  */
00044 /*@}*/
00045 
00046 #ifndef lint
00047 static const char RCSsuperell[] = "@(#)$Header: /cvsroot/brlcad/brlcad/src/librt/g_superell.c,v 14.16 2006/09/16 02:04:25 lbutler Exp $ (BRL)";
00048 #endif
00049 
00050 #include "common.h"
00051 
00052 #include <stddef.h>
00053 #include <stdio.h>
00054 #ifdef HAVE_STRING_H
00055 #  include <string.h>
00056 #else
00057 #  include <strings.h>
00058 #endif
00059 #include <math.h>
00060 
00061 #include "machine.h"
00062 #include "vmath.h"
00063 #include "db.h"
00064 #include "nmg.h"
00065 #include "raytrace.h"
00066 #include "nurb.h"
00067 #include "rtgeom.h"
00068 #include "./debug.h"
00069 
00070 const struct bu_structparse rt_superell_parse[] = {
00071     { "%f", 3, "V", bu_offsetof(struct rt_superell_internal, v[X]), BU_STRUCTPARSE_FUNC_NULL },
00072     { "%f", 3, "A", bu_offsetof(struct rt_superell_internal, a[X]), BU_STRUCTPARSE_FUNC_NULL },
00073     { "%f", 3, "B", bu_offsetof(struct rt_superell_internal, b[X]), BU_STRUCTPARSE_FUNC_NULL },
00074     { "%f", 3, "C", bu_offsetof(struct rt_superell_internal, c[X]), BU_STRUCTPARSE_FUNC_NULL },
00075     { "%f", 1, "n", bu_offsetof(struct rt_superell_internal, n), BU_STRUCTPARSE_FUNC_NULL },
00076     { "%f", 1, "e", bu_offsetof(struct rt_superell_internal, e), BU_STRUCTPARSE_FUNC_NULL },
00077     { {'\0','\0','\0','\0'}, 0, (char *)NULL, 0, BU_STRUCTPARSE_FUNC_NULL }
00078  };
00079 
00080 
00081 /*
00082  *  Algorithm:
00083  *
00084  *  Given V, A, B, and C, there is a set of points on this superellipsoid
00085  *
00086  *  { (x,y,z) | (x,y,z) is on superellipsoid defined by V, A, B, C }
00087  *
00088  *  Through a series of Affine Transformations, this set will be
00089  *  transformed into a set of points on a unit sphere at the origin
00090  *
00091  *  { (x',y',z') | (x',y',z') is on Sphere at origin }
00092  *
00093  *  The transformation from X to X' is accomplished by:
00094  *
00095  *  X' = S(R( X - V ))
00096  *
00097  *  where R(X) =  ( A/(|A|) )
00098  *               (  B/(|B|)  ) . X
00099  *                ( C/(|C|) )
00100  *
00101  *  and S(X) =   (  1/|A|   0     0   )
00102  *              (    0    1/|B|   0    ) . X
00103  *               (   0      0   1/|C| )
00104  *
00105  *  To find the intersection of a line with the superellipsoid, consider
00106  *  the parametric line L:
00107  *
00108  *      L : { P(n) | P + t(n) . D }
00109  *
00110  *  Call W the actual point of intersection between L and the superellipsoid.
00111  *  Let W' be the point of intersection between L' and the unit sphere.
00112  *
00113  *      L' : { P'(n) | P' + t(n) . D' }
00114  *
00115  *  W = invR( invS( W' ) ) + V
00116  *
00117  *  Where W' = k D' + P'.
00118  *
00119  *  Let dp = D' dot P'
00120  *  Let dd = D' dot D'
00121  *  Let pp = P' dot P'
00122  *
00123  *  and k = [ -dp +/- sqrt( dp*dp - dd * (pp - 1) ) ] / dd
00124  *  which is constant.
00125  *
00126  *  Now, D' = S( R( D ) )
00127  *  and  P' = S( R( P - V ) )
00128  *
00129  *  Substituting,
00130  *
00131  *  W = V + invR( invS[ k *( S( R( D ) ) ) + S( R( P - V ) ) ] )
00132  *    = V + invR( ( k * R( D ) ) + R( P - V ) )
00133  *    = V + k * D + P - V
00134  *    = k * D + P
00135  *
00136  *  Note that ``k'' is constant, and is the same in the formulations
00137  *  for both W and W'.
00138  *
00139  *  NORMALS.  Given the point W on the superellipsoid, what is the vector
00140  *  normal to the tangent plane at that point?
00141  *
00142  *  Map W onto the unit sphere, ie:  W' = S( R( W - V ) ).
00143  *
00144  *  Plane on unit sphere at W' has a normal vector of the same value(!).
00145  *  N' = W'
00146  *
00147  *  The plane transforms back to the tangent plane at W, and this
00148  *  new plane (on the superellipsoid) has a normal vector of N, viz:
00149  *
00150  *  N = inverse[ transpose( inverse[ S o R ] ) ] ( N' )
00151  *
00152  *  because if H is perpendicular to plane Q, and matrix M maps from
00153  *  Q to Q', then inverse[ transpose(M) ] (H) is perpendicular to Q'.
00154  *  Here, H and Q are in "prime space" with the unit sphere.
00155  *  [Somehow, the notation here is backwards].
00156  *  So, the mapping matrix M = inverse( S o R ), because
00157  *  S o R maps from normal space to the unit sphere.
00158  *
00159  *  N = inverse[ transpose( inverse[ S o R ] ) ] ( N' )
00160  *    = inverse[ transpose(invR o invS) ] ( N' )
00161  *    = inverse[ transpose(invS) o transpose(invR) ] ( N' )
00162  *    = inverse[ inverse(S) o R ] ( N' )
00163  *    = invR o S ( N' )
00164  *
00165  *    = invR o S ( W' )
00166  *    = invR( S( S( R( W - V ) ) ) )
00167  *
00168  *  because inverse(R) = transpose(R), so R = transpose( invR ),
00169  *  and S = transpose( S ).
00170  *
00171  *  Note that the normal vector N produced above will not have unit length.
00172  */
00173 
00174 struct superell_specific {
00175   vect_t superell_V; /* Vector to center of superellipsoid */
00176   vect_t superell_Au; /* unit-length A vector */
00177   vect_t superell_Bu;
00178   vect_t superell_Cu;
00179   double superell_n; /* north-south curvature power */
00180   double superell_e; /* east-west curvature power */
00181   double superell_invmsAu; /* 1.0 / |Au|^2 */
00182   double superell_invmsBu; /* 1.0 / |Bu|^2 */
00183   double superell_invmsCu; /* 1.0 / |Cu|^2 */
00184   vect_t superell_invsq;
00185   mat_t superell_SoR; /* matrix for local cordinate system, Scale(Rotate(V))*/
00186   mat_t superell_invRSSR; /* invR(Scale(Scale(Rot(V)))) */
00187   mat_t superell_invR; /* transposed rotation matrix */
00188 };
00189 #define SUPERELL_NULL   ((struct superell_specific *)0)
00190 
00191 /**
00192  *                      R T _ S U P E R E L L _ P R E P
00193  *
00194  *  Given a pointer to a GED database record, and a transformation matrix,
00195  *  determine if this is a valid superellipsoid, and if so, precompute various
00196  *  terms of the formula.
00197  *
00198  *  Returns -
00199  *      0       SUPERELL is OK
00200  *      !0      Error in description
00201  *
00202  *  Implicit return -
00203  *      A struct superell_specific is created, and it's address is stored in
00204  *      stp->st_specific for use by rt_superell_shot().
00205  */
00206 int
00207 rt_superell_prep(struct soltab *stp, struct rt_db_internal *ip, struct rt_i *rtip)
00208 {
00209 
00210   register struct superell_specific *superell;
00211   struct rt_superell_internal   *eip;
00212   LOCAL fastf_t magsq_a, magsq_b, magsq_c;
00213   LOCAL mat_t   R, TEMP;
00214   LOCAL vect_t  Au, Bu, Cu;     /* A,B,C with unit length */
00215   LOCAL vect_t  w1, w2, P;      /* used for bounding RPP */
00216   LOCAL fastf_t f;
00217 
00218   eip = (struct rt_superell_internal *)ip->idb_ptr;
00219   RT_SUPERELL_CK_MAGIC(eip);
00220 
00221   /* Validate that |A| > 0, |B| > 0, |C| > 0 */
00222   magsq_a = MAGSQ( eip->a );
00223   magsq_b = MAGSQ( eip->b );
00224   magsq_c = MAGSQ( eip->c );
00225 
00226   if( magsq_a < rtip->rti_tol.dist || magsq_b < rtip->rti_tol.dist || magsq_c < rtip->rti_tol.dist ) {
00227     bu_log("superell(%s):  near-zero length A(%g), B(%g), or C(%g) vector\n",
00228            stp->st_name, magsq_a, magsq_b, magsq_c );
00229     return(1);          /* BAD */
00230   }
00231   if (eip->n < rtip->rti_tol.dist || eip->e < rtip->rti_tol.dist) {
00232     bu_log("superell(%s):  near-zero length <n,e> curvature (%g, %g) causes problems\n",
00233            stp->st_name, eip->n, eip->e);
00234     /* BAD */
00235   }
00236   if (eip->n > 10000.0 || eip->e > 10000.0) {
00237     bu_log("superell(%s):  very large <n,e> curvature (%g, %g) causes problems\n",
00238            stp->st_name, eip->n, eip->e);
00239     /* BAD */
00240   }
00241 
00242   /* Create unit length versions of A,B,C */
00243   f = 1.0/sqrt(magsq_a);
00244   VSCALE( Au, eip->a, f );
00245   f = 1.0/sqrt(magsq_b);
00246   VSCALE( Bu, eip->b, f );
00247   f = 1.0/sqrt(magsq_c);
00248   VSCALE( Cu, eip->c, f );
00249 
00250   /* Validate that A.B == 0, B.C == 0, A.C == 0 (check dir only) */
00251   f = VDOT( Au, Bu );
00252   if( ! NEAR_ZERO(f, rtip->rti_tol.dist) )  {
00253     bu_log("superell(%s):  A not perpendicular to B, f=%f\n",stp->st_name, f);
00254     return(1);          /* BAD */
00255   }
00256   f = VDOT( Bu, Cu );
00257   if( ! NEAR_ZERO(f, rtip->rti_tol.dist) )  {
00258     bu_log("superell(%s):  B not perpendicular to C, f=%f\n",stp->st_name, f);
00259     return(1);          /* BAD */
00260   }
00261   f = VDOT( Au, Cu );
00262   if( ! NEAR_ZERO(f, rtip->rti_tol.dist) )  {
00263     bu_log("superell(%s):  A not perpendicular to C, f=%f\n",stp->st_name, f);
00264     return(1);          /* BAD */
00265   }
00266 
00267   /* Solid is OK, compute constant terms now */
00268 
00269   BU_GETSTRUCT( superell, superell_specific );
00270   stp->st_specific = (genptr_t)superell;
00271 
00272   superell->superell_n = eip->n;
00273   superell->superell_e = eip->e;
00274 
00275   VMOVE( superell->superell_V, eip->v );
00276 
00277   VSET( superell->superell_invsq, 1.0/magsq_a, 1.0/magsq_b, 1.0/magsq_c );
00278   VMOVE( superell->superell_Au, Au );
00279   VMOVE( superell->superell_Bu, Bu );
00280   VMOVE( superell->superell_Cu, Cu );
00281 
00282   /* compute the inverse magnitude square for equations during shot */
00283   superell->superell_invmsAu = 1.0 / magsq_a;
00284   superell->superell_invmsBu = 1.0 / magsq_b;
00285   superell->superell_invmsCu = 1.0 / magsq_c;
00286 
00287   /* compute the rotation matrix */
00288   MAT_IDN(R);
00289   VMOVE( &R[0], Au );
00290   VMOVE( &R[4], Bu );
00291   VMOVE( &R[8], Cu );
00292   bn_mat_trn( superell->superell_invR, R );
00293 
00294   /* computer invRSSR */
00295   MAT_IDN(superell->superell_invRSSR);
00296   MAT_IDN(TEMP);
00297   TEMP[0] = superell->superell_invsq[0];
00298   TEMP[5] = superell->superell_invsq[1];
00299   TEMP[10] = superell->superell_invsq[2];
00300   bn_mat_mul(TEMP, TEMP, R);
00301   bn_mat_mul(superell->superell_invRSSR, superell->superell_invR, TEMP);
00302 
00303   /* compute Scale(Rotate(vect)) */
00304   MAT_IDN(superell->superell_SoR);
00305   VSCALE( &superell->superell_SoR[0], eip->a, superell->superell_invsq[0]);
00306   VSCALE( &superell->superell_SoR[4], eip->b, superell->superell_invsq[1]);
00307   VSCALE( &superell->superell_SoR[8], eip->c, superell->superell_invsq[2]);
00308 
00309   /* Compute bounding sphere */
00310   VMOVE( stp->st_center, eip->v );
00311   f = magsq_a;
00312   if( magsq_b > f )
00313     f = magsq_b;
00314   if( magsq_c > f )
00315     f = magsq_c;
00316   stp->st_aradius = stp->st_bradius = sqrt(f);
00317 
00318   /* Compute bounding RPP */
00319   VSET( w1, magsq_a, magsq_b, magsq_c );
00320 
00321   /* X */
00322   VSET( P, 1.0, 0, 0 );         /* bounding plane normal */
00323   MAT3X3VEC( w2, R, P );                /* map plane to local coord syst */
00324   VELMUL( w2, w2, w2 );         /* square each term */
00325   f = VDOT( w1, w2 );
00326   f = sqrt(f);
00327   stp->st_min[X] = superell->superell_V[X] - f; /* V.P +/- f */
00328   stp->st_max[X] = superell->superell_V[X] + f;
00329 
00330   /* Y */
00331   VSET( P, 0, 1.0, 0 );         /* bounding plane normal */
00332   MAT3X3VEC( w2, R, P );                /* map plane to local coord syst */
00333   VELMUL( w2, w2, w2 );         /* square each term */
00334   f = VDOT( w1, w2 );
00335   f = sqrt(f);
00336   stp->st_min[Y] = superell->superell_V[Y] - f; /* V.P +/- f */
00337   stp->st_max[Y] = superell->superell_V[Y] + f;
00338 
00339   /* Z */
00340   VSET( P, 0, 0, 1.0 );         /* bounding plane normal */
00341   MAT3X3VEC( w2, R, P );                /* map plane to local coord syst */
00342   VELMUL( w2, w2, w2 );         /* square each term */
00343   f = VDOT( w1, w2 );
00344   f = sqrt(f);
00345   stp->st_min[Z] = superell->superell_V[Z] - f; /* V.P +/- f */
00346   stp->st_max[Z] = superell->superell_V[Z] + f;
00347 
00348   return(0);                    /* OK */
00349 }
00350 
00351 /**
00352  *                      R T _ S U P E R E L L _ P R I N T
00353  */
00354 void
00355 rt_superell_print(register const struct soltab *stp)
00356 {
00357         register struct superell_specific *superell =
00358                 (struct superell_specific *)stp->st_specific;
00359 
00360         VPRINT("V", superell->superell_V);
00361 }
00362 
00363 /* Equation of a superellipsoid:
00364  *
00365  * f(x) = [ (x^(2/e2) + y^(2/e2))^(e2/e1) + z^(2/e1) ]^(e1/2) - 1
00366  */
00367 
00368 /*
00369  *                      R T _ S U P E R E L L _ S H O T
00370  *
00371  *  Intersect a ray with an superellipsoid, where all constant terms have
00372  *  been precomputed by rt_superell_prep().  If an intersection occurs,
00373  *  a struct seg will be acquired and filled in.
00374  *
00375  *  Returns -
00376  *      0       MISS
00377  *      >0      HIT
00378  */
00379 int
00380 rt_superell_shot(struct soltab *stp, register struct xray *rp, struct application *ap, struct seg *seghead)
00381 {
00382   static int counter=10;
00383 
00384   register struct superell_specific *superell = (struct superell_specific *)stp->st_specific;
00385   LOCAL bn_poly_t equation; /* equation of superell to be solved */
00386   LOCAL vect_t translated;  /* translated shot vector */
00387   LOCAL vect_t newShotPoint; /* P' */
00388   LOCAL vect_t newShotDir; /* D' */
00389   LOCAL vect_t normalizedShotPoint; /* P' with normalized dist from superell */
00390   LOCAL bn_complex_t complexRoot[4]; /* roots returned from poly solver */
00391   LOCAL double realRoot[4];  /* real ray distance values */
00392   register int i,j;
00393   register struct seg *segp;
00394 
00395   /* translate ray point */
00396   /*  VSUB2(translated, rp->r_pt, superell->superell_V); */
00397   (translated)[X] = (rp->r_pt)[X] - (superell->superell_V)[X];
00398   (translated)[Y] = (rp->r_pt)[Y] - (superell->superell_V)[Y];
00399   (translated)[Z] = (rp->r_pt)[Z] - (superell->superell_V)[Z];
00400 
00401   /* scale and rotate point to get P' */
00402 
00403   /*  MAT4X3VEC(newShotPoint, superell->superell_SoR, translated); */
00404   newShotPoint[X] = (superell->superell_SoR[0]*translated[X] + superell->superell_SoR[1]*translated[Y] + superell->superell_SoR[ 2]*translated[Z]) * 1.0/(superell->superell_SoR[15]);
00405   newShotPoint[Y] = (superell->superell_SoR[4]*translated[X] + superell->superell_SoR[5]*translated[Y] + superell->superell_SoR[ 6]*translated[Z]) * 1.0/(superell->superell_SoR[15]);
00406   newShotPoint[Z] = (superell->superell_SoR[8]*translated[X] + superell->superell_SoR[9]*translated[Y] + superell->superell_SoR[10]*translated[Z]) * 1.0/(superell->superell_SoR[15]);
00407 
00408   /* translate ray direction vector */
00409   MAT4X3VEC(newShotDir, superell->superell_SoR, rp->r_dir);
00410   VUNITIZE(newShotDir);
00411 
00412   /* normalize distance from the superell.  substitues a corrected ray
00413    * point, which contains a translation along the ray direction to the
00414    * closest approach to vertex of the superell.  Translating the ray
00415    * along the direction of the ray to the closest point near the
00416    * primitives center vertex.  New ray origin is hence, normalized.
00417    */
00418   VSCALE( normalizedShotPoint, newShotDir,
00419           VDOT( newShotPoint, newShotDir ));
00420   VSUB2( normalizedShotPoint, newShotPoint, normalizedShotPoint );
00421 
00422   /* Now generate the polynomial equation for passing to the root finder */
00423 
00424   equation.dgr = 2;
00425 
00426   /* (x^2 / A) + (y^2 / B) + (z^2 / C) - 1 */
00427   equation.cf[0] = newShotPoint[X] * newShotPoint[X] * superell->superell_invmsAu + newShotPoint[Y] * newShotPoint[Y] * superell->superell_invmsBu + newShotPoint[Z] * newShotPoint[Z] * superell->superell_invmsCu - 1;
00428   /* (2xX / A) + (2yY / B) + (2zZ / C) */
00429   equation.cf[1] = 2 * newShotDir[X] * newShotPoint[X] * superell->superell_invmsAu + 2 * newShotDir[Y] * newShotPoint[Y] * superell->superell_invmsBu + 2 * newShotDir[Z] * newShotPoint[Z] * superell->superell_invmsCu;
00430   /* (X^2 / A) + (Y^2 / B) + (Z^2 / C) */
00431   equation.cf[2] = newShotDir[X] * newShotDir[X] * superell->superell_invmsAu  + newShotDir[Y] * newShotDir[Y] * superell->superell_invmsBu + newShotDir[Z] * newShotDir[Z] * superell->superell_invmsCu;
00432 
00433   if ( (i = rt_poly_roots( &equation, complexRoot, stp->st_dp->d_namep)) != 2 ) {
00434       if (i > 0) {
00435           bu_log("superell: rt_poly_roots() 2 != %d\n", i);
00436           bn_pr_roots(stp->st_name, complexRoot, i);
00437       } else if (i < 0) {
00438           static int reported=0;
00439           bu_log("The root solver failed to converge on a solution for %s\n", stp->st_dp->d_namep);
00440           if (!reported) {
00441               VPRINT("while shooting from:\t", rp->r_pt);
00442               VPRINT("while shooting at:\t", rp->r_dir);
00443               bu_log("Additional superellipsoid convergence failure details will be suppressed.\n");
00444               reported=1;
00445           }
00446       }
00447       return (0); /* MISS */
00448   }
00449 
00450   /* XXX BEGIN CUT */
00451   /*  Only real roots indicate an intersection in real space.
00452    *
00453    *  Look at each root returned; if the imaginary part is zero
00454    *  or sufficiently close, then use the real part as one value
00455    *  of 't' for the intersections
00456    */
00457   for ( j=0, i=0; j < 2; j++ ){
00458     if( NEAR_ZERO( complexRoot[j].im, 0.001 ) )
00459       realRoot[i++] = complexRoot[j].re;
00460   }
00461 
00462   /* reverse above translation by adding distance to all 'k' values. */
00463   /*  for( j = 0; j < i; ++j )
00464       realRoot[j] -= VDOT(newShotPoint, newShotDir);
00465   */
00466 
00467   /* Here, 'i' is number of points found */
00468   switch( i )  {
00469   case 0:
00470     return(0);          /* No hit */
00471 
00472   default:
00473     bu_log("rt_superell_shot: reduced 4 to %d roots\n",i);
00474     bn_pr_roots( stp->st_name, complexRoot, 4 );
00475     return(0);          /* No hit */
00476 
00477   case 2:
00478     {
00479       /* Sort most distant to least distant. */
00480       FAST fastf_t      u;
00481       if( (u=realRoot[0]) < realRoot[1] )  {
00482         /* bubble larger towards [0] */
00483         realRoot[0] = realRoot[1];
00484         realRoot[1] = u;
00485       }
00486     }
00487     break;
00488   case 4:
00489     {
00490       register short    n;
00491       register short    lim;
00492 
00493       /*  Inline rt_pt_sort().  Sorts realRoot[] into descending order. */
00494       for( lim = i-1; lim > 0; lim-- )  {
00495         for( n = 0; n < lim; n++ )  {
00496           FAST fastf_t  u;
00497           if( (u=realRoot[n]) < realRoot[n+1] )  {
00498             /* bubble larger towards [0] */
00499             realRoot[n] = realRoot[n+1];
00500             realRoot[n+1] = u;
00501           }
00502         }
00503       }
00504     }
00505     break;
00506   }
00507 
00508   if (counter > 0) {
00509     bu_log("realroot: in %d  out %d\n", realRoot[1], realRoot[0]);
00510     counter--;
00511   }
00512 
00513 
00514   /* Now, t[0] > t[npts-1] */
00515   /* realRoot[1] is entry point, and realRoot[0] is farthest exit point */
00516   RT_GET_SEG(segp, ap->a_resource);
00517   segp->seg_stp = stp;
00518   segp->seg_in.hit_dist = realRoot[1];
00519   segp->seg_out.hit_dist = realRoot[0];
00520   /*  segp->seg_in.hit_surfno = segp->seg_out.hit_surfno = 0; */
00521   /* Set aside vector for rt_superell_norm() later */
00522   /*  VJOIN1( segp->seg_in.hit_vpriv, newShotPoint, realRoot[1], newShotDir ); */
00523   /*  VJOIN1( segp->seg_out.hit_vpriv, newShotPoint, realRoot[0], newShotDir ); */
00524   BU_LIST_INSERT( &(seghead->l), &(segp->l) );
00525 
00526   if( i == 2 ) {
00527     return(2);                  /* HIT */
00528   }
00529 
00530   /* 4 points */
00531   /* realRoot[3] is entry point, and realRoot[2] is exit point */
00532   RT_GET_SEG(segp, ap->a_resource);
00533   segp->seg_stp = stp;
00534   segp->seg_in.hit_dist = realRoot[3]*superell->superell_e;
00535   segp->seg_out.hit_dist = realRoot[2]*superell->superell_e;
00536   segp->seg_in.hit_surfno = segp->seg_out.hit_surfno = 1;
00537   VJOIN1( segp->seg_in.hit_vpriv, newShotPoint, realRoot[3], newShotDir );
00538   VJOIN1( segp->seg_out.hit_vpriv, newShotPoint, realRoot[2], newShotDir );
00539   BU_LIST_INSERT( &(seghead->l), &(segp->l) );
00540   return(4);                    /* HIT */
00541   /* XXX END CUT */
00542 
00543   /* Is there any possibility of hitting another segment?  Only when there
00544    * is a concave curvature (<n,e> > <2.0, 2.0>).
00545    */
00546   if ( (superell->superell_n > 2.0) || (superell->superell_e > 2.0) ) {
00547 
00548   }
00549 
00550   return 1;
00551 }
00552 
00553 
00554 /**
00555  *                      R T _ S U P E R E L L _ V S H O T
00556  *
00557  *  This is the Becker vector version.
00558  */
00559 void
00560 rt_superell_vshot(struct soltab **stp, struct xray **rp, struct seg *segp, int n, struct application *ap)
00561                                /* An array of solid pointers */
00562                                /* An array of ray pointers */
00563                                /* array of segs (results returned) */
00564                                /* Number of ray/object pairs */
00565 
00566 {
00567   return;
00568 }
00569 
00570 
00571 /**
00572  *                      R T _ S U P E R E L L _ N O R M
00573  *
00574  *  Given ONE ray distance, return the normal and entry/exit point.
00575  */
00576 void
00577 rt_superell_norm(register struct hit *hitp, struct soltab *stp, register struct xray *rp)
00578 {
00579   register struct superell_specific *superell =
00580     (struct superell_specific *)stp->st_specific;
00581 
00582   LOCAL vect_t xlated;
00583   LOCAL fastf_t scale;
00584 
00585   VJOIN1( hitp->hit_point, rp->r_pt, hitp->hit_dist, rp->r_dir );
00586   VSUB2( xlated, hitp->hit_point, superell->superell_V );
00587   MAT4X3VEC( hitp->hit_normal, superell->superell_invRSSR, xlated );
00588   scale = 1.0 / MAGNITUDE( hitp->hit_normal );
00589   VSCALE( hitp->hit_normal, hitp->hit_normal, scale );
00590 
00591   return;
00592 }
00593 
00594 
00595 /**
00596  *                      R T _ S U P E R E L L _ C U R V E
00597  *
00598  *  Return the curvature of the superellipsoid.
00599  */
00600 void
00601 rt_superell_curve(register struct curvature *cvp, register struct hit *hitp, struct soltab *stp)
00602 {
00603   bu_log("called rt_superell_curve!\n");
00604   return;
00605 }
00606 
00607 
00608 /**
00609  *                      R T _ S U P E R E L L _ U V
00610  *
00611  *  For a hit on the surface of an SUPERELL, return the (u,v) coordinates
00612  *  of the hit point, 0 <= u,v <= 1.
00613  *  u = azimuth
00614  *  v = elevation
00615  */
00616 void
00617 rt_superell_uv(struct application *ap, struct soltab *stp, register struct hit *hitp, register struct uvcoord *uvp)
00618 {
00619   bu_log("called rt_superell_uv!\n");
00620   return;
00621 }
00622 
00623 /**
00624  *                      R T _ S U P E R E L L _ F R E E
00625  */
00626 void
00627 rt_superell_free(register struct soltab *stp)
00628 {
00629         register struct superell_specific *superell =
00630                 (struct superell_specific *)stp->st_specific;
00631 
00632         bu_free( (char *)superell, "superell_specific" );
00633 }
00634 
00635 
00636 int
00637 rt_superell_class(void)
00638 {
00639         return(0);
00640 }
00641 
00642 
00643 /**
00644  *                      R T _ S U P E R E L L _ 1 6 P T S
00645  *
00646  * Also used by the TGC code
00647  */
00648 #define SUPERELLOUT(n)  ov+(n-1)*3
00649 void
00650 rt_superell_16pts(register fastf_t *ov,
00651              register fastf_t *V,
00652              fastf_t *A,
00653              fastf_t *B)
00654 {
00655         static fastf_t c, d, e, f,g,h;
00656 
00657         e = h = .92388;                 /* cos(22.5) */
00658         c = d = .707107;                /* cos(45) */
00659         g = f = .382683;                /* cos(67.5) */
00660 
00661         /*
00662          * For angle theta, compute surface point as
00663          *
00664          *      V  +  cos(theta) * A  + sin(theta) * B
00665          *
00666          * note that sin(theta) is cos(90-theta).
00667          */
00668 
00669         VADD2( SUPERELLOUT(1), V, A );
00670         VJOIN2( SUPERELLOUT(2), V, e, A, f, B );
00671         VJOIN2( SUPERELLOUT(3), V, c, A, d, B );
00672         VJOIN2( SUPERELLOUT(4), V, g, A, h, B );
00673         VADD2( SUPERELLOUT(5), V, B );
00674         VJOIN2( SUPERELLOUT(6), V, -g, A, h, B );
00675         VJOIN2( SUPERELLOUT(7), V, -c, A, d, B );
00676         VJOIN2( SUPERELLOUT(8), V, -e, A, f, B );
00677         VSUB2( SUPERELLOUT(9), V, A );
00678         VJOIN2( SUPERELLOUT(10), V, -e, A, -f, B );
00679         VJOIN2( SUPERELLOUT(11), V, -c, A, -d, B );
00680         VJOIN2( SUPERELLOUT(12), V, -g, A, -h, B );
00681         VSUB2( SUPERELLOUT(13), V, B );
00682         VJOIN2( SUPERELLOUT(14), V, g, A, -h, B );
00683         VJOIN2( SUPERELLOUT(15), V, c, A, -d, B );
00684         VJOIN2( SUPERELLOUT(16), V, e, A, -f, B );
00685 }
00686 
00687 /**
00688  *                      R T _ S U P E R E L L _ P L O T
00689  */
00690 int
00691 rt_superell_plot(struct bu_list *vhead, struct rt_db_internal *ip, const struct rt_tess_tol *ttol, const struct bn_tol *tol)
00692 {
00693         register int            i;
00694         struct rt_superell_internal     *eip;
00695         fastf_t top[16*3];
00696         fastf_t middle[16*3];
00697         fastf_t bottom[16*3];
00698 
00699         RT_CK_DB_INTERNAL(ip);
00700         eip = (struct rt_superell_internal *)ip->idb_ptr;
00701         RT_SUPERELL_CK_MAGIC(eip);
00702 
00703         rt_superell_16pts( top, eip->v, eip->a, eip->b );
00704         rt_superell_16pts( bottom, eip->v, eip->b, eip->c );
00705         rt_superell_16pts( middle, eip->v, eip->a, eip->c );
00706 
00707         RT_ADD_VLIST( vhead, &top[15*ELEMENTS_PER_VECT], BN_VLIST_LINE_MOVE );
00708         for( i=0; i<16; i++ )  {
00709                 RT_ADD_VLIST( vhead, &top[i*ELEMENTS_PER_VECT], BN_VLIST_LINE_DRAW );
00710         }
00711 
00712         RT_ADD_VLIST( vhead, &bottom[15*ELEMENTS_PER_VECT], BN_VLIST_LINE_MOVE );
00713         for( i=0; i<16; i++ )  {
00714                 RT_ADD_VLIST( vhead, &bottom[i*ELEMENTS_PER_VECT], BN_VLIST_LINE_DRAW );
00715         }
00716 
00717         RT_ADD_VLIST( vhead, &middle[15*ELEMENTS_PER_VECT], BN_VLIST_LINE_MOVE );
00718         for( i=0; i<16; i++ )  {
00719                 RT_ADD_VLIST( vhead, &middle[i*ELEMENTS_PER_VECT], BN_VLIST_LINE_DRAW );
00720         }
00721         return(0);
00722 }
00723 
00724 #if 0
00725 static point_t  octa_verts[6] = {
00726         { 1, 0, 0 },    /* XPLUS */
00727         {-1, 0, 0 },    /* XMINUS */
00728         { 0, 1, 0 },    /* YPLUS */
00729         { 0,-1, 0 },    /* YMINUS */
00730         { 0, 0, 1 },    /* ZPLUS */
00731         { 0, 0,-1 }     /* ZMINUS */
00732 };
00733 
00734 #define XPLUS 0
00735 #define XMIN  1
00736 #define YPLUS 2
00737 #define YMIN  3
00738 #define ZPLUS 4
00739 #define ZMIN  5
00740 #endif
00741 
00742 /* Vertices of a unit octahedron */
00743 /* These need to be organized properly to give reasonable normals */
00744 /* static struct usvert {
00745  *      int     a;
00746  *      int     b;
00747  *      int     c;
00748  * } octahedron[8] = {
00749  *     { XPLUS, ZPLUS, YPLUS },
00750  *     { YPLUS, ZPLUS, XMIN  },
00751  *     { XMIN , ZPLUS, YMIN  },
00752  *     { YMIN , ZPLUS, XPLUS },
00753  *     { XPLUS, YPLUS, ZMIN  },
00754  *     { YPLUS, XMIN , ZMIN  },
00755  *     { XMIN , YMIN , ZMIN  },
00756  *     { YMIN , XPLUS, ZMIN  }
00757  * };
00758  */
00759 struct superell_state {
00760         struct shell    *s;
00761         struct rt_superell_internal     *eip;
00762         mat_t           invRinvS;
00763         mat_t           invRoS;
00764         fastf_t         theta_tol;
00765 };
00766 
00767 struct superell_vert_strip {
00768         int             nverts_per_strip;
00769         int             nverts;
00770         struct vertex   **vp;
00771         vect_t          *norms;
00772         int             nfaces;
00773         struct faceuse  **fu;
00774 };
00775 
00776 /**
00777  *                      R T _ S U P E R E L L _ T E S S
00778  *
00779  *  Tesssuperellate an superellipsoid.
00780  *
00781  *  The strategy is based upon the approach of Jon Leech 3/24/89,
00782  *  from program "sphere", which generates a polygon mesh
00783  *  approximating a sphere by
00784  *  recursive subdivision. First approximation is an octahedron;
00785  *  each level of refinement increases the number of polygons by
00786  *  a factor of 4.
00787  *  Level 3 (128 polygons) is a good tradeoff if gouraud
00788  *  shading is used to render the database.
00789  *
00790  *  At the start, points ABC lie on surface of the unit sphere.
00791  *  Pick DEF as the midpoints of the three edges of ABC.
00792  *  Normalize the new points to lie on surface of the unit sphere.
00793  *
00794  *        1
00795  *        B
00796  *       /\
00797  *    3 /  \ 4
00798  *    D/____\E
00799  *    /\    /\
00800  *   /  \  /  \
00801  *  /____\/____\
00802  * A      F     C
00803  * 0      5     2
00804  *
00805  *  Returns -
00806  *      -1      failure
00807  *       0      OK.  *r points to nmgregion that holds this tesssuperellation.
00808  */
00809 int
00810 rt_superell_tess(struct nmgregion **r, struct model *m, struct rt_db_internal *ip, const struct rt_tess_tol *ttol, const struct bn_tol *tol)
00811 {
00812   bu_log("rt_superell_tess called!\n");
00813   return -1;
00814 }
00815 
00816 /**
00817  *                      R T _ S U P E R E L L _ I M P O R T
00818  *
00819  *  Import an superellipsoid/sphere from the database format to
00820  *  the internal structure.
00821  *  Apply modeling transformations as wsuperell.
00822  */
00823 int
00824 rt_superell_import(struct rt_db_internal *ip, const struct bu_external *ep, register const fastf_t *mat, const struct db_i *dbip)
00825 {
00826         struct rt_superell_internal     *eip;
00827         union record            *rp;
00828         LOCAL fastf_t   vec[3*4 + 2];
00829 
00830         BU_CK_EXTERNAL( ep );
00831         rp = (union record *)ep->ext_buf;
00832         /* Check record type */
00833         if( rp->u_id != ID_SOLID )  {
00834                 bu_log("rt_superell_import: defective record\n");
00835                 return(-1);
00836         }
00837 
00838         RT_CK_DB_INTERNAL( ip );
00839         ip->idb_major_type = DB5_MAJORTYPE_BRLCAD;
00840         ip->idb_type = ID_SUPERELL;
00841         ip->idb_meth = &rt_functab[ID_SUPERELL];
00842         ip->idb_ptr = bu_malloc( sizeof(struct rt_superell_internal), "rt_superell_internal");
00843         eip = (struct rt_superell_internal *)ip->idb_ptr;
00844         eip->magic = RT_SUPERELL_INTERNAL_MAGIC;
00845 
00846         /* Convert from database to internal format */
00847         rt_fastf_float( vec, rp->s.s_values, 4 );
00848 
00849         /* Apply modeling transformations */
00850         MAT4X3PNT( eip->v, mat, &vec[0*3] );
00851         MAT4X3VEC( eip->a, mat, &vec[1*3] );
00852         MAT4X3VEC( eip->b, mat, &vec[2*3] );
00853         MAT4X3VEC( eip->c, mat, &vec[3*3] );
00854         eip->n = rp->s.s_values[12];
00855         eip->e = rp->s.s_values[13];
00856 
00857         return(0);              /* OK */
00858 }
00859 
00860 /**
00861  *                      R T _ S U P E R E L L _ E X P O R T
00862  */
00863 int
00864 rt_superell_export(struct bu_external *ep, const struct rt_db_internal *ip, double local2mm, const struct db_i *dbip)
00865 {
00866         struct rt_superell_internal     *tip;
00867         union record            *rec;
00868 
00869         RT_CK_DB_INTERNAL(ip);
00870         if( ip->idb_type != ID_SUPERELL )  return(-1);
00871         tip = (struct rt_superell_internal *)ip->idb_ptr;
00872         RT_SUPERELL_CK_MAGIC(tip);
00873 
00874         BU_CK_EXTERNAL(ep);
00875         ep->ext_nbytes = sizeof(union record);
00876         ep->ext_buf = (genptr_t)bu_calloc( 1, ep->ext_nbytes, "superell external");
00877         rec = (union record *)ep->ext_buf;
00878 
00879         rec->s.s_id = ID_SOLID;
00880         rec->s.s_type = SUPERELL;
00881 
00882         /* NOTE: This also converts to dbfloat_t */
00883         VSCALE( &rec->s.s_values[0], tip->v, local2mm );
00884         VSCALE( &rec->s.s_values[3], tip->a, local2mm );
00885         VSCALE( &rec->s.s_values[6], tip->b, local2mm );
00886         VSCALE( &rec->s.s_values[9], tip->c, local2mm );
00887 
00888         printf("SUPERELL: %g %g\n", tip->n, tip->e);
00889 
00890         rec->s.s_values[12] = tip->n;
00891         rec->s.s_values[13] = tip->e;
00892 
00893         return(0);
00894 }
00895 
00896 /**
00897  *                      R T _ S U P E R E L L _ I M P O R T 5
00898  *
00899  *  Import an superellipsoid/sphere from the database format to
00900  *  the internal structure.
00901  *  Apply modeling transformations as wsuperell.
00902  */
00903 int
00904 rt_superell_import5(struct rt_db_internal *ip, const struct bu_external *ep, register const fastf_t *mat, const struct db_i *dbip)
00905 {
00906         struct rt_superell_internal     *eip;
00907         fastf_t                 vec[ELEMENTS_PER_VECT*4 + 2];
00908 
00909         RT_CK_DB_INTERNAL( ip );
00910         BU_CK_EXTERNAL( ep );
00911 
00912         BU_ASSERT_LONG( ep->ext_nbytes, ==, SIZEOF_NETWORK_DOUBLE * (ELEMENTS_PER_VECT*4 + 2));
00913 
00914         ip->idb_major_type = DB5_MAJORTYPE_BRLCAD;
00915         ip->idb_type = ID_SUPERELL;
00916         ip->idb_meth = &rt_functab[ID_SUPERELL];
00917         ip->idb_ptr = bu_malloc( sizeof(struct rt_superell_internal), "rt_superell_internal");
00918 
00919         eip = (struct rt_superell_internal *)ip->idb_ptr;
00920         eip->magic = RT_SUPERELL_INTERNAL_MAGIC;
00921 
00922         /* Convert from database (network) to internal (host) format */
00923         ntohd( (unsigned char *)vec, ep->ext_buf, ELEMENTS_PER_VECT*4 + 2);
00924 
00925         /* Apply modeling transformations */
00926         MAT4X3PNT( eip->v, mat, &vec[0*ELEMENTS_PER_VECT] );
00927         MAT4X3VEC( eip->a, mat, &vec[1*ELEMENTS_PER_VECT] );
00928         MAT4X3VEC( eip->b, mat, &vec[2*ELEMENTS_PER_VECT] );
00929         MAT4X3VEC( eip->c, mat, &vec[3*ELEMENTS_PER_VECT] );
00930         eip->n = vec[4*ELEMENTS_PER_VECT];
00931         eip->e = vec[4*ELEMENTS_PER_VECT + 1];
00932 
00933         return 0;               /* OK */
00934 }
00935 
00936 /**
00937  *                      R T _ S U P E R E L L _ E X P O R T 5
00938  *
00939  *  The external format is:
00940  *      V point
00941  *      A vector
00942  *      B vector
00943  *      C vector
00944  */
00945 int
00946 rt_superell_export5(struct bu_external *ep, const struct rt_db_internal *ip, double local2mm, const struct db_i *dbip)
00947 {
00948         struct rt_superell_internal     *eip;
00949         fastf_t                 vec[ELEMENTS_PER_VECT*4 + 2];
00950 
00951         RT_CK_DB_INTERNAL(ip);
00952         if( ip->idb_type != ID_SUPERELL )  return(-1);
00953         eip = (struct rt_superell_internal *)ip->idb_ptr;
00954         RT_SUPERELL_CK_MAGIC(eip);
00955 
00956         BU_CK_EXTERNAL(ep);
00957         ep->ext_nbytes = SIZEOF_NETWORK_DOUBLE * (ELEMENTS_PER_VECT*4 + 2);
00958         ep->ext_buf = (genptr_t)bu_malloc( ep->ext_nbytes, "superell external");
00959 
00960         /* scale 'em into local buffer */
00961         VSCALE( &vec[0*ELEMENTS_PER_VECT], eip->v, local2mm );
00962         VSCALE( &vec[1*ELEMENTS_PER_VECT], eip->a, local2mm );
00963         VSCALE( &vec[2*ELEMENTS_PER_VECT], eip->b, local2mm );
00964         VSCALE( &vec[3*ELEMENTS_PER_VECT], eip->c, local2mm );
00965 
00966         vec[4*ELEMENTS_PER_VECT] = eip->n;
00967         vec[4*ELEMENTS_PER_VECT + 1] = eip->e;
00968 
00969         /* Convert from internal (host) to database (network) format */
00970         htond( ep->ext_buf, (unsigned char *)vec, ELEMENTS_PER_VECT*4 + 2 );
00971 
00972         return 0;
00973 }
00974 
00975 /**
00976  *                      R T _ S U P E R E L L _ D E S C R I B E
00977  *
00978  *  Make human-readable formatted presentation of this solid.
00979  *  First line describes type of solid.
00980  *  Additional lines are indented one tab, and give parameter values.
00981  */
00982 int
00983 rt_superell_describe(struct bu_vls *str, const struct rt_db_internal *ip, int verbose, double mm2local)
00984 {
00985         register struct rt_superell_internal    *tip =
00986                 (struct rt_superell_internal *)ip->idb_ptr;
00987         fastf_t mag_a, mag_b, mag_c;
00988         char    buf[256];
00989         double  angles[5];
00990         vect_t  unitv;
00991 
00992         RT_SUPERELL_CK_MAGIC(tip);
00993         bu_vls_strcat( str, "superellipsoid (SUPERELL)\n");
00994 
00995         sprintf(buf, "\tV (%g, %g, %g)\n",
00996                 INTCLAMP(tip->v[X] * mm2local),
00997                 INTCLAMP(tip->v[Y] * mm2local),
00998                 INTCLAMP(tip->v[Z] * mm2local) );
00999         bu_vls_strcat( str, buf );
01000 
01001         mag_a = MAGNITUDE(tip->a);
01002         mag_b = MAGNITUDE(tip->b);
01003         mag_c = MAGNITUDE(tip->c);
01004 
01005         sprintf(buf, "\tA (%g, %g, %g) mag=%g\n",
01006                 INTCLAMP(tip->a[X] * mm2local),
01007                 INTCLAMP(tip->a[Y] * mm2local),
01008                 INTCLAMP(tip->a[Z] * mm2local),
01009                 INTCLAMP(mag_a * mm2local) );
01010         bu_vls_strcat( str, buf );
01011 
01012         sprintf(buf, "\tB (%g, %g, %g) mag=%g\n",
01013                 INTCLAMP(tip->b[X] * mm2local),
01014                 INTCLAMP(tip->b[Y] * mm2local),
01015                 INTCLAMP(tip->b[Z] * mm2local),
01016                 INTCLAMP(mag_b * mm2local) );
01017         bu_vls_strcat( str, buf );
01018 
01019         sprintf(buf, "\tC (%g, %g, %g) mag=%g\n",
01020                 INTCLAMP(tip->c[X] * mm2local),
01021                 INTCLAMP(tip->c[Y] * mm2local),
01022                 INTCLAMP(tip->c[Z] * mm2local),
01023                 INTCLAMP(mag_c * mm2local) );
01024         bu_vls_strcat( str, buf );
01025 
01026         sprintf(buf, "\t<n,e> (%g, %g)\n", INTCLAMP(tip->n), INTCLAMP(tip->e));
01027         bu_vls_strcat(str, buf);
01028 
01029         if( !verbose )  return(0);
01030 
01031         VSCALE( unitv, tip->a, 1/mag_a );
01032         rt_find_fallback_angle( angles, unitv );
01033         rt_pr_fallback_angle( str, "\tA", angles );
01034 
01035         VSCALE( unitv, tip->b, 1/mag_b );
01036         rt_find_fallback_angle( angles, unitv );
01037         rt_pr_fallback_angle( str, "\tB", angles );
01038 
01039         VSCALE( unitv, tip->c, 1/mag_c );
01040         rt_find_fallback_angle( angles, unitv );
01041         rt_pr_fallback_angle( str, "\tC", angles );
01042 
01043         return(0);
01044 }
01045 
01046 /**
01047  *                      R T _ S U P E R E L L _ I F R E E
01048  *
01049  *  Free the storage associated with the rt_db_internal version of this solid.
01050  */
01051 void
01052 rt_superell_ifree(struct rt_db_internal *ip)
01053 {
01054         RT_CK_DB_INTERNAL(ip);
01055         bu_free( ip->idb_ptr, "superell ifree" );
01056         ip->idb_ptr = GENPTR_NULL;
01057 }
01058 
01059 /*  The U parameter runs south to north.
01060  *  In order to orient loop CCW, need to start with 0,1-->0,0 transition
01061  *  at the south pole.
01062  */
01063 static const fastf_t rt_superell_uvw[5*ELEMENTS_PER_VECT] = {
01064         0, 1, 0,
01065         0, 0, 0,
01066         1, 0, 0,
01067         1, 1, 0,
01068         0, 1, 0
01069 };
01070 
01071 /**
01072  *                      R T _ S U P E R E L L _ T N U R B
01073  */
01074 int
01075 rt_superell_tnurb(struct nmgregion **r, struct model *m, struct rt_db_internal *ip, const struct bn_tol *tol)
01076 {
01077   bu_log("rt_superell_tnurb called!\n");
01078         return 0;
01079 }
01080 
01081 
01082 /*
01083  * Local Variables:
01084  * mode: C
01085  * tab-width: 8
01086  * c-basic-offset: 4
01087  * indent-tabs-mode: t
01088  * End:
01089  * ex: shiftwidth=4 tabstop=8
01090  */

Generated on Mon Sep 18 01:24:51 2006 for BRL-CAD by  doxygen 1.4.6