BRL-CAD
gqa.c
Go to the documentation of this file.
1 /* G Q A . C
2  * BRL-CAD
3  *
4  * Copyright (c) 2008-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 /** @file libged/gqa.c
21  *
22  * performs a set of quantitative analyses on geometry.
23  *
24  * XXX need to look at gap computation
25  *
26  * plot the points where overlaps start/stop
27  *
28  * Designed to be a framework for 3d sampling of the geometry volume.
29  * TODO: Need to move the sample pattern logic into LIBRT.
30  *
31  */
32 
33 #include "common.h"
34 
35 #include <stdlib.h>
36 #include <string.h>
37 #include <errno.h>
38 #include <sys/types.h>
39 #include <sys/stat.h>
40 #include <math.h>
41 #include <limits.h> /* home of INT_MAX aka MAXINT */
42 
43 
44 #include "bu/parallel.h"
45 #include "bu/getopt.h"
46 #include "vmath.h"
47 #include "raytrace.h"
48 #include "plot3.h"
49 #include "sysv.h"
50 #include "analyze.h"
51 
52 #include "./ged_private.h"
53 
54 
55 /* bu_getopt() options */
56 char *options = "A:a:de:f:g:Gn:N:pP:qrS:s:t:U:u:vV:W:h?";
57 char *options_str = "[-A A|a|b|c|e|g|m|o|p|v|w] [-a az] [-d] [-e el] [-f densityFile] [-g spacing|upper,lower|upper-lower] [-G] [-n nhits] [-N nviews] [-p] [-P ncpus] [-q] [-r] [-S nsamples] [-t overlap_tol] [-U useair] [-u len_units vol_units wt_units] [-v] [-V volume_tol] [-W weight_tol]";
58 
59 #define ANALYSIS_VOLUME 1
60 #define ANALYSIS_WEIGHT 2
61 #define ANALYSIS_OVERLAPS 4
62 #define ANALYSIS_ADJ_AIR 8
63 #define ANALYSIS_GAP 16
64 #define ANALYSIS_EXP_AIR 32 /* exposed air */
65 #define ANALYSIS_BOX 64
66 #define ANALYSIS_INTERFACES 128
67 #define ANALYSIS_CENTROIDS 256
68 #define ANALYSIS_MOMENTS 512
69 #define ANALYSIS_PLOT_OVERLAPS 1024
70 
71 /* Note: struct parsing requires no space after the commas. take care
72  * when formatting this file. if the compile breaks here, it means
73  * that spaces got inserted incorrectly.
74  */
75 #define COMMA ','
76 #define STRCOMMA ","
77 
78 static int analysis_flags;
79 static int multiple_analyses;
80 
81 static double azimuth_deg;
82 static double elevation_deg;
83 static char *densityFileName;
84 static double gridSpacing;
85 static double gridSpacingLimit;
86 static char makeOverlapAssemblies;
87 static size_t require_num_hits;
88 static int ncpu;
89 static double Samples_per_model_axis;
90 static double overlap_tolerance;
91 static double volume_tolerance;
92 static double weight_tolerance;
93 static const fastf_t one_twelfth = 1.0 / 12.0;
94 static int aborted = 0;
95 
96 static int print_per_region_stats;
97 static int max_region_name_len;
98 static int use_air;
99 static int num_objects; /* number of objects specified on command line */
100 static int max_cpus;
101 static int num_views;
102 static int verbose;
103 static int quiet_missed_report;
104 
105 static int plot_files; /* Boolean: Should we produce plot files? */
106 static FILE *plot_weight;
107 static FILE *plot_volume;
108 static FILE *plot_overlaps;
109 static FILE *plot_adjair;
110 static FILE *plot_gaps;
111 static FILE *plot_expair;
112 
113 static int overlap_color[3] = { 255, 255, 0 }; /* yellow */
114 static int gap_color[3] = { 128, 192, 255 }; /* cyan */
115 static int adjAir_color[3] = { 128, 255, 192 }; /* pale green */
116 static int expAir_color[3] = { 255, 128, 255 }; /* magenta */
117 
118 static int debug = 0;
119 #define DLOG if (debug) bu_vls_printf
120 
121 /* Some defines for re-using the values from the application structure
122  * for other purposes
123  */
124 #define A_LENDEN a_color[0]
125 #define A_LEN a_color[1]
126 #define A_STATE a_uptr
127 
128 
129 struct cstate {
130  int curr_view; /* the "view" number we are shooting */
131  int u_axis; /* these 3 are in the range 0..2 inclusive and indicate which axis (X, Y, or Z) */
132  int v_axis; /* is being used for the U, V, or invariant vector direction */
133  int i_axis;
134 
135  /* GED_SEM_WORKER protects this */
136  int v; /* indicates how many "grid_size" steps in the v direction have been taken */
137 
138  /* GED_SEM_STATS protects this */
139  double *m_lenDensity;
140  double *m_len;
141  double *m_volume;
142  double *m_weight;
143  unsigned long *shots;
144  int first; /* this is the first time we've computed a set of views */
145 
146  vect_t u_dir; /* direction of U vector for "current view" */
147  vect_t v_dir; /* direction of V vector for "current view" */
148  struct rt_i *rtip;
149  long steps[3]; /* this is per-dimension, not per-view */
150  vect_t span; /* How much space does the geometry span in each of X, Y, Z directions */
151  vect_t area; /* area of the view for view with invariant at index */
152 
153  fastf_t *m_lenTorque; /* torque vector for each view */
154  fastf_t *m_moi; /* one vector per view for collecting the partial moments of inertia calculation */
155  fastf_t *m_poi; /* one vector per view for collecting the partial products of inertia calculation */
156 
157  struct resource *resp;
158 };
159 
160 
161 struct ged_gqa_plot {
162  struct bn_vlblock *vbp;
163  struct bu_list *vhead;
164 } ged_gqa_plot;
165 
166 /* the entries in the density table */
167 struct density_entry *densities = NULL;
168 static int num_densities;
169 
170 /* summary data structure for objects specified on command line */
171 static struct per_obj_data {
172  char *o_name;
173  double *o_len;
174  double *o_lenDensity;
175  double *o_volume;
176  double *o_weight;
177  fastf_t *o_lenTorque; /* torque vector for each view */
178  fastf_t *o_moi; /* one vector per view for collecting the partial moments of inertia calculation */
179  fastf_t *o_poi; /* one vector per view for collecting the partial products of inertia calculation */
180 } *obj_tbl;
181 
182 /**
183  * this is the data we track for each region
184  */
185 static struct per_region_data {
186  unsigned long hits;
187  double *r_lenDensity; /* for per-region per-view weight computation */
188  double *r_len; /* for per-region, per-view computation */
189  double *r_weight;
190  double *r_volume;
191  struct per_obj_data *optr;
192 } *reg_tbl;
193 
194 
195 /* Access to these lists should be in sections
196  * of code protected by GED_SEM_LIST
197  */
198 
199 /**
200  * list of gaps
201  */
202 static struct region_pair gapList = {
203  {
205  (struct bu_list *)&gapList,
206  (struct bu_list *)&gapList
207  },
208  { "Gaps" },
209  (struct region *)NULL,
210  (unsigned long)0,
211  (double)0.0,
212  {0.0, 0.0, 0.0, }
213 };
214 
215 
216 /**
217  * list of adjacent air
218  */
219 static struct region_pair adjAirList = {
220  {
222  (struct bu_list *)&adjAirList,
223  (struct bu_list *)&adjAirList
224  },
225  { (char *)"Adjacent Air" },
226  (struct region *)NULL,
227  (unsigned long)0,
228  (double)0.0,
229  {0.0, 0.0, 0.0, }
230 };
231 
232 
233 /**
234  * list of exposed air
235  */
236 static struct region_pair exposedAirList = {
237  {
239  (struct bu_list *)&exposedAirList,
240  (struct bu_list *)&exposedAirList
241  },
242  { "Exposed Air" },
243  (struct region *)NULL,
244  (unsigned long)0,
245  (double)0.0,
246  {0.0, 0.0, 0.0, }
247 };
248 
249 
250 /**
251  * list of overlaps
252  */
253 static struct region_pair overlapList = {
254  {
256  (struct bu_list *)&overlapList,
257  (struct bu_list *)&overlapList
258  },
259  { "Overlaps" },
260  (struct region *)NULL,
261  (unsigned long)0,
262  (double)0.0,
263  {0.0, 0.0, 0.0, }
264 };
265 
266 
267 /**
268  * This structure holds the name of a unit value, and the conversion
269  * factor necessary to convert from/to BRL-CAD standard units.
270  *
271  * The standard units are millimeters, cubic millimeters, and grams.
272  *
273  * XXX this section should be extracted to libbu/units.c
274  */
275 struct cvt_tab {
276  double val;
277  char name[32];
278 };
279 
280 
281 static const struct cvt_tab units_tab[3][40] = {
282  {
283  /* length, stolen from bu/units.c with the "none" value
284  * removed Values for converting from given units to mm
285  */
286  {1.0, "mm"}, /* default */
287  /* {0.0, "none"}, */ /* this is removed to force a certain
288  * amount of error checking for the user
289  */
290  {1.0e-7, "angstrom"},
291  {1.0e-7, "decinanometer"},
292  {1.0e-6, "nm"},
293  {1.0e-6, "nanometer"},
294  {1.0e-3, "um"},
295  {1.0e-3, "micrometer"},
296  {1.0e-3, "micron"},
297  {1.0, "millimeter"},
298  {10.0, "cm"},
299  {10.0, "centimeter"},
300  {1000.0, "m"},
301  {1000.0, "meter"},
302  {1000000.0, "km"},
303  {1000000.0, "kilometer"},
304  {25.4, "in"},
305  {25.4, "inch"},
306  {25.4, "inches"}, /* for plural */
307  {304.8, "ft"},
308  {304.8, "foot"},
309  {304.8, "feet"},
310  {456.2, "cubit"},
311  {914.4, "yd"},
312  {914.4, "yard"},
313  {5029.2, "rd"},
314  {5029.2, "rod"},
315  {1609344.0, "mi"},
316  {1609344.0, "mile"},
317  {1852000.0, "nmile"},
318  {1852000.0, "nautical mile"},
319  {1.495979e+14, "AU"},
320  {1.495979e+14, "astronomical unit"},
321  {9.460730e+18, "lightyear"},
322  {3.085678e+19, "pc"},
323  {3.085678e+19, "parsec"},
324  {0.0, ""} /* LAST ENTRY */
325  },
326  {
327  /* volume
328  * Values for converting from given units to mm^3
329  */
330  {1.0, "cu mm"}, /* default */
331 
332  {1.0, "mm"},
333  {1.0, "mm^3"},
334 
335  {1.0e3, "cm"},
336  {1.0e3, "cm^3"},
337  {1.0e3, "cu cm"},
338  {1.0e3, "cc"},
339 
340  {1.0e6, "l"},
341  {1.0e6, "liter"},
342  {1.0e6, "litre"},
343 
344  {1.0e9, "m"},
345  {1.0e9, "m^3"},
346  {1.0e9, "cu m"},
347 
348  {16387.064, "in"},
349  {16387.064, "in^3"},
350  {16387.064, "cu in"},
351 
352  {28316846.592, "ft"},
353 
354  {28316846.592, "ft^3"},
355  {28316846.592, "cu ft"},
356 
357  {764554857.984, "yds"},
358  {764554857.984, "yards"},
359  {764554857.984, "cu yards"},
360 
361  {0.0, ""} /* LAST ENTRY */
362  },
363  {
364  /* weight
365  * Values for converting given units to grams
366  */
367  {1.0, "grams"}, /* default */
368 
369  {1.0, "g"},
370  {0.0648, "gr"},
371  {0.0648, "grains"},
372 
373  {1.0e3, "kg"},
374  {1.0e3, "kilos"},
375  {1.0e3, "kilograms"},
376 
377  {28.35, "oz"},
378  {28.35, "ounce"},
379 
380  {453.6, "lb"},
381  {453.6, "lbs"},
382  {0.0, ""} /* LAST ENTRY */
383  }
384 };
385 
386 
387 /* this table keeps track of the "current" or "user selected units and
388  * the associated conversion values
389  */
390 #define LINE 0
391 #define VOL 1
392 #define WGT 2
393 static const struct cvt_tab *units[3] = {
394  &units_tab[0][0], /* linear */
395  &units_tab[1][0], /* volume */
396  &units_tab[2][0] /* weight */
397 };
398 
399 
400 /**
401  * read_units_double
402  *
403  * Read a non-negative floating point value with optional units
404  *
405  * Return
406  * 1 Failure
407  * 0 Success
408  */
409 int
410 read_units_double(double *val, char *buf, const struct cvt_tab *cvt)
411 {
412  double a;
413  char units_string[256] = {0};
414  int i;
415 
416 
417  i = sscanf(buf, "%lg%256s", &a, units_string);
418 
419  if (i < 0) return 1;
420 
421  if (i == 1) {
422  *val = a;
423 
424  return 0;
425  }
426  if (i == 2) {
427  *val = a;
428  for (; cvt->name[0] != '\0';) {
429  if (!bu_strncmp(cvt->name, units_string, 256)) {
430  goto found_units;
431  } else {
432  cvt++;
433  }
434  }
435  bu_vls_printf(_ged_current_gedp->ged_result_str, "Bad units specifier \"%s\" on value \"%s\"\n", units_string, buf);
436  return 1;
437 
438  found_units:
439  *val = a * cvt->val;
440  return 0;
441  }
442  bu_vls_printf(_ged_current_gedp->ged_result_str, "%s sscanf problem on \"%s\" got %d\n", BU_FLSTR, buf, i);
443  return 1;
444 }
445 
446 
447 /* the above should be extracted to libbu/units.c */
448 
449 
450 /**
451  * Parse through command line flags
452  */
453 static int
454 parse_args(int ac, char *av[])
455 {
456  int c;
457  int i;
458  double a;
459  char *p;
460 
461  /* Turn off getopt's error messages */
462  bu_opterr = 0;
463  bu_optind = 1;
464 
465  /* get all the option flags from the command line */
466  while ((c=bu_getopt(ac, av, options)) != -1) {
467  switch (c) {
468  case 'A':
469  {
470  analysis_flags = 0;
471  multiple_analyses = 0;
472  for (p = bu_optarg; *p; p++) {
473  switch (*p) {
474  case 'A' :
475  analysis_flags = ANALYSIS_VOLUME | ANALYSIS_WEIGHT | \
478  multiple_analyses = 1;
479  break;
480  case 'a' :
481  if (analysis_flags)
482  multiple_analyses = 1;
483 
484  analysis_flags |= ANALYSIS_ADJ_AIR;
485 
486  break;
487  case 'b' :
488  if (analysis_flags)
489  multiple_analyses = 1;
490 
491  analysis_flags |= ANALYSIS_BOX;
492 
493  break;
494  case 'c' :
495  if (analysis_flags)
496  multiple_analyses = 1;
497 
498  analysis_flags |= ANALYSIS_WEIGHT;
499  analysis_flags |= ANALYSIS_CENTROIDS;
500 
501  break;
502  case 'e' :
503  if (analysis_flags)
504  multiple_analyses = 1;
505 
506  analysis_flags |= ANALYSIS_EXP_AIR;
507  break;
508  case 'g' :
509  if (analysis_flags)
510  multiple_analyses = 1;
511 
512  analysis_flags |= ANALYSIS_GAP;
513  break;
514  case 'm' :
515  if (analysis_flags)
516  multiple_analyses = 1;
517 
518  analysis_flags |= ANALYSIS_WEIGHT;
519  analysis_flags |= ANALYSIS_CENTROIDS;
520  analysis_flags |= ANALYSIS_MOMENTS;
521 
522  break;
523  case 'o' :
524  if (analysis_flags)
525  multiple_analyses = 1;
526 
527  analysis_flags |= ANALYSIS_OVERLAPS;
528  break;
529  case 'p' :
530  if (analysis_flags)
531  multiple_analyses = 1;
532 
533  analysis_flags |= ANALYSIS_OVERLAPS;
534  analysis_flags |= ANALYSIS_PLOT_OVERLAPS;
535  break;
536  case 'v' :
537  if (analysis_flags)
538  multiple_analyses = 1;
539 
540  analysis_flags |= ANALYSIS_VOLUME;
541  break;
542  case 'w' :
543  if (analysis_flags)
544  multiple_analyses = 1;
545 
546  analysis_flags |= ANALYSIS_WEIGHT;
547  break;
548  default:
549  bu_vls_printf(_ged_current_gedp->ged_result_str, "Unknown analysis type \"%c\" requested.\n", *p);
550  return -1;
551  }
552  }
553  break;
554  }
555  case 'a':
556  bu_vls_printf(_ged_current_gedp->ged_result_str, "azimuth not implemented\n");
557  if (sscanf(bu_optarg, "%lg", &azimuth_deg) != 1) {
558  bu_vls_printf(_ged_current_gedp->ged_result_str, "error parsing azimuth \"%s\"\n", bu_optarg);
559  return -1;
560  }
561  break;
562  case 'e':
563  bu_vls_printf(_ged_current_gedp->ged_result_str, "elevation not implemented\n");
564  if (sscanf(bu_optarg, "%lg", &elevation_deg) != 1) {
565  bu_vls_printf(_ged_current_gedp->ged_result_str, "error parsing elevation \"%s\"\n", bu_optarg);
566  return -1;
567  }
568  break;
569  case 'd': debug = 1; break;
570 
571  case 'f': densityFileName = bu_optarg; break;
572 
573  case 'g':
574  {
575  double value1, value2;
576  i = 0;
577 
578  /* find out if we have two or one args; user can
579  * separate them with , or - delimiter
580  */
581  p = strchr(bu_optarg, COMMA);
582  if (p)
583  *p++ = '\0';
584  else {
585  p = strchr(bu_optarg, '-');
586  if (p)
587  *p++ = '\0';
588  }
589 
590 
591  if (read_units_double(&value1, bu_optarg, units_tab[0])) {
592  bu_vls_printf(_ged_current_gedp->ged_result_str, "error parsing grid spacing value \"%s\"\n", bu_optarg);
593  return -1;
594  }
595 
596  if (p) {
597  /* we've got 2 values, they are upper limit
598  * and lower limit.
599  */
600  if (read_units_double(&value2, p, units_tab[0])) {
601  bu_vls_printf(_ged_current_gedp->ged_result_str, "error parsing grid spacing limit value \"%s\"\n", p);
602  return -1;
603  }
604 
605  gridSpacing = value1;
606  gridSpacingLimit = value2;
607  } else {
608  gridSpacingLimit = value1;
609 
610  gridSpacing = 0.0; /* flag it */
611  }
612  break;
613  }
614  case 'G':
615  makeOverlapAssemblies = 1;
616  bu_vls_printf(_ged_current_gedp->ged_result_str, "-G option unimplemented\n");
617  return -1;
618  case 'n':
619  if (sscanf(bu_optarg, "%d", &c) != 1 || c < 0) {
620  bu_vls_printf(_ged_current_gedp->ged_result_str, "num_hits must be integer value >= 0, not \"%s\"\n", bu_optarg);
621  return -1;
622  }
623 
624  require_num_hits = (size_t)c;
625  break;
626 
627  case 'N':
628  num_views = atoi(bu_optarg);
629  break;
630  case 'p':
631  plot_files = ! plot_files;
632  break;
633  case 'P':
634  /* cannot ask for more cpu's than the machine has */
635  if ((c=atoi(bu_optarg)) > 0 && c <= max_cpus) ncpu = c;
636  break;
637  case 'q':
638  quiet_missed_report = 1;
639  break;
640  case 'r':
641  print_per_region_stats = 1;
642  break;
643  case 'S':
644  if (sscanf(bu_optarg, "%lg", &a) != 1 || a <= 1.0) {
645  bu_vls_printf(_ged_current_gedp->ged_result_str, "error in specifying minimum samples per model axis: \"%s\"\n", bu_optarg);
646  break;
647  }
648  Samples_per_model_axis = a + 1;
649  break;
650  case 't':
651  if (read_units_double(&overlap_tolerance, bu_optarg, units_tab[0])) {
652  bu_vls_printf(_ged_current_gedp->ged_result_str, "error in overlap tolerance distance \"%s\"\n", bu_optarg);
653  return -1;
654  }
655  break;
656  case 'v':
657  verbose = 1;
658  break;
659  case 'V':
660  if (read_units_double(&volume_tolerance, bu_optarg, units_tab[1])) {
661  bu_vls_printf(_ged_current_gedp->ged_result_str, "error in volume tolerance \"%s\"\n", bu_optarg);
662  return -1;
663  }
664  break;
665  case 'W':
666  if (read_units_double(&weight_tolerance, bu_optarg, units_tab[2])) {
667  bu_vls_printf(_ged_current_gedp->ged_result_str, "error in weight tolerance \"%s\"\n", bu_optarg);
668  return -1;
669  }
670  break;
671 
672  case 'U':
673  use_air = strtol(bu_optarg, (char **)NULL, 10);
674  if (errno == ERANGE || errno == EINVAL) {
675  bu_vls_printf(_ged_current_gedp->ged_result_str, "error in air argument %s\n", bu_optarg);
676  return -1;
677  }
678  break;
679  case 'u':
680  {
681  char *ptr = bu_optarg;
682  const struct cvt_tab *cv;
683  static const char *dim[3] = {"length", "volume", "weight"};
684  char *units_name[3] = {NULL, NULL, NULL};
685  char **units_ap;
686 
687  /* fill in units_name with the names we parse out */
688  units_ap = units_name;
689 
690  /* acquire unit names */
691  for (i = 0; i < 3 && ptr; i++) {
692  int found_unit;
693 
694  if (i == 0) {
695  *units_ap = strtok(ptr, STRCOMMA);
696  } else {
697  *units_ap = strtok(NULL, STRCOMMA);
698  }
699 
700  /* got something? */
701  if (*units_ap == NULL)
702  break;
703 
704  /* got something valid? */
705  found_unit = 0;
706  for (cv = &units_tab[i][0]; cv->name[0] != '\0'; cv++) {
707  if (units_name[i] && BU_STR_EQUAL(cv->name, units_name[i])) {
708  units[i] = cv;
709  found_unit = 1;
710  break;
711  }
712  }
713 
714  if (!found_unit) {
715  bu_vls_printf(_ged_current_gedp->ged_result_str, "Units \"%s\" not found in conversion table\n", units_name[i]);
716  return -1;
717  }
718 
719  ++units_ap;
720  }
721 
723  for (i = 0; i < 3; i++) {
724  bu_vls_printf(_ged_current_gedp->ged_result_str, " %s: %s", dim[i], units[i]->name);
725  }
727  }
728  break;
729 
730  default: /* '?' 'h' */
731  return -1;
732  }
733  }
734 
735  return bu_optind;
736 }
737 
738 
739 /**
740  * Returns
741  * 0 on success
742  * !0 on failure
743  */
744 int
746 {
747  struct stat sb;
748 
749  FILE *fp = (FILE *)NULL;
750  char *buf = NULL;
751  int ret = 0;
752  size_t sret = 0;
753 
754  fp = fopen(name, "rb");
755  if (fp == (FILE *)NULL) {
756  bu_vls_printf(_ged_current_gedp->ged_result_str, "Could not open file - %s\n", name);
757  return GED_ERROR;
758  }
759 
760  if (stat(name, &sb)) {
761  bu_vls_printf(_ged_current_gedp->ged_result_str, "Could not read file - %s\n", name);
762  fclose(fp);
763  return GED_ERROR;
764  }
765 
766  densities = (struct density_entry *)bu_calloc(128, sizeof(struct density_entry), "density entries");
767  num_densities = 128;
768 
769  /* a mapped file would make more sense here */
770  buf = (char *)bu_malloc(sb.st_size+1, "density buffer");
771  sret = fread(buf, sb.st_size, 1, fp);
772  if (sret != 1)
773  perror("fread");
774  ret = parse_densities_buffer(buf, (unsigned long)sb.st_size, densities, _ged_current_gedp->ged_result_str, &num_densities);
775  bu_free(buf, "density buffer");
776  fclose(fp);
777 
778  return ret;
779 }
780 
781 
782 /**
783  * Returns
784  * 0 on success
785  * !0 on failure
786  */
787 int
789 {
790  struct directory *dp;
791  struct rt_db_internal intern;
792  struct rt_binunif_internal *bu;
793  int ret;
794  char *buf;
795 
796  dp = db_lookup(rtip->rti_dbip, "_DENSITIES", LOOKUP_QUIET);
797  if (dp == (struct directory *)NULL) {
798  bu_vls_printf(_ged_current_gedp->ged_result_str, "No \"_DENSITIES\" density table object in database.");
799  bu_vls_printf(_ged_current_gedp->ged_result_str, " If you do not have density data you can still get adjacent air, bounding box, exposed air, gaps, volume or overlaps by using the -Aa, -Ab, -Ae, -Ag, -Av, or -Ao options respectively.\n");
800  return GED_ERROR;
801  }
802 
803  if (rt_db_get_internal(&intern, dp, rtip->rti_dbip, NULL, &rt_uniresource) < 0) {
804  bu_vls_printf(_ged_current_gedp->ged_result_str, "could not import %s\n", dp->d_namep);
805  return GED_ERROR;
806  }
807 
808  if ((intern.idb_major_type & DB5_MAJORTYPE_BINARY_MASK) == 0)
809  return GED_ERROR;
810 
811  bu = (struct rt_binunif_internal *)intern.idb_ptr;
812 
813  RT_CHECK_BINUNIF (bu);
814 
815  densities = (struct density_entry *)bu_calloc(128, sizeof(struct density_entry), "density entries");
816  num_densities = 128;
817 
818  /* Acquire one extra byte to accommodate parse_densities_buffer()
819  * (i.e. it wants to write an EOS in buf[bu->count]).
820  */
821  buf = (char *)bu_malloc(bu->count+1, "density buffer");
822  memcpy(buf, bu->u.int8, bu->count);
823  ret = parse_densities_buffer(buf, bu->count, densities, _ged_current_gedp->ged_result_str, &num_densities);
824  bu_free((void *)buf, "density buffer");
825 
826  return ret;
827 }
828 
829 
830 /**
831  * Write end points of partition to the standard output. If this
832  * routine return !0, this partition will be dropped from the boolean
833  * evaluation.
834  *
835  * Returns:
836  * 0 to eliminate partition with overlap entirely
837  * 1 to retain partition in output list, claimed by reg1
838  * 2 to retain partition in output list, claimed by reg2
839  *
840  * This routine must be prepared to run in parallel
841  */
842 int
843 overlap(struct application *ap,
844  struct partition *pp,
845  struct region *reg1,
846  struct region *reg2,
847  struct partition *hp)
848 {
849 
850  struct xray *rp = &ap->a_ray;
851  struct hit *ihitp = pp->pt_inhit;
852  struct hit *ohitp = pp->pt_outhit;
853  point_t ihit;
854  point_t ohit;
855  double depth;
856 
857  if (!hp) /* unexpected */
858  return 0;
859 
860  /* if one of the regions is air, let it loose */
861  if (reg1->reg_aircode && ! reg2->reg_aircode)
862  return 2;
863  if (reg2->reg_aircode && ! reg1->reg_aircode)
864  return 1;
865 
866  depth = ohitp->hit_dist - ihitp->hit_dist;
867 
868  if (depth < overlap_tolerance)
869  /* too small to matter, pick one or none */
870  return 1;
871 
872  VJOIN1(ihit, rp->r_pt, ihitp->hit_dist, rp->r_dir);
873  VJOIN1(ohit, rp->r_pt, ohitp->hit_dist, rp->r_dir);
874 
875  if (plot_overlaps) {
877  pl_color(plot_overlaps, V3ARGS(overlap_color));
878  pdv_3line(plot_overlaps, ihit, ohit);
880  }
881 
882  if (analysis_flags & ANALYSIS_PLOT_OVERLAPS) {
887  }
888 
889  if (analysis_flags & ANALYSIS_OVERLAPS) {
890  add_unique_pair(&overlapList, reg1, reg2, depth, ihit);
891 
892  if (plot_overlaps) {
894  pl_color(plot_overlaps, V3ARGS(overlap_color));
895  pdv_3line(plot_overlaps, ihit, ohit);
897  }
898  } else {
900  bu_vls_printf(_ged_current_gedp->ged_result_str, "overlap %s %s\n", reg1->reg_name, reg2->reg_name);
902  }
903 
904  /* XXX We should somehow flag the volume/weight calculations as invalid */
905 
906  /* since we have no basis to pick one over the other, just pick */
907  return 1; /* No further consideration to this partition */
908 }
909 
910 
911 /**
912  * Does nothing.
913  */
914 void
916  const struct partition *pp,
917  const struct bu_ptbl *regiontable,
918  const struct partition *InputHdp)
919 {
920  RT_CK_AP(ap);
921  RT_CK_PT(pp);
922  BU_CK_PTBL(regiontable);
923  if (!InputHdp)
924  return;
925 
926  /* do nothing */
927 
928  return;
929 }
930 
931 
932 void exposed_air(struct partition *pp,
933  point_t last_out_point,
934  point_t pt,
935  point_t opt)
936 {
937  /* this shouldn't be air */
938 
939  add_unique_pair(&exposedAirList,
940  pp->pt_regionp,
941  (struct region *)NULL,
942  pp->pt_outhit->hit_dist - pp->pt_inhit->hit_dist, /* thickness */
943  last_out_point); /* location */
944 
945  if (plot_expair) {
947  pl_color(plot_expair, V3ARGS(expAir_color));
948  pdv_3line(plot_expair, pt, opt);
950  }
951 }
952 
953 
954 /**
955  * rt_shootray() was told to call this on a hit. It passes the
956  * application structure which describes the state of the world (see
957  * raytrace.h), and a circular linked list of partitions, each one
958  * describing one in and out segment of one region.
959  *
960  * this routine must be prepared to run in parallel
961  */
962 int
963 hit(struct application *ap, struct partition *PartHeadp, struct seg *segs)
964 {
965  /* see raytrace.h for all of these guys */
966  struct partition *pp;
967  point_t pt, opt, last_out_point;
968  int last_air = 0; /* what was the aircode of the last item */
969  int air_first = 1; /* are we in an air before a solid */
970  double dist; /* the thickness of the partition */
971  double gap_dist;
972  double last_out_dist = -1.0;
973  double val;
974  struct cstate *state = (struct cstate *)ap->A_STATE;
975 
976  if (!segs) /* unexpected */
977  return 0;
978 
979  if (PartHeadp->pt_forw == PartHeadp) return 1;
980 
981 
982  /* examine each partition until we get back to the head */
983  for (pp=PartHeadp->pt_forw; pp != PartHeadp; pp = pp->pt_forw) {
984  struct density_entry *de;
985 
986  /* inhit info */
987  dist = pp->pt_outhit->hit_dist - pp->pt_inhit->hit_dist;
988  VJOIN1(pt, ap->a_ray.r_pt, pp->pt_inhit->hit_dist, ap->a_ray.r_dir);
989  VJOIN1(opt, ap->a_ray.r_pt, pp->pt_outhit->hit_dist, ap->a_ray.r_dir);
990 
991  if (debug) {
994  pp->pt_inhit->hit_dist, pp->pt_outhit->hit_dist);
996  }
997 
998  /* checking for air sticking out of the model. This is done
999  * here because there may be any number of air regions
1000  * sticking out of the model along the ray.
1001  */
1002  if (analysis_flags & ANALYSIS_EXP_AIR) {
1003 
1004  gap_dist = (pp->pt_inhit->hit_dist - last_out_dist);
1005 
1006  /* if air is first on the ray, or we're moving from void/gap to air
1007  * then this is exposed air
1008  */
1009  if (pp->pt_regionp->reg_aircode &&
1010  (air_first || gap_dist > overlap_tolerance)) {
1011  exposed_air(pp, last_out_point, pt, opt);
1012  } else {
1013  air_first = 0;
1014  }
1015  }
1016 
1017  /* looking for voids in the model */
1018  if (analysis_flags & ANALYSIS_GAP) {
1019  if (pp->pt_back != PartHeadp) {
1020  /* if this entry point is further than the previous
1021  * exit point then we have a void
1022  */
1023  gap_dist = pp->pt_inhit->hit_dist - last_out_dist;
1024 
1025  if (gap_dist > overlap_tolerance) {
1026 
1027  /* like overlaps, we only want to report unique pairs */
1028  add_unique_pair(&gapList,
1029  pp->pt_regionp,
1030  pp->pt_back->pt_regionp,
1031  gap_dist,
1032  pt);
1033 
1034  /* like overlaps, let's plot */
1035  if (plot_gaps) {
1036  vect_t gapEnd;
1037  VJOIN1(gapEnd, pt, -gap_dist, ap->a_ray.r_dir);
1038 
1040  pl_color(plot_gaps, V3ARGS(gap_color));
1041  pdv_3line(plot_gaps, pt, gapEnd);
1043  }
1044  }
1045  }
1046  }
1047 
1048  /* computing the weight of the objects */
1049  if (analysis_flags & ANALYSIS_WEIGHT) {
1050  if (debug) {
1052  bu_vls_printf(_ged_current_gedp->ged_result_str, "Hit %s doing weight\n", pp->pt_regionp->reg_name);
1054  }
1055 
1056  /* make sure mater index is within range of densities */
1057  if (pp->pt_regionp->reg_gmater >= num_densities) {
1059  bu_vls_printf(_ged_current_gedp->ged_result_str, "density index %d on region %s is outside of range of table [1..%d]\nSet GIFTmater on region or add entry to density table\n",
1060  pp->pt_regionp->reg_gmater,
1061  pp->pt_regionp->reg_name,
1062  num_densities); /* XXX this should do something else */
1064  return GED_ERROR;
1065  } else {
1066 
1067  /* make sure the density index has been set */
1068  de = &densities[pp->pt_regionp->reg_gmater];
1069  if (de->magic == DENSITY_MAGIC) {
1070  struct per_region_data *prd;
1071  vect_t cmass;
1072  vect_t lenTorque;
1073  fastf_t Lx = state->span[0]/state->steps[0];
1074  fastf_t Ly = state->span[1]/state->steps[1];
1075  fastf_t Lz = state->span[2]/state->steps[2];
1076  fastf_t Lx_sq;
1077  fastf_t Ly_sq;
1078  fastf_t Lz_sq;
1079  fastf_t cell_area;
1080  int los;
1081 
1082  switch (state->i_axis) {
1083  case 0:
1084  Lx_sq = dist*pp->pt_regionp->reg_los*0.01;
1085  Lx_sq *= Lx_sq;
1086  Ly_sq = Ly*Ly;
1087  Lz_sq = Lz*Lz;
1088  cell_area = Ly_sq;
1089  break;
1090  case 1:
1091  Lx_sq = Lx*Lx;
1092  Ly_sq = dist*pp->pt_regionp->reg_los*0.01;
1093  Ly_sq *= Ly_sq;
1094  Lz_sq = Lz*Lz;
1095  cell_area = Lx_sq;
1096  break;
1097  case 2:
1098  default:
1099  Lx_sq = Lx*Lx;
1100  Ly_sq = Ly*Ly;
1101  Lz_sq = dist*pp->pt_regionp->reg_los*0.01;
1102  Lz_sq *= Lz_sq;
1103  cell_area = Lx_sq;
1104  break;
1105  }
1106 
1107  /* factor in the density of this object weight
1108  * computation, factoring in the LOS percentage
1109  * material of the object
1110  */
1111  los = pp->pt_regionp->reg_los;
1112 
1113  if (los < 1) {
1115  bu_vls_printf(_ged_current_gedp->ged_result_str, "bad LOS (%d) on %s\n", los, pp->pt_regionp->reg_name);
1117  }
1118 
1119  /* accumulate the total weight values */
1120  val = de->grams_per_cu_mm * dist * (pp->pt_regionp->reg_los * 0.01);
1121  ap->A_LENDEN += val;
1122 
1123  prd = ((struct per_region_data *)pp->pt_regionp->reg_udata);
1124  /* accumulate the per-region per-view weight values */
1126  prd->r_lenDensity[state->i_axis] += val;
1127 
1128  /* accumulate the per-object per-view weight values */
1129  prd->optr->o_lenDensity[state->i_axis] += val;
1130 
1131  if (analysis_flags & ANALYSIS_CENTROIDS) {
1132  /* calculate the center of mass for this partition */
1133  VJOIN1(cmass, pt, dist*0.5, ap->a_ray.r_dir);
1134 
1135  /* calculate the lenTorque for this partition (i.e. centerOfMass * lenDensity) */
1136  VSCALE(lenTorque, cmass, val);
1137 
1138  /* accumulate per-object per-view torque values */
1139  VADD2(&prd->optr->o_lenTorque[state->i_axis*3], &prd->optr->o_lenTorque[state->i_axis*3], lenTorque);
1140 
1141  /* accumulate the total lenTorque */
1142  VADD2(&state->m_lenTorque[state->i_axis*3], &state->m_lenTorque[state->i_axis*3], lenTorque);
1143 
1144  if (analysis_flags & ANALYSIS_MOMENTS) {
1145  vectp_t moi;
1146  vectp_t poi = &state->m_poi[state->i_axis*3];
1147  fastf_t dx_sq = cmass[X]*cmass[X];
1148  fastf_t dy_sq = cmass[Y]*cmass[Y];
1149  fastf_t dz_sq = cmass[Z]*cmass[Z];
1150  fastf_t mass = val * cell_area;
1151 
1152  /* Collect moments and products of inertia for the current object */
1153  moi = &prd->optr->o_moi[state->i_axis*3];
1154  moi[X] += one_twelfth*mass*(Ly_sq + Lz_sq) + mass*(dy_sq + dz_sq);
1155  moi[Y] += one_twelfth*mass*(Lx_sq + Lz_sq) + mass*(dx_sq + dz_sq);
1156  moi[Z] += one_twelfth*mass*(Lx_sq + Ly_sq) + mass*(dx_sq + dy_sq);
1157  poi = &prd->optr->o_poi[state->i_axis*3];
1158  poi[X] -= mass*cmass[X]*cmass[Y];
1159  poi[Y] -= mass*cmass[X]*cmass[Z];
1160  poi[Z] -= mass*cmass[Y]*cmass[Z];
1161 
1162  /* Collect moments and products of inertia for all objects */
1163  moi = &state->m_moi[state->i_axis*3];
1164  moi[X] += one_twelfth*mass*(Ly_sq + Lz_sq) + mass*(dy_sq + dz_sq);
1165  moi[Y] += one_twelfth*mass*(Lx_sq + Lz_sq) + mass*(dx_sq + dz_sq);
1166  moi[Z] += one_twelfth*mass*(Lx_sq + Ly_sq) + mass*(dx_sq + dy_sq);
1167  poi = &state->m_poi[state->i_axis*3];
1168  poi[X] -= mass*cmass[X]*cmass[Y];
1169  poi[Y] -= mass*cmass[X]*cmass[Z];
1170  poi[Z] -= mass*cmass[Y]*cmass[Z];
1171  }
1172  }
1173 
1175 
1176  } else {
1178  bu_vls_printf(_ged_current_gedp->ged_result_str, "density index %d from region %s is not set.\nAdd entry to density table\n",
1181 
1182  aborted = 1;
1183  return GED_ERROR;
1184  }
1185  }
1186  }
1187 
1188  /* compute the volume of the object */
1189  if (analysis_flags & ANALYSIS_VOLUME) {
1190  struct per_region_data *prd = ((struct per_region_data *)pp->pt_regionp->reg_udata);
1191  ap->A_LEN += dist; /* add to total volume */
1192  {
1194 
1195  /* add to region volume */
1196  prd->r_len[state->curr_view] += dist;
1197 
1198  /* add to object volume */
1199  prd->optr->o_len[state->curr_view] += dist;
1200 
1202  }
1203  if (debug) {
1205  bu_vls_printf(_ged_current_gedp->ged_result_str, "\t\tvol hit %s oDist:%g objVol:%g %s\n",
1206  pp->pt_regionp->reg_name, dist, prd->optr->o_len[state->curr_view], prd->optr->o_name);
1208  }
1209 
1210  if (plot_volume) {
1211  VJOIN1(opt, ap->a_ray.r_pt, pp->pt_outhit->hit_dist, ap->a_ray.r_dir);
1212 
1214  if (ap->a_user & 1) {
1215  pl_color(plot_volume, V3ARGS(gap_color));
1216  } else {
1217  pl_color(plot_volume, V3ARGS(adjAir_color));
1218  }
1219 
1220  pdv_3line(plot_volume, pt, opt);
1222  }
1223  }
1224 
1225 
1226  /* look for two adjacent air regions */
1227  if (analysis_flags & ANALYSIS_ADJ_AIR) {
1228  if (last_air && pp->pt_regionp->reg_aircode &&
1229  pp->pt_regionp->reg_aircode != last_air) {
1230 
1231  double d = pp->pt_outhit->hit_dist - pp->pt_inhit->hit_dist;
1232  point_t aapt;
1233 
1234  add_unique_pair(&adjAirList, pp->pt_back->pt_regionp, pp->pt_regionp, 0.0, pt);
1235 
1236 
1237  d *= 0.25;
1238  VJOIN1(aapt, pt, d, ap->a_ray.r_dir);
1239 
1241  pl_color(plot_adjair, V3ARGS(adjAir_color));
1242  pdv_3line(plot_adjair, pt, aapt);
1244 
1245  }
1246  }
1247 
1248  /* note that this region has been seen */
1249  ((struct per_region_data *)pp->pt_regionp->reg_udata)->hits++;
1250 
1251  last_air = pp->pt_regionp->reg_aircode;
1252  last_out_dist = pp->pt_outhit->hit_dist;
1253  VJOIN1(last_out_point, ap->a_ray.r_pt, pp->pt_outhit->hit_dist, ap->a_ray.r_dir);
1254  }
1255 
1256 
1257  if (analysis_flags & ANALYSIS_EXP_AIR && last_air) {
1258  /* the last thing we hit was air. Make a note of that */
1259  pp = PartHeadp->pt_back;
1260 
1261  exposed_air(pp, last_out_point, pt, opt);
1262  }
1263 
1264 
1265  /* This value is returned by rt_shootray a hit usually returns 1,
1266  * miss 0.
1267  */
1268  return 1;
1269 }
1270 
1271 
1272 /**
1273  * rt_shootray() was told to call this on a miss.
1274  *
1275  * This routine must be prepared to run in parallel
1276  */
1277 int
1278 miss(struct application *ap)
1279 {
1280  RT_CK_APPLICATION(ap);
1281 
1282  return 0;
1283 }
1284 
1285 
1286 /**
1287  * This routine must be prepared to run in parallel
1288  */
1289 int
1290 get_next_row(struct cstate *state)
1291 {
1292  int v;
1293  /* look for more work */
1295 
1296  if (state->v < state->steps[state->v_axis])
1297  v = state->v++; /* get a row to work on */
1298  else
1299  v = 0; /* signal end of work */
1300 
1302 
1303  return v;
1304 }
1305 
1306 
1307 /**
1308  * This routine must be prepared to run in parallel
1309  */
1310 void
1311 plane_worker(int cpu, void *ptr)
1312 {
1313  struct application ap;
1314  int u, v;
1315  double v_coord;
1316  struct cstate *state = (struct cstate *)ptr;
1317  unsigned long shot_cnt;
1318 
1319  if (aborted)
1320  return;
1321 
1322  RT_APPLICATION_INIT(&ap);
1323  ap.a_rt_i = (struct rt_i *)state->rtip; /* application uses this instance */
1324  ap.a_hit = hit; /* where to go on a hit */
1325  ap.a_miss = miss; /* where to go on a miss */
1326  ap.a_logoverlap = logoverlap;
1327  ap.a_overlap = overlap;
1328  ap.a_resource = &state->resp[cpu];
1329  ap.A_LENDEN = 0.0; /* really the cumulative length*density for weight computation*/
1330  ap.A_LEN = 0.0; /* really the cumulative length for volume computation */
1331 
1332  /* gross hack */
1333  ap.a_ray.r_dir[state->u_axis] = ap.a_ray.r_dir[state->v_axis] = 0.0;
1334  ap.a_ray.r_dir[state->i_axis] = 1.0;
1335 
1336  ap.A_STATE = ptr; /* really copying the state ptr to the a_uptr */
1337 
1338  u = -1;
1339 
1340  v = get_next_row(state);
1341 
1342  shot_cnt = 0;
1343  while (v) {
1344 
1345  v_coord = v * gridSpacing;
1346  if (debug) {
1348  bu_vls_printf(_ged_current_gedp->ged_result_str, " v = %d v_coord=%g\n", v, v_coord);
1350  }
1351 
1352  if ((v&1) || state->first) {
1353  /* shoot all the rays in this row. This is either the
1354  * first time a view has been computed or it is an odd
1355  * numbered row in a grid refinement
1356  */
1357  for (u=1; u < state->steps[state->u_axis]; u++) {
1358  ap.a_ray.r_pt[state->u_axis] = ap.a_rt_i->mdl_min[state->u_axis] + u*gridSpacing;
1359  ap.a_ray.r_pt[state->v_axis] = ap.a_rt_i->mdl_min[state->v_axis] + v*gridSpacing;
1360  ap.a_ray.r_pt[state->i_axis] = ap.a_rt_i->mdl_min[state->i_axis];
1361 
1362  if (debug) {
1364  bu_vls_printf(_ged_current_gedp->ged_result_str, "%5g %5g %5g -> %g %g %g\n", V3ARGS(ap.a_ray.r_pt),
1365  V3ARGS(ap.a_ray.r_dir));
1367  }
1368  ap.a_user = v;
1369  (void)rt_shootray(&ap);
1370 
1371  if (aborted)
1372  return;
1373 
1374  shot_cnt++;
1375  }
1376  } else {
1377  /* shoot only the rays we need to on this row. Some of
1378  * them have been computed in a previous iteration.
1379  */
1380  for (u=1; u < state->steps[state->u_axis]; u+=2) {
1381  ap.a_ray.r_pt[state->u_axis] = ap.a_rt_i->mdl_min[state->u_axis] + u*gridSpacing;
1382  ap.a_ray.r_pt[state->v_axis] = ap.a_rt_i->mdl_min[state->v_axis] + v*gridSpacing;
1383  ap.a_ray.r_pt[state->i_axis] = ap.a_rt_i->mdl_min[state->i_axis];
1384 
1385  if (debug) {
1387  bu_vls_printf(_ged_current_gedp->ged_result_str, "%5g %5g %5g -> %g %g %g\n", V3ARGS(ap.a_ray.r_pt),
1388  V3ARGS(ap.a_ray.r_dir));
1390  }
1391  ap.a_user = v;
1392  (void)rt_shootray(&ap);
1393 
1394  if (aborted)
1395  return;
1396 
1397  shot_cnt++;
1398 
1399  if (debug)
1400  if (u+1 < state->steps[state->u_axis]) {
1404  }
1405  }
1406  }
1407 
1408  /* iterate */
1409  v = get_next_row(state);
1410  }
1411 
1412  if (debug && (u == -1)) {
1414  bu_vls_printf(_ged_current_gedp->ged_result_str, "didn't shoot any rays\n");
1416  }
1417 
1418  /* There's nothing else left to work on in this view. It's time
1419  * to add the values we have accumulated to the totals for the
1420  * view and return. When all threads have been through here,
1421  * we'll have returned to serial computation.
1422  */
1424  state->shots[state->curr_view] += shot_cnt;
1425  state->m_lenDensity[state->curr_view] += ap.A_LENDEN; /* add our length*density value */
1426  state->m_len[state->curr_view] += ap.A_LEN; /* add our volume value */
1428 }
1429 
1430 
1431 int
1432 find_cmd_line_obj(struct per_obj_data *obj_rpt, const char *name)
1433 {
1434  int i;
1435  char *str = bu_strdup(name);
1436  char *p;
1437 
1438  p=strchr(str, '/');
1439  if (p) {
1440  *p = '\0';
1441  }
1442 
1443  for (i = 0; i < num_objects; i++) {
1444  if (BU_STR_EQUAL(obj_rpt[i].o_name, str)) {
1445  bu_free(str, "");
1446  return i;
1447  }
1448  }
1449 
1450  bu_vls_printf(_ged_current_gedp->ged_result_str, "%s Didn't find object named \"%s\" in %d entries\n", BU_FLSTR, name, num_objects);
1451 
1452  return GED_ERROR;
1453 }
1454 
1455 
1456 /**
1457  * Allocate data structures for tracking statistics on a per-view
1458  * basis for each of the view, object and region levels.
1459  */
1460 void
1461 allocate_per_region_data(struct cstate *state, int start, int ac, const char *av[])
1462 {
1463  struct region *regp;
1464  struct rt_i *rtip = state->rtip;
1465  int i;
1466  int m;
1467 
1468  if (start > ac) {
1469  /* what? */
1470  bu_log("WARNING: Internal error (start:%d > ac:%d).\n", start, ac);
1471  return;
1472  }
1473 
1474  if (num_objects < 1) {
1475  /* what?? */
1476  bu_log("WARNING: No objects remaining.\n");
1477  return;
1478  }
1479 
1480  if (num_views == 0) {
1481  /* crap. */
1482  bu_log("WARNING: No views specified.\n");
1483  return;
1484  }
1485 
1486  if (rtip->nregions == 0) {
1487  /* dammit! */
1488  bu_log("WARNING: No regions remaining.\n");
1489  return;
1490  }
1491 
1492  state->m_lenDensity = (double *)bu_calloc(num_views, sizeof(double), "densityLen");
1493  state->m_len = (double *)bu_calloc(num_views, sizeof(double), "volume");
1494  state->m_volume = (double *)bu_calloc(num_views, sizeof(double), "volume");
1495  state->m_weight = (double *)bu_calloc(num_views, sizeof(double), "volume");
1496  state->shots = (unsigned long *)bu_calloc(num_views, sizeof(unsigned long), "volume");
1497  state->m_lenTorque = (fastf_t *)bu_calloc(num_views, sizeof(vect_t), "lenTorque");
1498  state->m_moi = (fastf_t *)bu_calloc(num_views, sizeof(vect_t), "moments of inertia");
1499  state->m_poi = (fastf_t *)bu_calloc(num_views, sizeof(vect_t), "products of inertia");
1500 
1501  /* build data structures for the list of objects the user
1502  * specified on the command line
1503  */
1504  obj_tbl = (struct per_obj_data *)bu_calloc(num_objects, sizeof(struct per_obj_data), "report tables");
1505  for (i = 0; i < num_objects; i++) {
1506  obj_tbl[i].o_name = (char *)av[start+i];
1507  obj_tbl[i].o_len = (double *)bu_calloc(num_views, sizeof(double), "o_len");
1508  obj_tbl[i].o_lenDensity = (double *)bu_calloc(num_views, sizeof(double), "o_lenDensity");
1509  obj_tbl[i].o_volume = (double *)bu_calloc(num_views, sizeof(double), "o_volume");
1510  obj_tbl[i].o_weight = (double *)bu_calloc(num_views, sizeof(double), "o_weight");
1511  obj_tbl[i].o_lenTorque = (fastf_t *)bu_calloc(num_views, sizeof(vect_t), "lenTorque");
1512  obj_tbl[i].o_moi = (fastf_t *)bu_calloc(num_views, sizeof(vect_t), "moments of inertia");
1513  obj_tbl[i].o_poi = (fastf_t *)bu_calloc(num_views, sizeof(vect_t), "products of inertia");
1514  }
1515 
1516  /* build objects for each region */
1517  reg_tbl = (struct per_region_data *)bu_calloc(rtip->nregions, sizeof(struct per_region_data), "per_region_data");
1518 
1519 
1520  for (i = 0, BU_LIST_FOR (regp, region, &(rtip->HeadRegion)), i++) {
1521  regp->reg_udata = &reg_tbl[i];
1522 
1523  reg_tbl[i].r_lenDensity = (double *)bu_calloc(num_views, sizeof(double), "r_lenDensity");
1524  reg_tbl[i].r_len = (double *)bu_calloc(num_views, sizeof(double), "r_len");
1525  reg_tbl[i].r_volume = (double *)bu_calloc(num_views, sizeof(double), "len");
1526  reg_tbl[i].r_weight = (double *)bu_calloc(num_views, sizeof(double), "len");
1527 
1528  m = (int)strlen(regp->reg_name);
1529  if (m > max_region_name_len) max_region_name_len = m;
1530  reg_tbl[i].optr = &obj_tbl[ find_cmd_line_obj(obj_tbl, &regp->reg_name[1]) ];
1531 
1532  }
1533 }
1534 
1535 
1536 /**
1537  * list_report
1538  */
1539 void
1541 {
1542  struct region_pair *rp;
1543 
1544  if (BU_LIST_IS_EMPTY(&list->l)) {
1545  bu_vls_printf(_ged_current_gedp->ged_result_str, "No %s\n", (char *)list->r.name);
1546 
1547  return;
1548  }
1549 
1550  bu_vls_printf(_ged_current_gedp->ged_result_str, "list %s:\n", (char *)list->r.name);
1551 
1552  for (BU_LIST_FOR (rp, region_pair, &(list->l))) {
1553  if (rp->r2) {
1554  bu_vls_printf(_ged_current_gedp->ged_result_str, "%s %s count:%lu dist:%g%s @ (%g %g %g)\n",
1555  rp->r.r1->reg_name, rp->r2->reg_name, rp->count,
1556  rp->max_dist / units[LINE]->val, units[LINE]->name, V3ARGS(rp->coord));
1557  } else {
1558  bu_vls_printf(_ged_current_gedp->ged_result_str, "%s count:%lu dist:%g%s @ (%g %g %g)\n",
1559  rp->r.r1->reg_name, rp->count,
1560  rp->max_dist / units[LINE]->val, units[LINE]->name, V3ARGS(rp->coord));
1561  }
1562  }
1563 }
1564 
1565 
1566 /**
1567  * Do some computations prior to raytracing based upon options the
1568  * user has specified
1569  *
1570  * Returns:
1571  * 0 continue, ready to go
1572  * !0 error encountered, terminate processing
1573  */
1574 int
1575 options_prep(struct rt_i *rtip, vect_t span)
1576 {
1577  double newGridSpacing = gridSpacing;
1578  int axis;
1579 
1580  /* figure out where the density values are coming from and get
1581  * them.
1582  */
1583  if (analysis_flags & ANALYSIS_WEIGHT) {
1584  if (densityFileName) {
1585  DLOG(_ged_current_gedp->ged_result_str, "density from file\n");
1586  if (get_densities_from_file(densityFileName) != GED_OK) {
1587  return GED_ERROR;
1588  }
1589  } else {
1590  DLOG(_ged_current_gedp->ged_result_str, "density from db\n");
1591  if (get_densities_from_database(rtip) != GED_OK) {
1592  return GED_ERROR;
1593  }
1594  }
1595  }
1596  /* refine the grid spacing if the user has set a lower bound on
1597  * the number of rays per model axis
1598  */
1599  for (axis=0; axis < 3; axis++) {
1600  if (span[axis] < newGridSpacing*Samples_per_model_axis) {
1601  /* along this axis, the gridSpacing is larger than the
1602  * model span. We need to refine.
1603  */
1604  newGridSpacing = span[axis] / Samples_per_model_axis;
1605  }
1606  }
1607 
1608  if (!ZERO(newGridSpacing - gridSpacing)) {
1609  bu_log("Initial grid spacing %g %s does not allow %g samples per axis.\n",
1610  gridSpacing / units[LINE]->val, units[LINE]->name, Samples_per_model_axis - 1);
1611 
1612  bu_log("Adjusted initial grid spacing to %g %s to get %g samples per model axis.\n",
1613  newGridSpacing / units[LINE]->val, units[LINE]->name, Samples_per_model_axis);
1614 
1615  gridSpacing = newGridSpacing;
1616  }
1617 
1618  /* if the vol/weight tolerances are not set, pick something */
1619  if (analysis_flags & ANALYSIS_VOLUME) {
1620  char *name = "volume.plot3";
1621  if (volume_tolerance < 0.0) {
1622  /* using 1/1000th the volume as a default tolerance, no particular reason */
1623  volume_tolerance = span[X] * span[Y] * span[Z] * 0.001;
1624  bu_log("Using estimated volume tolerance %g %s\n",
1625  volume_tolerance / units[VOL]->val, units[VOL]->name);
1626  } else
1627  bu_log("Using volume tolerance %g %s\n",
1628  volume_tolerance / units[VOL]->val, units[VOL]->name);
1629  if (plot_files)
1630  if ((plot_volume=fopen(name, "wb")) == (FILE *)NULL) {
1631  bu_vls_printf(_ged_current_gedp->ged_result_str, "cannot open plot file %s\n", name);
1632  }
1633  }
1634  if (analysis_flags & ANALYSIS_WEIGHT) {
1635  if (weight_tolerance < 0.0) {
1636  double max_den = 0.0;
1637  int i;
1638  for (i = 0; i < num_densities; i++) {
1639  if (densities[i].grams_per_cu_mm > max_den)
1640  max_den = densities[i].grams_per_cu_mm;
1641  }
1642  weight_tolerance = span[X] * span[Y] * span[Z] * 0.1 * max_den;
1643  bu_vls_printf(_ged_current_gedp->ged_result_str, "setting weight tolerance to %g %s\n",
1644  weight_tolerance / units[WGT]->val,
1645  units[WGT]->name);
1646  } else {
1647  bu_vls_printf(_ged_current_gedp->ged_result_str, "weight tolerance %g\n", weight_tolerance);
1648  }
1649  }
1650  if (analysis_flags & ANALYSIS_GAP) {
1651  char *name = "gaps.plot3";
1652  if (plot_files)
1653  if ((plot_gaps=fopen(name, "wb")) == (FILE *)NULL) {
1654  bu_vls_printf(_ged_current_gedp->ged_result_str, "cannot open plot file %s\n", name);
1655  return GED_ERROR;
1656  }
1657  }
1658  if (analysis_flags & ANALYSIS_OVERLAPS) {
1659  if (!ZERO(overlap_tolerance))
1660  bu_vls_printf(_ged_current_gedp->ged_result_str, "overlap tolerance to %g\n", overlap_tolerance);
1661  if (plot_files) {
1662  char *name = "overlaps.plot3";
1663  if ((plot_overlaps=fopen(name, "wb")) == (FILE *)NULL) {
1664  bu_vls_printf(_ged_current_gedp->ged_result_str, "cannot open plot file %s\n", name);
1665  return GED_ERROR;
1666  }
1667  }
1668  }
1669 
1670  if (print_per_region_stats)
1671  if ((analysis_flags & (ANALYSIS_VOLUME|ANALYSIS_WEIGHT)) == 0)
1672  bu_vls_printf(_ged_current_gedp->ged_result_str, "Note: -r option ignored: neither volume or weight options requested\n");
1673 
1674  if (analysis_flags & ANALYSIS_ADJ_AIR)
1675  if (plot_files) {
1676  char *name = "adj_air.plot3";
1677  if ((plot_adjair=fopen(name, "wb")) == (FILE *)NULL) {
1678  bu_vls_printf(_ged_current_gedp->ged_result_str, "cannot open plot file %s\n", name);
1679  return GED_ERROR;
1680  }
1681  }
1682 
1683  if (analysis_flags & ANALYSIS_EXP_AIR)
1684  if (plot_files) {
1685  char *name = "exp_air.plot3";
1686  if ((plot_expair=fopen(name, "wb")) == (FILE *)NULL) {
1687  bu_vls_printf(_ged_current_gedp->ged_result_str, "cannot open plot file %s\n", name);
1688  return GED_ERROR;
1689  }
1690  }
1691 
1692 
1693  if ((analysis_flags & (ANALYSIS_ADJ_AIR|ANALYSIS_EXP_AIR)) && ! use_air) {
1694  bu_vls_printf(_ged_current_gedp->ged_result_str, "Error: Air regions discarded but air analysis requested!\nSet use_air non-zero or eliminate air analysis\n");
1695  return GED_ERROR;
1696  }
1697 
1698  return GED_OK;
1699 }
1700 
1701 
1702 void
1703 view_reports(struct cstate *state)
1704 {
1705  if (analysis_flags & ANALYSIS_VOLUME) {
1706  int obj;
1707  int view;
1708 
1709  /* for each object, compute the volume for all views */
1710  for (obj = 0; obj < num_objects; obj++) {
1711  double val;
1712  /* compute volume of object for given view */
1713  view = state->curr_view;
1714 
1715  /* compute the per-view volume of this object */
1716 
1717  if (state->shots[view] > 0) {
1718  val = obj_tbl[obj].o_volume[view] =
1719  obj_tbl[obj].o_len[view] * (state->area[view] / state->shots[view]);
1720 
1721  if (verbose)
1722  bu_vls_printf(_ged_current_gedp->ged_result_str, "\t%s volume %g %s\n",
1723  obj_tbl[obj].o_name,
1724  val / units[VOL]->val,
1725  units[VOL]->name);
1726  }
1727  }
1728  }
1729  if (analysis_flags & ANALYSIS_WEIGHT) {
1730  int obj;
1731  int view = state->curr_view;
1732 
1733  for (obj = 0; obj < num_objects; obj++) {
1734  double grams_per_cu_mm = obj_tbl[obj].o_lenDensity[view] *
1735  (state->area[view] / state->shots[view]);
1736 
1737 
1738  if (verbose)
1740  obj_tbl[obj].o_name,
1741  grams_per_cu_mm / units[WGT]->val,
1742  units[WGT]->name);
1743  }
1744  }
1745 }
1746 
1747 
1748 /**
1749  * These checks are unique because they must both be completed. Early
1750  * termination before they are done is not an option. The results
1751  * computed here are used later.
1752  *
1753  * Returns:
1754  * 0 terminate
1755  * 1 continue processing
1756  */
1757 static int
1758 weight_volume_terminate(struct cstate *state)
1759 {
1760  /* Both weight and volume computations rely on this routine to
1761  * compute values that are printed in summaries. Hence, both
1762  * checks must always be done before this routine exits. So we
1763  * store the status (can we terminate processing?) in this
1764  * variable and act on it once both volume and weight computations
1765  * are done.
1766  */
1767  int can_terminate = 1;
1768 
1769  double low, hi, val, delta;
1770 
1771  if (analysis_flags & ANALYSIS_WEIGHT) {
1772  /* for each object, compute the weight for all views */
1773  int obj;
1774 
1775  for (obj = 0; obj < num_objects; obj++) {
1776  int view;
1777  double tmp;
1778 
1779  if (verbose)
1780  bu_vls_printf(_ged_current_gedp->ged_result_str, "object %d\n", obj);
1781 
1782  /* compute weight of object for given view */
1783  low = INFINITY;
1784  hi = -INFINITY;
1785  tmp = 0.0;
1786  for (view = 0; view < num_views; view++) {
1787  val = obj_tbl[obj].o_weight[view] =
1788  obj_tbl[obj].o_lenDensity[view] * (state->area[view] / state->shots[view]);
1789  V_MIN(low, val);
1790  V_MAX(hi, val);
1791  tmp += val;
1792  }
1793  delta = hi - low;
1794 
1795  if (verbose)
1797  "\t%s running avg weight %g %s hi=(%g) low=(%g)\n",
1798  obj_tbl[obj].o_name,
1799  (tmp / num_views) / units[WGT]->val,
1800  units[WGT]->name,
1801  hi / units[WGT]->val,
1802  low / units[WGT]->val);
1803 
1804  if (delta > weight_tolerance) {
1805  /* this object differs too much in each view, so we
1806  * need to refine the grid. signal that we cannot
1807  * terminate.
1808  */
1809  can_terminate = 0;
1810  if (verbose)
1811  bu_vls_printf(_ged_current_gedp->ged_result_str, "\t%s differs too much in weight per view.\n",
1812  obj_tbl[obj].o_name);
1813  }
1814  }
1815  if (can_terminate) {
1816  if (verbose)
1817  bu_vls_printf(_ged_current_gedp->ged_result_str, "all objects within tolerance on weight calculation\n");
1818  }
1819  }
1820 
1821  if (analysis_flags & ANALYSIS_VOLUME) {
1822  /* find the range of values for object volumes */
1823  int obj;
1824 
1825  /* for each object, compute the volume for all views */
1826  for (obj = 0; obj < num_objects; obj++) {
1827  int view;
1828  double tmp;
1829 
1830  /* compute volume of object for given view */
1831  low = INFINITY;
1832  hi = -INFINITY;
1833  tmp = 0.0;
1834  for (view = 0; view < num_views; view++) {
1835  val = obj_tbl[obj].o_volume[view] =
1836  obj_tbl[obj].o_len[view] * (state->area[view] / state->shots[view]);
1837  V_MIN(low, val);
1838  V_MAX(hi, val);
1839  tmp += val;
1840  }
1841  delta = hi - low;
1842 
1843  if (verbose)
1845  "\t%s running avg volume %g %s hi=(%g) low=(%g)\n",
1846  obj_tbl[obj].o_name,
1847  (tmp / num_views) / units[VOL]->val, units[VOL]->name,
1848  hi / units[VOL]->val,
1849  low / units[VOL]->val);
1850 
1851  if (delta > volume_tolerance) {
1852  /* this object differs too much in each view, so we
1853  * need to refine the grid.
1854  */
1855  can_terminate = 0;
1856  if (verbose)
1857  bu_vls_printf(_ged_current_gedp->ged_result_str, "\tvolume tol not met on %s. Refine grid\n",
1858  obj_tbl[obj].o_name);
1859  break;
1860  }
1861  }
1862  }
1863 
1864  if (can_terminate) {
1865  return 0; /* signal we don't want to go onward */
1866  }
1867  return 1;
1868 }
1869 
1870 
1871 /**
1872  * Check to see if we are done processing due to some user specified
1873  * limit being achieved.
1874  *
1875  * Returns:
1876  * 0 Terminate
1877  * 1 Continue processing
1878  */
1879 int
1880 terminate_check(struct cstate *state)
1881 {
1882  int wv_status;
1883  int view;
1884  int obj;
1885 
1886  DLOG(_ged_current_gedp->ged_result_str, "terminate_check\n");
1887  RT_CK_RTI(state->rtip);
1888 
1889  if (plot_overlaps) fflush(plot_overlaps);
1890  if (plot_weight) fflush(plot_weight);
1891  if (plot_volume) fflush(plot_volume);
1892  if (plot_adjair) fflush(plot_adjair);
1893  if (plot_gaps) fflush(plot_gaps);
1894  if (plot_expair) fflush(plot_expair);
1895 
1896  /* this computation is done first, because there are side effects
1897  * that must be obtained whether we terminate or not
1898  */
1899  wv_status = weight_volume_terminate(state);
1900 
1901 
1902  /* if we've reached the grid limit, we're done, no matter what */
1903  if (gridSpacing < gridSpacingLimit) {
1904  bu_vls_printf(_ged_current_gedp->ged_result_str, "NOTE: Stopped, grid spacing refined to %g (below lower limit %g).\n",
1905  gridSpacing, gridSpacingLimit);
1906  return 0;
1907  }
1908 
1909  /* if we are doing one of the "Error" checking operations:
1910  * Overlap, gap, adj_air, exp_air, then we ALWAYS go to the grid
1911  * spacing limit and we ALWAYS terminate on first error/list-entry
1912  */
1913  if ((analysis_flags & ANALYSIS_OVERLAPS)) {
1914  if (BU_LIST_NON_EMPTY(&overlapList.l)) {
1915  /* since we've found an overlap, we can quit */
1916  return 0;
1917  } else {
1918  bu_vls_printf(_ged_current_gedp->ged_result_str, "overlaps list at %gmm is empty\n", gridSpacing);
1919  }
1920  }
1921  if ((analysis_flags & ANALYSIS_GAP)) {
1922  if (BU_LIST_NON_EMPTY(&gapList.l)) {
1923  /* since we've found a gap, we can quit */
1924  return 0;
1925  }
1926  }
1927  if ((analysis_flags & ANALYSIS_ADJ_AIR)) {
1928  if (BU_LIST_NON_EMPTY(&adjAirList.l)) {
1929  /* since we've found adjacent air, we can quit */
1930  return 0;
1931  }
1932  }
1933  if ((analysis_flags & ANALYSIS_EXP_AIR)) {
1934  if (BU_LIST_NON_EMPTY(&exposedAirList.l)) {
1935  /* since we've found exposed air, we can quit */
1936  return 0;
1937  }
1938  }
1939 
1940 
1941  if (analysis_flags & (ANALYSIS_WEIGHT|ANALYSIS_VOLUME)) {
1942  /* volume/weight checks only get to terminate processing if
1943  * there are no "error" check computations being done
1944  */
1945  if (analysis_flags & (ANALYSIS_GAP|ANALYSIS_ADJ_AIR|ANALYSIS_OVERLAPS|ANALYSIS_EXP_AIR)) {
1946  if (verbose)
1947  bu_vls_printf(_ged_current_gedp->ged_result_str, "Volume/Weight tolerance met. Cannot terminate calculation due to error computations\n");
1948  } else {
1949  struct region *regp;
1950  int all_hit = 1;
1951  size_t hits;
1952 
1953  if (require_num_hits > 0) {
1954  /* check to make sure every region was hit at least once */
1955  for (BU_LIST_FOR (regp, region, &(state->rtip->HeadRegion))) {
1956  RT_CK_REGION(regp);
1957 
1958  hits = (size_t)((struct per_region_data *)regp->reg_udata)->hits;
1959  if (hits < require_num_hits) {
1960  all_hit = 0;
1961  if (verbose) {
1962  if (hits == 0 && !quiet_missed_report) {
1963  bu_vls_printf(_ged_current_gedp->ged_result_str, "%s was not hit\n", regp->reg_name);
1964  } else if (hits) {
1965  bu_vls_printf(_ged_current_gedp->ged_result_str, "%s hit only %zu times (< %zu)\n",
1966  regp->reg_name, hits, require_num_hits);
1967  }
1968  }
1969  }
1970  }
1971 
1972  if (all_hit && wv_status == 0) {
1973  if (verbose)
1974  bu_vls_printf(_ged_current_gedp->ged_result_str, "%s: Volume/Weight tolerance met. Terminate\n", BU_FLSTR);
1975  return 0; /* terminate */
1976  }
1977  } else {
1978  if (wv_status == 0) {
1979  if (verbose)
1980  bu_vls_printf(_ged_current_gedp->ged_result_str, "%s: Volume/Weight tolerance met. Terminate\n", BU_FLSTR);
1981  return 0; /* terminate */
1982  }
1983  }
1984  }
1985  }
1986 
1987  for (view=0; view < num_views; view++) {
1988  for (obj = 0; obj < num_objects; obj++) {
1989  VSCALE(&obj_tbl[obj].o_moi[view*3], &obj_tbl[obj].o_moi[view*3], 0.25);
1990  VSCALE(&obj_tbl[obj].o_poi[view*3], &obj_tbl[obj].o_poi[view*3], 0.25);
1991  }
1992 
1993  VSCALE(&state->m_moi[view*3], &state->m_moi[view*3], 0.25);
1994  VSCALE(&state->m_poi[view*3], &state->m_poi[view*3], 0.25);
1995  }
1996 
1997  return 1;
1998 }
1999 
2000 
2001 /**
2002  * summary_reports
2003  */
2004 void
2005 summary_reports(struct cstate *state)
2006 {
2007  int view;
2008  int obj;
2009  double avg_mass;
2010  struct region *regp;
2011 
2012  if (multiple_analyses)
2013  bu_vls_printf(_ged_current_gedp->ged_result_str, "Summaries (%gmm grid spacing):\n", gridSpacing*2);
2014  else
2015  bu_vls_printf(_ged_current_gedp->ged_result_str, "Summary (%gmm grid spacing):\n", gridSpacing*2);
2016 
2017  if (analysis_flags & ANALYSIS_WEIGHT) {
2019  for (obj = 0; obj < num_objects; obj++) {
2020  avg_mass = 0.0;
2021 
2022  for (view=0; view < num_views; view++) {
2023  /* computed in terminate_check() */
2024  avg_mass += obj_tbl[obj].o_weight[view];
2025  }
2026  avg_mass /= num_views;
2027  bu_vls_printf(_ged_current_gedp->ged_result_str, "\t%*s %g %s\n", -max_region_name_len, obj_tbl[obj].o_name,
2028  avg_mass / units[WGT]->val, units[WGT]->name);
2029 
2030  if (analysis_flags & ANALYSIS_CENTROIDS &&
2031  !ZERO(avg_mass)) {
2032  vect_t centroid;
2033  fastf_t Dx_sq, Dy_sq, Dz_sq;
2034  fastf_t inv_total_mass = 1.0/avg_mass;
2035 
2036  VSETALL(centroid, 0.0);
2037  for (view=0; view < num_views; view++) {
2038  vect_t torque;
2039  fastf_t cell_area = state->area[view] / state->shots[view];
2040 
2041  VSCALE(torque, &obj_tbl[obj].o_lenTorque[view*3], cell_area);
2042  VADD2(centroid, centroid, torque);
2043  }
2044 
2045  VSCALE(centroid, centroid, 1.0/(fastf_t)num_views);
2046  VSCALE(centroid, centroid, inv_total_mass);
2048  "\t\tcentroid: (%g %g %g) mm\n", V3ARGS(centroid));
2049 
2050  /* Do the final calculations for the moments of
2051  * inertia for the current object.
2052  */
2053  if (analysis_flags & ANALYSIS_MOMENTS) {
2054  struct bu_vls title = BU_VLS_INIT_ZERO;
2055  mat_t tmat; /* total mat */
2056 
2057  MAT_ZERO(tmat);
2058  for (view=0; view < num_views; view++) {
2059  vectp_t moi = &obj_tbl[obj].o_moi[view*3];
2060  vectp_t poi = &obj_tbl[obj].o_poi[view*3];
2061 
2062  tmat[MSX] += moi[X];
2063  tmat[MSY] += moi[Y];
2064  tmat[MSZ] += moi[Z];
2065  tmat[1] += poi[X];
2066  tmat[2] += poi[Y];
2067  tmat[6] += poi[Z];
2068  }
2069 
2070  tmat[MSX] /= (fastf_t)num_views;
2071  tmat[MSY] /= (fastf_t)num_views;
2072  tmat[MSZ] /= (fastf_t)num_views;
2073  tmat[1] /= (fastf_t)num_views;
2074  tmat[2] /= (fastf_t)num_views;
2075  tmat[6] /= (fastf_t)num_views;
2076 
2077  /* Lastly, apply the parallel axis theorem */
2078  Dx_sq = centroid[X]*centroid[X];
2079  Dy_sq = centroid[Y]*centroid[Y];
2080  Dz_sq = centroid[Z]*centroid[Z];
2081  tmat[MSX] -= avg_mass*(Dy_sq + Dz_sq);
2082  tmat[MSY] -= avg_mass*(Dx_sq + Dz_sq);
2083  tmat[MSZ] -= avg_mass*(Dx_sq + Dy_sq);
2084  tmat[1] += avg_mass*centroid[X]*centroid[Y];
2085  tmat[2] += avg_mass*centroid[X]*centroid[Z];
2086  tmat[6] += avg_mass*centroid[Y]*centroid[Z];
2087 
2088  tmat[4] = tmat[1];
2089  tmat[8] = tmat[2];
2090  tmat[9] = tmat[6];
2091 
2092  bu_vls_printf(&title, "For the Moments and Products of Inertia For %s", obj_tbl[obj].o_name);
2094  bu_vls_free(&title);
2095  }
2096  }
2097  }
2098 
2099 
2100  if (print_per_region_stats) {
2101  double *wv;
2103  for (BU_LIST_FOR (regp, region, &(state->rtip->HeadRegion))) {
2104  double low = INFINITY;
2105  double hi = -INFINITY;
2106 
2107  avg_mass = 0.0;
2108 
2109  for (view=0; view < num_views; view++) {
2110  wv = &((struct per_region_data *)regp->reg_udata)->r_weight[view];
2111 
2112  *wv = ((struct per_region_data *)regp->reg_udata)->r_lenDensity[view] *
2113  (state->area[view]/state->shots[view]);
2114 
2115  *wv /= units[WGT]->val;
2116 
2117  avg_mass += *wv;
2118 
2119  if (*wv < low) low = *wv;
2120  if (*wv > hi) hi = *wv;
2121  }
2122 
2123  avg_mass /= num_views;
2124  bu_vls_printf(_ged_current_gedp->ged_result_str, "\t%s %g %s +(%g) -(%g)\n",
2125  regp->reg_name,
2126  avg_mass,
2127  units[WGT]->name,
2128  hi - avg_mass,
2129  avg_mass - low);
2130  }
2131  }
2132 
2133  /* print grand totals */
2134  avg_mass = 0.0;
2135  for (view=0; view < num_views; view++) {
2136  avg_mass += state->m_weight[view] =
2137  state->m_lenDensity[view] *
2138  (state->area[view] / state->shots[view]);
2139  }
2140 
2141  avg_mass /= num_views;
2142  bu_vls_printf(_ged_current_gedp->ged_result_str, " Average total weight: %g %s\n", avg_mass / units[WGT]->val, units[WGT]->name);
2143 
2144  if (analysis_flags & ANALYSIS_CENTROIDS &&
2145  !ZERO(avg_mass)) {
2146  vect_t centroid;
2147  fastf_t Dx_sq, Dy_sq, Dz_sq;
2148  fastf_t inv_total_mass = 1.0/avg_mass;
2149 
2150  VSETALL(centroid, 0.0);
2151  for (view=0; view < num_views; view++) {
2152  vect_t torque;
2153  fastf_t cell_area = state->area[view] / state->shots[view];
2154 
2155  VSCALE(torque, &state->m_lenTorque[view*3], cell_area);
2156  VADD2(centroid, centroid, torque);
2157  }
2158 
2159  VSCALE(centroid, centroid, 1.0/(fastf_t)num_views);
2160  VSCALE(centroid, centroid, inv_total_mass);
2162  " Average centroid: (%g %g %g) mm\n", V3ARGS(centroid));
2163 
2164  /* Do the final calculations for the moments of inertia
2165  * for the current object.
2166  */
2167  if (analysis_flags & ANALYSIS_MOMENTS) {
2168  mat_t tmat; /* total mat */
2169 
2170  MAT_ZERO(tmat);
2171  for (view=0; view < num_views; view++) {
2172  vectp_t moi = &state->m_moi[view*3];
2173  vectp_t poi = &state->m_poi[view*3];
2174 
2175  tmat[MSX] += moi[X];
2176  tmat[MSY] += moi[Y];
2177  tmat[MSZ] += moi[Z];
2178  tmat[1] += poi[X];
2179  tmat[2] += poi[Y];
2180  tmat[6] += poi[Z];
2181  }
2182 
2183  tmat[MSX] /= (fastf_t)num_views;
2184  tmat[MSY] /= (fastf_t)num_views;
2185  tmat[MSZ] /= (fastf_t)num_views;
2186  tmat[1] /= (fastf_t)num_views;
2187  tmat[2] /= (fastf_t)num_views;
2188  tmat[6] /= (fastf_t)num_views;
2189 
2190  /* Lastly, apply the parallel axis theorem */
2191  Dx_sq = centroid[X]*centroid[X];
2192  Dy_sq = centroid[Y]*centroid[Y];
2193  Dz_sq = centroid[Z]*centroid[Z];
2194  tmat[MSX] -= avg_mass*(Dy_sq + Dz_sq);
2195  tmat[MSY] -= avg_mass*(Dx_sq + Dz_sq);
2196  tmat[MSZ] -= avg_mass*(Dx_sq + Dy_sq);
2197  tmat[1] += avg_mass*centroid[X]*centroid[Y];
2198  tmat[2] += avg_mass*centroid[X]*centroid[Z];
2199  tmat[6] += avg_mass*centroid[Y]*centroid[Z];
2200 
2201  tmat[4] = tmat[1];
2202  tmat[8] = tmat[2];
2203  tmat[9] = tmat[6];
2204 
2205  bn_mat_print_vls("For the Moments and Products of Inertia For\n\tAll Specified Objects",
2207  }
2208  }
2209  }
2210 
2211 
2212  if (analysis_flags & ANALYSIS_VOLUME) {
2214 
2215  /* print per-object */
2216  for (obj = 0; obj < num_objects; obj++) {
2217  avg_mass = 0.0;
2218 
2219  for (view=0; view < num_views; view++)
2220  avg_mass += obj_tbl[obj].o_volume[view];
2221 
2222  avg_mass /= num_views;
2223  bu_vls_printf(_ged_current_gedp->ged_result_str, "\t%*s %g %s\n", -max_region_name_len, obj_tbl[obj].o_name,
2224  avg_mass / units[VOL]->val, units[VOL]->name);
2225  }
2226 
2227  if (print_per_region_stats) {
2228  double *vv;
2229 
2231  for (BU_LIST_FOR (regp, region, &(state->rtip->HeadRegion))) {
2232  double low = INFINITY;
2233  double hi = -INFINITY;
2234  avg_mass = 0.0;
2235 
2236  for (view=0; view < num_views; view++) {
2237  vv = &((struct per_region_data *)regp->reg_udata)->r_volume[view];
2238 
2239  /* convert view length to a volume */
2240  *vv = ((struct per_region_data *)regp->reg_udata)->r_len[view] *
2241  (state->area[view] / state->shots[view]);
2242 
2243  /* convert to user's units */
2244  *vv /= units[VOL]->val;
2245 
2246  /* find limits of values */
2247  if (*vv < low) low = *vv;
2248  if (*vv > hi) hi = *vv;
2249 
2250  avg_mass += *vv;
2251  }
2252 
2253  avg_mass /= num_views;
2254 
2255  bu_vls_printf(_ged_current_gedp->ged_result_str, "\t%s volume:%g %s +(%g) -(%g)\n",
2256  regp->reg_name,
2257  avg_mass,
2258  units[VOL]->name,
2259  hi - avg_mass,
2260  avg_mass - low);
2261  }
2262  }
2263 
2264 
2265  /* print grand totals */
2266  avg_mass = 0.0;
2267  for (view=0; view < num_views; view++) {
2268  avg_mass += state->m_volume[view] =
2269  state->m_len[view] * (state->area[view] / state->shots[view]);
2270  }
2271 
2272  avg_mass /= num_views;
2273  bu_vls_printf(_ged_current_gedp->ged_result_str, " Average total volume: %g %s\n", avg_mass / units[VOL]->val, units[VOL]->name);
2274  }
2275  if (analysis_flags & ANALYSIS_OVERLAPS) list_report(&overlapList);
2276  if (analysis_flags & ANALYSIS_ADJ_AIR) list_report(&adjAirList);
2277  if (analysis_flags & ANALYSIS_GAP) list_report(&gapList);
2278  if (analysis_flags & ANALYSIS_EXP_AIR) list_report(&exposedAirList);
2279 
2280  for (BU_LIST_FOR (regp, region, &(state->rtip->HeadRegion))) {
2281  size_t hits;
2282  struct region_pair *rp;
2283  int is_overlap_only_hit;
2284 
2285  RT_CK_REGION(regp);
2286  hits = (size_t)((struct per_region_data *)regp->reg_udata)->hits;
2287  if (hits < require_num_hits) {
2288  if (hits == 0 && !quiet_missed_report) {
2289  is_overlap_only_hit = 0;
2290  if (analysis_flags & ANALYSIS_OVERLAPS) {
2291  /* If the region is in the overlap list, it has
2292  * been hit even though the hit count is zero.
2293  * Do not report zero hit regions if they are in
2294  * the overlap list.
2295  */
2296  for (BU_LIST_FOR (rp, region_pair, &(overlapList.l))) {
2297  if (rp->r.r1->reg_name == regp->reg_name) {
2298  is_overlap_only_hit = 1;
2299  break;
2300  } else if (rp->r2) {
2301  if (rp->r2->reg_name == regp->reg_name) {
2302  is_overlap_only_hit = 1;
2303  break;
2304  }
2305  }
2306  }
2307  }
2308  if (!is_overlap_only_hit) {
2309  bu_vls_printf(_ged_current_gedp->ged_result_str, "%s was not hit\n", regp->reg_name);
2310  }
2311  } else if (hits) {
2312  bu_vls_printf(_ged_current_gedp->ged_result_str, "%s hit only %zu times (< %zu)\n",
2313  regp->reg_name, hits, require_num_hits);
2314  }
2315  }
2316  }
2317 }
2318 
2319 
2320 int
2321 ged_gqa(struct ged *gedp, int argc, const char *argv[])
2322 {
2323  int arg_count;
2324  struct rt_i *rtip;
2325  int i;
2326  struct cstate state;
2327  int start_objs; /* index in command line args where geom object list starts */
2328  struct region_pair *rp;
2329  struct region *regp;
2330  static const char *usage = "object [object ...]";
2331  struct resource resp[MAX_PSW]; /* memory resources for multi-cpu processing */
2332 
2334  GED_CHECK_ARGC_GT_0(gedp, argc, GED_ERROR);
2335 
2336  /* initialize result */
2337  bu_vls_trunc(gedp->ged_result_str, 0);
2338 
2339  /* must be wanting help */
2340  if (argc == 1) {
2341  bu_vls_printf(gedp->ged_result_str, "Usage: %s %s %s", argv[0], options_str, usage);
2342  return GED_HELP;
2343  }
2344 
2345  _ged_current_gedp = gedp;
2346 
2347  analysis_flags = ANALYSIS_VOLUME | ANALYSIS_OVERLAPS | ANALYSIS_WEIGHT |
2349  multiple_analyses = 1;
2350  azimuth_deg = 0.0;
2351  elevation_deg = 0.0;
2352  densityFileName = (char *)0;
2353 
2354  /* FIXME: this is completely arbitrary, should probably be based
2355  * on the model size.
2356  */
2357  gridSpacing = 50.0;
2358 
2359  /* default grid spacing limit is based on the current distance
2360  * tolerance, one order of magnitude greater.
2361  *
2362  * FIXME: should probably be based on the model size.
2363  */
2364  gridSpacingLimit = 10.0 * gedp->ged_wdbp->wdb_tol.dist;
2365 
2366  makeOverlapAssemblies = 0;
2367  require_num_hits = 1;
2368  max_cpus = ncpu = bu_avail_cpus();
2369  Samples_per_model_axis = 2.0;
2370  overlap_tolerance = 0.0;
2371  volume_tolerance = -1.0;
2372  weight_tolerance = -1.0;
2373  print_per_region_stats = 0;
2374  max_region_name_len = 0;
2375  use_air = 1;
2376  num_objects = 0;
2377  num_views = 3;
2378  verbose = 0;
2379  quiet_missed_report = 0;
2380  plot_files = 0;
2381  plot_weight = (FILE *)0;
2382  plot_volume = (FILE *)0;
2383  plot_overlaps = (FILE *)0;
2384  plot_adjair = (FILE *)0;
2385  plot_gaps = (FILE *)0;
2386  plot_expair = (FILE *)0;
2387  debug = 0;
2388 
2389  /* parse command line arguments */
2390  arg_count = parse_args(argc, (char **)argv);
2391 
2392  if (arg_count < 0 || (argc-arg_count) < 1) {
2393  bu_vls_printf(gedp->ged_result_str, "Usage: %s %s %s", argv[0], options_str, usage);
2394  return GED_ERROR;
2395  }
2396 
2397  if (analysis_flags & ANALYSIS_PLOT_OVERLAPS) {
2399  ged_gqa_plot.vhead = rt_vlblock_find(ged_gqa_plot.vbp, 0xFF, 0xFF, 0x00);
2400  }
2401 
2402  rtip = rt_new_rti(gedp->ged_wdbp->dbip);
2403  rtip->useair = use_air;
2404 
2405  start_objs = arg_count;
2406  num_objects = argc - arg_count;
2407 
2408  /* Initialize all the per-CPU memory resources. The number of
2409  * processors can change at runtime, init them all.
2410  */
2411  memset(resp, 0, sizeof(resp));
2412  for (i = 0; i < MAX_PSW; i++) {
2413  rt_init_resource(&resp[i], i, rtip);
2414  }
2415  state.resp = resp;
2416 
2417  /* Walk trees. Here we identify any object trees in the database
2418  * that the user wants included in the ray trace.
2419  */
2420  for (; arg_count < argc; arg_count++) {
2421  if (rt_gettree(rtip, argv[arg_count]) < 0) {
2422  fprintf(stderr, "rt_gettree(%s) FAILED\n", argv[arg_count]);
2423  return GED_ERROR;
2424  }
2425  }
2426 
2427  /* This gets the database ready for ray tracing. (it precomputes
2428  * some values, sets up space partitioning, etc.)
2429  */
2430  rt_prep_parallel(rtip, ncpu);
2431 
2432  /* we now have to subdivide space */
2433  VSUB2(state.span, rtip->mdl_max, rtip->mdl_min);
2434  state.area[0] = state.span[1] * state.span[2];
2435  state.area[1] = state.span[2] * state.span[0];
2436  state.area[2] = state.span[0] * state.span[1];
2437 
2438  if (analysis_flags & ANALYSIS_BOX) {
2439  bu_vls_printf(gedp->ged_result_str, "bounding box: %g %g %g %g %g %g\n",
2440  V3ARGS(rtip->mdl_min), V3ARGS(rtip->mdl_max));
2441 
2442  bu_vls_printf(gedp->ged_result_str, "Area: (%g, %g, %g)\n", state.area[X], state.area[Y], state.area[Z]);
2443  }
2444  if (verbose) bu_vls_printf(gedp->ged_result_str, "ncpu: %d\n", ncpu);
2445 
2446  /* if the user did not specify the initial grid spacing limit, we
2447  * need to compute a reasonable one for them.
2448  */
2449  if (ZERO(gridSpacing)) {
2450  double min_span = MAX_FASTF;
2451  VPRINT("span", state.span);
2452 
2453  V_MIN(min_span, state.span[X]);
2454  V_MIN(min_span, state.span[Y]);
2455  V_MIN(min_span, state.span[Z]);
2456 
2457  gridSpacing = gridSpacingLimit;
2458  do {
2459  gridSpacing *= 2.0;
2460  } while (gridSpacing < min_span);
2461 
2462  gridSpacing *= 0.25;
2463  if (gridSpacing < gridSpacingLimit) gridSpacing = gridSpacingLimit;
2464 
2465  bu_log("Trying estimated initial grid spacing: %g %s\n",
2466  gridSpacing / units[LINE]->val, units[LINE]->name);
2467  } else {
2468  bu_log("Trying initial grid spacing: %g %s\n",
2469  gridSpacing / units[LINE]->val, units[LINE]->name);
2470  }
2471 
2472  bu_log("Using grid spacing lower limit: %g %s\n",
2473  gridSpacingLimit / units[LINE]->val, units[LINE]->name);
2474 
2475  if (options_prep(rtip, state.span) != GED_OK) return GED_ERROR;
2476 
2477  /* initialize some stuff */
2478  state.rtip = rtip;
2479  state.first = 1;
2480  allocate_per_region_data(&state, start_objs, argc, argv);
2481 
2482  /* compute */
2483  do {
2484  double inv_spacing = 1.0/gridSpacing;
2485  int view;
2486 
2487  VSCALE(state.steps, state.span, inv_spacing);
2488 
2489  bu_log("Processing with grid spacing %g %s %ld x %ld x %ld\n",
2490  gridSpacing / units[LINE]->val,
2491  units[LINE]->name,
2492  state.steps[0]-1,
2493  state.steps[1]-1,
2494  state.steps[2]-1);
2495 
2496 
2497  for (view=0; view < num_views; view++) {
2498 
2499  if (verbose)
2500  bu_vls_printf(gedp->ged_result_str, " view %d\n", view);
2501 
2502  /* gross hack. By assuming we have <= 3 views, we can let
2503  * the view # indicate a coordinate axis. Note this is
2504  * used as an index into state.area[]
2505  */
2506  state.i_axis = state.curr_view = view;
2507  state.u_axis = (state.curr_view+1) % 3;
2508  state.v_axis = (state.curr_view+2) % 3;
2509 
2510  state.u_dir[state.u_axis] = 1;
2511  state.u_dir[state.v_axis] = 0;
2512  state.u_dir[state.i_axis] = 0;
2513 
2514  state.v_dir[state.u_axis] = 0;
2515  state.v_dir[state.v_axis] = 1;
2516  state.v_dir[state.i_axis] = 0;
2517  state.v = 1;
2518 
2519  bu_parallel(plane_worker, ncpu, (void *)&state);
2520 
2521  if (aborted)
2522  goto aborted;
2523 
2524  view_reports(&state);
2525  }
2526 
2527  state.first = 0;
2528  gridSpacing *= 0.5;
2529 
2530  } while (terminate_check(&state));
2531 
2532 aborted:
2533  if (plot_overlaps) fclose(plot_overlaps);
2534  if (plot_weight) fclose(plot_weight);
2535  if (plot_volume) fclose(plot_volume);
2536  if (plot_adjair) fclose(plot_adjair);
2537  if (plot_gaps) fclose(plot_gaps);
2538  if (plot_expair) fclose(plot_expair);
2539 
2540 
2541  if (verbose)
2542  bu_vls_printf(gedp->ged_result_str, "Computation Done\n");
2543 
2544  if (!aborted) {
2545  summary_reports(&state);
2546 
2547  if (analysis_flags & ANALYSIS_PLOT_OVERLAPS)
2548  _ged_cvt_vlblock_to_solids(gedp, ged_gqa_plot.vbp, "OVERLAPS", 0);
2549  } else
2550  aborted = 0; /* reset flag */
2551 
2552  if (analysis_flags & ANALYSIS_PLOT_OVERLAPS)
2554 
2555  /* Clear out the lists */
2556  while (BU_LIST_WHILE (rp, region_pair, &overlapList.l)) {
2557  BU_LIST_DEQUEUE(&rp->l);
2558  bu_free(rp, "overlapList items");
2559  }
2560  while (BU_LIST_WHILE (rp, region_pair, &adjAirList.l)) {
2561  BU_LIST_DEQUEUE(&rp->l);
2562  bu_free(rp, "adjAirList items");
2563  }
2564  while (BU_LIST_WHILE (rp, region_pair, &gapList.l)) {
2565  BU_LIST_DEQUEUE(&rp->l);
2566  bu_free(rp, "gapList items");
2567  }
2568  while (BU_LIST_WHILE (rp, region_pair, &exposedAirList.l)) {
2569  BU_LIST_DEQUEUE(&rp->l);
2570  bu_free(rp, "exposedAirList items");
2571  }
2572 
2573  /* Free dynamically allocated state */
2574  bu_free(state.m_lenDensity, "m_lenDensity");
2575  bu_free(state.m_len, "m_len");
2576  bu_free(state.m_volume, "m_volume");
2577  bu_free(state.m_weight, "m_weight");
2578  bu_free(state.shots, "m_shots");
2579  bu_free(state.m_lenTorque, "m_lenTorque");
2580  bu_free(state.m_moi, "m_moi");
2581  bu_free(state.m_poi, "m_poi");
2582 
2583  for (i = 0; i < num_objects; i++) {
2584  bu_free(obj_tbl[i].o_len, "o_len");
2585  bu_free(obj_tbl[i].o_lenDensity, "o_lenDensity");
2586  bu_free(obj_tbl[i].o_volume, "o_volume");
2587  bu_free(obj_tbl[i].o_weight, "o_weight");
2588  bu_free(obj_tbl[i].o_lenTorque, "o_lenTorque");
2589  bu_free(obj_tbl[i].o_moi, "o_moi");
2590  bu_free(obj_tbl[i].o_poi, "o_poi");
2591  }
2592  bu_free(obj_tbl, "object table");
2593  obj_tbl = NULL;
2594 
2595  for (i = 0, BU_LIST_FOR (regp, region, &(rtip->HeadRegion)), i++) {
2596  bu_free(reg_tbl[i].r_lenDensity, "r_lenDensity");
2597  bu_free(reg_tbl[i].r_len, "r_len");
2598  bu_free(reg_tbl[i].r_volume, "r_volume");
2599  bu_free(reg_tbl[i].r_weight, "r_weight");
2600  }
2601  bu_free(reg_tbl, "object table");
2602  reg_tbl = NULL;
2603 
2604  if (densities != NULL) {
2605  bu_free(densities, "densities");
2606  densities = NULL;
2607  }
2608 
2609  rt_free_rti(rtip);
2610 
2611  return GED_OK;
2612 }
2613 
2614 
2615 /*
2616  * Local Variables:
2617  * tab-width: 8
2618  * mode: C
2619  * indent-tabs-mode: t
2620  * c-file-style: "stroustrup"
2621  * End:
2622  * ex: shiftwidth=4 tabstop=8
2623  */
#define COMMA
Definition: gqa.c:75
Definition: db_flip.c:35
void usage(struct ged *gedp)
Definition: coil.c:315
#define GED_OK
Definition: ged.h:55
char * d_namep
pointer to name string
Definition: raytrace.h:859
struct xray a_ray
Actual ray to be shot.
Definition: raytrace.h:1583
#define BU_LIST_FOR(p, structure, hp)
Definition: list.h:365
#define ANALYSIS_VOLUME
Definition: gqa.c:59
void bu_log(const char *,...) _BU_ATTR_PRINTF12
Definition: log.c:176
#define STRCOMMA
Definition: gqa.c:76
int rt_db_get_internal(struct rt_db_internal *ip, const struct directory *dp, const struct db_i *dbip, const mat_t mat, struct resource *resp)
Definition: dir.c:76
#define GED_SEM_WORKER
Definition: ged.h:89
#define DLOG
Definition: gqa.c:119
struct region * pt_regionp
ptr to containing region
Definition: raytrace.h:580
#define DENSITY_MAGIC
Definition: analyze.h:55
struct hit * pt_outhit
OUT hit ptr.
Definition: raytrace.h:579
void rt_vlblock_free(struct bn_vlblock *vbp)
Definition: vlist.c:78
int read_units_double(double *val, char *buf, const struct cvt_tab *cvt)
Definition: gqa.c:410
fastf_t * m_poi
Definition: gqa.c:155
struct resource * resp
Definition: gqa.c:157
fastf_t * m_moi
Definition: gqa.c:154
#define ANALYSIS_MOMENTS
Definition: gqa.c:68
Definition: list.h:118
point_t mdl_min
min corner of model bounding RPP
Definition: raytrace.h:1769
void exposed_air(struct partition *pp, point_t last_out_point, point_t pt, point_t opt)
Definition: gqa.c:932
#define VOL
Definition: gqa.c:391
struct rt_i * rtip
Definition: gqa.c:148
#define ANALYSIS_BOX
Definition: gqa.c:65
#define RT_CK_APPLICATION(_p)
Definition: raytrace.h:1675
#define RT_CK_RTI(_p)
Definition: raytrace.h:1833
double dist
>= 0
Definition: tol.h:73
Definition: ged.h:338
void rt_free_rti(struct rt_i *rtip)
Definition: prep.c:156
struct db_i * dbip
Definition: raytrace.h:1266
struct bu_list * vhead
Definition: gqa.c:163
Definition: clone.c:90
unsigned long count
Definition: analyze.h:75
int useair
1="air" regions are retained while prepping
Definition: raytrace.h:1756
#define VSETALL(a, s)
Definition: color.c:54
struct bu_list HeadRegion
ptr of list of regions in model
Definition: raytrace.h:1778
Definition: raytrace.h:215
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 GED_CHECK_ARGC_GT_0(_gedp, _argc, _flags)
Definition: ged.h:202
struct directory * db_lookup(const struct db_i *, const char *name, int noisy)
Definition: db_lookup.c:153
Definition: raytrace.h:368
size_t nregions
total # of regions participating
Definition: raytrace.h:1782
Definition: raytrace.h:248
vect_t span
Definition: gqa.c:150
void logoverlap(struct application *ap, const struct partition *pp, const struct bu_ptbl *regiontable, const struct partition *InputHdp)
Definition: gqa.c:915
vect_t area
Definition: gqa.c:151
int u_axis
Definition: gqa.c:131
struct rt_wdb * ged_wdbp
Definition: ged.h:340
char * bu_optarg
Definition: globals.c:91
struct bu_list l
Definition: analyze.h:69
char * options_str
Definition: gqa.c:57
Header file for the BRL-CAD common definitions.
int bu_optind
Definition: globals.c:89
#define BU_CK_PTBL(_p)
Definition: ptbl.h:74
const char * reg_name
Identifying string.
Definition: raytrace.h:539
void list_report(struct region_pair *list)
Definition: gqa.c:1540
double grams_per_cu_mm
Definition: analyze.h:60
#define ANALYSIS_PLOT_OVERLAPS
Definition: gqa.c:69
struct resource * a_resource
dynamic memory resources
Definition: raytrace.h:1591
fastf_t * m_lenTorque
Definition: gqa.c:153
struct bn_vlblock * rt_vlblock_init(void)
Definition: vlist.c:71
#define BU_LIST_NON_EMPTY(hp)
Definition: list.h:296
int first
Definition: gqa.c:144
int bu_getopt(int nargc, char *const nargv[], const char *ostr)
Definition: getopt.c:43
double val
Definition: units.c:42
#define MAX_FASTF
Definition: defines.h:340
#define RT_CK_REGION(_p)
Definition: raytrace.h:559
vect_t u_dir
Definition: gqa.c:146
union rt_binunif_internal::@9 u
int(* a_hit)(struct application *, struct partition *, struct seg *)
called when shot hits model
Definition: raytrace.h:1584
#define GED_ERROR
Definition: ged.h:61
struct region * r2
Definition: analyze.h:74
struct hit * pt_inhit
IN hit pointer.
Definition: raytrace.h:577
void * bu_malloc(size_t siz, const char *str)
Definition: malloc.c:314
double * m_weight
Definition: gqa.c:142
Definition: ptbl.h:62
vect_t v_dir
Definition: gqa.c:147
void rt_init_resource(struct resource *resp, int cpu_num, struct rt_i *rtip)
Definition: prep.c:585
#define ANALYSIS_ADJ_AIR
Definition: gqa.c:62
struct region * r1
Definition: analyze.h:72
if(share_geom)
Definition: nmg_mod.c:3829
char * strchr(const char *sp, int c)
struct bu_list * rt_vlblock_find(struct bn_vlblock *vbp, int r, int g, int b)
Definition: vlist.c:98
int idb_major_type
Definition: raytrace.h:192
void bu_vls_free(struct bu_vls *vp)
Definition: vls.c:248
int bu_strncmp(const char *string1, const char *string2, size_t n)
Definition: str.c:191
Definition: color.c:49
void plane_worker(int cpu, void *ptr)
Definition: gqa.c:1311
struct ged_gqa_plot ged_gqa_plot
char name[32]
Definition: units.c:43
struct rt_i * a_rt_i
this librt instance
Definition: raytrace.h:1588
void * memset(void *s, int c, size_t n)
struct resource rt_uniresource
default. Defined in librt/globals.c
Definition: globals.c:41
#define GED_CHECK_DATABASE_OPEN(_gedp, _flags)
Definition: ged.h:114
#define RT_CHECK_BINUNIF(_p)
Definition: raytrace.h:1015
int reg_los
approximate line-of-sight thickness equivalence
Definition: raytrace.h:545
void * bu_calloc(size_t nelem, size_t elsize, const char *str)
Definition: malloc.c:321
struct density_entry * densities
Definition: gqa.c:167
double * m_len
Definition: gqa.c:140
#define LOOKUP_QUIET
Definition: raytrace.h:893
#define BN_VLIST_LINE_MOVE
Definition: vlist.h:82
#define RT_CK_AP(_p)
Definition: raytrace.h:1674
int a_user
application-specific value
Definition: raytrace.h:1617
#define V3ARGS(a)
Definition: color.c:56
char * strtok(char *s, const char *delim)
#define ANALYSIS_OVERLAPS
Definition: gqa.c:61
char * options
Definition: gqa.c:56
#define BN_VLIST_LINE_DRAW
Definition: vlist.h:83
int reg_gmater
GIFT Material code.
Definition: raytrace.h:544
#define BN_ADD_VLIST(_free_hd, _dest_hd, pnt, draw)
Definition: vlist.h:123
#define BU_FLSTR
Definition: defines.h:142
#define ANALYSIS_EXP_AIR
Definition: gqa.c:64
Definition: gqa.c:129
#define WGT
Definition: gqa.c:392
Definition: units.c:41
struct bn_vlblock * vbp
Definition: gqa.c:162
void allocate_per_region_data(struct cstate *state, int start, int ac, const char *av[])
Definition: gqa.c:1461
vect_t coord
Definition: analyze.h:77
void view_reports(struct cstate *state)
Definition: gqa.c:1703
struct partition * pt_back
backwards link
Definition: raytrace.h:575
int(* a_overlap)(struct application *, struct partition *, struct region *, struct region *, struct partition *)
DEPRECATED.
Definition: raytrace.h:1592
int curr_view
Definition: gqa.c:130
char * bu_vls_addr(const struct bu_vls *vp)
Definition: vls.c:111
struct bu_vls * ged_result_str
Definition: ged.h:357
#define BU_LIST_WHILE(p, structure, hp)
Definition: list.h:410
struct ged * _ged_current_gedp
Definition: loadview.c:37
int find_cmd_line_obj(struct per_obj_data *obj_rpt, const char *name)
Definition: gqa.c:1432
void bu_semaphore_release(unsigned int i)
Definition: semaphore.c:218
void pl_color(register FILE *plotfp, int r, int g, int b)
Definition: plot3.c:325
size_t bu_avail_cpus(void)
Definition: parallel.c:193
vect_t r_dir
Direction of ray (UNIT Length)
Definition: raytrace.h:219
void pdv_3line(register FILE *plotfp, const fastf_t *a, const fastf_t *b)
Definition: plot3.c:642
char * name
Definition: analyze.h:71
#define MAX_PSW
Definition: parallel.h:48
long steps[3]
Definition: gqa.c:149
void bu_parallel(void(*func)(int func_ncpu, void *func_data), int ncpu, void *data)
struct rt_i * rt_new_rti(struct db_i *dbip)
Definition: prep.c:58
void bn_mat_print_vls(const char *title, const mat_t m, struct bu_vls *vls)
Definition: mat.c:91
void _ged_cvt_vlblock_to_solids(struct ged *gedp, struct bn_vlblock *vbp, const char *name, int copy)
Definition: draw.c:554
int rt_gettree(struct rt_i *rtip, const char *node)
Definition: tree.c:869
union region_pair::@16 r
#define ZERO(val)
Definition: units.c:38
int terminate_check(struct cstate *state)
Definition: gqa.c:1880
uint32_t magic
Definition: analyze.h:59
int parse_densities_buffer(char *buf, size_t len, struct density_entry *densities, struct bu_vls *result_str, int *num_densities)
Definition: density.c:28
void * idb_ptr
Definition: raytrace.h:195
void(* a_logoverlap)(struct application *, const struct partition *, const struct bu_ptbl *, const struct partition *)
called to log overlaps
Definition: raytrace.h:1594
int overlap(struct application *ap, struct partition *pp, struct region *reg1, struct region *reg2, struct partition *hp)
Definition: gqa.c:843
point_t r_pt
Point at which ray starts.
Definition: raytrace.h:218
int get_densities_from_database(struct rt_i *rtip)
Definition: gqa.c:788
struct bn_tol wdb_tol
Definition: raytrace.h:1269
int reg_aircode
Region ID AIR code.
Definition: raytrace.h:543
int i_axis
Definition: gqa.c:133
int hit(struct application *ap, struct partition *PartHeadp, struct seg *segs)
Definition: gqa.c:963
#define BU_SEM_SYSCALL
Definition: parallel.h:178
axis
Definition: color.c:48
struct db_i * rti_dbip
prt to Database instance struct
Definition: raytrace.h:1774
int bu_opterr
Definition: globals.c:88
point_t mdl_max
max corner of model bounding RPP
Definition: raytrace.h:1770
#define RT_CK_PT(_p)
Definition: raytrace.h:589
int get_densities_from_file(char *name)
Definition: gqa.c:745
int v
Definition: gqa.c:136
void bu_vls_printf(struct bu_vls *vls, const char *fmt,...) _BU_ATTR_PRINTF23
Definition: vls.c:694
void summary_reports(struct cstate *state)
Definition: gqa.c:2005
Definition: analyze.h:58
int ged_gqa(struct ged *gedp, int argc, const char *argv[])
Definition: gqa.c:2321
double max_dist
Definition: analyze.h:76
#define GED_SEM_STATS
Definition: ged.h:90
#define GED_HELP
Definition: ged.h:62
int miss(struct application *ap)
Definition: gqa.c:1278
Definition: color.c:51
void rt_prep_parallel(struct rt_i *rtip, int ncpu)
int options_prep(struct rt_i *rtip, vect_t span)
Definition: gqa.c:1575
int rt_shootray(struct application *ap)
void bu_free(void *ptr, const char *str)
Definition: malloc.c:328
int(* a_miss)(struct application *)
called when shot misses
Definition: raytrace.h:1585
void * reg_udata
User appl. data for material.
Definition: raytrace.h:548
int v_axis
Definition: gqa.c:132
#define ANALYSIS_WEIGHT
Definition: gqa.c:60
#define BU_VLS_INIT_ZERO
Definition: vls.h:84
#define ANALYSIS_GAP
Definition: gqa.c:63
#define BU_LIST_DEQUEUE(cur)
Definition: list.h:209
struct bu_list * free_vlist_hd
where to get/put free vlists
Definition: vlist.h:175
#define BU_LIST_HEAD_MAGIC
Definition: magic.h:56
struct partition * pt_forw
forwards link
Definition: raytrace.h:574
#define RT_APPLICATION_INIT(_p)
Definition: raytrace.h:1676
fastf_t hit_dist
dist from r_pt to hit_point
Definition: raytrace.h:250
HIDDEN void verbose(struct human_data_t *dude)
Definition: human.c:2008
Definition: vls.h:56
unsigned long * shots
Definition: gqa.c:143
HIDDEN const point_t delta
Definition: sh_prj.c:618
#define LINE
Definition: gqa.c:390
double * m_lenDensity
Definition: gqa.c:139
double fastf_t
Definition: defines.h:300
#define ANALYSIS_CENTROIDS
Definition: gqa.c:67
#define VPRINT(a, b)
Definition: raytrace.h:1881
double * m_volume
Definition: gqa.c:141
int get_next_row(struct cstate *state)
Definition: gqa.c:1290
struct region_pair * add_unique_pair(struct region_pair *list, struct region *r1, struct region *r2, double dist, point_t pt)
Definition: overlaps.c:31
Definition: color.c:50
#define bu_strdup(s)
Definition: str.h:71
#define BU_STR_EQUAL(s1, s2)
Definition: str.h:126