g_rhc.c

Go to the documentation of this file.
00001 /*                         G _ R H C . C
00002  * BRL-CAD
00003  *
00004  * Copyright (c) 1990-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_rhc.c
00026  *      Intersect a ray with a Right Hyperbolic Cylinder.
00027  *
00028  *  Algorithm -
00029  *
00030  *  Given V, H, R, and B, there is a set of points on this rhc
00031  *
00032  *  { (x,y,z) | (x,y,z) is on rhc }
00033  *
00034  *  Through a series of Affine Transformations, this set of points will be
00035  *  transformed into a set of points on an rhc located at the origin
00036  *  with a rectangular halfwidth R of 1 along the Y axis, a height H of +1
00037  *  along the -X axis, a distance B of 1 along the -Z axis between the
00038  *  vertex V and the tip of the hyperbola, and a distance c between the
00039  *  tip of the hyperbola and the vertex of the asymptotic cone.
00040  *
00041  *
00042  *  { (x',y',z') | (x',y',z') is on rhc at origin }
00043  *
00044  *  The transformation from X to X' is accomplished by:
00045  *
00046  *  X' = S(R( X - V ))
00047  *
00048  *  where R(X) =  ( H/(-|H|) )
00049  *               (  R/( |R|)  ) . X
00050  *                ( B/(-|B|) )
00051  *
00052  *  and S(X) =   (  1/|H|   0     0   )
00053  *              (    0    1/|R|   0    ) . X
00054  *               (   0      0   1/|B| )
00055  *
00056  *  To find the intersection of a line with the surface of the rhc,
00057  *  consider the parametric line L:
00058  *
00059  *      L : { P(n) | P + t(n) . D }
00060  *
00061  *  Call W the actual point of intersection between L and the rhc.
00062  *  Let W' be the point of intersection between L' and the unit rhc.
00063  *
00064  *      L' : { P'(n) | P' + t(n) . D' }
00065  *
00066  *  W = invR( invS( W' ) ) + V
00067  *
00068  *  Where W' = k D' + P'.
00069  *
00070  *  If Dy' and Dz' are both 0, then there is no hit on the rhc;
00071  *  but the end plates need checking.  If there is now only 1 hit
00072  *  point, the top plate needs to be checked as well.
00073  *
00074  *  Line L' hits the infinitely long canonical rhc at W' when
00075  *
00076  *      A * k**2 + B * k + C = 0
00077  *
00078  *  where
00079  *
00080  *  A = Dz'**2 - Dy'**2 * (1 + 2*c')
00081  *  B = 2 * ((1 + c' + Pz') * Dz' - (1 + 2*c') * Dy' * Py'
00082  *  C = (Pz' + c' + 1)**2 - (1 + 2*c') * Py'**2 - c'**2
00083  *  b = |Breadth| = 1.0
00084  *  h = |Height| = 1.0
00085  *  r = 1.0
00086  *  c' = c / |Breadth|
00087  *
00088  *  The quadratic formula yields k (which is constant):
00089  *
00090  *  k = [ -B +/- sqrt( B**2 - 4 * A * C )] / (2 * A)
00091  *
00092  *  Now, D' = S( R( D ) )
00093  *  and  P' = S( R( P - V ) )
00094  *
00095  *  Substituting,
00096  *
00097  *  W = V + invR( invS[ k *( S( R( D ) ) ) + S( R( P - V ) ) ] )
00098  *    = V + invR( ( k * R( D ) ) + R( P - V ) )
00099  *    = V + k * D + P - V
00100  *    = k * D + P
00101  *
00102  *  Note that ``k'' is constant, and is the same in the formulations
00103  *  for both W and W'.
00104  *
00105  *  The hit at ``k'' is a hit on the canonical rhc IFF
00106  *  -1 <= Wx' <= 0 and -1 <= Wz' <= 0.
00107  *
00108  *  NORMALS.  Given the point W on the surface of the rhc,
00109  *  what is the vector normal to the tangent plane at that point?
00110  *
00111  *  Map W onto the unit rhc, ie:  W' = S( R( W - V ) ).
00112  *
00113  *  Plane on unit rhc at W' has a normal vector N' where
00114  *
00115  *  N' = <0, Wy'*(1 + 2*c), -z-c-1>.
00116  *
00117  *  The plane transforms back to the tangent plane at W, and this
00118  *  new plane (on the original rhc) has a normal vector of N, viz:
00119  *
00120  *  N = inverse[ transpose( inverse[ S o R ] ) ] ( N' )
00121  *
00122  *  because if H is perpendicular to plane Q, and matrix M maps from
00123  *  Q to Q', then inverse[ transpose(M) ] (H) is perpendicular to Q'.
00124  *  Here, H and Q are in "prime space" with the unit sphere.
00125  *  [Somehow, the notation here is backwards].
00126  *  So, the mapping matrix M = inverse( S o R ), because
00127  *  S o R maps from normal space to the unit sphere.
00128  *
00129  *  N = inverse[ transpose( inverse[ S o R ] ) ] ( N' )
00130  *    = inverse[ transpose(invR o invS) ] ( N' )
00131  *    = inverse[ transpose(invS) o transpose(invR) ] ( N' )
00132  *    = inverse[ inverse(S) o R ] ( N' )
00133  *    = invR o S ( N' )
00134  *
00135  *  because inverse(R) = transpose(R), so R = transpose( invR ),
00136  *  and S = transpose( S ).
00137  *
00138  *  Note that the normal vector produced above will not have unit length.
00139  *
00140  *  THE TOP AND END PLATES.
00141  *
00142  *  If Dz' == 0, line L' is parallel to the top plate, so there is no
00143  *  hit on the top plate.  Otherwise, rays intersect the top plate
00144  *  with k = (0 - Pz')/Dz'.  The solution is within the top plate
00145  *  IFF  -1 <= Wx' <= 0 and -1 <= Wy' <= 1.
00146  *
00147  *  If Dx' == 0, line L' is parallel to the end plates, so there is no
00148  *  hit on the end plates.  Otherwise, rays intersect the front plate
00149  *  line L' hits the front plate with k = (0 - Px') / Dx', and
00150  *  and hits the back plate with k = (-1 - Px') / Dx'.
00151  *
00152  *  The solution W' is within an end plate IFF
00153  *
00154  *      (Wz' + c + 1)**2 - (Wy'**2 * (2*c + 1) >= c**2  and  Wz' <= 1.0
00155  *
00156  *  The normal for a hit on the top plate is -Bunit.
00157  *  The normal for a hit on the front plate is -Hunit, and
00158  *  the normal for a hit on the back plate is +Hunit.
00159  *
00160  *  Authors -
00161  *      Michael J. Markowski
00162  *
00163  *  Source -
00164  *      SECAD/VLD Computing Consortium, Bldg 394
00165  *      The U. S. Army Ballistic Research Laboratory
00166  *      Aberdeen Proving Ground, Maryland  21005-5066
00167  *
00168  */
00169 /*@}*/
00170 
00171 #ifndef lint
00172 static const char RCSrhc[] = "@(#)$Header: /cvsroot/brlcad/brlcad/src/librt/g_rhc.c,v 14.13 2006/09/16 02:04:24 lbutler Exp $ (BRL)";
00173 #endif
00174 
00175 #include "common.h"
00176 
00177 #include <stddef.h>
00178 #include <stdio.h>
00179 #ifdef HAVE_STRING_H
00180 #  include <string.h>
00181 #else
00182 #  include <strings.h>
00183 #endif
00184 #include <math.h>
00185 
00186 #include "machine.h"
00187 #include "vmath.h"
00188 #include "db.h"
00189 #include "nmg.h"
00190 #include "raytrace.h"
00191 #include "rtgeom.h"
00192 #include "./debug.h"
00193 
00194 struct rhc_specific {
00195         point_t rhc_V;          /* vector to rhc origin */
00196         vect_t  rhc_Bunit;      /* unit B vector */
00197         vect_t  rhc_Hunit;      /* unit H vector */
00198         vect_t  rhc_Runit;      /* unit vector, B x H */
00199         mat_t   rhc_SoR;        /* Scale(Rot(vect)) */
00200         mat_t   rhc_invRoS;     /* invRot(Scale(vect)) */
00201         fastf_t rhc_b;          /* |B| */
00202         fastf_t rhc_c;          /* c */
00203         fastf_t rhc_cprime;     /* c / |B| */
00204         fastf_t rhc_rsq;        /* r * r */
00205 };
00206 
00207 const struct bu_structparse rt_rhc_parse[] = {
00208     { "%f", 3, "V", bu_offsetof(struct rt_rhc_internal, rhc_V[X]), BU_STRUCTPARSE_FUNC_NULL },
00209     { "%f", 3, "H", bu_offsetof(struct rt_rhc_internal, rhc_H[X]), BU_STRUCTPARSE_FUNC_NULL },
00210     { "%f", 3, "B", bu_offsetof(struct rt_rhc_internal, rhc_B[X]), BU_STRUCTPARSE_FUNC_NULL },
00211     { "%f", 1, "r", bu_offsetof(struct rt_rhc_internal, rhc_r),    BU_STRUCTPARSE_FUNC_NULL },
00212     { "%f", 1, "c", bu_offsetof(struct rt_rhc_internal, rhc_c),    BU_STRUCTPARSE_FUNC_NULL },
00213     { {'\0','\0','\0','\0'}, 0, (char *)NULL, 0, BU_STRUCTPARSE_FUNC_NULL }
00214  };
00215 
00216 /**
00217  *                      R T _ R H C _ P R E P
00218  *
00219  *  Given a pointer to a GED database record, and a transformation matrix,
00220  *  determine if this is a valid RHC, and if so, precompute various
00221  *  terms of the formula.
00222  *
00223  *  Returns -
00224  *      0       RHC is OK
00225  *      !0      Error in description
00226  *
00227  *  Implicit return -
00228  *      A struct rhc_specific is created, and it's address is stored in
00229  *      stp->st_specific for use by rhc_shot().
00230  */
00231 int
00232 rt_rhc_prep(struct soltab *stp, struct rt_db_internal *ip, struct rt_i *rtip)
00233 {
00234         struct rt_rhc_internal          *xip;
00235         register struct rhc_specific    *rhc;
00236         LOCAL fastf_t   magsq_b, magsq_h, magsq_r;
00237         LOCAL fastf_t   mag_b, mag_h, mag_r;
00238         LOCAL fastf_t   f;
00239         LOCAL mat_t     R;
00240         LOCAL mat_t     Rinv;
00241         LOCAL mat_t     S;
00242         LOCAL vect_t    invsq;  /* [ 1/(|H|**2), 1/(|R|**2), 1/(|B|**2) ] */
00243 
00244         RT_CK_DB_INTERNAL(ip);
00245         RT_CK_RTI(rtip);
00246         xip = (struct rt_rhc_internal *)ip->idb_ptr;
00247         RT_RHC_CK_MAGIC(xip);
00248 
00249         /* compute |B| |H| */
00250         mag_b = sqrt( magsq_b = MAGSQ( xip->rhc_B ) );
00251         mag_h = sqrt( magsq_h = MAGSQ( xip->rhc_H ) );
00252         mag_r = xip->rhc_r;
00253         magsq_r = mag_r * mag_r;
00254 
00255         /* Check for |H| > 0, |B| > 0, |R| > 0 */
00256         if( NEAR_ZERO(mag_h, RT_LEN_TOL) || NEAR_ZERO(mag_b, RT_LEN_TOL)
00257          || NEAR_ZERO(mag_r, RT_LEN_TOL) )  {
00258                 return(1);              /* BAD, too small */
00259         }
00260 
00261         /* Check for B.H == 0 */
00262         f = VDOT( xip->rhc_B, xip->rhc_H ) / (mag_b * mag_h);
00263         if( ! NEAR_ZERO(f, RT_DOT_TOL) )  {
00264                 return(1);              /* BAD */
00265         }
00266 
00267         /*
00268          *  RHC is ok
00269          */
00270         stp->st_id = ID_RHC;            /* set soltab ID */
00271         stp->st_meth = &rt_functab[ID_RHC];
00272 
00273         BU_GETSTRUCT( rhc, rhc_specific );
00274         stp->st_specific = (genptr_t)rhc;
00275         rhc->rhc_b = mag_b;
00276         rhc->rhc_rsq = magsq_r;
00277         rhc->rhc_c = xip->rhc_c;
00278 
00279         /* make unit vectors in B, H, and BxH directions */
00280         VMOVE(    rhc->rhc_Hunit, xip->rhc_H );
00281         VUNITIZE( rhc->rhc_Hunit );
00282         VMOVE(    rhc->rhc_Bunit, xip->rhc_B );
00283         VUNITIZE( rhc->rhc_Bunit );
00284         VCROSS(   rhc->rhc_Runit, rhc->rhc_Bunit, rhc->rhc_Hunit );
00285 
00286         VMOVE( rhc->rhc_V, xip->rhc_V );
00287         rhc->rhc_cprime = xip->rhc_c / mag_b;
00288 
00289         /* Compute R and Rinv matrices */
00290         MAT_IDN( R );
00291         VREVERSE( &R[0], rhc->rhc_Hunit );
00292         VMOVE(    &R[4], rhc->rhc_Runit );
00293         VREVERSE( &R[8], rhc->rhc_Bunit );
00294         bn_mat_trn( Rinv, R );                  /* inv of rot mat is trn */
00295 
00296         /* Compute S */
00297         VSET( invsq, 1.0/magsq_h, 1.0/magsq_r, 1.0/magsq_b );
00298         MAT_IDN( S );
00299         S[ 0] = sqrt( invsq[0] );
00300         S[ 5] = sqrt( invsq[1] );
00301         S[10] = sqrt( invsq[2] );
00302 
00303         /* Compute SoR and invRoS */
00304         bn_mat_mul( rhc->rhc_SoR, S, R );
00305         bn_mat_mul( rhc->rhc_invRoS, Rinv, S );
00306 
00307         /* Compute bounding sphere and RPP */
00308         /* bounding sphere center */
00309         VJOIN2( stp->st_center, rhc->rhc_V,
00310                 mag_h / 2.0,    rhc->rhc_Hunit,
00311                 mag_b / 2.0,    rhc->rhc_Bunit );
00312         /* bounding radius */
00313         stp->st_bradius = 0.5 * sqrt(magsq_h + 4.0*magsq_r + magsq_b);
00314         /* approximate bounding radius */
00315         stp->st_aradius = stp->st_bradius;
00316 
00317         /* cheat, make bounding RPP by enclosing bounding sphere */
00318         stp->st_min[X] = stp->st_center[X] - stp->st_bradius;
00319         stp->st_max[X] = stp->st_center[X] + stp->st_bradius;
00320         stp->st_min[Y] = stp->st_center[Y] - stp->st_bradius;
00321         stp->st_max[Y] = stp->st_center[Y] + stp->st_bradius;
00322         stp->st_min[Z] = stp->st_center[Z] - stp->st_bradius;
00323         stp->st_max[Z] = stp->st_center[Z] + stp->st_bradius;
00324 
00325         return(0);                      /* OK */
00326 }
00327 
00328 /**
00329  *                      R T _ R H C _ P R I N T
00330  */
00331 void
00332 rt_rhc_print(register const struct soltab *stp)
00333 {
00334         register const struct rhc_specific *rhc =
00335                 (struct rhc_specific *)stp->st_specific;
00336 
00337         VPRINT("V", rhc->rhc_V);
00338         VPRINT("Bunit", rhc->rhc_Bunit);
00339         VPRINT("Hunit", rhc->rhc_Hunit);
00340         VPRINT("Runit", rhc->rhc_Runit);
00341         bn_mat_print("S o R", rhc->rhc_SoR );
00342         bn_mat_print("invR o S", rhc->rhc_invRoS );
00343 }
00344 
00345 /* hit_surfno is set to one of these */
00346 #define RHC_NORM_BODY   (1)             /* compute normal */
00347 #define RHC_NORM_TOP    (2)             /* copy rhc_N */
00348 #define RHC_NORM_FRT    (3)             /* copy reverse rhc_N */
00349 #define RHC_NORM_BACK   (4)
00350 
00351 /**
00352  *                      R T _ R H C _ S H O T
00353  *
00354  *  Intersect a ray with a rhc.
00355  *  If an intersection occurs, a struct seg will be acquired
00356  *  and filled in.
00357  *
00358  *  Returns -
00359  *      0       MISS
00360  *      >0      HIT
00361  */
00362 int
00363 rt_rhc_shot(struct soltab *stp, register struct xray *rp, struct application *ap, struct seg *seghead)
00364 {
00365         register struct rhc_specific *rhc =
00366                 (struct rhc_specific *)stp->st_specific;
00367         LOCAL vect_t    dprime;         /* D' */
00368         LOCAL vect_t    pprime;         /* P' */
00369         LOCAL fastf_t   k1, k2;         /* distance constants of solution */
00370         LOCAL fastf_t   x;
00371         LOCAL vect_t    xlated;         /* translated vector */
00372         LOCAL struct hit hits[3];       /* 2 potential hit points */
00373         register struct hit *hitp;      /* pointer to hit point */
00374 
00375         hitp = &hits[0];
00376 
00377         /* out, Mat, vect */
00378         MAT4X3VEC( dprime, rhc->rhc_SoR, rp->r_dir );
00379         VSUB2( xlated, rp->r_pt, rhc->rhc_V );
00380         MAT4X3VEC( pprime, rhc->rhc_SoR, xlated );
00381 
00382         x = rhc->rhc_cprime;
00383 
00384         if ( NEAR_ZERO(dprime[Y], SMALL) && NEAR_ZERO(dprime[Z], SMALL) )
00385                 goto check_plates;
00386 
00387         /* Find roots of the equation, using formula for quadratic */
00388         {
00389                 FAST fastf_t    a, b, c;        /* coeffs of polynomial */
00390                 FAST fastf_t    disc;           /* disc of radical */
00391 
00392                 a = dprime[Z] * dprime[Z] - dprime[Y] * dprime[Y] * (1 + 2*x);
00393                 b = 2*((pprime[Z] + x + 1) * dprime[Z]
00394                         - (2*x + 1) * dprime[Y] * pprime[Y]);
00395                 c = (pprime[Z]+x+1)*(pprime[Z]+x+1)
00396                         - (2*x + 1) * pprime[Y] * pprime[Y] - x*x;
00397                 if ( !NEAR_ZERO(a, RT_PCOEF_TOL) ) {
00398                         disc = b*b - 4 * a * c;
00399                         if (disc <= 0)
00400                                 goto check_plates;
00401                         disc = sqrt(disc);
00402 
00403                         k1 = (-b + disc) / (2.0 * a);
00404                         k2 = (-b - disc) / (2.0 * a);
00405 
00406                         /*
00407                          *  k1 and k2 are potential solutions to intersection with side.
00408                          *  See if they fall in range.
00409                          */
00410                         VJOIN1( hitp->hit_vpriv, pprime, k1, dprime );          /* hit' */
00411                         if( hitp->hit_vpriv[X] >= -1.0
00412                                 && hitp->hit_vpriv[X] <= 0.0
00413                                 && hitp->hit_vpriv[Z] >= -1.0
00414                                 && hitp->hit_vpriv[Z] <= 0.0 ) {
00415                                 hitp->hit_magic = RT_HIT_MAGIC;
00416                                 hitp->hit_dist = k1;
00417                                 hitp->hit_surfno = RHC_NORM_BODY;       /* compute N */
00418                                 hitp++;
00419                         }
00420 
00421                         VJOIN1( hitp->hit_vpriv, pprime, k2, dprime );          /* hit' */
00422                         if( hitp->hit_vpriv[X] >= -1.0
00423                                 && hitp->hit_vpriv[X] <= 0.0
00424                                 && hitp->hit_vpriv[Z] >= -1.0
00425                                 && hitp->hit_vpriv[Z] <= 0.0 ) {
00426                                 hitp->hit_magic = RT_HIT_MAGIC;
00427                                 hitp->hit_dist = k2;
00428                                 hitp->hit_surfno = RHC_NORM_BODY;       /* compute N */
00429                                 hitp++;
00430                         }
00431                 } else if ( !NEAR_ZERO(b, RT_PCOEF_TOL) ) {
00432                         k1 = -c/b;
00433                         VJOIN1( hitp->hit_vpriv, pprime, k1, dprime );          /* hit' */
00434                         if( hitp->hit_vpriv[X] >= -1.0
00435                                 && hitp->hit_vpriv[X] <= 0.0
00436                                 && hitp->hit_vpriv[Z] >= -1.0
00437                                 && hitp->hit_vpriv[Z] <= 0.0 ) {
00438                                 hitp->hit_magic = RT_HIT_MAGIC;
00439                                 hitp->hit_dist = k1;
00440                                 hitp->hit_surfno = RHC_NORM_BODY;       /* compute N */
00441                                 hitp++;
00442                         }
00443                 }
00444         }
00445 
00446 
00447         /*
00448          * Check for hitting the top and end plates.
00449          */
00450 check_plates:
00451         /* check front and back plates */
00452         if( hitp < &hits[2]  &&  !NEAR_ZERO(dprime[X], SMALL) )  {
00453                 /* 0 or 1 hits so far, this is worthwhile */
00454                 k1 = -pprime[X] / dprime[X];            /* front plate */
00455                 k2 = (-1.0 - pprime[X]) / dprime[X];    /* back plate */
00456 
00457                 VJOIN1( hitp->hit_vpriv, pprime, k1, dprime );  /* hit' */
00458                 if( (hitp->hit_vpriv[Z] + x + 1.0)
00459                         * (hitp->hit_vpriv[Z] + x + 1.0)
00460                         - hitp->hit_vpriv[Y] * hitp->hit_vpriv[Y]
00461                         * (1.0 + 2*x) >= x*x
00462                         && hitp->hit_vpriv[Z] >= -1.0
00463                         && hitp->hit_vpriv[Z] <= 0.0)  {
00464                         hitp->hit_magic = RT_HIT_MAGIC;
00465                         hitp->hit_dist = k1;
00466                         hitp->hit_surfno = RHC_NORM_FRT;        /* -H */
00467                         hitp++;
00468                 }
00469 
00470                 VJOIN1( hitp->hit_vpriv, pprime, k2, dprime );  /* hit' */
00471                 if( (hitp->hit_vpriv[Z] + x + 1.0)
00472                         * (hitp->hit_vpriv[Z] + x + 1.0)
00473                         - hitp->hit_vpriv[Y] * hitp->hit_vpriv[Y]
00474                         * (1.0 + 2*x) >= x*x
00475                         && hitp->hit_vpriv[Z] >= -1.0
00476                         && hitp->hit_vpriv[Z] <= 0.0)  {
00477                         hitp->hit_magic = RT_HIT_MAGIC;
00478                         hitp->hit_dist = k2;
00479                         hitp->hit_surfno = RHC_NORM_BACK;       /* +H */
00480                         hitp++;
00481                 }
00482         }
00483 
00484         /* check top plate */
00485         if( hitp == &hits[1]  &&  !NEAR_ZERO(dprime[Z], SMALL) )  {
00486                 /* 0 or 1 hits so far, this is worthwhile */
00487                 k1 = -pprime[Z] / dprime[Z];            /* top plate */
00488 
00489                 VJOIN1( hitp->hit_vpriv, pprime, k1, dprime );  /* hit' */
00490                 if( hitp->hit_vpriv[X] >= -1.0 &&  hitp->hit_vpriv[X] <= 0.0
00491                         && hitp->hit_vpriv[Y] >= -1.0
00492                         && hitp->hit_vpriv[Y] <= 1.0 ) {
00493                         hitp->hit_magic = RT_HIT_MAGIC;
00494                         hitp->hit_dist = k1;
00495                         hitp->hit_surfno = RHC_NORM_TOP;        /* -B */
00496                         hitp++;
00497                 }
00498         }
00499 
00500         if( hitp != &hits[2] )
00501                 return(0);      /* MISS */
00502 
00503         if( hits[0].hit_dist < hits[1].hit_dist )  {
00504                 /* entry is [0], exit is [1] */
00505                 register struct seg *segp;
00506 
00507                 RT_GET_SEG(segp, ap->a_resource);
00508                 segp->seg_stp = stp;
00509                 segp->seg_in = hits[0];         /* struct copy */
00510                 segp->seg_out = hits[1];        /* struct copy */
00511                 BU_LIST_INSERT( &(seghead->l), &(segp->l) );
00512         } else {
00513                 /* entry is [1], exit is [0] */
00514                 register struct seg *segp;
00515 
00516                 RT_GET_SEG(segp, ap->a_resource);
00517                 segp->seg_stp = stp;
00518                 segp->seg_in = hits[1];         /* struct copy */
00519                 segp->seg_out = hits[0];        /* struct copy */
00520                 BU_LIST_INSERT( &(seghead->l), &(segp->l) );
00521         }
00522         return(2);                      /* HIT */
00523 }
00524 
00525 #define RT_RHC_SEG_MISS(SEG)    (SEG).seg_stp=RT_SOLTAB_NULL
00526 
00527 /**
00528  *                      R T _ R H C _ V S H O T
00529  *
00530  *  Vectorized version.
00531  */
00532 void
00533 rt_rhc_vshot(struct soltab **stp, struct xray **rp, struct seg *segp, int n, struct application *ap)
00534                                /* An array of solid pointers */
00535                                /* An array of ray pointers */
00536                                /* array of segs (results returned) */
00537                                /* Number of ray/object pairs */
00538 
00539 {
00540         rt_vstub( stp, rp, segp, n, ap );
00541 }
00542 
00543 /**
00544  *                      R T _ R H C _ N O R M
00545  *
00546  *  Given ONE ray distance, return the normal and entry/exit point.
00547  */
00548 void
00549 rt_rhc_norm(register struct hit *hitp, struct soltab *stp, register struct xray *rp)
00550 {
00551         fastf_t c;
00552         vect_t  can_normal;     /* normal to canonical rhc */
00553         register struct rhc_specific *rhc =
00554                 (struct rhc_specific *)stp->st_specific;
00555 
00556         VJOIN1( hitp->hit_point, rp->r_pt, hitp->hit_dist, rp->r_dir );
00557         switch( hitp->hit_surfno )  {
00558         case RHC_NORM_BODY:
00559                 c = rhc->rhc_cprime;
00560                 VSET( can_normal,
00561                         0.0,
00562                         hitp->hit_vpriv[Y] * (1.0 + 2.0*c),
00563                         -hitp->hit_vpriv[Z] - c - 1.0 );
00564                 MAT4X3VEC( hitp->hit_normal, rhc->rhc_invRoS, can_normal );
00565                 VUNITIZE( hitp->hit_normal );
00566                 break;
00567         case RHC_NORM_TOP:
00568                 VREVERSE( hitp->hit_normal, rhc->rhc_Bunit );
00569                 break;
00570         case RHC_NORM_FRT:
00571                 VREVERSE( hitp->hit_normal, rhc->rhc_Hunit );
00572                 break;
00573         case RHC_NORM_BACK:
00574                 VMOVE( hitp->hit_normal, rhc->rhc_Hunit );
00575                 break;
00576         default:
00577                 bu_log("rt_rhc_norm: surfno=%d bad\n", hitp->hit_surfno);
00578                 break;
00579         }
00580 }
00581 
00582 /**
00583  *                      R T _ R H C _ C U R V E
00584  *
00585  *  Return the curvature of the rhc.
00586  */
00587 void
00588 rt_rhc_curve(register struct curvature *cvp, register struct hit *hitp, struct soltab *stp)
00589 {
00590         fastf_t b, c, rsq, y;
00591         fastf_t zp1_sq, zp2;    /* 1st deriv sqr, 2nd deriv */
00592         register struct rhc_specific *rhc =
00593                 (struct rhc_specific *)stp->st_specific;
00594 
00595         switch( hitp->hit_surfno )  {
00596         case RHC_NORM_BODY:
00597                 /* most nearly flat direction */
00598                 VMOVE( cvp->crv_pdir, rhc->rhc_Hunit );
00599                 cvp->crv_c1 = 0;
00600                 /* k = z'' / (1 + z'^2) ^ 3/2 */
00601                 b = rhc->rhc_b;
00602                 c = rhc->rhc_c;
00603                 y = hitp->hit_point[Y];
00604                 rsq = rhc->rhc_rsq;
00605                 zp1_sq = y * (b*b + 2*b*c) / rsq;
00606                 zp1_sq *= zp1_sq / (c*c + y*y*(b*b + 2*b*c)/rsq);
00607                 zp2 = c*c / (rsq*c*c + y*y*(b*b + 2*b*c));
00608                 cvp->crv_c2 = zp2 / pow( (1 + zp1_sq), 1.5);
00609                 break;
00610         case RHC_NORM_BACK:
00611         case RHC_NORM_FRT:
00612         case RHC_NORM_TOP:
00613                 /* any tangent direction */
00614                 bn_vec_ortho( cvp->crv_pdir, hitp->hit_normal );
00615                 cvp->crv_c1 = cvp->crv_c2 = 0;
00616                 break;
00617         }
00618 }
00619 
00620 /**
00621  *                      R T _ R H C _ U V
00622  *
00623  *  For a hit on the surface of an rhc, return the (u,v) coordinates
00624  *  of the hit point, 0 <= u,v <= 1.
00625  *  u = azimuth
00626  *  v = elevation
00627  */
00628 void
00629 rt_rhc_uv(struct application *ap, struct soltab *stp, register struct hit *hitp, register struct uvcoord *uvp)
00630 {
00631         register struct rhc_specific *rhc =
00632                 (struct rhc_specific *)stp->st_specific;
00633         LOCAL vect_t work;
00634         LOCAL vect_t pprime;
00635         FAST fastf_t len;
00636 
00637         /*
00638          * hit_point is on surface;  project back to unit rhc,
00639          * creating a vector from vertex to hit point.
00640          */
00641         VSUB2( work, hitp->hit_point, rhc->rhc_V );
00642         MAT4X3VEC( pprime, rhc->rhc_SoR, work );
00643 
00644         switch( hitp->hit_surfno )  {
00645         case RHC_NORM_BODY:
00646                 /* Skin.  x,y coordinates define rotation.  radius = 1 */
00647                 len = sqrt(pprime[Y]*pprime[Y] + pprime[Z]*pprime[Z]);
00648                 uvp->uv_u = acos(pprime[Y]/len) * bn_invpi;
00649                 uvp->uv_v = -pprime[X];         /* height */
00650                 break;
00651         case RHC_NORM_FRT:
00652         case RHC_NORM_BACK:
00653                 /* end plates - circular mapping, not seamless w/body, top */
00654                 len = sqrt(pprime[Y]*pprime[Y] + pprime[Z]*pprime[Z]);
00655                 uvp->uv_u = acos(pprime[Y]/len) * bn_invpi;
00656                 uvp->uv_v = len;        /* rim v = 1 for both plates */
00657                 break;
00658         case RHC_NORM_TOP:
00659                 uvp->uv_u = 1.0 - (pprime[Y] + 1.0)/2.0;
00660                 uvp->uv_v = -pprime[X];         /* height */
00661                 break;
00662         }
00663 
00664         /* uv_du should be relative to rotation, uv_dv relative to height */
00665         uvp->uv_du = uvp->uv_dv = 0;
00666 }
00667 
00668 /**
00669  *              R T _ R H C _ F R E E
00670  */
00671 void
00672 rt_rhc_free(register struct soltab *stp)
00673 {
00674         register struct rhc_specific *rhc =
00675                 (struct rhc_specific *)stp->st_specific;
00676 
00677         bu_free( (char *)rhc, "rhc_specific" );
00678 }
00679 
00680 /**
00681  *                      R T _ R H C _ C L A S S
00682  */
00683 int
00684 rt_rhc_class(void)
00685 {
00686         return(0);
00687 }
00688 
00689 
00690 /**
00691  *                      R T _ R H C _ P L O T
00692  */
00693 int
00694 rt_rhc_plot(struct bu_list *vhead, struct rt_db_internal *ip, const struct rt_tess_tol *ttol, const struct bn_tol *tol)
00695 {
00696         int             i, n;
00697         fastf_t         b, c, *back, f, *front, h, rh;
00698         fastf_t         dtol, ntol;
00699         vect_t          Bu, Hu, Ru;
00700         LOCAL mat_t     R;
00701         LOCAL mat_t     invR;
00702         LOCAL struct rt_rhc_internal    *xip;
00703         struct rt_pt_node       *old, *pos, *pts, *rt_ptalloc(void);
00704 
00705         RT_CK_DB_INTERNAL(ip);
00706         xip = (struct rt_rhc_internal *)ip->idb_ptr;
00707         RT_RHC_CK_MAGIC(xip);
00708 
00709         /* compute |B| |H| */
00710         b = MAGNITUDE( xip->rhc_B );    /* breadth */
00711         rh = xip->rhc_r;                /* rectangular halfwidth */
00712         h = MAGNITUDE( xip->rhc_H );    /* height */
00713         c = xip->rhc_c;                 /* dist to asympt origin */
00714 
00715         /* Check for |H| > 0, |B| > 0, rh > 0, c > 0 */
00716         if( NEAR_ZERO(h, RT_LEN_TOL) || NEAR_ZERO(b, RT_LEN_TOL)
00717          || NEAR_ZERO(rh, RT_LEN_TOL) || NEAR_ZERO(c, RT_LEN_TOL))  {
00718                 bu_log("rt_rhc_plot:  zero length H, B, c, or rh\n");
00719                 return(-2);             /* BAD */
00720         }
00721 
00722         /* Check for B.H == 0 */
00723         f = VDOT( xip->rhc_B, xip->rhc_H ) / (b * h);
00724         if( ! NEAR_ZERO(f, RT_DOT_TOL) )  {
00725                 bu_log("rt_rhc_plot: B not perpendicular to H, f=%f\n", f);
00726                 return(-3);             /* BAD */
00727         }
00728 
00729         /* make unit vectors in B, H, and BxH directions */
00730         VMOVE(    Hu, xip->rhc_H );
00731         VUNITIZE( Hu );
00732         VMOVE(    Bu, xip->rhc_B );
00733         VUNITIZE( Bu );
00734         VCROSS(   Ru, Bu, Hu );
00735 
00736         /* Compute R and Rinv matrices */
00737         MAT_IDN( R );
00738         VREVERSE( &R[0], Hu );
00739         VMOVE(    &R[4], Ru );
00740         VREVERSE( &R[8], Bu );
00741         bn_mat_trn( invR, R );                  /* inv of rot mat is trn */
00742 
00743         /*
00744          *  Establish tolerances
00745          */
00746         if( ttol->rel <= 0.0 || ttol->rel >= 1.0 )  {
00747                 dtol = 0.0;             /* none */
00748         } else {
00749                 /* Convert rel to absolute by scaling by smallest side */
00750                 if (rh < b)
00751                         dtol = ttol->rel * 2 * rh;
00752                 else
00753                         dtol = ttol->rel * 2 * b;
00754         }
00755         if( ttol->abs <= 0.0 )  {
00756                 if( dtol <= 0.0 )  {
00757                         /* No tolerance given, use a default */
00758                         if (rh < b)
00759                                 dtol = 2 * 0.10 * rh;   /* 10% */
00760                         else
00761                                 dtol = 2 * 0.10 * b;    /* 10% */
00762                 } else {
00763                         /* Use absolute-ized relative tolerance */
00764                 }
00765         } else {
00766                 /* Absolute tolerance was given, pick smaller */
00767                 if( ttol->rel <= 0.0 || dtol > ttol->abs )
00768                         dtol = ttol->abs;
00769         }
00770 
00771         /* To ensure normal tolerance, remain below this angle */
00772         if( ttol->norm > 0.0 )
00773                 ntol = ttol->norm;
00774         else
00775                 /* tolerate everything */
00776                 ntol = bn_pi;
00777 
00778         /* initial hyperbola approximation is a single segment */
00779         pts = rt_ptalloc();
00780         pts->next = rt_ptalloc();
00781         pts->next->next = NULL;
00782         VSET( pts->p,       0, -rh, 0);
00783         VSET( pts->next->p, 0,  rh, 0);
00784         /* 2 endpoints in 1st approximation */
00785         n = 2;
00786         /* recursively break segment 'til within error tolerances */
00787         n += rt_mk_hyperbola( pts, rh, b, c, dtol, ntol );
00788 
00789         /* get mem for arrays */
00790         front = (fastf_t *)bu_malloc(3*n * sizeof(fastf_t), "fast_t");
00791         back  = (fastf_t *)bu_malloc(3*n * sizeof(fastf_t), "fast_t");
00792 
00793         /* generate front & back plates in world coordinates */
00794         pos = pts;
00795         i = 0;
00796         while (pos) {
00797                 /* rotate back to original position */
00798                 MAT4X3VEC( &front[i], invR, pos->p );
00799                 /* move to origin vertex origin */
00800                 VADD2( &front[i], &front[i], xip->rhc_V );
00801                 /* extrude front to create back plate */
00802                 VADD2( &back[i], &front[i], xip->rhc_H );
00803                 i += 3;
00804                 old = pos;
00805                 pos = pos->next;
00806                 bu_free ( (char *)old, "rt_pt_node" );
00807         }
00808 
00809         /* Draw the front */
00810         RT_ADD_VLIST( vhead, &front[(n-1)*ELEMENTS_PER_VECT],
00811                 BN_VLIST_LINE_MOVE );
00812         for( i = 0; i < n; i++ )  {
00813                 RT_ADD_VLIST( vhead, &front[i*ELEMENTS_PER_VECT], BN_VLIST_LINE_DRAW );
00814         }
00815 
00816         /* Draw the back */
00817         RT_ADD_VLIST( vhead, &back[(n-1)*ELEMENTS_PER_VECT], BN_VLIST_LINE_MOVE );
00818         for( i = 0; i < n; i++ )  {
00819                 RT_ADD_VLIST( vhead, &back[i*ELEMENTS_PER_VECT], BN_VLIST_LINE_DRAW );
00820         }
00821 
00822         /* Draw connections */
00823         for( i = 0; i < n; i++ )  {
00824                 RT_ADD_VLIST( vhead, &front[i*ELEMENTS_PER_VECT], BN_VLIST_LINE_MOVE );
00825                 RT_ADD_VLIST( vhead, &back[i*ELEMENTS_PER_VECT], BN_VLIST_LINE_DRAW );
00826         }
00827 
00828         /* free mem */
00829         bu_free( (char *)front, "fastf_t");
00830         bu_free( (char *)back, "fastf_t");
00831 
00832         return(0);
00833 }
00834 
00835 /**
00836  *      R T _ M K _ H Y P E R B O L A
00837  */
00838 int
00839 rt_mk_hyperbola(struct rt_pt_node *pts, fastf_t r, fastf_t b, fastf_t c, fastf_t dtol, fastf_t ntol)
00840 {
00841         fastf_t A, B, C, discr, dist, intr, j, k, m, theta0, theta1, z0;
00842         int     n;
00843         point_t mpt, p0, p1;
00844         vect_t  norm_line, norm_hyperb;
00845         struct rt_pt_node *new, *rt_ptalloc(void);
00846 
00847 #define RHC_TOL .0001
00848         /* endpoints of segment approximating hyperbola */
00849         VMOVE( p0, pts->p );
00850         VMOVE( p1, pts->next->p );
00851         /* slope and intercept of segment */
00852         m = ( p1[Z] - p0[Z] ) / ( p1[Y] - p0[Y] );
00853         intr = p0[Z] - m * p0[Y];
00854         /* find point on hyperbola with max dist between hyperbola and line */
00855         j = b + c;
00856         k = 1 - m*m*r*r/(b*(b + 2*c));
00857         A = k;
00858         B = 2*j*k;
00859         C = j*j*k - c*c;
00860         discr = sqrt(B*B - 4*A*C);
00861         z0 = (-B + discr) / (2. * A);
00862         if ( z0+RHC_TOL >= -b ) /* use top sheet of hyperboloid */
00863                 mpt[Z] = z0;
00864         else
00865                 mpt[Z] = (-B - discr) / (2. * A);
00866         if (NEAR_ZERO( mpt[Z], RHC_TOL))
00867                 mpt[Z] = 0.;
00868         mpt[X] = 0;
00869         mpt[Y] = ((mpt[Z] + b + c) * (mpt[Z] + b + c) - c*c) / (b*(b + 2*c));
00870         if (NEAR_ZERO( mpt[Y], RHC_TOL ))
00871                 mpt[Y] = 0.;
00872         mpt[Y] = r * sqrt( mpt[Y] );
00873         if (p0[Y] < 0.)
00874                 mpt[Y] = -mpt[Y];
00875         /* max distance between that point and line */
00876         dist = fabs( m * mpt[Y] - mpt[Z] + intr ) / sqrt( m * m + 1 );
00877         /* angles between normal of line and of hyperbola at line endpoints */
00878         VSET( norm_line, m, -1., 0.);
00879         VSET( norm_hyperb, 0., (2*c + 1) / (p0[Z] + c + 1), -1.);
00880         VUNITIZE( norm_line );
00881         VUNITIZE( norm_hyperb );
00882         theta0 = fabs( acos( VDOT( norm_line, norm_hyperb )));
00883         VSET( norm_hyperb, 0., (2*c + 1) / (p1[Z] + c + 1), -1.);
00884         VUNITIZE( norm_hyperb );
00885         theta1 = fabs( acos( VDOT( norm_line, norm_hyperb )));
00886         /* split segment at widest point if not within error tolerances */
00887         if ( dist > dtol || theta0 > ntol || theta1 > ntol ) {
00888                 /* split segment */
00889                 new = rt_ptalloc();
00890                 VMOVE( new->p, mpt );
00891                 new->next = pts->next;
00892                 pts->next = new;
00893                 /* keep track of number of pts added */
00894                 n = 1;
00895                 /* recurse on first new segment */
00896                 n += rt_mk_hyperbola( pts, r, b, c, dtol, ntol );
00897                 /* recurse on second new segment */
00898                 n += rt_mk_hyperbola( new, r, b, c, dtol, ntol );
00899         } else
00900                 n  = 0;
00901         return( n );
00902 }
00903 
00904 /**
00905  *                      R T _ R H C _ T E S S
00906  *
00907  *  Returns -
00908  *      -1      failure
00909  *       0      OK.  *r points to nmgregion that holds this tessellation.
00910  */
00911 int
00912 rt_rhc_tess(struct nmgregion **r, struct model *m, struct rt_db_internal *ip, const struct rt_tess_tol *ttol, const struct bn_tol *tol)
00913 {
00914         int             i, j, n;
00915         fastf_t         b, c, *back, f, *front, h, rh;
00916         fastf_t         dtol, ntol;
00917         vect_t          Bu, Hu, Ru;
00918         LOCAL mat_t     R;
00919         LOCAL mat_t     invR;
00920         LOCAL struct rt_rhc_internal    *xip;
00921         struct rt_pt_node       *old, *pos, *pts, *rt_ptalloc(void);
00922         struct shell    *s;
00923         struct faceuse  **outfaceuses;
00924         struct vertex   **vfront, **vback, **vtemp, *vertlist[4];
00925         vect_t          *norms;
00926         fastf_t         bb_plus_2bc,b_plus_c,r_sq;
00927         int             failure=0;
00928 
00929         NMG_CK_MODEL( m );
00930         BN_CK_TOL( tol );
00931         RT_CK_TESS_TOL( ttol );
00932 
00933         RT_CK_DB_INTERNAL(ip);
00934         xip = (struct rt_rhc_internal *)ip->idb_ptr;
00935         RT_RHC_CK_MAGIC(xip);
00936 
00937         /* compute |B| |H| */
00938         b = MAGNITUDE( xip->rhc_B );    /* breadth */
00939         rh = xip->rhc_r;                /* rectangular halfwidth */
00940         h = MAGNITUDE( xip->rhc_H );    /* height */
00941         c = xip->rhc_c;                 /* dist to asympt origin */
00942 
00943         /* Check for |H| > 0, |B| > 0, rh > 0, c > 0 */
00944         if( NEAR_ZERO(h, RT_LEN_TOL) || NEAR_ZERO(b, RT_LEN_TOL)
00945          || NEAR_ZERO(rh, RT_LEN_TOL) || NEAR_ZERO(c, RT_LEN_TOL))  {
00946                 bu_log("rt_rhc_tess:  zero length H, B, c, or rh\n");
00947                 return(-2);             /* BAD */
00948         }
00949 
00950         /* Check for B.H == 0 */
00951         f = VDOT( xip->rhc_B, xip->rhc_H ) / (b * h);
00952         if( ! NEAR_ZERO(f, RT_DOT_TOL) )  {
00953                 bu_log("rt_rhc_tess: B not perpendicular to H, f=%f\n", f);
00954                 return(-3);             /* BAD */
00955         }
00956 
00957         /* make unit vectors in B, H, and BxH directions */
00958         VMOVE(    Hu, xip->rhc_H );
00959         VUNITIZE( Hu );
00960         VMOVE(    Bu, xip->rhc_B );
00961         VUNITIZE( Bu );
00962         VCROSS(   Ru, Bu, Hu );
00963 
00964         /* Compute R and Rinv matrices */
00965         MAT_IDN( R );
00966         VREVERSE( &R[0], Hu );
00967         VMOVE(    &R[4], Ru );
00968         VREVERSE( &R[8], Bu );
00969         bn_mat_trn( invR, R );                  /* inv of rot mat is trn */
00970 
00971         /*
00972          *  Establish tolerances
00973          */
00974         if( ttol->rel <= 0.0 || ttol->rel >= 1.0 )  {
00975                 dtol = 0.0;             /* none */
00976         } else {
00977                 /* Convert rel to absolute by scaling by smallest side */
00978                 if (rh < b)
00979                         dtol = ttol->rel * 2 * rh;
00980                 else
00981                         dtol = ttol->rel * 2 * b;
00982         }
00983         if( ttol->abs <= 0.0 )  {
00984                 if( dtol <= 0.0 )  {
00985                         /* No tolerance given, use a default */
00986                         if (rh < b)
00987                                 dtol = 2 * 0.10 * rh;   /* 10% */
00988                         else
00989                                 dtol = 2 * 0.10 * b;    /* 10% */
00990                 } else {
00991                         /* Use absolute-ized relative tolerance */
00992                 }
00993         } else {
00994                 /* Absolute tolerance was given, pick smaller */
00995                 if( ttol->rel <= 0.0 || dtol > ttol->abs )
00996                         dtol = ttol->abs;
00997         }
00998 
00999         /* To ensure normal tolerance, remain below this angle */
01000         if( ttol->norm > 0.0 )
01001                 ntol = ttol->norm;
01002         else
01003                 /* tolerate everything */
01004                 ntol = bn_pi;
01005 
01006         /* initial hyperbola approximation is a single segment */
01007         pts = rt_ptalloc();
01008         pts->next = rt_ptalloc();
01009         pts->next->next = NULL;
01010         VSET( pts->p,       0, -rh, 0);
01011         VSET( pts->next->p, 0,  rh, 0);
01012         /* 2 endpoints in 1st approximation */
01013         n = 2;
01014         /* recursively break segment 'til within error tolerances */
01015         n += rt_mk_hyperbola( pts, rh, b, c, dtol, ntol );
01016 
01017         /* get mem for arrays */
01018         front = (fastf_t *)bu_malloc(3*n * sizeof(fastf_t), "fastf_t");
01019         back  = (fastf_t *)bu_malloc(3*n * sizeof(fastf_t), "fastf_t");
01020         norms = (vect_t *)bu_calloc( n , sizeof( vect_t ) , "rt_rhc_tess: norms" );
01021         vfront = (struct vertex **)bu_malloc((n+1) * sizeof(struct vertex *), "vertex *");
01022         vback = (struct vertex **)bu_malloc((n+1) * sizeof(struct vertex *), "vertex *");
01023         vtemp = (struct vertex **)bu_malloc((n+1) * sizeof(struct vertex *), "vertex *");
01024         outfaceuses =
01025                 (struct faceuse **)bu_malloc((n+2) * sizeof(struct faceuse *), "faceuse *");
01026         if (!front || !back || !vfront || !vback || !vtemp || !outfaceuses) {
01027                 fprintf(stderr, "rt_rhc_tess: no memory!\n");
01028                 goto fail;
01029         }
01030 
01031         /* generate front & back plates in world coordinates */
01032         bb_plus_2bc = b*b + 2.0*b*c;
01033         b_plus_c = b + c;
01034         r_sq = rh*rh;
01035         pos = pts;
01036         i = 0;
01037         j = 0;
01038         while (pos) {
01039                 vect_t tmp_norm;
01040 
01041                 /* calculate normal for 2D hyperbola */
01042                 VSET( tmp_norm , 0.0 , pos->p[Y]*bb_plus_2bc , (-r_sq*(pos->p[Z]+b_plus_c)) );
01043                 MAT4X3VEC( norms[j] , invR , tmp_norm );
01044                 VUNITIZE( norms[j] );
01045                 /* rotate back to original position */
01046                 MAT4X3VEC( &front[i], invR, pos->p );
01047                 /* move to origin vertex origin */
01048                 VADD2( &front[i], &front[i], xip->rhc_V );
01049                 /* extrude front to create back plate */
01050                 VADD2( &back[i], &front[i], xip->rhc_H );
01051                 i += 3;
01052                 j++;
01053                 old = pos;
01054                 pos = pos->next;
01055                 bu_free ( (char *)old, "rt_pt_node" );
01056         }
01057 
01058         *r = nmg_mrsv( m );     /* Make region, empty shell, vertex */
01059         s = BU_LIST_FIRST(shell, &(*r)->s_hd);
01060 
01061         for( i=0; i<n; i++ )  {
01062                 vfront[i] = vtemp[i] = (struct vertex *)0;
01063         }
01064 
01065         /* Front face topology.  Verts are considered to go CCW */
01066         outfaceuses[0] = nmg_cface(s, vfront, n);
01067 
01068         (void)nmg_mark_edges_real( &outfaceuses[0]->l.magic );
01069 
01070         /* Back face topology.  Verts must go in opposite dir (CW) */
01071         outfaceuses[1] = nmg_cface(s, vtemp, n);
01072         for( i=0; i<n; i++ )  vback[i] = vtemp[n-1-i];
01073 
01074         (void)nmg_mark_edges_real( &outfaceuses[1]->l.magic );
01075 
01076         /* Duplicate [0] as [n] to handle loop end condition, below */
01077         vfront[n] = vfront[0];
01078         vback[n] = vback[0];
01079 
01080         /* Build topology for all the rectangular side faces (n of them)
01081          * connecting the front and back faces.
01082          * increasing indices go towards counter-clockwise (CCW).
01083          */
01084         for( i=0; i<n; i++ )  {
01085                 vertlist[0] = vfront[i];        /* from top, */
01086                 vertlist[1] = vback[i];         /* straight down, */
01087                 vertlist[2] = vback[i+1];       /* to left, */
01088                 vertlist[3] = vfront[i+1];      /* straight up. */
01089                 outfaceuses[2+i] = nmg_cface(s, vertlist, 4);
01090         }
01091 
01092         (void)nmg_mark_edges_real( &outfaceuses[n+1]->l.magic );
01093 
01094         for( i=0; i<n; i++ )  {
01095                 NMG_CK_VERTEX(vfront[i]);
01096                 NMG_CK_VERTEX(vback[i]);
01097         }
01098 
01099         /* Associate the vertex geometry, CCW */
01100         for( i=0; i<n; i++ )  {
01101                 nmg_vertex_gv( vfront[i], &front[3*(i)] );
01102         }
01103         for( i=0; i<n; i++ )  {
01104                 nmg_vertex_gv( vback[i], &back[3*(i)] );
01105         }
01106 
01107         /* Associate the face geometry */
01108         for (i=0 ; i < n+2 ; i++) {
01109                 if( nmg_fu_planeeqn( outfaceuses[i], tol ) < 0 )
01110                 {
01111                         failure = (-1);
01112                         goto fail;
01113                 }
01114         }
01115 
01116         /* Associate vertexuse normals */
01117         for( i=0 ; i<n ; i++ )
01118         {
01119                 struct vertexuse *vu;
01120                 struct faceuse *fu;
01121                 vect_t rev_norm;
01122 
01123                 VREVERSE( rev_norm , norms[i] );
01124 
01125                 /* do "front" vertices */
01126                 NMG_CK_VERTEX( vfront[i] );
01127                 for( BU_LIST_FOR( vu , vertexuse , &vfront[i]->vu_hd ) )
01128                 {
01129                         NMG_CK_VERTEXUSE( vu );
01130                         fu = nmg_find_fu_of_vu( vu );
01131                         NMG_CK_FACEUSE( fu );
01132                         if( fu->f_p == outfaceuses[0]->f_p ||
01133                             fu->f_p == outfaceuses[1]->f_p ||
01134                             fu->f_p == outfaceuses[n+1]->f_p )
01135                                         continue;       /* skip flat faces */
01136 
01137                         if( fu->orientation == OT_SAME )
01138                                 nmg_vertexuse_nv( vu , norms[i] );
01139                         else if( fu->orientation == OT_OPPOSITE )
01140                                 nmg_vertexuse_nv( vu , rev_norm );
01141                 }
01142 
01143                 /* and "back" vertices */
01144                 NMG_CK_VERTEX( vback[i] );
01145                 for( BU_LIST_FOR( vu , vertexuse , &vback[i]->vu_hd ) )
01146                 {
01147                         NMG_CK_VERTEXUSE( vu );
01148                         fu = nmg_find_fu_of_vu( vu );
01149                         NMG_CK_FACEUSE( fu );
01150                         if( fu->f_p == outfaceuses[0]->f_p ||
01151                             fu->f_p == outfaceuses[1]->f_p ||
01152                             fu->f_p == outfaceuses[n+1]->f_p )
01153                                         continue;       /* skip flat faces */
01154 
01155                         if( fu->orientation == OT_SAME )
01156                                 nmg_vertexuse_nv( vu , norms[i] );
01157                         else if( fu->orientation == OT_OPPOSITE )
01158                                 nmg_vertexuse_nv( vu , rev_norm );
01159                 }
01160         }
01161 
01162         /* Glue the edges of different outward pointing face uses together */
01163         nmg_gluefaces( outfaceuses, n+2, tol );
01164 
01165         /* Compute "geometry" for region and shell */
01166         nmg_region_a( *r, tol );
01167 
01168 fail:
01169         /* free mem */
01170         bu_free( (char *)front, "fastf_t");
01171         bu_free( (char *)back, "fastf_t");
01172         bu_free( (char*)vfront, "vertex *");
01173         bu_free( (char*)vback, "vertex *");
01174         bu_free( (char*)vtemp, "vertex *");
01175         bu_free( (char *)norms , "rt_rhc_tess: norms" );
01176         bu_free( (char*)outfaceuses, "faceuse *");
01177 
01178         return( failure );
01179 }
01180 
01181 /**
01182  *                      R T _ R H C _ I M P O R T
01183  *
01184  *  Import an RHC from the database format to the internal format.
01185  *  Apply modeling transformations as well.
01186  */
01187 int
01188 rt_rhc_import(struct rt_db_internal *ip, const struct bu_external *ep, register const fastf_t *mat, const struct db_i *dbip)
01189 {
01190         LOCAL struct rt_rhc_internal    *xip;
01191         union record                    *rp;
01192 
01193         BU_CK_EXTERNAL( ep );
01194         rp = (union record *)ep->ext_buf;
01195         /* Check record type */
01196         if( rp->u_id != ID_SOLID )  {
01197                 bu_log("rt_rhc_import: defective record\n");
01198                 return(-1);
01199         }
01200 
01201         RT_CK_DB_INTERNAL( ip );
01202         ip->idb_major_type = DB5_MAJORTYPE_BRLCAD;
01203         ip->idb_type = ID_RHC;
01204         ip->idb_meth = &rt_functab[ID_RHC];
01205         ip->idb_ptr = bu_malloc( sizeof(struct rt_rhc_internal), "rt_rhc_internal");
01206         xip = (struct rt_rhc_internal *)ip->idb_ptr;
01207         xip->rhc_magic = RT_RHC_INTERNAL_MAGIC;
01208 
01209         /* Warning:  type conversion */
01210         MAT4X3PNT( xip->rhc_V, mat, &rp->s.s_values[0*3] );
01211         MAT4X3VEC( xip->rhc_H, mat, &rp->s.s_values[1*3] );
01212         MAT4X3VEC( xip->rhc_B, mat, &rp->s.s_values[2*3] );
01213         xip->rhc_r = rp->s.s_values[3*3] / mat[15];
01214         xip->rhc_c = rp->s.s_values[3*3+1] / mat[15];
01215 
01216         if( xip->rhc_r < SMALL_FASTF || xip->rhc_c < SMALL_FASTF )
01217         {
01218                 bu_log( "rt_rhc_import: r or c are zero\n" );
01219                 bu_free( (char *)ip->idb_ptr , "rt_rhc_import: ip->idb_ptr" );
01220                 return( -1 );
01221         }
01222 
01223         return(0);                      /* OK */
01224 }
01225 
01226 /**
01227  *                      R T _ R H C _ E X P O R T
01228  *
01229  *  The name is added by the caller, in the usual place.
01230  */
01231 int
01232 rt_rhc_export(struct bu_external *ep, const struct rt_db_internal *ip, double local2mm, const struct db_i *dbip)
01233 {
01234         struct rt_rhc_internal  *xip;
01235         union record            *rhc;
01236 
01237         RT_CK_DB_INTERNAL(ip);
01238         if( ip->idb_type != ID_RHC )  return(-1);
01239         xip = (struct rt_rhc_internal *)ip->idb_ptr;
01240         RT_RHC_CK_MAGIC(xip);
01241 
01242         BU_CK_EXTERNAL(ep);
01243         ep->ext_nbytes = sizeof(union record);
01244         ep->ext_buf = (genptr_t)bu_calloc( 1, ep->ext_nbytes, "rhc external");
01245         rhc = (union record *)ep->ext_buf;
01246 
01247         rhc->s.s_id = ID_SOLID;
01248         rhc->s.s_type = RHC;
01249 
01250         if (MAGNITUDE(xip->rhc_B) < RT_LEN_TOL
01251                 || MAGNITUDE(xip->rhc_H) < RT_LEN_TOL
01252                 || xip->rhc_r < RT_LEN_TOL
01253                 || xip->rhc_c < RT_LEN_TOL) {
01254                 bu_log("rt_rhc_export: not all dimensions positive!\n");
01255                 return(-1);
01256         }
01257 
01258         {
01259                 vect_t ub, uh;
01260 
01261                 VMOVE(ub, xip->rhc_B);
01262                 VUNITIZE(ub);
01263                 VMOVE(uh, xip->rhc_H);
01264                 VUNITIZE(uh);
01265                 if ( !NEAR_ZERO( VDOT(ub, uh), RT_DOT_TOL) ) {
01266                         bu_log("rt_rhc_export: B and H are not perpendicular!\n");
01267                         return(-1);
01268                 }
01269         }
01270 
01271         /* Warning:  type conversion */
01272         VSCALE( &rhc->s.s_values[0*3], xip->rhc_V, local2mm );
01273         VSCALE( &rhc->s.s_values[1*3], xip->rhc_H, local2mm );
01274         VSCALE( &rhc->s.s_values[2*3], xip->rhc_B, local2mm );
01275         rhc->s.s_values[3*3] = xip->rhc_r * local2mm;
01276         rhc->s.s_values[3*3+1] = xip->rhc_c * local2mm;
01277 
01278         return(0);
01279 }
01280 
01281 /**
01282  *                      R T _ R H C _ I M P O R T 5
01283  *
01284  *  Import an RHC from the database format to the internal format.
01285  *  Apply modeling transformations as well.
01286  */
01287 int
01288 rt_rhc_import5(struct rt_db_internal *ip, const struct bu_external *ep, register const fastf_t *mat, const struct db_i *dbip)
01289 {
01290         LOCAL struct rt_rhc_internal    *xip;
01291         fastf_t                 vec[11];
01292 
01293         BU_CK_EXTERNAL( ep );
01294 
01295         BU_ASSERT_LONG( ep->ext_nbytes, ==, SIZEOF_NETWORK_DOUBLE * 11 );
01296 
01297         RT_CK_DB_INTERNAL( ip );
01298         ip->idb_major_type = DB5_MAJORTYPE_BRLCAD;
01299         ip->idb_type = ID_RHC;
01300         ip->idb_meth = &rt_functab[ID_RHC];
01301         ip->idb_ptr = bu_malloc( sizeof(struct rt_rhc_internal), "rt_rhc_internal");
01302 
01303         xip = (struct rt_rhc_internal *)ip->idb_ptr;
01304         xip->rhc_magic = RT_RHC_INTERNAL_MAGIC;
01305 
01306         /* Convert from database (network) to internal (host) format */
01307         ntohd( (unsigned char *)vec, ep->ext_buf, 11 );
01308 
01309         /* Apply modeling transformations */
01310         MAT4X3PNT( xip->rhc_V, mat, &vec[0*3] );
01311         MAT4X3VEC( xip->rhc_H, mat, &vec[1*3] );
01312         MAT4X3VEC( xip->rhc_B, mat, &vec[2*3] );
01313         xip->rhc_r = vec[3*3] / mat[15];
01314         xip->rhc_c = vec[3*3+1] / mat[15];
01315 
01316         if( xip->rhc_r < SMALL_FASTF || xip->rhc_c < SMALL_FASTF )
01317         {
01318                 bu_log( "rt_rhc_import: r or c are zero\n" );
01319                 bu_free( (char *)ip->idb_ptr , "rt_rhc_import: ip->idb_ptr" );
01320                 return( -1 );
01321         }
01322 
01323         return(0);                      /* OK */
01324 }
01325 
01326 /**
01327  *                      R T _ R H C _ E X P O R T 5
01328  *
01329  *  The name is added by the caller, in the usual place.
01330  */
01331 int
01332 rt_rhc_export5(struct bu_external *ep, const struct rt_db_internal *ip, double local2mm, const struct db_i *dbip)
01333 {
01334         struct rt_rhc_internal  *xip;
01335         fastf_t                 vec[11];
01336 
01337         RT_CK_DB_INTERNAL(ip);
01338         if( ip->idb_type != ID_RHC )  return(-1);
01339         xip = (struct rt_rhc_internal *)ip->idb_ptr;
01340         RT_RHC_CK_MAGIC(xip);
01341 
01342         BU_CK_EXTERNAL(ep);
01343         ep->ext_nbytes = SIZEOF_NETWORK_DOUBLE * 11;
01344         ep->ext_buf = (genptr_t)bu_malloc( ep->ext_nbytes, "rhc external");
01345 
01346         if (MAGNITUDE(xip->rhc_B) < RT_LEN_TOL
01347                 || MAGNITUDE(xip->rhc_H) < RT_LEN_TOL
01348                 || xip->rhc_r < RT_LEN_TOL
01349                 || xip->rhc_c < RT_LEN_TOL) {
01350                 bu_log("rt_rhc_export: not all dimensions positive!\n");
01351                 return(-1);
01352         }
01353 
01354         {
01355                 vect_t ub, uh;
01356 
01357                 VMOVE(ub, xip->rhc_B);
01358                 VUNITIZE(ub);
01359                 VMOVE(uh, xip->rhc_H);
01360                 VUNITIZE(uh);
01361                 if ( !NEAR_ZERO( VDOT(ub, uh), RT_DOT_TOL) ) {
01362                         bu_log("rt_rhc_export: B and H are not perpendicular!\n");
01363                         return(-1);
01364                 }
01365         }
01366 
01367         /* scale 'em into local buffer */
01368         VSCALE( &vec[0*3], xip->rhc_V, local2mm );
01369         VSCALE( &vec[1*3], xip->rhc_H, local2mm );
01370         VSCALE( &vec[2*3], xip->rhc_B, local2mm );
01371         vec[3*3] = xip->rhc_r * local2mm;
01372         vec[3*3+1] = xip->rhc_c * local2mm;
01373 
01374         /* Convert from internal (host) to database (network) format */
01375         htond( ep->ext_buf, (unsigned char *)vec, 11 );
01376 
01377         return(0);
01378 }
01379 
01380 /**
01381  *                      R T _ R H C _ D E S C R I B E
01382  *
01383  *  Make human-readable formatted presentation of this solid.
01384  *  First line describes type of solid.
01385  *  Additional lines are indented one tab, and give parameter values.
01386  */
01387 int
01388 rt_rhc_describe(struct bu_vls *str, const struct rt_db_internal *ip, int verbose, double mm2local)
01389 {
01390         register struct rt_rhc_internal *xip =
01391                 (struct rt_rhc_internal *)ip->idb_ptr;
01392         char    buf[256];
01393 
01394         RT_RHC_CK_MAGIC(xip);
01395         bu_vls_strcat( str, "Right Hyperbolic Cylinder (RHC)\n");
01396 
01397         sprintf(buf, "\tV (%g, %g, %g)\n",
01398                 INTCLAMP(xip->rhc_V[X] * mm2local),
01399                 INTCLAMP(xip->rhc_V[Y] * mm2local),
01400                 INTCLAMP(xip->rhc_V[Z] * mm2local) );
01401         bu_vls_strcat( str, buf );
01402 
01403         sprintf(buf, "\tB (%g, %g, %g) mag=%g\n",
01404                 INTCLAMP(xip->rhc_B[X] * mm2local),
01405                 INTCLAMP(xip->rhc_B[Y] * mm2local),
01406                 INTCLAMP(xip->rhc_B[Z] * mm2local),
01407                 INTCLAMP(MAGNITUDE(xip->rhc_B) * mm2local));
01408         bu_vls_strcat( str, buf );
01409 
01410         sprintf(buf, "\tH (%g, %g, %g) mag=%g\n",
01411                 INTCLAMP(xip->rhc_H[X] * mm2local),
01412                 INTCLAMP(xip->rhc_H[Y] * mm2local),
01413                 INTCLAMP(xip->rhc_H[Z] * mm2local),
01414                 INTCLAMP(MAGNITUDE(xip->rhc_H) * mm2local));
01415         bu_vls_strcat( str, buf );
01416 
01417         sprintf(buf, "\tr=%g\n", INTCLAMP(xip->rhc_r * mm2local));
01418         bu_vls_strcat( str, buf );
01419 
01420         sprintf(buf, "\tc=%g\n", INTCLAMP(xip->rhc_c * mm2local));
01421         bu_vls_strcat( str, buf );
01422 
01423         return(0);
01424 }
01425 
01426 /**
01427  *                      R T _ R H C _ I F R E E
01428  *
01429  *  Free the storage associated with the rt_db_internal version of this solid.
01430  */
01431 void
01432 rt_rhc_ifree(struct rt_db_internal *ip)
01433 {
01434         register struct rt_rhc_internal *xip;
01435 
01436         RT_CK_DB_INTERNAL(ip);
01437         xip = (struct rt_rhc_internal *)ip->idb_ptr;
01438         RT_RHC_CK_MAGIC(xip);
01439         xip->rhc_magic = 0;             /* sanity */
01440 
01441         bu_free( (char *)xip, "rhc ifree" );
01442         ip->idb_ptr = GENPTR_NULL;      /* sanity */
01443 }
01444 
01445 /*
01446  * Local Variables:
01447  * mode: C
01448  * tab-width: 8
01449  * c-basic-offset: 4
01450  * indent-tabs-mode: t
01451  * End:
01452  * ex: shiftwidth=4 tabstop=8
01453  */

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