BRL-CAD
hf.c
Go to the documentation of this file.
1 /* H F . C
2  * BRL-CAD
3  *
4  * Copyright (c) 1994-2014 United States Government as represented by
5  * the U.S. Army Research Laboratory.
6  *
7  * This library is free software; you can redistribute it and/or
8  * modify it under the terms of the GNU Lesser General Public License
9  * version 2.1 as published by the Free Software Foundation.
10  *
11  * This library is distributed in the hope that it will be useful, but
12  * WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14  * Lesser General Public License for more details.
15  *
16  * You should have received a copy of the GNU Lesser General Public
17  * License along with this file; see the file named COPYING for more
18  * information.
19  */
20 /** @addtogroup primitives */
21 /** @{ */
22 /** @file primitives/hf/hf.c
23  *
24  * Intersect a ray with a height field, where the heights are imported
25  * from an external data file, and where some (or all) of the
26  * parameters of that data file may be read in from an external
27  * control file.
28  *
29  * Description of the external string description of the HF.
30  *
31  * There are two versions of this parse table. The string solid in
32  * the .g file can set any parameter. The indirect control file
33  * (cfile) can only set parameters relating to dfile parameters, and
34  * not to the geometric position, orientation, and scale of the HF's
35  * bounding RPP.
36  *
37  * In general, the cfile should be thought of as describing the data
38  * arrangement of the dfile, and the string solid should be thought of
39  * as describing the "geometry" of the height field's bounding RPP.
40  *
41  * The string solid is parsed first. If a cfile is present, it is
42  * parsed second, and any parameters specified in the cfile override
43  * the values taken from the string solid.
44  */
45 
46 #include "common.h"
47 
48 #include <stdlib.h>
49 #include <stddef.h>
50 #include <math.h>
51 #include <string.h>
52 #include "bio.h"
53 
54 #include "bu/cv.h"
55 #include "bu/parallel.h"
56 #include "vmath.h"
57 #include "db.h"
58 #include "nmg.h"
59 #include "rtgeom.h"
60 #include "raytrace.h"
61 
62 
63 #define HF_O(m) bu_offsetof(struct rt_hf_internal, m)
64 
65 /* All fields valid in string solid */
66 const struct bu_structparse rt_hf_parse[] = {
67  {"%s", 128, "cfile", HF_O(cfile), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL},
68  {"%s", 128, "dfile", HF_O(dfile), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL},
69  {"%s", 8, "fmt", HF_O(fmt), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL},
70  {"%d", 1, "w", HF_O(w), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL },
71  {"%d", 1, "n", HF_O(n), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL },
72  {"%d", 1, "shorts", HF_O(shorts), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL },
73  {"%f", 1, "file2mm", HF_O(file2mm), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL },
74  {"%f", 3, "v", HF_O(v), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL },
75  {"%f", 3, "x", HF_O(x), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL },
76  {"%f", 3, "y", HF_O(y), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL },
77  {"%f", 1, "xlen", HF_O(xlen), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL },
78  {"%f", 1, "ylen", HF_O(ylen), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL },
79  {"%f", 1, "zscale", HF_O(zscale), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL },
80  {"", 0, (char *)0, 0, BU_STRUCTPARSE_FUNC_NULL, NULL, NULL }
81 };
82 /* Subset of fields found in cfile */
83 const struct bu_structparse rt_hf_cparse[] = {
84  {"%s", 128, "dfile", HF_O(dfile), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL},
85  {"%s", 8, "fmt", HF_O(fmt), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL},
86  {"%d", 1, "w", HF_O(w), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL },
87  {"%d", 1, "n", HF_O(n), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL },
88  {"%d", 1, "shorts", HF_O(shorts), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL },
89  {"%f", 1, "file2mm", HF_O(file2mm), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL },
90  {"", 0, (char *)0, 0, BU_STRUCTPARSE_FUNC_NULL, NULL, NULL }
91 };
92 
93 
94 struct hf_specific {
95  vect_t hf_V; /* min vertex/origin of HF */
96  vect_t hf_VO; /* max vertex of HF */
97  vect_t hf_X; /* X direction vector */
98  fastf_t hf_Xlen; /* magnitude of HF in X direction */
99  vect_t hf_Y; /* Y Direction vector */
100  fastf_t hf_Ylen; /* magnitude of HF in Y direction */
101  vect_t hf_N; /* dir of elevation */
102  fastf_t hf_min; /* bounding box of hf solid */
104  fastf_t hf_file2mm; /* scale file elevation units to mm */
105  int hf_w; /* X dimension of file */
106  int hf_n; /* Y dimension of file */
107  int hf_shorts; /* Boolean: use shorts instead of double */
109 };
110 
111 
112 /**
113  * Convert in-memory form of a height-field (HF) to a displacement-map
114  * solid (DSP) in internal representation. There is no record in the
115  * V5 database for an HF.
116  */
117 int
118 rt_hf_to_dsp(struct rt_db_internal *db_intern)
119 {
120  struct rt_hf_internal *hip = (struct rt_hf_internal *)db_intern->idb_ptr;
121  struct rt_dsp_internal *dsp;
122  vect_t tmp;
123 
124  RT_CK_DB_INTERNAL(db_intern);
125  RT_HF_CK_MAGIC(hip);
126 
127  if (! hip->shorts) {
128  bu_log("cannot convert float HF to DSP\n");
129  return -1;
130  }
131 
132  BU_ALLOC(dsp, struct rt_dsp_internal);
133  bu_vls_init(&dsp->dsp_name);
134  bu_vls_strcat(&dsp->dsp_name, hip->dfile);
135  dsp->dsp_xcnt = hip->w;
136  dsp->dsp_ycnt = hip->n;
137  dsp->dsp_smooth = 0;
138  dsp->dsp_cuttype = DSP_CUT_DIR_ADAPT;
139  if (RT_G_DEBUG & DEBUG_HF) {
140  bu_log("Converting from cut-style lower-left/upper-right to adaptive\n");
141  }
142  dsp->dsp_datasrc = RT_DSP_SRC_FILE;
143 
144 
145  MAT_IDN(dsp->dsp_stom);
146  MAT_DELTAS_VEC(dsp->dsp_stom, hip->v); /* translate */
147 
148  /* Apply modeling transformations */
149  VUNITIZE(hip->x);
150  VSCALE(tmp, hip->x, hip->xlen/(fastf_t)(hip->w - 1));
151  dsp->dsp_stom[0] = tmp[0];
152  dsp->dsp_stom[4] = tmp[1];
153  dsp->dsp_stom[8] = tmp[2];
154  VUNITIZE(hip->y);
155  VSCALE(tmp, hip->y, hip->ylen/(fastf_t)(hip->n - 1));
156  dsp->dsp_stom[1] = tmp[0];
157  dsp->dsp_stom[5] = tmp[1];
158  dsp->dsp_stom[9] = tmp[2];
159  VCROSS(tmp, hip->x, hip->y);
160  VUNITIZE(tmp);
161 
162  /* The next line should be:
163  *
164  * VSCALE(tmp, tmp, hip->zscale * hip->file2mm);
165  *
166  * This will make the converted DSP plot in MGED agree with what
167  * he HF looks like, but the HF ignores 'zscale' in the shot
168  * routine. So we choose to duplicate the raytrace behavior
169  * (ignore zscale)
170  */
171  VSCALE(tmp, tmp, hip->file2mm);
172 
173  dsp->dsp_stom[2] = tmp[0];
174  dsp->dsp_stom[6] = tmp[1];
175  dsp->dsp_stom[10] = tmp[2];
176 
177  bn_mat_inv(dsp->dsp_mtos, dsp->dsp_stom);
178 
179  dsp->magic = RT_DSP_INTERNAL_MAGIC;
180 
181  rt_db_free_internal(db_intern);
182 
183  db_intern->idb_ptr = (void *)dsp;
184  db_intern->idb_major_type = DB5_MAJORTYPE_BRLCAD;
185  db_intern->idb_type = ID_DSP;
186  db_intern->idb_meth = &OBJ[ID_DSP];
187 
188  return 0;
189 
190 }
191 
192 /**
193  * Calculate the bounding RPP for an hf
194  */
195 int
196 rt_hf_bbox(struct rt_db_internal *ip, point_t *min_pt, point_t *max_pt, const struct bn_tol *UNUSED(tol)) {
197  struct rt_hf_internal *hip;
198  vect_t height, work;
199  vect_t hf_N, hf_X, hf_Y;
200  fastf_t hf_max;
201 
202  RT_CK_DB_INTERNAL(ip);
203  hip = (struct rt_hf_internal *)ip->idb_ptr;
204  RT_HF_CK_MAGIC(hip);
205 
206  VMOVE(hf_X, hip->x);
207  VUNITIZE(hf_X);
208  VMOVE(hf_Y, hip->y);
209  VUNITIZE(hf_Y);
210 
211  VCROSS(hf_N, hf_X, hf_Y);
212 
213  /*
214  * Locate the min-max of the HF.
215  */
216  if (hip->shorts) {
217  int max, min;
218  int len;
219  unsigned short *sp;
220  int i;
221 
222  sp = (unsigned short *)hip->mp->apbuf;
223  min = max = *sp++;
224  len = hip->w * hip->n;
225  for (i = 1; i < len; i++, sp++) {
226  if ((int)*sp > max) max=*sp;
227  if ((int)*sp < min) min=*sp;
228  }
229  hf_max = max * hip->file2mm;
230  } else {
231  fastf_t max, min;
232  int len;
233  int i;
234  fastf_t *fp;
235 
236  fp = (fastf_t *) hip->mp->apbuf;
237  min = max = *fp++;
238  len = hip->w * hip->n;
239  for (i = 1; i < len; i++, fp++) {
240  if (*fp > max) max = *fp;
241  if (*fp < min) min = *fp;
242  }
243  hf_max = max * hip->file2mm;
244  }
245 
246  VSCALE(height, hf_N, hf_max);
247 
248  VMOVE((*min_pt), hip->v);
249  VMOVE((*max_pt), hip->v);
250  VADD2(work, hip->v, height);
251  VMINMAX((*min_pt), (*max_pt), work);
252  VJOIN1(work, hip->v, hip->xlen, hf_X);
253  VMINMAX((*min_pt), (*max_pt), work);
254  VADD2(work, work, height);
255  VMINMAX((*min_pt), (*max_pt), work);
256  VJOIN1(work, hip->v, hip->ylen, hf_Y);
257  VMINMAX((*min_pt), (*max_pt), work);
258  VADD2(work, work, height);
259  VMINMAX((*min_pt), (*max_pt), work);
260  VJOIN2(work, hip->v, hip->xlen, hf_X, hip->ylen, hf_Y);
261  VMINMAX((*min_pt), (*max_pt), work);
262  VADD2(work, work, height);
263  VMINMAX((*min_pt), (*max_pt), work);
264 
265  return 0;
266 }
267 
268 
269 /**
270  * Given a pointer to a GED database record, and a transformation
271  * matrix, determine if this is a valid HF, and if so, precompute
272  * various terms of the formula.
273  *
274  * Returns -
275  * 0 HF is OK
276  * !0 Error in description
277  *
278  * Implicit return -
279  * A struct hf_specific is created, and its address is stored in
280  * stp->st_specific for use by hf_shot().
281  */
282 int
283 rt_hf_prep(struct soltab *stp, struct rt_db_internal *ip, struct rt_i *rtip)
284 {
285  struct rt_hf_internal *hip;
286  struct hf_specific *hf;
287  const struct bn_tol *tol = &rtip->rti_tol;
288  double dot;
289  vect_t height, work;
290  static int first_time = 1;
291 
292  if (first_time) {
293  first_time = 0;
294  }
295  RT_CK_DB_INTERNAL(ip);
296  hip = (struct rt_hf_internal *)ip->idb_ptr;
297  RT_HF_CK_MAGIC(hip);
298 
299  BU_GET(hf, struct hf_specific);
300  stp->st_specific = (void *) hf;
301  /*
302  * The stuff that is given to us.
303  */
304  VMOVE(hf->hf_V, hip->v);
305  VMOVE(hf->hf_X, hip->x);
306  VUNITIZE(hf->hf_X);
307  hf->hf_Xlen = hip->xlen;
308  VMOVE(hf->hf_Y, hip->y);
309  VUNITIZE(hf->hf_Y);
310  hf->hf_Ylen = hip->ylen;
311  hf->hf_file2mm = hip->file2mm;
312  hf->hf_w = hip->w;
313  hf->hf_n = hip->n;
314  hf->hf_mp = hip->mp;
316  ++hf->hf_mp->uses;
318  hf->hf_shorts = hip->shorts;
319  /*
320  * From here down, we are calculating new values on a one time
321  * basis.
322  */
323 
324  /*
325  * Start finding the location of the opposite vertex to V
326  */
327  VJOIN2(hf->hf_VO, hip->v, hip->xlen, hip->x, hip->ylen, hip->y);
328 
329  /*
330  * get the normal.
331  */
332  dot = VDOT(hf->hf_X, hf->hf_Y);
333  if (fabs(dot) >tol->perp) {
334  /* not perpendicular, bad hf */
335  bu_log("Hf(%s): X not perpendicular to Y.\n", stp->st_name);
336  BU_PUT(hf, struct hf_specific);
337  stp->st_specific = (struct hf_specific *)NULL;
338  return 1; /* BAD */
339  }
340  VCROSS(hf->hf_N, hf->hf_X, hf->hf_Y);
341  VUNITIZE(hf->hf_N); /* Not needed (?) */
342 
343  /*
344  * Locate the min-max of the HF for use in determining VO and
345  * bounding boxes et so forth.
346  */
347  if (hf->hf_shorts) {
348  int max, min;
349  int len;
350  unsigned short *sp;
351  int i;
352 
353  sp = (unsigned short *)hf->hf_mp->apbuf;
354  min = max = *sp++;
355  len = hf->hf_w * hf->hf_n;
356  for (i = 1; i < len; i++, sp++) {
357  if ((int)*sp > max) max=*sp;
358  if ((int)*sp < min) min=*sp;
359  }
360  hf->hf_min = min * hf->hf_file2mm;
361  hf->hf_max = max * hf->hf_file2mm;
362  } else {
363  fastf_t max, min;
364  int len;
365  int i;
366  fastf_t *fp;
367 
368  fp = (fastf_t *) hf->hf_mp->apbuf;
369  min = max = *fp++;
370  len = hf->hf_w * hf->hf_n;
371  for (i = 1; i < len; i++, fp++) {
372  if (*fp > max) max = *fp;
373  if (*fp < min) min = *fp;
374  }
375  hf->hf_min = min * hf->hf_file2mm;
376  hf->hf_max = max * hf->hf_file2mm;
377  }
378 
379  VSCALE(height, hf->hf_N, hf->hf_max);
380  VADD2(hf->hf_VO, hf->hf_VO, height);
381 
382  /*
383  * Now we compute the bounding box and sphere.
384  */
385  VMOVE(stp->st_min, hf->hf_V);
386  VMOVE(stp->st_max, hf->hf_V);
387  VADD2(work, hf->hf_V, height);
388  VMINMAX(stp->st_min, stp->st_max, work);
389  VJOIN1(work, hf->hf_V, hf->hf_Xlen, hf->hf_X);
390  VMINMAX(stp->st_min, stp->st_max, work);
391  VADD2(work, work, height);
392  VMINMAX(stp->st_min, stp->st_max, work);
393  VJOIN1(work, hf->hf_V, hf->hf_Ylen, hf->hf_Y);
394  VMINMAX(stp->st_min, stp->st_max, work);
395  VADD2(work, work, height);
396  VMINMAX(stp->st_min, stp->st_max, work);
397  VJOIN2(work, hf->hf_V, hf->hf_Xlen, hf->hf_X, hf->hf_Ylen, hf->hf_Y);
398  VMINMAX(stp->st_min, stp->st_max, work);
399  VADD2(work, work, height);
400  VMINMAX(stp->st_min, stp->st_max, work);
401  /* Now find the center and radius for a bounding sphere. */
402  {
403  fastf_t dx, dy, dz;
404  fastf_t f;
405 
406  VADD2SCALE(stp->st_center, stp->st_max, stp->st_min, 0.5);
407 
408  dx = (stp->st_max[X] - stp->st_min[X])/2;
409  dy = (stp->st_max[Y] - stp->st_min[Y])/2;
410  dz = (stp->st_max[Z] - stp->st_min[Z])/2;
411  f = dx;
412  if (dy > f) f = dy;
413  if (dz > f) f = dz;
414  stp->st_aradius = f;
415  stp->st_bradius = sqrt(dx*dx + dy*dy + dz*dz);
416  }
417  return 0;
418 }
419 
420 
421 void
422 rt_hf_print(const struct soltab *stp)
423 {
424  const struct hf_specific *hf =
425  (struct hf_specific *)stp->st_specific;
426  VPRINT("V", hf->hf_V);
427  VPRINT("X", hf->hf_X);
428  VPRINT("Y", hf->hf_Y);
429  VPRINT("N", hf->hf_N);
430  bu_log("XL %g\n", hf->hf_Xlen);
431  bu_log("YL %g\n", hf->hf_Ylen);
432 }
433 
434 
435 static int
436 hf_cell_shot(struct soltab *stp, struct xray *rp, struct application *ap, struct hit *hitp, int xCell, int yCell)
437 {
438  struct hf_specific *hfp =
439  (struct hf_specific *)stp->st_specific;
440 
441  fastf_t dn, abs_dn, k1st = 0, k2nd = 0, alpha, beta;
442  int dir1st, dir2nd;
443  vect_t wxb, xp;
444  vect_t tri_wn1st, tri_wn2nd, tri_BA1st, tri_BA2nd;
445  vect_t tri_CA1st, tri_CA2nd;
446  vect_t xvect, yvect, tri_A, tri_B, tri_C;
447  int fnd1, fnd2;
448  double hf2mm = hfp->hf_file2mm;
449 
450  if (ap) RT_CK_APPLICATION(ap);
451 
452  if (RT_G_DEBUG & DEBUG_HF) {
453  bu_log("hf_cell_shot(%s): %d, %d\n", stp->st_name,
454  xCell, yCell);
455  }
456  {
457  fastf_t scale;
458  scale = hfp->hf_Xlen/((double)hfp->hf_w-1);
459  VSCALE(xvect, hfp->hf_X, scale);
460  scale = hfp->hf_Ylen/((double)hfp->hf_n-1);
461  VSCALE(yvect, hfp->hf_Y, scale);
462  }
463  if (hfp->hf_shorts) {
464  unsigned short *sp;
465  sp = (unsigned short *)hfp->hf_mp->apbuf +
466  yCell * hfp->hf_w + xCell;
467 
468  /* Get the points of one of the triangles */
469 
470  /* 0, 0 -> tri_A */
471  VJOIN3(tri_A, hfp->hf_V, *sp*hf2mm, hfp->hf_N, xCell+0, xvect,
472  yCell+0, yvect);
473  sp++;
474  /* 1, 0 -> tri_B */
475  VJOIN3(tri_B, hfp->hf_V, *sp*hf2mm, hfp->hf_N, xCell+1, xvect,
476  yCell+0, yvect);
477  sp += hfp->hf_w;
478  /* 1, 1 -> tri_C */
479  VJOIN3(tri_C, hfp->hf_V, *sp*hf2mm, hfp->hf_N, xCell+1, xvect,
480  yCell+1, yvect);
481 
482  VSUB2(tri_CA1st, tri_C, tri_A);
483  VSUB2(tri_BA1st, tri_B, tri_A);
484  VCROSS(tri_wn1st, tri_BA1st, tri_CA1st);
485  VMOVE(tri_B, tri_C); /* This can optimize down. */
486  --sp;
487 
488  /* 0, 1 */
489  VJOIN3(tri_C, hfp->hf_V, *sp*hf2mm, hfp->hf_N, xCell+0, xvect,
490  yCell+1, yvect);
491  VSUB2(tri_CA2nd, tri_C, tri_A);
492 
493  VSUB2(tri_BA2nd, tri_B, tri_A);
494  VCROSS(tri_wn2nd, tri_BA2nd, tri_CA2nd);
495  } else {
496  double *fp;
497  fp = (double *)hfp->hf_mp->apbuf +
498  yCell * hfp->hf_w + xCell;
499  /* 0, 0 -> A */
500  VJOIN3(tri_A, hfp->hf_V, *fp*hf2mm, hfp->hf_N, xCell+0, xvect,
501  yCell+0, yvect);
502  fp++;
503  /* 1, 0 */
504  VJOIN3(tri_B, hfp->hf_V, *fp*hf2mm, hfp->hf_N, xCell+1, xvect,
505  yCell+0, yvect);
506  fp += hfp->hf_w;
507  /* 1, 1 */
508  VJOIN3(tri_C, hfp->hf_V, *fp*hf2mm, hfp->hf_N, xCell+1, xvect,
509  yCell+1, yvect);
510  VSUB2(tri_CA1st, tri_C, tri_A);
511  VSUB2(tri_BA1st, tri_B, tri_A);
512  VCROSS(tri_wn1st, tri_BA1st, tri_CA1st);
513  VMOVE(tri_B, tri_C); /* This can optimize down. */
514  --fp;
515  /* 0, 1 */
516  VJOIN3(tri_C, hfp->hf_V, *fp*hf2mm, hfp->hf_N, xCell+0, xvect,
517  yCell+1, yvect);
518  VSUB2(tri_CA2nd, tri_C, tri_A);
519  VSUB2(tri_BA2nd, tri_B, tri_A);
520  VCROSS(tri_wn2nd, tri_BA2nd, tri_CA2nd);
521  }
522 
523  /* 0, 1 1, 1
524  * o o
525  * _
526  * ^ //|
527  * CA2nd| //
528  * | //
529  * | BA2nd//
530  * | //
531  * | // CA1st
532  * | //
533  * | //
534  * | //
535  * | //
536  *
537  * o --------> o
538  * 0, 0 BA1st 1, 0
539  *
540  * wn1st and wn2nd are non-unit normal vectors pointing out of screen
541  */
542 
543  fnd1 = fnd2 = 0;
544 
545  /* Ray triangle intersection.
546  * See: "Graphics Gems" An Efficient Ray-Polygon Intersection P:390
547  */
548 
549  dn = VDOT(tri_wn1st, rp->r_dir); /* wn1st points out */
550  if (dn >= 0.0) {
551  dir1st = 1;
552  abs_dn=dn;
553  } else {
554  dir1st = 0;
555  abs_dn = -dn;
556  }
557 
558  /* make sure ray not parallel to plane of triangle */
559  if (abs_dn >= SQRT_SMALL_FASTF) {
560  VSUB2(wxb, tri_A, rp->r_pt);
561  VCROSS(xp, wxb, rp->r_dir);
562 
563  /* alpha = dist along CA1 to isect pt */
564  alpha = VDOT(tri_CA1st, xp);
565  if (dn < 0.0) alpha = -alpha;
566 
567  /* if pt before CA1st or beyond end of CA1st pt is
568  * outside triangle.
569  *
570  * XXX Can someone explain the "alpha <= abs_dn" part? -- Lee
571  * I know it's supposed to be determining if the point
572  * is beyond the end of CA1st, but I don't see how the math
573  * here does that.
574  */
575  if (alpha >= 0.0 && alpha <= abs_dn) {
576 
577  /* beta = dist along BA1 to isect pt */
578  beta = VDOT(tri_BA1st, xp);
579  if (dn > 0.0) beta = -beta;
580 
581  if (beta >= 0.0 && beta <= abs_dn) {
582  if (alpha + beta <= abs_dn) {
583  k1st = VDOT(wxb, tri_wn1st) / dn;
584  fnd1 = 1;
585  }
586  }
587  }
588  }
589 
590  /* XXX This is really hard to read. Need to fix this like above */
591  dn = VDOT(tri_wn2nd, rp->r_dir);
592  if (dn >= 0.0) {
593  dir2nd = 1;
594  abs_dn = dn;
595  } else {
596  dir2nd = 0;
597  abs_dn = -dn;
598  }
599  if (abs_dn < SQRT_SMALL_FASTF) goto leave;
600  VSUB2(wxb, tri_A, rp->r_pt);
601  VCROSS(xp, wxb, rp->r_dir);
602  alpha = VDOT(tri_CA2nd, xp);
603  if (dn <0.0) alpha = -alpha;
604  if (alpha < 0.0 || alpha > abs_dn) goto leave;
605  beta = VDOT(tri_BA2nd, xp);
606  if (dn > 0.0) beta = -beta;
607  if (beta < 0.0 || beta > abs_dn) goto leave;
608  if (alpha+beta > abs_dn) goto leave;
609  k2nd = VDOT(wxb, tri_wn2nd)/ dn;
610  fnd2 = 1;
611  leave:
612  if (!fnd1 && !fnd2) return 0;
613 
614  if (RT_G_DEBUG & DEBUG_HF) {
615  bu_log("hf_cell_shot: hit(%d).\n", fnd1+fnd2);
616  }
617 
618  /*
619  * XXX - This is now wrong.
620  *
621  * We have now done the ray-triangle intersection. dn1st and dn
622  * tell us the direction of the normal, <0 is in and >0 is out.
623  * k1st and k2nd tell us the distance from the start point.
624  *
625  * We are only interested in the closest hit. and that will
626  * replace the out if dn>0 or the in if dn<0.
627  */
628 
629  if (!fnd2) {
630  hitp->hit_magic = RT_HIT_MAGIC;
631  hitp->hit_dist = k1st;
632  VMOVE(hitp->hit_normal, tri_wn1st);
633  VUNITIZE(hitp->hit_normal);
634  hitp->hit_surfno = yCell*hfp->hf_w+xCell;
635  return 1;
636  }
637  if (!fnd1) {
638  hitp->hit_magic = RT_HIT_MAGIC;
639  hitp->hit_dist = k2nd;
640  VMOVE(hitp->hit_normal, tri_wn2nd);
641  VUNITIZE(hitp->hit_normal);
642  hitp->hit_surfno = yCell*hfp->hf_w+xCell;
643  return 1;
644  }
645  /*
646  * This is the two hit situation which can cause interesting
647  * problems. There are basically five different cases that must be
648  * dealt with and each one requires that the ray be classified
649  *
650  * 1) The ray has hit two different planes at two different
651  * locations (k1st != k2nd). Return both hit points.
652  * 2) The ray is going from inside to outside, return one hit point.
653  * 3) The ray is going from outside to inside, return one hit point.
654  * 4) The ray is going from inside to inside, return two hit points.
655  * 5) The ray is going from outside to outside, return two hit points.
656  */
657  hitp->hit_magic = RT_HIT_MAGIC;
658  hitp->hit_dist = k1st;
659  VMOVE(hitp->hit_normal, tri_wn1st);
660  VUNITIZE(hitp->hit_normal);
661  hitp->hit_surfno = yCell*hfp->hf_w+xCell;
662  hitp++;
663 
664  if ((fabs(k1st-k2nd) <= SMALL_FASTF) &&
665  (dir1st + dir2nd != 1)) return 1;
666 
667  hitp->hit_magic = RT_HIT_MAGIC;
668  hitp->hit_dist = k2nd;
669  VMOVE(hitp->hit_normal, tri_wn2nd);
670  VUNITIZE(hitp->hit_normal);
671  hitp->hit_surfno = yCell*hfp->hf_w+xCell;
672  return 2;
673 }
674 
675 
676 #define MAXHITS 128 /* # of surfaces hit, must be even */
677 
678 /**
679  * For the given plane of the bounding box of the hf solid, compute
680  * the hit distance and add it to the list of hits. If the plane
681  * happens to be the "Zmax" face, then the hit is really in the
682  * elevation data, and we skip it. That will be handled elsewhere.
683  */
684 static void
685 axis_plane_isect(int plane, fastf_t inout, struct xray *rp, struct hf_specific *hf, double xWidth, double yWidth, struct hit **hp, int *nhits)
686 {
687  double left, right, xx = 0, xright = 0, answer;
688  vect_t loc;
689  int CellX = 0, CellY = 0;
690 
691  if (plane == -6) return;
692 
693  if (plane == -3) {
694  (*hp)->hit_magic = RT_HIT_MAGIC;
695  (*hp)->hit_dist = inout;
696  (*hp)->hit_surfno = plane;
697  (*hp)++;
698  if ((*nhits)++>=MAXHITS) bu_bomb("g_hf.c: too many hits.\n");
699  return;
700  }
701 
702  VJOIN1(loc, rp->r_pt, inout, rp->r_dir);
703  VSUB2(loc, loc, hf->hf_V);
704 
705  /* find the left, right and xx */
706  switch (plane) {
707  case -1:
708  CellY = loc[Y]/ yWidth;
709  CellX = 0;
710  xright = yWidth;
711  xx = loc[Y] - (CellY* yWidth);
712  break;
713  case -2:
714  CellY = 0;
715  CellX = loc[X]/ xWidth;
716  xright = xWidth;
717  xx = loc[X] - CellX * xWidth;
718  break;
719  case -4:
720  CellY = loc[Y]/ yWidth;
721  CellX = hf->hf_n-1;
722  xright = yWidth;
723  xx = loc[Y] - (CellY* yWidth);
724  break;
725  case -5:
726  CellY = hf->hf_w-1;
727  CellX = loc[X]/ xWidth;
728  xright = xWidth;
729  xx = loc[X] - CellX* xWidth;
730  break;
731  }
732 
733  /* What does this indicate that it generates so much noise? */
734  if (xx < 0) {
735  bu_log("hf: xx < 0, plane = %d\n", plane);
736  }
737 
738  if (hf->hf_shorts) {
739  unsigned short *sp;
740  sp = (unsigned short *)hf->hf_mp->apbuf +
741  CellY * hf->hf_w + CellX;
742  left = *sp;
743  if (plane == -2 || plane == -5) {
744  sp++;
745  } else {
746  sp += hf->hf_w;
747  }
748  right = *sp;
749  } else {
750  double *fp;
751  fp = (double *) hf->hf_mp->apbuf +
752  CellY * hf->hf_w + CellX;
753  left = *fp;
754  if (plane == -2 || plane == -5) {
755  fp++;
756  } else {
757  fp += hf->hf_w;
758  }
759  right = *fp;
760  }
761  left *= hf->hf_file2mm;
762  right *= hf->hf_file2mm;
763  answer = (right-left)/xright*xx+left;
764 
765  if (loc[Z]-SQRT_SMALL_FASTF < answer) {
766  (*hp)->hit_magic = RT_HIT_MAGIC;
767  (*hp)->hit_dist = inout;
768  (*hp)->hit_surfno = plane;
769  VJOIN1((*hp)->hit_point, rp->r_pt, inout, rp->r_dir);
770  (*hp)++;
771  if ((*nhits)++>=MAXHITS) bu_bomb("g_hf.c: too many hits.\n");
772  }
773 }
774 
775 
776 /**
777  * Intersect a ray with a height field. If an intersection occurs, a
778  * struct seg will be acquired and filled in.
779  *
780  * Returns -
781  * 0 MISS
782  * >0 HIT
783  */
784 int
785 rt_hf_shot(struct soltab *stp, struct xray *rp, struct application *ap, struct seg *seghead)
786 {
787  struct hf_specific *hf =
788  (struct hf_specific *)stp->st_specific;
789 
790  struct hit hits[MAXHITS];
791  struct hit *hp;
792  int nhits;
793  double xWidth, yWidth;
794 
795  vect_t peqn;
796  fastf_t pdist = 0;
797  fastf_t allDist[6]; /* The hit point for all rays. */
798  fastf_t cosine;
799 
800  int iplane, oplane, j;
801  fastf_t in, out;
802  vect_t aray, curloc;
803 
804  VSETALL(peqn, 0);
805 
806  memset(hits, 0, sizeof(hits));
807 
808  in = -INFINITY;
809  out = INFINITY;
810  iplane = oplane = 0;
811 
812  nhits = 0;
813  hp = &hits[0];
814 
815 
816  /*
817  * First step in raytracing the HF is to find the intersection of
818  * the ray with the bounding box. Since the ray might not even
819  * hit the box.
820  *
821  * The results of this intercept can be used to find which cell of
822  * the dda is the start cell.
823  */
824  for (j=-1; j>-7; j--) {
825  fastf_t dn; /* Direction dot Normal */
826  fastf_t dxbdn; /* distance between d and b * dn */
827  fastf_t s; /* actual distance in mm */
828  int allIndex;
829 
830  switch (j) {
831  case -1:
832  /* Xmin plane */
833  VREVERSE(peqn, hf->hf_X);
834  pdist = VDOT(peqn, hf->hf_V);
835  break;
836  case -2:
837  /* Ymin plane */
838  VREVERSE(peqn, hf->hf_Y);
839  pdist = VDOT(peqn, hf->hf_V);
840  break;
841  case -3:
842  /* Zmin plane */
843  VREVERSE(peqn, hf->hf_N);
844  pdist = VDOT(peqn, hf->hf_V);
845  break;
846  case -4:
847  /* Xmax plane */
848  VMOVE(peqn, hf->hf_X);
849  pdist = VDOT(peqn, hf->hf_VO);
850  break;
851  case -5:
852  /* Ymax plane */
853  VMOVE(peqn, hf->hf_Y);
854  pdist = VDOT(peqn, hf->hf_VO);
855  break;
856  case -6:
857  /* Zmax plane */
858  VMOVE(peqn, hf->hf_N);
859  pdist = VDOT(peqn, hf->hf_VO);
860  break;
861  }
862  allIndex = abs(j)-1;
863 
864  dxbdn = VDOT(peqn, rp->r_pt) - pdist;
865  dn = -VDOT(peqn, rp->r_dir);
866  if (RT_G_DEBUG & DEBUG_HF) {
867  VPRINT("hf: Plane Equation", peqn);
868  bu_log("hf: dn=%g, dxbdn=%g, ", dn, dxbdn);
869  }
870 
871  if (dn < -SQRT_SMALL_FASTF) {
872  /* Leaving */
873  allDist[allIndex] = s = dxbdn/dn;
874  if (out > s) {
875  out = s;
876  oplane = j;
877  }
878  if (RT_G_DEBUG & DEBUG_HF) {
879  bu_log("s=%g out=%g\n", s, out);
880  }
881  } else if (dn > SQRT_SMALL_FASTF) {
882  /* entering */
883  allDist[allIndex] = s = dxbdn/dn;
884  if (in < s) {
885  in = s;
886  iplane = j;
887  }
888  if (RT_G_DEBUG & DEBUG_HF) {
889  bu_log("s=%g in=%g\n", s, in);
890  }
891  } else {
892  /*
893  * if the ray is outside the solid, then this
894  * is a miss.
895  */
896  if (RT_G_DEBUG & DEBUG_HF) {
897  bu_log("s=DIVIDE_BY_ZERO\n");
898  }
899  if (dxbdn > SQRT_SMALL_FASTF) {
900  return 0; /* MISS */
901  }
902  allDist[allIndex] = INFINITY;
903  }
904  if (in > out) {
905  if (RT_G_DEBUG & DEBUG_HF) {
906  bu_log("rt_hf_shoot(%s): in(%g) > out(%g)\n",
907  stp->st_name, in, out);
908  }
909  return 0; /* MISS */
910  }
911  }
912 
913  if (iplane >= 0 || oplane >= 0) {
914  bu_log("rt_hf_shoot(%s): 1 hit => MISS\n",
915  stp->st_name);
916  return 0; /* MISS */
917  }
918 
919  if (fabs(in-out) <= SMALL_FASTF || out >= INFINITY) {
920  if (RT_G_DEBUG & DEBUG_HF) {
921  bu_log("rt_hf_shoot(%s): in(%g) >= out(%g) || out >= INFINITY\n",
922  stp->st_name, in, out);
923  }
924  return 0;
925  }
926 
927  /*
928  * Once that translation is done, we start a DDA to walk across
929  * the field. Checking each "cell" and when changing cells in the
930  * off direction, the 2 cells filling the corner are also checked.
931  */
932 
933  xWidth = hf->hf_Xlen/((double)(hf->hf_w-1));
934  yWidth = hf->hf_Ylen/((double)(hf->hf_n-1));
935 
936  if (RT_G_DEBUG & DEBUG_HF) {
937  bu_log("hf: xWidth=%g, yWidth=%g, in=%g, out=%g\n", xWidth,
938  yWidth, in, out);
939  }
940 
941 
942  /*
943  * add the sides, and bottom to the hit list.
944  */
945  axis_plane_isect(iplane, in, rp, hf, xWidth, yWidth, &hp, &nhits);
946  axis_plane_isect(oplane, out, rp, hf, xWidth, yWidth, &hp, &nhits);
947 
948  /*
949  * Gee, we've gotten much closer, we know that we hit the
950  * solid. Now it's time to see which cell we hit. The Key here is
951  * to use a fast DDA to check ONLY the cells we are interested in.
952  * The basic idea and some of the pseudo code comes from:
953  *
954  * Grid Tracing: Fast Ray Tracing for Height Fields By: F. Kenton
955  * Musgrave.
956  */
957 
958  /*
959  * Now we figure out which direction we are going to be moving, X,
960  * Y or Z.
961  */
962  {
963  vect_t tmp;
964  VMOVE(tmp, rp->r_dir);
965  tmp[Z] = 0.0; /* XXX Bogus? Assumes X, Y in XY plane */
966  VUNITIZE(tmp);
967  cosine = VDOT(tmp, hf->hf_X);
968  }
969  if (fabs(cosine) <= SMALL_FASTF) {
970  /* near enough to Z */
971  vect_t tmp;
972  int xCell, yCell, r;
973  if (RT_G_DEBUG & DEBUG_HF) {
974  bu_log("hf: Vertical shoot\n");
975  }
976  VSUB2(tmp, rp->r_pt, hf->hf_V);
977  xCell = tmp[X]/hf->hf_Xlen*hf->hf_w;
978  yCell = tmp[Y]/hf->hf_Ylen*hf->hf_n;
979  r = hf_cell_shot(stp, rp, ap, hp, xCell, yCell);
980  if (r) {
981  nhits += r;
982  if (nhits > MAXHITS)
983  bu_bomb("g_hf.c: too many hits.\n");
984  hp += r;
985  }
986  } else if (cosine*cosine > 0.5) {
987  double tmp;
988  double farZ, minZ, maxZ;
989  int xCell, yCell, signX, signY;
990  double highest, lowest, error, delta;
991  double deltaZ;
992 
993  vect_t goesIN, goesOUT;
994 
995  VJOIN1(goesIN, rp->r_pt, allDist[3], rp->r_dir); /* Xmax plane */
996  VJOIN1(goesOUT, rp->r_pt, allDist[0], rp->r_dir); /* Xmin plane */
997  VSUB2(aray, goesOUT, goesIN);
998  VSUB2(curloc, goesIN, hf->hf_V);
999 
1000 
1001  /*
1002  * We will be stepping one "cell width" in the X direction
1003  * each time through the loop. In simple case we have???
1004  * cell_width = htfp->htf_Xlen/(htfp->htf_i-1); deltaX =
1005  * (Xdist*cell_width)/Xdist deltaY = (Ydist*cell_width)/Xdist;
1006  */
1007  tmp = xWidth/fabs(aray[X]);
1008 
1009  VSCALE(aray, aray, tmp);
1010 
1011  /*
1012  * Some fudges here. First, the size of the array of samples
1013  * is iXj, but the size of the array of CELLS is (i-1)X(j-1),
1014  * therefore number of CELLS is one less than number of
1015  * samples and the scaling factor is used. Second, the math
1016  * is nice to us. IF we are entering at the far end
1017  * (curloc[X] == Xlen || curloc[Y] == Ylen) then the result we
1018  * will get back is of the cell following this (out of bounds)
1019  * So we add a check for that problem.
1020  */
1021  xCell = curloc[X]/xWidth;
1022  tmp = curloc[Y];
1023  if (tmp < 0) {
1024  yCell = tmp/yWidth - 1;
1025  } else {
1026  yCell = tmp/yWidth;
1027  }
1028 
1029  signX = (aray[X] < 0.0) ? -1 : 1;
1030  signY = (aray[Y] < 0.0) ? -1 : 1;
1031 
1032  if (RT_G_DEBUG & DEBUG_HF) {
1033  bu_log("hf: curloc=(%g, %g, %g) aray=(%g, %g, %g)\n", curloc[X], curloc[Y],
1034  curloc[Z], aray[X], aray[Y], aray[Z]);
1035  bu_log("hf: from=(%g, %g) to=(%g, %g)\n",
1036  goesIN[X]/xWidth,
1037  goesIN[Y]/yWidth,
1038  goesOUT[X]/xWidth,
1039  goesOUT[Y]/yWidth);
1040  }
1041  error = curloc[Y]-yCell*yWidth;
1042  error /= yWidth;
1043 
1044  delta = aray[Y]/fabs(aray[X]);
1045 
1046  if (delta < 0.0) {
1047  delta = -delta;
1048  error = -error;
1049  } else {
1050  error -= 1.0;
1051  }
1052 
1053  /*
1054  * if at the far end (Xlen) then we need to move one step
1055  * forward along aray.
1056  */
1057  if (xCell >= hf->hf_w-1) {
1058  xCell+=signX;
1059  error += delta;
1060  }
1061  if (RT_G_DEBUG & DEBUG_HF) {
1062  bu_log("hf: delta=%g, error=%g, %d, %d\n",
1063  delta, error, xCell, yCell);
1064  }
1065 
1066 
1067  deltaZ = (aray[Z] < 0) ? -aray[Z] : aray[Z];
1068  do {
1069  farZ = curloc[Z] + aray[Z];
1070  maxZ = (curloc[Z] > farZ) ? curloc[Z] : farZ;
1071  minZ = (curloc[Z] < farZ) ? curloc[Z] : farZ;
1072  if (RT_G_DEBUG & DEBUG_HF) {
1073  bu_log("hf: cell %d, %d [%g -- %g]",
1074  xCell, yCell, minZ, maxZ);
1075  }
1076  /*
1077  * Are we on the grid yet? If not, then we will check for
1078  * a side step and inc.
1079  */
1080  /* CTJ - Or maxZ < hf->hf_min then no chance to hit */
1081  if (yCell < 0 || yCell > hf->hf_n-2) {
1082  if (error > -SQRT_SMALL_FASTF) {
1083  if (yCell >= -1) goto skip_first;
1084  yCell += signY;
1085  error -= 1.0;
1086  }
1087  xCell += signX;
1088  error += delta;
1089  VADD2(curloc, curloc, aray);
1090  continue;
1091  }
1092 
1093  /*
1094  * Get the min/max of the four corners of a given cell.
1095  * Since the data in memory can be in one of two forms,
1096  * unsigned short and double, we have this simple if
1097  * statement around the find min/max to reference the data
1098  * correctly.
1099  */
1100  if (hf->hf_shorts) {
1101  unsigned short *sp;
1102  sp = (unsigned short *)hf->hf_mp->apbuf +
1103  yCell * hf->hf_w + xCell;
1104  /* 0, 0 */
1105  highest = lowest = *sp++;
1106  /* 1, 0 */
1107  if (lowest > (double)*sp) lowest=*sp;
1108  if (highest < (double)*sp) highest=*sp;
1109  sp+=hf->hf_w;
1110  /* 1, 1 */
1111  if (lowest > (double)*sp) lowest=*sp;
1112  if (highest < (double)*sp) highest=*sp;
1113  sp--;
1114  /* 0, 1 */
1115  if (lowest > (double)*sp) lowest = *sp;
1116  if (highest < (double)*sp) highest = *sp;
1117  lowest *= hf->hf_file2mm;
1118  highest *= hf->hf_file2mm;
1119  } else {
1120  double *fp;
1121  fp = (double *)hf->hf_mp->apbuf +
1122  yCell * hf->hf_w + xCell;
1123  /* 0, 0 */
1124  highest = lowest = *fp++;
1125  /* 1, 0 */
1126  if (lowest > *fp) lowest=*fp;
1127  if (highest < *fp) highest=*fp;
1128  fp+=hf->hf_w;
1129  /* 1, 1 */
1130  if (lowest > *fp) lowest=*fp;
1131  if (highest < *fp) highest=*fp;
1132  fp--;
1133  /* 0, 1 */
1134  if (lowest > *fp) lowest = *fp;
1135  if (highest < *fp) highest = *fp;
1136  lowest *= hf->hf_file2mm;
1137  highest *= hf->hf_file2mm;
1138  }
1139 
1140  if (RT_G_DEBUG & DEBUG_HF) {
1141  bu_log("lowest=%g, highest=%g\n",
1142  lowest, highest);
1143  }
1144 
1145  /* This is the primary test. It is designed to get all
1146  * cells that the ray passes through.
1147  */
1148  if (maxZ+deltaZ > lowest &&
1149  minZ-deltaZ < highest) {
1150  int r = hf_cell_shot(stp, rp, ap, hp, xCell, yCell);
1151  if (r) {
1152  nhits += r;
1153  if (nhits >= MAXHITS)
1154  bu_bomb("g_hf.c: too many hits.\n");
1155  hp += r;
1156  }
1157  }
1158  /* This is the DDA trying to fill in the corners as it
1159  * walks the path.
1160  */
1161  skip_first:
1162  if (error > SQRT_SMALL_FASTF) {
1163  yCell += signY;
1164  if (RT_G_DEBUG & DEBUG_HF) {
1165  bu_log("hf: cell %d, %d ", xCell, yCell);
1166  }
1167  if ((yCell < 0) || yCell > hf->hf_n-2) {
1168  error -= 1.0;
1169  xCell += signX;
1170  error += delta;
1171  VADD2(curloc, curloc, aray);
1172  continue;
1173  }
1174  if (hf->hf_shorts) {
1175  unsigned short *sp;
1176  sp = (unsigned short *)hf->hf_mp->apbuf +
1177  yCell * hf->hf_w + xCell;
1178  /* 0, 0 */
1179  highest = lowest = *sp++;
1180  /* 1, 0 */
1181  if (lowest > (double)*sp) lowest=*sp;
1182  if (highest < (double)*sp) highest=*sp;
1183  sp+=hf->hf_w;
1184  /* 1, 1 */
1185  if (lowest > (double)*sp) lowest=*sp;
1186  if (highest < (double)*sp) highest=*sp;
1187  sp--;
1188  /* 0, 1 */
1189  if (lowest > (double)*sp) lowest = *sp;
1190  if (highest < (double)*sp) highest = *sp;
1191  lowest *= hf->hf_file2mm;
1192  highest *= hf->hf_file2mm;
1193  } else {
1194  double *fp;
1195  fp = (double *)hf->hf_mp->apbuf +
1196  yCell * hf->hf_w + xCell;
1197  /* 0, 0 */
1198  highest = lowest = *fp++;
1199  /* 1, 0 */
1200  if (lowest > *fp) lowest=*fp;
1201  if (highest < *fp) highest=*fp;
1202  fp+=hf->hf_w;
1203  /* 1, 1 */
1204  if (lowest > *fp) lowest=*fp;
1205  if (highest < *fp) highest=*fp;
1206  fp--;
1207  /* 0, 1 */
1208  if (lowest > *fp) lowest = *fp;
1209  if (highest < *fp) highest = *fp;
1210  lowest *= hf->hf_file2mm;
1211  highest *= hf->hf_file2mm;
1212  }
1213  if (maxZ+deltaZ > lowest &&
1214  minZ-deltaZ < highest) {
1215  int r = hf_cell_shot(stp, rp, ap, hp, xCell, yCell);
1216  /* DO HIT */
1217  if (r) {
1218  nhits += r;
1219  if (nhits >= MAXHITS)
1220  bu_bomb("g_hf.c: too many hits.\n");
1221  hp+=r;
1222  }
1223  }
1224  error -= 1.0;
1225  } else if (error > -SQRT_SMALL_FASTF) {
1226  yCell += signY;
1227  error -= 1.0;
1228  }
1229  xCell += signX;
1230  error += delta;
1231  VADD2(curloc, curloc, aray);
1232  } while (xCell >= 0 && xCell < hf->hf_w-1);
1233  if (RT_G_DEBUG & DEBUG_HF) {
1234  bu_log("htf: leaving loop, %d, %d, %g vs. 0--%d, 0--%d, 0.0--%g\n",
1235  xCell, yCell, curloc[Z], hf->hf_w-1, hf->hf_n-1, hf->hf_max);
1236  }
1237 
1238  } else {
1239  /* OTHER HALF */
1240 
1241  double tmp;
1242  double farZ, minZ, maxZ;
1243  double deltaZ;
1244  int xCell, yCell, signX, signY;
1245  double highest, lowest, error, delta;
1246 
1247  vect_t goesIN, goesOUT;
1248 
1249  VJOIN1(goesIN, rp->r_pt, allDist[4], rp->r_dir);
1250  VJOIN1(goesOUT, rp->r_pt, allDist[1], rp->r_dir);
1251  VSUB2(aray, goesOUT, goesIN);
1252  VSUB2(curloc, goesIN, hf->hf_V);
1253 
1254 
1255  /*
1256  * We will be stepping one "cell width" in the X direction
1257  * each time through the loop. In simple case we have???
1258  * cell_width = htfp->htf_Xlen/(htfp->htf_i-1);
1259  * deltaX = (Xdist*cell_width)/Xdist
1260  * deltaY = (Ydist*cell_width)/Xdist;
1261  */
1262  tmp = yWidth/fabs(aray[Y]);
1263 
1264  VSCALE(aray, aray, tmp);
1265 
1266  /*
1267  * Some fudges here. First, the size of the array of samples
1268  * is iXj, but the size of the array of CELLS is (i-1)X(j-1),
1269  * therefore number of CELLS is one less than number of
1270  * samples and the scaling factor is used. Second, the math
1271  * is nice to us. IF we are entering at the far end
1272  * (curloc[X] == Xlen || curloc[Y] == Ylen) then the result we
1273  * will get back is of the cell following this (out of bounds)
1274  * So we add a check for that problem.
1275  */
1276  yCell = curloc[Y]/yWidth;
1277  tmp = curloc[X];
1278  if (tmp < 0) {
1279  xCell = tmp/xWidth - 1;
1280  } else {
1281  xCell = tmp/xWidth;
1282  }
1283 
1284  signX = (aray[X] < 0.0) ? -1 : 1;
1285  signY = (aray[Y] < 0.0) ? -1 : 1;
1286 
1287  if (RT_G_DEBUG & DEBUG_HF) {
1288  bu_log("hf: curloc=(%g, %g, %g) aray=(%g, %g, %g)\n", curloc[X], curloc[Y],
1289  curloc[Z], aray[X], aray[Y], aray[Z]);
1290  bu_log("hf: from=(%g, %g) to=(%g, %g)\n",
1291  goesIN[X]/xWidth,
1292  goesIN[Y]/yWidth,
1293  goesOUT[X]/xWidth,
1294  goesOUT[Y]/yWidth);
1295  }
1296  error = curloc[X]-xCell*xWidth;
1297  error /= xWidth;
1298 
1299 /* delta = aray[X]/xWidth; */
1300  delta = aray[X]/fabs(aray[Y]);
1301 
1302  if (delta < 0.0) {
1303  delta = -delta;
1304  error = -error;
1305  } else {
1306  error -= 1.0;
1307  }
1308 
1309  /*
1310  * if at the far end (Ylen) then we need to move one step
1311  * forward along aray.
1312  */
1313  if (yCell >= hf->hf_n-1) {
1314  yCell+=signY;
1315  error += delta;
1316  }
1317  if (RT_G_DEBUG & DEBUG_HF) {
1318  bu_log("hf: delta=%g, error=%g, %d, %d\n",
1319  delta, error, xCell, yCell);
1320  }
1321 
1322  deltaZ = (aray[Z] < 0) ? -aray[Z] : aray[Z];
1323  do {
1324  farZ = curloc[Z] + aray[Z];
1325  maxZ = (curloc[Z] > farZ) ? curloc[Z] : farZ;
1326  minZ = (curloc[Z] < farZ) ? curloc[Z] : farZ;
1327  if (RT_G_DEBUG & DEBUG_HF) {
1328  bu_log("hf: cell %d, %d [%g -- %g] ",
1329  xCell, yCell, minZ, maxZ);
1330  }
1331  /* CTJ - Or maxZ < hf->hf_min */
1332  if (xCell < 0 || xCell > hf->hf_w-2) {
1333  if (error > -SQRT_SMALL_FASTF) {
1334  if (xCell >= -1) goto skip_2nd;
1335  xCell += signX;
1336  error -= 1.0;
1337  }
1338  yCell += signY;
1339  error += delta;
1340  VADD2(curloc, curloc, aray);
1341  continue;
1342  }
1343 
1344  if (hf->hf_shorts) {
1345  unsigned short *sp;
1346  sp = (unsigned short *)hf->hf_mp->apbuf +
1347  yCell * hf->hf_w + xCell;
1348  /* 0, 0 */
1349  highest = lowest = *sp++;
1350  /* 1, 0 */
1351  if (lowest > (double)*sp) lowest=*sp;
1352  if (highest < (double)*sp) highest=*sp;
1353  sp+=hf->hf_w;
1354  /* 1, 1 */
1355  if (lowest > (double)*sp) lowest=*sp;
1356  if (highest < (double)*sp) highest=*sp;
1357  sp--;
1358  /* 0, 1 */
1359  if (lowest > (double)*sp) lowest = *sp;
1360  if (highest < (double)*sp) highest = *sp;
1361  lowest *= hf->hf_file2mm;
1362  highest *= hf->hf_file2mm;
1363  } else {
1364  double *fp;
1365  fp = (double *)hf->hf_mp->apbuf +
1366  yCell * hf->hf_w + xCell;
1367  /* 0, 0 */
1368  highest = lowest = *fp++;
1369  /* 1, 0 */
1370  if (lowest > *fp) lowest=*fp;
1371  if (highest < *fp) highest=*fp;
1372  fp+=hf->hf_w;
1373  /* 1, 1 */
1374  if (lowest > *fp) lowest=*fp;
1375  if (highest < *fp) highest=*fp;
1376  fp--;
1377  /* 0, 1 */
1378  if (lowest > *fp) lowest = *fp;
1379  if (highest < *fp) highest = *fp;
1380  lowest *= hf->hf_file2mm;
1381  highest *= hf->hf_file2mm;
1382  }
1383 
1384 
1385  if (RT_G_DEBUG & DEBUG_HF) {
1386  bu_log("lowest=%g, highest=%g\n",
1387  lowest, highest);
1388  }
1389 
1390  /* This is the primary test. It is designed to get all
1391  * cells that the ray passes through.
1392  */
1393  if (maxZ+deltaZ > lowest &&
1394  minZ-deltaZ < highest) {
1395  int r = hf_cell_shot(stp, rp, ap, hp, xCell, yCell);
1396  if (r) {
1397  nhits += r;
1398  if (nhits >= MAXHITS)
1399  bu_bomb("g_hf.c: too many hits.\n");
1400  hp+=r;
1401  }
1402  }
1403  /* This is the DDA trying to fill in the corners as it
1404  * walks the path.
1405  */
1406  skip_2nd:
1407  if (error > SQRT_SMALL_FASTF) {
1408  xCell += signX;
1409  if (RT_G_DEBUG & DEBUG_HF) {
1410  bu_log("hf: cell %d, %d\n", xCell, yCell);
1411  }
1412  if ((xCell < 0) || xCell > hf->hf_w-2) {
1413  error -= 1.0;
1414  yCell += signY;
1415  error += delta;
1416  VADD2(curloc, curloc, aray);
1417  continue;
1418  }
1419  if (hf->hf_shorts) {
1420  unsigned short *sp;
1421  sp = (unsigned short *)hf->hf_mp->apbuf +
1422  yCell * hf->hf_w + xCell;
1423  /* 0, 0 */
1424  highest = lowest = *sp++;
1425  /* 1, 0 */
1426  if (lowest > (double)*sp) lowest=*sp;
1427  if (highest < (double)*sp) highest=*sp;
1428  sp+=hf->hf_w;
1429  /* 1, 1 */
1430  if (lowest > (double)*sp) lowest=*sp;
1431  if (highest < (double)*sp) highest=*sp;
1432  sp--;
1433  /* 0, 1 */
1434  if (lowest > (double)*sp) lowest = *sp;
1435  if (highest < (double)*sp) highest = *sp;
1436  lowest *= hf->hf_file2mm;
1437  highest *= hf->hf_file2mm;
1438  } else {
1439  double *fp;
1440  fp = (double *)hf->hf_mp->apbuf +
1441  yCell * hf->hf_w + xCell;
1442  /* 0, 0 */
1443  highest = lowest = *fp++;
1444  /* 1, 0 */
1445  if (lowest > *fp) lowest=*fp;
1446  if (highest < *fp) highest=*fp;
1447  fp+=hf->hf_w;
1448  /* 1, 1 */
1449  if (lowest > *fp) lowest=*fp;
1450  if (highest < *fp) highest=*fp;
1451  fp--;
1452  /* 0, 1 */
1453  if (lowest > *fp) lowest = *fp;
1454  if (highest < *fp) highest = *fp;
1455  lowest *= hf->hf_file2mm;
1456  highest *= hf->hf_file2mm;
1457  }
1458  if (maxZ+deltaZ > lowest &&
1459  minZ-deltaZ < highest) {
1460  int r = hf_cell_shot(stp, rp, ap, hp, xCell, yCell);
1461  /* DO HIT */
1462  if (r) {
1463  nhits += r;
1464  if (nhits >= MAXHITS)
1465  bu_bomb("g_hf.c: too many hits.\n");
1466  hp += r;
1467  }
1468  }
1469  error -= 1.0;
1470  } else if (error > -SQRT_SMALL_FASTF) {
1471  xCell += signX;
1472  error -= 1.0;
1473  }
1474  yCell += signY;
1475  error += delta;
1476  VADD2(curloc, curloc, aray);
1477  } while (yCell >= 0 && yCell < hf->hf_n-1);
1478  if (RT_G_DEBUG & DEBUG_HF) {
1479  bu_log("htf: leaving loop, %d, %d, %g vs. 0--%d, 0--%d, 0.0--%g\n",
1480  xCell, yCell, curloc[Z], hf->hf_w-1, hf->hf_n-1, hf->hf_max);
1481  }
1482  }
1483 
1484  /* Sort hits, near to Far */
1485  {
1486  int i;
1487  struct hit tmp;
1488  for (i = 0; i < nhits-1; i++) {
1489  for (j = i+1; j < nhits; j++) {
1490  if (hits[i].hit_dist <= hits[j].hit_dist) continue;
1491  tmp = hits[j];
1492  hits[j]=hits[i];
1493  hits[i]=tmp;
1494  }
1495  }
1496  }
1497 
1498  if (nhits & 1) {
1499  int i;
1500  static int nerrors = 0;
1501 
1502  if (nhits >= MAXHITS) {
1503  bu_bomb("g_hf.c: too many hits.\n");
1504  }
1505  hits[nhits] = hits[nhits-1]; /* struct copy*/
1506  VREVERSE(hits[nhits].hit_normal, hits[nhits-1].hit_normal);
1507  nhits++;
1508  if (nerrors++ < 300) {
1509  bu_log("rt_hf_shot(%s): %d hit(s)@ %d %d: ", stp->st_name, nhits-1, ap->a_x, ap->a_y);
1510  for (i = 0; i < nhits; i++) {
1511  bu_log("%f(%d), ", hits[i].hit_dist, hits[i].hit_surfno);
1512  }
1513  bu_log("\n");
1514  }
1515  }
1516 
1517  /* nhits is even, build segments */
1518  {
1519  struct seg *segp;
1520  int i;
1521  for (i = 0; i < nhits; i += 2) {
1522  RT_GET_SEG(segp, ap->a_resource);
1523  segp->seg_stp = stp;
1524  segp->seg_in = hits[i];
1525  segp->seg_out= hits[i+1];
1526  BU_LIST_INSERT(&(seghead->l), &(segp->l));
1527  }
1528  }
1529  return nhits; /* hits or misses */
1530 }
1531 
1532 
1533 /**
1534  * Given ONE ray distance, return the normal and entry/exit point.
1535  */
1536 void
1537 rt_hf_norm(struct hit *hitp, struct soltab *stp, struct xray *rp)
1538 {
1539  struct hf_specific *hf =
1540  (struct hf_specific *)stp->st_specific;
1541  int j;
1542 
1543  j = hitp->hit_surfno;
1544 
1545  VJOIN1(hitp->hit_point, rp->r_pt, hitp->hit_dist, rp->r_dir);
1546  if (j >= 0) {
1547  /* Normals computed in rt_htf_shot, nothing to do here. */
1548  return;
1549  }
1550 
1551  switch (j) {
1552  case -1:
1553  VREVERSE(hitp->hit_normal, hf->hf_X);
1554  break;
1555  case -2:
1556  VREVERSE(hitp->hit_normal, hf->hf_Y);
1557  break;
1558  case -3:
1559  VREVERSE(hitp->hit_normal, hf->hf_N);
1560  break;
1561  case -4:
1562  VMOVE(hitp->hit_normal, hf->hf_X);
1563  break;
1564  case -5:
1565  VMOVE(hitp->hit_normal, hf->hf_Y);
1566  break;
1567  case -6:
1568  VMOVE(hitp->hit_normal, hf->hf_N);
1569  break;
1570  }
1571  VUNITIZE(hitp->hit_normal);
1572 
1573 }
1574 
1575 
1576 /**
1577  * Return the curvature of the hf.
1578  */
1579 void
1580 rt_hf_curve(struct curvature *cvp, struct hit *hitp, struct soltab *stp)
1581 {
1582  struct hf_specific *hf = (struct hf_specific *)stp->st_specific;
1583  if (!hf || !cvp || !hitp)
1584  return;
1585  RT_CK_HIT(hitp);
1586 
1587  cvp->crv_c1 = cvp->crv_c2 = 0;
1588 
1589  /* any tangent direction */
1590  bn_vec_ortho(cvp->crv_pdir, hitp->hit_normal);
1591 }
1592 
1593 
1594 /**
1595  * For a hit on the surface of an hf, return the (u, v) coordinates of
1596  * the hit point, 0 <= u, v <= 1.
1597  *
1598  * u = azimuth
1599  * v = elevation
1600  */
1601 void
1602 rt_hf_uv(struct application *ap, struct soltab *stp, struct hit *hitp, struct uvcoord *uvp)
1603 {
1604  struct hf_specific *hf;
1605  vect_t delta;
1606  fastf_t r = 0;
1607 
1608  if (ap) RT_CK_APPLICATION(ap);
1609  if (!hitp || !uvp || !stp)
1610  return;
1611  hf = (struct hf_specific *)stp->st_specific;
1612  if (hf == NULL)
1613  return;
1614  RT_CK_HIT(hitp);
1615 
1616  VSUB2(delta, hitp->hit_point, hf->hf_V);
1617  uvp->uv_u = delta[X] / hf->hf_Xlen;
1618  uvp->uv_v = delta[Y] / hf->hf_Ylen;
1619  r = 0.0;
1620  if (uvp->uv_u < 0.0) uvp->uv_u = 0.0;
1621  if (uvp->uv_u > 1.0) uvp->uv_u = 1.0;
1622  if (uvp->uv_v < 0.0) uvp->uv_v = 0.0;
1623  if (uvp->uv_v > 1.0) uvp->uv_v = 1.0;
1624 
1625  uvp->uv_du = r;
1626  uvp->uv_dv = r;
1627 }
1628 
1629 
1630 void
1631 rt_hf_free(struct soltab *stp)
1632 {
1633  struct hf_specific *hf =
1634  (struct hf_specific *)stp->st_specific;
1635 
1636  if (hf->hf_mp) {
1638  hf->hf_mp = (struct bu_mapped_file *)0;
1639  }
1640  BU_PUT(hf, struct hf_specific);
1641 }
1642 
1643 
1644 int
1645 rt_hf_plot(struct bu_list *vhead, struct rt_db_internal *ip, const struct rt_tess_tol *ttol, const struct bn_tol *UNUSED(tol), const struct rt_view_info *UNUSED(info))
1646 {
1647  struct rt_hf_internal *xip;
1648  unsigned short *sp = (unsigned short *)NULL;
1649  double *dp;
1650  vect_t xbasis;
1651  vect_t ybasis;
1652  vect_t zbasis;
1653  point_t start;
1654  point_t cur;
1655  size_t x;
1656  size_t y;
1657  int cmd;
1658  int step;
1659  int half_step;
1660  int goal;
1661 
1662  BU_CK_LIST_HEAD(vhead);
1663  RT_CK_DB_INTERNAL(ip);
1664  xip = (struct rt_hf_internal *)ip->idb_ptr;
1665  RT_HF_CK_MAGIC(xip);
1666 
1667  VSCALE(xbasis, xip->x, xip->xlen / (xip->w - 1));
1668  VSCALE(ybasis, xip->y, xip->ylen / (xip->n - 1));
1669  VCROSS(zbasis, xip->x, xip->y);
1670  VSCALE(zbasis, zbasis, xip->zscale * xip->file2mm);
1671 
1672  /* XXX This should be set from the tessellation tolerance */
1673  goal = 20000;
1674 
1675  /* Draw the 4 corners of the base plate */
1676  RT_ADD_VLIST(vhead, xip->v, BN_VLIST_LINE_MOVE);
1677 
1678  VJOIN1(start, xip->v, xip->xlen, xip->x);
1679  RT_ADD_VLIST(vhead, start, BN_VLIST_LINE_DRAW);
1680 
1681  VJOIN2(start, xip->v, xip->xlen, xip->x, xip->ylen, xip->y);
1682  RT_ADD_VLIST(vhead, start, BN_VLIST_LINE_DRAW);
1683 
1684  VJOIN1(start, xip->v, xip->ylen, xip->y);
1685  RT_ADD_VLIST(vhead, start, BN_VLIST_LINE_DRAW);
1686 
1687  RT_ADD_VLIST(vhead, xip->v, BN_VLIST_LINE_DRAW);
1688  goal -= 5;
1689 
1690 #define HF_GET(_p, _x, _y) ((_p)[(_y)*xip->w+(_x)])
1691  /*
1692  * Draw the four "ridge lines" at full resolution, for edge matching.
1693  */
1694  if (xip->shorts) {
1695  /* X direction, Y=0, with edges down to base */
1696  RT_ADD_VLIST(vhead, xip->v, BN_VLIST_LINE_MOVE);
1697  sp = &HF_GET((unsigned short *)xip->mp->apbuf, 0, 0);
1698  for (x = 0; x < xip->w; x++) {
1699  VJOIN2(cur, xip->v, x, xbasis, *sp, zbasis);
1700  RT_ADD_VLIST(vhead, cur, BN_VLIST_LINE_DRAW);
1701  sp++;
1702  }
1703  VJOIN1(cur, xip->v, xip->xlen, xip->x);
1704  RT_ADD_VLIST(vhead, cur, BN_VLIST_LINE_DRAW);
1705 
1706  /* X direction, Y=n-1, with edges down to base */
1707  VJOIN1(start, xip->v, xip->ylen, xip->y);
1708  RT_ADD_VLIST(vhead, start, BN_VLIST_LINE_MOVE);
1709  sp = &HF_GET((unsigned short *)xip->mp->apbuf, 0, xip->n - 1);
1710  VJOIN1(start, xip->v, xip->ylen, xip->y);
1711  for (x = 0; x < xip->w; x++) {
1712  VJOIN2(cur, start, x, xbasis, *sp, zbasis);
1713  RT_ADD_VLIST(vhead, cur, BN_VLIST_LINE_DRAW);
1714  sp++;
1715  }
1716  VJOIN2(cur, xip->v, xip->xlen, xip->x, xip->ylen, xip->y);
1717  RT_ADD_VLIST(vhead, cur, BN_VLIST_LINE_DRAW);
1718 
1719  /* Y direction, X=0 */
1720  cmd = BN_VLIST_LINE_MOVE;
1721  sp = &HF_GET((unsigned short *)xip->mp->apbuf, 0, 0);
1722  for (y = 0; y < xip->n; y++) {
1723  VJOIN2(cur, xip->v, y, ybasis, *sp, zbasis);
1724  RT_ADD_VLIST(vhead, cur, cmd);
1725  cmd = BN_VLIST_LINE_DRAW;
1726  sp += xip->w;
1727  }
1728 
1729  /* Y direction, X=w-1 */
1730  cmd = BN_VLIST_LINE_MOVE;
1731  sp = &HF_GET((unsigned short *)xip->mp->apbuf, xip->w - 1, 0);
1732  VJOIN1(start, xip->v, xip->xlen, xip->x);
1733  for (y = 0; y < xip->n; y++) {
1734  VJOIN2(cur, start, y, ybasis, *sp, zbasis);
1735  RT_ADD_VLIST(vhead, cur, cmd);
1736  cmd = BN_VLIST_LINE_DRAW;
1737  sp += xip->w;
1738  }
1739  } else {
1740  /* X direction, Y=0, with edges down to base */
1741  RT_ADD_VLIST(vhead, xip->v, BN_VLIST_LINE_MOVE);
1742  dp = &HF_GET((double *)xip->mp->apbuf, 0, 0);
1743  for (x = 0; x < xip->w; x++) {
1744  VJOIN2(cur, xip->v, x, xbasis, *dp, zbasis);
1745  RT_ADD_VLIST(vhead, cur, BN_VLIST_LINE_DRAW);
1746  dp++;
1747  }
1748  VJOIN1(cur, xip->v, xip->xlen, xip->x);
1749  RT_ADD_VLIST(vhead, cur, BN_VLIST_LINE_DRAW);
1750 
1751  /* X direction, Y=n-1, with edges down to base */
1752  VJOIN1(start, xip->v, xip->ylen, xip->y);
1753  RT_ADD_VLIST(vhead, start, BN_VLIST_LINE_MOVE);
1754  dp = &HF_GET((double *)xip->mp->apbuf, 0, xip->n - 1);
1755  VJOIN1(start, xip->v, xip->ylen, xip->y);
1756  for (x = 0; x < xip->w; x++) {
1757  VJOIN2(cur, start, x, xbasis, *dp, zbasis);
1758  RT_ADD_VLIST(vhead, cur, BN_VLIST_LINE_DRAW);
1759  dp++;
1760  }
1761  VJOIN2(cur, xip->v, xip->xlen, xip->x, xip->ylen, xip->y);
1762  RT_ADD_VLIST(vhead, cur, BN_VLIST_LINE_DRAW);
1763 
1764  /* Y direction, X=0 */
1765  cmd = BN_VLIST_LINE_MOVE;
1766  dp = &HF_GET((double *)xip->mp->apbuf, 0, 0);
1767  for (y = 0; y < xip->n; y++) {
1768  VJOIN2(cur, xip->v, y, ybasis, *dp, zbasis);
1769  RT_ADD_VLIST(vhead, cur, cmd);
1770  cmd = BN_VLIST_LINE_DRAW;
1771  sp += xip->w;
1772  }
1773 
1774  /* Y direction, X=w-1 */
1775  cmd = BN_VLIST_LINE_MOVE;
1776  dp = &HF_GET((double *)xip->mp->apbuf, xip->w - 1, 0);
1777  VJOIN1(start, xip->v, xip->xlen, xip->x);
1778  for (y = 0; y < xip->n; y++) {
1779  VJOIN2(cur, start, y, ybasis, *dp, zbasis);
1780  RT_ADD_VLIST(vhead, cur, cmd);
1781  cmd = BN_VLIST_LINE_DRAW;
1782  dp += xip->w;
1783  }
1784  }
1785  goal -= 4 + 2 * (xip->w + xip->n);
1786 
1787  /* Apply relative tolerance, if specified */
1788  if (!ZERO(ttol->rel)) {
1789  size_t rstep;
1790  rstep = xip->w;
1791  V_MAX(rstep, xip->n);
1792  step = ttol->rel * rstep;
1793  } else {
1794  /* No relative tol specified, limit drawing to 'goal' # of vectors */
1795  if (goal <= 0) return 0; /* no vectors for interior */
1796 
1797  /* Compute data stride based upon producing no more than 'goal' vectors */
1798  step = ceil(sqrt(2*(xip->w-1)*(xip->n-1) / (double)goal));
1799  }
1800  if (step < 1) step = 1;
1801  if ((half_step = step/2) < 1) half_step = 1;
1802 
1803  /* Draw the contour lines in W (x) direction. Don't redo ridges. */
1804  for (y = half_step; y < xip->n-half_step; y += step) {
1805  VJOIN1(start, xip->v, y, ybasis);
1806  x = 0;
1807  if (xip->shorts) {
1808  sp = &HF_GET((unsigned short *)xip->mp->apbuf, x, y);
1809  VJOIN2(cur, start, x, xbasis, *sp, zbasis);
1810  RT_ADD_VLIST(vhead, cur, BN_VLIST_LINE_MOVE);
1811  x += half_step;
1812  sp = &HF_GET((unsigned short *)xip->mp->apbuf, x, y);
1813  for (; x < xip->w; x += step) {
1814  VJOIN2(cur, start, x, xbasis, *sp, zbasis);
1815  RT_ADD_VLIST(vhead, cur, BN_VLIST_LINE_DRAW);
1816  sp += step;
1817  }
1818  if (x != step+xip->w-1+step) {
1819  x = xip->w - 1;
1820  sp = &HF_GET((unsigned short *)xip->mp->apbuf, x, y);
1821  VJOIN2(cur, start, x, xbasis, *sp, zbasis);
1822  RT_ADD_VLIST(vhead, cur, BN_VLIST_LINE_DRAW);
1823  }
1824  } else {
1825  /* doubles */
1826  dp = &HF_GET((double *)xip->mp->apbuf, x, y);
1827  VJOIN2(cur, start, x, xbasis, *dp, zbasis);
1828  RT_ADD_VLIST(vhead, cur, BN_VLIST_LINE_MOVE);
1829  x += half_step;
1830  dp = &HF_GET((double *)xip->mp->apbuf, x, y);
1831  for (; x < xip->w; x+=step) {
1832  VJOIN2(cur, start, x, xbasis, *dp, zbasis);
1833  RT_ADD_VLIST(vhead, cur, BN_VLIST_LINE_DRAW);
1834  dp += step;
1835  }
1836  if (x != step+xip->w-1+step) {
1837  x = xip->w - 1;
1838  dp = &HF_GET((double *)xip->mp->apbuf, x, y);
1839  VJOIN2(cur, start, x, xbasis, *dp, zbasis);
1840  RT_ADD_VLIST(vhead, cur, BN_VLIST_LINE_DRAW);
1841  }
1842  }
1843  }
1844 
1845  /* Draw the contour lines in the N (y) direction */
1846  if (xip->shorts) {
1847  for (x = half_step; x < xip->w-half_step; x += step) {
1848  VJOIN1(start, xip->v, x, xbasis);
1849  y = 0;
1850  sp = &HF_GET((unsigned short *)xip->mp->apbuf, x, y);
1851  VJOIN2(cur, start, y, ybasis, *sp, zbasis);
1852  RT_ADD_VLIST(vhead, cur, BN_VLIST_LINE_MOVE);
1853  y += half_step;
1854  for (; y < xip->n; y += step) {
1855  sp = &HF_GET((unsigned short *)xip->mp->apbuf, x, y);
1856  VJOIN2(cur, start, y, ybasis, *sp, zbasis);
1857  RT_ADD_VLIST(vhead, cur, BN_VLIST_LINE_DRAW);
1858  }
1859  if (y != step+xip->n-1+step) {
1860  y = xip->n - 1;
1861  sp = &HF_GET((unsigned short *)xip->mp->apbuf, x, y);
1862  VJOIN2(cur, start, y, ybasis, *sp, zbasis);
1863  RT_ADD_VLIST(vhead, cur, BN_VLIST_LINE_DRAW);
1864  }
1865  }
1866  } else {
1867  /* doubles */
1868  for (x = half_step; x < xip->w-half_step; x += step) {
1869  VJOIN1(start, xip->v, x, xbasis);
1870  y = 0;
1871  dp = &HF_GET((double *)xip->mp->apbuf, x, y);
1872  VJOIN2(cur, start, y, ybasis, *dp, zbasis);
1873  RT_ADD_VLIST(vhead, cur, BN_VLIST_LINE_MOVE);
1874  y += half_step;
1875  for (; y < xip->n; y += step) {
1876  dp = &HF_GET((double *)xip->mp->apbuf, x, y);
1877  VJOIN2(cur, start, y, ybasis, *dp, zbasis);
1878  RT_ADD_VLIST(vhead, cur, BN_VLIST_LINE_DRAW);
1879  }
1880  if (y != step+xip->n-1+step) {
1881  y = xip->n - 1;
1882  dp = &HF_GET((double *)xip->mp->apbuf, x, y);
1883  VJOIN2(cur, start, y, ybasis, *dp, zbasis);
1884  RT_ADD_VLIST(vhead, cur, BN_VLIST_LINE_DRAW);
1885  }
1886  }
1887  }
1888  return 0;
1889 }
1890 
1891 
1892 /**
1893  * Returns -
1894  * -1 failure
1895  * 0 OK. *r points to nmgregion that holds this tessellation.
1896  */
1897 int
1898 rt_hf_tess(struct nmgregion **r, struct model *m, struct rt_db_internal *ip, const struct rt_tess_tol *UNUSED(ttol), const struct bn_tol *UNUSED(tol))
1899 {
1900  struct rt_hf_internal *xip;
1901 
1902  RT_CK_DB_INTERNAL(ip);
1903  xip = (struct rt_hf_internal *)ip->idb_ptr;
1904  RT_HF_CK_MAGIC(xip);
1905 
1906  if (r) *r = NULL;
1907  if (m) NMG_CK_MODEL(m);
1908 
1909  return -1;
1910 }
1911 
1912 
1913 /**
1914  * Import an HF from the database format to the internal format.
1915  * Apply modeling transformations as well.
1916  */
1917 int
1918 rt_hf_import4(struct rt_db_internal *ip, const struct bu_external *ep, const fastf_t *mat, const struct db_i *dbip)
1919 {
1920  struct rt_hf_internal *xip;
1921  union record *rp;
1922  struct bu_vls str = BU_VLS_INIT_ZERO;
1923  struct bu_mapped_file *mp;
1924  vect_t tmp;
1925  int in_cookie; /* format cookie */
1926  int in_len;
1927  int out_cookie;
1928  size_t count;
1929  size_t got;
1930 
1931  if (dbip) RT_CK_DBI(dbip);
1932 
1933  BU_CK_EXTERNAL(ep);
1934  rp = (union record *)ep->ext_buf;
1935  /* Check record type */
1936  if (rp->u_id != DBID_STRSOL) {
1937  bu_log("rt_hf_import4: defective record\n");
1938  return -1;
1939  }
1940 
1941  RT_CK_DB_INTERNAL(ip);
1942  ip->idb_major_type = DB5_MAJORTYPE_BRLCAD;
1943  ip->idb_type = ID_HF;
1944  ip->idb_meth = &OBJ[ID_HF];
1945  BU_ALLOC(ip->idb_ptr, struct rt_hf_internal);
1946 
1947  xip = (struct rt_hf_internal *)ip->idb_ptr;
1948  xip->magic = RT_HF_INTERNAL_MAGIC;
1949 
1950  /* Provide defaults. Only non-defaulted fields are dfile, w, n */
1951  xip->shorts = 1; /* for now */
1952  xip->file2mm = 1.0;
1953  VSETALL(xip->v, 0);
1954  VSET(xip->x, 1, 0, 0);
1955  VSET(xip->y, 0, 1, 0);
1956  xip->xlen = 1000;
1957  xip->ylen = 1000;
1958  xip->zscale = 1;
1959  bu_strlcpy(xip->fmt, "nd", sizeof(xip->fmt));
1960 
1961  /* Process parameters found in .g file */
1962  bu_vls_strcpy(&str, rp->ss.ss_args);
1963  if (bu_struct_parse(&str, rt_hf_parse, (char *)xip, NULL) < 0) {
1964  bu_vls_free(&str);
1965  err1:
1966  bu_free((char *)xip, "rt_hf_import4: xip");
1967  ip->idb_type = ID_NULL;
1968  ip->idb_ptr = (void *)NULL;
1969  return -2;
1970  }
1971  bu_vls_free(&str);
1972 
1973  /* If "cfile" was specified, process parameters from there */
1974  if (xip->cfile[0]) {
1975  FILE *fp;
1976 
1978  fp = fopen(xip->cfile, "rb");
1980  if (!fp) {
1981  perror(xip->cfile);
1982  bu_log("rt_hf_import4() unable to open cfile=%s\n", xip->cfile);
1983  goto err1;
1984  }
1985  while (bu_vls_gets(&str, fp) >= 0)
1986  bu_vls_strcat(&str, " ");
1988  fclose(fp);
1990  if (bu_struct_parse(&str, rt_hf_cparse, (char *)xip, NULL) < 0) {
1991  bu_log("rt_hf_import4() parse error in cfile input '%s'\n",
1992  bu_vls_addr(&str));
1993  bu_vls_free(&str);
1994  goto err1;
1995  }
1996  }
1997 
1998  /* Check for reasonable values */
1999  if (!xip->dfile[0]) {
2000  /* XXX Should create 2x2 data file instead, for positioning use (FPO) */
2001  bu_log("rt_hf_import4() no dfile specified\n");
2002  goto err1;
2003  }
2004  if (xip->w < 2 || xip->n < 2) {
2005  bu_log("rt_hf_import4() w=%u, n=%u too small\n", xip->w, xip->n);
2006  goto err1;
2007  }
2008  if (xip->xlen <= 0 || xip->ylen <= 0) {
2009  bu_log("rt_hf_import4() xlen=%g, ylen=%g too small\n", xip->xlen, xip->ylen);
2010  goto err1;
2011  }
2012 
2013  /* Apply modeling transformations */
2014  if (mat == NULL) mat = bn_mat_identity;
2015  MAT4X3PNT(tmp, mat, xip->v);
2016  VMOVE(xip->v, tmp);
2017  MAT4X3VEC(tmp, mat, xip->x);
2018  VMOVE(xip->x, tmp);
2019  MAT4X3VEC(tmp, mat, xip->y);
2020  VMOVE(xip->y, tmp);
2021  xip->xlen /= mat[15];
2022  xip->ylen /= mat[15];
2023  xip->zscale /= mat[15];
2024 
2025  VUNITIZE(xip->x);
2026  VUNITIZE(xip->y);
2027 
2028  /* Prepare for cracking input file format */
2029  if ((in_cookie = bu_cv_cookie(xip->fmt)) == 0) {
2030  bu_log("rt_hf_import4() fmt='%s' unknown\n", xip->fmt);
2031  goto err1;
2032  }
2033  in_len = bu_cv_itemlen(in_cookie);
2034 
2035  /*
2036  * Load data file, and transform to internal format
2037  */
2038  mp = bu_open_mapped_file(xip->dfile, "hf");
2039  if (!mp) {
2040  bu_log("rt_hf_import4() unable to open '%s'\n", xip->dfile);
2041  goto err1;
2042  }
2043  xip->mp = mp;
2044  count = mp->buflen / in_len;
2045 
2046  /* If this data has already been mapped, all done */
2047  if (mp->apbuf) return 0; /* OK */
2048 
2049  /* Transform external data to internal format -- short or double */
2050  if (xip->shorts) {
2051  mp->apbuflen = sizeof(unsigned short) * count;
2052  out_cookie = bu_cv_cookie("hus");
2053  } else {
2054  mp->apbuflen = sizeof(double) * count;
2055  out_cookie = bu_cv_cookie("hd");
2056  }
2057 
2058  if (bu_cv_optimize(in_cookie) == bu_cv_optimize(out_cookie)) {
2059  /* Don't replicate the data, just re-use the pointer */
2060  mp->apbuf = mp->buf;
2061  return 0; /* OK */
2062  }
2063 
2064  mp->apbuf = (void *)bu_malloc(mp->apbuflen, "rt_hf_import4 apbuf[]");
2065  got = bu_cv_w_cookie(mp->apbuf, out_cookie, mp->apbuflen,
2066  mp->buf, in_cookie, count);
2067  if (got != count) {
2068  bu_log("rt_hf_import4(%s) bu_cv_w_cookie count=%zu, got=%zu\n",
2069  xip->dfile, count, got);
2070  }
2071 
2072  return 0; /* OK */
2073 }
2074 
2075 
2076 /**
2077  * The name is added by the caller, in the usual place.
2078  *
2079  * The meaning of the export here is slightly different than that of
2080  * most other solids. The cfile and dfile are not modified, only
2081  * changes to the string solid parameters are placed back into the .g
2082  * file. Note that any parameters taken from a cfile are included in
2083  * the new string solid. This isn't a problem, because if the cfile
2084  * is changed (perhaps to substitute a different resolution height
2085  * field of the same location in space), its new parameters will
2086  * override those stored in the string solid (including the dfile
2087  * name).
2088  */
2089 int
2090 rt_hf_export4(struct bu_external *ep, const struct rt_db_internal *ip, double local2mm, const struct db_i *dbip)
2091 {
2092  struct rt_hf_internal *xip;
2093  union record *rec;
2094  struct bu_vls str = BU_VLS_INIT_ZERO;
2095 
2096  if (dbip) RT_CK_DBI(dbip);
2097 
2098  RT_CK_DB_INTERNAL(ip);
2099  if (ip->idb_type != ID_HF) return -1;
2100  xip = (struct rt_hf_internal *)ip->idb_ptr;
2101  RT_HF_CK_MAGIC(xip);
2102 
2103  /* Apply any scale transformation */
2104  xip->xlen /= local2mm;
2105  xip->ylen /= local2mm;
2106  xip->zscale /= local2mm;
2107 
2108  BU_CK_EXTERNAL(ep);
2109  ep->ext_nbytes = sizeof(union record) * DB_SS_NGRAN;
2110  ep->ext_buf = (uint8_t *)bu_calloc(1, ep->ext_nbytes, "hf external");
2111  rec = (union record *)ep->ext_buf;
2112 
2113  bu_vls_struct_print(&str, rt_hf_parse, (char *)xip);
2114 
2115  /* Any changes made by solid editing affect .g file only, and not
2116  * the cfile, if specified.
2117  */
2118 
2119  rec->s.s_id = DBID_STRSOL;
2120  bu_strlcpy(rec->ss.ss_keyword, "hf", sizeof(rec->ss.ss_keyword));
2121  bu_strlcpy(rec->ss.ss_args, bu_vls_addr(&str), DB_SS_LEN);
2122  bu_vls_free(&str);
2123 
2124  bu_log("DEPRECATED: The 'hf' height field primitive is no longer supported. Use the 'dsp' displacement map instead.\n");
2125 
2126  return 0;
2127 }
2128 
2129 
2130 int
2131 rt_hf_import5(struct rt_db_internal *ip, const struct bu_external *ep, const fastf_t *mat, const struct db_i *dbip)
2132 {
2133  if (ip) RT_CK_DB_INTERNAL(ip);
2134  if (!ep || !mat) return -1;
2135  if (dbip) RT_CK_DBI(dbip);
2136 
2137  bu_log("DEPRECATED: The 'hf' height field primitive is no longer supported. Use the 'dsp' displacement map instead.\n");
2138  bu_log("\tTo convert HF primitives to DSP, use 'dbupgrade'.\n");
2139  /* The rt_hf_to_dsp() routine can also be used */
2140  return -1;
2141 }
2142 
2143 
2144 int
2145 rt_hf_export5(struct bu_external *ep, const struct rt_db_internal *ip, double UNUSED(local2mm), const struct db_i *dbip)
2146 {
2147  if (!ep) return -1;
2148  if (ip) RT_CK_DB_INTERNAL(ip);
2149  if (dbip) RT_CK_DBI(dbip);
2150 
2151  bu_log("DEPRECATED: The 'hf' height field primitive is no longer supported. Use the 'dsp' displacement map instead.\n");
2152  bu_log("\tTo convert HF primitives to DSP, use 'dbupgrade'.\n");
2153  /* The rt_hf_to_dsp() routine can also be used */
2154  return -1;
2155 }
2156 
2157 
2158 /**
2159  * Make human-readable formatted presentation of this solid. First
2160  * line describes type of solid. Additional lines are indented one
2161  * tab, and give parameter values.
2162  */
2163 int
2164 rt_hf_describe(struct bu_vls *str, const struct rt_db_internal *ip, int verbose, double mm2local)
2165 {
2166  struct rt_hf_internal *xip;
2167 
2168  BU_CK_VLS(str);
2169 
2170  xip = (struct rt_hf_internal *)ip->idb_ptr;
2171  RT_HF_CK_MAGIC(xip);
2172 
2173  bu_vls_printf(str, "Height Field (HF) mm2local=%g\n", mm2local);
2174  if (!verbose)
2175  return 0;
2176 
2177  bu_vls_struct_print(str, rt_hf_parse, (const char *)ip->idb_ptr);
2178  bu_vls_strcat(str, "\n");
2179 
2180  return 0;
2181 }
2182 
2183 
2184 /**
2185  * Free the storage associated with the rt_db_internal version of this
2186  * solid.
2187  */
2188 void
2190 {
2191  struct rt_hf_internal *xip;
2192 
2193  RT_CK_DB_INTERNAL(ip);
2194  xip = (struct rt_hf_internal *)ip->idb_ptr;
2195  RT_HF_CK_MAGIC(xip);
2196  xip->magic = 0; /* sanity */
2197 
2198  if (xip->mp) {
2199  BU_CK_MAPPED_FILE(xip->mp);
2200  bu_close_mapped_file(xip->mp);
2201  }
2202 
2203  bu_free((char *)xip, "hf ifree");
2204  ip->idb_ptr = ((void *)0); /* sanity */
2205 }
2206 int
2207 rt_hf_params(struct pc_pc_set *UNUSED(ps), const struct rt_db_internal *ip)
2208 {
2209  if (ip) RT_CK_DB_INTERNAL(ip);
2210 
2211  return 0; /* OK */
2212 }
2213 
2214 
2215 /** @} */
2216 /*
2217  * Local Variables:
2218  * mode: C
2219  * tab-width: 8
2220  * indent-tabs-mode: t
2221  * c-file-style: "stroustrup"
2222  * End:
2223  * ex: shiftwidth=4 tabstop=8
2224  */
void bu_vls_init(struct bu_vls *vp)
Definition: vls.c:56
Definition: raytrace.h:800
uint32_t hit_magic
Definition: raytrace.h:249
void bu_log(const char *,...) _BU_ATTR_PRINTF12
Definition: log.c:176
#define BU_LIST_INSERT(old, new)
Definition: list.h:183
int rt_hf_import5(struct rt_db_internal *ip, const struct bu_external *ep, const fastf_t *mat, const struct db_i *dbip)
Definition: hf.c:2131
size_t bu_cv_itemlen(int cookie)
struct hit seg_in
IN information.
Definition: raytrace.h:370
#define RT_DSP_INTERNAL_MAGIC
Definition: magic.h:88
Definition: list.h:118
int hf_w
Definition: hf.c:105
void rt_hf_print(const struct soltab *stp)
Definition: hf.c:422
int rt_hf_shot(struct soltab *stp, struct xray *rp, struct application *ap, struct seg *seghead)
Definition: hf.c:785
#define RT_CK_APPLICATION(_p)
Definition: raytrace.h:1675
vect_t crv_pdir
Principle direction.
Definition: raytrace.h:307
const mat_t bn_mat_identity
Matrix and vector functionality.
Definition: mat.c:46
fastf_t uv_u
Range 0..1.
Definition: raytrace.h:341
struct soltab * seg_stp
pointer back to soltab
Definition: raytrace.h:372
if lu s
Definition: nmg_mod.c:3860
void bu_vls_strcat(struct bu_vls *vp, const char *s)
Definition: vls.c:368
#define VSET(a, b, c, d)
Definition: color.c:53
#define VSETALL(a, s)
Definition: color.c:54
Definition: raytrace.h:215
void bu_semaphore_acquire(unsigned int i)
Definition: semaphore.c:180
const struct bu_structparse rt_hf_cparse[]
Definition: hf.c:83
int rt_hf_to_dsp(struct rt_db_internal *db_intern)
Definition: hf.c:118
Definition: pc.h:108
Definition: raytrace.h:368
fastf_t hf_file2mm
Definition: hf.c:104
Definition: raytrace.h:248
#define SMALL_FASTF
Definition: defines.h:342
fastf_t st_aradius
Radius of APPROXIMATING sphere.
Definition: raytrace.h:433
#define ID_NULL
Unused.
Definition: raytrace.h:458
#define ID_HF
Height Field.
Definition: raytrace.h:482
Header file for the BRL-CAD common definitions.
#define HF_O(m)
Definition: hf.c:63
int rt_hf_describe(struct bu_vls *str, const struct rt_db_internal *ip, int verbose, double mm2local)
Definition: hf.c:2164
int rt_hf_plot(struct bu_list *vhead, struct rt_db_internal *ip, const struct rt_tess_tol *ttol, const struct bn_tol *tol, const struct rt_view_info *info)
Definition: hf.c:1645
struct resource * a_resource
dynamic memory resources
Definition: raytrace.h:1591
int rt_hf_params(struct pc_pc_set *ps, const struct rt_db_internal *ip)
Definition: hf.c:2207
fastf_t hf_max
Definition: hf.c:103
#define DEBUG_HF
27 Height Field solids
Definition: raytrace.h:114
fastf_t hf_Xlen
Definition: hf.c:98
void rt_hf_norm(struct hit *hitp, struct soltab *stp, struct xray *rp)
Definition: hf.c:1537
int hf_n
Definition: hf.c:106
vect_t hf_X
Definition: hf.c:97
struct bu_list l
Definition: raytrace.h:369
void * bu_malloc(size_t siz, const char *str)
Definition: malloc.c:314
void rt_hf_free(struct soltab *stp)
Definition: hf.c:1631
double rel
rel dist tol
Definition: raytrace.h:181
fastf_t hf_min
Definition: hf.c:102
if(share_geom)
Definition: nmg_mod.c:3829
int rt_hf_export5(struct bu_external *ep, const struct rt_db_internal *ip, double local2mm, const struct db_i *dbip)
Definition: hf.c:2145
size_t apbuflen
Definition: mapped_file.h:90
int idb_major_type
Definition: raytrace.h:192
int rt_hf_export4(struct bu_external *ep, const struct rt_db_internal *ip, double local2mm, const struct db_i *dbip)
Definition: hf.c:2090
void bu_vls_free(struct bu_vls *vp)
Definition: vls.c:248
Definition: color.c:49
#define BU_CK_VLS(_vp)
Definition: vls.h:69
#define ID_DSP
Displacement map.
Definition: raytrace.h:483
void * memset(void *s, int c, size_t n)
const struct bu_structparse rt_hf_parse[]
Definition: hf.c:66
#define RT_G_DEBUG
Definition: raytrace.h:1718
#define RT_ADD_VLIST(hd, pnt, draw)
Definition: raytrace.h:1865
size_t bu_cv_w_cookie(void *out, int outcookie, size_t size, void *in, int incookie, size_t count)
Definition: convert.c:474
#define RT_CK_DB_INTERNAL(_p)
Definition: raytrace.h:207
#define BU_ALLOC(_ptr, _type)
Definition: malloc.h:223
void * bu_calloc(size_t nelem, size_t elsize, const char *str)
Definition: malloc.c:321
vect_t hf_N
Definition: hf.c:101
fastf_t st_bradius
Radius of BOUNDING sphere.
Definition: raytrace.h:434
fastf_t crv_c2
curvature in other direction
Definition: raytrace.h:309
#define RT_CK_HIT(_p)
Definition: raytrace.h:259
#define BN_VLIST_LINE_MOVE
Definition: vlist.h:82
const struct rt_functab * idb_meth
for ft_ifree(), etc.
Definition: raytrace.h:194
struct bu_mapped_file * hf_mp
Definition: hf.c:108
fastf_t uv_dv
delta in v
Definition: raytrace.h:344
#define bu_strlcpy(dst, src, size)
Definition: str.h:60
void bn_mat_inv(mat_t output, const mat_t input)
int a_x
Screen X of ray, if applicable.
Definition: raytrace.h:1596
uint8_t * ext_buf
Definition: parse.h:216
#define BU_GET(_ptr, _type)
Definition: malloc.h:201
point_t hit_point
DEPRECATED: Intersection point, use VJOIN1 hit_dist.
Definition: raytrace.h:251
point_t st_max
max X, Y, Z of bounding RPP
Definition: raytrace.h:438
#define BU_CK_MAPPED_FILE(_mf)
Definition: mapped_file.h:101
#define SQRT_SMALL_FASTF
Definition: defines.h:346
struct hit seg_out
OUT information.
Definition: raytrace.h:371
#define BN_VLIST_LINE_DRAW
Definition: vlist.h:83
void rt_hf_curve(struct curvature *cvp, struct hit *hitp, struct soltab *stp)
Definition: hf.c:1580
int bu_cv_cookie(const char *in)
Definition: convert.c:33
#define UNUSED(parameter)
Definition: common.h:239
#define BU_PUT(_ptr, _type)
Definition: malloc.h:215
struct bu_mapped_file * bu_open_mapped_file(const char *name, const char *appl)
Definition: mappedfile.c:56
goto out
Definition: nmg_mod.c:3846
Support for uniform tolerances.
Definition: tol.h:71
int rt_hf_bbox(struct rt_db_internal *ip, point_t *min_pt, point_t *max_pt, const struct bn_tol *tol)
Definition: hf.c:196
char * bu_vls_addr(const struct bu_vls *vp)
Definition: vls.c:111
void bu_vls_struct_print(struct bu_vls *vls, const struct bu_structparse *sdp, const char *base)
#define BU_STRUCTPARSE_FUNC_NULL
Definition: parse.h:153
void bu_semaphore_release(unsigned int i)
Definition: semaphore.c:218
ustring alpha
vect_t r_dir
Direction of ray (UNIT Length)
Definition: raytrace.h:219
Definition: hf.c:94
int rt_hf_prep(struct soltab *stp, struct rt_db_internal *ip, struct rt_i *rtip)
Definition: hf.c:283
int bu_vls_gets(struct bu_vls *vp, FILE *fp)
Definition: vls.c:621
struct bn_tol rti_tol
Math tolerances for this model.
Definition: raytrace.h:1765
#define RT_CK_DBI(_p)
Definition: raytrace.h:829
#define RT_GET_SEG(p, res)
Definition: raytrace.h:379
double perp
nearly 0
Definition: tol.h:75
vect_t hf_VO
Definition: hf.c:96
#define ZERO(val)
Definition: units.c:38
void * idb_ptr
Definition: raytrace.h:195
point_t r_pt
Point at which ray starts.
Definition: raytrace.h:218
point_t st_min
min X, Y, Z of bounding RPP
Definition: raytrace.h:437
const struct rt_functab OBJ[]
Definition: table.c:159
#define BU_SEM_SYSCALL
Definition: parallel.h:178
void bn_vec_ortho(vect_t out, const vect_t in)
int bu_struct_parse(const struct bu_vls *in_vls, const struct bu_structparse *desc, const char *base, void *data)
Definition: parse.c:878
int rt_hf_import4(struct rt_db_internal *ip, const struct bu_external *ep, const fastf_t *mat, const struct db_i *dbip)
Definition: hf.c:1918
int a_y
Screen Y of ray, if applicable.
Definition: raytrace.h:1597
#define HF_GET(_p, _x, _y)
void bu_vls_printf(struct bu_vls *vls, const char *fmt,...) _BU_ATTR_PRINTF23
Definition: vls.c:694
vect_t hf_Y
Definition: hf.c:99
void * st_specific
-> ID-specific (private) struct
Definition: raytrace.h:435
fastf_t hf_Ylen
Definition: hf.c:100
fastf_t uv_du
delta in u
Definition: raytrace.h:343
void bu_close_mapped_file(struct bu_mapped_file *mp)
Definition: mappedfile.c:339
#define MAXHITS
Definition: hf.c:676
fastf_t crv_c1
curvature in principle dir
Definition: raytrace.h:308
Definition: color.c:51
void bu_vls_strcpy(struct bu_vls *vp, const char *s)
Definition: vls.c:310
int hit_surfno
solid-specific surface indicator
Definition: raytrace.h:255
#define RT_HF_INTERNAL_MAGIC
Definition: magic.h:97
vect_t hf_V
Definition: hf.c:95
#define RT_SEM_MODEL
Definition: raytrace.h:1733
void bu_free(void *ptr, const char *str)
Definition: malloc.c:328
vect_t hit_normal
DEPRECATED: Surface Normal at hit_point, use RT_HIT_NORMAL.
Definition: raytrace.h:252
#define BU_CK_LIST_HEAD(_p)
Definition: list.h:142
#define BU_CK_EXTERNAL(_p)
Definition: parse.h:224
#define BU_VLS_INIT_ZERO
Definition: vls.h:84
void rt_hf_ifree(struct rt_db_internal *ip)
Definition: hf.c:2189
size_t ext_nbytes
Definition: parse.h:210
void rt_hf_uv(struct application *ap, struct soltab *stp, struct hit *hitp, struct uvcoord *uvp)
Definition: hf.c:1602
fastf_t hit_dist
dist from r_pt to hit_point
Definition: raytrace.h:250
HIDDEN void verbose(struct human_data_t *dude)
Definition: human.c:2008
Definition: vls.h:56
void bu_bomb(const char *str) _BU_ATTR_NORETURN
Definition: bomb.c:91
HIDDEN const point_t delta
Definition: sh_prj.c:618
double fastf_t
Definition: defines.h:300
#define VPRINT(a, b)
Definition: raytrace.h:1881
int bu_cv_optimize(int cookie)
void rt_db_free_internal(struct rt_db_internal *ip)
Definition: dir.c:216
int rt_hf_tess(struct nmgregion **r, struct model *m, struct rt_db_internal *ip, const struct rt_tess_tol *ttol, const struct bn_tol *tol)
Definition: hf.c:1898
Definition: color.c:50
int hf_shorts
Definition: hf.c:107
fastf_t uv_v
Range 0..1.
Definition: raytrace.h:342
point_t st_center
Centroid of solid.
Definition: raytrace.h:432
#define RT_HIT_MAGIC
Definition: magic.h:161