BRL-CAD
tabdata.c
Go to the documentation of this file.
1 /* T A B D A T A . C
2  * BRL-CAD
3  *
4  * Copyright (c) 2004-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 
21 /** @addtogroup anim */
22 /** @{ */
23 /** @file libbn/tabdata.c
24  *
25  * @brief
26  * Routines for processing tables (curves) of data with one independent
27  * parameter which is common to many sets of dependent data values.
28  *
29  * Operates on bn_table (independent var) and
30  * bn_tabdata (dependent variable) structures.
31  *
32  * One application is for storing spectral curves, see spectrum.c
33  *
34  * @par Inspired by -
35  * Roy Hall and his book "Illumination and Color in Computer
36  *@n Generated Imagery", Springer Verlag, New York, 1989.
37  *@n ISBN 0-387-96774-5
38  *
39  * With thanks to Russ Moulton Jr, EOSoft Inc. for his "rad.c" module.
40  */
41 /** @} */
42 
43 #include "common.h"
44 
45 #include <math.h>
46 #include <string.h>
47 #include "bio.h"
48 
49 #include "bu/debug.h"
50 #include "bu/log.h"
51 #include "bu/malloc.h"
52 #include "bu/parallel.h"
53 #include "vmath.h"
54 #include "bn/tabdata.h"
55 
56 void
57 bn_table_free(struct bn_table *tabp)
58 {
59  if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_table_free(%p)\n", (void *)tabp);
60  BN_CK_TABLE(tabp);
61 
62  tabp->nx = 0; /* sanity */
63  bu_free(tabp, "struct bn_table");
64 }
65 
66 void
68 {
69  if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_tabdata_free(%p)\n", (void *)data);
70 
71  BN_CK_TABDATA(data);
72  BN_CK_TABLE(data->table);
73 
74  data->ny = 0; /* sanity */
75  data->table = NULL; /* sanity */
76  bu_free(data, "struct bn_tabdata");
77 }
78 
79 void
80 bn_ck_table(const struct bn_table *tabp)
81 {
82  register size_t i;
83 
84  if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_ck_table(%p)\n", (void *)tabp);
85 
86  BN_CK_TABLE(tabp);
87 
88  if (tabp->nx < 2)
89  bu_bomb("bn_ck_table() less than 2 wavelengths\n");
90 
91  for (i=0; i < tabp->nx; i++) {
92  if (tabp->x[i] >= tabp->x[i+1])
93  bu_bomb("bn_ck_table() wavelengths not in strictly ascending order\n");
94  }
95 }
96 
97 struct bn_table *
98 bn_table_make_uniform(size_t num, double first, double last) {
99  struct bn_table *tabp;
100  fastf_t *fp;
101  fastf_t delta;
102  int j;
103 
104  if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_table_make_uniform( num=%zu, %g, %g )\n", num, first, last);
105 
106  if (first >= last) bu_bomb("bn_table_make_uniform() first >= last\n");
107 
108  BN_GET_TABLE(tabp, num);
109 
110  delta = (last - first) / (double)num;
111 
112  fp = &tabp->x[0];
113  for (j = num; j > 0; j--) {
114  *fp++ = first;
115  first += delta;
116  }
117  tabp->x[num] = last;
118 
119  return tabp;
120 }
121 
122 void
123 bn_tabdata_add(struct bn_tabdata *out, const struct bn_tabdata *in1, const struct bn_tabdata *in2)
124 {
125  register int j;
126  register fastf_t *op;
127  register const fastf_t *i1, *i2;
128 
129  if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_tabdata_add(%p, %p, %p)\n", (void *)out, (void *)in1, (void *)in2);
130 
131  BN_CK_TABDATA(out);
132  BN_CK_TABDATA(in1);
133  BN_CK_TABDATA(in2);
134 
135  if (in1->table != in2->table || in1->table != out->table)
136  bu_bomb("bn_tabdata_add(): samples drawn from different tables\n");
137  if (in1->ny != in2->ny || in1->ny != out->ny)
138  bu_bomb("bn_tabdata_add(): different tabdata lengths?\n");
139 
140  op = out->y;
141  i1 = in1->y;
142  i2 = in2->y;
143  for (j = in1->ny; j > 0; j--)
144  *op++ = *i1++ + *i2++;
145  /* VADD2N( out->y, i1->y, i2->y, in1->ny ); */
146 }
147 
148 void
149 bn_tabdata_mul(struct bn_tabdata *out, const struct bn_tabdata *in1, const struct bn_tabdata *in2)
150 {
151  register int j;
152  register fastf_t *op;
153  register const fastf_t *i1, *i2;
154 
155  if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_tabdata_mul(%p, %p, %p)\n", (void *)out, (void *)in1, (void *)in2);
156 
157  BN_CK_TABDATA(out);
158  BN_CK_TABDATA(in1);
159  BN_CK_TABDATA(in2);
160 
161  if (in1->table != in2->table || in1->table != out->table)
162  bu_bomb("bn_tabdata_mul(): samples drawn from different tables\n");
163  if (in1->ny != in2->ny || in1->ny != out->ny)
164  bu_bomb("bn_tabdata_mul(): different tabdata lengths?\n");
165 
166  op = out->y;
167  i1 = in1->y;
168  i2 = in2->y;
169  for (j = in1->ny; j > 0; j--)
170  *op++ = *i1++ * *i2++;
171  /* VELMUL2N( out->y, i1->y, i2->y, in1->ny ); */
172 }
173 
174 void
175 bn_tabdata_mul3(struct bn_tabdata *out, const struct bn_tabdata *in1, const struct bn_tabdata *in2, const struct bn_tabdata *in3)
176 {
177  register int j;
178  register fastf_t *op;
179  register const fastf_t *i1, *i2, *i3;
180 
181  if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_tabdata_mul3(%p, %p, %p, %p)\n", (void *)out, (void *)in1, (void *)in2, (void *)in3);
182 
183  BN_CK_TABDATA(out);
184  BN_CK_TABDATA(in1);
185  BN_CK_TABDATA(in2);
186  BN_CK_TABDATA(in3);
187 
188  if (in1->table != in2->table || in1->table != out->table || in1->table != in2->table)
189  bu_bomb("bn_tabdata_mul(): samples drawn from different tables\n");
190  if (in1->ny != in2->ny || in1->ny != out->ny)
191  bu_bomb("bn_tabdata_mul(): different tabdata lengths?\n");
192 
193  op = out->y;
194  i1 = in1->y;
195  i2 = in2->y;
196  i3 = in3->y;
197  for (j = in1->ny; j > 0; j--)
198  *op++ = *i1++ * *i2++ * *i3++;
199  /* VELMUL3N( out->y, i1->y, i2->y, i3->y, in1->ny ); */
200 }
201 
202 void
203 bn_tabdata_incr_mul3_scale(struct bn_tabdata *out, const struct bn_tabdata *in1, const struct bn_tabdata *in2, const struct bn_tabdata *in3, register double scale)
204 {
205  register int j;
206  register fastf_t *op;
207  register const fastf_t *i1, *i2, *i3;
208 
209  if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_tabdata_incr_mul3_scale(%p, %p, %p, %p, %g)\n", (void *)out, (void *)in1, (void *)in2, (void *)in3, scale);
210 
211  BN_CK_TABDATA(out);
212  BN_CK_TABDATA(in1);
213  BN_CK_TABDATA(in2);
214  BN_CK_TABDATA(in3);
215 
216  if (in1->table != in2->table || in1->table != out->table || in1->table != in3->table)
217  bu_bomb("bn_tabdata_mul(): samples drawn from different tables\n");
218  if (in1->ny != in2->ny || in1->ny != out->ny)
219  bu_bomb("bn_tabdata_mul(): different tabdata lengths?\n");
220 
221  op = out->y;
222  i1 = in1->y;
223  i2 = in2->y;
224  i3 = in3->y;
225  for (j = in1->ny; j > 0; j--)
226  *op++ += *i1++ * *i2++ * *i3++ * scale;
227 }
228 
229 void
230 bn_tabdata_incr_mul2_scale(struct bn_tabdata *out, const struct bn_tabdata *in1, const struct bn_tabdata *in2, register double scale)
231 {
232  register int j;
233  register fastf_t *op;
234  register const fastf_t *i1, *i2;
235 
236  if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_tabdata_incr_mul2_scale(%p, %p, %p, %g)\n", (void *)out, (void *)in1, (void *)in2, scale);
237 
238  BN_CK_TABDATA(out);
239  BN_CK_TABDATA(in1);
240  BN_CK_TABDATA(in2);
241 
242  if (in1->table != in2->table || in1->table != out->table)
243  bu_bomb("bn_tabdata_mul(): samples drawn from different tables\n");
244  if (in1->ny != in2->ny || in1->ny != out->ny)
245  bu_bomb("bn_tabdata_mul(): different tabdata lengths?\n");
246 
247  op = out->y;
248  i1 = in1->y;
249  i2 = in2->y;
250  for (j = in1->ny; j > 0; j--)
251  *op++ += *i1++ * *i2++ * scale;
252 }
253 
254 void
255 bn_tabdata_scale(struct bn_tabdata *out, const struct bn_tabdata *in1, register double scale)
256 {
257  register int j;
258  register fastf_t *op;
259  register const fastf_t *i1;
260 
261  if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_tabdata_scale(%p, %p, %g)\n", (void *)out, (void *)in1, scale);
262 
263  BN_CK_TABDATA(out);
264  BN_CK_TABDATA(in1);
265 
266  if (in1->table != out->table)
267  bu_bomb("bn_tabdata_scale(): samples drawn from different tables\n");
268  if (in1->ny != out->ny)
269  bu_bomb("bn_tabdata_scale(): different tabdata lengths?\n");
270 
271  op = out->y;
272  i1 = in1->y;
273  for (j = in1->ny; j > 0; j--)
274  *op++ = *i1++ * scale;
275  /* VSCALEN( out->y, in->y, scale ); */
276 }
277 
278 void
279 bn_table_scale(struct bn_table *tabp, register double scale)
280 {
281  register int j;
282  register fastf_t *op;
283 
284  if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_table_scale(%p, %g)\n", (void *)tabp, scale);
285 
286  BN_CK_TABLE(tabp);
287 
288  op = tabp->x;
289  for (j = tabp->nx+1; j > 0; j--)
290  *op++ *= scale;
291  /* VSCALEN( tabp->x, tabp->x, scale, tabp->nx+1 ); */
292 }
293 
294 void
295 bn_tabdata_join1(struct bn_tabdata *out, const struct bn_tabdata *in1, register double scale, const struct bn_tabdata *in2)
296 {
297  register int j;
298  register fastf_t *op;
299  register const fastf_t *i1, *i2;
300 
301  if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_tabdata_join1(%p, %p, %g, %p)\n", (void *)out, (void *)in1, scale, (void *)in2);
302 
303  BN_CK_TABDATA(out);
304  BN_CK_TABDATA(in1);
305  BN_CK_TABDATA(in2);
306 
307  if (in1->table != out->table)
308  bu_bomb("bn_tabdata_join1(): samples drawn from different tables\n");
309  if (in1->table != in2->table)
310  bu_bomb("bn_tabdata_join1(): samples drawn from different tables\n");
311  if (in1->ny != out->ny)
312  bu_bomb("bn_tabdata_join1(): different tabdata lengths?\n");
313 
314  op = out->y;
315  i1 = in1->y;
316  i2 = in2->y;
317  for (j = in1->ny; j > 0; j--)
318  *op++ = *i1++ + scale * *i2++;
319  /* VJOIN1N( out->y, in1->y, scale, in2->y ); */
320 }
321 
322 void
323 bn_tabdata_join2(struct bn_tabdata *out, const struct bn_tabdata *in1, register double scale2, const struct bn_tabdata *in2, register double scale3, const struct bn_tabdata *in3)
324 {
325  register int j;
326  register fastf_t *op;
327  register const fastf_t *i1, *i2, *i3;
328 
329  if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_tabdata_join2(%p, %p, %g, %p, %g, %p)\n", (void *)out, (void *)in1, scale2, (void *)in2, scale3, (void *)in3);
330 
331  BN_CK_TABDATA(out);
332  BN_CK_TABDATA(in1);
333  BN_CK_TABDATA(in2);
334 
335  if (in1->table != out->table)
336  bu_bomb("bn_tabdata_join1(): samples drawn from different tables\n");
337  if (in1->table != in2->table)
338  bu_bomb("bn_tabdata_join1(): samples drawn from different tables\n");
339  if (in1->table != in3->table)
340  bu_bomb("bn_tabdata_join1(): samples drawn from different tables\n");
341  if (in1->ny != out->ny)
342  bu_bomb("bn_tabdata_join1(): different tabdata lengths?\n");
343 
344  op = out->y;
345  i1 = in1->y;
346  i2 = in2->y;
347  i3 = in3->y;
348  for (j = in1->ny; j > 0; j--)
349  *op++ = *i1++ + scale2 * *i2++ + scale3 * *i3++;
350  /* VJOIN2N( out->y, in1->y, scale2, in2->y, scale3, in3->y ); */
351 }
352 
353 void
354 bn_tabdata_blend2(struct bn_tabdata *out, register double scale1, const struct bn_tabdata *in1, register double scale2, const struct bn_tabdata *in2)
355 {
356  register int j;
357  register fastf_t *op;
358  register const fastf_t *i1, *i2;
359 
360  if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_tabdata_blend2(%p, %g, %p, %g, %p)\n", (void *)out, scale1, (void *)in1, scale2, (void *)in2);
361 
362  BN_CK_TABDATA(out);
363  BN_CK_TABDATA(in1);
364  BN_CK_TABDATA(in2);
365 
366  if (in1->table != out->table)
367  bu_bomb("bn_tabdata_blend2(): samples drawn from different tables\n");
368  if (in1->table != in2->table)
369  bu_bomb("bn_tabdata_blend2(): samples drawn from different tables\n");
370  if (in1->ny != out->ny)
371  bu_bomb("bn_tabdata_blend2(): different tabdata lengths?\n");
372 
373  op = out->y;
374  i1 = in1->y;
375  i2 = in2->y;
376  for (j = in1->ny; j > 0; j--)
377  *op++ = scale1 * *i1++ + scale2 * *i2++;
378  /* VBLEND2N( out->y, scale1, in1->y, scale2, in2->y ); */
379 }
380 
381 void
382 bn_tabdata_blend3(struct bn_tabdata *out, register double scale1, const struct bn_tabdata *in1, register double scale2, const struct bn_tabdata *in2, register double scale3, const struct bn_tabdata *in3)
383 {
384  register int j;
385  register fastf_t *op;
386  register const fastf_t *i1, *i2, *i3;
387 
388  if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_tabdata_blend3(%p, %g, %p, %g, %p, %g, %p)\n", (void *)out, scale1, (void *)in1, scale2, (void *)in2, scale3, (void *)in3);
389 
390  BN_CK_TABDATA(out);
391  BN_CK_TABDATA(in1);
392  BN_CK_TABDATA(in2);
393  BN_CK_TABDATA(in3);
394 
395  if (in1->table != out->table)
396  bu_bomb("bn_tabdata_blend3(): samples drawn from different tables\n");
397  if (in1->table != in2->table)
398  bu_bomb("bn_tabdata_blend3(): samples drawn from different tables\n");
399  if (in1->table != in3->table)
400  bu_bomb("bn_tabdata_blend3(): samples drawn from different tables\n");
401  if (in1->ny != out->ny)
402  bu_bomb("bn_tabdata_blend3(): different tabdata lengths?\n");
403 
404  op = out->y;
405  i1 = in1->y;
406  i2 = in2->y;
407  i3 = in3->y;
408  for (j = in1->ny; j > 0; j--)
409  *op++ = scale1 * *i1++ + scale2 * *i2++ + scale3 * *i3++;
410  /* VBLEND3N( out->y, scale1, in1->y, scale2, in2->y, scale3, in3->y ); */
411 }
412 
413 double
414 bn_tabdata_area1(const struct bn_tabdata *in)
415 {
416  register fastf_t area;
417  register const fastf_t *ip;
418  register int j;
419 
420  BN_CK_TABDATA(in);
421 
422  area = 0;
423  ip = in->y;
424  for (j = in->ny; j > 0; j--)
425  area += *ip++;
426 
427  if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_tabdata_area(%p) = %g\n", (void *)in, area);
428 
429  return area;
430 }
431 
432 double
433 bn_tabdata_area2(const struct bn_tabdata *in)
434 {
435  const struct bn_table *tabp;
436  register fastf_t area;
437  fastf_t width;
438  register int j;
439 
440  BN_CK_TABDATA(in);
441  tabp = in->table;
442  BN_CK_TABLE(tabp);
443 
444  area = 0;
445  for (j = in->ny-1; j >= 0; j--) {
446  width = tabp->x[j+1] - tabp->x[j];
447  area += in->y[j] * width;
448  }
449 
450  if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_tabdata_area2(%p) = %g\n", (void *)in, area);
451  return area;
452 }
453 
454 double
455 bn_tabdata_mul_area1(const struct bn_tabdata *in1, const struct bn_tabdata *in2)
456 {
457  register fastf_t area;
458  register const fastf_t *i1, *i2;
459  register int j;
460 
461  BN_CK_TABDATA(in1);
462  BN_CK_TABDATA(in2);
463 
464  area = 0;
465  i1 = in1->y;
466  i2 = in2->y;
467  for (j = in1->ny; j > 0; j--)
468  area += *i1++ * *i2++;
469 
470  if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_tabdata_mul_area1(%p, %p) = %g\n", (void *)in1, (void *)in2, area);
471  return area;
472 }
473 
474 double
475 bn_tabdata_mul_area2(const struct bn_tabdata *in1, const struct bn_tabdata *in2)
476 {
477  const struct bn_table *tabp;
478  register fastf_t area;
479  fastf_t width;
480  register int j;
481 
482  BN_CK_TABDATA(in1);
483  BN_CK_TABDATA(in2);
484  tabp = in1->table;
485  BN_CK_TABLE(tabp);
486 
487  area = 0;
488  for (j = in1->ny-1; j >= 0; j--) {
489  width = tabp->x[j+1] - tabp->x[j];
490  area += in1->y[j] * in2->y[j] * width;
491  }
492 
493  if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_tabdata_mul_area2(%p, %p) = %g\n", (void *)in1, (void *)in2, area);
494  return area;
495 }
496 
497 /*
498  *@brief
499  * Return the index in the table's x[] array of the interval which
500  * contains 'xval'.
501  *
502  *
503  * @return -1 if less than start value
504  * @return -2 if greater than end value.
505  * @return 0..nx-1 otherwise.
506  * Never returns a value of nx.
507  *
508  * A binary search would be more efficient, as the wavelengths (x values)
509  * are known to be sorted in ascending order.
510  */
511 long
512 bn_table_find_x(const struct bn_table *tabp, double xval)
513 {
514  register int i;
515 
516  if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_table_find_x(%p, %g)\n", (void *)tabp, xval);
517  BN_CK_TABLE(tabp);
518 
519  if (xval > tabp->x[tabp->nx]) return -2;
520  if (xval >= tabp->x[tabp->nx-1]) return tabp->nx-1;
521 
522  /* Search for proper interval in input spectrum */
523  for (i = tabp->nx-2; i >=0; i--) {
524  if (xval >= tabp->x[i]) return i;
525  }
526  /* if ( xval < tabp->x[0] ) return -1; */
527  return -1;
528 }
529 
530 fastf_t
531 bn_table_lin_interp(const struct bn_tabdata *samp, register double wl)
532 {
533  const struct bn_table *tabp;
534  size_t i;
535  int idx;
536  fastf_t fract;
537  fastf_t ret;
538 
539  if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_table_lin_interp(%p, %g)\n", (void *)samp, wl);
540 
541  BN_CK_TABDATA(samp);
542  tabp = samp->table;
543  BN_CK_TABLE(tabp);
544 
545  idx = bn_table_find_x(tabp, wl);
546  if (idx < 0) {
547  if (bu_debug&BU_DEBUG_TABDATA)
548  bu_log("bn_table_lin_interp(%g) out of range %g to %g\n", wl, tabp->x[0], tabp->x[tabp->nx]);
549  return 0;
550  }
551  i = (size_t)idx;
552 
553  if (wl < tabp->x[i] || wl >= tabp->x[i+1]) {
554  bu_log("bn_table_lin_interp(%g) assertion1 failed at %g\n", wl, tabp->x[i]);
555  bu_bomb("bn_table_lin_interp() assertion1 failed\n");
556  }
557 
558  if (i >= tabp->nx-2) {
559  /* Assume value is constant in final interval. */
560  if (bu_debug&BU_DEBUG_TABDATA)bu_log("bn_table_lin_interp(%g)=%g off end of range %g to %g\n", wl, samp->y[tabp->nx-1], tabp->x[0], tabp->x[tabp->nx]);
561  return samp->y[tabp->nx-1];
562  }
563 
564  /* The interval has been found */
565  fract = (wl - tabp->x[i]) / (tabp->x[i+1] - tabp->x[i]);
566  if (fract < 0 || fract > 1) bu_bomb("bn_table_lin_interp() assertion2 failed\n");
567  ret = (1-fract) * samp->y[i] + fract * samp->y[i+1];
568  if (bu_debug&BU_DEBUG_TABDATA)bu_log("bn_table_lin_interp(%g)=%g in range %g to %g\n",
569  wl, ret, tabp->x[i], tabp->x[i+1]);
570  return ret;
571 }
572 
573 struct bn_tabdata *
574 bn_tabdata_resample_max(const struct bn_table *newtable, const struct bn_tabdata *olddata) {
575  const struct bn_table *oldtable;
576  struct bn_tabdata *newsamp;
577  size_t i;
578  int j, k;
579 
580  if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_tabdata_resample_max(%p, %p)\n", (void *)newtable, (void *)olddata);
581 
582  BN_CK_TABLE(newtable);
583  BN_CK_TABDATA(olddata);
584  oldtable = olddata->table;
585  BN_CK_TABLE(oldtable);
586 
587  if (oldtable == newtable) bu_log("bn_tabdata_resample_max() NOTICE old and new bn_table structs are the same\n");
588 
589  BN_GET_TABDATA(newsamp, newtable);
590 
591  for (i = 0; i < newtable->nx; i++) {
592  /*
593  * Find good value(s) in olddata to represent the span from
594  * newtable->x[i] to newtable->x[i+1].
595  */
596  j = bn_table_find_x(oldtable, newtable->x[i]);
597  k = bn_table_find_x(oldtable, newtable->x[i+1]);
598  if (k == -1) {
599  /* whole new span is off left side of old table */
600  newsamp->y[i] = 0;
601  continue;
602  }
603  if (j == -2) {
604  /* whole new span is off right side of old table */
605  newsamp->y[i] = 0;
606  continue;
607  }
608 
609  if (j == k && j > 0) {
610  register fastf_t tmp;
611  /*
612  * Simple case, ends of output span are completely
613  * contained within one input span.
614  * Interpolate for both ends, take max.
615  */
616  newsamp->y[i] = bn_table_lin_interp(olddata, newtable->x[i]);
617  tmp = bn_table_lin_interp(olddata, newtable->x[i+1]);
618  V_MAX(newsamp->y[i], tmp);
619  } else {
620  register fastf_t tmp, n;
621  register int s;
622  /*
623  * Complex case: find good representative value.
624  * Interpolate both ends, and consider all
625  * intermediate old samples in span. Take max.
626  * One (but not both) new ends may be off old table.
627  */
628  n = bn_table_lin_interp(olddata, newtable->x[i]);
629  tmp = bn_table_lin_interp(olddata, newtable->x[i+1]);
630  V_MAX(n, tmp);
631 
632  for (s = j+1; s <= k; s++) {
633  if ((tmp = olddata->y[s]) > n)
634  n = tmp;
635  }
636  newsamp->y[i] = n;
637  }
638  }
639  return newsamp;
640 }
641 
642 struct bn_tabdata *
643 bn_tabdata_resample_avg(const struct bn_table *newtable, const struct bn_tabdata *olddata) {
644  const struct bn_table *oldtable;
645  struct bn_tabdata *newsamp;
646  size_t i;
647  int j, k;
648 
649  if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_tabdata_resample_avg(%p, %p)\n", (void *)newtable, (void *)olddata);
650 
651  BN_CK_TABLE(newtable);
652  BN_CK_TABDATA(olddata);
653  oldtable = olddata->table;
654  BN_CK_TABLE(oldtable);
655 
656  if (oldtable == newtable) bu_log("bn_tabdata_resample_avg() NOTICE old and new bn_table structs are the same\n");
657 
658  BN_GET_TABDATA(newsamp, newtable);
659 
660  for (i = 0; i < newtable->nx; i++) {
661  /*
662  * Find good value(s) in olddata to represent the span from
663  * newtable->x[i] to newtable->x[i+1].
664  */
665  j = bn_table_find_x(oldtable, newtable->x[i]);
666  k = bn_table_find_x(oldtable, newtable->x[i+1]);
667 
668  if (j < 0 || k < 0 || j == k) {
669  /*
670  * Simple case, ends of output span are completely
671  * contained within one input span.
672  * Interpolate for both ends, take average.
673  */
674  newsamp->y[i] = 0.5 * (
675  bn_table_lin_interp(olddata, newtable->x[i]) +
676  bn_table_lin_interp(olddata, newtable->x[i+1]));
677  } else {
678  /*
679  * Complex case: find average value.
680  * Interpolate both ends, and consider all
681  * intermediate old spans.
682  * There are three parts to sum:
683  * Partial interval from newx[i] to j+1
684  * Full intervals from j+1 to k
685  * Partial interval from k to newx[i+1]
686  */
687  fastf_t wsum; /* weighted sum */
688  fastf_t a, b; /* values being averaged */
689  int s;
690 
691  /* Partial interval from newx[i] to j+1 */
692  a = bn_table_lin_interp(olddata, newtable->x[i]); /* in "j" bin */
693  b = olddata->y[j+1];
694  wsum = 0.5 * (a+b) * (oldtable->x[j+1] - newtable->x[i]);
695 
696  /* Full intervals from j+1 to k */
697  for (s = j+1; s < k; s++) {
698  a = olddata->y[s];
699  b = olddata->y[s+1];
700  wsum += 0.5 * (a+b) * (oldtable->x[s+1] - oldtable->x[s]);
701  }
702 
703  /* Partial interval from k to newx[i+1] */
704  a = olddata->y[k];
705  b = bn_table_lin_interp(olddata, newtable->x[i+1]); /* in "k" bin */
706  wsum += 0.5 * (a+b) * (newtable->x[i+1] - oldtable->x[k]);
707 
708  /* Adjust the weighted sum by the total width */
709  newsamp->y[i] =
710  wsum / (newtable->x[i+1] - newtable->x[i]);
711  }
712  }
713  return newsamp;
714 }
715 
716 int
717 bn_table_write(const char *filename, const struct bn_table *tabp)
718 {
719  FILE *fp;
720  size_t j;
721 
722  if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_table_write(%s, %p)\n", filename, (void *)tabp);
723 
724  BN_CK_TABLE(tabp);
725 
727  fp = fopen(filename, "wb");
729 
730  if (fp == NULL) {
731  perror(filename);
732  bu_log("bn_table_write(%s, %p) FAILED\n", filename, (void *)tabp);
733  return -1;
734  }
735 
737  fprintf(fp, " %lu sample starts, and one end.\n", (long unsigned int)tabp->nx);
738  for (j=0; j <= tabp->nx; j++) {
739  fprintf(fp, "%g\n", tabp->x[j]);
740  }
741  fclose(fp);
743  return 0;
744 }
745 
746 struct bn_table *
747 bn_table_read(const char *filename) {
748  int ret;
749  struct bn_table *tabp;
750  struct bu_vls line = BU_VLS_INIT_ZERO;
751  FILE *fp;
752  size_t nw;
753  size_t j;
754 
755  if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_table_read(%s)\n", filename);
756 
758  fp = fopen(filename, "rb");
760 
761  if (fp == NULL) {
762  perror(filename);
763  bu_log("bn_table_read(%s) FAILED\n", filename);
764  return NULL;
765  }
766 
767  if (bu_vls_gets(&line, fp) < 0) {
768  perror(filename);
769  bu_log("Failed to read line\n");
770  fclose(fp);
771  return NULL;
772  }
773  nw = 0;
774  /* TODO: %lu to size_t isn't right. We may need a bu_sscanf() that can cope
775  * with native pointer width correctly, as %p is not ubiquitous and %z is a
776  * C99-ism */
777  sscanf(bu_vls_addr(&line), "%lu", (long unsigned int *)(&nw));
778  bu_vls_free(&line);
779 
780  if (nw <= 0) bu_bomb("bn_table_read() bad nw value\n");
781 
782  BN_GET_TABLE(tabp, nw);
783 
785  for (j=0; j <= tabp->nx; j++) {
786  double val;
787  ret = fscanf(fp, "%lf", &val);
788  if (ret != 1) {
789  bu_log("bn_table_read(%s) READ FAILURE. Abort\n", filename);
790  break;
791  }
792  tabp->x[j] = val;
793  }
794  fclose(fp);
796 
797  bn_ck_table(tabp);
798 
799  return tabp;
800 }
801 
802 void
803 bn_pr_table(const char *title, const struct bn_table *tabp)
804 {
805  size_t j;
806 
807  bu_log("%s\n", title);
808  BN_CK_TABLE(tabp);
809 
810  for (j=0; j <= tabp->nx; j++) {
811  bu_log("%3zu: %g\n", j, tabp->x[j]);
812  }
813 }
814 
815 void
816 bn_pr_tabdata(const char *title, const struct bn_tabdata *data)
817 {
818  size_t j;
819 
820  bu_log("%s: ", title);
821  BN_CK_TABDATA(data);
822 
823  for (j=0; j < data->ny; j++) {
824  bu_log("%g, ", data->y[j]);
825  }
826  bu_log("\n");
827 }
828 
829 int
831 {
832  FILE *fp;
833  const struct bn_table *tabp;
834  size_t j;
835 
836  if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_print_table_and_tabdata(%s, %p)\n", filename, (void *)data);
837 
838  BN_CK_TABDATA(data);
839  tabp = data->table;
840  BN_CK_TABLE(tabp);
841 
843  fp = fopen(filename, "wb");
845 
846  if (fp == NULL) {
847  perror(filename);
848  bu_log("bn_print_table_and_tabdata(%s, %p) FAILED\n", filename, (void *)data);
849  return -1;
850  }
851 
853  for (j=0; j < tabp->nx; j++) {
854  fprintf(fp, "%g %g\n", tabp->x[j], data->y[j]);
855  }
856  fprintf(fp, "%g (novalue)\n", tabp->x[tabp->nx]);
857  fclose(fp);
859  return 0;
860 }
861 
862 struct bn_tabdata *
864  struct bn_table *tabp;
865  struct bn_tabdata *data;
866  FILE *fp;
867  char buf[128];
868  size_t count = 0;
869  size_t i;
870 
871  if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_read_table_and_tabdata(%s)\n", filename);
872 
874  fp = fopen(filename, "rb");
876 
877  if (fp == NULL) {
878  perror(filename);
879  bu_log("bn_read_table_and_tabdata(%s) FAILED. Couldn't open file.\n", filename);
880  return NULL;
881  }
882 
883  /* First pass: Count number of lines */
885  for (;;) {
886  if (bu_fgets(buf, sizeof(buf), fp) == NULL) break;
887  count++;
888  }
889  fclose(fp);
891 
892  if (count < 2) {
893  perror(filename);
894  bu_log("bn_read_table_and_tabdata(%s) FAILED. Expected at least 2 lines in file.\n", filename);
895  return NULL;
896  }
897 
898  /* Allocate storage */
899  BN_GET_TABLE(tabp, count);
900  BN_GET_TABDATA(data, tabp);
901 
902  /* Second pass: Read only as much data as storage was allocated for */
904  fp = fopen(filename, "rb");
905  for (i=0; i < count; i++) {
906  double xval, yval;
907  int ret;
908 
909  buf[0] = '\0';
910  if (bu_fgets(buf, sizeof(buf), fp) == NULL) {
911  bu_log("bn_read_table_and_tabdata(%s) unexpected EOF on line %zu\n", filename, i);
912  break;
913  }
914  ret = sscanf(buf, "%lf %lf", &xval, &yval);
915  if (ret != 2) {
916  bu_log("Malformatted table data encountered in %s\n", filename);
917  break;
918  }
919  xval = tabp->x[i];
920  yval = data->y[i];
921  }
922  fclose(fp);
924 
925  /* Complete final interval */
926  tabp->x[count] = 2 * tabp->x[count-1] - tabp->x[count-2];
927 
928  bn_ck_table(tabp);
929 
930  return data;
931 }
932 
933 struct bn_tabdata *
934 bn_tabdata_binary_read(const char *filename, size_t num, const struct bn_table *tabp) {
935  struct bn_tabdata *data;
936  char *cp;
937  int nbytes;
938  int len;
939  int got;
940  int fd;
941  long i;
942 
943  if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_tabdata_binary_read(%s, num=%zu, %p)\n", filename, num, (void *)tabp);
944 
945  BN_CK_TABLE(tabp);
946 
947  nbytes = BN_SIZEOF_TABDATA(tabp);
948  len = num * nbytes;
949 
951  fd = open(filename, 0);
953  if (fd <= 0) {
954  perror(filename);
955  bu_log("bn_tabdata_binary_read open failed on \"%s\"\n", filename);
956  return (struct bn_tabdata *)NULL;
957  }
958 
959  /* Get a big array of structures for reading all at once */
960  data = (struct bn_tabdata *)bu_malloc(len+8, "bn_tabdata[]");
961 
963  got = read(fd, (char *)data, len);
965  if (got != len) {
966  if (got < 0) {
967  perror(filename);
968  bu_log("bn_tabdata_binary_read read error on \"%s\"\n", filename);
969  } else {
970  bu_log("bn_tabdata_binary_read(%s) expected %d got %d\n", filename, len, got);
971  }
972  bu_free(data, "bn_tabdata[]");
974  close(fd);
976  return (struct bn_tabdata *)NULL;
977  }
979  close(fd);
981 
982  /* Connect data[i].table pointer to tabp */
983  cp = (char *)data;
984  for (i = num-1; i >= 0; i--, cp += nbytes) {
985  register struct bn_tabdata *sp;
986  sp = (struct bn_tabdata *)cp;
987  BN_CK_TABDATA(sp);
988  sp->table = tabp;
989  }
990 
991  return data;
992 }
993 
994 struct bn_tabdata *
995 bn_tabdata_malloc_array(const struct bn_table *tabp, size_t num) {
996  struct bn_tabdata *data;
997  char *cp;
998  size_t i;
999  size_t nw;
1000  size_t nbytes;
1001 
1002  if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_tabdata_malloc_array(%p, num=%zu)\n", (void *)tabp, num);
1003 
1004  BN_CK_TABLE(tabp);
1005  nw = tabp->nx;
1006  nbytes = BN_SIZEOF_TABDATA(tabp);
1007 
1008  data = (struct bn_tabdata *)bu_calloc(num,
1009  nbytes, "struct bn_tabdata[]");
1010 
1011  cp = (char *)data;
1012  for (i = 0; i < num; i++) {
1013  register struct bn_tabdata *sp;
1014 
1015  sp = (struct bn_tabdata *)cp;
1016  sp->magic = BN_TABDATA_MAGIC;
1017  sp->ny = nw;
1018  sp->table = tabp;
1019  cp += nbytes;
1020  }
1021  return data;
1022 }
1023 
1024 void
1025 bn_tabdata_copy(struct bn_tabdata *out, const struct bn_tabdata *in)
1026 {
1027  if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_tabdata_copy(%p, %p)\n", (void *)out, (void *)in);
1028 
1029  BN_CK_TABDATA(out);
1030  BN_CK_TABDATA(in);
1031 
1032  if (in->table != out->table)
1033  bu_bomb("bn_tabdata_copy(): samples drawn from different tables\n");
1034  if (in->ny != out->ny)
1035  bu_bomb("bn_tabdata_copy(): different tabdata lengths?\n");
1036 
1037  memcpy((char *)out->y, (const char *)in->y, BN_SIZEOF_TABDATA_Y(in));
1038 }
1039 
1040 struct bn_tabdata *
1041 bn_tabdata_dup(const struct bn_tabdata *in) {
1042  struct bn_tabdata *data;
1043 
1044  BN_CK_TABDATA(in);
1045  BN_GET_TABDATA(data, in->table);
1046 
1047  memcpy((char *)data->y, (const char *)in->y, BN_SIZEOF_TABDATA_Y(in));
1048 
1049  if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_tabdata_dup(%p) = %p\n", (void *)in, (void *)data);
1050  return data;
1051 }
1052 
1053 struct bn_tabdata *
1054 bn_tabdata_get_constval(double val, const struct bn_table *tabp) {
1055  struct bn_tabdata *data;
1056  int todo;
1057  register fastf_t *op;
1058 
1059  BN_CK_TABLE(tabp);
1060  BN_GET_TABDATA(data, tabp);
1061 
1062  op = data->y;
1063  for (todo = data->ny-1; todo >= 0; todo--)
1064  *op++ = val;
1065 
1066  if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_tabdata_get_constval(val=%g, %p)=%p\n", val, (void *)tabp, (void *)data);
1067 
1068  return data;
1069 }
1070 
1071 void
1073 {
1074  int todo;
1075  register fastf_t *op;
1076 
1077  if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_tabdata_constval(%p, val=%g)\n", (void *)data, val);
1078 
1079  BN_CK_TABDATA(data);
1080 
1081  op = data->y;
1082  for (todo = data->ny-1; todo >= 0; todo--)
1083  *op++ = val;
1084 }
1085 
1086 void
1087 bn_tabdata_to_tcl(struct bu_vls *vp, const struct bn_tabdata *data)
1088 {
1089  const struct bn_table *tabp;
1090  register size_t i;
1091  register fastf_t minval = MAX_FASTF, maxval = -MAX_FASTF;
1092 
1093  if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_tabdata_to_tcl(%p, %p)\n", (void *)vp, (void *)data);
1094 
1095  BU_CK_VLS(vp);
1096  BN_CK_TABDATA(data);
1097  tabp = data->table;
1098  BN_CK_TABLE(tabp);
1099 
1100  bu_vls_strcat(vp, "x {");
1101  for (i=0; i < tabp->nx; i++) {
1102  bu_vls_printf(vp, "%g ", tabp->x[i]);
1103  }
1104  bu_vls_strcat(vp, "} y {");
1105  for (i=0; i < data->ny; i++) {
1106  register fastf_t val = data->y[i];
1107  bu_vls_printf(vp, "%g ", val);
1108  V_MIN(minval, val);
1109  V_MAX(maxval, val);
1110  }
1111  bu_vls_printf(vp, "} nx %zu ymin %g ymax %g",
1112  tabp->nx, minval, maxval);
1113 }
1114 
1115 struct bn_tabdata *
1116 bn_tabdata_from_array(const double *array) {
1117  register const double *dp;
1118  size_t len = 0;
1119  struct bn_table *tabp;
1120  struct bn_tabdata *data;
1121  register size_t i;
1122 
1123  /* First, find len */
1124  for (dp = array; *dp > 0; dp += 2)
1125  ; /* NIL */
1126  len = (dp - array) >> 1;
1127 
1128  /* Second, build bn_table */
1129  BN_GET_TABLE(tabp, len);
1130  for (i = 0; i < len; i++) {
1131  tabp->x[i] = array[i<<1];
1132  }
1133  tabp->x[len] = tabp->x[len-1] + 1; /* invent span end */
1134 
1135  /* Third, build bn_tabdata (last input "y" is ignored) */
1136  BN_GET_TABDATA(data, tabp);
1137  for (i = 0; i < len-1; i++) {
1138  data->y[i] = array[(i<<1)+1];
1139  }
1140 
1141  if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_tabdata_from_array(%p) = %p\n", (void *)array, (void *)data);
1142  return data;
1143 }
1144 
1145 void
1146 bn_tabdata_freq_shift(struct bn_tabdata *out, const struct bn_tabdata *in, double offset)
1147 {
1148  const struct bn_table *tabp;
1149  register size_t i;
1150 
1151  if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_tabdata_freq_shift(%p, %p, offset=%g)\n", (void *)out, (void *)in, offset);
1152 
1153  BN_CK_TABDATA(out);
1154  BN_CK_TABDATA(in);
1155  tabp = in->table;
1156 
1157  if (tabp != out->table)
1158  bu_bomb("bn_tabdata_freq_shift(): samples drawn from different tables\n");
1159  if (in->ny != out->ny)
1160  bu_bomb("bn_tabdata_freq_shift(): different tabdata lengths?\n");
1161 
1162  for (i=0; i < out->ny; i++) {
1163  out->y[i] = bn_table_lin_interp(in, tabp->x[i]+offset);
1164  }
1165 }
1166 
1167 int
1168 bn_table_interval_num_samples(const struct bn_table *tabp, double low, double hi)
1169 {
1170  register size_t i;
1171  register size_t count = 0;
1172 
1173  BN_CK_TABLE(tabp);
1174 
1175  for (i=0; i < tabp->nx-1; i++) {
1176  if (tabp->x[i] >= low && tabp->x[i] <= hi) count++;
1177  }
1178 
1179  if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_table_interval_num_samples(%p, low=%g, hi=%g) = %zu\n", (void *)tabp, low, hi, count);
1180  return count;
1181 }
1182 
1183 int
1184 bn_table_delete_sample_pts(struct bn_table *tabp, unsigned int i, unsigned int j)
1185 {
1186  size_t tokill;
1187  unsigned int k;
1188 
1189  if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_table_delete_samples(%p, %d, %d)\n", (void *)tabp, i, j);
1190 
1191  BN_CK_TABLE(tabp);
1192 
1193  if (i >= tabp->nx || j >= tabp->nx)
1194  bu_bomb("bn_table_delete_sample_pts() index out of range\n");
1195 
1196  tokill = j - i + 1;
1197  if (tokill < 1) bu_bomb("bn_table_delete_sample_pts(): nothing to delete\n");
1198  if (tokill >= tabp->nx) bu_bomb("bn_table_delete_sample_pts(): you can't kill 'em all!\n");
1199 
1200  tabp->nx -= tokill;
1201 
1202  for (k = i; k < tabp->nx; k++) {
1203  tabp->x[k] = tabp->x[k+tokill];
1204  }
1205  return tokill;
1206 }
1207 
1208 struct bn_table *
1209 bn_table_merge2(const struct bn_table *a, const struct bn_table *b) {
1210  struct bn_table *newtab;
1211  register size_t i, j, k;
1212 
1213  BN_CK_TABLE(a);
1214  BN_CK_TABLE(b);
1215 
1216  BN_GET_TABLE(newtab, a->nx + b->nx + 2);
1217 
1218  i = j = 0; /* input subscripts */
1219  k = 0; /* output subscript */
1220  while (i <= a->nx || j <= b->nx) {
1221  if (i > a->nx) {
1222  while (j <= b->nx)
1223  newtab->x[k++] = b->x[j++];
1224  break;
1225  }
1226  if (j > b->nx) {
1227  while (i <= a->nx)
1228  newtab->x[k++] = a->x[i++];
1229  break;
1230  }
1231  /* Both have remaining elements, take lower one */
1232  if (ZERO(a->x[i] - b->x[j])) {
1233  newtab->x[k++] = a->x[i++];
1234  j++; /* compress out duplicate */
1235  continue;
1236  }
1237  if (a->x[i] <= b->x[j]) {
1238  newtab->x[k++] = a->x[i++];
1239  } else {
1240  newtab->x[k++] = b->x[j++];
1241  }
1242  }
1243  if (k > newtab->nx)
1244  bu_bomb("bn_table_merge2() assertion failed, k>nx?\n");
1245  newtab->nx = k-1;
1246 
1247  if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_table_merge2(%p, %p) = %p\n", (void *)a, (void *)b, (void *)newtab);
1248  return newtab;
1249 }
1250 
1251 struct bn_tabdata *
1252 bn_tabdata_mk_linear_filter(const struct bn_table *spectrum, double lower_wavelen, double upper_wavelen) {
1253  struct bn_tabdata *filt;
1254  int first;
1255  int last;
1256  fastf_t filt_range;
1257  fastf_t cell_range;
1258  fastf_t frac;
1259  size_t i;
1260 
1261  BN_CK_TABLE(spectrum);
1262 
1263  if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_tabdata_mk_linear_filter(%p, low=%g, up=%g)\n", (void *)spectrum, lower_wavelen, upper_wavelen);
1264 
1265  if (lower_wavelen < spectrum->x[0])
1266  if (upper_wavelen > spectrum->x[spectrum->nx])
1267  bu_log("bn_tabdata_mk_linear_filter() warning, upper_wavelen %g > highest sampled wavelen %g\n",
1268  upper_wavelen, spectrum->x[spectrum->nx]);
1269 
1270  /* First, find first (possibly partial) sample */
1271  first = bn_table_find_x(spectrum, lower_wavelen);
1272  if (first == -1) {
1273  first = 0;
1274  bu_log("bn_tabdata_mk_linear_filter() warning, lower_wavelen %g < lowest sampled wavelen %g\n",
1275  lower_wavelen, spectrum->x[0]);
1276  } else if (first <= -2) {
1277  bu_log("bn_tabdata_mk_linear_filter() ERROR, lower_wavelen %g > highest sampled wavelen %g\n",
1278  lower_wavelen, spectrum->x[spectrum->nx]);
1279  return NULL;
1280  }
1281 
1282  /* Second, find last (possibly partial) sample */
1283  last = bn_table_find_x(spectrum, upper_wavelen);
1284  if (last == -1) {
1285  bu_log("bn_tabdata_mk_linear_filter() ERROR, upper_wavelen %g < lowest sampled wavelen %g\n",
1286  upper_wavelen, spectrum->x[0]);
1287  return NULL;
1288  } else if (last <= -2) {
1289  last = spectrum->nx-1;
1290  bu_log("bn_tabdata_mk_linear_filter() warning, upper_wavelen %g > highest sampled wavelen %g\n",
1291  upper_wavelen, spectrum->x[spectrum->nx]);
1292  }
1293 
1294  /* 'filt' is filled with zeros by default */
1295  BN_GET_TABDATA(filt, spectrum);
1296 
1297  /* Special case: first and last are in same sample cell */
1298  if (first == last) {
1299  filt_range = upper_wavelen - lower_wavelen;
1300  cell_range = spectrum->x[first+1] - spectrum->x[first];
1301  frac = filt_range / cell_range;
1302 
1303  /* Could use a BU_ASSERT_RANGE_FLOAT */
1304  BU_ASSERT((frac >= 0.0) && (frac <= 1.0));
1305 
1306  filt->y[first] = frac;
1307  return filt;
1308  }
1309 
1310  /* Calculate fraction 0..1.0 for first and last samples */
1311  filt_range = spectrum->x[first+1] - lower_wavelen;
1312  cell_range = spectrum->x[first+1] - spectrum->x[first];
1313  frac = filt_range / cell_range;
1314  V_MIN(frac, 1);
1315 
1316  BU_ASSERT((frac >= 0.0) && (frac <= 1.0));
1317  filt->y[first] = frac;
1318 
1319  filt_range = upper_wavelen - spectrum->x[last];
1320  cell_range = spectrum->x[last+1] - spectrum->x[last];
1321  frac = filt_range / cell_range;
1322  V_MIN(frac, 1);
1323 
1324  BU_ASSERT((frac >= 0.0) && (frac <= 1.0));
1325  filt->y[last] = frac;
1326 
1327  /* Fill in range between with 1.0 values */
1328  for (i = first+1; i < (size_t)last; i++)
1329  filt->y[i] = 1.0;
1330 
1331  return filt;
1332 }
1333 /** @} */
1334 /*
1335  * Local Variables:
1336  * mode: C
1337  * tab-width: 8
1338  * indent-tabs-mode: t
1339  * c-file-style: "stroustrup"
1340  * End:
1341  * ex: shiftwidth=4 tabstop=8
1342  */
void bn_table_free(struct bn_table *tabp)
Definition: tabdata.c:57
Definition: db_flip.c:35
char filename[MAXLENGTH]
Definition: human.c:105
struct bn_tabdata * bn_tabdata_resample_max(const struct bn_table *newtable, const struct bn_tabdata *olddata)
Definition: tabdata.c:574
void bu_log(const char *,...) _BU_ATTR_PRINTF12
Definition: log.c:176
int bn_table_interval_num_samples(const struct bn_table *tabp, double low, double hi)
Definition: tabdata.c:1168
struct bn_tabdata * bn_tabdata_malloc_array(const struct bn_table *tabp, size_t num)
Definition: tabdata.c:995
#define BN_GET_TABDATA(_data, _table)
Definition: tabdata.h:114
void bn_tabdata_free(struct bn_tabdata *data)
Definition: tabdata.c:67
void bn_tabdata_incr_mul2_scale(struct bn_tabdata *out, const struct bn_tabdata *in1, const struct bn_tabdata *in2, register double scale)
Definition: tabdata.c:230
void bn_tabdata_incr_mul3_scale(struct bn_tabdata *out, const struct bn_tabdata *in1, const struct bn_tabdata *in2, const struct bn_tabdata *in3, register double scale)
Definition: tabdata.c:203
if lu s
Definition: nmg_mod.c:3860
size_t ny
Definition: tabdata.h:102
void bu_vls_strcat(struct bu_vls *vp, const char *s)
Definition: vls.c:368
void bu_semaphore_acquire(unsigned int i)
Definition: semaphore.c:180
int bn_table_write(const char *filename, const struct bn_table *tabp)
Definition: tabdata.c:717
Header file for the BRL-CAD common definitions.
#define BU_ASSERT(_equation)
Definition: defines.h:216
double bn_tabdata_area2(const struct bn_tabdata *in)
Definition: tabdata.c:433
#define MAX_FASTF
Definition: defines.h:340
ustring width
#define BN_CK_TABDATA(_p)
Definition: tabdata.h:106
void * bu_malloc(size_t siz, const char *str)
Definition: malloc.c:314
struct bn_table * bn_table_read(const char *filename)
Definition: tabdata.c:747
#define BN_TABDATA_MAGIC
Definition: magic.h:72
void bn_tabdata_copy(struct bn_tabdata *out, const struct bn_tabdata *in)
Definition: tabdata.c:1025
if(share_geom)
Definition: nmg_mod.c:3829
void bn_tabdata_to_tcl(struct bu_vls *vp, const struct bn_tabdata *data)
Definition: tabdata.c:1087
void bu_vls_free(struct bu_vls *vp)
Definition: vls.c:248
COMPLEX data[64]
Definition: fftest.c:34
size_t nx
Definition: tabdata.h:75
#define BU_CK_VLS(_vp)
Definition: vls.h:69
#define BN_SIZEOF_TABDATA(_table)
Definition: tabdata.h:110
#define BN_CK_TABLE(_p)
Definition: tabdata.h:79
void bn_tabdata_freq_shift(struct bn_tabdata *out, const struct bn_tabdata *in, double offset)
Definition: tabdata.c:1146
fastf_t bn_table_lin_interp(const struct bn_tabdata *samp, register double wl)
Definition: tabdata.c:531
void * bu_calloc(size_t nelem, size_t elsize, const char *str)
Definition: malloc.c:321
fastf_t y[1]
array of ny samples, dynamically sized
Definition: tabdata.h:104
struct bn_table * bn_table_make_uniform(size_t num, double first, double last)
Definition: tabdata.c:98
#define BN_SIZEOF_TABDATA_Y(_tabdata)
Definition: tabdata.h:109
void bn_tabdata_mul(struct bn_tabdata *out, const struct bn_tabdata *in1, const struct bn_tabdata *in2)
Definition: tabdata.c:149
struct bn_tabdata * bn_tabdata_resample_avg(const struct bn_table *newtable, const struct bn_tabdata *olddata)
Definition: tabdata.c:643
long bn_table_find_x(const struct bn_table *tabp, double xval)
Definition: tabdata.c:512
void bn_table_scale(struct bn_table *tabp, register double scale)
Definition: tabdata.c:279
struct bn_tabdata * bn_tabdata_get_constval(double val, const struct bn_table *tabp)
Definition: tabdata.c:1054
void bn_tabdata_add(struct bn_tabdata *out, const struct bn_tabdata *in1, const struct bn_tabdata *in2)
Definition: tabdata.c:123
int bn_print_table_and_tabdata(const char *filename, const struct bn_tabdata *data)
Definition: tabdata.c:830
struct bn_tabdata * bn_read_table_and_tabdata(const char *filename)
Definition: tabdata.c:863
goto out
Definition: nmg_mod.c:3846
char * bu_vls_addr(const struct bu_vls *vp)
Definition: vls.c:111
double bn_tabdata_mul_area2(const struct bn_tabdata *in1, const struct bn_tabdata *in2)
Definition: tabdata.c:475
struct bn_tabdata * bn_tabdata_mk_linear_filter(const struct bn_table *spectrum, double lower_wavelen, double upper_wavelen)
Definition: tabdata.c:1252
void bn_tabdata_constval(struct bn_tabdata *data, double val)
Definition: tabdata.c:1072
void bu_semaphore_release(unsigned int i)
Definition: semaphore.c:218
struct bn_tabdata * bn_tabdata_binary_read(const char *filename, size_t num, const struct bn_table *tabp)
Definition: tabdata.c:934
void bn_pr_tabdata(const char *title, const struct bn_tabdata *data)
Definition: tabdata.c:816
struct bn_table * bn_table_merge2(const struct bn_table *a, const struct bn_table *b)
Definition: tabdata.c:1209
int bu_vls_gets(struct bu_vls *vp, FILE *fp)
Definition: vls.c:621
double bn_tabdata_mul_area1(const struct bn_tabdata *in1, const struct bn_tabdata *in2)
Definition: tabdata.c:455
void bn_pr_table(const char *title, const struct bn_table *tabp)
Definition: tabdata.c:803
#define ZERO(val)
Definition: units.c:38
void bn_tabdata_join2(struct bn_tabdata *out, const struct bn_tabdata *in1, register double scale2, const struct bn_tabdata *in2, register double scale3, const struct bn_tabdata *in3)
Definition: tabdata.c:323
int bu_debug
Definition: globals.c:87
#define BU_SEM_SYSCALL
Definition: parallel.h:178
fastf_t x[1]
array of nx+1 wavelengths, dynamically sized
Definition: tabdata.h:76
uint32_t magic
Definition: tabdata.h:101
void bu_vls_printf(struct bu_vls *vls, const char *fmt,...) _BU_ATTR_PRINTF23
Definition: vls.c:694
void bn_tabdata_mul3(struct bn_tabdata *out, const struct bn_tabdata *in1, const struct bn_tabdata *in2, const struct bn_tabdata *in3)
Definition: tabdata.c:175
struct bn_tabdata * bn_tabdata_dup(const struct bn_tabdata *in)
Definition: tabdata.c:1041
void bu_free(void *ptr, const char *str)
Definition: malloc.c:328
void bn_tabdata_join1(struct bn_tabdata *out, const struct bn_tabdata *in1, register double scale, const struct bn_tabdata *in2)
Definition: tabdata.c:295
struct bn_tabdata * bn_tabdata_from_array(const double *array)
Definition: tabdata.c:1116
struct bn_table * spectrum
Definition: init.c:41
char * bu_fgets(char *s, int size, FILE *stream)
Definition: fgets.c:31
int bn_table_delete_sample_pts(struct bn_table *tabp, unsigned int i, unsigned int j)
Definition: tabdata.c:1184
#define BU_VLS_INIT_ZERO
Definition: vls.h:84
#define BU_DEBUG_TABDATA
Definition: debug.h:73
double bn_tabdata_area1(const struct bn_tabdata *in)
Definition: tabdata.c:414
Definition: vls.h:56
void bu_bomb(const char *str) _BU_ATTR_NORETURN
Definition: bomb.c:91
HIDDEN const point_t delta
Definition: sh_prj.c:618
double fastf_t
Definition: defines.h:300
void bn_tabdata_blend2(struct bn_tabdata *out, register double scale1, const struct bn_tabdata *in1, register double scale2, const struct bn_tabdata *in2)
Definition: tabdata.c:354
const struct bn_table * table
Up pointer to definition of X axis.
Definition: tabdata.h:103
void bn_tabdata_scale(struct bn_tabdata *out, const struct bn_tabdata *in1, register double scale)
Definition: tabdata.c:255
void bn_ck_table(const struct bn_table *tabp)
Definition: tabdata.c:80
void bn_tabdata_blend3(struct bn_tabdata *out, register double scale1, const struct bn_tabdata *in1, register double scale2, const struct bn_tabdata *in2, register double scale3, const struct bn_tabdata *in3)
Definition: tabdata.c:382
#define BN_GET_TABLE(_table, _nx)
Definition: tabdata.h:84