g_ebm.c

Go to the documentation of this file.
00001 /*                         G _ E B M . C
00002  * BRL-CAD
00003  *
00004  * Copyright (c) 1988-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_ebm.c
00026  *      Intersect a ray with an Extruded Bitmap,
00027  *      where the bitmap is taken from a bw(5) file.
00028  *
00029  *  Author -
00030  *      Michael John Muuss
00031  *
00032  *  Source -
00033  *      SECAD/VLD Computing Consortium, Bldg 394
00034  *      The U. S. Army Ballistic Research Laboratory
00035  *      Aberdeen Proving Ground, Maryland  21005
00036  *
00037  */
00038 
00039 #ifndef lint
00040 static const char RCSebm[] = "@(#)$Header: /cvsroot/brlcad/brlcad/src/librt/g_ebm.c,v 14.16 2006/09/16 02:04:24 lbutler Exp $ (BRL)";
00041 #endif
00042 
00043 #include "common.h"
00044 
00045 #include <stdlib.h>
00046 #include <stddef.h>
00047 #include <stdio.h>
00048 #include <ctype.h>
00049 #include <math.h>
00050 #ifdef HAVE_STRING_H
00051 #  include <string.h>
00052 #else
00053 #  include <strings.h>
00054 #endif
00055 #include <ctype.h>
00056 
00057 #include "tcl.h"
00058 #include "machine.h"
00059 #include "vmath.h"
00060 #include "db.h"
00061 #include "nmg.h"
00062 #include "rtgeom.h"
00063 #include "raytrace.h"
00064 #include "./debug.h"
00065 #include "./fixpt.h"
00066 
00067 /*
00068 NOTES:
00069         Changed small to small11 for win32 compatibility
00070 */
00071 
00072 
00073 struct rt_ebm_specific {
00074         struct rt_ebm_internal  ebm_i;
00075         vect_t          ebm_xnorm;      /* local +X norm in model coords */
00076         vect_t          ebm_ynorm;
00077         vect_t          ebm_znorm;
00078         vect_t          ebm_cellsize;   /* ideal coords: size of each cell */
00079         vect_t          ebm_origin;     /* local coords of grid origin (0,0,0) for now */
00080         vect_t          ebm_large;      /* local coords of XYZ max */
00081         mat_t           ebm_mat;        /* model to ideal space */
00082 };
00083 #define RT_EBM_NULL     ((struct rt_ebm_specific *)0)
00084 
00085 #define RT_EBM_O(m)     bu_offsetof(struct rt_ebm_internal, m)
00086 
00087 const struct bu_structparse rt_ebm_parse[] = {
00088 #if CRAY && !__STDC__
00089         {"%s",  RT_EBM_NAME_LEN, "file",        1,      BU_STRUCTPARSE_FUNC_NULL },
00090 #else
00091         {"%s",  RT_EBM_NAME_LEN, "file", bu_offsetofarray(struct rt_ebm_internal, file), BU_STRUCTPARSE_FUNC_NULL },
00092 #endif
00093         {"%d",  1, "w",         RT_EBM_O(xdim),         BU_STRUCTPARSE_FUNC_NULL },
00094         {"%d",  1, "n",         RT_EBM_O(ydim),         BU_STRUCTPARSE_FUNC_NULL },
00095         {"%f",  1, "d",         RT_EBM_O(tallness),     BU_STRUCTPARSE_FUNC_NULL },
00096         {"%f",  16, "mat", bu_offsetofarray(struct rt_ebm_internal, mat), BU_STRUCTPARSE_FUNC_NULL },
00097         {"",    0, (char *)0, 0,                        BU_STRUCTPARSE_FUNC_NULL }
00098 };
00099 
00100 struct ebm_hit_private {
00101         int     x_cell;
00102         int     y_cell;
00103 };
00104 
00105 
00106 BU_EXTERN(int rt_ebm_dda,(struct xray *rp, struct soltab *stp,
00107         struct application *ap, struct seg *seghead));
00108 BU_EXTERN(int rt_seg_planeclip,(struct seg *out_hd, struct seg *in_hd,
00109         vect_t out_norm, fastf_t in, fastf_t out,
00110         struct xray *rp, struct application *ap));
00111 BU_EXTERN( void rt_ebm_plate, ( int x1, int y1, int x2, int y2,
00112         double t, mat_t mat, struct bu_list *vhead ) );
00113 
00114 /*
00115  *  Codes to represent surface normals.
00116  *  In a bitmap, there are only 4 possible normals.
00117  *  With this code, reverse the sign to reverse the direction.
00118  *  As always, the normal is expected to point outwards.
00119  */
00120 #define NORM_ZPOS       3
00121 #define NORM_YPOS       2
00122 #define NORM_XPOS       1
00123 #define NORM_XNEG       (-1)
00124 #define NORM_YNEG       (-2)
00125 #define NORM_ZNEG       (-3)
00126 
00127 /*
00128  *  Regular bit addressing is used:  (0..W-1, 0..N-1),
00129  *  but the bitmap is stored with two cells of zeros all around,
00130  *  so permissible subscripts run (-2..W+1, -2..N+1).
00131  *  This eliminates special-case code for the boundary conditions.
00132  */
00133 #define BIT_XWIDEN      2
00134 #define BIT_YWIDEN      2
00135 #define BIT(_eip,_xx,_yy)       \
00136         ((unsigned char *)((_eip)->mp->apbuf))[ \
00137                 ((_yy)+BIT_YWIDEN)*((_eip)->xdim + \
00138                 BIT_XWIDEN*2)+(_xx)+BIT_XWIDEN ]
00139 
00140 /**
00141  *                      R T _ S E G _ P L A N E C L I P
00142  *
00143  *  Take a segment chain, in sorted order (ascending hit_dist),
00144  *  and clip to the range (in, out) along the normal "out_norm".
00145  *  For the particular ray "rp", find the parametric distances:
00146  *      kmin is the minimum permissible parameter, "in" units away
00147  *      kmax is the maximum permissible parameter, "out" units away
00148  *
00149  *  Returns -
00150  *      1       OK: trimmed segment chain, still in sorted order
00151  *      0       ERROR
00152  */
00153 int
00154 rt_seg_planeclip(struct seg *out_hd, struct seg *in_hd, fastf_t *out_norm, fastf_t in, fastf_t out, struct xray *rp, struct application *ap)
00155 {
00156         fastf_t         norm_dist_min, norm_dist_max;
00157         fastf_t         slant_factor;
00158         fastf_t         kmin, kmax;
00159         vect_t          in_norm;
00160         register struct seg     *curr;
00161         int             out_norm_code;
00162         int             count;
00163 
00164         norm_dist_min = in - VDOT( rp->r_pt, out_norm );
00165         slant_factor = VDOT( rp->r_dir, out_norm );     /* always abs < 1 */
00166         if( NEAR_ZERO( slant_factor, SQRT_SMALL_FASTF ) )  {
00167                 if( norm_dist_min < 0.0 )  {
00168                         bu_log("rt_seg_planeclip ERROR -- ray parallel to baseplane, outside \n");
00169                         /* XXX Free segp chain */
00170                         return(0);
00171                 }
00172                 kmin = -INFINITY;
00173         } else
00174                 kmin =  norm_dist_min / slant_factor;
00175 
00176         VREVERSE( in_norm, out_norm );
00177         norm_dist_max = out - VDOT( rp->r_pt, in_norm );
00178         slant_factor = VDOT( rp->r_dir, in_norm );      /* always abs < 1 */
00179         if( NEAR_ZERO( slant_factor, SQRT_SMALL_FASTF ) )  {
00180                 if( norm_dist_max < 0.0 )  {
00181                         bu_log("rt_seg_planeclip ERROR -- ray parallel to baseplane, outside \n");
00182                         /* XXX Free segp chain */
00183                         return(0);
00184                 }
00185                 kmax = INFINITY;
00186         } else
00187                 kmax =  norm_dist_max / slant_factor;
00188 
00189         if( kmin > kmax )  {
00190                 /* If r_dir[Z] < 0, will need to swap min & max */
00191                 slant_factor = kmax;
00192                 kmax = kmin;
00193                 kmin = slant_factor;
00194                 out_norm_code = NORM_ZPOS;
00195         } else {
00196                 out_norm_code = NORM_ZNEG;
00197         }
00198         if(RT_G_DEBUG&DEBUG_EBM)bu_log("kmin=%g, kmax=%g, out_norm_code=%d\n", kmin, kmax, out_norm_code );
00199 
00200         count = 0;
00201         while( BU_LIST_WHILE( curr, seg, &(in_hd->l) ) )  {
00202                 BU_LIST_DEQUEUE( &(curr->l) );
00203                 if(RT_G_DEBUG&DEBUG_EBM)bu_log(" rt_seg_planeclip seg( %g, %g )\n", curr->seg_in.hit_dist, curr->seg_out.hit_dist );
00204                 if( curr->seg_out.hit_dist <= kmin )  {
00205                         if(RT_G_DEBUG&DEBUG_EBM)bu_log("seg_out %g <= kmin %g, freeing\n", curr->seg_out.hit_dist, kmin );
00206                         RT_FREE_SEG(curr, ap->a_resource);
00207                         continue;
00208                 }
00209                 if( curr->seg_in.hit_dist >= kmax )  {
00210                         if(RT_G_DEBUG&DEBUG_EBM)bu_log("seg_in  %g >= kmax %g, freeing\n", curr->seg_in.hit_dist, kmax );
00211                         RT_FREE_SEG(curr, ap->a_resource);
00212                         continue;
00213                 }
00214                 if( curr->seg_in.hit_dist <= kmin )  {
00215                         if(RT_G_DEBUG&DEBUG_EBM)bu_log("seg_in = kmin %g\n", kmin );
00216                         curr->seg_in.hit_dist = kmin;
00217                         curr->seg_in.hit_surfno = out_norm_code;
00218                 }
00219                 if( curr->seg_out.hit_dist >= kmax )  {
00220                         if(RT_G_DEBUG&DEBUG_EBM)bu_log("seg_out= kmax %g\n", kmax );
00221                         curr->seg_out.hit_dist = kmax;
00222                         curr->seg_out.hit_surfno = (-out_norm_code);
00223                 }
00224                 BU_LIST_INSERT( &(out_hd->l), &(curr->l) );
00225                 count += 2;
00226         }
00227         return( count );
00228 }
00229 
00230 static int rt_ebm_normtab[3] = { NORM_XPOS, NORM_YPOS, NORM_ZPOS };
00231 
00232 
00233 /**
00234  *                      R T _ E B M _ D D A
00235  *
00236  *  Step through the 2-D array, in local coordinates ("ideal space").
00237  */
00238 int
00239 rt_ebm_dda(register struct xray *rp, struct soltab *stp, struct application *ap, struct seg *seghead)
00240 {
00241         register struct rt_ebm_specific *ebmp =
00242                 (struct rt_ebm_specific *)stp->st_specific;
00243         vect_t  invdir;
00244         double  t0;     /* in point of cell */
00245         double  t1;     /* out point of cell */
00246         double  tmax;   /* out point of entire grid */
00247         vect_t  t;      /* next t value for XYZ cell plane intersect */
00248         vect_t  delta;  /* spacing of XYZ cell planes along ray */
00249         int     igrid[3];/* Grid cell coordinates of cell (integerized) */
00250         vect_t  P;      /* hit point */
00251         int     inside; /* inside/outside a solid flag */
00252         int     in_index;
00253         int     out_index;
00254         int     j;
00255 
00256         /* Compute the inverse of the direction cosines */
00257         if( !NEAR_ZERO( rp->r_dir[X], SQRT_SMALL_FASTF ) )  {
00258                 invdir[X]=1.0/rp->r_dir[X];
00259         } else {
00260                 invdir[X] = INFINITY;
00261                 rp->r_dir[X] = 0.0;
00262         }
00263         if( !NEAR_ZERO( rp->r_dir[Y], SQRT_SMALL_FASTF ) )  {
00264                 invdir[Y]=1.0/rp->r_dir[Y];
00265         } else {
00266                 invdir[Y] = INFINITY;
00267                 rp->r_dir[Y] = 0.0;
00268         }
00269         if( !NEAR_ZERO( rp->r_dir[Z], SQRT_SMALL_FASTF ) )  {
00270                 invdir[Z]=1.0/rp->r_dir[Z];
00271         } else {
00272                 invdir[Z] = INFINITY;
00273                 rp->r_dir[Z] = 0.0;
00274         }
00275 
00276         /* intersect ray with ideal grid rpp */
00277         VSETALL( P, 0 );
00278         if( ! rt_in_rpp(rp, invdir, P, ebmp->ebm_large ) )
00279                 return(0);      /* MISS */
00280 
00281         VJOIN1( P, rp->r_pt, rp->r_min, rp->r_dir );
00282         /* P is hit point (on RPP?) */
00283 
00284 if(RT_G_DEBUG&DEBUG_EBM)VPRINT("ebm_origin", ebmp->ebm_origin);
00285 if(RT_G_DEBUG&DEBUG_EBM)VPRINT("r_pt", rp->r_pt);
00286 if(RT_G_DEBUG&DEBUG_EBM)VPRINT("P", P);
00287 if(RT_G_DEBUG&DEBUG_EBM)VPRINT("cellsize", ebmp->ebm_cellsize);
00288         t0 = rp->r_min;
00289         tmax = rp->r_max;
00290 if(RT_G_DEBUG&DEBUG_EBM)bu_log("[shoot: r_min=%g, r_max=%g]\n", rp->r_min, rp->r_max);
00291 
00292         /* find grid cell where ray first hits ideal space bounding RPP */
00293         igrid[X] = (P[X] - ebmp->ebm_origin[X]) / ebmp->ebm_cellsize[X];
00294         igrid[Y] = (P[Y] - ebmp->ebm_origin[Y]) / ebmp->ebm_cellsize[Y];
00295         if( igrid[X] < 0 )  {
00296                 igrid[X] = 0;
00297         } else if( igrid[X] >= ebmp->ebm_i.xdim ) {
00298                 igrid[X] = ebmp->ebm_i.xdim-1;
00299         }
00300         if( igrid[Y] < 0 )  {
00301                 igrid[Y] = 0;
00302         } else if( igrid[Y] >= ebmp->ebm_i.ydim ) {
00303                 igrid[Y] = ebmp->ebm_i.ydim-1;
00304         }
00305 if(RT_G_DEBUG&DEBUG_EBM)bu_log("g[X] = %d, g[Y] = %d\n", igrid[X], igrid[Y]);
00306 
00307         if( rp->r_dir[X] == 0.0 && rp->r_dir[Y] == 0.0 )  {
00308                 register struct seg     *segp;
00309 
00310                 /*  Ray is traveling exactly along Z axis.
00311                  *  Just check the one cell hit.
00312                  *  Depend on higher level to clip ray to Z extent.
00313                  */
00314 if(RT_G_DEBUG&DEBUG_EBM)bu_log("ray on local Z axis\n");
00315                 if( BIT( &ebmp->ebm_i, igrid[X], igrid[Y] ) == 0 )
00316                         return(0);      /* MISS */
00317                 RT_GET_SEG(segp, ap->a_resource);
00318                 segp->seg_stp = stp;
00319                 segp->seg_in.hit_dist = 0;
00320                 segp->seg_out.hit_dist = INFINITY;
00321 
00322                 segp->seg_in.hit_vpriv[X] =
00323                         (double) igrid[X] / ebmp->ebm_i.xdim;
00324                 segp->seg_in.hit_vpriv[Y] =
00325                         (double) igrid[Y] / ebmp->ebm_i.ydim;
00326 
00327                 segp->seg_out.hit_vpriv[X] =
00328                         (double) igrid[X] / ebmp->ebm_i.xdim;
00329                 segp->seg_out.hit_vpriv[Y] =
00330                         (double) igrid[Y] / ebmp->ebm_i.ydim;
00331 
00332                 if( rp->r_dir[Z] < 0 )  {
00333                         segp->seg_in.hit_surfno = NORM_ZPOS;
00334                         segp->seg_out.hit_surfno = NORM_ZNEG;
00335                 } else {
00336                         segp->seg_in.hit_surfno = NORM_ZNEG;
00337                         segp->seg_out.hit_surfno = NORM_ZPOS;
00338                 }
00339                 BU_LIST_INSERT( &(seghead->l), &(segp->l) );
00340                 return(2);                      /* HIT */
00341         }
00342 
00343         /* X setup */
00344         if( rp->r_dir[X] == 0.0 )  {
00345                 t[X] = INFINITY;
00346                 delta[X] = 0;
00347         } else {
00348                 j = igrid[X];
00349                 if( rp->r_dir[X] < 0 ) j++;
00350                 t[X] = (ebmp->ebm_origin[X] + j*ebmp->ebm_cellsize[X] -
00351                         rp->r_pt[X]) * invdir[X];
00352                 delta[X] = ebmp->ebm_cellsize[X] * fabs(invdir[X]);
00353         }
00354         /* Y setup */
00355         if( rp->r_dir[Y] == 0.0 )  {
00356                 t[Y] = INFINITY;
00357                 delta[Y] = 0;
00358         } else {
00359                 j = igrid[Y];
00360                 if( rp->r_dir[Y] < 0 ) j++;
00361                 t[Y] = (ebmp->ebm_origin[Y] + j*ebmp->ebm_cellsize[Y] -
00362                         rp->r_pt[Y]) * invdir[Y];
00363                 delta[Y] = ebmp->ebm_cellsize[Y] * fabs(invdir[Y]);
00364         }
00365 #if 0
00366         /* Z setup */
00367         if( rp->r_dir[Z] == 0.0 )  {
00368                 t[Z] = INFINITY;
00369         } else {
00370                 /* Consider igrid[Z] to be either 0 or 1 */
00371                 if( rp->r_dir[Z] < 0 )  {
00372                         t[Z] = (ebmp->ebm_origin[Z] + ebmp->ebm_cellsize[Z] -
00373                                 rp->r_pt[Z]) * invdir[Z];
00374                 } else {
00375                         t[Z] = (ebmp->ebm_origin[Z] - rp->r_pt[Z]) * invdir[Z];
00376                 }
00377         }
00378 #endif
00379 
00380         /* The delta[] elements *must* be positive, as t must increase */
00381 if(RT_G_DEBUG&DEBUG_EBM)bu_log("t[X] = %g, delta[X] = %g\n", t[X], delta[X] );
00382 if(RT_G_DEBUG&DEBUG_EBM)bu_log("t[Y] = %g, delta[Y] = %g\n", t[Y], delta[Y] );
00383 
00384         /* Find face of entry into first cell -- max initial t value */
00385         if( t[X] == INFINITY ) {
00386                 in_index = Y;
00387                 t0 = t[Y];
00388         }
00389         else if( t[Y] == INFINITY ) {
00390                 in_index = X;
00391                 t0 = t[X];
00392         }
00393         else if( t[X] >= t[Y]  )  {
00394                 in_index = X;
00395                 t0 = t[X];
00396         } else {
00397                 in_index = Y;
00398                 t0 = t[Y];
00399         }
00400 if(RT_G_DEBUG&DEBUG_EBM)bu_log("Entry index is %s, t0=%g\n", in_index==X ? "X" : "Y", t0);
00401 
00402         /* Advance to next exits */
00403         t[X] += delta[X];
00404         t[Y] += delta[Y];
00405 
00406         /* Ensure that next exit is after first entrance */
00407         if( t[X] < t0 )  {
00408                 bu_log("*** advancing t[X]\n");
00409                 t[X] += delta[X];
00410         }
00411         if( t[Y] < t0 )  {
00412                 bu_log("*** advancing t[Y]\n");
00413                 t[Y] += delta[Y];
00414         }
00415 if(RT_G_DEBUG&DEBUG_EBM)bu_log("Exit t[X]=%g, t[Y]=%g\n", t[X], t[Y] );
00416 
00417         inside = 0;
00418 
00419         while( t0 < tmax ) {
00420                 int     val;
00421                 struct seg      *segp;
00422 
00423                 /* find minimum exit t value */
00424                 out_index = t[X] < t[Y] ? X : Y;
00425 
00426                 t1 = t[out_index];
00427 
00428                 /* Ray passes through cell igrid[XY] from t0 to t1 */
00429                 val = BIT( &ebmp->ebm_i, igrid[X], igrid[Y] );
00430 if(RT_G_DEBUG&DEBUG_EBM)bu_log("igrid [%d %d] from %g to %g, val=%d\n",
00431                         igrid[X], igrid[Y],
00432                         t0, t1, val );
00433 if(RT_G_DEBUG&DEBUG_EBM)bu_log("Exit index is %s, t[X]=%g, t[Y]=%g\n",
00434                         out_index==X ? "X" : "Y", t[X], t[Y] );
00435 
00436 
00437                 if( t1 <= t0 )  bu_log("ERROR ebm t1=%g < t0=%g\n", t1, t0 );
00438                 if( !inside )  {
00439                         if( val > 0 )  {
00440                                 /* Handle the transition from vacuum to solid */
00441                                 /* Start of segment (entering a full voxel) */
00442                                 inside = 1;
00443 
00444                                 RT_GET_SEG(segp, ap->a_resource);
00445                                 segp->seg_stp = stp;
00446                                 segp->seg_in.hit_dist = t0;
00447 
00448                                 segp->seg_in.hit_vpriv[X] =
00449                                         (double) igrid[X] / ebmp->ebm_i.xdim;
00450                                 segp->seg_in.hit_vpriv[Y] =
00451                                         (double) igrid[Y] / ebmp->ebm_i.ydim;
00452 
00453                                 /* Compute entry normal */
00454                                 if( rp->r_dir[in_index] < 0 )  {
00455                                         /* Go left, entry norm goes right */
00456                                         segp->seg_in.hit_surfno =
00457                                                 rt_ebm_normtab[in_index];
00458                                 }  else  {
00459                                         /* go right, entry norm goes left */
00460                                         segp->seg_in.hit_surfno =
00461                                                 (-rt_ebm_normtab[in_index]);
00462                                 }
00463                                 BU_LIST_INSERT( &(seghead->l), &(segp->l) );
00464 
00465                                 if(RT_G_DEBUG&DEBUG_EBM) bu_log("START t=%g, surfno=%d\n",
00466                                         t0, segp->seg_in.hit_surfno);
00467                         } else {
00468                                 /* Do nothing, marching through void */
00469                         }
00470                 } else {
00471                         register struct seg     *tail;
00472                         if( val > 0 )  {
00473                                 /* Do nothing, marching through solid */
00474                         } else {
00475                                 /* End of segment (now in an empty voxel) */
00476                                 /* Handle transition from solid to vacuum */
00477                                 inside = 0;
00478 
00479                                 tail = BU_LIST_LAST( seg, &(seghead->l) );
00480                                 tail->seg_out.hit_dist = t0;
00481 
00482                                 tail->seg_out.hit_vpriv[X] =
00483                                         (double) igrid[X] / ebmp->ebm_i.xdim;
00484                                 tail->seg_out.hit_vpriv[Y] =
00485                                         (double) igrid[Y] / ebmp->ebm_i.ydim;
00486 
00487                                 /* Compute exit normal */
00488                                 if( rp->r_dir[in_index] < 0 )  {
00489                                         /* Go left, exit normal goes left */
00490                                         tail->seg_out.hit_surfno =
00491                                                 (-rt_ebm_normtab[in_index]);
00492                                 }  else  {
00493                                         /* go right, exit norm goes right */
00494                                         tail->seg_out.hit_surfno =
00495                                                 rt_ebm_normtab[in_index];
00496                                 }
00497                                 if(RT_G_DEBUG&DEBUG_EBM) bu_log("END t=%g, surfno=%d\n",
00498                                         t0, tail->seg_out.hit_surfno );
00499                         }
00500                 }
00501 
00502                 /* Take next step */
00503                 t0 = t1;
00504                 in_index = out_index;
00505                 t[out_index] += delta[out_index];
00506                 if( rp->r_dir[out_index] > 0 ) {
00507                         igrid[out_index]++;
00508                 } else {
00509                         igrid[out_index]--;
00510                 }
00511         }
00512 
00513         if( inside )  {
00514                 register struct seg     *tail;
00515 
00516                 /* Close off the final segment */
00517                 tail = BU_LIST_LAST( seg, &(seghead->l) );
00518                 tail->seg_out.hit_dist = tmax;
00519                 VSETALL(tail->seg_out.hit_vpriv, 0.0);
00520 
00521                 /* Compute exit normal.  Previous out_index is now in_index */
00522                 if( rp->r_dir[in_index] < 0 )  {
00523                         /* Go left, exit normal goes left */
00524                         tail->seg_out.hit_surfno = (-rt_ebm_normtab[in_index]);
00525                 }  else  {
00526                         /* go right, exit norm goes right */
00527                         tail->seg_out.hit_surfno = rt_ebm_normtab[in_index];
00528                 }
00529                 if(RT_G_DEBUG&DEBUG_EBM) bu_log("closed END t=%g, surfno=%d\n",
00530                         tmax, tail->seg_out.hit_surfno );
00531         }
00532 
00533         if( BU_LIST_IS_EMPTY( &(seghead->l) ) )
00534                 return(0);
00535         return(2);
00536 }
00537 
00538 /**
00539  *                      R T _ E B M _ I M P O R T
00540  *
00541  *  Read in the information from the string solid record.
00542  *  Then, as a service to the application, read in the bitmap
00543  *  and set up some of the associated internal variables.
00544  */
00545 int
00546 rt_ebm_import(struct rt_db_internal *ip, const struct bu_external *ep, const fastf_t *mat, const struct db_i *dbip)
00547 {
00548         union record    *rp;
00549         register struct rt_ebm_internal *eip;
00550         struct bu_vls   str;
00551         int             nbytes;
00552         mat_t           tmat;
00553         struct bu_mapped_file   *mp;
00554 
00555         BU_CK_EXTERNAL( ep );
00556         rp = (union record *)ep->ext_buf;
00557         if( rp->u_id != DBID_STRSOL )  {
00558                 bu_log("rt_ebm_import: defective strsol record\n");
00559                 return(-1);
00560         }
00561 
00562         RT_CK_DB_INTERNAL( ip );
00563         ip->idb_major_type = DB5_MAJORTYPE_BRLCAD;
00564         ip->idb_type = ID_EBM;
00565         ip->idb_meth = &rt_functab[ID_EBM];
00566         ip->idb_ptr = bu_calloc(1, sizeof(struct rt_ebm_internal), "rt_ebm_internal");
00567         eip = (struct rt_ebm_internal *)ip->idb_ptr;
00568         eip->magic = RT_EBM_INTERNAL_MAGIC;
00569 
00570         /* Provide default orientation info */
00571         MAT_IDN( eip->mat );
00572 
00573         bu_vls_init( &str );
00574         bu_vls_strcpy( &str, rp->ss.ss_args );
00575         if( bu_struct_parse( &str, rt_ebm_parse, (char *)eip ) < 0 )  {
00576                 bu_vls_free( &str );
00577                 bu_free( (char *)eip , "rt_ebm_import: eip" );
00578                 ip->idb_type = ID_NULL;
00579                 ip->idb_ptr = (genptr_t)NULL;
00580                 return -2;
00581         }
00582         bu_vls_free( &str );
00583 
00584         /* Check for reasonable values */
00585         if( eip->file[0] == '\0' || eip->xdim < 1 ||
00586             eip->ydim < 1 || eip->mat[15] <= 0.0 ||
00587             eip->tallness <= 0.0 )  {
00588                 bu_struct_print( "Unreasonable EBM parameters", rt_ebm_parse,
00589                         (char *)eip );
00590                 bu_free( (char *)eip , "rt_ebm_import: eip" );
00591                 ip->idb_type = ID_NULL;
00592                 ip->idb_ptr = (genptr_t)NULL;
00593                 return -1;
00594         }
00595 
00596         /* Apply any modeling transforms to get final matrix */
00597         bn_mat_mul( tmat, mat, eip->mat );
00598         MAT_COPY( eip->mat, tmat );
00599 
00600         /* Get bit map from .bw(5) file */
00601         if( !(mp = bu_open_mapped_file_with_path( dbip->dbi_filepath, eip->file, "ebm" )) )  {
00602                 bu_log("rt_ebm_import() unable to open '%s'\n", eip->file);
00603                 bu_free( (char *)eip , "rt_ebm_import: eip" );
00604 fail:
00605                 ip->idb_type = ID_NULL;
00606                 ip->idb_ptr = (genptr_t)NULL;
00607                 return -1;
00608         }
00609         eip->mp = mp;
00610         if( mp->buflen < eip->xdim*eip->ydim )  {
00611                 bu_log("rt_ebm_import() file '%s' is too short %d < %d\n",
00612                         eip->file, mp->buflen, eip->xdim*eip->ydim );
00613                 goto fail;
00614         }
00615 
00616         nbytes = (eip->xdim+BIT_XWIDEN*2)*(eip->ydim+BIT_YWIDEN*2);
00617 
00618         /* If first use of this file, prepare in-memory buffer */
00619         if( !mp->apbuf )  {
00620                 register int    y;
00621                 unsigned char   *cp;
00622 
00623                 /* Prevent a multi-processor race */
00624                 bu_semaphore_acquire(RT_SEM_MODEL);
00625                 if( mp->apbuf )  {
00626                         /* someone else beat us, nothing more to do */
00627                         bu_semaphore_release(RT_SEM_MODEL);
00628                         return 0;
00629                 }
00630                 mp->apbuf = (genptr_t)bu_calloc(
00631                         1, nbytes, "rt_ebm_import bitmap" );
00632                 mp->apbuflen = nbytes;
00633 
00634                 bu_semaphore_release(RT_SEM_MODEL);
00635 
00636                 /* Because of in-memory padding, read each scanline separately */
00637                 cp = (unsigned char *)mp->buf;
00638                 for( y=0; y < eip->ydim; y++ )  {
00639                         /* BIT() addresses into mp->apbuf */
00640                         bcopy( cp, &BIT( eip, 0, y), eip->xdim );
00641                         cp += eip->xdim;
00642                 }
00643         }
00644         return( 0 );
00645 }
00646 
00647 /**
00648  *                      R T _ E B M _ E X P O R T
00649  *
00650  *  The name will be added by the caller.
00651  */
00652 int
00653 rt_ebm_export(struct bu_external *ep, const struct rt_db_internal *ip, double local2mm, const struct db_i *dbip)
00654 {
00655         struct rt_ebm_internal  *eip;
00656         struct rt_ebm_internal  ebm;    /* scaled version */
00657         union record            *rec;
00658         struct bu_vls           str;
00659 
00660         RT_CK_DB_INTERNAL(ip);
00661         if( ip->idb_type != ID_EBM )  return(-1);
00662         eip = (struct rt_ebm_internal *)ip->idb_ptr;
00663         RT_EBM_CK_MAGIC(eip);
00664         ebm = *eip;                     /* struct copy */
00665 
00666         /* Apply scale factor */
00667         ebm.mat[15] /= local2mm;
00668 
00669         BU_CK_EXTERNAL(ep);
00670         ep->ext_nbytes = sizeof(union record)*DB_SS_NGRAN;
00671         ep->ext_buf = (genptr_t)bu_calloc( 1, ep->ext_nbytes, "ebm external");
00672         rec = (union record *)ep->ext_buf;
00673 
00674         bu_vls_init( &str );
00675         bu_vls_struct_print( &str, rt_ebm_parse, (char *)&ebm );
00676 
00677         rec->ss.ss_id = DBID_STRSOL;
00678         strncpy( rec->ss.ss_keyword, "ebm", NAMESIZE-1 );
00679         strncpy( rec->ss.ss_args, bu_vls_addr(&str), DB_SS_LEN-1 );
00680         bu_vls_free( &str );
00681 
00682         return(0);
00683 }
00684 
00685 
00686 /**
00687  *                      R T _ E B M _ I M P O R T 5
00688  *
00689  *  Read in the information from the string solid record.
00690  *  Then, as a service to the application, read in the bitmap
00691  *  and set up some of the associated internal variables.
00692  */
00693 int
00694 rt_ebm_import5(struct rt_db_internal *ip, const struct bu_external *ep, const fastf_t *mat, const struct db_i *dbip)
00695 {
00696         register struct rt_ebm_internal *eip;
00697         struct bu_vls   str;
00698         int             nbytes;
00699         mat_t           tmat;
00700         struct bu_mapped_file   *mp;
00701 
00702         BU_CK_EXTERNAL( ep );
00703         RT_CK_DB_INTERNAL( ip );
00704 
00705         ip->idb_major_type = DB5_MAJORTYPE_BRLCAD;
00706         ip->idb_type = ID_EBM;
00707         ip->idb_meth = &rt_functab[ID_EBM];
00708         ip->idb_ptr = bu_calloc(1, sizeof(struct rt_ebm_internal), "rt_ebm_internal");
00709         eip = (struct rt_ebm_internal *)ip->idb_ptr;
00710         eip->magic = RT_EBM_INTERNAL_MAGIC;
00711 
00712         /* Provide default orientation info */
00713         MAT_IDN( eip->mat );
00714 
00715         bu_vls_init( &str );
00716         bu_vls_strcpy( &str, ep->ext_buf );
00717         if( bu_struct_parse( &str, rt_ebm_parse, (char *)eip ) < 0 )  {
00718                 bu_vls_free( &str );
00719                 bu_free( (char *)eip , "rt_ebm_import: eip" );
00720                 ip->idb_type = ID_NULL;
00721                 ip->idb_ptr = (genptr_t)NULL;
00722                 return -2;
00723         }
00724         bu_vls_free( &str );
00725 
00726         /* Check for reasonable values */
00727         if( eip->file[0] == '\0' || eip->xdim < 1 ||
00728             eip->ydim < 1 || eip->mat[15] <= 0.0 ||
00729             eip->tallness <= 0.0 )  {
00730                 bu_struct_print( "Unreasonable EBM parameters", rt_ebm_parse,
00731                         (char *)eip );
00732                 bu_free( (char *)eip , "rt_ebm_import: eip" );
00733                 ip->idb_type = ID_NULL;
00734                 ip->idb_ptr = (genptr_t)NULL;
00735                 return -1;
00736         }
00737 
00738         /* Apply any modeling transforms to get final matrix */
00739         bn_mat_mul( tmat, mat, eip->mat );
00740         MAT_COPY( eip->mat, tmat );
00741 
00742         /* Get bit map from .bw(5) file */
00743         if( !(mp = bu_open_mapped_file_with_path( dbip->dbi_filepath, eip->file, "ebm" )) )  {
00744                 bu_log("rt_ebm_import() unable to open '%s'\n", eip->file);
00745                 bu_free( (char *)eip , "rt_ebm_import: eip" );
00746 fail:
00747                 ip->idb_type = ID_NULL;
00748                 ip->idb_ptr = (genptr_t)NULL;
00749                 return -1;
00750         }
00751         eip->mp = mp;
00752         if( mp->buflen < eip->xdim*eip->ydim )  {
00753                 bu_log("rt_ebm_import() file '%s' is too short %d < %d\n",
00754                         eip->file, mp->buflen, eip->xdim*eip->ydim );
00755                 goto fail;
00756         }
00757 
00758         nbytes = (eip->xdim+BIT_XWIDEN*2)*(eip->ydim+BIT_YWIDEN*2);
00759 
00760         /* If first use of this file, prepare in-memory buffer */
00761         if( !mp->apbuf )  {
00762                 register int    y;
00763                 unsigned char   *cp;
00764 
00765                 /* Prevent a multi-processor race */
00766                 bu_semaphore_acquire(RT_SEM_MODEL);
00767                 if( mp->apbuf )  {
00768                         /* someone else beat us, nothing more to do */
00769                         bu_semaphore_release(RT_SEM_MODEL);
00770                         return 0;
00771                 }
00772                 mp->apbuf = (genptr_t)bu_calloc(
00773                         1, nbytes, "rt_ebm_import bitmap" );
00774                 mp->apbuflen = nbytes;
00775 
00776                 bu_semaphore_release(RT_SEM_MODEL);
00777 
00778                 /* Because of in-memory padding, read each scanline separately */
00779                 cp = (unsigned char *)mp->buf;
00780                 for( y=0; y < eip->ydim; y++ )  {
00781                         /* BIT() addresses into mp->apbuf */
00782                         bcopy( cp, &BIT( eip, 0, y), eip->xdim );
00783                         cp += eip->xdim;
00784                 }
00785         }
00786         return( 0 );
00787 }
00788 
00789 /**
00790  *                      R T _ E B M _ E X P O R T 5
00791  *
00792  *  The name will be added by the caller.
00793  */
00794 int
00795 rt_ebm_export5(struct bu_external *ep, const struct rt_db_internal *ip, double local2mm, const struct db_i *dbip)
00796 {
00797         struct rt_ebm_internal  *eip;
00798         struct rt_ebm_internal  ebm;    /* scaled version */
00799         struct bu_vls           str;
00800 
00801         RT_CK_DB_INTERNAL(ip);
00802         if( ip->idb_type != ID_EBM )  return(-1);
00803         eip = (struct rt_ebm_internal *)ip->idb_ptr;
00804         RT_EBM_CK_MAGIC(eip);
00805         ebm = *eip;                     /* struct copy */
00806 
00807         /* Apply scale factor */
00808         ebm.mat[15] /= local2mm;
00809 
00810         BU_CK_EXTERNAL(ep);
00811 
00812         bu_vls_init( &str );
00813         bu_vls_struct_print( &str, rt_ebm_parse, (char *)&ebm );
00814 
00815         ep->ext_nbytes = bu_vls_strlen( &str ) + 1;
00816         ep->ext_buf = (genptr_t)bu_calloc( 1, ep->ext_nbytes, "ebm external");
00817 
00818         strcpy( ep->ext_buf , bu_vls_addr(&str) );
00819         bu_vls_free( &str );
00820 
00821         return(0);
00822 }
00823 
00824 /**
00825  *                      R T _ E B M _ D E S C R I B E
00826  *
00827  *  Make human-readable formatted presentation of this solid.
00828  *  First line describes type of solid.
00829  *  Additional lines are indented one tab, and give parameter values.
00830  */
00831 int
00832 rt_ebm_describe(struct bu_vls *str, const struct rt_db_internal *ip, int verbose, double mm2local)
00833 {
00834         register struct rt_ebm_internal *eip =
00835                 (struct rt_ebm_internal *)ip->idb_ptr;
00836         int i;
00837         struct bu_vls substr;
00838 
00839         RT_EBM_CK_MAGIC(eip);
00840 
00841         bu_vls_init( &substr );
00842         bu_vls_strcat( str, "extruded bitmap (EBM)\n\t");
00843 
00844 /*      bu_vls_struct_print( str, rt_ebm_parse, (char *)eip );
00845         bu_vls_strcat( str, "\n" );     */
00846 
00847         bu_vls_printf( &substr, "  file=\"%s\" w=%d n=%d depth=%g\n   mat=",
00848                 eip->file, eip->xdim, eip->ydim, INTCLAMP(eip->tallness*mm2local) );
00849         bu_vls_vlscat( str, &substr );
00850         for( i=0 ; i<15 ; i++ )
00851         {
00852                 bu_vls_trunc2( &substr, 0 );
00853                 bu_vls_printf( &substr, "%g,", INTCLAMP(eip->mat[i]) );
00854                 bu_vls_vlscat( str, &substr );
00855         }
00856         bu_vls_trunc2( &substr, 0 );
00857         bu_vls_printf( &substr, "%g\n", INTCLAMP(eip->mat[15]) );
00858         bu_vls_vlscat( str, &substr );
00859 
00860         bu_vls_free( &substr );
00861 
00862         return(0);
00863 }
00864 
00865 /**
00866  *                      R T _ E B M _ I F R E E
00867  *
00868  *  Free the storage associated with the rt_db_internal version of this solid.
00869  */
00870 void
00871 rt_ebm_ifree(struct rt_db_internal *ip)
00872 {
00873         register struct rt_ebm_internal *eip;
00874 
00875         RT_CK_DB_INTERNAL(ip);
00876         eip = (struct rt_ebm_internal *)ip->idb_ptr;
00877         RT_EBM_CK_MAGIC(eip);
00878 
00879         if(eip->mp)  {
00880                 BU_CK_MAPPED_FILE(eip->mp);
00881                 bu_close_mapped_file(eip->mp);
00882         }
00883 
00884         eip->magic = 0;                 /* sanity */
00885         eip->mp = (struct bu_mapped_file *)0;
00886         bu_free( (char *)eip, "ebm ifree" );
00887         ip->idb_ptr = GENPTR_NULL;      /* sanity */
00888         ip->idb_type = ID_NULL;         /* sanity */
00889 }
00890 
00891 /**
00892  *                      R T _ E B M _ P R E P
00893  *
00894  *  Returns -
00895  *      0       OK
00896  *      !0      Failure
00897  *
00898  *  Implicit return -
00899  *      A struct rt_ebm_specific is created, and it's address is stored
00900  *      in stp->st_specific for use by rt_ebm_shot().
00901  */
00902 int
00903 rt_ebm_prep(struct soltab *stp, struct rt_db_internal *ip, struct rt_i *rtip)
00904 {
00905         struct rt_ebm_internal  *eip;
00906         register struct rt_ebm_specific *ebmp;
00907         vect_t  norm;
00908         vect_t  radvec;
00909         vect_t  diam;
00910         vect_t  small1;
00911 
00912         eip = (struct rt_ebm_internal *)ip->idb_ptr;
00913         RT_EBM_CK_MAGIC(eip);
00914 
00915         BU_GETSTRUCT( ebmp, rt_ebm_specific );
00916         ebmp->ebm_i = *eip;             /* struct copy */
00917 
00918         /* "steal" the bitmap storage */
00919         eip->mp = (struct bu_mapped_file *)0;   /* "steal" the mapped file */
00920 
00921         /* build Xform matrix from model(world) to ideal(local) space */
00922         bn_mat_inv( ebmp->ebm_mat, eip->mat );
00923 
00924         /* Pre-compute the necessary normals.  Rotate only. */
00925         VSET( norm, 1, 0 , 0 );
00926         MAT3X3VEC( ebmp->ebm_xnorm, eip->mat, norm );
00927         VSET( norm, 0, 1, 0 );
00928         MAT3X3VEC( ebmp->ebm_ynorm, eip->mat, norm );
00929         VSET( norm, 0, 0, 1 );
00930         MAT3X3VEC( ebmp->ebm_znorm, eip->mat, norm );
00931 
00932         stp->st_specific = (genptr_t)ebmp;
00933 
00934         /* Find bounding RPP of rotated local RPP */
00935         VSETALL( small1, 0 );
00936         VSET( ebmp->ebm_large, ebmp->ebm_i.xdim, ebmp->ebm_i.ydim, ebmp->ebm_i.tallness );
00937         bn_rotate_bbox( stp->st_min, stp->st_max, eip->mat,
00938                 small1, ebmp->ebm_large );
00939 
00940         /* for now, EBM origin in ideal coordinates is at origin */
00941         VSETALL( ebmp->ebm_origin, 0 );
00942         VADD2( ebmp->ebm_large, ebmp->ebm_large, ebmp->ebm_origin );
00943 
00944         /* for now, EBM cell size in ideal coordinates is one unit/cell */
00945         VSETALL( ebmp->ebm_cellsize, 1 );
00946 
00947         VSUB2( diam, stp->st_max, stp->st_min );
00948         VADD2SCALE( stp->st_center, stp->st_min, stp->st_max, 0.5 );
00949         VSCALE( radvec, diam, 0.5 );
00950         stp->st_aradius = stp->st_bradius = MAGNITUDE( radvec );
00951 
00952         return(0);              /* OK */
00953 }
00954 
00955 /**
00956  *                      R T _ E B M _ P R I N T
00957  */
00958 void
00959 rt_ebm_print(register const struct soltab *stp)
00960 {
00961         register const struct rt_ebm_specific *ebmp =
00962                 (struct rt_ebm_specific *)stp->st_specific;
00963 
00964         bu_log("ebm file = %s\n", ebmp->ebm_i.file );
00965         bu_log("dimensions = (%d, %d, %g)\n",
00966                 ebmp->ebm_i.xdim, ebmp->ebm_i.ydim,
00967                 ebmp->ebm_i.tallness );
00968         VPRINT("model cellsize", ebmp->ebm_cellsize);
00969         VPRINT("model grid origin", ebmp->ebm_origin);
00970 }
00971 
00972 /**
00973  *                      R T _ E B M _ S H O T
00974  *
00975  *  Intersect a ray with an extruded bitmap.
00976  *  If intersection occurs, a pointer to a sorted linked list of
00977  *  "struct seg"s will be returned.
00978  *
00979  *  Returns -
00980  *      0       MISS
00981  *      >0      HIT
00982  */
00983 int
00984 rt_ebm_shot(struct soltab *stp, register struct xray *rp, struct application *ap, struct seg *seghead)
00985 {
00986         register struct rt_ebm_specific *ebmp =
00987                 (struct rt_ebm_specific *)stp->st_specific;
00988         vect_t          norm;
00989         struct xray     ideal_ray;
00990         struct seg      myhead;
00991         int             i;
00992 
00993         BU_LIST_INIT( &(myhead.l) );
00994 
00995         /* Transform actual ray into ideal space at origin in X-Y plane */
00996         MAT4X3PNT( ideal_ray.r_pt, ebmp->ebm_mat, rp->r_pt );
00997         MAT4X3VEC( ideal_ray.r_dir, ebmp->ebm_mat, rp->r_dir );
00998 
00999 #if 0
01000 bu_log("%g %g %g %g %g %g\n",
01001 ideal_ray.r_pt[X], ideal_ray.r_pt[Y], ideal_ray.r_pt[Z],
01002 ideal_ray.r_dir[X], ideal_ray.r_dir[Y], ideal_ray.r_dir[Z] );
01003 #endif
01004         if( rt_ebm_dda( &ideal_ray, stp, ap, &myhead ) <= 0 )
01005                 return(0);
01006 
01007         VSET( norm, 0, 0, -1 );         /* letters grow in +z, which is "inside" the halfspace */
01008         i = rt_seg_planeclip( seghead, &myhead, norm, 0.0, ebmp->ebm_i.tallness,
01009                 &ideal_ray, ap );
01010 #if 0
01011         if( segp )  {
01012                 vect_t  a, b;
01013                 /* Plot where WE think the ray goes */
01014                 VJOIN1( a, rp->r_pt, segp->seg_in.hit_dist, rp->r_dir );
01015                 VJOIN1( b, rp->r_pt, segp->seg_out.hit_dist, rp->r_dir );
01016                 pl_color( stdout, 0, 0, 255 );  /* B */
01017                 pdv_3line( stdout, a, b );
01018         }
01019 #endif
01020         return(i);
01021 }
01022 
01023 /**
01024  *                      R T _ E B M _ N O R M
01025  *
01026  *  Given one ray distance, return the normal and
01027  *  entry/exit point.
01028  *  This is mostly a matter of translating the stored
01029  *  code into the proper normal.
01030  */
01031 void
01032 rt_ebm_norm(register struct hit *hitp, struct soltab *stp, register struct xray *rp)
01033 {
01034         register struct rt_ebm_specific *ebmp =
01035                 (struct rt_ebm_specific *)stp->st_specific;
01036 
01037         VJOIN1( hitp->hit_point, rp->r_pt, hitp->hit_dist, rp->r_dir );
01038 
01039         switch( hitp->hit_surfno )  {
01040         case NORM_XPOS:
01041                 VMOVE( hitp->hit_normal, ebmp->ebm_xnorm );
01042                 break;
01043         case NORM_XNEG:
01044                 VREVERSE( hitp->hit_normal, ebmp->ebm_xnorm );
01045                 break;
01046 
01047         case NORM_YPOS:
01048                 VMOVE( hitp->hit_normal, ebmp->ebm_ynorm );
01049                 break;
01050         case NORM_YNEG:
01051                 VREVERSE( hitp->hit_normal, ebmp->ebm_ynorm );
01052                 break;
01053 
01054         case NORM_ZPOS:
01055                 VMOVE( hitp->hit_normal, ebmp->ebm_znorm );
01056                 break;
01057         case NORM_ZNEG:
01058                 VREVERSE( hitp->hit_normal, ebmp->ebm_znorm );
01059                 break;
01060 
01061         default:
01062                 bu_log("ebm_norm(%s): surfno=%d bad\n",
01063                         stp->st_name, hitp->hit_surfno );
01064                 VSETALL( hitp->hit_normal, 0 );
01065                 break;
01066         }
01067 }
01068 
01069 /**
01070  *                      R T _ E B M _ C U R V E
01071  *
01072  *  Everything has sharp edges.  This makes things easy.
01073  */
01074 void
01075 rt_ebm_curve(register struct curvature *cvp, register struct hit *hitp, struct soltab *stp)
01076 {
01077 /*      register struct rt_ebm_specific *ebmp =
01078                 (struct rt_ebm_specific *)stp->st_specific; */
01079 
01080         bn_vec_ortho( cvp->crv_pdir, hitp->hit_normal );
01081         cvp->crv_c1 = cvp->crv_c2 = 0;
01082 }
01083 
01084 /**
01085  *                      R T _ E B M _ U V
01086  *
01087  *  Map the hit point in 2-D into the range 0..1
01088  *  untransformed X becomes U, and Y becomes V.
01089  */
01090 void
01091 rt_ebm_uv(struct application *ap, struct soltab *stp, register struct hit *hitp, register struct uvcoord *uvp)
01092 {
01093         uvp->uv_u = hitp->hit_vpriv[X];
01094         uvp->uv_v = hitp->hit_vpriv[Y];
01095 
01096         /* XXX should compute this based upon footprint of ray in ebm space */
01097         uvp->uv_du = 0.0;
01098         uvp->uv_dv = 0.0;
01099 }
01100 
01101 /**
01102  *                      R T _ E B M _ F R E E
01103  */
01104 void
01105 rt_ebm_free(struct soltab *stp)
01106 {
01107         register struct rt_ebm_specific *ebmp =
01108                 (struct rt_ebm_specific *)stp->st_specific;
01109 
01110         BU_CK_MAPPED_FILE(ebmp->ebm_i.mp);
01111         bu_close_mapped_file(ebmp->ebm_i.mp);
01112 
01113         bu_free( (char *)ebmp, "rt_ebm_specific" );
01114 }
01115 
01116 int
01117 rt_ebm_class(void)
01118 {
01119         return(0);
01120 }
01121 
01122 /**
01123  *                      R T _ E B M _ P L O T
01124  */
01125 int
01126 rt_ebm_plot(struct bu_list *vhead, struct rt_db_internal *ip, const struct rt_tess_tol *ttol, const struct bn_tol *tol)
01127 {
01128         register struct rt_ebm_internal *eip;
01129         register int    x,y;
01130         register int    following;
01131         register int    base;
01132 /*      int             i; */
01133 
01134         RT_CK_DB_INTERNAL(ip);
01135         eip = (struct rt_ebm_internal *)ip->idb_ptr;
01136         RT_EBM_CK_MAGIC(eip);
01137 
01138         /* Find vertical lines */
01139         base = 0;       /* lint */
01140         for( x=0; x <= eip->xdim; x++ )  {
01141                 following = 0;
01142                 for( y=0; y <= eip->ydim; y++ )  {
01143                         if( following )  {
01144                                 if( (BIT( eip, x-1, y )==0) != (BIT( eip, x, y )==0) )
01145                                         continue;
01146                                 rt_ebm_plate( x, base, x, y, eip->tallness,
01147                                         eip->mat, vhead );
01148                                 following = 0;
01149                         } else {
01150                                 if( (BIT( eip, x-1, y )==0) == (BIT( eip, x, y )==0) )
01151                                         continue;
01152                                 following = 1;
01153                                 base = y;
01154                         }
01155                 }
01156         }
01157 
01158         /* Find horizontal lines */
01159         for( y=0; y <= eip->ydim; y++ )  {
01160                 following = 0;
01161                 for( x=0; x <= eip->xdim; x++ )  {
01162                         if( following )  {
01163                                 if( (BIT( eip, x, y-1 )==0) != (BIT( eip, x, y )==0) )
01164                                         continue;
01165                                 rt_ebm_plate( base, y, x, y, eip->tallness,
01166                                         eip->mat, vhead );
01167                                 following = 0;
01168                         } else {
01169                                 if( (BIT( eip, x, y-1 )==0) == (BIT( eip, x, y )==0) )
01170                                         continue;
01171                                 following = 1;
01172                                 base = x;
01173                         }
01174                 }
01175         }
01176         return(0);
01177 }
01178 
01179 /* either x1==x2, or y1==y2 */
01180 void
01181 rt_ebm_plate(int x1, int y1, int x2, int y2, double t, register fastf_t *mat, register struct bu_list *vhead)
01182 {
01183         LOCAL point_t   s, p;
01184         LOCAL point_t   srot, prot;
01185 
01186         VSET( s, x1, y1, 0.0 );
01187         MAT4X3PNT( srot, mat, s );
01188         RT_ADD_VLIST( vhead, srot, BN_VLIST_LINE_MOVE );
01189 
01190         VSET( p, x1, y1, t );
01191         MAT4X3PNT( prot, mat, p );
01192         RT_ADD_VLIST( vhead, prot, BN_VLIST_LINE_DRAW );
01193 
01194         VSET( p, x2, y2, t );
01195         MAT4X3PNT( prot, mat, p );
01196         RT_ADD_VLIST( vhead, prot, BN_VLIST_LINE_DRAW );
01197 
01198         p[Z] = 0;
01199         MAT4X3PNT( prot, mat, p );
01200         RT_ADD_VLIST( vhead, prot, BN_VLIST_LINE_DRAW );
01201 
01202         RT_ADD_VLIST( vhead, srot, BN_VLIST_LINE_DRAW );
01203 }
01204 
01205 struct ebm_edge
01206 {
01207         struct bu_list  l;
01208         int             x1,y1;
01209         int             x2,y2;
01210         int             left;   /* 1=>material to left, 0=>material to right */
01211         struct vertex   *v;     /* vertex at x1,y1 */
01212 };
01213 
01214 /* either x1==x2, or y1==y2 */
01215 static void
01216 rt_ebm_edge(int x1, int y1, int x2, int y2, int left, struct ebm_edge *edges)
01217 {
01218         struct ebm_edge *new_edge;
01219 
01220         new_edge = (struct ebm_edge *)bu_malloc( sizeof( struct ebm_edge ) , "rt_ebm_tess: new_edge" );
01221 
01222         /* make all edges go from lower values to larger */
01223         if( y1 < y2 || x1 < x2 )
01224         {
01225                 new_edge->x1 = x1;
01226                 new_edge->y1 = y1;
01227                 new_edge->x2 = x2;
01228                 new_edge->y2 = y2;
01229                 new_edge->left = left;
01230         }
01231         else
01232         {
01233                 new_edge->x1 = x2;
01234                 new_edge->y1 = y2;
01235                 new_edge->x2 = x1;
01236                 new_edge->y2 = y1;
01237                 new_edge->left = (!left);
01238         }
01239         new_edge->v = (struct vertex *)NULL;
01240         BU_LIST_APPEND( &edges->l , &new_edge->l );
01241 }
01242 
01243 static int
01244 rt_ebm_sort_edges(struct ebm_edge *edges)
01245 {
01246         struct ebm_edge loops;
01247         int vertical;
01248         int done;
01249         int from_x,from_y,to_x,to_y;
01250         int start_x,start_y;
01251         int max_loop_length=0;
01252         int loop_length;
01253 
01254         /* create another list to hold the edges as they are sorted */
01255         BU_LIST_INIT( &loops.l );
01256 
01257         while( BU_LIST_NON_EMPTY( &edges->l ) )
01258         {
01259                 struct ebm_edge *start,*next;
01260 
01261                 /* look for a vertical edge starting in lower left (smallest x and y ) */
01262                 start = (struct ebm_edge *)NULL;
01263                 next = BU_LIST_FIRST( ebm_edge , &edges->l );
01264                 while( BU_LIST_NOT_HEAD( &next->l , &edges->l ) )
01265                 {
01266                         if( next->x1 != next->x2 )
01267                         {
01268                                 next = BU_LIST_PNEXT( ebm_edge , &next->l );
01269                                 continue;       /* not a vertical edge */
01270                         }
01271 
01272                         if( !start )
01273                                 start = next;
01274                         else if( next->x1 < start->x1 || next->y1 < start->y1 )
01275                                 start = next;
01276 
01277                         next = BU_LIST_PNEXT( ebm_edge , &next->l );
01278                 }
01279 
01280                 if( !start )
01281                         rt_bomb( "rt_ebm_tess: rt_ebm_sort_edges: no vertical edges left!\n" );
01282 
01283                 /* put starting edge on the loop list */
01284                 BU_LIST_DEQUEUE( &start->l );
01285                 BU_LIST_INSERT( &loops.l , &start->l );
01286 
01287                 next = (struct ebm_edge *)NULL;
01288                 vertical = 0;   /* look for horizontal edge */
01289                 done = 0;
01290                 to_x = start->x2;
01291                 to_y = start->y2;
01292                 from_x = start->x1;
01293                 from_y = start->y1;
01294                 start_x = from_x;
01295                 start_y = from_y;
01296                 loop_length = 1;
01297                 while( !done )
01298                 {
01299                         struct ebm_edge *e,*e_poss[2];
01300                         int poss;
01301 
01302                         /* now find an edge that starts where this one stops (at to_x,to_y) */
01303                         poss = 0;
01304                         for( BU_LIST_FOR( e , ebm_edge , &edges->l ) )
01305                         {
01306                                 if( (vertical && e->y1 == e->y2) ||
01307                                    (!vertical && e->x1 == e->x2) )
01308                                         continue;
01309 
01310                                 if( (e->x1 == to_x && e->y1 == to_y) ||
01311                                     (e->x2 == to_x && e->y2 == to_y) )
01312                                         e_poss[poss++] = e;
01313                                 if( poss > 2 )
01314                                         rt_bomb( "rt_ebm_tess: rt_ebm_sort_edges: too many edges at one point\n" );
01315                         }
01316 
01317                         if( poss == 0 )
01318                                 rt_bomb( "rt_ebm_tess: rt_ebm_sort_edges: no edge to continue loop\n" );
01319                         if( poss == 1 )
01320                         {
01321                                 next = e_poss[0];
01322                         }
01323                         else
01324                         {
01325                                 /* must choose between two possibilities */
01326                                 if( vertical )
01327                                 {
01328                                         if( to_x < from_x )
01329                                         {
01330                                                 if( e_poss[0]->y1 > to_y || e_poss[0]->y2 > to_y )
01331                                                         next = e_poss[0];
01332                                                 else
01333                                                         next = e_poss[1];
01334                                         }
01335                                         else
01336                                         {
01337                                                 if( e_poss[0]->y1 < to_y || e_poss[0]->y2 < to_y )
01338                                                         next = e_poss[0];
01339                                                 else
01340                                                         next = e_poss[1];
01341                                         }
01342                                 }
01343                                 else
01344                                 {
01345                                         if( to_y < from_y )
01346                                         {
01347                                                 if( e_poss[0]->x1 < to_x || e_poss[0]->x2 < to_x )
01348                                                         next = e_poss[0];
01349                                                 else
01350                                                         next = e_poss[1];
01351                                         }
01352                                         else
01353                                         {
01354                                                 if( e_poss[0]->x1 > to_x || e_poss[0]->x2 > to_x )
01355                                                         next = e_poss[0];
01356                                                 else
01357                                                         next = e_poss[1];
01358                                         }
01359                                 }
01360                         }
01361 
01362                         if( next->x2 == to_x && next->y2 == to_y )
01363                         {
01364                                 /* reverse direction of edge */
01365                                 next->x2 = next->x1;
01366                                 next->y2 = next->y1;
01367                                 next->x1 = to_x;
01368                                 next->y1 = to_y;
01369                                 next->left = (!next->left);
01370                         }
01371                         to_x = next->x2;
01372                         to_y = next->y2;
01373                         from_x = next->x1;
01374                         from_y = next->y1;
01375                         loop_length++;
01376 
01377                         BU_LIST_DEQUEUE( &next->l );
01378                         BU_LIST_INSERT( &loops.l , &next->l );
01379 
01380                         if( to_x == start_x && to_y == start_y )
01381                         {
01382                                 if( loop_length > max_loop_length )
01383                                         max_loop_length = loop_length;
01384                                 done = 1;       /* complete loop */
01385                         }
01386 
01387                         if( vertical )
01388                                 vertical = 0;
01389                         else
01390                                 vertical = 1;
01391                 }
01392         }
01393 
01394         /* move sorted list back to "edges" */
01395         BU_LIST_INSERT_LIST( &edges->l , &loops.l );
01396 
01397         return( max_loop_length );
01398 }
01399 
01400 /**
01401  *                      R T _ E B M _ T E S S
01402  */
01403 int
01404 rt_ebm_tess(struct nmgregion **r, struct model *m, struct rt_db_internal *ip, const struct rt_tess_tol *ttol, const struct bn_tol *tol)
01405 {
01406         struct rt_ebm_internal  *eip;
01407         struct shell    *s;
01408         struct faceuse  *fu=(struct faceuse*)NULL;
01409         register int    i;
01410         struct vertex   ***vertp;       /* dynam array of ptrs to pointers */
01411         struct vertex   **loop_verts;
01412         struct ebm_edge edges;          /* list of edges */
01413         struct ebm_edge *e,*start_loop;
01414         int             start,x,y,left;
01415         int             max_loop_length;
01416         int             loop_length;
01417         vect_t          height,h;
01418 
01419         BN_CK_TOL( tol );
01420         NMG_CK_MODEL( m );
01421 
01422         RT_CK_DB_INTERNAL(ip);
01423         eip = (struct rt_ebm_internal *)ip->idb_ptr;
01424         RT_EBM_CK_MAGIC(eip);
01425 
01426         BU_LIST_INIT( &edges.l );
01427 
01428         x = 0;
01429         while( x <= eip->xdim )
01430         {
01431                 y = 0;
01432                 while( y <= eip->ydim )
01433                 {
01434                         if( (BIT( eip , x-1 , y ) == 0 ) != (BIT( eip , x , y ) == 0 ) )
01435                         {
01436                                 /* a vertical edge starts here */
01437                                 start = y;
01438                                 left = (BIT( eip , x , y ) != 0 );
01439 
01440                                 /* find other end */
01441                                 while( (BIT( eip , x-1 , y ) == 0 ) != (BIT( eip , x , y ) == 0 ) &&
01442                                         (BIT( eip , x , y ) != 0 ) == left )
01443                                                 y++;
01444                                 rt_ebm_edge( x , start , x , y , left , &edges );
01445                         }
01446                         else
01447                                 y++;
01448                 }
01449                 x++;
01450         }
01451 
01452         y = 0;
01453         while( y <= eip->ydim )
01454         {
01455                 x = 0;
01456                 while( x <= eip->xdim )
01457                 {
01458                         if( (BIT( eip, x, y-1 )==0) != (BIT( eip, x, y )==0) )
01459                         {
01460                                 /* a horizontal edge starts here */
01461                                 start = x;
01462                                 left = (BIT( eip , x , y-1 ) != 0 );
01463 
01464                                 /* find other end */
01465                                 while( (BIT( eip, x, y-1 )==0) != (BIT( eip, x, y )==0) &&
01466                                         (BIT( eip , x , y-1 ) != 0 ) == left )
01467                                                 x++;
01468                                 rt_ebm_edge( start , y , x , y , left , &edges );
01469                         }
01470                         else
01471                                 x++;
01472                 }
01473                 y++;
01474         }
01475 
01476         /* Sort edges into loops */
01477         max_loop_length = rt_ebm_sort_edges( &edges );
01478 
01479         /* Make NMG structures */
01480 
01481         /* make region, shell, vertex */
01482         *r = nmg_mrsv( m );
01483         s = BU_LIST_FIRST(shell, &(*r)->s_hd);
01484 
01485 
01486         vertp = (struct vertex ***)bu_calloc( max_loop_length , sizeof( struct vertex **) ,
01487                 "rt_ebm_tess: vertp" );
01488         loop_verts = (struct vertex **)bu_calloc( max_loop_length , sizeof( struct vertex *),
01489                 "rt_ebm_tess: loop_verts" );
01490 
01491         e = BU_LIST_FIRST( ebm_edge , &edges.l );
01492         while( BU_LIST_NOT_HEAD( &e->l , &edges.l ) )
01493         {
01494                 start_loop = e;
01495                 loop_length = 0;
01496                 vertp[loop_length++] = &start_loop->v;
01497 
01498                 e = BU_LIST_PNEXT( ebm_edge , &start_loop->l );
01499                 while( BU_LIST_NOT_HEAD( &e->l , &edges.l ) )
01500                 {
01501                         vertp[loop_length++] = &e->v;
01502                         if( e->x2 == start_loop->x1 && e->y2 == start_loop->y1 )
01503                         {
01504                                 struct faceuse *fu1;
01505                                 struct ebm_edge *e1;
01506                                 point_t pt_ebm,pt_model;
01507                                 int done=0;
01508 
01509                                 if( e->left )
01510                                 {
01511                                         /* make a face */
01512                                         fu1 = nmg_cmface( s , vertp , loop_length );
01513                                         NMG_CK_FACEUSE( fu1 );
01514 
01515                                         /* assign geometry to the vertices used in this face */
01516                                         e1 = start_loop;
01517                                         while( !done )
01518                                         {
01519                                                 if( e1->v )
01520                                                 {
01521                                                         VSET( pt_ebm , e1->x1 , e1->y1 , 0.0 );
01522                                                         MAT4X3PNT( pt_model , eip->mat , pt_ebm );
01523                                                         nmg_vertex_gv( e1->v , pt_model );
01524                                                 }
01525                                                 if( e1 == e )
01526                                                         done = 1;
01527                                                 else
01528                                                         e1 = BU_LIST_PNEXT( ebm_edge , &e1->l );
01529                                         }
01530 
01531                                         /* assign face geometry */
01532                                         if( nmg_fu_planeeqn( fu1 , tol ) )
01533                                                 goto fail;
01534 
01535                                         if( !fu )
01536                                                 fu = fu1;
01537                                 }
01538                                 else
01539                                 {
01540                                         /* make a hole */
01541                                         for( i=0 ; i<loop_length ; i++ )
01542                                         {
01543                                                 if( *vertp[loop_length-i-1] )
01544                                                         loop_verts[i] = (*vertp[loop_length-i-1]);
01545                                                 else
01546                                                         loop_verts[i] = (struct vertex *)NULL;
01547                                         }
01548 
01549                                         (void) nmg_add_loop_to_face( s , fu , loop_verts ,
01550                                                         loop_length , OT_OPPOSITE );
01551 
01552                                         /* Assign geometry to new vertices */
01553                                         done = 0;
01554                                         e1 = start_loop;
01555                                         i = loop_length - 1;
01556                                         while( !done )
01557                                         {
01558                                                 if( !loop_verts[i]->vg_p )
01559                                                 {
01560                                                         VSET( pt_ebm , e1->x1 , e1->y1 , 0.0 );
01561                                                         MAT4X3PNT( pt_model , eip->mat , pt_ebm );
01562                                                         nmg_vertex_gv( loop_verts[i] , pt_model );
01563                                                         e1->v = loop_verts[i];
01564                                                 }
01565                                                 if( e1 == e )
01566                                                         done = 1;
01567                                                 else
01568                                                 {
01569                                                         e1 = BU_LIST_PNEXT( ebm_edge , &e1->l );
01570                                                         i--;
01571                                                 }
01572                                         }
01573                                 }
01574                                 break;
01575                         }
01576                         e = BU_LIST_PNEXT( ebm_edge , &e->l );
01577                 }
01578                 e = BU_LIST_PNEXT( ebm_edge , &e->l );
01579         }
01580 
01581         /* all faces should merge into one */
01582         nmg_shell_coplanar_face_merge( s , tol , 1 );
01583 
01584         fu = BU_LIST_FIRST( faceuse , &s->fu_hd );
01585         NMG_CK_FACEUSE( fu );
01586 
01587         VSET( h , 0.0 , 0.0 , eip->tallness );
01588         MAT4X3VEC( height , eip->mat , h );
01589 
01590         nmg_extrude_face( fu , height , tol );
01591 
01592         nmg_region_a( *r , tol );
01593 
01594         (void)nmg_mark_edges_real( &s->l.magic );
01595 
01596         bu_free( (char *)vertp , "rt_ebm_tess: vertp" );
01597         bu_free( (char *)loop_verts , "rt_ebm_tess: loop_verts" );
01598         while( BU_LIST_NON_EMPTY( &edges.l ) )
01599         {
01600                 e = BU_LIST_FIRST( ebm_edge , &edges.l );
01601                 BU_LIST_DEQUEUE( &e->l );
01602                 bu_free( (char *)e , "rt_ebm_tess: e" );
01603         }
01604         return( 0 );
01605 
01606 fail:
01607         bu_free( (char *)vertp , "rt_ebm_tess: vertp" );
01608         bu_free( (char *)loop_verts , "rt_ebm_tess: loop_verts" );
01609         while( BU_LIST_NON_EMPTY( &edges.l ) )
01610         {
01611                 e = BU_LIST_FIRST( ebm_edge , &edges.l );
01612                 BU_LIST_DEQUEUE( &e->l );
01613                 bu_free( (char *)e , "rt_ebm_tess: e" );
01614         }
01615 
01616         return(-1);
01617 }
01618 
01619 /******** Test Driver *********/
01620 #ifdef TEST_DRIVER
01621 
01622 FILE                    *plotfp;
01623 
01624 struct soltab           Tsolid;
01625 struct directory        Tdir;
01626 struct rt_ebm_specific  *bmsp;
01627 struct resource         resource;
01628 mat_t                   Tmat;
01629 struct application      Tappl;
01630 
01631 main( int argc, char * *argv )
01632 {
01633         point_t pt1;
01634         point_t pt2;
01635         register int    x, y;
01636         fastf_t         xx, yy;
01637         mat_t           mat;
01638         register struct rt_ebm_specific *ebmp;
01639         int             arg;
01640         FILE            *fp;
01641         union record    rec;
01642 
01643 
01644         if( argc > 1 )  {
01645                 rt_g.debug |= DEBUG_EBM;
01646                 arg = atoi(argv[1]);
01647         }
01648 
01649         plotfp = fopen( "ebm.pl", "w" );
01650 
01651         RT_DIR_SET_NAMEP(&Tdir, "Tsolid");
01652         Tsolid.st_dp = &Tdir;
01653 
01654         RT_APPLICATION_INIT(&Tappl);
01655 
01656         Tappl.a_purpose = "testing";
01657         Tappl.a_resource = &resource;
01658         Tsolid.st_matp = &Tmat;
01659         MAT_IDN( Tsolid.st_matp );
01660 
01661         strcpy( rec.ss.ss_keyword, "ebm" );
01662         strcpy( rec.ss.ss_args, "file=bm.bw w=6 n=6 d=6.0" );
01663 
01664         if( rt_ebm_prep( &Tsolid, &rec, 0 ) != 0 )  {
01665                 printf("prep failed\n");
01666                 exit(1);
01667         }
01668         rt_ebm_print( &Tsolid );
01669         ebmp = bmsp = (struct rt_ebm_specific *)Tsolid.st_specific;
01670 
01671         outline( Tsolid.st_matp, &rec );
01672 
01673 #if 1
01674         if( (fp = fopen("ebm.rays", "r")) == NULL )  {
01675                 perror("ebm.rays");
01676                 exit(1);
01677         }
01678         for(;;)  {
01679                 x = fscanf( fp, "%lf %lf %lf %lf %lf %lf\n",
01680                         &pt1[X], &pt1[Y], &pt1[Z],
01681                         &pt2[X], &pt2[Y], &pt2[Z] );
01682                 if( x < 6 )  break;
01683                 VADD2( pt2, pt2, pt1 );
01684                 trial( pt1, pt2 );
01685         }
01686 #endif
01687 #if 0
01688         y = arg;
01689         for( x=0; x<=ebmp->ebm_i.xdim; x++ )  {
01690                 VSET( pt1, 0, y, 1 );
01691                 VSET( pt2, x, 0, 1 );
01692                 trial( pt1, pt2 );
01693         }
01694 #endif
01695 #if 0
01696         y = arg;
01697         for( x=0; x<=ebmp->ebm_i.xdim; x++ )  {
01698                 VSET( pt1, 0, y, 2 );
01699                 VSET( pt2, x, 0, 4 );
01700                 trial( pt1, pt2 );
01701         }
01702 #endif
01703 #if 0
01704         for( y=0; y<=ebmp->ebm_i.ydim; y++ )  {
01705                 for( x=0; x<=ebmp->ebm_i.xdim; x++ )  {
01706                         VSET( pt1, 0, y, 2 );
01707                         VSET( pt2, x, 0, 4 );
01708                         trial( pt1, pt2 );
01709                 }
01710         }
01711 #endif
01712 #if 0
01713         for( y= -1; y<=ebmp->ebm_i.ydim; y++ )  {
01714                 for( x= -1; x<=ebmp->ebm_i.xdim; x++ )  {
01715                         VSET( pt1, x, y, 10 );
01716                         VSET( pt2, x+2, y+3, -4 );
01717                         trial( pt1, pt2 );
01718                 }
01719         }
01720 #endif
01721 #if 0
01722         for( y=0; y<=ebmp->ebm_i.ydim; y++ )  {
01723                 for( x=0; x<=ebmp->ebm_i.xdim; x++ )  {
01724                         VSET( pt1, ebmp->ebm_i.xdim, y, 2 );
01725                         VSET( pt2, x, ebmp->ebm_i.ydim, 4 );
01726                         trial( pt1, pt2 );
01727                 }
01728         }
01729 #endif
01730 
01731 #if 0
01732         for( yy = 2.0; yy < 6.0; yy += 0.3 )  {
01733                 VSET( pt1, 0, yy, 2 );
01734                 VSET( pt2, 6, 0, 4 );
01735                 trial( pt1, pt2 );
01736         }
01737 #endif
01738 #if 0
01739         for( yy=0; yy<=ebmp->ebm_i.ydim; yy += 0.3 )  {
01740                 for( xx=0; xx<=ebmp->ebm_i.xdim; xx += 0.3 )  {
01741                         VSET( pt1, 0, yy, 2 );
01742                         VSET( pt2, xx, 0, 4 );
01743                         trial( pt1, pt2 );
01744                 }
01745         }
01746 #endif
01747 #if 0
01748         for( yy=0; yy<=ebmp->ebm_i.ydim; yy += 0.3 )  {
01749                 for( xx=0; xx<=ebmp->ebm_i.xdim; xx += 0.3 )  {
01750                         VSET( pt1, ebmp->ebm_i.xdim, yy, 2 );
01751                         VSET( pt2, xx, ebmp->ebm_i.ydim, 4 );
01752                         trial( pt1, pt2 );
01753                 }
01754         }
01755 #endif
01756 
01757 #if 0
01758         /* (6,5) (2,5) (3,4) (1.2,1.2)
01759          * were especially troublesome */
01760         xx=6;
01761         yy=1.5;
01762         VSET( pt1, 0, yy, 2 );
01763         VSET( pt2, xx, 0, 4 );
01764         trial( pt1, pt2 );
01765 #endif
01766 
01767 #if 0
01768         /* (1,2) (3,2) (0,2) (0,0.3)
01769          * were especially troublesome */
01770         xx=0;
01771         yy=0.3;
01772         VSET( pt1, ebmp->ebm_i.xdim, yy, 2 );
01773         VSET( pt2, xx, ebmp->ebm_i.ydim, 4 );
01774         trial( pt1, pt2 );
01775 #endif
01776 
01777 #if 0
01778         VSET( pt1, 0, 1.5, 2 );
01779         VSET( pt2, 4.75, 6, 4 );
01780         trial( pt1, pt2 );
01781 
01782         VSET( pt1, 0, 2.5, 2 );
01783         VSET( pt2, 4.5, 0, 4 );
01784         trial( pt1, pt2 );
01785 #endif
01786 
01787 
01788 #if 0
01789         /* With Z=-10, it works, but looks like trouble, due to
01790          * the Z-clipping causing lots of green vectors on the 2d proj.
01791          * VSET( pt1, 0.75, 1.1, -10 );
01792          */
01793         VSET( pt1, 0.75, 1.1, 0 );
01794         {
01795                 for( yy=0; yy<=ebmp->ebm_i.ydim; yy += 0.3 )  {
01796                         for( xx=0; xx<=ebmp->ebm_i.xdim; xx += 0.3 )  {
01797                                 VSET( pt2, xx, yy, 4 );
01798                                 trial( pt1, pt2 );
01799                         }
01800                 }
01801         }
01802 #endif
01803 #if 0
01804         for( x=0; x<ebmp->ebm_i.xdim; x++ )  {
01805                 VSET( pt1, x+0.75, 1.1, -10 );
01806                 VSET( pt2, x+0.75, 1.1, 4 );
01807                 trial( pt1, pt2 );
01808         }
01809 #endif
01810 
01811         exit(0);
01812 }
01813 
01814 trial(p1, p2)
01815 vect_t  p1, p2;
01816 {
01817         struct seg      *segp, *next;
01818         fastf_t         lastk;
01819         struct xray     ray;
01820         register struct rt_ebm_specific *ebmp = bmsp;
01821 
01822         VMOVE( ray.r_pt, p1 );
01823         VSUB2( ray.r_dir, p2, p1 );
01824         if( MAGNITUDE( ray.r_dir ) < 1.0e-10 )
01825                 return;
01826         VUNITIZE( ray.r_dir );
01827 
01828         printf("------- (%g, %g, %g) to (%g, %g, %g), dir=(%g, %g, %g)\n",
01829                 ray.r_pt[X], ray.r_pt[Y], ray.r_pt[Z],
01830                 p2[X], p2[Y], p2[Z],
01831                 ray.r_dir[X], ray.r_dir[Y], ray.r_dir[Z] );
01832 
01833 
01834         segp = rt_ebm_shot( &Tsolid, &ray, &Tappl );
01835 
01836         lastk = 0;
01837         while( segp != SEG_NULL )  {
01838                 /* Draw 2-D segments */
01839                 draw2seg( ray.r_pt, ray.r_dir, lastk, segp->seg_in.hit_dist, 0 );
01840                 draw2seg( ray.r_pt, ray.r_dir, segp->seg_in.hit_dist, segp->seg_out.hit_dist, 1 );
01841                 lastk = segp->seg_out.hit_dist;
01842 
01843                 draw3seg( segp, ray.r_pt, ray.r_dir );
01844 
01845                 next = segp->seg_next;
01846                 FREE_SEG(segp, Tappl.a_resource);
01847                 segp = next;
01848         }
01849 }
01850 
01851 outline(mat, rp)
01852 mat_t   mat;
01853 union record    *rp;
01854 {
01855         register struct rt_ebm_specific *ebmp = bmsp;
01856         register struct vlist   *vp;
01857         struct vlhead   vhead;
01858 
01859         vhead.vh_first = vhead.vh_last = VL_NULL;
01860 
01861         pl_3space( plotfp, -BIT_XWIDEN,-BIT_YWIDEN,-BIT_XWIDEN,
01862                  ebmp->ebm_i.xdim+BIT_XWIDEN, ebmp->ebm_i.ydim+BIT_YWIDEN, (int)(ebmp->ebm_i.tallness+1.99) );
01863         pl_3box( plotfp, -BIT_XWIDEN,-BIT_YWIDEN,-BIT_XWIDEN,
01864                  ebmp->ebm_i.xdim+BIT_XWIDEN, ebmp->ebm_i.ydim+BIT_YWIDEN, (int)(ebmp->ebm_i.tallness+1.99) );
01865 
01866         /* Get vlist, then just draw the vlist */
01867         rt_ebm_plot( rp, mat, &vhead, 0 );
01868 
01869         for( vp = vhead.vh_first; vp != VL_NULL; vp = vp->vl_forw )  {
01870                 if( vp->vl_draw == 0 )
01871                         pdv_3move( plotfp, vp->vl_pnt );
01872                 else
01873                         pdv_3cont( plotfp, vp->vl_pnt );
01874         }
01875         FREE_VL( vhead.vh_first );
01876 }
01877 
01878 
01879 draw2seg( pt, dir, k1, k2, inside )
01880 vect_t  pt, dir;
01881 fastf_t k1, k2;
01882 int     inside;
01883 {
01884         vect_t  a, b;
01885 
01886         a[0] = pt[0] + k1 * dir[0];
01887         a[1] = pt[1] + k1 * dir[1];
01888         a[2] = 0;
01889         b[0] = pt[0] + k2 * dir[0];
01890         b[1] = pt[1] + k2 * dir[1];
01891         b[2] = 0;
01892 
01893         if( inside )
01894                 pl_color( plotfp, 255, 0, 0 );  /* R */
01895         else
01896                 pl_color( plotfp, 0, 255, 0 );  /* G */
01897         pdv_3line( plotfp, a, b );
01898 }
01899 
01900 draw3seg( segp, pt, dir )
01901 register struct seg     *segp;
01902 vect_t  pt, dir;
01903 {
01904         vect_t  a, b;
01905 
01906         VJOIN1( a, pt, segp->seg_in.hit_dist, dir );
01907         VJOIN1( b, pt, segp->seg_out.hit_dist, dir );
01908         pl_color( plotfp, 0, 0, 255 );  /* B */
01909         pdv_3line( plotfp, a, b );
01910 }
01911 #endif /* test driver */
01912 
01913 /**
01914  *              R T _ E B M _ T C L G E T
01915  *
01916  *      Routine to format the parameters of an EBM for "db get"
01917  *
01918  *      Legal requested parameters are:
01919  *              "F" - bitmap file to extrude
01920  *              "W" - number of cells in X direction
01921  *              "N" - number of cells in Y direction
01922  *              "H" - height of each cell (mm)
01923  *              "M" - matrix to transform EBM solid into model coordinates
01924  *
01925  *      no paramaters requested returns all
01926  */
01927 int
01928 rt_ebm_tclget(Tcl_Interp *interp, const struct rt_db_internal *intern, const char *attr)
01929 {
01930         register struct rt_ebm_internal *ebm=(struct rt_ebm_internal *)intern->idb_ptr;
01931         Tcl_DString     ds;
01932         struct bu_vls   vls;
01933         int             i;
01934 
01935         RT_EBM_CK_MAGIC( ebm );
01936 
01937         Tcl_DStringInit( &ds );
01938         bu_vls_init( &vls );
01939 
01940         if( attr == (char *)NULL ) {
01941                 bu_vls_strcpy( &vls, "ebm" );
01942                 bu_vls_printf( &vls, " F %s W %d N %d H %.25g",
01943                                ebm->file, ebm->xdim, ebm->ydim, ebm->tallness );
01944                 bu_vls_printf( &vls, " M {" );
01945                 for( i=0 ; i<16 ; i++ )
01946                         bu_vls_printf( &vls, " %.25g", ebm->mat[i] );
01947                 bu_vls_printf( &vls, " }" );
01948         }
01949         else if( !strcmp( attr, "F" ) )
01950                 bu_vls_printf( &vls, "%s", ebm->file );
01951         else if( !strcmp( attr, "W" ) )
01952                 bu_vls_printf( &vls, "%d", ebm->xdim );
01953         else if( !strcmp( attr, "N" ) )
01954                 bu_vls_printf( &vls, "%d", ebm->ydim );
01955         else if( !strcmp( attr, "H" ) )
01956                 bu_vls_printf( &vls, "%.25g", ebm->tallness );
01957         else if( !strcmp( attr, "M" ) ) {
01958                 for( i=0 ; i<16 ; i++ )
01959                         bu_vls_printf( &vls, "%.25g ", ebm->mat[i] );
01960         }
01961         else {
01962                 Tcl_SetResult( interp,"ERROR: Unknown attribute, choices are F, W, N, or H\n",
01963                 TCL_STATIC );
01964                 bu_vls_free( &vls );
01965                 return( TCL_ERROR );
01966         }
01967 
01968         Tcl_DStringAppend( &ds, bu_vls_addr( &vls ), -1 );
01969         Tcl_DStringResult( interp, &ds );
01970         Tcl_DStringFree( &ds );
01971         bu_vls_free( &vls );
01972         return( TCL_OK );
01973 }
01974 
01975 
01976 /**
01977  *              R T _ E B M _ T C L A D J U S T
01978  *
01979  *      Routine to adjust the parameters of an EBM
01980  *
01981  *      Legal parameters are:
01982  *              "F" - bitmap file to extrude
01983  *              "W" - number of cells in X direction
01984  *              "N" - number of cells in Y direction
01985  *              "H" - height of each cell (mm)
01986  *              "M" - matrix to transform EBM solid into model coordinates
01987  */
01988 
01989 int
01990 rt_ebm_tcladjust(Tcl_Interp *interp, struct rt_db_internal *intern, int argc, char **argv)
01991 {
01992         struct rt_ebm_internal *ebm;
01993 
01994         RT_CK_DB_INTERNAL( intern );
01995 
01996         ebm = (struct rt_ebm_internal *)intern->idb_ptr;
01997         RT_EBM_CK_MAGIC( ebm );
01998 
01999         while( argc >= 2 ) {
02000                 if( !strcmp( argv[0], "F" ) ) {
02001                         if( strlen( argv[1] ) >= RT_EBM_NAME_LEN ) {
02002                                 Tcl_SetResult( interp,
02003                                                "ERROR: File name too long",
02004                                                TCL_STATIC );
02005                                 return( TCL_ERROR );
02006                         }
02007                         strcpy( ebm->file, argv[1] );
02008                 }
02009                 else if( !strcmp( argv[0], "W" ) ) {
02010                         ebm->xdim = atoi( argv[1] );
02011                 }
02012                 else if( !strcmp( argv[0], "N" ) ) {
02013                         ebm->ydim = atoi( argv[1] );
02014                 }
02015                 else if( !strcmp( argv[0], "H" ) ) {
02016                         ebm->tallness = atof( argv[1] );
02017                 }
02018                 else if( !strcmp( argv[0], "M" ) ) {
02019                         int len=16;
02020                         fastf_t array[16];
02021                         fastf_t *ar_ptr;
02022 
02023                         ar_ptr = array;
02024 
02025                         if( tcl_list_to_fastf_array( interp, argv[1], &ar_ptr, &len) !=
02026                             len ) {
02027                                 Tcl_SetResult( interp,
02028                                       "ERROR: incorrect number of coefficents for matrix\n",
02029                                       TCL_STATIC );
02030                                 return( TCL_ERROR );
02031                         }
02032                         MAT_COPY( ebm->mat, array )
02033                 }
02034                 else {
02035                         Tcl_SetResult( interp,
02036                               "ERROR: illegal argument, choices are F, W, N, or H\n",
02037                               TCL_STATIC );
02038                         return( TCL_ERROR );
02039                 }
02040                 argc -= 2;
02041                 argv += 2;
02042         }
02043         return( TCL_OK );
02044 }
02045 
02046 int
02047 rt_ebm_tclform( const struct rt_functab *ftp, Tcl_Interp *interp )
02048 {
02049         RT_CK_FUNCTAB(ftp);
02050 
02051         Tcl_AppendResult( interp,
02052                           "F %s W %d N %d H %f M { %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f", (char *)NULL );
02053 
02054         return TCL_OK;
02055 
02056 }
02057 
02058 
02059 /**
02060  *              R T _ E B M _ M A K E
02061  *
02062  *      Routine to make a new EBM solid. The only purpose of this routine is
02063  *      to initialize the matrix and height to legal values.
02064  */
02065 void
02066 rt_ebm_make(const struct rt_functab *ftp, struct rt_db_internal *intern, double diameter)
02067 {
02068         struct rt_ebm_internal *ebm;
02069 
02070         intern->idb_major_type = DB5_MAJORTYPE_BRLCAD;
02071         intern->idb_type = ID_EBM;
02072         BU_ASSERT(&rt_functab[intern->idb_type] == ftp);
02073 
02074         intern->idb_meth = ftp;
02075         ebm = (struct rt_ebm_internal *)bu_calloc( sizeof( struct rt_ebm_internal ), 1,
02076                                      "rt_ebm_internal");
02077         intern->idb_ptr = (genptr_t)ebm;
02078         ebm->magic = RT_EBM_INTERNAL_MAGIC;
02079         MAT_IDN( ebm->mat );
02080         ebm->tallness = 1.0;
02081 }
02082 
02083 /*@}*/
02084 /*
02085  * Local Variables:
02086  * mode: C
02087  * tab-width: 8
02088  * c-basic-offset: 4
02089  * indent-tabs-mode: t
02090  * End:
02091  * ex: shiftwidth=4 tabstop=8
02092  */

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