BRL-CAD
dsp.c
Go to the documentation of this file.
1 /* D S P . C
2  * BRL-CAD
3  *
4  * Copyright (c) 1999-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/dsp/dsp.c
23  *
24  * Intersect a ray with a displacement map.
25  *
26  * The bounding box planes (in dsp coordinates) are numbered 0 .. 5
27  *
28  * For purposes of the "struct hit" surface number, the
29  * "non-elevation" surfaces are numbered 0 .. 7 where:
30  *
31  * Plane # Name plane dist
32  * --------------------------------------------------
33  * 0 XMIN (dist = 0)
34  * 1 XMAX (dist = xsiz)
35  * 2 YMIN (dist = 0)
36  * 3 YMAX (dist = ysiz)
37  * 4 ZMIN (dist = 0)
38  * 5 ZMAX (dsp_max)
39  *
40  * 6 ZMID (dsp_min)
41  * 7 ZTOP (computed)
42  *
43  * if the "struct hit" surfno surface is ZMAX, then hit_vpriv[X, Y]
44  * holds the cell that was hit. hit_vpriv[Z] is 0 if this was an
45  * in-hit. 1 if an out-hit.
46  *
47  */
48 
49 #include "common.h"
50 
51 #include <stddef.h>
52 #include <stdio.h>
53 #include <string.h>
54 #include <math.h>
55 #include <setjmp.h>
56 #include "bnetwork.h"
57 
58 #include "bu/cv.h"
59 #include "bu/parallel.h"
60 #include "vmath.h"
61 #include "raytrace.h"
62 #include "rtgeom.h"
63 #include "db.h"
64 #include "plot3.h"
65 
66 /* private header */
67 #include "./dsp.h"
68 
69 
70 #define FULL_DSP_DEBUGGING 1
71 
72 #define ORDERED_ISECT 1
73 
74 #define DIM_BB_CHILDREN 4
75 #define NUM_BB_CHILDREN (DIM_BB_CHILDREN*DIM_BB_CHILDREN)
76 
77 #define IMPORT_FAIL(_s) \
78  if (dsp_ip) { \
79  bu_log("rt_dsp_import4(%d) '%s' %s\n", __LINE__, bu_vls_addr(&dsp_ip->dsp_name), _s); \
80  bu_free((char *)dsp_ip, "rt_dsp_import4: dsp_ip"); \
81  } \
82  ip->idb_type = ID_NULL; \
83  ip->idb_ptr = (void *)NULL; \
84  return -2
85 
86 
87 struct dsp_rpp {
88  unsigned short dsp_min[3];
89  unsigned short dsp_max[3];
90 };
91 
92 
93 /**
94  * This structure contains a bounding box for a portion of the DSP
95  * along with information about sub-bounding boxes, and what layer
96  * (resolution) of the DSP this box bounds
97  */
98 struct dsp_bb {
99  uint32_t magic;
100  struct dsp_rpp dspb_rpp; /* our bounding box */
101  /*
102  * the next two elements indicate the number and locations of
103  * sub-bounding rpps.
104  *
105  * dsp_b_ch_dim is typically DIM_BB_CHILDREN, DIM_BB_CHILDREN
106  * except for "border" areas of the array
107  */
108  unsigned short dspb_subcell_size;/* XXX This is not yet computed */
109  unsigned short dspb_ch_dim[2]; /* dimensions of children[] */
111 };
112 
113 
114 #define MAGIC_dsp_bb 234
115 #define DSP_BB_CK(_p) BU_CKMAG(_p, MAGIC_dsp_bb, "dsp_bb")
116 /*
117  * This structure provides a handle to all of the bounding boxes for
118  * the DSP at a particular resolution.
119  */
120 #define LAYER(l, x, y) l->p[l->dim[1]*y+x]
121 struct dsp_bb_layer {
122  unsigned int dim[2]; /* the dimensions of the array at element p */
123  struct dsp_bb *p; /* array of dsp_bb's for this level */
124 };
125 
126 # define XCNT(_p) (((struct rt_dsp_internal *)_p)->dsp_xcnt)
127 # define YCNT(_p) (((struct rt_dsp_internal *)_p)->dsp_ycnt)
128 # define XSIZ(_p) (_p->dsp_i.dsp_xcnt - 1)
129 # define YSIZ(_p) (_p->dsp_i.dsp_ycnt - 1)
130 
131 
132 /* FIXME: rename? */
133 extern int rt_retrieve_binunif(struct rt_db_internal *intern,
134  const struct db_i *dbip,
135  const char *name);
136 
137 
138 #define dlog if (RT_G_DEBUG & DEBUG_HF) bu_log
139 
140 
141 #define BBOX_PLANES 7 /* 2 tops & 5 other sides */
142 #define XMIN 0
143 #define XMAX 1
144 #define YMIN 2
145 #define YMAX 3
146 #define ZMIN 4
147 #define ZMAX 5
148 #define ZMID 6
149 #define ZTOP 7
150 
151 
152 /**
153  * per-solid ray tracing form of solid, including precomputed terms
154  *
155  * The dsp_i element MUST BE FIRST so that we can cast a pointer to a
156  * dsp_specific to a rt_dsp_internal.
157  */
158 struct dsp_specific {
159  struct rt_dsp_internal dsp_i; /* MUST BE FIRST */
161  int xsiz;
162  int ysiz;
163  int layers;
165  struct dsp_bb *bb_array;
166 };
167 
168 
169 struct bbox_isect {
172  int in_surf;
173  int out_surf;
174 };
175 
176 
177 /**
178  * per-ray ray tracing information
179  */
180 struct isect_stuff {
181  struct dsp_specific *dsp;
182  struct bu_list seglist; /**< list of segments */
183  struct xray r; /**< solid space ray */
184  vect_t inv_dir; /**< inverses of ray direction */
185  struct application *ap;
186  struct soltab *stp;
187  struct bn_tol *tol;
188 
189  struct bbox_isect bbox;
191 
192  int num_segs;
193  int dmin, dmax; /* for dsp_in_rpp, {X, Y, Z}MIN/MAX */
194 };
195 
196 
197 /**
198  * plane equations (minus offset distances) for bounding RPP
199  */
200 static const vect_t dsp_pl[BBOX_PLANES] = {
201  {-1.0, 0.0, 0.0},
202  { 1.0, 0.0, 0.0},
203 
204  {0.0, -1.0, 0.0},
205  {0.0, 1.0, 0.0},
206 
207  {0.0, 0.0, -1.0},
208  {0.0, 0.0, 1.0},
209  {0.0, 0.0, 1.0},
210 };
211 
212 
213 /**
214  * This function computes the DSP mtos matrix from the stom matrix
215  * whenever the stom matrix is parsed using a bu_structparse.
216  */
217 HIDDEN void
219  const struct bu_structparse *sp,
220  const char *sp_name,
221  void *base,
222  const char *UNUSED(p),
223  void *UNUSED(data))
224 {
225  struct rt_dsp_internal *dsp_ip = (struct rt_dsp_internal *)base;
226 
227  if (!sp) return;
228  if (!sp_name) return;
229  if (!base) return;
230 
231  bn_mat_inv(dsp_ip->dsp_mtos, dsp_ip->dsp_stom);
232 }
233 
234 
235 HIDDEN void
237  const struct bu_structparse *sp,
238  const char *sp_name,
239  void *base,
240  const char *UNUSED(p),
241  void *UNUSED(data))
242 {
243  struct rt_dsp_internal *dsp_ip = (struct rt_dsp_internal *)base;
244 
245  if (!sp) return;
246  if (!sp_name) return;
247  if (!base) return;
248 
249  dsp_ip->dsp_datasrc = RT_DSP_SRC_V4_FILE;
250  dsp_ip->dsp_bip = (struct rt_db_internal *)NULL;
251 }
252 
253 
254 #define DSP_O(m) bu_offsetof(struct rt_dsp_internal, m)
255 
256 /** only used when editing a v4 database */
257 const struct bu_structparse rt_dsp_parse[] = {
258  {"%V", 1, "file", DSP_O(dsp_name), hook_file, NULL, NULL },
259  {"%i", 1, "sm", DSP_O(dsp_smooth), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL },
260  {"%d", 1, "w", DSP_O(dsp_xcnt), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL },
261  {"%d", 1, "n", DSP_O(dsp_ycnt), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL },
262  {"%f", 16, "stom", DSP_O(dsp_stom), hook_mtos_from_stom, NULL, NULL },
263  {"", 0, (char *)0, 0, BU_STRUCTPARSE_FUNC_NULL, NULL, NULL }
264 };
265 
266 
267 /** only used when editing a v4 database */
268 const struct bu_structparse rt_dsp_ptab[] = {
269  {"%V", 1, "file", DSP_O(dsp_name), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL },
270  {"%i", 1, "sm", DSP_O(dsp_smooth), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL },
271  {"%d", 1, "w", DSP_O(dsp_xcnt), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL },
272  {"%d", 1, "n", DSP_O(dsp_ycnt), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL },
273  {"%f", 16, "stom", DSP_O(dsp_stom), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL },
274  {"", 0, (char *)0, 0, BU_STRUCTPARSE_FUNC_NULL, NULL, NULL }
275 };
276 
277 
278 /*
279  * This one is used by the rt_dsp_tclget()
280  */
281 
282 
283 static int plot_file_num = 0;
284 
285 
286 /**
287  * Plot an RPP to a file in the given color
288  */
289 HIDDEN void
290 plot_rpp(FILE *fp, struct bound_rpp *rpp, int r, int g, int b)
291 {
292  pl_color(fp, r, g, b);
293 
294  pd_3move(fp, rpp->min[X], rpp->min[Y], rpp->min[Z]);
295  pd_3cont(fp, rpp->max[X], rpp->min[Y], rpp->min[Z]);
296  pd_3cont(fp, rpp->max[X], rpp->max[Y], rpp->min[Z]);
297  pd_3cont(fp, rpp->min[X], rpp->max[Y], rpp->min[Z]);
298  pd_3cont(fp, rpp->min[X], rpp->min[Y], rpp->min[Z]);
299 
300  pd_3cont(fp, rpp->min[X], rpp->min[Y], rpp->max[Z]);
301  pd_3cont(fp, rpp->max[X], rpp->min[Y], rpp->max[Z]);
302  pd_3cont(fp, rpp->max[X], rpp->max[Y], rpp->max[Z]);
303  pd_3cont(fp, rpp->min[X], rpp->max[Y], rpp->max[Z]);
304  pd_3cont(fp, rpp->min[X], rpp->min[Y], rpp->max[Z]);
305 }
306 
307 
308 /**
309  * Plot a dsp_bb structure
310  */
311 HIDDEN void
312 plot_dsp_bb(FILE *fp, struct dsp_bb *dsp_bb,
313  struct dsp_specific *dsp,
314  int r, int g, int b, int blather)
315 {
316  fastf_t *stom = &dsp->dsp_i.dsp_stom[0];
317  struct bound_rpp rpp;
318  point_t pt;
319 
320  DSP_BB_CK(dsp_bb);
321 
322  VMOVE(pt, dsp_bb->dspb_rpp.dsp_min); /* int->float conversion */
323  MAT4X3PNT(rpp.min, stom, pt);
324 
325  VMOVE(pt, dsp_bb->dspb_rpp.dsp_max); /* int->float conversion */
326  MAT4X3PNT(rpp.max, stom, pt);
327 
328  if (blather)
329  bu_log(" (%g %g %g) (%g %g %g)\n",
330  V3ARGS(rpp.min), V3ARGS(rpp.max));
331 
332  plot_rpp(fp, &rpp, r, g, b);
333 }
334 
335 
336 /*
337  * drawing support for isect_ray_dsp_bb()
338  */
339 HIDDEN FILE *
340 draw_dsp_bb(int *plotnum,
341  struct dsp_bb *dsp_bb,
342  struct dsp_specific *dsp,
343  int r, int g, int b)
344 {
345  char buf[64];
346  FILE *fp;
347  struct dsp_bb bb;
348 
349  sprintf(buf, "dsp_bb%03d.plot3", (*plotnum)++);
350  if ((fp=fopen(buf, "wb")) == (FILE *)NULL) {
351  perror(buf);
352  bu_bomb("");
353  }
354 
355  bu_log("plotting %s", buf);
356  bb = *dsp_bb; /* struct copy */
357  bb.dspb_rpp.dsp_min[Z] = 0.0;
358  plot_dsp_bb(fp, &bb, dsp, r, g, b, 1);
359 
360  return fp;
361 }
362 
363 
364 #define PLOT_LAYERS
365 #ifdef PLOT_LAYERS
366 
367 
368 /**
369  * Plot the bounding box layers for a dsp
370  */
371 HIDDEN void
372 plot_layers(struct dsp_specific *dsp_sp)
373 {
374  FILE *fp;
375  int l, n;
376  unsigned int x, y;
377  char buf[32];
378  static int colors[7][3] = {
379  {255, 0, 0},
380  {0, 255, 0},
381  {0, 0, 255},
382  {255, 255, 0},
383  {255, 0, 255},
384  {0, 255, 255},
385  {255, 255, 255}
386  };
387  int r, g, b, c;
388  struct dsp_bb *d_bb;
389 
390  for (l = 0; l < dsp_sp->layers; l++) {
392  sprintf(buf, "Dsp_layer%d.plot3", l);
393  fp=fopen(buf, "wb");
395  if (fp == (FILE *)NULL) {
396  bu_log("%s:%d error opening %s\n", __FILE__, __LINE__,
397  buf);
398  return;
399  } else
400  bu_log("plotting \"%s\" dim:%d, %d\n", buf,
401  dsp_sp->layer[l].dim[X],
402  dsp_sp->layer[l].dim[Y]);
403  c = l % 6;
404  r = colors[c][0];
405  g = colors[c][1];
406  b = colors[c][2];
407 
408  for (y = 0; y < dsp_sp->layer[l].dim[Y]; y+= 2) {
409  for (x = 0; x < dsp_sp->layer[l].dim[X]; x+= 2) {
410  n = y * dsp_sp->layer[l].dim[X] + x;
411  d_bb = &dsp_sp->layer[l].p[n];
412  plot_dsp_bb(fp, d_bb, dsp_sp, r, g, b, 0);
413 
414  }
415  }
416  fclose(fp);
417  }
418 }
419 
420 
421 #endif
422 
423 
424 /**
425  * Plot the results of intersecting a ray with the top of a cell
426  */
427 HIDDEN void
429  struct dsp_bb *dsp_bb,
430  point_t A,
431  point_t B,
432  point_t C,
433  point_t D,
434  struct hit hitlist[],
435  int hitflags,
436  int style) /* plot diagonal */
437 {
438  fastf_t *stom = &isect->dsp->dsp_i.dsp_stom[0];
439  char buf[64];
440  static int plotcnt = 0;
441  static int cnt = 0;
442  FILE *fp;
443  point_t p1, p2, p3, p4;
444  int i;
445  int in_seg;
446  static unsigned char colors[4][3] = {
447  {255, 255, 128},
448  {255, 128, 255},
449  {128, 255, 128},
450  {128, 255, 255},
451  };
452 
453  DSP_BB_CK(dsp_bb);
454 
456  if (style)
457  sprintf(buf, "dsp_cell_isect%04d.plot3", cnt++);
458  else
459  sprintf(buf, "dsp_cell_top%04d.plot3", plotcnt++);
460 
461  fp=fopen(buf, "wb");
462 
464 
465  if (fp == (FILE *)NULL) {
466  bu_log("error opening \"%s\"\n", buf);
467  return;
468  } else {
469  bu_log("plotting %s flags 0x%x\n\t", buf, hitflags);
470  }
471 
472  plot_dsp_bb(fp, dsp_bb, isect->dsp, 128, 128, 128, 1);
473 
474  /* plot the triangulation */
475  pl_color(fp, 255, 255, 255);
476  MAT4X3PNT(p1, stom, A);
477  MAT4X3PNT(p2, stom, B);
478  MAT4X3PNT(p3, stom, C);
479  MAT4X3PNT(p4, stom, D);
480 
481  pdv_3move(fp, p1);
482  if (style) {
483  pdv_3cont(fp, p2);
484  pdv_3cont(fp, p4);
485  pdv_3cont(fp, p1);
486  pdv_3cont(fp, p3);
487  pdv_3cont(fp, p4);
488  } else {
489  pdv_3cont(fp, p2);
490  pdv_3cont(fp, p4);
491  pdv_3cont(fp, p3);
492  pdv_3cont(fp, p1);
493  }
494 
495  /* plot the hit points */
496 
497  for (in_seg = 0, i = 0; i < 4; i++) {
498  if (hitflags & (1<<i)) {
499  if (in_seg) {
500  in_seg = 0;
501  MAT4X3PNT(p1, stom, hitlist[i].hit_point);
502  pdv_3cont(fp, p1);
503  } else {
504  in_seg = 1;
505  pl_color(fp, colors[i][0], colors[i][1], colors[i][2]);
506  MAT4X3PNT(p1, stom, hitlist[i].hit_point);
507  pdv_3move(fp, p1);
508  }
509  }
510  }
511  fclose(fp);
512 }
513 
514 
515 HIDDEN void
516 dsp_print_v4(struct bu_vls *vls, const struct rt_dsp_internal *dsp_ip)
517 {
518  point_t pt, v;
519  RT_DSP_CK_MAGIC(dsp_ip);
520  BU_CK_VLS(&dsp_ip->dsp_name);
521  BU_CK_VLS(vls);
522 
523  bu_vls_printf(vls, "Displacement Map\n file='%s' w=%u n=%u sm=%d",
524  bu_vls_addr(&dsp_ip->dsp_name),
525  dsp_ip->dsp_xcnt,
526  dsp_ip->dsp_ycnt,
527  dsp_ip->dsp_smooth);
528 
529  VSETALL(pt, 0.0);
530 
531  MAT4X3PNT(v, dsp_ip->dsp_stom, pt);
532 
533  bu_vls_printf(vls, " (origin at %g %g %g)mm\n", V3INTCLAMPARGS(v));
534 
535  bu_vls_printf(vls, " stom=\n");
536  bu_vls_printf(vls, " %8.3f %8.3f %8.3f %8.3f\n",
537  V4INTCLAMPARGS(dsp_ip->dsp_stom));
538 
539  bu_vls_printf(vls, " %8.3f %8.3f %8.3f %8.3f\n",
540  V4INTCLAMPARGS(&dsp_ip->dsp_stom[4]));
541 
542  bu_vls_printf(vls, " %8.3f %8.3f %8.3f %8.3f\n",
543  V4INTCLAMPARGS(&dsp_ip->dsp_stom[8]));
544 
545  bu_vls_printf(vls, " %8.3f %8.3f %8.3f %8.3f\n",
546  V4INTCLAMPARGS(&dsp_ip->dsp_stom[12]));
547 }
548 
549 
550 HIDDEN void
551 dsp_print_v5(struct bu_vls *vls,
552  const struct rt_dsp_internal *dsp_ip)
553 {
554  point_t pt, v;
555 
556  BU_CK_VLS(vls);
557 
558  RT_DSP_CK_MAGIC(dsp_ip);
559  BU_CK_VLS(&dsp_ip->dsp_name);
560 
561 
562  bu_vls_printf(vls, "Displacement Map\n");
563 
564  switch (dsp_ip->dsp_datasrc) {
565  case RT_DSP_SRC_V4_FILE:
566  bu_vls_printf(vls, " Error Error Error");
567  break;
568  case RT_DSP_SRC_FILE:
569  bu_vls_printf(vls, " file");
570  break;
571  case RT_DSP_SRC_OBJ:
572  bu_vls_printf(vls, " obj");
573  break;
574  default:
575  bu_vls_printf(vls, "unk src type'%c'", dsp_ip->dsp_datasrc);
576  break;
577  }
578 
579  bu_vls_printf(vls, "='%s'\n w=%u n=%u sm=%d ",
580  bu_vls_addr(&dsp_ip->dsp_name),
581  dsp_ip->dsp_xcnt,
582  dsp_ip->dsp_ycnt,
583  dsp_ip->dsp_smooth);
584 
585  switch (dsp_ip->dsp_cuttype) {
586  case DSP_CUT_DIR_ADAPT:
587  bu_vls_printf(vls, "cut=ad"); break;
588  case DSP_CUT_DIR_llUR:
589  bu_vls_printf(vls, "cut=lR"); break;
590  case DSP_CUT_DIR_ULlr:
591  bu_vls_printf(vls, "cut=Lr"); break;
592  default:
593  bu_vls_printf(vls, "cut bogus('%c'/%d)",
594  dsp_ip->dsp_cuttype,
595  dsp_ip->dsp_cuttype); break;
596  }
597 
598 
599  VSETALL(pt, 0.0);
600 
601  MAT4X3PNT(v, dsp_ip->dsp_stom, pt);
602 
603  bu_vls_printf(vls, " (origin at %g %g %g)mm\n", V3INTCLAMPARGS(v));
604 
605  bu_vls_printf(vls, " stom=\n");
606  bu_vls_printf(vls, " %8.3f %8.3f %8.3f %8.3f\n", V4INTCLAMPARGS(dsp_ip->dsp_stom));
607  bu_vls_printf(vls, " %8.3f %8.3f %8.3f %8.3f\n", V4INTCLAMPARGS(&dsp_ip->dsp_stom[4]));
608  bu_vls_printf(vls, " %8.3f %8.3f %8.3f %8.3f\n", V4INTCLAMPARGS(&dsp_ip->dsp_stom[8]));
609  bu_vls_printf(vls, " %8.3f %8.3f %8.3f %8.3f\n", V4INTCLAMPARGS(&dsp_ip->dsp_stom[12]));
610 }
611 
612 
613 void
614 rt_dsp_print(register const struct soltab *stp)
615 {
616  register const struct dsp_specific *dsp =
617  (struct dsp_specific *)stp->st_specific;
618  struct bu_vls vls = BU_VLS_INIT_ZERO;
619 
620 
621  RT_DSP_CK_MAGIC(dsp);
622  BU_CK_VLS(&dsp->dsp_i.dsp_name);
623 
624  bu_vls_printf(&vls, "\n---------db version: %d----------\n",
625  db_version(stp->st_rtip->rti_dbip));
626 
627  switch (db_version(stp->st_rtip->rti_dbip)) {
628  case 4:
629  BU_CK_VLS(&vls);
630  dsp_print_v4(&vls, &(dsp->dsp_i));
631  break;
632  case 5:
633  BU_CK_VLS(&vls);
634  dsp_print_v5(&vls, &(dsp->dsp_i));
635  break;
636  }
637 
638  bu_log("%s", bu_vls_addr(&vls));
639 
640  if (BU_VLS_IS_INITIALIZED(&vls)) bu_vls_free(&vls);
641 
642 }
643 
644 
645 /**
646  * compute bounding boxes for each cell, then compute bounding boxes
647  * for collections of bounding boxes
648  */
649 HIDDEN void
650 dsp_layers(struct dsp_specific *dsp, unsigned short *d_min, unsigned short *d_max)
651 {
652  int idx, curr_layer, xs, ys, xv, yv, tot;
653  unsigned int x, y, i, j, k;
654  unsigned short dsp_min, dsp_max, cell_min, cell_max;
655  unsigned short elev;
656  struct dsp_bb *dsp_bb;
657  struct dsp_rpp *t;
658  struct dsp_bb_layer *curr, *prev;
659  unsigned short subcell_size;
660 
661  /* First we compute the total number of struct dsp_bb's we will need */
662  xs = dsp->xsiz;
663  ys = dsp->ysiz;
664  tot = xs * ys;
665  /* bu_log("layer %d %dx%d\n", 0, xs, ys); */
666  dsp->layers = 1;
667  while (xs > 1 || ys > 1) {
668  xv = xs / DIM_BB_CHILDREN;
669  yv = ys / DIM_BB_CHILDREN;
670  if (xs % DIM_BB_CHILDREN) xv++;
671  if (ys % DIM_BB_CHILDREN) yv++;
672 
673 #ifdef FULL_DSP_DEBUGGING
674  if (RT_G_DEBUG & DEBUG_HF)
675  bu_log("layer %d %dx%d\n", dsp->layers, xv, yv);
676 #endif
677  tot += xv * yv;
678 
679  if (xv > 0) xs = xv;
680  else xs = 1;
681 
682  if (yv > 0) ys = yv;
683  else ys = 1;
684  dsp->layers++;
685  }
686 
687 
688 #ifdef FULL_DSP_DEBUGGING
689  if (RT_G_DEBUG & DEBUG_HF)
690  bu_log("%d layers total\n", dsp->layers);
691 #endif
692 
693  /* allocate the struct dsp_bb's we will need */
694  dsp->layer = (struct dsp_bb_layer *)bu_malloc(dsp->layers * sizeof(struct dsp_bb_layer),
695  "dsp_bb_layers array");
696  dsp->bb_array = (struct dsp_bb *)bu_malloc(tot * sizeof(struct dsp_bb), "dsp_bb array");
697 
698  /* now we fill in the "lowest" layer of struct dsp_bb's from the
699  * raw data
700  */
701  dsp->layer[0].dim[X] = dsp->xsiz;
702  dsp->layer[0].dim[Y] = dsp->ysiz;
703  dsp->layer[0].p = dsp->bb_array;
704 
705  xs = dsp->xsiz;
706  ys = dsp->ysiz;
707 
708  dsp_min = 0xffff;
709  dsp_max = 0;
710 
711  for (y = 0; y < YSIZ(dsp); y++) {
712 
713  cell_min = 0xffff;
714  cell_max = 0;
715 
716  for (x = 0; x < XSIZ(dsp); x++) {
717 
718  elev = DSP(&dsp->dsp_i, x, y);
719  cell_min = cell_max = elev;
720 
721  elev = DSP(&dsp->dsp_i, x+1, y);
722  V_MIN(cell_min, elev);
723  V_MAX(cell_max, elev);
724 
725  elev = DSP(&dsp->dsp_i, x, y+1);
726  V_MIN(cell_min, elev);
727  V_MAX(cell_max, elev);
728 
729  elev = DSP(&dsp->dsp_i, x+1, y+1);
730  V_MIN(cell_min, elev);
731  V_MAX(cell_max, elev);
732 
733  /* factor the cell min/max into the overall min/max */
734  V_MIN(dsp_min, cell_min);
735  V_MAX(dsp_max, cell_max);
736 
737  /* fill in the dsp_rpp cell min/max */
738  i = y*XSIZ(dsp) + x;
739  dsp_bb = &dsp->layer[0].p[i];
740  VSET(dsp_bb->dspb_rpp.dsp_min, x, y, cell_min);
741  VSET(dsp_bb->dspb_rpp.dsp_max, x+1, y+1, cell_max);
742 
743  dsp_bb->dspb_subcell_size = 0;
744 
745  /* There are no "children" of a layer 0 element */
746  dsp_bb->dspb_ch_dim[X] = 0;
747  dsp_bb->dspb_ch_dim[Y] = 0;
748  for (k = 0; k < NUM_BB_CHILDREN; k++) {
749  dsp_bb->dspb_children[k] =
750  (struct dsp_bb *)NULL;
751  }
752  dsp_bb->magic = MAGIC_dsp_bb;
753 
754  /* XXX should we compute the triangle orientation and
755  * save it here too?
756  */
757  }
758  }
759 
760  *d_min = dsp_min;
761  *d_max = dsp_max;
762 
763 
764  if (RT_G_DEBUG & DEBUG_HF)
765  bu_log("layer 0 filled\n");
766 
767  subcell_size = 1;
768 
769  /* now we compute successive layers from the initial layer */
770  for (curr_layer = 1; curr_layer < dsp->layers; curr_layer++) {
771  /* compute the number of cells in each direction for this layer */
772 
773  xs = dsp->layer[curr_layer-1].dim[X];
774  if (xs % DIM_BB_CHILDREN)
775  dsp->layer[curr_layer].dim[X] =
776  xs / DIM_BB_CHILDREN + 1;
777  else
778  dsp->layer[curr_layer].dim[X] =
779  xs / DIM_BB_CHILDREN;
780 
781  ys = dsp->layer[curr_layer-1].dim[Y];
782  if (ys % DIM_BB_CHILDREN)
783  dsp->layer[curr_layer].dim[Y] =
784  ys / DIM_BB_CHILDREN + 1;
785  else
786  dsp->layer[curr_layer].dim[Y] =
787  ys / DIM_BB_CHILDREN;
788 
789  /* set the start of the array for this layer */
790  dsp->layer[curr_layer].p =
791  &dsp->layer[curr_layer-1].p[dsp->layer[curr_layer-1].dim[X] * dsp->layer[curr_layer-1].dim[Y] ];
792 
793  curr = &dsp->layer[curr_layer];
794  prev = &dsp->layer[curr_layer-1];
795 
796  if (RT_G_DEBUG & DEBUG_HF)
797  bu_log("layer %d subcell size %d\n", curr_layer, subcell_size);
798 
799  /* walk the grid and fill in the values for this layer */
800  for (y = 0; y < curr->dim[Y]; y++) {
801  for (x = 0; x < curr->dim[X]; x++) {
802  int n, xp, yp;
803  /* x, y are in the coordinates in the current
804  * layer. xp, yp are the coordinates of the
805  * same area in the previous (lower) layer.
806  */
807  xp = x * DIM_BB_CHILDREN;
808  yp = y * DIM_BB_CHILDREN;
809 
810  /* initialize the current dsp_bb cell */
811  dsp_bb = &curr->p[y*curr->dim[X]+x];
812  dsp_bb->magic = MAGIC_dsp_bb;
813  n = lrint(pow((double)DIM_BB_CHILDREN, (double)curr_layer));
814  VSET(dsp_bb->dspb_rpp.dsp_min,
815  x * n, y * n, 0x0ffff);
816  VSET(dsp_bb->dspb_rpp.dsp_max,
817  x * n, y * n, 0);
818 
819  /* record the dimensions of our children */
820  dsp_bb->dspb_subcell_size = subcell_size;
821 
822 
823  tot = 0;
824  i = 0;
825  for (j = 0; j < DIM_BB_CHILDREN && (yp+j)<prev->dim[Y]; j++) {
826  for (i = 0; i < DIM_BB_CHILDREN && (xp+i)<prev->dim[X]; i++) {
827 
828  idx = (yp+j) * prev->dim[X] + xp+i;
829 
830  t = &prev->p[ idx ].dspb_rpp;
831 
832  VMINMAX(dsp_bb->dspb_rpp.dsp_min,
833  dsp_bb->dspb_rpp.dsp_max, t->dsp_min);
834  VMINMAX(dsp_bb->dspb_rpp.dsp_min,
835  dsp_bb->dspb_rpp.dsp_max, t->dsp_max);
836 
837  dsp_bb->dspb_children[tot++] = &prev->p[ idx ];
838 
839  }
840  }
841 
842  dsp_bb->dspb_ch_dim[X] = i;
843  dsp_bb->dspb_ch_dim[Y] = j;
844  }
845  }
846  subcell_size *= DIM_BB_CHILDREN;
847  }
848 
849 #ifdef PLOT_LAYERS
850  if (RT_G_DEBUG & DEBUG_HF) {
851  plot_layers(dsp);
852  bu_log("_ x:%u y:%u min %d max %d\n",
853  XCNT(dsp), YCNT(dsp), dsp_min, dsp_max);
854  }
855 #endif
856 }
857 
858 /**
859  * Calculate the bounding box for a dsp.
860  */
861 int
862 rt_dsp_bbox(struct rt_db_internal *ip, point_t *min, point_t *max, const struct bn_tol *UNUSED(tol)) {
863  struct rt_dsp_internal *dsp_ip;
864  struct dsp_specific ds;
865  unsigned short dsp_min, dsp_max;
866  point_t pt, bbpt;
867 
868  RT_CK_DB_INTERNAL(ip);
869  dsp_ip = (struct rt_dsp_internal *)ip->idb_ptr;
870  RT_DSP_CK_MAGIC(dsp_ip);
871  BU_CK_VLS(&dsp_ip->dsp_name);
872 
873  switch (dsp_ip->dsp_datasrc) {
874  case RT_DSP_SRC_V4_FILE:
875  case RT_DSP_SRC_FILE:
876  if (!dsp_ip->dsp_mp) {
877  bu_log("dsp(%s): no data file or data file empty\n", bu_vls_addr(&dsp_ip->dsp_name));
878  return 1; /* BAD */
879  }
880  BU_CK_MAPPED_FILE(dsp_ip->dsp_mp);
881 
882  /* we do this here and now because we will need it for the
883  * dsp_specific structure in a few lines
884  */
886  ++dsp_ip->dsp_mp->uses;
888  break;
889  case RT_DSP_SRC_OBJ:
890  if (!dsp_ip->dsp_bip) {
891  bu_log("dsp(%s): no data\n", bu_vls_addr(&dsp_ip->dsp_name));
892  return 1; /* BAD */
893  }
894  RT_CK_DB_INTERNAL(dsp_ip->dsp_bip);
895  RT_CK_BINUNIF(dsp_ip->dsp_bip->idb_ptr);
896  break;
897  }
898 
899  memset(&ds, 0, sizeof(struct dsp_specific));
900 
901  /* this works ok, because the mapped file keeps track of the
902  * number of uses. However, the binunif interface does not.
903  * We'll have to copy the data for that one.
904  */
905  ds.dsp_i = *dsp_ip; /* struct copy */
906 
907  /* this keeps the binary internal object from being freed */
908  dsp_ip->dsp_bip = (struct rt_db_internal *)NULL;
909 
910 
911  ds.xsiz = dsp_ip->dsp_xcnt-1; /* size is # cells or values-1 */
912  ds.ysiz = dsp_ip->dsp_ycnt-1; /* size is # cells or values-1 */
913 
914 
915  /* compute the multi-resolution bounding boxes */
916  dsp_layers(&ds, &dsp_min, &dsp_max);
917 
918 
919  /* record the distance to each of the bounding planes */
920  ds.dsp_pl_dist[XMIN] = 0.0;
921  ds.dsp_pl_dist[XMAX] = (fastf_t)ds.xsiz;
922  ds.dsp_pl_dist[YMIN] = 0.0;
923  ds.dsp_pl_dist[YMAX] = (fastf_t)ds.ysiz;
924  ds.dsp_pl_dist[ZMIN] = 0.0;
925  ds.dsp_pl_dist[ZMAX] = (fastf_t)dsp_max;
926  ds.dsp_pl_dist[ZMID] = (fastf_t)dsp_min;
927 
928  /* compute enlarged bounding box and sphere */
929  VSETALL((*min), INFINITY);
930  VSETALL((*max), -INFINITY);
931 
932 #define BBOX_PT(_x, _y, _z) \
933  VSET(pt, (fastf_t)_x, (fastf_t)_y, (fastf_t)_z); \
934  MAT4X3PNT(bbpt, dsp_ip->dsp_stom, pt); \
935  VMINMAX((*min), (*max), bbpt)
936 
937  BBOX_PT(-.1, -.1, -.1);
938  BBOX_PT(dsp_ip->dsp_xcnt+.1, -.1, -.1);
939  BBOX_PT(dsp_ip->dsp_xcnt+.1, dsp_ip->dsp_ycnt+.1, -1);
940  BBOX_PT(-.1, dsp_ip->dsp_ycnt+.1, -1);
941  BBOX_PT(-.1, -.1, dsp_max+.1);
942  BBOX_PT(dsp_ip->dsp_xcnt+.1, -.1, dsp_max+.1);
943  BBOX_PT(dsp_ip->dsp_xcnt+.1, dsp_ip->dsp_ycnt+.1, dsp_max+.1);
944  BBOX_PT(-.1, dsp_ip->dsp_ycnt+.1, dsp_max+.1);
945 
946 #undef BBOX_PT
947 
948  switch (ds.dsp_i.dsp_datasrc) {
949  case RT_DSP_SRC_V4_FILE:
950  case RT_DSP_SRC_FILE:
951  if (ds.dsp_i.dsp_mp) {
952  bu_close_mapped_file(ds.dsp_i.dsp_mp);
953  } else if (ds.dsp_i.dsp_buf) {
954  bu_free(ds.dsp_i.dsp_buf, "dsp fake data");
955  }
956  break;
957  case RT_DSP_SRC_OBJ:
958  break;
959  }
960 
961  return 0;
962 }
963 
964 /**
965  * Given a pointer to a GED database record, and a transformation
966  * matrix, determine if this is a valid DSP, and if so, precompute
967  * various terms of the formula.
968  *
969  * Returns -
970  * 0 DSP is OK
971  * !0 Error in description
972  *
973  * Implicit return -
974  * A struct dsp_specific is created, and its address is stored in
975  * stp->st_specific for use by dsp_shot().
976  *
977  * Note: because the stand-along bbox calculation requires much
978  * of the prep logic, the in-prep bbox calculations are left
979  * in to avoid duplication rather than calling rt_dsp_bbox.
980  */
981 int
982 rt_dsp_prep(struct soltab *stp, struct rt_db_internal *ip, struct rt_i *rtip)
983 {
984  struct rt_dsp_internal *dsp_ip;
985  register struct dsp_specific *dsp;
986  unsigned short dsp_min, dsp_max;
987  point_t pt, bbpt;
988  vect_t work;
989  fastf_t f;
990 
991  if (RT_G_DEBUG & DEBUG_HF)
992  bu_log("rt_dsp_prep()\n");
993 
994  if (rtip) RT_CK_RTI(rtip);
995 
996  RT_CK_DB_INTERNAL(ip);
997  dsp_ip = (struct rt_dsp_internal *)ip->idb_ptr;
998  RT_DSP_CK_MAGIC(dsp_ip);
999  BU_CK_VLS(&dsp_ip->dsp_name);
1000 
1001  switch (dsp_ip->dsp_datasrc) {
1002  case RT_DSP_SRC_V4_FILE:
1003  case RT_DSP_SRC_FILE:
1004  if (!dsp_ip->dsp_mp) {
1005  bu_log("dsp(%s): no data file or data file empty\n", bu_vls_addr(&dsp_ip->dsp_name));
1006  return 1; /* BAD */
1007  }
1008  BU_CK_MAPPED_FILE(dsp_ip->dsp_mp);
1009 
1010  /* we do this here and now because we will need it for the
1011  * dsp_specific structure in a few lines
1012  */
1014  ++dsp_ip->dsp_mp->uses;
1016  break;
1017  case RT_DSP_SRC_OBJ:
1018  if (!dsp_ip->dsp_bip) {
1019  bu_log("dsp(%s): no data\n", bu_vls_addr(&dsp_ip->dsp_name));
1020  return 1; /* BAD */
1021  }
1022  RT_CK_DB_INTERNAL(dsp_ip->dsp_bip);
1023  RT_CK_BINUNIF(dsp_ip->dsp_bip->idb_ptr);
1024  break;
1025  }
1026 
1027 
1028  BU_GET(dsp, struct dsp_specific);
1029  stp->st_specific = (void *) dsp;
1030 
1031  /* this works ok, because the mapped file keeps track of the
1032  * number of uses. However, the binunif interface does not.
1033  * We'll have to copy the data for that one.
1034  */
1035  dsp->dsp_i = *dsp_ip; /* struct copy */
1036 
1037  /* this keeps the binary internal object from being freed */
1038  dsp_ip->dsp_bip = (struct rt_db_internal *)NULL;
1039 
1040 
1041  dsp->xsiz = dsp_ip->dsp_xcnt-1; /* size is # cells or values-1 */
1042  dsp->ysiz = dsp_ip->dsp_ycnt-1; /* size is # cells or values-1 */
1043 
1044 
1045  /* compute the multi-resolution bounding boxes */
1046  dsp_layers(dsp, &dsp_min, &dsp_max);
1047 
1048 
1049  /* record the distance to each of the bounding planes */
1050  dsp->dsp_pl_dist[XMIN] = 0.0;
1051  dsp->dsp_pl_dist[XMAX] = (fastf_t)dsp->xsiz;
1052  dsp->dsp_pl_dist[YMIN] = 0.0;
1053  dsp->dsp_pl_dist[YMAX] = (fastf_t)dsp->ysiz;
1054  dsp->dsp_pl_dist[ZMIN] = 0.0;
1055  dsp->dsp_pl_dist[ZMAX] = (fastf_t)dsp_max;
1056  dsp->dsp_pl_dist[ZMID] = (fastf_t)dsp_min;
1057 
1058  /* compute enlarged bounding box and sphere */
1059 
1060 #define BBOX_PT(_x, _y, _z) \
1061  VSET(pt, (fastf_t)_x, (fastf_t)_y, (fastf_t)_z); \
1062  MAT4X3PNT(bbpt, dsp_ip->dsp_stom, pt); \
1063  VMINMAX(stp->st_min, stp->st_max, bbpt)
1064 
1065  BBOX_PT(-.1, -.1, -.1);
1066  BBOX_PT(dsp_ip->dsp_xcnt+.1, -.1, -.1);
1067  BBOX_PT(dsp_ip->dsp_xcnt+.1, dsp_ip->dsp_ycnt+.1, -1);
1068  BBOX_PT(-.1, dsp_ip->dsp_ycnt+.1, -1);
1069  BBOX_PT(-.1, -.1, dsp_max+.1);
1070  BBOX_PT(dsp_ip->dsp_xcnt+.1, -.1, dsp_max+.1);
1071  BBOX_PT(dsp_ip->dsp_xcnt+.1, dsp_ip->dsp_ycnt+.1, dsp_max+.1);
1072  BBOX_PT(-.1, dsp_ip->dsp_ycnt+.1, dsp_max+.1);
1073 
1074 #undef BBOX_PT
1075 
1076  VADD2SCALE(stp->st_center, stp->st_min, stp->st_max, 0.5);
1077  VSUB2SCALE(work, stp->st_max, stp->st_min, 0.5);
1078 
1079  f = work[X];
1080  if (work[Y] > f) f = work[Y];
1081  if (work[Z] > f) f = work[Z];
1082  stp->st_aradius = f;
1083  stp->st_bradius = MAGNITUDE(work);
1084 
1085  if (RT_G_DEBUG & DEBUG_HF) {
1086  bu_log(" model space bbox (%g %g %g) (%g %g %g)\n",
1087  V3ARGS(stp->st_min),
1088  V3ARGS(stp->st_max));
1089  }
1090 
1091  return 0;
1092 }
1093 
1094 
1095 HIDDEN void
1096 plot_seg(struct isect_stuff *isect,
1097  struct hit *in_hit,
1098  struct hit *out_hit,
1099  const point_t bbmin, /* The bounding box of what you are adding ... */
1100  const point_t bbmax, /* ... */
1101  int r, int g, int b)/* ... this is strictly for debug plot purposes */
1102 {
1103  fastf_t *stom = &isect->dsp->dsp_i.dsp_stom[0];
1104  struct bound_rpp rpp;
1105  char fname[32];
1106  FILE *fp;
1107  static int segnum = 0;
1108 
1109  /* plot the bounding box and the seg */
1111  sprintf(fname, "dsp_seg%04d.plot3", segnum++);
1112  fp=fopen(fname, "wb");
1114 
1115  if (fp != (FILE *)NULL) {
1116  bu_log("plotting %s\n", fname);
1117 
1118  MAT4X3PNT(rpp.min, stom, bbmin);
1119  MAT4X3PNT(rpp.max, stom, bbmax);
1120  plot_rpp(fp, &rpp, r/2, g/2, b/2);
1121 
1122  /* re-use the rpp for the points for the segment */
1123  MAT4X3PNT(rpp.min, stom, in_hit->hit_point);
1124  MAT4X3PNT(rpp.max, stom, out_hit->hit_point);
1125 
1126  pl_color(fp, r, g, b);
1127  pdv_3line(fp, rpp.min, rpp.max);
1128 
1130  fclose(fp);
1132  }
1133 }
1134 
1135 
1136 /**
1137  * Add a segment to the list of intersections in DSP space
1138  *
1139  * Return:
1140  * 0 continue to intersect
1141  * 1 All intersections computed, terminate intersection processing
1142  */
1143 HIDDEN int
1144 add_seg(struct isect_stuff *isect,
1145  struct hit *in_hit,
1146  struct hit *out_hit,
1147  const point_t bbmin, /* The bounding box of what you are adding ... */
1148  const point_t bbmax, /* ... */
1149  int r, int g, int b /* ... this is strictly for debug plot purposes */
1150  )
1151 
1152 {
1153  struct seg *seg;
1154  fastf_t tt = isect->tol->dist;
1155  fastf_t delta;
1156 #ifndef ORDERED_ISECT
1157  struct bu_list *spot;
1158 #endif
1159 
1160  tt *= tt;
1161 
1162 #ifdef ORDERED_ISECT
1163  if (BU_LIST_NON_EMPTY(&isect->seglist)) {
1164  /* if the new in-distance equals the old out distance
1165  * we just extend the old segment
1166  */
1167  seg = BU_LIST_LAST(seg, &isect->seglist);
1168  if (fabs(seg->seg_out.hit_dist - in_hit->hit_dist) <= tt) {
1169 
1170  seg->seg_out.hit_dist = out_hit->hit_dist;
1171  seg->seg_out.hit_surfno = out_hit->hit_surfno;
1172  VMOVE(seg->seg_out.hit_normal, out_hit->hit_normal);
1173  if (out_hit->hit_surfno == ZTOP) {
1174  seg->seg_out.hit_vpriv[X] = out_hit->hit_vpriv[X];
1175  seg->seg_out.hit_vpriv[Y] = out_hit->hit_vpriv[Y];
1176  }
1177  seg->seg_out.hit_vpriv[Z] = 1.0; /* flag as out-hit */
1178 
1179  if (RT_G_DEBUG & DEBUG_HF) {
1180  bu_log("extending previous seg to %g\n", out_hit->hit_dist);
1181  plot_seg(isect, in_hit, out_hit, bbmin, bbmax, r, g, b);
1182  }
1183  return 0;
1184  }
1185  }
1186 #else
1187  /* insert the segment in the list by in-hit distance */
1188  dlog("searching for insertion point for seg w in/out dist %g %g\n",
1189  in_hit->hit_dist, out_hit->hit_dist);
1190 
1191  for (BU_LIST_FOR(seg, seg, &isect->seglist)) {
1192  dlog("checking %g->%g seg\n", seg->seg_in.hit_dist,
1193  seg->seg_out.hit_dist);
1194  /* found the spot for this one */
1195  if (fabs(seg->seg_out.hit_dist - in_hit->hit_dist) <= tt) {
1196  seg->seg_out = *out_hit; /* struct copy */
1198  seg->seg_out.hit_vpriv[Z] = 1.0; /* flag as out-hit */
1199 
1200  if (RT_G_DEBUG & DEBUG_HF) {
1201  bu_log("extending previous seg to %g\n",
1202  out_hit->hit_dist);
1203  plot_seg(isect, in_hit, out_hit, bbmin, bbmax, r, g, b);
1204  }
1205  return 0;
1206  }
1207  if (in_hit->hit_dist < seg->seg_in.hit_dist) {
1208  spot = &seg->l;
1209  dlog("insert before this one\n");
1210  goto found_spot;
1211  }
1212  }
1213  spot = &isect->seglist;
1214  seg = BU_LIST_LAST(seg, &isect->seglist);
1215  dlog("insert at end\n");
1216  found_spot:
1217 #endif
1218 
1219 
1220  /* if both points are on the "floor" of the DSP, then we
1221  * don't have a hit segment
1222  */
1223  if (NEAR_ZERO(in_hit->hit_point[Z], isect->tol->dist) &&
1224  NEAR_ZERO(out_hit->hit_point[Z], isect->tol->dist)) {
1225  return 0;
1226  }
1227 
1228  /* throw away any zero length segments,
1229  * mostly to avoid seeing inside-out segments
1230  */
1231  delta = out_hit->hit_dist - in_hit->hit_dist;
1232 
1233  if (ZERO(delta)) {
1234  return 0;
1235  }
1236 
1237  /* if the segment is inside-out, we need to say something about it */
1238  if (delta < 0.0 && !NEAR_ZERO(delta, isect->tol->dist)) {
1239  bu_log(" %s:%dDSP: Adding inside-out seg in:%g out:%g\n",
1240  __FILE__, __LINE__,
1241  in_hit->hit_dist, out_hit->hit_dist);
1242 
1243  VPRINT("\tin_pt", in_hit->hit_point);
1244  VPRINT("\tout_pt", out_hit->hit_point);
1245  }
1246 
1247 
1248  RT_GET_SEG(seg, isect->ap->a_resource);
1249 
1250  seg->seg_in.hit_dist = in_hit->hit_dist;
1251  seg->seg_out.hit_dist = out_hit->hit_dist;
1252 
1253  seg->seg_in.hit_surfno = in_hit->hit_surfno;
1254  seg->seg_out.hit_surfno = out_hit->hit_surfno;
1255 
1256  VMOVE(seg->seg_in.hit_normal, in_hit->hit_normal);
1257  VMOVE(seg->seg_out.hit_normal, out_hit->hit_normal);
1258 
1259  if (in_hit->hit_surfno >= ZMAX) {
1260  seg->seg_in.hit_vpriv[X] = in_hit->hit_vpriv[X];
1261  seg->seg_in.hit_vpriv[Y] = in_hit->hit_vpriv[Y];
1262  }
1263 
1264  if (out_hit->hit_surfno >= ZMAX) {
1265  seg->seg_out.hit_vpriv[X] = out_hit->hit_vpriv[X];
1266  seg->seg_out.hit_vpriv[Y] = out_hit->hit_vpriv[Y];
1267  }
1268 
1269  seg->seg_in.hit_vpriv[Z] = 0.0; /* flag as in-hit */
1270  seg->seg_out.hit_vpriv[Z] = 1.0; /* flag as out-hit */
1271 
1272  seg->seg_stp = isect->stp;
1273 
1274 #ifdef FULL_DSP_DEBUGGING
1275  if (VDOT(seg->seg_in.hit_normal, isect->r.r_dir) > 0) {
1276  bu_log("----------------------------------------------------------\n");
1277  bu_log("pixel %d, %d bogus seg in, inward normal\nN: %g %g %g\nd: %g %g %g\n",
1278  isect->ap->a_x, isect->ap->a_y,
1279  V3ARGS(seg->seg_in.hit_normal), V3ARGS(isect->r.r_dir));
1280  bu_bomb("");
1281  }
1282  if (VDOT(seg->seg_out.hit_normal, isect->r.r_dir) < 0) {
1283  bu_log("----------------------------------------------------------\n");
1284  bu_log("pixel %d, %d bogus seg out, inward normal\nN: %g %g %g\nd: %g %g %g\n",
1285  isect->ap->a_x, isect->ap->a_y,
1286  V3ARGS(seg->seg_out.hit_normal), V3ARGS(isect->r.r_dir));
1287  bu_bomb("");
1288  }
1289 #endif
1290 
1291 #ifdef ORDERED_ISECT
1292  BU_LIST_INSERT(&isect->seglist, &seg->l);
1293 #else
1294  BU_LIST_INSERT(spot, &seg->l);
1295 #endif
1296 
1297 
1298  if (RT_G_DEBUG & DEBUG_HF)
1299  plot_seg(isect, in_hit, out_hit, bbmin, bbmax, r, g, b);
1300 
1301 
1302  if (seg->seg_in.hit_dist > 0.0 || seg->seg_out.hit_dist > 0.0) {
1303  return ++isect->num_segs > isect->ap->a_onehit;
1304  }
1305  return 0;
1306 }
1307 
1308 
1309 /**
1310  * Side Effects:
1311  * dist and P may be set
1312  *
1313  * Return:
1314  * 1 Ray intersects triangle
1315  * 0 Ray misses triangle
1316  * -1 Ray/plane parallel
1317  */
1318 HIDDEN int
1320  point_t A,
1321  point_t B,
1322  point_t C,
1323  struct hit *hitp,
1324  fastf_t alphabbeta[])
1325 {
1326  point_t P; /* plane intercept point */
1327  vect_t AB, AC, AP;
1328  plane_t N; /* Normal for plane of triangle */
1329  fastf_t NdotDir;
1330  fastf_t alpha, beta; /* barycentric distances */
1331  fastf_t hitdist; /* distance to ray/triangle intercept */
1332  fastf_t toldist; /* distance tolerance from isect->tol */
1333 
1334 #ifdef FULL_DSP_DEBUGGING
1335  if (RT_G_DEBUG & DEBUG_HF) {
1336  fastf_t *stom = &isect->dsp->dsp_i.dsp_stom[0];
1337 
1338  MAT4X3PNT(P, stom, A);
1339  bu_log("isect_ray_triangle...\n A %g %g %g (%g %g %g)\n",
1340  V3ARGS(A), V3ARGS(P));
1341 
1342  MAT4X3PNT(P, stom, B);
1343  bu_log(" B %g %g %g (%g %g %g)\n",
1344  V3ARGS(B), V3ARGS(P));
1345 
1346  MAT4X3PNT(P, stom, C);
1347  bu_log(" C %g %g %g (%g %g %g)\n",
1348  V3ARGS(C), V3ARGS(P));
1349 
1350  }
1351 #endif
1352  VSUB2(AB, B, A);
1353  VSUB2(AC, C, A);
1354 
1355  /* Compute the plane equation of the triangle */
1356  VCROSS(N, AB, AC);
1357  VUNITIZE(N);
1358  N[H] = VDOT(N, A);
1359 
1360 
1361  /* intersect ray with plane */
1362  NdotDir = VDOT(N, isect->r.r_dir);
1363  if (BN_VECT_ARE_PERP(NdotDir, isect->tol)) {
1364  /* Ray perpendicular to plane of triangle */
1365  return -1;
1366  }
1367 
1368  /* dist to plane icept */
1369  hitdist = (N[H]-VDOT(N, isect->r.r_pt)) / NdotDir;
1370 
1371  VJOIN1(P, isect->r.r_pt, hitdist, isect->r.r_dir);
1372 
1373 #ifdef FULL_DSP_DEBUGGING
1374  if (RT_G_DEBUG & DEBUG_HF) {
1375  FILE *fp;
1376  char buf[32];
1377  static int plotnum = 0;
1378  fastf_t *stom = &isect->dsp->dsp_i.dsp_stom[0];
1379  point_t p1, p2;
1380 
1382  sprintf(buf, "dsp_tri%03d.plot3", plotnum++);
1383  fp=fopen(buf, "wb");
1385 
1386  if (fp == (FILE *)NULL) {
1387  bu_log("%s:%d error opening \"%s\"\n", __FILE__, __LINE__, buf);
1388  bu_bomb("");
1389  } else
1390  bu_log(" plotting %s\n", buf);
1391 
1392 
1393  pl_color(fp, 255, 255, 255);
1394  MAT4X3PNT(p1, stom, A); pdv_3move(fp, p1);
1395  MAT4X3PNT(p2, stom, B); pdv_3cont(fp, p2);
1396  MAT4X3PNT(p2, stom, C); pdv_3cont(fp, p2);
1397  pdv_3cont(fp, p1);
1398 
1399  /* plot the ray */
1400  pl_color(fp, 255, 255, 0);
1401  MAT4X3PNT(p1, stom, P);
1402  MAT4X3PNT(p2, stom, isect->r.r_pt);
1403  pdv_3line(fp, p1, p2);
1404 
1406  fclose(fp);
1408 
1409  bu_log(" dist:%g plane point: %g %g %g\n", hitdist, V3ARGS(p2));
1410  }
1411 #endif
1412  /* We project the system into the XY plane at this point to
1413  * determine if the ray_plane_isect_pt is within the bounds of the
1414  * triangle
1415  *
1416  * The general idea here is to project the vector AP onto both
1417  * sides of the triangle. The distances along each side can be
1418  * used to determine if we are inside/outside the triangle.
1419  *
1420  * VSUB2(AP, P, A);
1421  * alpha = VDOT(AP, AB);
1422  * beta = VDOT(AP, AC);
1423  *
1424  * C
1425  * |
1426  * |
1427  * |
1428  *---- |--- P
1429  * | | /
1430  * | | / |
1431  * alpha | / |
1432  * | |/ |
1433  *---- A---------B
1434  *
1435  * b
1436  * |-e--|
1437  * t
1438  * a
1439  *
1440  * To save on computation, we do this calculation in 2D.
1441  *
1442  * See: "Graphics Gems", Andrew S. Glassner ed. PP390-393 for details
1443  */
1444 
1445  toldist = isect->tol->dist;
1446 
1447 
1448  VSUB2(AP, P, A);
1449 #ifdef FULL_DSP_DEBUGGING
1450  if (RT_G_DEBUG & DEBUG_HF) {
1451  VPRINT(" AP", AP);
1452  VPRINT(" AB", AB);
1453  VPRINT(" AC", AC);
1454  }
1455 #endif
1456 
1457  /* Ordinarily, in 2D we would say:
1458  *
1459  * beta = AB[X] * AP[X] + AB[Y] * AP[Y];
1460  * alpha = AC[X] * AP[X] + AC[Y] * AP[Y];
1461  *
1462  * However, in this case, we know that AB and AC will always be
1463  * axis-aligned, so we can short-cut. XXX consider: we know that
1464  * AB and AC will be unit length, so only the sign counts.
1465  */
1466 
1467  if (ZERO(AB[X])) {
1468  beta = AB[Y] * AP[Y];
1469  } else {
1470  beta = AB[X] * AP[X];
1471  }
1472 
1473  if (ZERO(AC[X])) {
1474  alpha = AC[Y] * AP[Y];
1475  } else {
1476  alpha = AC[X] * AP[X];
1477  }
1478 
1479  /* return 1 if we hit the triangle */
1480  alphabbeta[0] = alpha;
1481  alphabbeta[1] = beta;
1482 
1483 #ifdef FULL_DSP_DEBUGGING
1484  if (alpha < -toldist) {
1485  dlog("alpha < 0\n");
1486  return 0;
1487  }
1488  if (beta < -toldist) {
1489  dlog("beta < 0\n");
1490  return 0;
1491  }
1492  if ((alpha+beta) > (1.0 + toldist)) {
1493  dlog("alpha+beta > 1\n");
1494  return 0;
1495  }
1496 #else
1497  if (alpha < -toldist || beta < -toldist || (alpha+beta) > (1.0 + toldist))
1498  return 0;
1499 #endif
1500 
1501  hitp->hit_dist = hitdist;
1502  VMOVE(hitp->hit_normal, N);
1503  VMOVE(hitp->hit_point, P);
1504  return 1;
1505 }
1506 
1507 
1508 /**
1509  * For adaptive diagonal selection or for Upper-Left to lower right
1510  * cell cut, we must permute the vertices of the cell before handing
1511  * them to the intersection algorithm. That's what this function
1512  * does.
1513  */
1514 HIDDEN int
1515 permute_cell(point_t A,
1516  point_t B,
1517  point_t C,
1518  point_t D,
1519  struct dsp_specific *dsp,
1520  struct dsp_rpp *dsp_rpp)
1521 {
1522  int x, y;
1523 
1524 
1525 #ifdef FULL_DSP_DEBUGGING
1526  if (RT_G_DEBUG & DEBUG_HF) {
1527  VPRINT("\tA", A);
1528  VPRINT("\tB", B);
1529  VPRINT("\tC", C);
1530  VPRINT("\tD", D);
1531  }
1532 #endif
1533 
1534  switch (dsp->dsp_i.dsp_cuttype) {
1535  case DSP_CUT_DIR_llUR:
1536  return DSP_CUT_DIR_llUR;
1537  break;
1538 
1539  case DSP_CUT_DIR_ADAPT: {
1540  int lo[2], hi[2];
1541  point_t tmp;
1542  fastf_t h1, h2, h3, h4;
1543  fastf_t cAD, cBC; /* curvature in direction AD, and BC */
1544 
1545  if (RT_G_DEBUG & DEBUG_HF)
1546  bu_log("cell %d, %d adaptive triangulation... ",
1547  dsp_rpp->dsp_min[X],
1548  dsp_rpp->dsp_min[Y]);
1549 
1550  /*
1551  * We look at the points in the diagonal next cells to
1552  * determine the curvature along each diagonal of this
1553  * cell. This cell is divided into two triangles by
1554  * cutting across the cell in the direction of least
1555  * curvature.
1556  *
1557  * * * * *
1558  * \ /
1559  * \C D/
1560  * * *--* *
1561  * |\/|
1562  * |/\|
1563  * * *--* *
1564  * /A B\
1565  * / \
1566  * * * * *
1567  */
1568 
1569  lo[X] = dsp_rpp->dsp_min[X] - 1;
1570  lo[Y] = dsp_rpp->dsp_min[Y] - 1;
1571  hi[X] = dsp_rpp->dsp_max[X] + 1;
1572  hi[Y] = dsp_rpp->dsp_max[Y] + 1;
1573 
1574  /* a little bounds checking */
1575  if (lo[X] < 0) lo[X] = 0;
1576  if (lo[Y] < 0) lo[Y] = 0;
1577  if (hi[X] > dsp->xsiz)
1578  hi[X] = dsp->xsiz;
1579 
1580  if (hi[Y] > dsp->ysiz)
1581  hi[Y] = dsp->ysiz;
1582 
1583  /* compute curvature along the A->D direction */
1584  h1 = DSP(&dsp->dsp_i, lo[X], lo[Y]);
1585  h2 = A[Z];
1586  h3 = D[Z];
1587  h4 = DSP(&dsp->dsp_i, hi[X], hi[Y]);
1588 
1589  cAD = fabs(h3 + h1 - 2*h2) + fabs(h4 + h2 - 2*h3);
1590 
1591 
1592  /* compute curvature along the B->C direction */
1593  h1 = DSP(&dsp->dsp_i, hi[X], lo[Y]);
1594  h2 = B[Z];
1595  h3 = C[Z];
1596  h4 = DSP(&dsp->dsp_i, lo[X], hi[Y]);
1597 
1598  cBC = fabs(h3 + h1 - 2*h2) + fabs(h4 + h2 - 2*h3);
1599 
1600  if (cAD < cBC) {
1601  /* A-D cut is fine, no need to permute */
1602  if (RT_G_DEBUG & DEBUG_HF)
1603  bu_log("A-D cut\n");
1604 
1605  return DSP_CUT_DIR_llUR;
1606 
1607  }
1608 
1609  /* prefer the B-C cut */
1610  VMOVE(tmp, A);
1611  VMOVE(A, B);
1612  VMOVE(B, D);
1613  VMOVE(D, C);
1614  VMOVE(C, tmp);
1615  if (RT_G_DEBUG & DEBUG_HF)
1616  bu_log("B-C cut\n");
1617 
1618  return DSP_CUT_DIR_ULlr;
1619 
1620  break;
1621  }
1622  case DSP_CUT_DIR_ULlr:
1623  /* assign the values for the corner points
1624  *
1625  * D----C
1626  * | |
1627  * | |
1628  * | |
1629  * B----A
1630  */
1631  x = dsp_rpp->dsp_min[X];
1632  y = dsp_rpp->dsp_min[Y];
1633  VSET(B, x, y, DSP(&dsp->dsp_i, x, y));
1634 
1635  x = dsp_rpp->dsp_max[X];
1636  VSET(A, x, y, DSP(&dsp->dsp_i, x, y));
1637 
1638  y = dsp_rpp->dsp_max[Y];
1639  VSET(C, x, y, DSP(&dsp->dsp_i, x, y));
1640 
1641  x = dsp_rpp->dsp_min[X];
1642  VSET(D, x, y, DSP(&dsp->dsp_i, x, y));
1643 
1644  return DSP_CUT_DIR_ULlr;
1645  break;
1646  }
1647  bu_log("%s:%d Unknown DSP cut direction: %d\n",
1648  __FILE__, __LINE__, dsp->dsp_i.dsp_cuttype);
1649  bu_bomb("");
1650  /* not reached */
1651  return -1;
1652 }
1653 
1654 
1655 /**
1656  * determine if a point P is above/below the slope line on the
1657  * bounding box. e.g.:
1658  *
1659  * Bounding box side view:
1660  *
1661  * +-------+
1662  * | |
1663  * | * Determine if P (or Q) is above the
1664  * | /| diagonal from the two * points at the corners
1665  * | P. / | of the bounding box.
1666  * | / |
1667  * | / |
1668  * | / . |
1669  * | / Q |
1670  * |/ |
1671  * * |
1672  * | |
1673  * +-------+
1674  *
1675  * Return
1676  * 0 if pt above line (such as P in diagram)
1677  * 1 if pt at/below line (such as Q in diagram)
1678  */
1679 HIDDEN int
1680 check_bbpt_hit_elev(int i, /* indicates face of cell */
1681  point_t A,
1682  point_t B,
1683  point_t C,
1684  point_t D,
1685  point_t P)
1686 {
1687  fastf_t slope = 0.0;
1688  fastf_t delta = 0.0;
1689  fastf_t origin = 0.0;
1690 
1691 #ifdef FULL_DSP_DEBUGGING
1692  dlog("check_bbpt_hit_elev(");
1693 #endif
1694  switch (i) {
1695  case XMIN:
1696  /* the minimal YZ plane. Top view: * *
1697  * |
1698  * * *
1699  */
1700 #ifdef FULL_DSP_DEBUGGING
1701  dlog("XMIN)\n");
1702 #endif
1703  slope = C[Z] - A[Z];
1704  delta = P[Y] - A[Y];
1705  origin = A[Z];
1706  break;
1707  case XMAX:
1708  /* the maximal YZ plane. Top view: * *
1709  * |
1710  * * *
1711  */
1712 #ifdef FULL_DSP_DEBUGGING
1713  dlog("XMAX)\n");
1714 #endif
1715  slope = D[Z] - B[Z];
1716  delta = P[Y] - B[Y];
1717  origin = B[Z];
1718  break;
1719  case YMIN:
1720  /* the minimal XZ plane. Top view: * *
1721  *
1722  * * - *
1723  */
1724 #ifdef FULL_DSP_DEBUGGING
1725  dlog("YMIN)\n");
1726 #endif
1727  slope = B[Z] - A[Z];
1728  delta = P[X] - A[X];
1729  origin = A[Z];
1730  break;
1731  case YMAX:
1732  /* the maximal XZ plane. Top view: * - *
1733  *
1734  * * *
1735  */
1736 #ifdef FULL_DSP_DEBUGGING
1737  dlog("YMAX)\n");
1738 #endif
1739  slope = D[Z] - C[Z];
1740  delta = P[X] - C[X];
1741  origin = C[Z];
1742  break;
1743  case ZMIN:
1744 #ifdef FULL_DSP_DEBUGGING
1745  dlog("ZMIN)\n");
1746 #endif
1747  return 1;
1748  break;
1749  case ZMAX:
1750 #ifdef FULL_DSP_DEBUGGING
1751  dlog("ZMAX)\n");
1752 #endif
1753  return 0;
1754  break;
1755  default:
1756  bu_log("%s:%d Coding error, bad face %d\n", __FILE__, __LINE__, i);
1757  bu_bomb("");
1758  break;
1759  }
1760 
1761  if ((origin + slope * delta) < P[Z]) return 0;
1762 
1763  return 1;
1764 }
1765 
1766 
1767 /*
1768  * Return
1769  * 0 continue intersection calculations
1770  * 1 Terminate intersection computation
1771  */
1772 HIDDEN int
1774 {
1775  point_t A, B, C, D, P;
1776  int x, y;
1777  vect2d_t ab_first = V2INIT_ZERO;
1778  vect2d_t ab_second = V2INIT_ZERO;
1779  struct hit hits[4]; /* list of hits that are valid */
1780  struct hit *hitp;
1781  int hitf = 0; /* bit flags for valid hits in hits */
1782  int cond, i;
1783  int hitcount = 0;
1784  point_t bbmin, bbmax;
1785  fastf_t dot, dot2;
1786 
1787  for (x = 0; x < 4; x++)
1788  memset(hits+x, 0, sizeof(struct hit));
1789  x = 0;
1790 
1791  dlog("isect_ray_cell_top\n");
1792  DSP_BB_CK(dsp_bb);
1793 
1794  /* assign the values for the corner points
1795  *
1796  * C----D
1797  * | |
1798  * | |
1799  * | |
1800  * A----B
1801  */
1802  x = dsp_bb->dspb_rpp.dsp_min[X];
1803  y = dsp_bb->dspb_rpp.dsp_min[Y];
1804  VSET(A, x, y, DSP(&isect->dsp->dsp_i, x, y));
1805 
1806  x = dsp_bb->dspb_rpp.dsp_max[X];
1807  VSET(B, x, y, DSP(&isect->dsp->dsp_i, x, y));
1808 
1809  y = dsp_bb->dspb_rpp.dsp_max[Y];
1810  VSET(D, x, y, DSP(&isect->dsp->dsp_i, x, y));
1811 
1812  x = dsp_bb->dspb_rpp.dsp_min[X];
1813  VSET(C, x, y, DSP(&isect->dsp->dsp_i, x, y));
1814 
1815 
1816 #ifdef DEBUG_FULL
1817  if (RT_G_DEBUG & DEBUG_HF) {
1818  point_t p1, p2;
1819 
1820  VJOIN1(p1, isect->r.r_pt, isect->r.r_min, isect->r.r_dir);
1821  VMOVE(hits[0].hit_point, p1);
1822  hits[0].hit_dist = isect->r.r_min;
1823 
1824  VJOIN1(p2, isect->r.r_pt, isect->r.r_max, isect->r.r_dir);
1825  VMOVE(hits[1].hit_point, p2);
1826  hits[1].hit_dist = isect->r.r_max;
1827 
1828  plot_cell_top(isect, dsp_bb, A, B, C, D, hits, 3, 0);
1829  }
1830 #endif
1831 
1832 
1833  /* first order of business is to discard any "fake" hits on the
1834  * bounding box, and fill in any "real" hits in our list
1835  */
1836  VJOIN1(P, isect->r.r_pt, isect->r.r_min, isect->r.r_dir);
1837  if (check_bbpt_hit_elev(isect->dmin, A, B, C, D, P)) {
1838  hits[0].hit_dist = isect->r.r_min;
1839  VMOVE(hits[0].hit_point, P);
1840  VMOVE(hits[0].hit_normal, dsp_pl[isect->dmin]);
1841  /* vpriv */
1842  hits[0].hit_vpriv[X] = dsp_bb->dspb_rpp.dsp_min[X];
1843  hits[0].hit_vpriv[Y] = dsp_bb->dspb_rpp.dsp_min[Y];
1844  /* private */
1845  hits[0].hit_surfno = isect->dmin;
1846 
1847  hitcount++;
1848 
1849  hitf = 1;
1850  if (RT_G_DEBUG & DEBUG_HF) {
1851  dot = VDOT(hits[0].hit_normal, isect->r.r_dir);
1852  bu_log("hit ray/bb min Normal: %g %g %g %s\n",
1853  V3ARGS(hits[0].hit_normal),
1854  ((dot > 0.0) ? "outbound" : "inbound"));
1855  }
1856  } else {
1857  dlog("miss ray/bb min\n");
1858  }
1859 
1860 
1861  /* make sure the point P is below the cell top */
1862  VJOIN1(P, isect->r.r_pt, isect->r.r_max, isect->r.r_dir);
1863  if (check_bbpt_hit_elev(isect->dmax, A, B, C, D, P)) {
1864  /* P is at or below the top surface */
1865  hits[3].hit_dist = isect->r.r_max;
1866  VMOVE(hits[3].hit_point, P);
1867  VMOVE(hits[3].hit_normal, dsp_pl[isect->dmax]);
1868  /* vpriv */
1869  hits[3].hit_vpriv[X] = dsp_bb->dspb_rpp.dsp_min[X];
1870  hits[3].hit_vpriv[Y] = dsp_bb->dspb_rpp.dsp_min[Y];
1871  /* private */
1872  hits[3].hit_surfno = isect->dmax;
1873 
1874  hitcount++;
1875 
1876  hitf |= 8;
1877  if (RT_G_DEBUG & DEBUG_HF) {
1878  dot = VDOT(hits[3].hit_normal, isect->r.r_dir);
1879  bu_log("hit ray/bb max Normal: %g %g %g %s\n",
1880  V3ARGS(hits[3].hit_normal),
1881  ((dot > 0.0) ? "outbound" : "inbound"));
1882  }
1883  } else {
1884  dlog("miss ray/bb max\n");
1885  }
1886 
1887 
1888  (void)permute_cell(A, B, C, D, isect->dsp, &dsp_bb->dspb_rpp);
1889 
1890  if ((cond=isect_ray_triangle(isect, B, D, A, &hits[1], ab_first)) > 0.0) {
1891  /* hit triangle */
1892 
1893  /* record cell */
1894  hits[1].hit_vpriv[X] = dsp_bb->dspb_rpp.dsp_min[X];
1895  hits[1].hit_vpriv[Y] = dsp_bb->dspb_rpp.dsp_min[Y];
1896  hits[1].hit_surfno = ZTOP; /* indicate we hit the top */
1897 
1898  hitcount++;
1899  hitf |= 2;
1900  dlog(" hit triangle 1 (alpha: %g beta:%g alpha+beta: %g) vpriv %g %g\n",
1901  ab_first[0], ab_first[1], ab_first[0] + ab_first[1],
1902  hits[1].hit_vpriv[X], hits[1].hit_vpriv[Y]);
1903  } else {
1904  dlog(" miss triangle 1 (alpha: %g beta:%g a+b: %g) cond:%d\n",
1905  ab_first[0], ab_first[1], ab_first[0] + ab_first[1], cond);
1906  }
1907  if ((cond=isect_ray_triangle(isect, C, A, D, &hits[2], ab_second)) > 0.0) {
1908  /* hit triangle */
1909 
1910  /* record cell */
1911  hits[2].hit_vpriv[X] = dsp_bb->dspb_rpp.dsp_min[X];
1912  hits[2].hit_vpriv[Y] = dsp_bb->dspb_rpp.dsp_min[Y];
1913  hits[2].hit_surfno = ZTOP; /* indicate we hit the top */
1914 
1915  hitcount++;
1916 
1917  hitf |= 4;
1918  dlog(" hit triangle 2 (alpha: %g beta:%g alpha+beta: %g) vpriv %g %g\n",
1919  ab_second[0], ab_second[1], ab_second[0] + ab_second[1],
1920  hits[2].hit_vpriv[X], hits[2].hit_vpriv[Y]);
1921 
1922  if (hitf & 2) {
1923  /* if this hit occurs before the hit on the other triangle
1924  * swap the order
1925  */
1926  if (hits[1].hit_dist > hits[2].hit_dist) {
1927  struct hit tmp;
1928  tmp = hits[1]; /* struct copy */
1929  hits[1] = hits[2]; /* struct copy */
1930  hits[2] = tmp; /* struct copy */
1931  dlog("re-ordered triangle hits\n");
1932 
1933  } else
1934  dlog("triangle hits in order\n");
1935  }
1936 
1937 
1938  } else {
1939  dlog(" miss triangle 2 (alpha: %g beta:%g alpha+beta: %g) cond:%d\n",
1940  ab_second[0], ab_second[1], ab_second[0] + ab_second[1], cond);
1941  }
1942 
1943  if (RT_G_DEBUG & DEBUG_HF) {
1944  bu_log("hitcount: %d flags: 0x%0x\n", hitcount, hitf);
1945 
1946  plot_cell_top(isect, dsp_bb, A, B, C, D, hits, hitf, 1);
1947  for (i = 0; i < 4; i++) {
1948  if (hitf & (1<<i)) {
1949  fastf_t v = VDOT(isect->r.r_dir, hits[i].hit_normal);
1950 
1951  bu_log("%d dist:%g N:%g %g %g ",
1952  i, hits[i].hit_dist, V3ARGS(hits[i].hit_normal));
1953 
1954  if (v > 0.0) bu_log("outbound\n");
1955  else if (v < 0.0) bu_log("inbound\n");
1956  else bu_log("perp\n");
1957  }
1958  }
1959  bu_log("assembling segs\n");
1960  }
1961 
1962 
1963  /* fill out the segment structures */
1964 
1965  hitp = 0;
1966  for (i = 0; i < 4; i++) {
1967  if (hitf & (1<<i)) {
1968  if (hitp) {
1969 
1970  dot2 = VDOT(isect->r.r_dir, hits[i].hit_normal);
1971 
1972  /* if we have two entry points then pick the first one */
1973  if (dot2 < 0.0) {
1974  dlog("dot2(%g) < 0.0\n", dot2);
1975  if (hitp->hit_dist > hits[i].hit_dist) {
1976  dlog("skipping duplicate entry point at dist %g\n",
1977  hitp->hit_dist);
1978 
1979  hitp = &hits[i];
1980  } else {
1981  dlog("skipping duplicate entry point at dist %g\n",
1982  hits[i].hit_dist);
1983  }
1984 
1985  continue;
1986  }
1987 
1988  /* int/float conv */
1989  VMOVE(bbmin, dsp_bb->dspb_rpp.dsp_min);
1990  VMOVE(bbmax, dsp_bb->dspb_rpp.dsp_max);
1991 
1992  /* create seg with hits[i].hit_point as out point */
1993  if (add_seg(isect, hitp, &hits[i], bbmin, bbmax, 255, 255, 255))
1994  return 1;
1995 
1996  hitp = 0;
1997  } else {
1998  dot = VDOT(isect->r.r_dir, hits[i].hit_normal);
1999  if (dot >= 0.0)
2000  continue;
2001 
2002  /* remember hits[i].hit_point); as in point */
2003  if (RT_G_DEBUG & DEBUG_HF) {
2004  bu_log("in-hit at dist %g\n", hits[i].hit_dist);
2005  }
2006  hitp = &hits[i];
2007  }
2008  }
2009  }
2010 
2011  if (hitp && hitcount > 1) {
2012  point_t p1, p2;
2013 
2014  bu_log("----------------ERROR incomplete segment-------------\n");
2015  bu_log(" pixel %d %d\n", isect->ap->a_x, isect->ap->a_y);
2016 
2017  VJOIN1(p1, isect->r.r_pt, isect->r.r_min, isect->r.r_dir);
2018  VMOVE(hits[0].hit_point, p1);
2019  hits[0].hit_dist = isect->r.r_min;
2020 
2021  VJOIN1(p2, isect->r.r_pt, isect->r.r_max, isect->r.r_dir);
2022  VMOVE(hits[1].hit_point, p2);
2023  hits[1].hit_dist = isect->r.r_max;
2024 
2025  if (RT_G_DEBUG & DEBUG_HF)
2026  plot_cell_top(isect, dsp_bb, A, B, C, D, hits, 3, 0);
2027  }
2028  return 0;
2029 }
2030 
2031 
2032 /**
2033  * Compute the intersections of a ray with a rectangular parallelepiped
2034  * (RPP) that has faces parallel to the coordinate planes
2035  *
2036  * The algorithm here was developed by Gary Kuehl for GIFT. A good
2037  * description of the approach used can be found in "??" by XYZZY and
2038  * Barsky, ACM Transactions on Graphics, Vol 3 No 1, January 1984.
2039  *
2040  * Note - The computation of entry and exit distance is mandatory, as
2041  * the final test catches the majority of misses.
2042  *
2043  * Note - A hit is returned if the intersect is behind the start
2044  * point.
2045  *
2046  * Returns -
2047  * 0 if ray does not hit RPP,
2048  * !0 if ray hits RPP.
2049  *
2050  * Implicit return -
2051  * rp->r_min = dist from start of ray to point at which ray ENTERS solid
2052  * rp->r_max = dist from start of ray to point at which ray LEAVES solid
2053  * isect->dmin
2054  * isect->dmax
2055  */
2056 int
2057 dsp_in_rpp(struct isect_stuff *isect,
2058  register const fastf_t *min,
2059  register const fastf_t *max)
2060 {
2061  struct xray *rp = &isect->r;
2062  /* inverses of rp->r_dir[] */
2063  register const fastf_t *invdir = isect->inv_dir;
2064  register const fastf_t *pt = &rp->r_pt[0];
2065  register fastf_t sv;
2066  register fastf_t rmin = -INFINITY;
2067  register fastf_t rmax = INFINITY;
2068  int dmin = -1;
2069  int dmax = -1;
2070 
2071  /* Start with infinite ray, and trim it down */
2072 
2073  /* X axis */
2074  if (*invdir < 0.0) {
2075  /* Heading towards smaller numbers */
2076  /* if (*min > *pt) miss */
2077  if (rmax > (sv = (*min - *pt) * *invdir)) {
2078  rmax = sv;
2079  dmax = XMIN;
2080  }
2081  if (rmin < (sv = (*max - *pt) * *invdir)) {
2082  rmin = sv;
2083  dmin = XMAX;
2084  }
2085  } else if (*invdir > 0.0) {
2086  /* Heading towards larger numbers */
2087  /* if (*max < *pt) miss */
2088  if (rmax > (sv = (*max - *pt) * *invdir)) {
2089  rmax = sv;
2090  dmax = XMAX;
2091  }
2092  if (rmin < ((sv = (*min - *pt) * *invdir))) {
2093  rmin = sv;
2094  dmin = XMIN;
2095  }
2096  } else {
2097  /*
2098  * Direction cosines along this axis is NEAR 0, which implies
2099  * that the ray is perpendicular to the axis, so merely check
2100  * position against the boundaries.
2101  */
2102  if ((*min > *pt) || (*max < *pt))
2103  return 0; /* MISS */
2104  }
2105 
2106  /* Y axis */
2107  pt++; invdir++; max++; min++;
2108  if (*invdir < 0.0) {
2109  if (rmax > (sv = (*min - *pt) * *invdir)) {
2110  /* towards smaller */
2111  rmax = sv;
2112  dmax = YMIN;
2113  }
2114  if (rmin < (sv = (*max - *pt) * *invdir)) {
2115  rmin = sv;
2116  dmin = YMAX;
2117  }
2118  } else if (*invdir > 0.0) {
2119  /* towards larger */
2120  if (rmax > (sv = (*max - *pt) * *invdir)) {
2121  rmax = sv;
2122  dmax = YMAX;
2123  }
2124  if (rmin < ((sv = (*min - *pt) * *invdir))) {
2125  rmin = sv;
2126  dmin = YMIN;
2127  }
2128  } else {
2129  if ((*min > *pt) || (*max < *pt))
2130  return 0; /* MISS */
2131  }
2132 
2133  /* Z axis */
2134  pt++; invdir++; max++; min++;
2135  if (*invdir < 0.0) {
2136  /* towards smaller */
2137  if (rmax > (sv = (*min - *pt) * *invdir)) {
2138  rmax = sv;
2139  dmax = ZMIN;
2140  }
2141  if (rmin < (sv = (*max - *pt) * *invdir)) {
2142  rmin = sv;
2143  dmin = ZMAX;
2144  }
2145  } else if (*invdir > 0.0) {
2146  /* towards larger */
2147  if (rmax > (sv = (*max - *pt) * *invdir)) {
2148  rmax = sv;
2149  dmax = ZMAX;
2150  }
2151  if (rmin < ((sv = (*min - *pt) * *invdir))) {
2152  rmin = sv;
2153  dmin = ZMIN;
2154  }
2155  } else {
2156  if ((*min > *pt) || (*max < *pt))
2157  return 0; /* MISS */
2158  }
2159 
2160  /* If equal, RPP is actually a plane */
2161  if (rmin > rmax)
2162  return 0; /* MISS */
2163 
2164  /* HIT. Only now do rp->r_min and rp->r_max have to be written */
2165  rp->r_min = rmin;
2166  rp->r_max = rmax;
2167 
2168  isect->dmin = dmin;
2169  isect->dmax = dmax;
2170  return 1; /* HIT */
2171 }
2172 
2173 
2174 #ifdef ORDERED_ISECT
2175 HIDDEN int
2176 isect_ray_dsp_bb(struct isect_stuff *isect, struct dsp_bb *dsp_bb);
2177 
2178 
2179 /**
2180  * Return
2181  * 0 continue intersection calculations
2182  * 1 Terminate intersection computation
2183  */
2184 HIDDEN int
2186  struct dsp_bb *dsp_bb,
2187  point_t minpt, /* entry point of dsp_bb */
2188  point_t UNUSED(maxpt), /* exit point of dsp_bb */
2189  point_t bbmin, /* min point of bb (Z=0) */
2190  point_t UNUSED(bbmax)) /* max point of bb */
2191 {
2192  fastf_t tDX; /* dist along ray to span 1 cell in X dir */
2193  fastf_t tDY; /* dist along ray to span 1 cell in Y dir */
2194  fastf_t tX, tY; /* dist from hit pt. to next cell boundary */
2195  fastf_t curr_dist;
2196  short cX, cY; /* coordinates of current cell */
2197  short cs; /* cell X, Y dimension */
2198  short stepX, stepY; /* dist to step in child array for each dir */
2199  short stepPX, stepPY;
2200  fastf_t out_dist;
2201  struct dsp_bb **p;
2202  fastf_t *stom = &isect->dsp->dsp_i.dsp_stom[0];
2203  point_t pt, v;
2204  int loop = 0;
2205 
2206  DSP_BB_CK(dsp_bb);
2207 
2208  /* compute the size of a cell in each direction */
2209  cs = dsp_bb->dspb_subcell_size;
2210 
2211  /* compute current cell */
2212  cX = (minpt[X] - bbmin[X]) / cs;
2213  cY = (minpt[Y] - bbmin[Y]) / cs;
2214 
2215  /* a little bounds checking because a hit on XMAX or YMAX looks
2216  * like it should be in the next cell outside the box
2217  */
2218  if (cX >= dsp_bb->dspb_ch_dim[X]) cX = dsp_bb->dspb_ch_dim[X] - 1;
2219  if (cY >= dsp_bb->dspb_ch_dim[Y]) cY = dsp_bb->dspb_ch_dim[Y] - 1;
2220 
2221 #ifdef FULL_DSP_DEBUGGING
2222  dlog("recurse_dsp_bb cell size: %d current cell: %d %d\n",
2223  cs, cX, cY);
2224  dlog("dspb_ch_dim x:%d y:%d\n",
2225  dsp_bb->dspb_ch_dim[X], dsp_bb->dspb_ch_dim[Y]);
2226 #endif
2227 
2228  tX = tY = curr_dist = isect->r.r_min;
2229 
2230  if (isect->r.r_dir[X] < 0.0) {
2231  stepPX = stepX = -1;
2232  /* tDX is the distance along the ray we have to travel to
2233  * traverse a cell (travel a unit distance) along the X axis
2234  * of the grid
2235  */
2236  tDX = -cs / isect->r.r_dir[X];
2237 
2238  /* tX is the distance along the ray to the first cell boundary
2239  * in the X direction beyond our hit point (minpt)
2240  */
2241  tX += ((bbmin[X] + (cX * cs)) - minpt[X]) / isect->r.r_dir[X];
2242  } else {
2243  stepPX = stepX = 1;
2244  tDX = cs / isect->r.r_dir[X];
2245 
2246  if (isect->r.r_dir[X] > 0.0)
2247  tX += ((bbmin[X] + ((cX+1) * cs)) - minpt[X]) / isect->r.r_dir[X];
2248  else
2249  tX = MAX_FASTF; /* infinite distance to next X boundary */
2250  }
2251 
2252  if (isect->r.r_dir[Y] < 0) {
2253  /* distance in dspb_children we have to move to step in Y dir
2254  */
2255  stepY = -1;
2256  stepPY = -dsp_bb->dspb_ch_dim[X];
2257  tDY = -cs / isect->r.r_dir[Y];
2258  tY += ((bbmin[Y] + (cY * cs)) - minpt[Y]) / isect->r.r_dir[Y];
2259  } else {
2260  stepY = 1;
2261  stepPY = dsp_bb->dspb_ch_dim[X];
2262  tDY = cs / isect->r.r_dir[Y];
2263 
2264  if (isect->r.r_dir[Y] > 0.0)
2265  tY += ((bbmin[Y] + ((cY+1) * cs)) - minpt[Y]) / isect->r.r_dir[Y];
2266  else
2267  tY = MAX_FASTF;
2268  }
2269 
2270  /* factor in the tolerance to the out-distance */
2271  out_dist = isect->r.r_max - isect->tol->dist;
2272 
2273  p = &dsp_bb->dspb_children[dsp_bb->dspb_ch_dim[X] * cY + cX];
2274 #ifdef FULL_DSP_DEBUGGING
2275  dlog("tX:%g tY:%g\n", tX, tY);
2276 
2277 #endif
2278 
2279  do {
2280  /* intersect with the current cell */
2281  if (RT_G_DEBUG & DEBUG_HF) {
2282  if (loop)
2283  bu_log("\nisect sub-cell %d %d curr_dist:%g out_dist: %g",
2284  cX, cY, curr_dist, out_dist);
2285  else {
2286  bu_log("isect sub-cell %d %d curr_dist:%g out_dist %g",
2287  cX, cY, curr_dist, out_dist);
2288  loop = 1;
2289  }
2290  VJOIN1(pt, isect->r.r_pt, curr_dist, isect->r.r_dir);
2291  MAT4X3PNT(v, stom, pt);
2292  bu_log("pt %g %g %g\n", V3ARGS(v));
2293  }
2294 
2295  /* get pointer to the current child cell. Note the extra
2296  * level of indirection. We want to march p through
2297  * dspb_children using pointer addition, but dspb_children
2298  * contains pointers. We need to be sure to increment the
2299  * offset in the array, not the child pointer.
2300  */
2301 
2302  if (RT_G_DEBUG & DEBUG_HF)
2304 
2305  if (isect_ray_dsp_bb(isect, *p)) return 1;
2306 
2307  if (RT_G_DEBUG & DEBUG_HF)
2308  bu_log_indent_delta(-4);
2309 
2310  /* figure out which cell is next */
2311  if (tX < tY) {
2312  cX += stepX; /* track cell offset for debugging */
2313  p += stepPX;
2314 #ifdef FULL_DSP_DEBUGGING
2315  dlog("stepping X to %d because %g < %g\n", cX, tX, tY);
2316 #endif
2317  curr_dist = tX;
2318  tX += tDX;
2319  } else {
2320  cY += stepY; /* track cell offset for debugging */
2321  p += stepPY;
2322 #ifdef FULL_DSP_DEBUGGING
2323  dlog("stepping Y to %d because %g >= %g\n", cY, tX, tY);
2324 #endif
2325  curr_dist = tY;
2326  tY += tDY;
2327  }
2328 #ifdef FULL_DSP_DEBUGGING
2329  dlog("curr_dist %g, out_dist %g\n", curr_dist, out_dist);
2330 #endif
2331  } while (curr_dist < out_dist &&
2332  cX < dsp_bb->dspb_ch_dim[X] && cX >= 0 &&
2333  cY < dsp_bb->dspb_ch_dim[Y] && cY >= 0);
2334 
2335  return 0;
2336 }
2337 
2338 
2339 #endif
2340 
2341 /**
2342  * Intersect a ray with a DSP bounding box. This is the primary child
2343  * of rt_dsp_shot()
2344  *
2345  * Return
2346  * 0 continue intersection calculations
2347  * 1 Terminate intersection computation
2348  */
2349 HIDDEN int
2350 isect_ray_dsp_bb(struct isect_stuff *isect, struct dsp_bb *dsp_bb)
2351 {
2352  point_t bbmin, bbmax;
2353  point_t minpt, maxpt;
2354  fastf_t min_z;
2355  /* the rest of these support debugging output */
2356  FILE *fp;
2357  static int plotnum;
2358  fastf_t *stom;
2359  point_t pt;
2360  struct xray *r = &isect->r; /* Does this buy us anything? */
2361 
2362  if (dsp_bb) {
2363  DSP_BB_CK(dsp_bb);
2364  } else {
2365  bu_log("%s:%d null ptr pixel %d %d\n",
2366  __FILE__, __LINE__, isect->ap->a_x, isect->ap->a_y);
2367  bu_bomb("");
2368  }
2369 
2370  if (RT_G_DEBUG & DEBUG_HF) {
2371  bu_log("\nisect_ray_dsp_bb((%d, %d, %d) (%d, %d, %d))\n",
2372  V3ARGS(dsp_bb->dspb_rpp.dsp_min),
2373  V3ARGS(dsp_bb->dspb_rpp.dsp_max));
2374  }
2375 
2376  /* check to see if we miss the RPP for this area entirely */
2377  VMOVE(bbmax, dsp_bb->dspb_rpp.dsp_max);
2378  VSET(bbmin,
2379  dsp_bb->dspb_rpp.dsp_min[X],
2380  dsp_bb->dspb_rpp.dsp_min[Y], 0.0);
2381 
2382 
2383  if (! dsp_in_rpp(isect, bbmin, bbmax)) {
2384  /* missed it all, just return */
2385 
2386  if (RT_G_DEBUG & DEBUG_HF) {
2387  bu_log("missed... ");
2388  fclose(draw_dsp_bb(&plotnum, dsp_bb, isect->dsp, 0, 150, 0));
2389  }
2390 
2391  return 0;
2392  }
2393 
2394  /* At this point we know that we've hit the overall bounding box
2395  */
2396 
2397  VJOIN1(minpt, r->r_pt, r->r_min, r->r_dir);
2398  VJOIN1(maxpt, r->r_pt, r->r_max, r->r_dir);
2399  /* We could do the following:
2400  *
2401  * however, we only need the Z values now, and benchmarking show
2402  * that this is an expensive
2403  */
2404  if (RT_G_DEBUG & DEBUG_HF) {
2405 
2406  stom = &isect->dsp->dsp_i.dsp_stom[0];
2407 
2408  bu_log("hit b-box ");
2409  fp = draw_dsp_bb(&plotnum, dsp_bb, isect->dsp, 200, 200, 100);
2410 
2411  pl_color(fp, 150, 150, 255);
2412  MAT4X3PNT(pt, stom, minpt);
2413  pdv_3move(fp, pt);
2414  MAT4X3PNT(pt, stom, maxpt);
2415  pdv_3cont(fp, pt);
2416 
2417  fclose(fp);
2418  }
2419 
2420 
2421  /* if both hits are UNDER the top of the "foundation" pillar, we
2422  * can just add a segment for that range and return
2423  */
2424  min_z = dsp_bb->dspb_rpp.dsp_min[Z];
2425 
2426  if (minpt[Z] < min_z && maxpt[Z] < min_z) {
2427  /* add hit segment */
2428  struct hit seg_in, seg_out;
2429  VSETALL(seg_in.hit_vpriv, 0.0);
2430  VSETALL(seg_out.hit_vpriv, 0.0);
2431 
2432  seg_in.hit_magic = RT_HIT_MAGIC;
2433  seg_in.hit_dist = r->r_min;
2434  VMOVE(seg_in.hit_point, minpt);
2435  VMOVE(seg_in.hit_normal, dsp_pl[isect->dmin]);
2436  /* hit_priv */
2437  /* hit_private */
2438  seg_in.hit_surfno = isect->dmin;
2439  /* hit_rayp */
2440 
2441  seg_out.hit_dist = r->r_max;
2442 
2443  VMOVE(seg_out.hit_point, maxpt);
2444  VMOVE(seg_out.hit_normal, dsp_pl[isect->dmax]);
2445  /* hit_priv */
2446  /* hit_private */
2447  seg_out.hit_surfno = isect->dmax;
2448  /* hit_rayp */
2449 
2450  if (RT_G_DEBUG & DEBUG_HF) {
2451  /* we need these for debug output
2452  * VMOVE(seg_in.hit_point, minpt);
2453  * VMOVE(seg_out.hit_point, maxpt);
2454  */
2455  /* create a special bounding box for plotting purposes */
2456  VMOVE(bbmax, dsp_bb->dspb_rpp.dsp_max);
2457  VMOVE(bbmin, dsp_bb->dspb_rpp.dsp_min);
2458  bbmax[Z] = bbmin[Z];
2459  bbmin[Z] = 0.0;
2460  }
2461 
2462  /* outta here */
2463  return add_seg(isect, &seg_in, &seg_out, bbmin, bbmax, 0, 255, 255);
2464  }
2465 
2466 
2467  /* We've hit something where we might be going through the
2468  * boundary. We've got to intersect the children
2469  */
2470  if (dsp_bb->dspb_ch_dim[0]) {
2471 #ifdef ORDERED_ISECT
2472  return recurse_dsp_bb(isect, dsp_bb, minpt, maxpt, bbmin, bbmax);
2473 #else
2474  int i;
2475  /* there are children, so we recurse */
2476  i = dsp_bb->dspb_ch_dim[X] * dsp_bb->dspb_ch_dim[Y] - 1;
2477  if (RT_G_DEBUG & DEBUG_HF)
2479 
2480  for (; i >= 0; i--)
2481  isect_ray_dsp_bb(isect, dsp_bb->dspb_children[i]);
2482 
2483  if (RT_G_DEBUG & DEBUG_HF)
2484  bu_log_indent_delta(-4);
2485 
2486  return 0;
2487 
2488 #endif
2489  }
2490 
2491  /***********************************************************************
2492  *
2493  * This section is for level 0 intersections only
2494  *
2495  ***********************************************************************/
2496 
2497  /* intersect the DSP grid surface geometry */
2498 
2499  /* Check for a hit on the triangulated zone on top. This gives us
2500  * intersections on the triangulated top, and the sides and bottom
2501  * of the bounding box for the triangles.
2502  *
2503  * We do this first because we already know that the ray does NOT
2504  * just pass through the "foundation " pillar underneath (see test
2505  * above)
2506  */
2507  bbmin[Z] = dsp_bb->dspb_rpp.dsp_min[Z];
2508  if (dsp_in_rpp(isect, bbmin, bbmax)) {
2509  /* hit rpp */
2510 
2511  isect_ray_cell_top(isect, dsp_bb);
2512  }
2513 
2514 
2515  /* check for hits on the "foundation" pillar under the top. The
2516  * ray may have entered through the top of the pillar, possibly
2517  * after having come down through the triangles above
2518  */
2519  bbmax[Z] = dsp_bb->dspb_rpp.dsp_min[Z];
2520  bbmin[Z] = 0.0;
2521  if (dsp_in_rpp(isect, bbmin, bbmax)) {
2522  /* hit rpp */
2523  struct hit in_hit, out_hit;
2524  VSETALL(in_hit.hit_vpriv, 0.0);
2525  VSETALL(out_hit.hit_vpriv, 0.0);
2526 
2527  VJOIN1(minpt, r->r_pt, r->r_min, r->r_dir);
2528  VJOIN1(maxpt, r->r_pt, r->r_max, r->r_dir);
2529 
2530  in_hit.hit_dist = r->r_min;
2531  in_hit.hit_surfno = isect->dmin;
2532  VMOVE(in_hit.hit_point, minpt);
2533  VMOVE(in_hit.hit_normal, dsp_pl[isect->dmin]);
2534 
2535  out_hit.hit_dist = r->r_max;
2536  out_hit.hit_surfno = isect->dmax;
2537  VMOVE(out_hit.hit_point, maxpt);
2538  VMOVE(out_hit.hit_normal, dsp_pl[isect->dmax]);
2539 
2540  /* add a segment to the list */
2541  return add_seg(isect, &in_hit, &out_hit, bbmin, bbmax, 255, 255, 0);
2542  }
2543 
2544  return 0;
2545 }
2546 
2547 
2548 /**
2549  * Intersect a ray with a dsp.
2550  * If an intersection occurs, a struct seg will be acquired
2551  * and filled in.
2552  *
2553  * Returns -
2554  * 0 MISS
2555  * >0 HIT
2556  */
2557 int
2558 rt_dsp_shot(struct soltab *stp, register struct xray *rp, struct application *ap, struct seg *seghead)
2559 {
2560  register struct dsp_specific *dsp =
2561  (struct dsp_specific *)stp->st_specific;
2562  register struct seg *segp;
2563  int i;
2564  vect_t dir; /* temp storage */
2565  vect_t v;
2566  struct isect_stuff isect;
2567  fastf_t delta;
2568 
2569  RT_DSP_CK_MAGIC(dsp);
2570  BU_CK_VLS(&dsp->dsp_i.dsp_name);
2571 
2572  switch (dsp->dsp_i.dsp_datasrc) {
2573  case RT_DSP_SRC_V4_FILE:
2574  case RT_DSP_SRC_FILE:
2575  BU_CK_MAPPED_FILE(dsp->dsp_i.dsp_mp);
2576  break;
2577  case RT_DSP_SRC_OBJ:
2578  RT_CK_DB_INTERNAL(dsp->dsp_i.dsp_bip);
2579  RT_CK_BINUNIF(dsp->dsp_i.dsp_bip->idb_ptr);
2580  break;
2581  }
2582 
2583  /*
2584  * map ray into the coordinate system of the dsp
2585  */
2586  MAT4X3PNT(isect.r.r_pt, dsp->dsp_i.dsp_mtos, rp->r_pt);
2587  MAT4X3VEC(dir, dsp->dsp_i.dsp_mtos, rp->r_dir);
2588  VMOVE(isect.r.r_dir, dir);
2589  VUNITIZE(isect.r.r_dir);
2590 
2591  /* wrap a bunch of things together */
2592  isect.ap = ap;
2593  isect.stp = stp;
2594  isect.dsp = (struct dsp_specific *)stp->st_specific;
2595  isect.tol = &ap->a_rt_i->rti_tol;
2596  isect.num_segs = 0;
2597 
2598  VINVDIR(isect.inv_dir, isect.r.r_dir);
2599  BU_LIST_INIT(&isect.seglist);
2600 
2601 
2602  if (RT_G_DEBUG & DEBUG_HF) {
2603  bu_log("rt_dsp_shot(pt:(%g %g %g)\n\tdir[%g]:(%g %g %g))\n pixel(%d, %d)\n",
2604  V3ARGS(rp->r_pt),
2605  MAGNITUDE(rp->r_dir),
2606  V3ARGS(rp->r_dir),
2607  ap->a_x, ap->a_y);
2608 
2609  bn_mat_print("mtos", dsp->dsp_i.dsp_mtos);
2610  bu_log("Solid space ray pt:(%g %g %g)\n", V3ARGS(isect.r.r_pt));
2611  bu_log("\tdir[%g]: [%g %g %g]\n\tunit_dir(%g %g %g)\n",
2612  MAGNITUDE(dir),
2613  V3ARGS(dir),
2614  V3ARGS(isect.r.r_dir));
2615  }
2616 
2617  /* We look at the topmost layer of the bounding-box tree and make
2618  * sure that it has dimension 1. Otherwise, something is wrong
2619  */
2620  if (isect.dsp->layer[isect.dsp->layers-1].dim[X] != 1 ||
2621  isect.dsp->layer[isect.dsp->layers-1].dim[Y] != 1) {
2622  bu_log("%s:%d how do i find the topmost layer?\n",
2623  __FILE__, __LINE__);
2624  bu_bomb("");
2625  }
2626 
2627 
2628  /* intersect the ray with the bounding rpps */
2629  (void)isect_ray_dsp_bb(&isect, isect.dsp->layer[isect.dsp->layers-1].p);
2630 
2631  /* if we missed it all, give up now */
2632  if (BU_LIST_IS_EMPTY(&isect.seglist))
2633  return 0;
2634 
2635  /* map hit distances back to model space */
2636  i = 0;
2637  for (BU_LIST_FOR(segp, seg, &isect.seglist)) {
2638  i += 2;
2639  if (RT_G_DEBUG & DEBUG_HF) {
2640  bu_log("\nsolid in:%6g out:%6g\t",
2641  segp->seg_in.hit_dist,
2642  segp->seg_out.hit_dist);
2643  }
2644 
2645  /* form vector for hit point from start in solid space */
2646  VSCALE(dir, isect.r.r_dir, segp->seg_in.hit_dist);
2647  /* transform vector into model space */
2648  MAT4X3VEC(v, dsp->dsp_i.dsp_stom, dir);
2649  /* get magnitude */
2650  segp->seg_in.hit_dist = MAGNITUDE(v);
2651  /* XXX why is this necessary? */
2652  if (VDOT(v, rp->r_dir) < 0.0) segp->seg_in.hit_dist *= -1.0;
2653 
2654  VSCALE(dir, isect.r.r_dir, segp->seg_out.hit_dist);
2655  MAT4X3VEC(v, dsp->dsp_i.dsp_stom, dir);
2656  segp->seg_out.hit_dist = MAGNITUDE(v);
2657  if (VDOT(v, rp->r_dir) < 0.0) segp->seg_out.hit_dist *= -1.0;
2658 
2659 
2660  delta = segp->seg_out.hit_dist - segp->seg_in.hit_dist;
2661 
2662  if (delta < 0.0 && !NEAR_ZERO(delta, ap->a_rt_i->rti_tol.dist)) {
2663  bu_log("Pixel %d %d seg inside out in:%g out:%g seg_len:%g\n",
2664  ap->a_x, ap->a_y, segp->seg_in.hit_dist, segp->seg_out.hit_dist,
2665  delta);
2666  }
2667 
2668  if (RT_G_DEBUG & DEBUG_HF) {
2669  bu_log("model in:%6g out:%6g\t",
2670  segp->seg_in.hit_dist,
2671  segp->seg_out.hit_dist);
2672  }
2673  }
2674 
2675  if (RT_G_DEBUG & DEBUG_HF) {
2676  fastf_t NdotD;
2677  fastf_t d;
2678  static const plane_t plane = {0.0, 0.0, -1.0, 0.0};
2679 
2680  NdotD = VDOT(plane, rp->r_dir);
2681  d = - ((VDOT(plane, rp->r_pt) - plane[H]) / NdotD);
2682  bu_log("rp -> Z=0 dist: %g\n", d);
2683  }
2684 
2685  /* transfer list of hitpoints */
2686  BU_LIST_APPEND_LIST(&(seghead->l), &isect.seglist);
2687 
2688  return i;
2689 }
2690 
2691 
2692 /***********************************************************************
2693  *
2694  * Compute the model-space normal at a gridpoint
2695  *
2696  */
2697 HIDDEN void
2699  struct dsp_specific *dsp,
2700  unsigned int x,
2701  unsigned int y,
2702  FILE *fd,
2703  fastf_t len)
2704 {
2705  /* Gridpoint specified is "B" we compute normal by taking the
2706  * cross product of the vectors A->C, D->E
2707  *
2708  * E
2709  *
2710  * |
2711  *
2712  * A - B - C
2713  *
2714  * |
2715  *
2716  * D
2717  */
2718 
2719  point_t A, C, D, E, tmp, pt, endpt;
2720  vect_t Vac, Vde;
2721 
2722  if (RT_G_DEBUG & DEBUG_HF) {
2723  bu_log("normal at %d %d\n", x, y);
2724  bn_mat_print("\tstom", dsp->dsp_i.dsp_stom);
2725  }
2726  VSET(tmp, x, y, DSP(&dsp->dsp_i, x, y));
2727 
2728  if (x == 0) {
2729  VMOVE(A, tmp);
2730  } else {
2731  VSET(A, x-1, y, DSP(&dsp->dsp_i, x-1, y));
2732  }
2733 
2734  if (x >= XSIZ(dsp)) {
2735  VMOVE(C, tmp);
2736  } else {
2737  VSET(C, x+1, y, DSP(&dsp->dsp_i, x+1, y));
2738  }
2739 
2740  if (y == 0) {
2741  VMOVE(D, tmp);
2742  } else {
2743  VSET(D, x, y-1, DSP(&dsp->dsp_i, x, y-1));
2744  }
2745 
2746  if (y >= YSIZ(dsp)) {
2747  VMOVE(E, tmp);
2748  } else {
2749  VSET(E, x, y+1, DSP(&dsp->dsp_i, x, y+1));
2750  }
2751 
2752  MAT4X3PNT(pt, dsp->dsp_i.dsp_stom, tmp);
2753 
2754 
2755  /* Computing in world coordinates */
2756  VMOVE(tmp, A);
2757  MAT4X3PNT(A, dsp->dsp_i.dsp_stom, tmp);
2758 
2759  VMOVE(tmp, C);
2760  MAT4X3PNT(C, dsp->dsp_i.dsp_stom, tmp);
2761 
2762  VMOVE(tmp, D);
2763  MAT4X3PNT(D, dsp->dsp_i.dsp_stom, tmp);
2764 
2765  VMOVE(tmp, E);
2766  MAT4X3PNT(E, dsp->dsp_i.dsp_stom, tmp);
2767 
2768  VSUB2(Vac, C, A);
2769  VSUB2(Vde, E, D);
2770 
2771  VUNITIZE(Vac);
2772  VUNITIZE(Vde);
2773  VCROSS(N, Vac, Vde);
2774 
2775  if (RT_G_DEBUG & DEBUG_HF) {
2776  VPRINT("\tA", A);
2777  VPRINT("\tC", C);
2778  VPRINT("\tD", D);
2779  VPRINT("\tE", E);
2780  VPRINT("\tVac", Vac);
2781  VPRINT("\tVde", Vde);
2782  VPRINT("\tModel Cross N", N);
2783  }
2784  VUNITIZE(N);
2785 
2786  if (RT_G_DEBUG & DEBUG_HF) {
2787  VPRINT("\tModel Unit N", N);
2788  }
2789  if (fd) {
2790  VJOIN1(endpt, pt, len, N);
2791 
2792  pl_color(fd, 220, 220, 90);
2793  pdv_3line(fd, A, C);
2794  pdv_3line(fd, D, E);
2795 
2796  pdv_3line(fd, pt, endpt);
2797  }
2798 
2799 
2800 }
2801 
2802 
2803 /**
2804  * Given ONE ray distance, return the normal and entry/exit point.
2805  */
2806 void
2807 rt_dsp_norm(register struct hit *hitp, struct soltab *stp, register struct xray *rp)
2808 {
2809  vect_t N, t, tmp, A;
2810  struct dsp_specific *dsp = (struct dsp_specific *)stp->st_specific;
2811  vect_t Anorm, Bnorm, Dnorm, Cnorm, ABnorm, CDnorm;
2812  fastf_t Xfrac, Yfrac;
2813  int x, y;
2814  point_t pt;
2815  fastf_t dot;
2816  fastf_t len;
2817  FILE *fd = (FILE *)NULL;
2818 
2819 
2820  RT_DSP_CK_MAGIC(dsp);
2821  BU_CK_VLS(&dsp->dsp_i.dsp_name);
2822 
2823  switch (dsp->dsp_i.dsp_datasrc) {
2824  case RT_DSP_SRC_V4_FILE:
2825  case RT_DSP_SRC_FILE:
2826  BU_CK_MAPPED_FILE(dsp->dsp_i.dsp_mp);
2827  break;
2828  case RT_DSP_SRC_OBJ:
2829  RT_CK_DB_INTERNAL(dsp->dsp_i.dsp_bip);
2830  RT_CK_BINUNIF(dsp->dsp_i.dsp_bip->idb_ptr);
2831  break;
2832  }
2833 
2834  if (RT_G_DEBUG & DEBUG_HF) {
2835  bu_log("rt_dsp_norm(%g %g %g)\n", V3ARGS(hitp->hit_normal));
2836  VJOIN1(hitp->hit_point, rp->r_pt, hitp->hit_dist, rp->r_dir);
2837  VPRINT("\thit point", hitp->hit_point);
2838  bu_log("hit dist: %g\n", hitp->hit_dist);
2839  bu_log("%s:%d vpriv: %g, %g %g\n",
2840  __FILE__, __LINE__, V3ARGS(hitp->hit_vpriv));
2841  }
2842 
2843  if (hitp->hit_surfno < XMIN || hitp->hit_surfno > ZTOP) {
2844  bu_log("%s:%d bogus surface of DSP %d\n",
2845  __FILE__, __LINE__, hitp->hit_surfno);
2846  bu_bomb("");
2847  }
2848 
2849  /* compute hit point */
2850  VJOIN1(hitp->hit_point, rp->r_pt, hitp->hit_dist, rp->r_dir);
2851 
2852 
2853  if (hitp->hit_surfno != ZTOP || !dsp->dsp_i.dsp_smooth) {
2854  /* We've hit one of the sides or bottom, or the user didn't
2855  * ask for smoothing of the elevation data, so there's no
2856  * interpolation to do. Just transform the normal to model
2857  * space, and compute the actual hit point
2858  */
2859 
2860  /* transform normal into model space */
2861  MAT4X3VEC(tmp, dsp->dsp_i.dsp_mtos, hitp->hit_normal);
2862  VUNITIZE(tmp);
2863  VMOVE(hitp->hit_normal, tmp);
2864 
2865  if (RT_G_DEBUG & DEBUG_HF)
2866  bu_log("\tno Interpolation needed. Normal: %g, %g, %g\n",
2867  V3ARGS(hitp->hit_normal));
2868  return;
2869  }
2870 
2871  if (RT_G_DEBUG & DEBUG_HF)
2872  bu_log("\tNormal Interpolation flag: %d\n", dsp->dsp_i.dsp_smooth);
2873 
2874 
2875  /* compute the distance between grid points in model space */
2876  VSET(tmp, 1.0, 0.0, 0.0);
2877  MAT4X3VEC(t, dsp->dsp_i.dsp_stom, tmp);
2878  len = MAGNITUDE(t);
2879 
2880  if (RT_G_DEBUG & DEBUG_HF) {
2881  struct bu_vls str = BU_VLS_INIT_ZERO;
2882 
2883  bu_vls_printf(&str, "dsp_gourand%02d.plot3", plot_file_num++);
2884  bu_log("plotting normals in %s", bu_vls_addr(&str));
2885  fd = fopen(bu_vls_addr(&str), "w");
2886  bu_vls_free(&str);
2887 
2888  /* plot the ray */
2889  pl_color(fd, 255, 0, 0);
2890  pdv_3line(fd, rp->r_pt, hitp->hit_point);
2891 
2892  /* plot the normal we started with */
2893  pl_color(fd, 0, 255, 0);
2894  VJOIN1(tmp, hitp->hit_point, len, hitp->hit_normal);
2895  pdv_3line(fd, hitp->hit_point, tmp);
2896 
2897  }
2898 
2899  /* get the cell we hit */
2900  x = hitp->hit_vpriv[X];
2901  y = hitp->hit_vpriv[Y];
2902 
2903  compute_normal_at_gridpoint(Anorm, dsp, x, y, fd, len);
2904  compute_normal_at_gridpoint(Bnorm, dsp, x+1, y, fd, len);
2905  compute_normal_at_gridpoint(Dnorm, dsp, x+1, y+1, fd, len);
2906  compute_normal_at_gridpoint(Cnorm, dsp, x, y+1, fd, len);
2907 
2908  /* transform the hit point into DSP space for determining
2909  * interpolation
2910  */
2911  MAT4X3PNT(pt, dsp->dsp_i.dsp_mtos, hitp->hit_point);
2912 
2913  Xfrac = (pt[X] - x);
2914  Yfrac = (pt[Y] - y);
2915  if (RT_G_DEBUG & DEBUG_HF)
2916  bu_log("Xfract:%g Yfract:%g\n", Xfrac, Yfrac);
2917 
2918  if (Xfrac < 0.0) Xfrac = 0.0;
2919  else if (Xfrac > 1.0) Xfrac = 1.0;
2920 
2921  if (Yfrac < 0.0) Yfrac = 0.0;
2922  else if (Yfrac > 1.0) Yfrac = 1.0;
2923 
2924 
2925  if (dsp->dsp_i.dsp_smooth == 2) {
2926  /* This is an experiment to "flatten" the curvature of the dsp
2927  * near the grid points
2928  */
2929 #define SMOOTHSTEP(x) ((x)*(x)*(3 - 2*(x)))
2930  Xfrac = SMOOTHSTEP(Xfrac);
2931  Yfrac = SMOOTHSTEP(Yfrac);
2932 #undef SMOOTHSTEP
2933  }
2934 
2935  /* we compute the normal along the "X edges" of the cell */
2936  VSCALE(Anorm, Anorm, (1.0-Xfrac));
2937  VSCALE(Bnorm, Bnorm, Xfrac);
2938  VADD2(ABnorm, Anorm, Bnorm);
2939  VUNITIZE(ABnorm);
2940 
2941  VSCALE(Cnorm, Cnorm, (1.0-Xfrac));
2942  VSCALE(Dnorm, Dnorm, Xfrac);
2943  VADD2(CDnorm, Dnorm, Cnorm);
2944  VUNITIZE(CDnorm);
2945 
2946  /* now we interpolate the two X edge normals to get the final one
2947  */
2948  VSCALE(ABnorm, ABnorm, (1.0-Yfrac));
2949  VSCALE(CDnorm, CDnorm, Yfrac);
2950  VADD2(N, ABnorm, CDnorm);
2951 
2952  VUNITIZE(N);
2953 
2954  dot = VDOT(N, rp->r_dir);
2955  if (RT_G_DEBUG & DEBUG_HF)
2956  bu_log("interpolated %g %g %g dot:%g\n", V3ARGS(N), dot);
2957 
2958  if ((ZERO(hitp->hit_vpriv[Z]) && dot > 0.0)/* in-hit needs fix */ ||
2959  (ZERO(hitp->hit_vpriv[Z] - 1.0) && dot < 0.0)/* out-hit needs fix */) {
2960  /* bring the normal back to being perpendicular to the ray to
2961  * avoid "flipped normal" warnings
2962  */
2963  VCROSS(A, rp->r_dir, N);
2964  VCROSS(N, A, rp->r_dir);
2965 
2966  VUNITIZE(N);
2967 
2968  dot = VDOT(N, rp->r_dir);
2969 
2970 
2971  if (RT_G_DEBUG & DEBUG_HF)
2972  bu_log("corrected: %g %g %g dot:%g\n", V3ARGS(N), dot);
2973  }
2974  VMOVE(hitp->hit_normal, N);
2975 
2976  if (RT_G_DEBUG & DEBUG_HF) {
2977  pl_color(fd, 255, 255, 255);
2978  VJOIN1(tmp, hitp->hit_point, len, hitp->hit_normal);
2979  pdv_3line(fd, hitp->hit_point, tmp);
2980 
2982  fclose(fd);
2984  }
2985 }
2986 
2987 
2988 /**
2989  * Return the curvature of the dsp.
2990  */
2991 void
2992 rt_dsp_curve(register struct curvature *cvp, register struct hit *hitp, struct soltab *stp)
2993 {
2994  if (stp) RT_CK_SOLTAB(stp);
2995 
2996  if (RT_G_DEBUG & DEBUG_HF)
2997  bu_log("rt_dsp_curve()\n");
2998 
2999  cvp->crv_c1 = cvp->crv_c2 = 0;
3000 
3001  /* any tangent direction */
3002  bn_vec_ortho(cvp->crv_pdir, hitp->hit_normal);
3003 }
3004 
3005 
3006 /**
3007  * For a hit on the surface of a dsp, return the (u, v) coordinates
3008  * of the hit point, 0 <= u, v <= 1.
3009  * u = azimuth
3010  * v = elevation
3011  */
3012 void
3013 rt_dsp_uv(struct application *ap, struct soltab *stp, register struct hit *hitp, register struct uvcoord *uvp)
3014 {
3015  register struct dsp_specific *dsp =
3016  (struct dsp_specific *)stp->st_specific;
3017  point_t pt;
3018  vect_t tmp;
3019  fastf_t r;
3020  fastf_t min_r_U, min_r_V;
3021  vect_t norm;
3022  vect_t rev_dir;
3023  fastf_t dot_N;
3024  vect_t UV_dir;
3025  vect_t U_dir, V_dir;
3026  fastf_t U_len, V_len;
3027  fastf_t one_over_len;
3028 
3029  MAT4X3PNT(pt, dsp->dsp_i.dsp_mtos, hitp->hit_point);
3030 
3031  /* compute U, V */
3032  uvp->uv_u = pt[X] / (double)XSIZ(dsp);
3033  CLAMP(uvp->uv_u, 0.0, 1.0);
3034 
3035  uvp->uv_v = pt[Y] / (double)YSIZ(dsp);
3036  CLAMP(uvp->uv_v, 0.0, 1.0);
3037 
3038 
3039  /* du, dv indicate the extent of the ray radius in UV coordinates.
3040  * To compute this, transform unit vectors from solid space to
3041  * model space. We remember the length of the resultant vectors
3042  * and then unitize them to get u, v directions in model
3043  * coordinate space.
3044  */
3045  VSET(tmp, XSIZ(dsp), 0.0, 0.0); /* X direction vector */
3046  MAT4X3VEC(U_dir, dsp->dsp_i.dsp_stom, tmp); /* into model space */
3047  U_len = MAGNITUDE(U_dir); /* remember model space length */
3048  one_over_len = 1.0/U_len;
3049  VSCALE(U_dir, U_dir, one_over_len); /* scale to unit length in model space */
3050 
3051  VSET(tmp, 0.0, YSIZ(dsp), 0.0); /* Y direction vector */
3052  MAT4X3VEC(V_dir, dsp->dsp_i.dsp_stom, tmp);
3053  V_len = MAGNITUDE(V_dir);
3054  one_over_len = 1.0/V_len;
3055  VSCALE(V_dir, V_dir, one_over_len); /* scale to unit length in model space */
3056 
3057  /* divide the hit-point radius by the U/V unit length distance */
3058  r = ap->a_rbeam + ap->a_diverge * hitp->hit_dist;
3059  min_r_U = r / U_len;
3060  min_r_V = r / V_len;
3061 
3062  /* compute UV_dir, a vector in the plane of the hit point
3063  * (surface) which points in the anti-rayward direction
3064  */
3065  VREVERSE(rev_dir, ap->a_ray.r_dir);
3066  VMOVE(norm, hitp->hit_normal);
3067  dot_N = VDOT(rev_dir, norm);
3068  VJOIN1(UV_dir, rev_dir, -dot_N, norm);
3069  VUNITIZE(UV_dir);
3070 
3071  if (ZERO(dot_N)) {
3072  /* ray almost perfectly 90 degrees to surface */
3073  uvp->uv_du = min_r_U;
3074  uvp->uv_dv = min_r_V;
3075  } else {
3076  /* somehow this computes the extent of U and V in the radius */
3077  uvp->uv_du = (r / U_len) * VDOT(UV_dir, U_dir) / dot_N;
3078  uvp->uv_dv = (r / V_len) * VDOT(UV_dir, V_dir) / dot_N;
3079  }
3080 
3081 
3082  if (uvp->uv_du < 0.0)
3083  uvp->uv_du = -uvp->uv_du;
3084  if (uvp->uv_du < min_r_U)
3085  uvp->uv_du = min_r_U;
3086 
3087  if (uvp->uv_dv < 0.0)
3088  uvp->uv_dv = -uvp->uv_dv;
3089  if (uvp->uv_dv < min_r_V)
3090  uvp->uv_dv = min_r_V;
3091 
3092  if (RT_G_DEBUG & DEBUG_HF)
3093  bu_log("rt_dsp_uv(pt:%g, %g siz:%u, %u)\n U_len=%g V_len=%g\n r=%g rbeam=%g diverge=%g dist=%g\n u=%g v=%g du=%g dv=%g\n",
3094  pt[X], pt[Y], XSIZ(dsp), YSIZ(dsp),
3095  U_len, V_len,
3096  r, ap->a_rbeam, ap->a_diverge, hitp->hit_dist,
3097  uvp->uv_u, uvp->uv_v,
3098  uvp->uv_du, uvp->uv_dv);
3099 }
3100 
3101 
3102 void
3103 rt_dsp_free(register struct soltab *stp)
3104 {
3105  register struct dsp_specific *dsp =
3106  (struct dsp_specific *)stp->st_specific;
3107 
3108  if (RT_G_DEBUG & DEBUG_HF)
3109  bu_log("rt_dsp_free()\n");
3110 
3111  switch (dsp->dsp_i.dsp_datasrc) {
3112  case RT_DSP_SRC_V4_FILE:
3113  case RT_DSP_SRC_FILE:
3114  if (dsp->dsp_i.dsp_mp) {
3115  bu_close_mapped_file(dsp->dsp_i.dsp_mp);
3116  } else if (dsp->dsp_i.dsp_buf) {
3117  bu_free(dsp->dsp_i.dsp_buf, "dsp fake data");
3118  }
3119  break;
3120  case RT_DSP_SRC_OBJ:
3121  break;
3122  }
3123 
3124  BU_PUT(dsp, struct dsp_specific);
3125 }
3126 
3127 
3128 /* FIXME: can't 'rt_dsp_class' be replaced by 'rt_generic_class'? */
3129 int
3130 rt_dsp_class(const struct soltab *UNUSED(s), const vect_t UNUSED(v0), const vect_t UNUSED(v2), const struct bn_tol *UNUSED(b))
3131 {
3132  if (RT_G_DEBUG & DEBUG_HF)
3133  bu_log("rt_dsp_class()\n");
3134 
3135  return 0;
3136 }
3137 
3138 
3139 int
3140 rt_dsp_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))
3141 {
3142  struct rt_dsp_internal *dsp_ip =
3143  (struct rt_dsp_internal *)ip->idb_ptr;
3144  point_t m_pt;
3145  point_t s_pt;
3146  point_t o_pt;
3147  unsigned int x, y;
3148  int step;
3149  unsigned int xlim = dsp_ip->dsp_xcnt - 1;
3150  unsigned int ylim = dsp_ip->dsp_ycnt - 1;
3151  int xfudge, yfudge;
3152  int drawing;
3153 
3154  if (RT_G_DEBUG & DEBUG_HF)
3155  bu_log("rt_dsp_plot()\n");
3156 
3157  BU_CK_LIST_HEAD(vhead);
3158  RT_CK_DB_INTERNAL(ip);
3159  RT_DSP_CK_MAGIC(dsp_ip);
3160 
3161  switch (dsp_ip->dsp_datasrc) {
3162  case RT_DSP_SRC_V4_FILE:
3163  case RT_DSP_SRC_FILE:
3164  if (!dsp_ip->dsp_mp) {
3165  bu_log("WARNING: Cannot find data file for displacement map (DSP)\n");
3166  if (bu_vls_addr(&dsp_ip->dsp_name)) {
3167  bu_log(" DSP data file [%s] not found or empty\n", bu_vls_addr(&dsp_ip->dsp_name));
3168  } else {
3169  bu_log(" DSP data file not found or not specified\n");
3170  }
3171  return 0;
3172  }
3173  BU_CK_MAPPED_FILE(dsp_ip->dsp_mp);
3174  break;
3175  case RT_DSP_SRC_OBJ:
3176  if (!dsp_ip->dsp_bip) {
3177  bu_log("WARNING: Cannot find data object for displacement map (DSP)\n");
3178  if (bu_vls_addr(&dsp_ip->dsp_name)) {
3179  bu_log(" DSP data object [%s] not found or empty\n", bu_vls_addr(&dsp_ip->dsp_name));
3180  } else {
3181  bu_log(" DSP data object not found or not specified\n");
3182  }
3183  return 0;
3184  }
3185  RT_CK_DB_INTERNAL(dsp_ip->dsp_bip);
3186  RT_CK_BINUNIF(dsp_ip->dsp_bip->idb_ptr);
3187  break;
3188  }
3189 
3190 
3191 #define MOVE(_pt) \
3192  MAT4X3PNT(m_pt, dsp_ip->dsp_stom, _pt); \
3193  RT_ADD_VLIST(vhead, m_pt, BN_VLIST_LINE_MOVE)
3194 
3195 #define DRAW(_pt) \
3196  MAT4X3PNT(m_pt, dsp_ip->dsp_stom, _pt); \
3197  RT_ADD_VLIST(vhead, m_pt, BN_VLIST_LINE_DRAW)
3198 
3199 
3200  /* Draw the Bottom */
3201  VSETALL(s_pt, 0.0);
3202  MOVE(s_pt);
3203 
3204  s_pt[X] = xlim;
3205  DRAW(s_pt);
3206 
3207  s_pt[Y] = ylim;
3208  DRAW(s_pt);
3209 
3210  s_pt[X] = 0.0;
3211  DRAW(s_pt);
3212 
3213  s_pt[Y] = 0.0;
3214  DRAW(s_pt);
3215 
3216 
3217  /* Draw the corners */
3218  s_pt[Z] = DSP(dsp_ip, 0, 0);
3219  DRAW(s_pt);
3220 
3221  VSET(s_pt, xlim, 0.0, 0.0);
3222  MOVE(s_pt);
3223  s_pt[Z] = DSP(dsp_ip, xlim, 0);
3224  DRAW(s_pt);
3225 
3226 
3227  VSET(s_pt, xlim, ylim, 0.0);
3228  MOVE(s_pt);
3229  s_pt[Z] = DSP(dsp_ip, xlim, ylim);
3230  DRAW(s_pt);
3231 
3232  VSET(s_pt, 0.0, ylim, 0.0);
3233  MOVE(s_pt);
3234  s_pt[Z] = DSP(dsp_ip, 0, ylim);
3235  DRAW(s_pt);
3236 
3237 
3238  /* Draw the outside line of the top We draw the four sides of the
3239  * top at full resolution. This helps edge matching. The inside
3240  * of the top, we draw somewhat coarser
3241  */
3242  for (y = 0; y < dsp_ip->dsp_ycnt; y += ylim) {
3243  VSET(s_pt, 0.0, y, DSP(dsp_ip, 0, y));
3244  MOVE(s_pt);
3245 
3246  for (x = 0; x < dsp_ip->dsp_xcnt; x++) {
3247  s_pt[X] = x;
3248  s_pt[Z] = DSP(dsp_ip, x, y);
3249  DRAW(s_pt);
3250  }
3251  }
3252 
3253 
3254  for (x = 0; x < dsp_ip->dsp_xcnt; x += xlim) {
3255  VSET(s_pt, x, 0.0, DSP(dsp_ip, x, 0));
3256  MOVE(s_pt);
3257 
3258  for (y = 0; y < dsp_ip->dsp_ycnt; y++) {
3259  s_pt[Y] = y;
3260  s_pt[Z] = DSP(dsp_ip, x, y);
3261  DRAW(s_pt);
3262  }
3263  }
3264 
3265  /* now draw the body of the top */
3266  if (!ZERO(ttol->rel)) {
3267  unsigned int rstep;
3268  rstep = dsp_ip->dsp_xcnt;
3269  V_MAX(rstep, dsp_ip->dsp_ycnt);
3270  step = (int)(ttol->rel * rstep);
3271  } else {
3272  int goal = 10000;
3273  goal -= 5;
3274  goal -= 8 + 2 * (dsp_ip->dsp_xcnt+dsp_ip->dsp_ycnt);
3275 
3276  if (goal <= 0) return 0;
3277 
3278  /* Compute data stride based upon producing no more than 'goal' vectors */
3279  step = ceil(
3280  sqrt(2*(xlim)*(ylim) /
3281  (double)goal)
3282  );
3283  }
3284  if (step < 1) step = 1;
3285 
3286 
3287  xfudge = (dsp_ip->dsp_xcnt % step + step) / 2;
3288  yfudge = (dsp_ip->dsp_ycnt % step + step) / 2;
3289 
3290  if (xfudge < 1) xfudge = 1;
3291  if (yfudge < 1) yfudge = 1;
3292 
3293  /* draw the horizontal (y==const) lines */
3294  for (y = yfudge; y < ylim; y += step) {
3295  VSET(s_pt, 0.0, y, DSP(dsp_ip, 0, y));
3296  VMOVE(o_pt, s_pt);
3297  if (!ZERO(o_pt[Z])) {
3298  drawing = 1;
3299  MOVE(o_pt);
3300  } else {
3301  drawing = 0;
3302  }
3303 
3304  for (x = xfudge; x < xlim; x+=step) {
3305  s_pt[X] = x;
3306 
3307  s_pt[Z] = DSP(dsp_ip, x, y);
3308  if (!ZERO(s_pt[Z])) {
3309  if (drawing) {
3310  DRAW(s_pt);
3311  } else {
3312  MOVE(o_pt);
3313  DRAW(s_pt);
3314  drawing = 1;
3315  }
3316  } else {
3317  if (drawing) {
3318  DRAW(s_pt);
3319  drawing = 0;
3320  }
3321  }
3322 
3323  VMOVE(o_pt, s_pt);
3324  }
3325 
3326  s_pt[X] = xlim;
3327  s_pt[Z] = DSP(dsp_ip, xlim, y);
3328  if (!ZERO(s_pt[Z])) {
3329  if (drawing) {
3330  DRAW(s_pt);
3331  } else {
3332  MOVE(o_pt);
3333  DRAW(s_pt);
3334  drawing = 1;
3335  }
3336  } else {
3337  if (drawing) {
3338  DRAW(s_pt);
3339  drawing = 0;
3340  }
3341  }
3342 
3343  }
3344 
3345  for (x = xfudge; x < xlim; x += step) {
3346  VSET(s_pt, x, 0.0, DSP(dsp_ip, x, 0));
3347  VMOVE(o_pt, s_pt);
3348  if (!ZERO(o_pt[Z])) {
3349  drawing = 1;
3350  MOVE(o_pt);
3351  } else {
3352  drawing = 0;
3353  }
3354 
3355 
3356  for (y = yfudge; y < ylim; y+=step) {
3357  s_pt[Y] = y;
3358 
3359  s_pt[Z] = DSP(dsp_ip, x, y);
3360  if (!ZERO(s_pt[Z])) {
3361  if (drawing) {
3362  DRAW(s_pt);
3363  } else {
3364  MOVE(o_pt);
3365  DRAW(s_pt);
3366  drawing = 1;
3367  }
3368  } else {
3369  if (drawing) {
3370  DRAW(s_pt);
3371  drawing = 0;
3372  }
3373  }
3374 
3375  VMOVE(o_pt, s_pt);
3376  }
3377 
3378  s_pt[Y] = ylim;
3379  s_pt[Z] = DSP(dsp_ip, x, ylim);
3380  if (!ZERO(s_pt[Z])) {
3381  if (drawing) {
3382  DRAW(s_pt);
3383  } else {
3384  MOVE(o_pt);
3385  DRAW(s_pt);
3386  drawing = 1;
3387  }
3388  } else {
3389  if (drawing) {
3390  DRAW(s_pt);
3391  drawing = 0;
3392  }
3393  }
3394  }
3395 
3396 #undef MOVE
3397 #undef DRAW
3398  return 0;
3399 }
3400 
3401 
3402 /*
3403  * Determine the cut direction for a DSP cell. This routine is used by
3404  * rt_dsp_tess(). It somewhat duplicates code from permute_cell(),
3405  * which is a bad thing, but permute_cell() had other side effects not
3406  * desired here.
3407  *
3408  * inputs:
3409  * dsp_ip - pointer to the rt_dsp_internal struct fro thi DSP
3410  * x - the DSP cell x-coord of the lower left corner of the cell
3411  * y - the DSP cell y-coord of the lower left corner of the cell
3412  * xlim - the maximum value of the DSP x coordinates
3413  * ylim - the maximum value of the DSP y coordinates
3414  * return:
3415  * the direction for cutting this cell.
3416  * DSP_CUT_DIR_llUR - lower left to upper right
3417  * DSP_CUT_DIR_ULlr - upper right to lower left
3418  */
3419 HIDDEN int
3420 get_cut_dir(struct rt_dsp_internal *dsp_ip, int x, int y, int xlim, int ylim)
3421 {
3422 /*
3423  * height array contains DSP values:
3424  * height[0] is at (x, y)
3425  * height[1] is at (x+1, y)
3426  * height[2] is at (x+1, y+1)
3427  * height[3] is at (x, y+1)
3428  * height[4] is at (x-1, y-1)
3429  * height[5] is at (x+2, y-1)
3430  * height[6] is at (x+2, y+2)
3431  * height[7] is at (x-1, y+2)
3432  *
3433  * 7 6
3434  * \ /
3435  * \ /
3436  * 3---2
3437  * | |
3438  * | |
3439  * 0---1
3440  * / \
3441  * / \
3442  * 4 5
3443  *
3444  * (0, 1, 2, 3) is the cell of interest.
3445  */
3446  int height[8];
3447  int xx, yy;
3448  fastf_t c02, c13; /* curvature in direction 0<->2, and 1<->3 */
3449 
3450  if (dsp_ip->dsp_cuttype != DSP_CUT_DIR_ADAPT) {
3451  /* not using adaptive cut type, so just return the cut type */
3452  return dsp_ip->dsp_cuttype;
3453  }
3454 
3455  /* fill in the height array */
3456  xx = x;
3457  yy = y;
3458  height[0] = DSP(dsp_ip, xx, yy);
3459  xx = x+1;
3460  if (xx > xlim) xx = xlim;
3461  height[1] = DSP(dsp_ip, xx, yy);
3462  yy = y+1;
3463  if (yy > ylim) yy = ylim;
3464  height[2] = DSP(dsp_ip, xx, yy);
3465  xx = x;
3466  height[3] = DSP(dsp_ip, xx, yy);
3467  xx = x-1;
3468  if (xx < 0) xx = 0;
3469  yy = y-1;
3470  if (yy < 0) yy = 0;
3471  height[4] = DSP(dsp_ip, xx, yy);
3472  xx = x+2;
3473  if (xx > xlim) xx = xlim;
3474  height[5] = DSP(dsp_ip, xx, yy);
3475  yy = y+2;
3476  if (yy > ylim) yy = ylim;
3477  height[6] = DSP(dsp_ip, xx, yy);
3478  xx = x-1;
3479  if (xx < 0) xx = 0;
3480  height[7] = DSP(dsp_ip, xx, yy);
3481 
3482  /* compute curvature along the 0<->2 direction */
3483  c02 = fabs(height[2] + height[4] - 2*height[0]) + fabs(height[6] + height[0] - 2*height[2]);
3484 
3485 
3486  /* compute curvature along the 1<->3 direction */
3487  c13 = fabs(height[3] + height[5] - 2*height[1]) + fabs(height[7] + height[1] - 2*height[3]);
3488 
3489  if (c02 < c13) {
3490  /* choose the 0<->2 direction */
3491  return DSP_CUT_DIR_llUR;
3492  } else {
3493  /* choose the other direction */
3494  return DSP_CUT_DIR_ULlr;
3495  }
3496 }
3497 
3498 
3499 /* a macro to check if these vertices can produce a valid face */
3500 #define CHECK_VERTS(_v) \
3501  ((*_v[0] != *_v[1]) || (*_v[0] == NULL)) && \
3502  ((*_v[0] != *_v[2]) || (*_v[2] == NULL)) && \
3503  ((*_v[1] != *_v[2]) || (*_v[1] == NULL))
3504 
3505 /**
3506  * Returns -
3507  * -1 failure
3508  * 0 OK. *r points to nmgregion that holds this tessellation.
3509  */
3510 int
3511 rt_dsp_tess(struct nmgregion **r, struct model *m, struct rt_db_internal *ip, const struct rt_tess_tol *UNUSED(ttol), const struct bn_tol *tol)
3512 {
3513  struct rt_dsp_internal *dsp_ip;
3514  struct shell *s;
3515  int xlim;
3516  int ylim;
3517  int x, y;
3518  point_t pt[4];
3519  point_t tmp_pt;
3520  int base_vert_count;
3521  struct vertex **base_verts;
3522  struct vertex **verts[3];
3523  struct vertex *hole_verts[3];
3524  struct faceuse *fu;
3525  struct faceuse *base_fu;
3526  struct vertex **strip1Verts;
3527  struct vertex **strip2Verts;
3528  int base_vert_no1;
3529  int base_vert_no2;
3530  int has_holes = 0;
3531 
3532  if (RT_G_DEBUG & DEBUG_HF)
3533  bu_log("rt_dsp_tess()\n");
3534 
3535  /* do a bunch of checks to make sure all is well */
3536 
3537  RT_CK_DB_INTERNAL(ip);
3538  dsp_ip = (struct rt_dsp_internal *)ip->idb_ptr;
3539  RT_DSP_CK_MAGIC(dsp_ip);
3540 
3541  switch (dsp_ip->dsp_datasrc) {
3542  case RT_DSP_SRC_FILE:
3543  case RT_DSP_SRC_V4_FILE:
3544  if (!dsp_ip->dsp_mp) {
3545  bu_log("WARNING: Cannot find data file for displacement map (DSP)\n");
3546  if (bu_vls_addr(&dsp_ip->dsp_name)) {
3547  bu_log(" DSP data file [%s] not found or empty\n", bu_vls_addr(&dsp_ip->dsp_name));
3548  } else {
3549  bu_log(" DSP data file not found or not specified\n");
3550  }
3551  return -1;
3552  }
3553  BU_CK_MAPPED_FILE(dsp_ip->dsp_mp);
3554  break;
3555  case RT_DSP_SRC_OBJ:
3556  if (!dsp_ip->dsp_bip) {
3557  bu_log("WARNING: Cannot find data object for displacement map (DSP)\n");
3558  if (bu_vls_addr(&dsp_ip->dsp_name)) {
3559  bu_log(" DSP data object [%s] not found or empty\n", bu_vls_addr(&dsp_ip->dsp_name));
3560  } else {
3561  bu_log(" DSP data object not found or not specified\n");
3562  }
3563  return -1;
3564  }
3565  RT_CK_DB_INTERNAL(dsp_ip->dsp_bip);
3566  RT_CK_BINUNIF(dsp_ip->dsp_bip->idb_ptr);
3567  break;
3568  }
3569 
3570 
3571  xlim = dsp_ip->dsp_xcnt - 1;
3572  ylim = dsp_ip->dsp_ycnt - 1;
3573 
3574  /* Base_verts will contain the vertices for the base face ordered
3575  * correctly to create the base face.
3576  * Cannot simply use the four corners, because that would not
3577  * create a valid NMG.
3578  * base_verts[0] is at (0, 0)
3579  * base_verts[ylim] is at (0, ylim)
3580  * base_verts[ylim+xlim] is at (xlim, ylim)
3581  * base_verts[2*ylim+xlim] is at (xlim, 0)
3582  * base_verts[2*ylim+2*xlim-x] is at (x, 0)
3583  *
3584  * strip1Verts and strip2Verts are temporary storage for vertices
3585  * along the top of the dsp. For each strip of triangles at a
3586  * given x value, strip1Verts[y] is the vertex at (x, y, h) and
3587  * strip2Verts[y] is the vertex at (x+1, y, h), where h is the DSP
3588  * value at that point. After each strip of faces is created,
3589  * strip2Verts is copied to strip1Verts and strip2Verts is set to
3590  * all NULLs.
3591  */
3592 
3593  /* malloc space for the vertices */
3594  base_vert_count = 2*xlim + 2*ylim;
3595  base_verts = (struct vertex **)bu_calloc(base_vert_count, sizeof(struct vertex *), "base verts");
3596  strip1Verts = (struct vertex **)bu_calloc(ylim+1, sizeof(struct vertex *), "strip1Verts");
3597  strip2Verts = (struct vertex **)bu_calloc(ylim+1, sizeof(struct vertex *), "strip2Verts");
3598 
3599  /* Make region, empty shell, vertex */
3600  *r = nmg_mrsv(m);
3601  s = BU_LIST_FIRST(shell, &(*r)->s_hd);
3602 
3603  /* make the base face */
3604  base_fu = nmg_cface(s, base_verts, base_vert_count);
3605 
3606  /* assign geometry to the base_verts */
3607  /* start with x=0 edge */
3608  for (y = 0; y <= ylim; y++) {
3609  VSET(tmp_pt, 0, y, 0);
3610  MAT4X3PNT(pt[0], dsp_ip->dsp_stom, tmp_pt);
3611  nmg_vertex_gv(base_verts[y], pt[0]);
3612  }
3613  /* now do y=ylim edge */
3614  for (x = 1; x <= xlim; x++) {
3615  VSET(tmp_pt, x, ylim, 0);
3616  MAT4X3PNT(pt[0], dsp_ip->dsp_stom, tmp_pt);
3617  nmg_vertex_gv(base_verts[ylim+x], pt[0]);
3618  }
3619  /* now do x=xlim edge */
3620  for (y = 0; y < ylim; y++) {
3621  VSET(tmp_pt, xlim, y, 0);
3622  MAT4X3PNT(pt[0], dsp_ip->dsp_stom, tmp_pt);
3623  nmg_vertex_gv(base_verts[2*ylim+xlim-y], pt[0]);
3624  }
3625  /* now do y=0 edge */
3626  for (x = 1; x < xlim; x++) {
3627  VSET(tmp_pt, x, 0, 0);
3628  MAT4X3PNT(pt[0], dsp_ip->dsp_stom, tmp_pt);
3629  nmg_vertex_gv(base_verts[2*(xlim+ylim)-x], pt[0]);
3630  }
3631  if (nmg_fu_planeeqn(base_fu, tol) < 0) {
3632  bu_log("Failed to make base face\n");
3633  bu_free(base_verts, "base verts");
3634  bu_free(strip1Verts, "strip 1 verts");
3635  bu_free(strip2Verts, "strip 2 verts");
3636  return -1; /* FAIL */
3637  }
3638 
3639  /* if a displacement on the edge (x=0) is zero then strip1Verts at
3640  * that point is the corresponding base_vert
3641  */
3642  for (y = 0; y <= ylim; y++) {
3643  if (DSP(dsp_ip, 0, y) == 0) {
3644  strip1Verts[y] = base_verts[y];
3645  }
3646  }
3647 
3648  /* make faces along x=0 plane */
3649  for (y= 1; y <= ylim; y++) {
3650  verts[0] = &base_verts[y-1];
3651  verts[1] = &strip1Verts[y-1];
3652  verts[2] = &strip1Verts[y];
3653  if (CHECK_VERTS(verts)) {
3654  fu = nmg_cmface(s, verts, 3);
3655  if (y == 1 && strip1Verts[0]->vg_p == NULL) {
3656  VSET(tmp_pt, 0, 0, DSP(dsp_ip, 0, 0));
3657  MAT4X3PNT(pt[0], dsp_ip->dsp_stom, tmp_pt);
3658  nmg_vertex_gv(strip1Verts[0], pt[0]);
3659  }
3660  if (strip1Verts[y]->vg_p == NULL) {
3661  VSET(tmp_pt, 0, y, DSP(dsp_ip, 0, y));
3662  MAT4X3PNT(pt[0], dsp_ip->dsp_stom, tmp_pt);
3663  nmg_vertex_gv(strip1Verts[y], pt[0]);
3664  }
3665  if (nmg_fu_planeeqn(fu, tol) < 0) {
3666  bu_log("Failed to make x=0 face at y=%d\n", y);
3667  bu_free(base_verts, "base verts");
3668  bu_free(strip1Verts, "strip 1 verts");
3669  bu_free(strip2Verts, "strip 2 verts");
3670  return -1; /* FAIL */
3671  }
3672  }
3673  verts[0] = &base_verts[y-1];
3674  verts[1] = &strip1Verts[y];
3675  verts[2] = &base_verts[y];
3676  if (CHECK_VERTS(verts)) {
3677  fu = nmg_cmface(s, verts, 3);
3678  if (strip1Verts[y]->vg_p == NULL) {
3679  VSET(tmp_pt, 0, y, DSP(dsp_ip, 0, y));
3680  MAT4X3PNT(pt[0], dsp_ip->dsp_stom, tmp_pt);
3681  nmg_vertex_gv(strip1Verts[y], pt[0]);
3682  }
3683  if (nmg_fu_planeeqn(fu, tol) < 0) {
3684  bu_log("Failed to make x=0 face at y=%d\n", y);
3685  bu_free(base_verts, "base verts");
3686  bu_free(strip1Verts, "strip 1 verts");
3687  bu_free(strip2Verts, "strip 2 verts");
3688  return -1; /* FAIL */
3689  }
3690  }
3691  }
3692 
3693  /* make each strip of triangles. Make two triangles for each cell
3694  * (x, y)<->(x+1, y+1). Also make the vertical faces at y=0 and
3695  * y=ylim.
3696  */
3697  for (x = 0; x < xlim; x++) {
3698  /* set strip2Verts to base_verts values where the strip2Vert
3699  * is above a base_vert and DSP value is 0
3700  */
3701  if ((x+1) == xlim) {
3702  for (y = 0; y <= ylim; y++) {
3703  if (DSP(dsp_ip, xlim, y) == 0) {
3704  strip2Verts[y] = base_verts[2*ylim + xlim - y];
3705  }
3706  }
3707  } else {
3708  if (DSP(dsp_ip, x+1, 0) == 0) {
3709  strip2Verts[0] = base_verts[2*(ylim+xlim)-(x+1)];
3710  }
3711  if (DSP(dsp_ip, x+1, ylim) == 0) {
3712  strip2Verts[ylim] = base_verts[ylim+x+1];
3713  }
3714  }
3715 
3716  /* make the faces at y=0 for this strip */
3717  if (x == 0) {
3718  base_vert_no1 = 0;
3719  base_vert_no2 = 2*(ylim+xlim)-1;
3720  } else {
3721  base_vert_no1 = 2*(ylim + xlim) - x;
3722  base_vert_no2 = base_vert_no1 - 1;
3723  }
3724 
3725  verts[0] = &base_verts[base_vert_no1];
3726  verts[1] = &strip2Verts[0];
3727  verts[2] = &strip1Verts[0];
3728  if (CHECK_VERTS(verts)) {
3729  fu = nmg_cmface(s, verts, 3);
3730  VSET(tmp_pt, x+1, 0, DSP(dsp_ip, x+1, 0));
3731  MAT4X3PNT(pt[0], dsp_ip->dsp_stom, tmp_pt);
3732  nmg_vertex_gv(strip2Verts[0], pt[0]);
3733  if (nmg_fu_planeeqn(fu, tol) < 0) {
3734  bu_log("Failed to make first face at x=%d, y=%d\n", x, 0);
3735  bu_free(base_verts, "base verts");
3736  bu_free(strip1Verts, "strip 1 verts");
3737  bu_free(strip2Verts, "strip 2 verts");
3738  return -1; /* FAIL */
3739  }
3740  }
3741 
3742  verts[0] = &base_verts[base_vert_no1];
3743  verts[1] = &base_verts[base_vert_no2];
3744  verts[2] = &strip2Verts[0];
3745  if (CHECK_VERTS(verts)) {
3746  fu = nmg_cmface(s, verts, 3);
3747  if (nmg_fu_planeeqn(fu, tol) < 0) {
3748  bu_log("Failed to make second face at x=%d, y=%d\n", x, 0);
3749  bu_free(base_verts, "base verts");
3750  bu_free(strip1Verts, "strip 1 verts");
3751  bu_free(strip2Verts, "strip 2 verts");
3752  return -1; /* FAIL */
3753  }
3754  if (strip2Verts[0]->vg_p == NULL) {
3755  VSET(tmp_pt, x+1, 0, DSP(dsp_ip, x+1, 0));
3756  MAT4X3PNT(pt[0], dsp_ip->dsp_stom, tmp_pt);
3757  nmg_vertex_gv(strip2Verts[0], pt[0]);
3758  }
3759  }
3760 
3761  /* make the top faces for this strip */
3762  for (y = 0; y < ylim; y++) {
3763  int dir;
3764  /* get the cut direction for this cell */
3765  dir = get_cut_dir(dsp_ip, x, y, xlim, ylim);
3766 
3767  if (dir == DSP_CUT_DIR_llUR) {
3768  verts[0] = &strip1Verts[y];
3769  verts[1] = &strip2Verts[y];
3770  verts[2] = &strip2Verts[y+1];
3771  if (CHECK_VERTS(verts)) {
3772  if (DSP(dsp_ip, x, y) == 0 && DSP(dsp_ip, x+1, y) == 0 && DSP(dsp_ip, x+1, y+1) == 0) {
3773  /* make a hole, instead of a face */
3774  hole_verts[0] = strip1Verts[y];
3775  hole_verts[1] = strip2Verts[y];
3776  hole_verts[2] = strip2Verts[y+1];
3777  nmg_add_loop_to_face(s, base_fu, hole_verts, 3, OT_OPPOSITE);
3778  has_holes = 1;
3779  VSET(tmp_pt, x+1, y+1, DSP(dsp_ip, x+1, y+1));
3780  MAT4X3PNT(pt[0], dsp_ip->dsp_stom, tmp_pt);
3781  nmg_vertex_gv(hole_verts[2], pt[0]);
3782  strip2Verts[y+1] = hole_verts[2];
3783  } else {
3784  fu = nmg_cmface(s, verts, 3);
3785  VSET(tmp_pt, x+1, y+1, DSP(dsp_ip, x+1, y+1));
3786  MAT4X3PNT(pt[0], dsp_ip->dsp_stom, tmp_pt);
3787  nmg_vertex_gv(strip2Verts[y+1], pt[0]);
3788  if (nmg_fu_planeeqn(fu, tol) < 0) {
3789  bu_log("Failed to make first top face at x=%d, y=%d\n", x, y);
3790  bu_free(base_verts, "base verts");
3791  bu_free(strip1Verts, "strip 1 verts");
3792  bu_free(strip2Verts, "strip 2 verts");
3793  return -1; /* FAIL */
3794  }
3795  }
3796  }
3797  } else {
3798  verts[0] = &strip1Verts[y+1];
3799  verts[1] = &strip1Verts[y];
3800  verts[2] = &strip2Verts[y];
3801  if (CHECK_VERTS(verts)) {
3802  if (DSP(dsp_ip, x, y+1) == 0 && DSP(dsp_ip, x, y) == 0 && DSP(dsp_ip, x+1, y) == 0) {
3803  /* make a hole, instead of a face */
3804  hole_verts[0] = strip1Verts[y+1];
3805  hole_verts[1] = strip1Verts[y];
3806  hole_verts[2] = strip2Verts[y];
3807  nmg_add_loop_to_face(s, base_fu, hole_verts, 3, OT_OPPOSITE);
3808  has_holes = 1;
3809  } else {
3810  fu = nmg_cmface(s, verts, 3);
3811  if (nmg_fu_planeeqn(fu, tol) < 0) {
3812  bu_log("Failed to make first top face at x=%d, y=%d\n", x, y);
3813  bu_free(base_verts, "base verts");
3814  bu_free(strip1Verts, "strip 1 verts");
3815  bu_free(strip2Verts, "strip 2 verts");
3816  return -1; /* FAIL */
3817  }
3818  }
3819  }
3820  }
3821 
3822 
3823  if (dir == DSP_CUT_DIR_llUR) {
3824  verts[0] = &strip1Verts[y];
3825  verts[1] = &strip2Verts[y+1];
3826  verts[2] = &strip1Verts[y+1];
3827  if (CHECK_VERTS(verts)) {
3828  if (DSP(dsp_ip, x, y) == 0 && DSP(dsp_ip, x+1, y+1) == 0 && DSP(dsp_ip, x, y+1) == 0) {
3829  /* make a hole, instead of a face */
3830  hole_verts[0] = strip1Verts[y];
3831  hole_verts[1] = strip2Verts[y+1];
3832  hole_verts[2] = strip1Verts[y+1];
3833  nmg_add_loop_to_face(s, base_fu, hole_verts, 3, OT_OPPOSITE);
3834  has_holes = 1;
3835  VSET(tmp_pt, x+1, y+1, DSP(dsp_ip, x+1, y+1));
3836  MAT4X3PNT(pt[0], dsp_ip->dsp_stom, tmp_pt);
3837  nmg_vertex_gv(hole_verts[1], pt[0]);
3838  strip2Verts[y+1] = hole_verts[1];
3839  } else {
3840  fu = nmg_cmface(s, verts, 3);
3841  VSET(tmp_pt, x+1, y+1, DSP(dsp_ip, x+1, y+1));
3842  MAT4X3PNT(pt[0], dsp_ip->dsp_stom, tmp_pt);
3843  nmg_vertex_gv(strip2Verts[y+1], pt[0]);
3844  if (nmg_fu_planeeqn(fu, tol) < 0) {
3845  bu_log("Failed to make second top face at x=%d, y=%d\n", x, y);
3846  bu_free(base_verts, "base verts");
3847  bu_free(strip1Verts, "strip 1 verts");
3848  bu_free(strip2Verts, "strip 2 verts");
3849  return -1; /* FAIL */
3850  }
3851  }
3852  }
3853  } else {
3854  verts[0] = &strip2Verts[y];
3855  verts[1] = &strip2Verts[y+1];
3856  verts[2] = &strip1Verts[y+1];
3857  if (CHECK_VERTS(verts)) {
3858  if (DSP(dsp_ip, x+1, y) == 0 && DSP(dsp_ip, x+1, y+1) == 0 && DSP(dsp_ip, x, y+1) == 0) {
3859  /* make a hole, instead of a face */
3860  hole_verts[0] = strip2Verts[y];
3861  hole_verts[1] = strip2Verts[y+1];
3862  hole_verts[2] = strip1Verts[y+1];
3863  nmg_add_loop_to_face(s, base_fu, hole_verts, 3, OT_OPPOSITE);
3864  has_holes = 1;
3865  VSET(tmp_pt, x+1, y+1, DSP(dsp_ip, x+1, y+1));
3866  MAT4X3PNT(pt[0], dsp_ip->dsp_stom, tmp_pt);
3867  nmg_vertex_gv(hole_verts[1], pt[0]);
3868  strip2Verts[y+1] = hole_verts[1];
3869  } else {
3870  fu = nmg_cmface(s, verts, 3);
3871  VSET(tmp_pt, x+1, y+1, DSP(dsp_ip, x+1, y+1));
3872  MAT4X3PNT(pt[0], dsp_ip->dsp_stom, tmp_pt);
3873  nmg_vertex_gv(strip2Verts[y+1], pt[0]);
3874  if (nmg_fu_planeeqn(fu, tol) < 0) {
3875  bu_log("Failed to make second top face at x=%d, y=%d\n", x, y);
3876  bu_free(base_verts, "base verts");
3877  bu_free(strip1Verts, "strip 1 verts");
3878  bu_free(strip2Verts, "strip 2 verts");
3879  return -1; /* FAIL */
3880  }
3881  }
3882  }
3883  }
3884  }
3885 
3886  /* make the faces at the y=ylim plane for this strip */
3887  verts[0] = &strip1Verts[ylim];
3888  verts[1] = &strip2Verts[ylim];
3889  verts[2] = &base_verts[ylim+x+1];
3890  if (CHECK_VERTS(verts)) {
3891  fu = nmg_cmface(s, verts, 3);
3892  if (nmg_fu_planeeqn(fu, tol) < 0) {
3893  bu_log("Failed to make first face at x=%d, y=ylim\n", x);
3894  bu_free(base_verts, "base verts");
3895  bu_free(strip1Verts, "strip 1 verts");
3896  bu_free(strip2Verts, "strip 2 verts");
3897  return -1; /* FAIL */
3898  }
3899  }
3900 
3901  verts[0] = &base_verts[ylim+x+1];
3902  verts[1] = &base_verts[ylim+x];
3903  verts[2] = &strip1Verts[ylim];
3904  if (CHECK_VERTS(verts)) {
3905  fu = nmg_cmface(s, verts, 3);
3906  if (nmg_fu_planeeqn(fu, tol) < 0) {
3907  bu_log("Failed to make second face at x=%d, y=ylim\n", x);
3908  bu_free(base_verts, "base verts");
3909  bu_free(strip1Verts, "strip 1 verts");
3910  bu_free(strip2Verts, "strip 2 verts");
3911  return -1; /* FAIL */
3912  }
3913  }
3914 
3915  /* copy strip2 to strip1, set strip2 to all NULLs */
3916  for (y = 0; y <= ylim; y++) {
3917  strip1Verts[y] = strip2Verts[y];
3918  strip2Verts[y] = (struct vertex *)NULL;
3919  }
3920  }
3921 
3922  /* make faces at x=xlim plane */
3923  for (y = 0; y < ylim; y++) {
3924  base_vert_no1 = 2*ylim+xlim-y;
3925  verts[0] = &base_verts[base_vert_no1];
3926  verts[1] = &base_verts[base_vert_no1-1];
3927  verts[2] = &strip1Verts[y];
3928  if (CHECK_VERTS(verts)) {
3929  fu = nmg_cmface(s, verts, 3);
3930  if (nmg_fu_planeeqn(fu, tol) < 0) {
3931  bu_log("Failed to make first face at x=xlim, y=%d\n", y);
3932  bu_free(base_verts, "base verts");
3933  bu_free(strip1Verts, "strip 1 verts");
3934  bu_free(strip2Verts, "strip 2 verts");
3935  return -1; /* FAIL */
3936  }
3937  }
3938 
3939  verts[0] = &strip1Verts[y];
3940  verts[1] = &base_verts[base_vert_no1-1];
3941  verts[2] = &strip1Verts[y+1];
3942  if (CHECK_VERTS(verts)) {
3943  fu = nmg_cmface(s, verts, 3);
3944  if (nmg_fu_planeeqn(fu, tol) < 0) {
3945  bu_log("Failed to make second face at x=xlim, y=%d\n", y);
3946  bu_free(base_verts, "base verts");
3947  bu_free(strip1Verts, "strip 1 verts");
3948  bu_free(strip2Verts, "strip 2 verts");
3949  return -1; /* FAIL */
3950  }
3951  }
3952  }
3953 
3954 
3955  bu_free(base_verts, "base verts");
3956  bu_free(strip1Verts, "strip 1 verts");
3957  bu_free(strip2Verts, "strip 2 verts");
3958 
3959  if (has_holes) {
3960  /* do a bunch of joining and splitting of touching loops to
3961  * get the base face in a state that the triangulator can
3962  * handle
3963  */
3964  struct loopuse *lu;
3965  for (BU_LIST_FOR(lu, loopuse, &base_fu->lu_hd)) {
3967  }
3968  for (BU_LIST_FOR(lu, loopuse, &base_fu->lu_hd)) {
3969  if (lu->orientation != OT_UNSPEC) continue;
3970  nmg_lu_reorient(lu);
3971  }
3972  if (!nmg_kill_cracks(s))
3973  return -1;
3974 
3975  for (BU_LIST_FOR(lu, loopuse, &base_fu->lu_hd)) {
3976  nmg_split_touchingloops(lu, tol);
3977  }
3978  for (BU_LIST_FOR(lu, loopuse, &base_fu->lu_hd)) {
3979  if (lu->orientation != OT_UNSPEC) continue;
3980  nmg_lu_reorient(lu);
3981  }
3982  if (!nmg_kill_cracks(s))
3983  return -1;
3984 
3985  for (BU_LIST_FOR(lu, loopuse, &base_fu->lu_hd)) {
3987  }
3988  for (BU_LIST_FOR(lu, loopuse, &base_fu->lu_hd)) {
3989  if (lu->orientation != OT_UNSPEC) continue;
3990  nmg_lu_reorient(lu);
3991  }
3992  if (!nmg_kill_cracks(s))
3993  return -1;
3994  }
3995 
3996  /* Mark edges as real */
3997  (void)nmg_mark_edges_real(&s->l.magic);
3998 
3999  /* Compute "geometry" for region and shell */
4000  nmg_region_a(*r, tol);
4001 
4002  /* sanity check */
4003  nmg_make_faces_within_tol(s, tol);
4004 
4005  return 0;
4006 }
4007 
4008 
4009 /**
4010  * Retrieve data for DSP from external file.
4011  * Returns:
4012  * 0 Success
4013  * !0 Failure
4014  */
4015 HIDDEN int
4016 get_file_data(struct rt_dsp_internal *dsp_ip, const struct db_i *dbip)
4017 {
4018  struct bu_mapped_file *mf;
4019  int count, in_cookie, out_cookie;
4020 
4021  RT_DSP_CK_MAGIC(dsp_ip);
4022  RT_CK_DBI(dbip);
4023 
4024  /* get file */
4025  mf = dsp_ip->dsp_mp =
4027  bu_vls_addr(&dsp_ip->dsp_name), "dsp");
4028  if (!mf) {
4029  bu_log("mapped file open failure: %s%c%s\n",
4030  *dbip->dbi_filepath,BU_DIR_SEPARATOR,bu_vls_addr(&dsp_ip->dsp_name));
4031  return 0;
4032  }
4033 
4034  if ((size_t)dsp_ip->dsp_mp->buflen != (size_t)(dsp_ip->dsp_xcnt*dsp_ip->dsp_ycnt*2)) {
4035  bu_log("DSP buffer wrong size: %lu s/b %u ",
4036  (long unsigned int)dsp_ip->dsp_mp->buflen, dsp_ip->dsp_xcnt*dsp_ip->dsp_ycnt*2);
4037  return -1;
4038  }
4039 
4040  in_cookie = bu_cv_cookie("nus"); /* data is network unsigned short */
4041  out_cookie = bu_cv_cookie("hus");
4042 
4043  if (bu_cv_optimize(in_cookie) != bu_cv_optimize(out_cookie)) {
4044  size_t got;
4045  /* if we're on a little-endian machine we convert the input
4046  * file from network to host format
4047  */
4048  count = dsp_ip->dsp_xcnt * dsp_ip->dsp_ycnt;
4049  mf->apbuflen = count * sizeof(unsigned short);
4050  mf->apbuf = bu_malloc(mf->apbuflen, "apbuf");
4051 
4052  got = bu_cv_w_cookie(mf->apbuf, out_cookie, mf->apbuflen,
4053  mf->buf, in_cookie, count);
4054  if (got != (size_t)count) {
4055  bu_log("got %zu != count %d", got, count);
4056  bu_bomb("\n");
4057  }
4058  dsp_ip->dsp_buf = (short unsigned int *)dsp_ip->dsp_mp->apbuf;
4059  } else {
4060  dsp_ip->dsp_buf = (short unsigned int *)dsp_ip->dsp_mp->buf;
4061  }
4062  return 0;
4063 }
4064 
4065 
4066 /**
4067  * Retrieve data for DSP from a database object.
4068  */
4069 HIDDEN int
4070 get_obj_data(struct rt_dsp_internal *dsp_ip, const struct db_i *dbip)
4071 {
4072  struct rt_binunif_internal *bip;
4073  int in_cookie, out_cookie;
4074  size_t got;
4075  int ret;
4076 
4077  BU_ALLOC(dsp_ip->dsp_bip, struct rt_db_internal);
4078 
4079  ret = rt_retrieve_binunif (dsp_ip->dsp_bip, dbip, bu_vls_addr(&dsp_ip->dsp_name));
4080  if (ret)
4081  return -1;
4082 
4083  if (RT_G_DEBUG & DEBUG_HF) {
4084  bu_log("db_internal magic: 0x%08x major: %d minor: %d\n",
4085  dsp_ip->dsp_bip->idb_magic,
4086  dsp_ip->dsp_bip->idb_major_type,
4087  dsp_ip->dsp_bip->idb_minor_type);
4088  }
4089 
4090  bip = (struct rt_binunif_internal *)dsp_ip->dsp_bip->idb_ptr;
4091 
4092  if (RT_G_DEBUG & DEBUG_HF)
4093  bu_log("binunif magic: 0x%08x type: %d count:%zu data[0]:%u\n",
4094  bip->magic, bip->type, bip->count, bip->u.uint16[0]);
4095 
4096 
4097  in_cookie = bu_cv_cookie("nus"); /* data is network unsigned short */
4098  out_cookie = bu_cv_cookie("hus");
4099 
4100  if (bu_cv_optimize(in_cookie) != bu_cv_optimize(out_cookie)) {
4101  /* if we're on a little-endian machine we convert the input
4102  * file from network to host format
4103  */
4104 
4105  got = bu_cv_w_cookie(bip->u.uint16, out_cookie,
4106  bip->count * sizeof(unsigned short),
4107  bip->u.uint16, in_cookie, bip->count);
4108 
4109  if (got != bip->count) {
4110  bu_log("got %zu != count %zu", got, bip->count);
4111  bu_bomb("\n");
4112  }
4113  }
4114 
4115  dsp_ip->dsp_buf = bip->u.uint16;
4116  return 0;
4117 }
4118 
4119 
4120 /**
4121  * Handle things common to both the v4 and v5 database.
4122  *
4123  * This includes applying the modelling transform, and fetching the
4124  * actual data.
4125  *
4126  * Return:
4127  * 0 success
4128  * !0 failure
4129  */
4130 HIDDEN int
4131 dsp_get_data(struct rt_dsp_internal *dsp_ip, const mat_t mat, const struct db_i *dbip)
4132 {
4133  mat_t tmp;
4134  char *p;
4135 
4136  /* Apply Modeling transform */
4137  MAT_COPY(tmp, dsp_ip->dsp_stom);
4138  bn_mat_mul(dsp_ip->dsp_stom, mat, tmp);
4139 
4140  bn_mat_inv(dsp_ip->dsp_mtos, dsp_ip->dsp_stom);
4141 
4142  p = bu_vls_addr(&dsp_ip->dsp_name);
4143 
4144  switch (dsp_ip->dsp_datasrc) {
4145  case RT_DSP_SRC_FILE:
4146  case RT_DSP_SRC_V4_FILE:
4147  /* Retrieve the data from an external file */
4148  if (RT_G_DEBUG & DEBUG_HF)
4149  bu_log("getting data from file \"%s\"\n", p);
4150 
4151  if (get_file_data(dsp_ip, dbip) != 0) {
4152  p = "file";
4153  } else {
4154  return 0;
4155  }
4156 
4157  break;
4158  case RT_DSP_SRC_OBJ:
4159  /* Retrieve the data from an internal db object */
4160  if (RT_G_DEBUG & DEBUG_HF)
4161  bu_log("getting data from object \"%s\"\n", p);
4162 
4163  if (get_obj_data(dsp_ip, dbip) != 0) {
4164  p = "object";
4165  } else {
4166  RT_CK_DB_INTERNAL(dsp_ip->dsp_bip);
4167  RT_CK_BINUNIF(dsp_ip->dsp_bip->idb_ptr);
4168  return 0;
4169  }
4170  break;
4171  default:
4172  bu_log("%s:%d Odd dsp data src '%c' s/b '%c' or '%c'\n",
4173  __FILE__, __LINE__, dsp_ip->dsp_datasrc,
4174  RT_DSP_SRC_FILE, RT_DSP_SRC_OBJ);
4175  return -1;
4176  }
4177 
4178  bu_log("Cannot retrieve DSP data from %s \"%s\"\n", p,
4179  bu_vls_addr(&dsp_ip->dsp_name));
4180 
4181  dsp_ip->dsp_mp = (struct bu_mapped_file *)NULL;
4182  dsp_ip->dsp_buf = NULL;
4183 
4184  return 1;
4185 }
4186 
4187 
4188 /**
4189  * Import an DSP from the database format to the internal format.
4190  * Apply modeling transformations as well.
4191  */
4192 int
4193 rt_dsp_import4(struct rt_db_internal *ip, const struct bu_external *ep, register const fastf_t *mat, const struct db_i *dbip)
4194 {
4195  struct rt_dsp_internal *dsp_ip;
4196  union record *rp;
4197  struct bu_vls str = BU_VLS_INIT_ZERO;
4198 
4199  if (RT_G_DEBUG & DEBUG_HF)
4200  bu_log("rt_dsp_import4_v4()\n");
4201 
4202  BU_CK_EXTERNAL(ep);
4203  rp = (union record *)ep->ext_buf;
4204 
4205  if (RT_G_DEBUG & DEBUG_HF)
4206  bu_log("rt_dsp_import4(%s)\n", rp->ss.ss_args);
4207  /*----------------------------------------------------------------------*/
4208 
4209 
4210  /* Check record type */
4211  if (rp->u_id != DBID_STRSOL) {
4212  bu_log("rt_dsp_import4: defective record\n");
4213  return -1;
4214  }
4215 
4216  RT_CK_DB_INTERNAL(ip);
4217  ip->idb_major_type = DB5_MAJORTYPE_BRLCAD;
4218  ip->idb_type = ID_DSP;
4219  ip->idb_meth = &OBJ[ID_DSP];
4220  BU_ALLOC(ip->idb_ptr, struct rt_dsp_internal);
4221 
4222  dsp_ip = (struct rt_dsp_internal *)ip->idb_ptr;
4223  dsp_ip->magic = RT_DSP_INTERNAL_MAGIC;
4224 
4225  /* set defaults */
4226  /* XXX bu_struct_parse does not set the null?
4227  * memset(&dsp_ip->dsp_name[0], 0, DSP_NAME_LEN);
4228  */
4229  dsp_ip->dsp_xcnt = dsp_ip->dsp_ycnt = 0;
4230 
4231  dsp_ip->dsp_cuttype = DSP_CUT_DIR_ADAPT;
4232  dsp_ip->dsp_smooth = 1;
4233  MAT_IDN(dsp_ip->dsp_stom);
4234  MAT_IDN(dsp_ip->dsp_mtos);
4235 
4236  bu_vls_strcpy(&str, rp->ss.ss_args);
4237  if (bu_struct_parse(&str, rt_dsp_parse, (char *)dsp_ip, NULL) < 0) {
4238  if (BU_VLS_IS_INITIALIZED(&str)) bu_vls_free(&str);
4239  IMPORT_FAIL("parse error");
4240  }
4241 
4242  /* validate size */
4243  if (dsp_ip->dsp_xcnt == 0 || dsp_ip->dsp_ycnt == 0) {
4244  IMPORT_FAIL("zero dimension on map");
4245  }
4246 
4247  if (mat == NULL) mat = bn_mat_identity;
4248  if (dsp_get_data(dsp_ip, mat, dbip)!=0) {
4249  IMPORT_FAIL("unable to load displacement map data");
4250  }
4251 
4252  if (RT_G_DEBUG & DEBUG_HF) {
4253  bu_vls_trunc(&str, 0);
4254  bu_vls_struct_print(&str, rt_dsp_ptab, (char *)dsp_ip);
4255  bu_log(" imported as(%s)\n", bu_vls_addr(&str));
4256 
4257  }
4258 
4259  if (BU_VLS_IS_INITIALIZED(&str)) bu_vls_free(&str);
4260 
4261  RT_CK_DB_INTERNAL(dsp_ip->dsp_bip);
4262  RT_CK_BINUNIF(dsp_ip->dsp_bip->idb_ptr);
4263 
4264  return 0; /* OK */
4265 }
4266 
4267 
4268 /**
4269  * The name is added by the caller, in the usual place.
4270  */
4271 int
4272 rt_dsp_export4(struct bu_external *ep, const struct rt_db_internal *ip, double local2mm, const struct db_i *dbip)
4273 {
4274  struct rt_dsp_internal *dsp_ip;
4275  struct rt_dsp_internal dsp;
4276  union record *rec;
4277  struct bu_vls str = BU_VLS_INIT_ZERO;
4278 
4279  if (dbip) RT_CK_DBI(dbip);
4280 
4281  RT_CK_DB_INTERNAL(ip);
4282  if (ip->idb_type != ID_DSP) return -1;
4283  dsp_ip = (struct rt_dsp_internal *)ip->idb_ptr;
4284  RT_DSP_CK_MAGIC(dsp_ip);
4285  BU_CK_VLS(&dsp_ip->dsp_name);
4286 
4287  BU_CK_EXTERNAL(ep);
4288  ep->ext_nbytes = sizeof(union record)*DB_SS_NGRAN;
4289  ep->ext_buf = (uint8_t *)bu_calloc(1, ep->ext_nbytes, "dsp external");
4290  rec = (union record *)ep->ext_buf;
4291 
4292  dsp = *dsp_ip; /* struct copy */
4293 
4294  /* Since libwdb users may want to operate in units other than mm,
4295  * we offer the opportunity to scale the solid (to get it into mm)
4296  * on the way out.
4297  */
4298  dsp.dsp_stom[15] *= local2mm;
4299 
4300  bu_vls_struct_print(&str, rt_dsp_ptab, (char *)&dsp);
4301  if (RT_G_DEBUG & DEBUG_HF)
4302  bu_log("rt_dsp_export4_v4(%s)\n", bu_vls_addr(&str));
4303 
4304  rec->ss.ss_id = DBID_STRSOL;
4305  bu_strlcpy(rec->ss.ss_keyword, "dsp", sizeof(rec->ss.ss_keyword));
4306  bu_strlcpy(rec->ss.ss_args, bu_vls_addr(&str), DB_SS_LEN);
4307 
4308 
4309  if (BU_VLS_IS_INITIALIZED(&str)) bu_vls_free(&str);
4310 
4311  return 0;
4312 }
4313 
4314 
4315 /**
4316  * Import an DSP from the database format to the internal format.
4317  * Apply modeling transformations as well.
4318  */
4319 int
4320 rt_dsp_import5(struct rt_db_internal *ip, const struct bu_external *ep, register const fastf_t *mat, const struct db_i *dbip, struct resource *resp)
4321 {
4322  struct rt_dsp_internal *dsp_ip;
4323  unsigned char *cp;
4324 
4325  /* must be double for import and export */
4326  double scanmat[16];
4327 
4328  if (resp) RT_CK_RESOURCE(resp);
4329 
4330  if (RT_G_DEBUG & DEBUG_HF)
4331  bu_log("rt_dsp_import4_v5()\n");
4332 
4333 
4334  BU_CK_EXTERNAL(ep);
4335 
4336  BU_ASSERT_LONG(ep->ext_nbytes, >, 141);
4337 
4338  RT_CK_DB_INTERNAL(ip);
4339 
4340  ip->idb_major_type = DB5_MAJORTYPE_BRLCAD;
4341  ip->idb_type = ID_DSP;
4342  ip->idb_meth = &OBJ[ID_DSP];
4343  BU_ALLOC(ip->idb_ptr, struct rt_dsp_internal);
4344 
4345  dsp_ip = (struct rt_dsp_internal *)ip->idb_ptr;
4346  BU_VLS_INIT(&dsp_ip->dsp_name);
4347 
4348  dsp_ip->magic = RT_DSP_INTERNAL_MAGIC;
4349 
4350  /* get x, y counts */
4351  cp = (unsigned char *)ep->ext_buf;
4352 
4353  dsp_ip->dsp_xcnt = ntohl(*(uint32_t *)cp);
4354  cp += SIZEOF_NETWORK_LONG;
4355  if (dsp_ip->dsp_xcnt < 1) {
4356  bu_log("%s:%d DSP X dimension (%u) < 1 \n",
4357  __FILE__, __LINE__,
4358  dsp_ip->dsp_xcnt);
4359  }
4360 
4361  dsp_ip->dsp_ycnt = ntohl(*(uint32_t *)cp);
4362  cp += SIZEOF_NETWORK_LONG;
4363  if (dsp_ip->dsp_ycnt < 1) {
4364  bu_log("%s:%d DSP Y dimension (%u) < 1 \n",
4365  __FILE__, __LINE__,
4366  dsp_ip->dsp_ycnt);
4367  }
4368 
4369  /* validate size */
4370  if (dsp_ip->dsp_xcnt == 0 || dsp_ip->dsp_ycnt == 0) {
4371  IMPORT_FAIL("zero dimension on map");
4372  }
4373 
4374  /* convert matrix */
4375  bu_cv_ntohd((unsigned char *)scanmat, cp, ELEMENTS_PER_MAT);
4376  MAT_COPY(dsp_ip->dsp_stom, scanmat); /* double to fastf_t */
4377 
4378  cp += SIZEOF_NETWORK_DOUBLE * ELEMENTS_PER_MAT;
4379  bn_mat_inv(dsp_ip->dsp_mtos, dsp_ip->dsp_stom);
4380 
4381  /* convert smooth flag */
4382  dsp_ip->dsp_smooth = ntohs(*(uint16_t *)cp);
4383  cp += SIZEOF_NETWORK_SHORT;
4384 
4385  dsp_ip->dsp_datasrc = *cp;
4386  cp++;
4387  switch (dsp_ip->dsp_datasrc) {
4388  case RT_DSP_SRC_V4_FILE:
4389  case RT_DSP_SRC_FILE:
4390  case RT_DSP_SRC_OBJ:
4391  break;
4392  default:
4393  bu_log("%s:%d bogus DSP datasrc '%c' (%d)\n",
4394  __FILE__, __LINE__,
4395  dsp_ip->dsp_datasrc, dsp_ip->dsp_datasrc);
4396  break;
4397  }
4398 
4399  dsp_ip->dsp_cuttype = *cp;
4400  cp++;
4401  switch (dsp_ip->dsp_cuttype) {
4402  case DSP_CUT_DIR_ADAPT:
4403  case DSP_CUT_DIR_llUR:
4404  case DSP_CUT_DIR_ULlr:
4405  break;
4406  default:
4407  bu_log("%s:%d bogus DSP cut type '%c' (%d)\n",
4408  __FILE__, __LINE__,
4409  dsp_ip->dsp_cuttype, dsp_ip->dsp_cuttype);
4410  break;
4411  }
4412 
4413  /* convert name of data location */
4414  bu_vls_init(&dsp_ip->dsp_name);
4415  bu_vls_strncpy(&dsp_ip->dsp_name, (char *)cp,
4416  ep->ext_nbytes - (cp - (unsigned char *)ep->ext_buf));
4417 
4418  if (mat == NULL) mat = bn_mat_identity;
4419  if (dsp_get_data(dsp_ip, mat, dbip) != 0) {
4420  IMPORT_FAIL("unable to load displacement map data");
4421  }
4422 
4423  return 0; /* OK */
4424 }
4425 
4426 
4427 /**
4428  * The name is added by the caller, in the usual place.
4429  */
4430 int
4431 rt_dsp_export5(struct bu_external *ep, const struct rt_db_internal *ip, double local2mm, const struct db_i *dbip, struct resource *resp)
4432 {
4433  struct rt_dsp_internal *dsp_ip;
4434  unsigned long name_len;
4435  unsigned char *cp;
4436  size_t rem;
4437 
4438  /* must be double for import and export */
4439  double scanmat[16];
4440 
4441  if (resp) RT_CK_RESOURCE(resp);
4442  if (dbip) RT_CK_DBI(dbip);
4443 
4444  RT_CK_DB_INTERNAL(ip);
4445  if (ip->idb_type != ID_DSP) return -1;
4446  dsp_ip = (struct rt_dsp_internal *)ip->idb_ptr;
4447  RT_DSP_CK_MAGIC(dsp_ip);
4448 
4449  if (RT_G_DEBUG & DEBUG_HF)
4450  bu_log("rt_dsp_export4_v5()\n");
4451 
4452  name_len = bu_vls_strlen(&dsp_ip->dsp_name) + 1;
4453 
4454  BU_CK_EXTERNAL(ep);
4455 
4456  ep->ext_nbytes =
4457  SIZEOF_NETWORK_LONG * 2 +
4458  SIZEOF_NETWORK_DOUBLE * ELEMENTS_PER_MAT +
4460  2 + name_len;
4461 
4462  ep->ext_buf = (uint8_t *)bu_malloc(ep->ext_nbytes, "dsp external");
4463  cp = (unsigned char *)ep->ext_buf;
4464  rem = ep->ext_nbytes;
4465 
4466  memset(ep->ext_buf, 0, ep->ext_nbytes);
4467 
4468  /* Now we fill the buffer with the data, making sure everything is
4469  * converted to Big-Endian IEEE
4470  */
4471 
4472  *(uint32_t *)cp = htonl((uint32_t)dsp_ip->dsp_xcnt);
4473  cp += SIZEOF_NETWORK_LONG;
4474  rem -= SIZEOF_NETWORK_LONG;
4475 
4476  *(uint32_t *)cp = htonl((uint32_t)dsp_ip->dsp_ycnt);
4477  cp += SIZEOF_NETWORK_LONG;
4478  rem -= SIZEOF_NETWORK_LONG;
4479 
4480  /* Since libwdb users may want to operate in units other than mm,
4481  * we offer the opportunity to scale the solid (to get it into mm)
4482  * on the way out.
4483  */
4484  dsp_ip->dsp_stom[15] *= local2mm;
4485 
4486  MAT_COPY(scanmat, dsp_ip->dsp_stom); /* convert fastf_t to double */
4487  bu_cv_htond(cp, (unsigned char *)scanmat, ELEMENTS_PER_MAT);
4488 
4489  cp += SIZEOF_NETWORK_DOUBLE * ELEMENTS_PER_MAT;
4490  rem -= SIZEOF_NETWORK_DOUBLE * ELEMENTS_PER_MAT;
4491 
4492  *(uint16_t *)cp = htons((uint16_t)dsp_ip->dsp_smooth);
4493  cp += SIZEOF_NETWORK_SHORT;
4494  rem -= SIZEOF_NETWORK_SHORT;
4495 
4496  switch (dsp_ip->dsp_datasrc) {
4497  case RT_DSP_SRC_V4_FILE:
4498  case RT_DSP_SRC_FILE:
4499  case RT_DSP_SRC_OBJ:
4500  *cp = dsp_ip->dsp_datasrc;
4501  break;
4502  default:
4503  *cp = RT_DSP_SRC_FILE;
4504  break;
4505  }
4506  cp++;
4507  rem--;
4508 
4509  switch (dsp_ip->dsp_cuttype) {
4510  case DSP_CUT_DIR_ADAPT:
4511  case DSP_CUT_DIR_llUR:
4512  case DSP_CUT_DIR_ULlr:
4513  *cp = dsp_ip->dsp_cuttype;
4514  break;
4515  default:
4516  *cp = DSP_CUT_DIR_ADAPT;
4517  break;
4518  }
4519  cp++;
4520  rem--;
4521 
4522  bu_strlcpy((char *)cp, bu_vls_addr(&dsp_ip->dsp_name), rem);
4523 
4524  return 0; /* OK */
4525 }
4526 
4527 
4528 /**
4529  * Make human-readable formatted presentation of this solid. First
4530  * line describes type of solid. Additional lines are indented one
4531  * tab, and give parameter values.
4532  */
4533 int
4535  const struct rt_db_internal *ip,
4536  int UNUSED(verbose),
4537  double UNUSED(mm2local),
4538  struct resource *resp,
4539  struct db_i *db_ip)
4540 {
4541  register struct rt_dsp_internal *dsp_ip =
4542  (struct rt_dsp_internal *)ip->idb_ptr;
4543  struct bu_vls vls = BU_VLS_INIT_ZERO;
4544 
4545  if (resp) RT_CK_RESOURCE(resp);
4546  if (db_ip) RT_CK_DBI(db_ip);
4547 
4548  if (RT_G_DEBUG & DEBUG_HF)
4549  bu_log("rt_dsp_describe()\n");
4550 
4551  RT_DSP_CK_MAGIC(dsp_ip);
4552 
4553  switch (db_version(db_ip)) {
4554  case 4:
4555  dsp_print_v4(&vls, dsp_ip);
4556  break;
4557  case 5:
4558  dsp_print_v5(&vls, dsp_ip);
4559  break;
4560  }
4561 
4562  bu_vls_vlscat(str, &vls);
4563 
4564  if (BU_VLS_IS_INITIALIZED(&vls)) bu_vls_free(&vls);
4565 
4566  return 0;
4567 }
4568 
4569 
4570 /**
4571  * Free the storage associated with the rt_db_internal version of this
4572  * solid.
4573  */
4574 void
4576 {
4577  register struct rt_dsp_internal *dsp_ip;
4578 
4579  if (RT_G_DEBUG & DEBUG_HF)
4580  bu_log("rt_dsp_ifree()\n");
4581 
4582  RT_CK_DB_INTERNAL(ip);
4583 
4584  dsp_ip = (struct rt_dsp_internal *)ip->idb_ptr;
4585  RT_DSP_CK_MAGIC(dsp_ip);
4586 
4587  if (dsp_ip->dsp_mp) {
4588  bu_close_mapped_file(dsp_ip->dsp_mp);
4589  }
4590 
4591  if (dsp_ip->dsp_bip) {
4592  dsp_ip->dsp_bip->idb_meth->ft_ifree((struct rt_db_internal *) dsp_ip->dsp_bip);
4593  }
4594 
4595  dsp_ip->magic = 0; /* sanity */
4596  dsp_ip->dsp_mp = (struct bu_mapped_file *)0;
4597 
4598  if (BU_VLS_IS_INITIALIZED(&dsp_ip->dsp_name))
4599  bu_vls_free(&dsp_ip->dsp_name);
4600  else
4601  bu_log("Freeing Bogus DSP, VLS string not initialized\n");
4602 
4603 
4604  bu_free((char *)dsp_ip, "dsp ifree");
4605  ip->idb_ptr = ((void *)0); /* sanity */
4606 }
4607 
4608 
4609 HIDDEN void
4610 hook_verify(const struct bu_structparse *sp,
4611  const char *sp_name,
4612  void *base,
4613  const char *UNUSED(p),
4614  void *UNUSED(data))
4615 {
4616  struct rt_dsp_internal *dsp_ip = (struct rt_dsp_internal *)base;
4617 
4618  if (!sp || !sp_name || !base) return;
4619 
4620  if (BU_STR_EQUAL(sp_name, "src")) {
4621  switch (dsp_ip->dsp_datasrc) {
4622  case RT_DSP_SRC_V4_FILE:
4623  case RT_DSP_SRC_FILE:
4624  case RT_DSP_SRC_OBJ:
4625  break;
4626  default:
4627  bu_log("Error in DSP data source field s/b one of [%c%c%c]\n",
4628  RT_DSP_SRC_V4_FILE,
4629  RT_DSP_SRC_FILE,
4630  RT_DSP_SRC_OBJ);
4631  break;
4632  }
4633 
4634  } else if (BU_STR_EQUAL(sp_name, "w")) {
4635  if (dsp_ip->dsp_xcnt == 0)
4636  bu_log("Error in DSP width dimension (0)\n");
4637  } else if (BU_STR_EQUAL(sp_name, "n")) {
4638  if (dsp_ip->dsp_ycnt == 0)
4639  bu_log("Error in DSP width dimension (0)\n");
4640  } else if (BU_STR_EQUAL(sp_name, "cut")) {
4641  switch (dsp_ip->dsp_cuttype) {
4642  case DSP_CUT_DIR_ADAPT:
4643  case DSP_CUT_DIR_llUR:
4644  case DSP_CUT_DIR_ULlr:
4645  break;
4646  default:
4647  bu_log("Error in DSP cut type: %c s/b one of [%c%c%c]\n",
4648  dsp_ip->dsp_cuttype,
4649  DSP_CUT_DIR_ADAPT,
4650  DSP_CUT_DIR_llUR,
4651  DSP_CUT_DIR_ULlr);
4652  break;
4653  }
4654  }
4655 }
4656 
4657 
4659  {"%V", 1, "file", DSP_O(dsp_name), hook_file, NULL, NULL },
4660  {"%c", 1, "src", DSP_O(dsp_datasrc), hook_verify, NULL, NULL },
4661  {"%V", 1, "name", DSP_O(dsp_name), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL },
4662  {"%d", 1, "w", DSP_O(dsp_xcnt), hook_verify, NULL, NULL },
4663  {"%d", 1, "n", DSP_O(dsp_ycnt), hook_verify, NULL, NULL },
4664  {"%i", 1, "sm", DSP_O(dsp_smooth), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL },
4665  {"%c", 1, "cut", DSP_O(dsp_cuttype), hook_verify, NULL, NULL },
4666  {"%f", 16, "stom", DSP_O(dsp_stom), hook_verify, NULL, NULL },
4667  {"", 0, (char *)0, 0, BU_STRUCTPARSE_FUNC_NULL, NULL, NULL }
4668 };
4669 
4670 
4671 /**
4672  * This is the generic routine to be listed in OBJ[].ft_get for
4673  * those solid types which are fully described by their ft_parsetab
4674  * entry.
4675  *
4676  * 'attr' is specified to retrieve only one attribute, rather than
4677  * all.
4678  *
4679  * Example: "db get ell.s B" to get only the B vector.
4680  */
4681 int
4682 rt_dsp_get(struct bu_vls *logstr, const struct rt_db_internal *intern, const char *attr)
4683 {
4684  register const struct bu_structparse *sp = NULL;
4685  const struct rt_dsp_internal *dsp_ip;
4686 
4687  /* XXX if dsp_datasrc == RT_DSP_SRC_V4_FILE we have a V4 dsp
4688  * otherwise, a V5 dsp. Take advantage of this.
4689  */
4690 
4691  RT_CK_DB_INTERNAL(intern);
4692  dsp_ip = (struct rt_dsp_internal *)intern->idb_ptr;
4693 
4694  if (attr == (char *)0) {
4695  /* Print out solid type and all attributes */
4696 
4697  bu_vls_printf(logstr, "dsp");
4698 
4699  switch (dsp_ip->dsp_datasrc) {
4700  case RT_DSP_SRC_V4_FILE:
4701  sp = rt_dsp_ptab;
4702  break;
4703  case RT_DSP_SRC_FILE:
4704  case RT_DSP_SRC_OBJ:
4705  sp = fake_dsp_printab;
4706  break;
4707  }
4708 
4709  while (sp && sp->sp_name != NULL) {
4710  bu_vls_printf(logstr, " %s {", sp->sp_name);
4711  bu_vls_struct_item(logstr, sp, (char *)dsp_ip, ' ');
4712  bu_vls_printf(logstr, "}");
4713  ++sp;
4714  }
4715 
4716  return BRLCAD_OK;
4717  }
4718 
4719  switch (dsp_ip->dsp_datasrc) {
4720  case RT_DSP_SRC_V4_FILE:
4721  sp = rt_dsp_ptab;
4722  break;
4723  case RT_DSP_SRC_FILE:
4724  case RT_DSP_SRC_OBJ:
4725  sp = fake_dsp_printab;
4726  break;
4727  }
4728 
4729  if (bu_vls_struct_item_named(logstr, sp, attr,
4730  (char *)dsp_ip, ' ') < 0) {
4731  bu_vls_printf(logstr,
4732  "Objects of type %s do not have a %s attribute.",
4733  "dsp", attr);
4734  return BRLCAD_ERROR;
4735  }
4736 
4737  return BRLCAD_OK;
4738 }
4739 
4740 
4741 /**
4742  * For those solids entirely defined by their parsetab. Invoked via
4743  * OBJ[].ft_adjust()
4744  */
4745 int
4746 rt_dsp_adjust(struct bu_vls *logstr, struct rt_db_internal *intern, int argc, const char **argv)
4747 {
4748  register const struct bu_structparse *sp = NULL;
4749  const struct rt_dsp_internal *dsp_ip;
4750 
4751  RT_CK_DB_INTERNAL(intern);
4752  dsp_ip = (struct rt_dsp_internal *)intern->idb_ptr;
4753  RT_DSP_CK_MAGIC(dsp_ip);
4754  BU_CK_VLS(&dsp_ip->dsp_name);
4755 
4756  switch (dsp_ip->dsp_datasrc) {
4757  case RT_DSP_SRC_V4_FILE:
4758  sp = rt_dsp_ptab;
4759  break;
4760  case RT_DSP_SRC_FILE:
4761  case RT_DSP_SRC_OBJ:
4762  sp = fake_dsp_printab;
4763  break;
4764  }
4765 
4766  if (! sp) return BRLCAD_ERROR;
4767 
4768  return bu_structparse_argv(logstr, argc, argv, sp, (char *)intern->idb_ptr, NULL);
4769 }
4770 
4771 
4772 void
4773 rt_dsp_make(const struct rt_functab *ftp, struct rt_db_internal *intern)
4774 {
4775  struct rt_dsp_internal *dsp;
4776 
4777  intern->idb_major_type = DB5_MAJORTYPE_BRLCAD;
4778  intern->idb_type = ID_DSP;
4779 
4780  BU_ASSERT(&OBJ[intern->idb_type] == ftp);
4781  intern->idb_meth = ftp;
4782 
4783  BU_ALLOC(dsp, struct rt_dsp_internal);
4784 
4785  intern->idb_ptr = (void *)dsp;
4786  dsp->magic = RT_DSP_INTERNAL_MAGIC;
4787  bu_vls_init(&dsp->dsp_name);
4788  bu_vls_strcpy(&dsp->dsp_name, "/dev/null");
4789  dsp->dsp_cuttype = DSP_CUT_DIR_ADAPT;
4790  MAT_IDN(dsp->dsp_mtos);
4791  MAT_IDN(dsp->dsp_stom);
4792  dsp->dsp_datasrc = RT_DSP_SRC_FILE;
4793 
4794 }
4795 
4796 
4797 int
4798 rt_dsp_params(struct pc_pc_set *ps, const struct rt_db_internal *ip)
4799 {
4800  if (!ps) return 0;
4801  RT_CK_DB_INTERNAL(ip);
4802 
4803  return 0; /* OK */
4804 }
4805 
4806 
4807 HIDDEN int
4809  int B[3],
4810  int C[3],
4811  int D[3],
4812  struct dsp_specific *dsp)
4813 
4814 {
4815  switch (dsp->dsp_i.dsp_cuttype) {
4816  case DSP_CUT_DIR_llUR:
4817  return 0;
4818  break;
4819 
4820  case DSP_CUT_DIR_ADAPT: {
4821  int lo[2], hi[2];
4822  fastf_t h1, h2, h3, h4;
4823  fastf_t cAD, cBC; /* curvature in direction AD, and BC */
4824 
4825 
4826  /*
4827  * We look at the points in the diagonal next cells to
4828  * determine the curvature along each diagonal of this
4829  * cell. This cell is divided into two triangles by
4830  * cutting across the cell in the direction of least
4831  * curvature.
4832  *
4833  * * * * *
4834  * \ /
4835  * \C D/
4836  * * *--* *
4837  * |\/|
4838  * |/\|
4839  * * *--* *
4840  * /A B\
4841  * / \
4842  * * * * *
4843  */
4844 
4845  lo[X] = A[X] - 1;
4846  lo[Y] = A[Y] - 1;
4847  hi[X] = D[X] + 1;
4848  hi[Y] = D[Y] + 1;
4849 
4850  /* a little bounds checking */
4851  if (lo[X] < 0) lo[X] = 0;
4852  if (lo[Y] < 0) lo[Y] = 0;
4853  if (hi[X] > dsp->xsiz)
4854  hi[X] = dsp->xsiz;
4855 
4856  if (hi[Y] > dsp->ysiz)
4857  hi[Y] = dsp->ysiz;
4858 
4859  /* compute curvature along the A->D direction */
4860  h1 = DSP(&dsp->dsp_i, lo[X], lo[Y]);
4861  h2 = A[Z];
4862  h3 = D[Z];
4863  h4 = DSP(&dsp->dsp_i, hi[X], hi[Y]);
4864 
4865  cAD = fabs(h3 + h1 - 2*h2) + fabs(h4 + h2 - 2*h3);
4866 
4867 
4868  /* compute curvature along the B->C direction */
4869  h1 = DSP(&dsp->dsp_i, hi[X], lo[Y]);
4870  h2 = B[Z];
4871  h3 = C[Z];
4872  h4 = DSP(&dsp->dsp_i, lo[X], hi[Y]);
4873 
4874  cBC = fabs(h3 + h1 - 2*h2) + fabs(h4 + h2 - 2*h3);
4875 
4876  if (cAD < cBC) {
4877  /* A-D cut is fine, no need to permute */
4878  if (RT_G_DEBUG & DEBUG_HF)
4879  bu_log("A-D cut (no swap)\n");
4880 
4881  return 0;
4882  }
4883 
4884  /* fallthrough */
4885 
4886  }
4887  case DSP_CUT_DIR_ULlr: {
4888  /* prefer the B-C cut */
4889  int tmp[3];
4890 
4891  VMOVE(tmp, A);
4892  VMOVE(A, B);
4893  VMOVE(B, tmp);
4894 
4895  VMOVE(tmp, D);
4896  VMOVE(D, C);
4897  VMOVE(C, tmp);
4898 
4899  if (RT_G_DEBUG & DEBUG_HF)
4900  bu_log("B-C cut (swap)\n");
4901  }
4902  return 0;
4903  break;
4904  }
4905  bu_log("%s:%d Unknown DSP cut direction: %d\n",
4906  __FILE__, __LINE__, dsp->dsp_i.dsp_cuttype);
4907  bu_bomb("");
4908  /* not reached */
4909  return -1;
4910 }
4911 
4912 
4913 /**
4914  * triangle with origin at A, X axis along AB, Y axis along AC
4915  *
4916  * B--A C A--B C
4917  * | | | |
4918  * C A--B C B--A
4919  */
4920 int
4921 project_pt(point_t out,
4922  int A[3],
4923  int B[3],
4924  int C[3],
4925  point_t pt)
4926 {
4927  int dx, dy;
4928  fastf_t alpha, beta, x, y;
4929  vect_t AB, AC;
4930 
4931  if (RT_G_DEBUG & DEBUG_HF) {
4932  bu_log(" pt %g %g %g\n", V3ARGS(pt));
4933  bu_log(" origin %d %d %d\n", V3ARGS(A));
4934  bu_log(" Bpt %d %d %d\n", V3ARGS(B));
4935  bu_log(" Cpt %d %d %d\n", V3ARGS(C));
4936  }
4937  if ((B[X] - A[X]) < 0) dx = -1;
4938  else dx = 1;
4939  if ((C[Y] - A[Y]) < 0) dy = -1;
4940  else dy = 1;
4941 
4942  x = pt[X] - A[X];
4943  y = pt[Y] - A[Y];
4944 
4945  alpha = x * dx;
4946  beta = y * dy;
4947  if (RT_G_DEBUG & DEBUG_HF) {
4948  bu_log("alpha:%g beta:%g\n", alpha, beta);
4949  }
4950  if (alpha < -SMALL_FASTF) {
4951  bu_log("Alpha negative: %g x:%g dx:%d\n", alpha, x, dx);
4952  }
4953  if (beta < -SMALL_FASTF) {
4954  bu_log("Beta negative: %g x:%g dx:%d\n", beta, y, dy);
4955  }
4956  CLAMP(alpha, 0.0, 1.0);
4957  CLAMP(beta, 0.0, 1.0);
4958 
4959  if (RT_G_DEBUG & DEBUG_HF) {
4960  bu_log("x:%g y:%g dx:%d dy:%d alpha:%g beta:%g\n",
4961  x, y, dx, dy, alpha, beta);
4962  }
4963  if ((alpha+beta) > (1.0+SMALL_FASTF)) {
4964  if (RT_G_DEBUG & DEBUG_HF) bu_log("Not this triangle\n");
4965  return 1;
4966  }
4967 
4968  VSUB2(AB, B, A);
4969  VSUB2(AC, C, A);
4970 
4971  VJOIN2(out, A, alpha, AB, beta, AC);
4972  if (RT_G_DEBUG & DEBUG_HF) bu_log("out: %g %g %g\n", V3ARGS(out));
4973 
4974  return 0;
4975 }
4976 
4977 
4978 /**
4979  * Given an arbitrary point return the projection of that point onto
4980  * the surface of the DSP. If the point is outside the bounds of the
4981  * DSP then it will be projected to the nearest edge of the DSP.
4982  * Return the cell number and the height above/below the plane
4983  */
4984 int
4985 dsp_pos(point_t out, /* return value */
4986  struct soltab *stp, /* pointer to soltab for dsp */
4987  point_t p) /* model space point */
4988 {
4989  struct dsp_specific *dsp;
4990  point_t pt, tri_pt;
4991  unsigned int x, y;
4992  int A[3], B[3], C[3], D[3];
4993 
4994  /* init points */
4995  VSET(pt, 0, 0, 0);
4996  VSET(tri_pt, 0, 0, 0);
4997 
4998  RT_CK_SOLTAB(stp);
4999 
5000  if (stp->st_id != ID_DSP) return 1;
5001 
5002  dsp = (struct dsp_specific *)stp->st_specific;
5003 
5004  MAT4X3PNT(pt, dsp->dsp_i.dsp_mtos, p);
5005 
5006 
5007  if (RT_G_DEBUG & DEBUG_HF) {
5008  VPRINT("user_point", p);
5009  VPRINT("dsp_point", pt);
5010  }
5011  x = pt[X];
5012  y = pt[Y];
5013 
5014  V_MIN(x, XSIZ(dsp)-1);
5015  V_MIN(y, YSIZ(dsp)-1);
5016 
5017  if (RT_G_DEBUG & DEBUG_HF)
5018  bu_log("x:%d y:%d\n", x, y);
5019 
5020  VSET(A, x, y, DSP(&dsp->dsp_i, x, y));
5021  x += 1;
5022  VSET(B, x, y, DSP(&dsp->dsp_i, x, y));
5023  y += 1;
5024  VSET(D, x, y, DSP(&dsp->dsp_i, x, y));
5025  x -= 1;
5026  VSET(C, x, y, DSP(&dsp->dsp_i, x, y));
5027  y -= 1;
5028 
5029  if (RT_G_DEBUG & DEBUG_HF) {
5030  bu_log(" A: %d %d %d\n", V3ARGS(A));
5031  bu_log(" B: %d %d %d\n", V3ARGS(B));
5032  bu_log(" C: %d %d %d\n", V3ARGS(C));
5033  bu_log(" D: %d %d %d\n", V3ARGS(D));
5034  }
5035 
5036  swap_cell_pts(A, B, C, D, dsp);
5037  if (RT_G_DEBUG & DEBUG_HF) {
5038  bu_log(" A: %d %d %d\n", V3ARGS(A));
5039  bu_log(" B: %d %d %d\n", V3ARGS(B));
5040  bu_log(" C: %d %d %d\n", V3ARGS(C));
5041  bu_log(" D: %d %d %d\n", V3ARGS(D));
5042  }
5043 
5044  if (project_pt(tri_pt, B, A, D, pt)) {
5045  if (project_pt(tri_pt, C, D, A, pt)) {
5046  bu_log("Now what???\n");
5047  }
5048  }
5049 
5050  MAT4X3PNT(out, dsp->dsp_i.dsp_stom, tri_pt);
5051  if (RT_G_DEBUG & DEBUG_HF) {
5052  VPRINT("user_pt", p);
5053  VPRINT("tri_pt", tri_pt);
5054  VPRINT("model_space", out);
5055  bu_log("X: %d Y:%d\n", x, y);
5056  }
5057 
5058  return 0;
5059 }
5060 
5061 
5062 /* Important when concatenating source files together */
5063 #undef dlog
5064 #undef XMIN
5065 #undef XMAX
5066 #undef YMIN
5067 #undef YMAX
5068 #undef ZMIN
5069 #undef ZMAX
5070 #undef ZMID
5071 #undef DSP
5072 #undef XCNT
5073 #undef YCNT
5074 #undef XSIZ
5075 #undef YSIZ
5076 
5077 /** @} */
5078 /*
5079  * Local Variables:
5080  * mode: C
5081  * tab-width: 8
5082  * indent-tabs-mode: t
5083  * c-file-style: "stroustrup"
5084  * End:
5085  * ex: shiftwidth=4 tabstop=8
5086  */
int num_segs
Definition: dsp.c:192
void bu_vls_init(struct bu_vls *vp)
Definition: vls.c:56
struct xray a_ray
Actual ray to be shot.
Definition: raytrace.h:1583
#define BU_LIST_FOR(p, structure, hp)
Definition: list.h:365
void nmg_lu_reorient(struct loopuse *lu)
Definition: nmg_mod.c:3649
Definition: raytrace.h:800
uint32_t hit_magic
Definition: raytrace.h:249
#define ZMIN
Definition: dsp.c:146
void bu_log(const char *,...) _BU_ATTR_PRINTF12
Definition: log.c:176
#define BU_LIST_INSERT(old, new)
Definition: list.h:183
#define SIZEOF_NETWORK_DOUBLE
Definition: cv.h:48
struct faceuse * nmg_cface(struct shell *s, struct vertex **verts, int n)
Definition: nmg_mod.c:1130
const char * sp_name
Definition: parse.h:140
struct dsp_bb * dspb_children[NUM_BB_CHILDREN]
Definition: dsp.c:110
uint32_t magic
Definition: dsp.c:99
struct hit seg_in
IN information.
Definition: raytrace.h:370
#define RT_DSP_INTERNAL_MAGIC
Definition: magic.h:88
fastf_t C[2 *MAX_CNT+1][2 *MAX_CNT+1]
Definition: dsp_brep.cpp:38
HIDDEN void plot_seg(struct isect_stuff *isect, struct hit *in_hit, struct hit *out_hit, const point_t bbmin, const point_t bbmax, int r, int g, int b)
Definition: dsp.c:1096
Definition: list.h:118
int nmg_fu_planeeqn(struct faceuse *fu, const struct bn_tol *tol)
Definition: nmg_mod.c:1311
struct dsp_specific * dsp
Definition: dsp.c:181
HIDDEN void dsp_print_v4(struct bu_vls *vls, const struct rt_dsp_internal *dsp_ip)
Definition: dsp.c:516
#define YMAX
Definition: dsp.c:145
#define BU_LIST_LAST(structure, hp)
Definition: list.h:306
struct dsp_rpp dspb_rpp
Definition: dsp.c:100
HIDDEN int get_file_data(struct rt_dsp_internal *dsp_ip, const struct db_i *dbip)
Definition: dsp.c:4016
int dmin
Definition: dsp.c:193
#define NUM_BB_CHILDREN
Definition: dsp.c:75
#define RT_CK_RTI(_p)
Definition: raytrace.h:1833
#define BU_VLS_IS_INITIALIZED(_vp)
Definition: vls.h:92
void pdv_3move(register FILE *plotfp, const fastf_t *pt)
Definition: plot3.c:618
double dist
>= 0
Definition: tol.h:73
void nmg_make_faces_within_tol(struct shell *s, const struct bn_tol *tol)
Definition: nmg_misc.c:8463
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 xray r
Definition: dsp.c:183
struct rt_dsp_internal dsp_i
Definition: dsp.c:159
struct soltab * seg_stp
pointer back to soltab
Definition: raytrace.h:372
HIDDEN void plot_dsp_bb(FILE *fp, struct dsp_bb *dsp_bb, struct dsp_specific *dsp, int r, int g, int b, int blather)
Definition: dsp.c:312
if lu s
Definition: nmg_mod.c:3860
Definition: clone.c:90
HIDDEN int get_obj_data(struct rt_dsp_internal *dsp_ip, const struct db_i *dbip)
Definition: dsp.c:4070
int nmg_join_touchingloops(struct loopuse *lu)
Definition: nmg_mod.c:2678
lu
Definition: nmg_mod.c:3855
#define VSET(a, b, c, d)
Definition: color.c:53
unsigned short * uint16
Definition: raytrace.h:1010
#define VSETALL(a, s)
Definition: color.c:54
Definition: raytrace.h:215
#define N
Definition: randmt.c:39
#define IMPORT_FAIL(_s)
Definition: dsp.c:77
int rt_dsp_import5(struct rt_db_internal *ip, const struct bu_external *ep, register const fastf_t *mat, const struct db_i *dbip, struct resource *resp)
Definition: dsp.c:4320
void bu_semaphore_acquire(unsigned int i)
Definition: semaphore.c:180
void bu_vls_trunc(struct bu_vls *vp, int len)
Definition: vls.c:198
#define BU_LIST_IS_EMPTY(hp)
Definition: list.h:295
#define SIZEOF_NETWORK_LONG
Definition: cv.h:46
Definition: pc.h:108
Definition: raytrace.h:368
HIDDEN int check_bbpt_hit_elev(int i, point_t A, point_t B, point_t C, point_t D, point_t P)
Definition: dsp.c:1680
#define BU_ASSERT_LONG(_lhs, _relation, _rhs)
Definition: defines.h:240
Definition: raytrace.h:248
void nmg_vertex_gv(struct vertex *v, const fastf_t *pt)
Definition: nmg_mk.c:1668
int db_version(struct db_i *dbip)
Definition: db5_scan.c:414
HIDDEN void hook_mtos_from_stom(const struct bu_structparse *sp, const char *sp_name, void *base, const char *p, void *data)
Definition: dsp.c:218
char ** dbi_filepath
search path for aux file opens (convenience var)
Definition: raytrace.h:810
#define SMALL_FASTF
Definition: defines.h:342
fastf_t st_aradius
Radius of APPROXIMATING sphere.
Definition: raytrace.h:433
void bu_vls_strncpy(struct bu_vls *vp, const char *s, size_t n)
Definition: vls.c:339
HIDDEN int dsp_get_data(struct rt_dsp_internal *dsp_ip, const mat_t mat, const struct db_i *dbip)
Definition: dsp.c:4131
unsigned short dsp_min[3]
Definition: dsp.c:88
Header file for the BRL-CAD common definitions.
#define DRAW(_pt)
#define ZMID
Definition: dsp.c:148
#define XMAX
Definition: dsp.c:143
HIDDEN void dsp_layers(struct dsp_specific *dsp, unsigned short *d_min, unsigned short *d_max)
Definition: dsp.c:650
struct resource * a_resource
dynamic memory resources
Definition: raytrace.h:1591
unsigned short dspb_ch_dim[2]
Definition: dsp.c:109
HIDDEN void plot_cell_top(struct isect_stuff *isect, struct dsp_bb *dsp_bb, point_t A, point_t B, point_t C, point_t D, struct hit hitlist[], int hitflags, int style)
Definition: dsp.c:428
#define BU_ASSERT(_equation)
Definition: defines.h:216
void bu_cv_htond(unsigned char *out, const unsigned char *in, size_t count)
int dsp_in_rpp(struct isect_stuff *isect, register const fastf_t *min, register const fastf_t *max)
Definition: dsp.c:2057
#define BU_LIST_NON_EMPTY(hp)
Definition: list.h:296
void bu_log_indent_delta(int delta)
Definition: log.c:54
int rt_dsp_prep(struct soltab *stp, struct rt_db_internal *ip, struct rt_i *rtip)
Definition: dsp.c:982
int a_onehit
flag to stop on first hit
Definition: raytrace.h:1586
#define MAX_FASTF
Definition: defines.h:340
int rt_dsp_import4(struct rt_db_internal *ip, const struct bu_external *ep, register const fastf_t *mat, const struct db_i *dbip)
Definition: dsp.c:4193
point_t min
Definition: raytrace.h:415
#define DEBUG_HF
27 Height Field solids
Definition: raytrace.h:114
union rt_binunif_internal::@9 u
#define HIDDEN
Definition: common.h:86
#define RT_CK_BINUNIF(_p)
Definition: raytrace.h:1016
void pd_3move(register FILE *plotfp, double x, double y, double z)
Definition: plot3.c:624
struct bu_list l
Definition: raytrace.h:369
HIDDEN int add_seg(struct isect_stuff *isect, struct hit *in_hit, struct hit *out_hit, const point_t bbmin, const point_t bbmax, int r, int g, int b)
Definition: dsp.c:1144
int rt_dsp_adjust(struct bu_vls *logstr, struct rt_db_internal *intern, int argc, const char **argv)
Definition: dsp.c:4746
void * bu_malloc(size_t siz, const char *str)
Definition: malloc.c:314
#define YSIZ(_p)
Definition: dsp.c:129
double rel
rel dist tol
Definition: raytrace.h:181
int out_surf
Definition: dsp.c:173
Definition: dsp.c:87
#define XCNT(_p)
Definition: dsp.c:126
if(share_geom)
Definition: nmg_mod.c:3829
size_t apbuflen
Definition: mapped_file.h:90
HIDDEN int swap_cell_pts(int A[3], int B[3], int C[3], int D[3], struct dsp_specific *dsp)
Definition: dsp.c:4808
int idb_major_type
Definition: raytrace.h:192
double dsp_pl_dist[BBOX_PLANES]
Definition: dsp.c:160
void nmg_split_touchingloops(struct loopuse *lu, const struct bn_tol *tol)
Definition: nmg_mod.c:2585
void rt_dsp_norm(register struct hit *hitp, struct soltab *stp, register struct xray *rp)
Definition: dsp.c:2807
void bu_vls_free(struct bu_vls *vp)
Definition: vls.c:248
const struct bu_structparse fake_dsp_printab[]
Definition: dsp.c:4658
Definition: color.c:49
int rt_dsp_export4(struct bu_external *ep, const struct rt_db_internal *ip, double local2mm, const struct db_i *dbip)
Definition: dsp.c:4272
struct bbox_isect minbox
Definition: dsp.c:190
COMPLEX data[64]
Definition: fftest.c:34
struct rt_i * a_rt_i
this librt instance
Definition: raytrace.h:1588
#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)
#define RT_G_DEBUG
Definition: raytrace.h:1718
int rt_dsp_describe(struct bu_vls *str, const struct rt_db_internal *ip, int verbose, double mm2local, struct resource *resp, struct db_i *db_ip)
Definition: dsp.c:4534
int rt_dsp_bbox(struct rt_db_internal *ip, point_t *min, point_t *max, const struct bn_tol *tol)
Definition: dsp.c:862
#define CHECK_VERTS(_v)
Definition: dsp.c:3500
vect_t hit_vpriv
PRIVATE vector for xxx_*()
Definition: raytrace.h:253
struct dsp_bb * p
Definition: dsp.c:123
HIDDEN int isect_ray_dsp_bb(struct isect_stuff *isect, struct dsp_bb *dsp_bb)
Definition: dsp.c:2350
size_t bu_cv_w_cookie(void *out, int outcookie, size_t size, void *in, int incookie, size_t count)
Definition: convert.c:474
struct application * ap
Definition: dsp.c:185
void bn_mat_print(const char *title, const mat_t m)
Definition: mat.c:81
#define RT_CK_DB_INTERNAL(_p)
Definition: raytrace.h:207
struct bn_tol * tol
Definition: dsp.c:187
int dmax
Definition: dsp.c:193
#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
fastf_t st_bradius
Radius of BOUNDING sphere.
Definition: raytrace.h:434
void pd_3cont(register FILE *plotfp, double x, double y, double z)
Definition: plot3.c:636
int rt_dsp_get(struct bu_vls *logstr, const struct rt_db_internal *intern, const char *attr)
Definition: dsp.c:4682
const struct bu_structparse rt_dsp_parse[]
Definition: dsp.c:257
#define YMIN
Definition: dsp.c:144
fastf_t crv_c2
curvature in other direction
Definition: raytrace.h:309
int bu_vls_struct_item_named(struct bu_vls *vp, const struct bu_structparse *sdp, const char *name, const char *base, int sep_char)
Definition: parse.c:1203
#define ZTOP
Definition: dsp.c:149
fastf_t a_diverge
slope of beam divergence/mm
Definition: raytrace.h:1600
const struct rt_functab * idb_meth
for ft_ifree(), etc.
Definition: raytrace.h:194
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)
#define V3ARGS(a)
Definition: color.c:56
int a_x
Screen X of ray, if applicable.
Definition: raytrace.h:1596
void rt_dsp_uv(struct application *ap, struct soltab *stp, register struct hit *hitp, register struct uvcoord *uvp)
Definition: dsp.c:3013
const struct bu_structparse rt_dsp_ptab[]
Definition: dsp.c:268
#define NEAR_ZERO(val, epsilon)
Definition: color.c:55
#define BRLCAD_OK
Definition: defines.h:71
uint8_t * ext_buf
Definition: parse.h:216
#define BU_GET(_ptr, _type)
Definition: malloc.h:201
void rt_dsp_curve(register struct curvature *cvp, register struct hit *hitp, struct soltab *stp)
Definition: dsp.c:2992
#define SMOOTHSTEP(x)
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
struct faceuse * nmg_add_loop_to_face(struct shell *s, struct faceuse *fu, struct vertex **verts, int n, int dir)
Definition: nmg_mod.c:1211
struct hit seg_out
OUT information.
Definition: raytrace.h:371
void rt_dsp_ifree(struct rt_db_internal *ip)
Definition: dsp.c:4575
unsigned short dsp_max[3]
Definition: dsp.c:89
#define BBOX_PT(_x, _y, _z)
struct bbox_isect bbox
Definition: dsp.c:189
fastf_t r_max
exit dist from bounding sphere
Definition: raytrace.h:221
size_t bu_vls_strlen(const struct bu_vls *vp)
Definition: vls.c:189
void pdv_3cont(register FILE *plotfp, const fastf_t *pt)
Definition: plot3.c:630
int xsiz
Definition: dsp.c:161
struct dsp_bb * bb_array
Definition: dsp.c:165
int bu_cv_cookie(const char *in)
Definition: convert.c:33
HIDDEN void plot_layers(struct dsp_specific *dsp_sp)
Definition: dsp.c:372
#define UNUSED(parameter)
Definition: common.h:239
HIDDEN int isect_ray_triangle(struct isect_stuff *isect, point_t A, point_t B, point_t C, struct hit *hitp, fastf_t alphabbeta[])
Definition: dsp.c:1319
HIDDEN int recurse_dsp_bb(struct isect_stuff *isect, struct dsp_bb *dsp_bb, point_t minpt, point_t maxpt, point_t bbmin, point_t bbmax)
Definition: dsp.c:2185
int nmg_mark_edges_real(const uint32_t *magic_p)
Definition: nmg_misc.c:846
#define BU_PUT(_ptr, _type)
Definition: malloc.h:215
goto out
Definition: nmg_mod.c:3846
#define DSP_O(m)
Definition: dsp.c:254
void bn_mat_mul(mat_t o, const mat_t a, const mat_t b)
int layers
Definition: dsp.c:163
Support for uniform tolerances.
Definition: tol.h:71
void rt_dsp_make(const struct rt_functab *ftp, struct rt_db_internal *intern)
Definition: dsp.c:4773