tabdata.c

Go to the documentation of this file.
00001 /*                       T A B D A T A . C
00002  * BRL-CAD
00003  *
00004  * Copyright (c) 2004-2012 United States Government as represented by
00005  * the U.S. Army Research Laboratory.
00006  *
00007  * This library is free software; you can redistribute it and/or
00008  * modify it under the terms of the GNU Lesser General Public License
00009  * version 2.1 as published by the Free Software Foundation.
00010  *
00011  * This library is distributed in the hope that it will be useful, but
00012  * WITHOUT ANY WARRANTY; without even the implied warranty of
00013  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00014  * Lesser General Public License for more details.
00015  *
00016  * You should have received a copy of the GNU Lesser General Public
00017  * License along with this file; see the file named COPYING for more
00018  * information.
00019  */
00020 /** @addtogroup anim */
00021 /** @{ */
00022 /** @file libbn/tabdata.c
00023  *
00024  * @brief
00025  *  Routines for processing tables (curves) of data with one independent
00026  *  parameter which is common to many sets of dependent data values.
00027  *
00028  *  Operates on bn_table (independent var) and
00029  *  bn_tabdata (dependent variable) structures.
00030  *
00031  *  One application is for storing spectral curves, see spectrum.c
00032  *
00033  *  @par Inspired by -
00034  *      Roy Hall and his book "Illumination and Color in Computer
00035  *@n    Generated Imagery", Springer Verlag, New York, 1989.
00036  *@n    ISBN 0-387-96774-5
00037  *
00038  *  With thanks to Russ Moulton Jr, EOSoft Inc. for his "rad.c" module.
00039  */
00040 /** @} */
00041 
00042 #include "common.h"
00043 
00044 #include <math.h>
00045 #include <string.h>
00046 #include "bio.h"
00047 
00048 #include "vmath.h"
00049 #include "bu.h"
00050 #include "bn.h"
00051 
00052 /*
00053  *                      B N _ T A B L E _ F R E E
00054  */
00055 void
00056 bn_table_free(struct bn_table *tabp)
00057 {
00058     if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_table_free(%p)\n", (void *)tabp);
00059     BN_CK_TABLE(tabp);
00060 
00061     tabp->nx = 0;                       /* sanity */
00062     bu_free(tabp, "struct bn_table");
00063 }
00064 
00065 /*
00066  *                      B N _ T A B D A T A _ F R E E
00067  */
00068 void
00069 bn_tabdata_free(struct bn_tabdata *data)
00070 {
00071     if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_tabdata_free(%p)\n", (void *)data);
00072 
00073     BN_CK_TABDATA(data);
00074     BN_CK_TABLE(data->table);
00075 
00076     data->ny = 0;                       /* sanity */
00077     data->table = NULL;         /* sanity */
00078     bu_free(data, "struct bn_tabdata");
00079 }
00080 
00081 /*
00082  *                      B N _ C K _ T A B L E
00083  */
00084 void
00085 bn_ck_table(const struct bn_table *tabp)
00086 {
00087     register size_t i;
00088 
00089     if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_ck_table(%p)\n", (void *)tabp);
00090 
00091     BN_CK_TABLE(tabp);
00092 
00093     if (tabp->nx < 2)
00094         bu_bomb("bn_ck_table() less than 2 wavelengths\n");
00095 
00096     for (i=0; i < tabp->nx; i++)  {
00097         if (tabp->x[i] >= tabp->x[i+1])
00098             bu_bomb("bn_ck_table() wavelengths not in strictly ascending order\n");
00099     }
00100 }
00101 
00102 /*
00103  *                      B N _ T A B L E _ M A K E _ U N I F O R M
00104  *@brief
00105  *  Set up an independent "table margin" from 'first' to 'last',
00106  *  inclusive, using 'num' uniformly spaced samples.  Num >= 1.
00107  */
00108 struct bn_table *
00109 bn_table_make_uniform(size_t num, double first, double last) {
00110     struct bn_table     *tabp;
00111     fastf_t                     *fp;
00112     fastf_t                     delta;
00113     int                 j;
00114 
00115     if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_table_make_uniform( num=%zu, %g, %g )\n", num, first, last);
00116 
00117     if (first >= last)  bu_bomb("bn_table_make_uniform() first >= last\n");
00118 
00119     BN_GET_TABLE(tabp, num);
00120 
00121     delta = (last - first) / (double)num;
00122 
00123     fp = &tabp->x[0];
00124     for (j = num; j > 0; j--)  {
00125         *fp++ = first;
00126         first += delta;
00127     }
00128     tabp->x[num] = last;
00129 
00130     return tabp;
00131 }
00132 
00133 /*
00134  *                      B N _ T A B D A T A _ A D D
00135  *@brief
00136  *  Sum the values from two data tables.
00137  */
00138 void
00139 bn_tabdata_add(struct bn_tabdata *out, const struct bn_tabdata *in1, const struct bn_tabdata *in2)
00140 {
00141     register int                j;
00142     register fastf_t    *op;
00143     register const fastf_t      *i1, *i2;
00144 
00145     if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_tabdata_add(%p, %p, %p)\n", (void *)out, (void *)in1, (void *)in2);
00146 
00147     BN_CK_TABDATA(out);
00148     BN_CK_TABDATA(in1);
00149     BN_CK_TABDATA(in2);
00150 
00151     if (in1->table != in2->table || in1->table != out->table)
00152         bu_bomb("bn_tabdata_add(): samples drawn from different tables\n");
00153     if (in1->ny != in2->ny || in1->ny != out->ny)
00154         bu_bomb("bn_tabdata_add(): different tabdata lengths?\n");
00155 
00156     op = out->y;
00157     i1 = in1->y;
00158     i2 = in2->y;
00159     for (j = in1->ny; j > 0; j--)
00160         *op++ = *i1++ + *i2++;
00161     /* VADD2N( out->y, i1->y, i2->y, in1->ny ); */
00162 }
00163 
00164 /*
00165  *                      B N _ T A B D A T A _ M U L
00166  *@brief
00167  *  Element-by-element multiply the values from two data tables.
00168  */
00169 void
00170 bn_tabdata_mul(struct bn_tabdata *out, const struct bn_tabdata *in1, const struct bn_tabdata *in2)
00171 {
00172     register int                j;
00173     register fastf_t    *op;
00174     register const fastf_t      *i1, *i2;
00175 
00176     if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_tabdata_mul(%p, %p, %p)\n", (void *)out, (void *)in1, (void *)in2);
00177 
00178     BN_CK_TABDATA(out);
00179     BN_CK_TABDATA(in1);
00180     BN_CK_TABDATA(in2);
00181 
00182     if (in1->table != in2->table || in1->table != out->table)
00183         bu_bomb("bn_tabdata_mul(): samples drawn from different tables\n");
00184     if (in1->ny != in2->ny || in1->ny != out->ny)
00185         bu_bomb("bn_tabdata_mul(): different tabdata lengths?\n");
00186 
00187     op = out->y;
00188     i1 = in1->y;
00189     i2 = in2->y;
00190     for (j = in1->ny; j > 0; j--)
00191         *op++ = *i1++ * *i2++;
00192     /* VELMUL2N( out->y, i1->y, i2->y, in1->ny ); */
00193 }
00194 
00195 /*
00196  *                      B N _ T A B D A T A _ M U L 3
00197  *@brief
00198  *  Element-by-element multiply the values from three data tables.
00199  */
00200 void
00201 bn_tabdata_mul3(struct bn_tabdata *out, const struct bn_tabdata *in1, const struct bn_tabdata *in2, const struct bn_tabdata *in3)
00202 {
00203     register int                j;
00204     register fastf_t    *op;
00205     register const fastf_t      *i1, *i2, *i3;
00206 
00207     if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_tabdata_mul3(%p, %p, %p, %p)\n", (void *)out, (void *)in1, (void *)in2, (void *)in3);
00208 
00209     BN_CK_TABDATA(out);
00210     BN_CK_TABDATA(in1);
00211     BN_CK_TABDATA(in2);
00212     BN_CK_TABDATA(in3);
00213 
00214     if (in1->table != in2->table || in1->table != out->table || in1->table != in2->table)
00215         bu_bomb("bn_tabdata_mul(): samples drawn from different tables\n");
00216     if (in1->ny != in2->ny || in1->ny != out->ny)
00217         bu_bomb("bn_tabdata_mul(): different tabdata lengths?\n");
00218 
00219     op = out->y;
00220     i1 = in1->y;
00221     i2 = in2->y;
00222     i3 = in3->y;
00223     for (j = in1->ny; j > 0; j--)
00224         *op++ = *i1++ * *i2++ * *i3++;
00225     /* VELMUL3N( out->y, i1->y, i2->y, i3->y, in1->ny ); */
00226 }
00227 
00228 /*
00229  *                      B N _ T A B D A T A _ I N C R _ M U L 3 _ S C A L E
00230  *@brief
00231  *  Element-by-element multiply the values from three data tables and a scalor.
00232  *
00233  *      out += in1 * in2 * in3 * scale
00234  */
00235 void
00236 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)
00237 {
00238     register int                j;
00239     register fastf_t    *op;
00240     register const fastf_t      *i1, *i2, *i3;
00241 
00242     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);
00243 
00244     BN_CK_TABDATA(out);
00245     BN_CK_TABDATA(in1);
00246     BN_CK_TABDATA(in2);
00247     BN_CK_TABDATA(in3);
00248 
00249     if (in1->table != in2->table || in1->table != out->table || in1->table != in3->table)
00250         bu_bomb("bn_tabdata_mul(): samples drawn from different tables\n");
00251     if (in1->ny != in2->ny || in1->ny != out->ny)
00252         bu_bomb("bn_tabdata_mul(): different tabdata lengths?\n");
00253 
00254     op = out->y;
00255     i1 = in1->y;
00256     i2 = in2->y;
00257     i3 = in3->y;
00258     for (j = in1->ny; j > 0; j--)
00259         *op++ += *i1++ * *i2++ * *i3++ * scale;
00260 }
00261 
00262 /*
00263  *                      B N _ T A B D A T A _ I N C R _ M U L 2 _ S C A L E
00264  *@brief
00265  *  Element-by-element multiply the values from two data tables and a scalor.
00266  *
00267  *      out += in1 * in2 * scale
00268  */
00269 void
00270 bn_tabdata_incr_mul2_scale(struct bn_tabdata *out, const struct bn_tabdata *in1, const struct bn_tabdata *in2, register double scale)
00271 {
00272     register int                j;
00273     register fastf_t    *op;
00274     register const fastf_t      *i1, *i2;
00275 
00276     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);
00277 
00278     BN_CK_TABDATA(out);
00279     BN_CK_TABDATA(in1);
00280     BN_CK_TABDATA(in2);
00281 
00282     if (in1->table != in2->table || in1->table != out->table)
00283         bu_bomb("bn_tabdata_mul(): samples drawn from different tables\n");
00284     if (in1->ny != in2->ny || in1->ny != out->ny)
00285         bu_bomb("bn_tabdata_mul(): different tabdata lengths?\n");
00286 
00287     op = out->y;
00288     i1 = in1->y;
00289     i2 = in2->y;
00290     for (j = in1->ny; j > 0; j--)
00291         *op++ += *i1++ * *i2++ * scale;
00292 }
00293 
00294 /*
00295  *                      B N _ T A B D A T A _ S C A L E
00296  *@brief
00297  *  Multiply every element in a data table by a scalar value 'scale'.
00298  */
00299 void
00300 bn_tabdata_scale(struct bn_tabdata *out, const struct bn_tabdata *in1, register double scale)
00301 {
00302     register int                j;
00303     register fastf_t    *op;
00304     register const fastf_t      *i1;
00305 
00306     if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_tabdata_scale(%p, %p, %g)\n", (void *)out, (void *)in1, scale);
00307 
00308     BN_CK_TABDATA(out);
00309     BN_CK_TABDATA(in1);
00310 
00311     if (in1->table != out->table)
00312         bu_bomb("bn_tabdata_scale(): samples drawn from different tables\n");
00313     if (in1->ny != out->ny)
00314         bu_bomb("bn_tabdata_scale(): different tabdata lengths?\n");
00315 
00316     op = out->y;
00317     i1 = in1->y;
00318     for (j = in1->ny; j > 0; j--)
00319         *op++ = *i1++ * scale;
00320     /* VSCALEN( out->y, in->y, scale ); */
00321 }
00322 
00323 /*
00324  *                      B N _ T A B L E _ S C A L E
00325  *@brief
00326  *  Scale the indepentent axis of a table by 'scale'.
00327  */
00328 void
00329 bn_table_scale(struct bn_table *tabp, register double scale)
00330 {
00331     register int                j;
00332     register fastf_t    *op;
00333 
00334     if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_table_scale(%p, %g)\n", (void *)tabp, scale);
00335 
00336     BN_CK_TABLE(tabp);
00337 
00338     op = tabp->x;
00339     for (j = tabp->nx+1; j > 0; j--)
00340         *op++ *= scale;
00341     /* VSCALEN( tabp->x, tabp->x, scale, tabp->nx+1 ); */
00342 }
00343 
00344 /*
00345  *                      B N _ T A B D A T A _ J O I N 1
00346  *@brief
00347  *  Multiply every element in data table in2 by a scalar value 'scale',
00348  *  add it to the element in in1, and store in 'out'.
00349  *  'out' may overlap in1 or in2.
00350  */
00351 void
00352 bn_tabdata_join1(struct bn_tabdata *out, const struct bn_tabdata *in1, register double scale, const struct bn_tabdata *in2)
00353 {
00354     register int                j;
00355     register fastf_t    *op;
00356     register const fastf_t      *i1, *i2;
00357 
00358     if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_tabdata_join1(%p, %p, %g, %p)\n", (void *)out, (void *)in1, scale, (void *)in2);
00359 
00360     BN_CK_TABDATA(out);
00361     BN_CK_TABDATA(in1);
00362     BN_CK_TABDATA(in2);
00363 
00364     if (in1->table != out->table)
00365         bu_bomb("bn_tabdata_join1(): samples drawn from different tables\n");
00366     if (in1->table != in2->table)
00367         bu_bomb("bn_tabdata_join1(): samples drawn from different tables\n");
00368     if (in1->ny != out->ny)
00369         bu_bomb("bn_tabdata_join1(): different tabdata lengths?\n");
00370 
00371     op = out->y;
00372     i1 = in1->y;
00373     i2 = in2->y;
00374     for (j = in1->ny; j > 0; j--)
00375         *op++ = *i1++ + scale * *i2++;
00376     /* VJOIN1N( out->y, in1->y, scale, in2->y ); */
00377 }
00378 
00379 /*
00380  *                      B N _ T A B D A T A _ J O I N 2
00381  *@brief
00382  *  Multiply every element in data table in2 by a scalar value 'scale2',
00383  *  plus in3 * scale3, and
00384  *  add it to the element in in1, and store in 'out'.
00385  *  'out' may overlap in1 or in2.
00386  */
00387 void
00388 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)
00389 {
00390     register int                j;
00391     register fastf_t    *op;
00392     register const fastf_t      *i1, *i2, *i3;
00393 
00394     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);
00395 
00396     BN_CK_TABDATA(out);
00397     BN_CK_TABDATA(in1);
00398     BN_CK_TABDATA(in2);
00399 
00400     if (in1->table != out->table)
00401         bu_bomb("bn_tabdata_join1(): samples drawn from different tables\n");
00402     if (in1->table != in2->table)
00403         bu_bomb("bn_tabdata_join1(): samples drawn from different tables\n");
00404     if (in1->table != in3->table)
00405         bu_bomb("bn_tabdata_join1(): samples drawn from different tables\n");
00406     if (in1->ny != out->ny)
00407         bu_bomb("bn_tabdata_join1(): different tabdata lengths?\n");
00408 
00409     op = out->y;
00410     i1 = in1->y;
00411     i2 = in2->y;
00412     i3 = in3->y;
00413     for (j = in1->ny; j > 0; j--)
00414         *op++ = *i1++ + scale2 * *i2++ + scale3 * *i3++;
00415     /* VJOIN2N( out->y, in1->y, scale2, in2->y, scale3, in3->y ); */
00416 }
00417 
00418 /*
00419  *                      B N _ T A B D A T A _ B L E N D 2
00420  */
00421 void
00422 bn_tabdata_blend2(struct bn_tabdata *out, register double scale1, const struct bn_tabdata *in1, register double scale2, const struct bn_tabdata *in2)
00423 {
00424     register int                j;
00425     register fastf_t    *op;
00426     register const fastf_t      *i1, *i2;
00427 
00428     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);
00429 
00430     BN_CK_TABDATA(out);
00431     BN_CK_TABDATA(in1);
00432     BN_CK_TABDATA(in2);
00433 
00434     if (in1->table != out->table)
00435         bu_bomb("bn_tabdata_blend2(): samples drawn from different tables\n");
00436     if (in1->table != in2->table)
00437         bu_bomb("bn_tabdata_blend2(): samples drawn from different tables\n");
00438     if (in1->ny != out->ny)
00439         bu_bomb("bn_tabdata_blend2(): different tabdata lengths?\n");
00440 
00441     op = out->y;
00442     i1 = in1->y;
00443     i2 = in2->y;
00444     for (j = in1->ny; j > 0; j--)
00445         *op++ = scale1 * *i1++ + scale2 * *i2++;
00446     /* VBLEND2N( out->y, scale1, in1->y, scale2, in2->y ); */
00447 }
00448 
00449 /*
00450  *                      B N _ T A B D A T A _ B L E N D 3
00451  */
00452 void
00453 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)
00454 {
00455     register int                j;
00456     register fastf_t    *op;
00457     register const fastf_t      *i1, *i2, *i3;
00458 
00459     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);
00460 
00461     BN_CK_TABDATA(out);
00462     BN_CK_TABDATA(in1);
00463     BN_CK_TABDATA(in2);
00464     BN_CK_TABDATA(in3);
00465 
00466     if (in1->table != out->table)
00467         bu_bomb("bn_tabdata_blend3(): samples drawn from different tables\n");
00468     if (in1->table != in2->table)
00469         bu_bomb("bn_tabdata_blend3(): samples drawn from different tables\n");
00470     if (in1->table != in3->table)
00471         bu_bomb("bn_tabdata_blend3(): samples drawn from different tables\n");
00472     if (in1->ny != out->ny)
00473         bu_bomb("bn_tabdata_blend3(): different tabdata lengths?\n");
00474 
00475     op = out->y;
00476     i1 = in1->y;
00477     i2 = in2->y;
00478     i3 = in3->y;
00479     for (j = in1->ny; j > 0; j--)
00480         *op++ = scale1 * *i1++ + scale2 * *i2++ + scale3 * *i3++;
00481     /* VBLEND3N( out->y, scale1, in1->y, scale2, in2->y, scale3, in3->y ); */
00482 }
00483 
00484 /*
00485  *                      B N _ T A B D A T A _ A R E A 1
00486  *@brief
00487  *  Following interpretation #1, where y[j] stores the total (integral
00488  *  or area) value within the interval, return the area under the whole curve.
00489  *  This is simply totaling up the areas from each of the intervals.
00490  */
00491 double
00492 bn_tabdata_area1(const struct bn_tabdata *in)
00493 {
00494     register fastf_t            area;
00495     register const fastf_t      *ip;
00496     register int                j;
00497 
00498     BN_CK_TABDATA(in);
00499 
00500     area = 0;
00501     ip = in->y;
00502     for (j = in->ny; j > 0; j--)
00503         area += *ip++;
00504 
00505     if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_tabdata_area(%p) = %g\n", (void *)in, area);
00506 
00507     return area;
00508 }
00509 
00510 /*
00511  *                      B N _ T A B D A T A _ A R E A 2
00512  *@brief
00513  *  Following interpretation #2, where y[j] stores the average
00514  *  value for the interval, return the area under
00515  *  the whole curve.  Since the iterval spacing need not be uniform,
00516  *  sum the areas of the rectangles.
00517  */
00518 double
00519 bn_tabdata_area2(const struct bn_tabdata *in)
00520 {
00521     const struct bn_table       *tabp;
00522     register fastf_t            area;
00523     fastf_t                     width;
00524     register int                j;
00525 
00526     BN_CK_TABDATA(in);
00527     tabp = in->table;
00528     BN_CK_TABLE(tabp);
00529 
00530     area = 0;
00531     for (j = in->ny-1; j >= 0; j--)  {
00532         width = tabp->x[j+1] - tabp->x[j];
00533         area += in->y[j] * width;
00534     }
00535 
00536     if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_tabdata_area2(%p) = %g\n", (void *)in, area);
00537     return area;
00538 }
00539 
00540 /*
00541  *                      B N _ T A B D A T A _ M U L _ A R E A 1
00542  *@brief
00543  *  Following interpretation #1, where y[j] stores the total (integral
00544  *  or area) value within the interval, return the area under the whole curve.
00545  *  This is simply totaling up the areas from each of the intervals.
00546  *  The curve value is found by multiplying corresponding entries from
00547  *  in1 and in2.
00548  */
00549 double
00550 bn_tabdata_mul_area1(const struct bn_tabdata *in1, const struct bn_tabdata *in2)
00551 {
00552     register fastf_t            area;
00553     register const fastf_t      *i1, *i2;
00554     register int                j;
00555 
00556     BN_CK_TABDATA(in1);
00557     BN_CK_TABDATA(in2);
00558 
00559     area = 0;
00560     i1 = in1->y;
00561     i2 = in2->y;
00562     for (j = in1->ny; j > 0; j--)
00563         area += *i1++ * *i2++;
00564 
00565     if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_tabdata_mul_area1(%p, %p) = %g\n", (void *)in1, (void *)in2, area);
00566     return area;
00567 }
00568 
00569 /*
00570  *                      B N _ T A B D A T A _ M U L _ A R E A 2
00571  *@brief
00572  *  Following interpretation #2,
00573  *  return the area under the whole curve.
00574  *  The curve value is found by multiplying corresponding entries from
00575  *  in1 and in2.
00576  */
00577 double
00578 bn_tabdata_mul_area2(const struct bn_tabdata *in1, const struct bn_tabdata *in2)
00579 {
00580     const struct bn_table       *tabp;
00581     register fastf_t            area;
00582     fastf_t                     width;
00583     register int                j;
00584 
00585     BN_CK_TABDATA(in1);
00586     BN_CK_TABDATA(in2);
00587     tabp = in1->table;
00588     BN_CK_TABLE(tabp);
00589 
00590     area = 0;
00591     for (j = in1->ny-1; j >= 0; j--)  {
00592         width = tabp->x[j+1] - tabp->x[j];
00593         area += in1->y[j] * in2->y[j] * width;
00594     }
00595 
00596     if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_tabdata_mul_area2(%p, %p) = %g\n", (void *)in1, (void *)in2, area);
00597     return area;
00598 }
00599 
00600 /*
00601  *                      B N _ T A B L E _ F I N D _ X
00602  *@brief
00603  *  Return the index in the table's x[] array of the interval which
00604  *  contains 'xval'.
00605  *
00606  *
00607  * @return      -1      if less than start value
00608  * @return      -2      if greater than end value.
00609  * @return      0..nx-1 otherwise.
00610  *              Never returns a value of nx.
00611  *
00612  *  A binary search would be more efficient, as the wavelengths (x values)
00613  *  are known to be sorted in ascending order.
00614  */
00615 long
00616 bn_table_find_x(const struct bn_table *tabp, double xval)
00617 {
00618     register int        i;
00619 
00620     if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_table_find_x(%p, %g)\n", (void *)tabp, xval);
00621     BN_CK_TABLE(tabp);
00622 
00623     if (xval > tabp->x[tabp->nx])  return -2;
00624     if (xval >= tabp->x[tabp->nx-1])  return tabp->nx-1;
00625 
00626     /* Search for proper interval in input spectrum */
00627     for (i = tabp->nx-2; i >=0; i--)  {
00628         if (xval >= tabp->x[i])  return i;
00629     }
00630     /* if ( xval < tabp->x[0] )  return -1; */
00631     return -1;
00632 }
00633 
00634 /*
00635  *                      B N _ T A B L E _ L I N _ I N T E R P
00636  *@brief
00637  *  Return the value of the curve at independent parameter value 'wl'.
00638  *  Linearly interpolate between values in the input table.
00639  *  Zero is returned for values outside the sampled range.
00640  */
00641 fastf_t
00642 bn_table_lin_interp(const struct bn_tabdata *samp, register double wl)
00643 {
00644     const struct bn_table *tabp;
00645     size_t i;
00646     int idx;
00647     fastf_t fract;
00648     fastf_t ret;
00649 
00650     if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_table_lin_interp(%p, %g)\n", (void *)samp, wl);
00651 
00652     BN_CK_TABDATA(samp);
00653     tabp = samp->table;
00654     BN_CK_TABLE(tabp);
00655 
00656     idx = bn_table_find_x(tabp, wl);
00657     if (idx < 0)  {
00658         if (bu_debug&BU_DEBUG_TABDATA)
00659             bu_log("bn_table_lin_interp(%g) out of range %g to %g\n", wl, tabp->x[0], tabp->x[tabp->nx]);
00660         return 0;
00661     }
00662     i = (size_t)idx;
00663 
00664     if (wl < tabp->x[i] || wl >= tabp->x[i+1])  {
00665         bu_log("bn_table_lin_interp(%g) assertion1 failed at %g\n", wl, tabp->x[i]);
00666         bu_bomb("bn_table_lin_interp() assertion1 failed\n");
00667     }
00668 
00669     if (i >= tabp->nx-2)  {
00670         /* Assume value is constant in final interval. */
00671         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]);
00672         return samp->y[tabp->nx-1];
00673     }
00674 
00675     /* The interval has been found */
00676     fract = (wl - tabp->x[i]) / (tabp->x[i+1] - tabp->x[i]);
00677     if (fract < 0 || fract > 1)  bu_bomb("bn_table_lin_interp() assertion2 failed\n");
00678     ret = (1-fract) * samp->y[i] + fract * samp->y[i+1];
00679     if (bu_debug&BU_DEBUG_TABDATA)bu_log("bn_table_lin_interp(%g)=%g in range %g to %g\n",
00680                                              wl, ret, tabp->x[i], tabp->x[i+1]);
00681     return ret;
00682 }
00683 
00684 /*
00685  *                      B N _ T A B D A T A _ R E S A M P L E _ M A X
00686  *@brief
00687  *  Given a set of sampled data 'olddata', resample it for different
00688  *  spacing, by linearly interpolating the values when an output span
00689  *  is entirely contained within an input span, and by taking the
00690  *  maximum when an output span covers more than one input span.
00691  *
00692  *  This assumes interpretation (2) of the data, i.e. that the values
00693  *  are the average value across the interval.
00694  */
00695 struct bn_tabdata *
00696 bn_tabdata_resample_max(const struct bn_table *newtable, const struct bn_tabdata *olddata) {
00697     const struct bn_table       *oldtable;
00698     struct bn_tabdata   *newsamp;
00699     size_t i;
00700     int                 j, k;
00701 
00702     if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_tabdata_resample_max(%p, %p)\n", (void *)newtable, (void *)olddata);
00703 
00704     BN_CK_TABLE(newtable);
00705     BN_CK_TABDATA(olddata);
00706     oldtable = olddata->table;
00707     BN_CK_TABLE(oldtable);
00708 
00709     if (oldtable == newtable)  bu_log("bn_tabdata_resample_max() NOTICE old and new bn_table structs are the same\n");
00710 
00711     BN_GET_TABDATA(newsamp, newtable);
00712 
00713     for (i = 0; i < newtable->nx; i++)  {
00714         /*
00715          *  Find good value(s) in olddata to represent the span from
00716          *  newtable->x[i] to newtable->x[i+1].
00717          */
00718         j = bn_table_find_x(oldtable, newtable->x[i]);
00719         k = bn_table_find_x(oldtable, newtable->x[i+1]);
00720         if (k == -1)  {
00721             /* whole new span is off left side of old table */
00722             newsamp->y[i] = 0;
00723             continue;
00724         }
00725         if (j == -2)  {
00726             /* whole new span is off right side of old table */
00727             newsamp->y[i] = 0;
00728             continue;
00729         }
00730 
00731         if (j == k && j > 0)  {
00732             register fastf_t tmp;
00733             /*
00734              *  Simple case, ends of output span are completely
00735              *  contained within one input span.
00736              *  Interpolate for both ends, take max.
00737              *  XXX this could be more efficiently written inline here.
00738              */
00739             newsamp->y[i] = bn_table_lin_interp(olddata, newtable->x[i]);
00740             tmp = bn_table_lin_interp(olddata, newtable->x[i+1]);
00741             if (tmp > newsamp->y[i])  newsamp->y[i] = tmp;
00742         } else {
00743             register fastf_t tmp, n;
00744             register int        s;
00745             /*
00746              *  Complex case: find good representative value.
00747              *  Interpolate both ends, and consider all
00748              *  intermediate old samples in span.  Take max.
00749              *  One (but not both) new ends may be off old table.
00750              */
00751             n = bn_table_lin_interp(olddata, newtable->x[i]);
00752             tmp = bn_table_lin_interp(olddata, newtable->x[i+1]);
00753             if (tmp > n)  n = tmp;
00754             for (s = j+1; s <= k; s++)  {
00755                 if ((tmp = olddata->y[s]) > n)
00756                     n = tmp;
00757             }
00758             newsamp->y[i] = n;
00759         }
00760     }
00761     return newsamp;
00762 }
00763 
00764 /*
00765  *                      B N _ T A B D A T A _ R E S A M P L E _ A V G
00766  *@brief
00767  *  Given a set of sampled data 'olddata', resample it for different
00768  *  spacing, by linearly interpolating the values when an output span
00769  *  is entirely contained within an input span, and by taking the
00770  *  average when an output span covers more than one input span.
00771  *
00772  *  This assumes interpretation (2) of the data, i.e. that the values
00773  *  are the average value across the interval.
00774  */
00775 struct bn_tabdata *
00776 bn_tabdata_resample_avg(const struct bn_table *newtable, const struct bn_tabdata *olddata) {
00777     const struct bn_table       *oldtable;
00778     struct bn_tabdata   *newsamp;
00779     size_t i;
00780     int                 j, k;
00781 
00782     if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_tabdata_resample_avg(%p, %p)\n", (void *)newtable, (void *)olddata);
00783 
00784     BN_CK_TABLE(newtable);
00785     BN_CK_TABDATA(olddata);
00786     oldtable = olddata->table;
00787     BN_CK_TABLE(oldtable);
00788 
00789     if (oldtable == newtable)  bu_log("bn_tabdata_resample_avg() NOTICE old and new bn_table structs are the same\n");
00790 
00791     BN_GET_TABDATA(newsamp, newtable);
00792 
00793     for (i = 0; i < newtable->nx; i++)  {
00794         /*
00795          *  Find good value(s) in olddata to represent the span from
00796          *  newtable->x[i] to newtable->x[i+1].
00797          */
00798         j = bn_table_find_x(oldtable, newtable->x[i]);
00799         k = bn_table_find_x(oldtable, newtable->x[i+1]);
00800 
00801         if (j < 0 || k < 0 || j == k)  {
00802             /*
00803              *  Simple case, ends of output span are completely
00804              *  contained within one input span.
00805              *  Interpolate for both ends, take average.
00806              *  XXX this could be more efficiently written inline here.
00807              */
00808             newsamp->y[i] = 0.5 * (
00809                                 bn_table_lin_interp(olddata, newtable->x[i]) +
00810                                 bn_table_lin_interp(olddata, newtable->x[i+1]));
00811         } else {
00812             /*
00813              *  Complex case: find average value.
00814              *  Interpolate both end, and consider all
00815              *  intermediate old spans.
00816              *  There are three parts to sum:
00817              *  Partial interval from newx[i] to j+1
00818              *  Full intervals from j+1 to k
00819              *  Partial interval from k to newx[i+1]
00820              */
00821             fastf_t wsum;               /* weighted sum */
00822             fastf_t     a, b;           /* values being averaged */
00823             int s;
00824 
00825             /* Partial interval from newx[i] to j+1 */
00826             a = bn_table_lin_interp(olddata, newtable->x[i]);   /* in "j" bin */
00827             b = olddata->y[j+1];
00828             wsum = 0.5 * (a+b) * (oldtable->x[j+1] - newtable->x[i]);
00829 
00830             /* Full intervals from j+1 to k */
00831             for (s = j+1; s < k; s++)  {
00832                 a = olddata->y[s];
00833                 b = olddata->y[s+1];
00834                 wsum += 0.5 * (a+b) * (oldtable->x[s+1] - oldtable->x[s]);
00835             }
00836 
00837             /* Partial interval from k to newx[i+1] */
00838             a = olddata->y[k];
00839             b = bn_table_lin_interp(olddata, newtable->x[i+1]); /* in "k" bin */
00840             wsum += 0.5 * (a+b) * (newtable->x[i+1] - oldtable->x[k]);
00841 
00842             /* Adjust the weighted sum by the total width */
00843             newsamp->y[i] =
00844                 wsum / (newtable->x[i+1] - newtable->x[i]);
00845         }
00846     }
00847     return newsamp;
00848 }
00849 
00850 /*
00851  *                      B N _ T A B L E _ W R I T E
00852  *@brief
00853  *  Write out the table structure in an ASCII file,
00854  *  giving the number of values (minus 1), and the
00855  *  actual values.
00856  */
00857 int
00858 bn_table_write(const char *filename, const struct bn_table *tabp)
00859 {
00860     FILE *fp;
00861     size_t j;
00862 
00863     if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_table_write(%s, %p)\n", filename, (void *)tabp);
00864 
00865     BN_CK_TABLE(tabp);
00866 
00867     bu_semaphore_acquire(BU_SEM_SYSCALL);
00868     fp = fopen(filename, "wb");
00869     bu_semaphore_release(BU_SEM_SYSCALL);
00870 
00871     if (fp == NULL)  {
00872         perror(filename);
00873         bu_log("bn_table_write(%s, %p) FAILED\n", filename, (void *)tabp);
00874         return -1;
00875     }
00876 
00877     bu_semaphore_acquire(BU_SEM_SYSCALL);
00878     fprintf(fp, "  %lu sample starts, and one end.\n", (long unsigned int)tabp->nx);
00879     for (j=0; j <= tabp->nx; j++)  {
00880         fprintf(fp, "%g\n", tabp->x[j]);
00881     }
00882     fclose(fp);
00883     bu_semaphore_release(BU_SEM_SYSCALL);
00884     return 0;
00885 }
00886 
00887 /*
00888  *                      B N _ T A B L E _ R E A D
00889  *@brief
00890  *  Allocate and read in the independent variable values from an ASCII file,
00891  *  giving the number of samples (minus 1), and the
00892  *  actual values.
00893  */
00894 struct bn_table *
00895 bn_table_read(const char *filename) {
00896     int ret;
00897     struct bn_table *tabp;
00898     struct bu_vls line = BU_VLS_INIT_ZERO;
00899     FILE *fp;
00900     size_t nw;
00901     size_t j;
00902 
00903     if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_table_read(%s)\n", filename);
00904 
00905     bu_semaphore_acquire(BU_SEM_SYSCALL);
00906     fp = fopen(filename, "rb");
00907     bu_semaphore_release(BU_SEM_SYSCALL);
00908 
00909     if (fp == NULL)  {
00910         perror(filename);
00911         bu_log("bn_table_read(%s) FAILED\n", filename);
00912         return NULL;
00913     }
00914 
00915     if (bu_vls_gets(&line, fp) < 0) {
00916         perror(filename);
00917         bu_log("Failed to read line\n");
00918         fclose(fp);
00919         return NULL;
00920     }
00921     nw = 0;
00922     /* TODO: %lu to size_t isn't right. We may need a bu_sscanf() that can cope
00923      * with native pointer width correctly, as %p is not ubiquitous and %z is a
00924      * C99-ism */
00925     sscanf(bu_vls_addr(&line), "%lu", (long unsigned int *)(&nw));
00926     bu_vls_free(&line);
00927 
00928     if (nw <= 0) bu_bomb("bn_table_read() bad nw value\n");
00929 
00930     BN_GET_TABLE(tabp, nw);
00931 
00932     bu_semaphore_acquire(BU_SEM_SYSCALL);
00933     for (j=0; j <= tabp->nx; j++)  {
00934         /* XXX assumes fastf_t == double */
00935         ret = fscanf(fp, "%lf", &tabp->x[j]);
00936         if (ret != 1) {
00937             bu_log("bn_table_read(%s) READ FAILURE. Abort\n", filename);
00938             break;
00939         }
00940     }
00941     fclose(fp);
00942     bu_semaphore_release(BU_SEM_SYSCALL);
00943 
00944     bn_ck_table(tabp);
00945 
00946     return tabp;
00947 }
00948 
00949 /*
00950  *                      B N _ P R _ T A B L E
00951  */
00952 void
00953 bn_pr_table(const char *title, const struct bn_table *tabp)
00954 {
00955     size_t j;
00956 
00957     bu_log("%s\n", title);
00958     BN_CK_TABLE(tabp);
00959 
00960     for (j=0; j <= tabp->nx; j++)  {
00961         bu_log("%3zu: %g\n", j, tabp->x[j]);
00962     }
00963 }
00964 
00965 /*
00966  *                      B N _ P R _ T A B D A T A
00967  */
00968 void
00969 bn_pr_tabdata(const char *title, const struct bn_tabdata *data)
00970 {
00971     size_t j;
00972 
00973     bu_log("%s: ", title);
00974     BN_CK_TABDATA(data);
00975 
00976     for (j=0; j < data->ny; j++)  {
00977         bu_log("%g, ", data->y[j]);
00978     }
00979     bu_log("\n");
00980 }
00981 
00982 /*
00983  *                      B N _ P R I N T _ T A B L E _ A N D _ T A B D A T A
00984  *@brief
00985  *  Write out a given data table into an ASCII file,
00986  *  suitable for input to GNUPLOT.
00987  *
00988  *      (set term postscript)
00989  *      (set output "|print-postscript")
00990  *      (plot "filename" with lines)
00991  */
00992 int
00993 bn_print_table_and_tabdata(const char *filename, const struct bn_tabdata *data)
00994 {
00995     FILE        *fp;
00996     const struct bn_table       *tabp;
00997     size_t j;
00998 
00999     if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_print_table_and_tabdata(%s, %p)\n", filename, (void *)data);
01000 
01001     BN_CK_TABDATA(data);
01002     tabp = data->table;
01003     BN_CK_TABLE(tabp);
01004 
01005     bu_semaphore_acquire(BU_SEM_SYSCALL);
01006     fp = fopen(filename, "wb");
01007     bu_semaphore_release(BU_SEM_SYSCALL);
01008 
01009     if (fp == NULL)  {
01010         perror(filename);
01011         bu_log("bn_print_table_and_tabdata(%s, %p) FAILED\n", filename, (void *)data);
01012         return -1;
01013     }
01014 
01015     bu_semaphore_acquire(BU_SEM_SYSCALL);
01016     for (j=0; j < tabp->nx; j++)  {
01017         fprintf(fp, "%g %g\n", tabp->x[j], data->y[j]);
01018     }
01019     fprintf(fp, "%g (novalue)\n", tabp->x[tabp->nx]);
01020     fclose(fp);
01021     bu_semaphore_release(BU_SEM_SYSCALL);
01022     return 0;
01023 }
01024 
01025 /*
01026  *                      B N _ R E A D _ T A B L E _ A N D _ T A B D A T A
01027  *@brief
01028  *  Read in a file which contains two columns of numbers, the first
01029  *  column being the wavelength, the second column being the sample value
01030  *  at that wavelength.
01031  *
01032  *  A new bn_table structure and one bn_tabdata structure
01033  *  are created, a pointer to the bn_tabdata structure is returned.
01034  *  The final wavelength is guessed at.
01035  */
01036 struct bn_tabdata *
01037 bn_read_table_and_tabdata(const char *filename) {
01038     struct bn_table     *tabp;
01039     struct bn_tabdata   *data;
01040     FILE        *fp;
01041     char        buf[128];
01042     size_t count = 0;
01043     size_t i;
01044 
01045     if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_read_table_and_tabdata(%s)\n", filename);
01046 
01047     bu_semaphore_acquire(BU_SEM_SYSCALL);
01048     fp = fopen(filename, "rb");
01049     bu_semaphore_release(BU_SEM_SYSCALL);
01050 
01051     if (fp == NULL)  {
01052         perror(filename);
01053         bu_log("bn_read_table_and_tabdata(%s) FAILED. Couldn't open file.\n", filename);
01054         return NULL;
01055     }
01056 
01057     /* First pass:  Count number of lines */
01058     bu_semaphore_acquire(BU_SEM_SYSCALL);
01059     for (;;)  {
01060         if (bu_fgets(buf, sizeof(buf), fp) == NULL)  break;
01061         count++;
01062     }
01063     fclose(fp);
01064     bu_semaphore_release(BU_SEM_SYSCALL);
01065 
01066     if (count < 2)  {
01067         perror(filename);
01068         bu_log("bn_read_table_and_tabdata(%s) FAILED. Expected at least 2 lines in file.\n", filename);
01069         return NULL;
01070     }
01071 
01072     /* Allocate storage */
01073     BN_GET_TABLE(tabp, count);
01074     BN_GET_TABDATA(data, tabp);
01075 
01076     /* Second pass:  Read only as much data as storage was allocated for */
01077     bu_semaphore_acquire(BU_SEM_SYSCALL);
01078     fp = fopen(filename, "rb");
01079     for (i=0; i < count; i++) {
01080         buf[0] = '\0';
01081         if (bu_fgets(buf, sizeof(buf), fp) == NULL)  {
01082             bu_log("bn_read_table_and_tabdata(%s) unexpected EOF on line %zu\n", filename, i);
01083             break;
01084         }
01085         sscanf(buf, "%lf %lf", &tabp->x[i], &data->y[i]);
01086     }
01087     fclose(fp);
01088     bu_semaphore_release(BU_SEM_SYSCALL);
01089 
01090     /* Complete final interval */
01091     tabp->x[count] = 2 * tabp->x[count-1] - tabp->x[count-2];
01092 
01093     bn_ck_table(tabp);
01094 
01095     return data;
01096 }
01097 
01098 /*
01099  *                      B N _ T A B D A T A _ B I N A R Y _ R E A D
01100  */
01101 struct bn_tabdata *
01102 bn_tabdata_binary_read(const char *filename, size_t num, const struct bn_table *tabp) {
01103     struct bn_tabdata   *data;
01104     char        *cp;
01105     int nbytes;
01106     int len;
01107     int got;
01108     int fd;
01109     long i;
01110 
01111     if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_tabdata_binary_read(%s, num=%zu, %p)\n", filename, num, (void *)tabp);
01112 
01113     BN_CK_TABLE(tabp);
01114 
01115     nbytes = BN_SIZEOF_TABDATA(tabp);
01116     len = num * nbytes;
01117 
01118     bu_semaphore_acquire(BU_SEM_SYSCALL);
01119     fd = open(filename, 0);
01120     bu_semaphore_release(BU_SEM_SYSCALL);
01121     if (fd <= 0)  {
01122         perror(filename);
01123         bu_log("bn_tabdata_binary_read open failed on \"%s\"\n", filename);
01124         return (struct bn_tabdata *)NULL;
01125     }
01126 
01127     /* Get a big array of structures for reading all at once */
01128     data = (struct bn_tabdata *)bu_malloc(len+8, "bn_tabdata[]");
01129 
01130     bu_semaphore_acquire(BU_SEM_SYSCALL);
01131     got = read(fd, (char *)data, len);
01132     bu_semaphore_release(BU_SEM_SYSCALL);
01133     if (got != len)  {
01134         if (got < 0) {
01135             perror(filename);
01136             bu_log("bn_tabdata_binary_read read error on \"%s\"\n", filename);
01137         } else {
01138             bu_log("bn_tabdata_binary_read(%s) expected %d got %d\n", filename, len, got);
01139         }
01140         bu_free(data, "bn_tabdata[]");
01141         bu_semaphore_acquire(BU_SEM_SYSCALL);
01142         close(fd);
01143         bu_semaphore_release(BU_SEM_SYSCALL);
01144         return (struct bn_tabdata *)NULL;
01145     }
01146     bu_semaphore_acquire(BU_SEM_SYSCALL);
01147     close(fd);
01148     bu_semaphore_release(BU_SEM_SYSCALL);
01149 
01150     /* Connect data[i].table pointer to tabp */
01151     cp = (char *)data;
01152     for (i = num-1; i >= 0; i--, cp += nbytes)  {
01153         register struct bn_tabdata *sp;
01154         sp = (struct bn_tabdata *)cp;
01155         BN_CK_TABDATA(sp);
01156         sp->table = tabp;
01157     }
01158 
01159     return data;
01160 }
01161 
01162 /*
01163  *                      B N _ T A B D A T A _ M A L L O C _ A R R A Y
01164  *@brief
01165  *  Allocate storage for, and initialize, an array of 'num' data table
01166  *  structures.
01167  *  This subroutine is provided because the bn_tabdata structures
01168  *  are variable length.
01169  */
01170 struct bn_tabdata *
01171 bn_tabdata_malloc_array(const struct bn_table *tabp, size_t num) {
01172     struct bn_tabdata *data;
01173     char *cp;
01174     size_t i;
01175     size_t nw;
01176     size_t nbytes;
01177 
01178     if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_tabdata_malloc_array(%p, num=%zu)\n", (void *)tabp, num);
01179 
01180     BN_CK_TABLE(tabp);
01181     nw = tabp->nx;
01182     nbytes = BN_SIZEOF_TABDATA(tabp);
01183 
01184     data = (struct bn_tabdata *)bu_calloc(num,
01185                                           nbytes, "struct bn_tabdata[]");
01186 
01187     cp = (char *)data;
01188     for (i = 0; i < num; i++) {
01189         register struct bn_tabdata      *sp;
01190 
01191         sp = (struct bn_tabdata *)cp;
01192         sp->magic = BN_TABDATA_MAGIC;
01193         sp->ny = nw;
01194         sp->table = tabp;
01195         cp += nbytes;
01196     }
01197     return data;
01198 }
01199 
01200 /*
01201  *                      B N _ T A B D A T A _ C O P Y
01202  */
01203 void
01204 bn_tabdata_copy(struct bn_tabdata *out, const struct bn_tabdata *in)
01205 {
01206     if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_tabdata_copy(%p, %p)\n", (void *)out, (void *)in);
01207 
01208     BN_CK_TABDATA(out);
01209     BN_CK_TABDATA(in);
01210 
01211     if (in->table != out->table)
01212         bu_bomb("bn_tabdata_copy(): samples drawn from different tables\n");
01213     if (in->ny != out->ny)
01214         bu_bomb("bn_tabdata_copy(): different tabdata lengths?\n");
01215 
01216     memcpy((char *)out->y, (const char *)in->y, BN_SIZEOF_TABDATA_Y(in));
01217 }
01218 
01219 /*
01220  *                      B N _ T A B D A T A _ D U P
01221  */
01222 struct bn_tabdata *
01223 bn_tabdata_dup(const struct bn_tabdata *in) {
01224     struct bn_tabdata *data;
01225 
01226     BN_CK_TABDATA(in);
01227     BN_GET_TABDATA(data, in->table);
01228 
01229     memcpy((char *)data->y, (const char *)in->y, BN_SIZEOF_TABDATA_Y(in));
01230 
01231     if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_tabdata_dup(%p) = %p\n", (void *)in, (void *)data);
01232     return data;
01233 }
01234 
01235 /*
01236  *                      B N _ T A B D A T A _ G E T _ C O N S T V A L
01237  *@brief
01238  *  For a given table, allocate and return a tabdata structure
01239  *  with all elements initialized to 'val'.
01240  */
01241 struct bn_tabdata *
01242 bn_tabdata_get_constval(double val, const struct bn_table *tabp) {
01243     struct bn_tabdata   *data;
01244     int                 todo;
01245     register fastf_t    *op;
01246 
01247     BN_CK_TABLE(tabp);
01248     BN_GET_TABDATA(data, tabp);
01249 
01250     op = data->y;
01251     for (todo = data->ny-1; todo >= 0; todo--)
01252         *op++ = val;
01253 
01254     if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_tabdata_get_constval(val=%g, %p)=%p\n", val, (void *)tabp, (void *)data);
01255 
01256     return data;
01257 }
01258 
01259 /*
01260  *                      B N _ T A B D A T A _ C O N S T V A L
01261  *@brief
01262  *  Set all the tabdata elements to 'val'
01263  */
01264 void
01265 bn_tabdata_constval(struct bn_tabdata *data, double val)
01266 {
01267     int                 todo;
01268     register fastf_t    *op;
01269 
01270     if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_tabdata_constval(%p, val=%g)\n", (void *)data, val);
01271 
01272     BN_CK_TABDATA(data);
01273 
01274     op = data->y;
01275     for (todo = data->ny-1; todo >= 0; todo--)
01276         *op++ = val;
01277 }
01278 
01279 /*
01280  *                      B N _ T A B D A T A _ T O _ T C L
01281  *@brief
01282  *  Convert an bn_tabdata/bn_table pair into a Tcl compatible string
01283  *  appended to a VLS.  It will have form:
01284  *      x {...} y {...} nx # ymin # ymax #
01285  */
01286 void
01287 bn_tabdata_to_tcl(struct bu_vls *vp, const struct bn_tabdata *data)
01288 {
01289     const struct bn_table       *tabp;
01290     register size_t i;
01291     register fastf_t    minval = MAX_FASTF, maxval = -MAX_FASTF;
01292 
01293     if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_tabdata_to_tcl(%p, %p)\n", (void *)vp, (void *)data);
01294 
01295     BU_CK_VLS(vp);
01296     BN_CK_TABDATA(data);
01297     tabp = data->table;
01298     BN_CK_TABLE(tabp);
01299 
01300     bu_vls_strcat(vp, "x {");
01301     for (i=0; i < tabp->nx; i++)  {
01302         bu_vls_printf(vp, "%g ", tabp->x[i]);
01303     }
01304     bu_vls_strcat(vp, "} y {");
01305     for (i=0; i < data->ny; i++)  {
01306         register fastf_t val = data->y[i];
01307         bu_vls_printf(vp, "%g ", val);
01308         if (val < minval)  minval = val;
01309         if (val > maxval)  maxval = val;
01310     }
01311     bu_vls_printf(vp, "} nx %zu ymin %g ymax %g",
01312                   tabp->nx, minval, maxval);
01313 }
01314 
01315 /*
01316  *                      B N _ T A B D A T A _ F R O M _ A R R A Y
01317  *@brief
01318  *  Given an array of (x, y) pairs, build the relevant bn_table and
01319  *  bn_tabdata structures.
01320  *
01321  *  The table is terminated by an x value <= 0.
01322  *  Consistent with the interpretation of the spans,
01323  *  invent a final span ending x value.
01324  */
01325 struct bn_tabdata *
01326 bn_tabdata_from_array(const double *array) {
01327     register const double *dp;
01328     size_t len = 0;
01329     struct bn_table *tabp;
01330     struct bn_tabdata *data;
01331     register size_t i;
01332 
01333     /* First, find len */
01334     for (dp = array; *dp > 0; dp += 2)
01335         ; /* NIL */
01336     len = (dp - array) >> 1;
01337 
01338     /* Second, build bn_table */
01339     BN_GET_TABLE(tabp, len);
01340     for (i = 0; i < len; i++)  {
01341         tabp->x[i] = array[i<<1];
01342     }
01343     tabp->x[len] = tabp->x[len-1] + 1;  /* invent span end */
01344 
01345     /* Third, build bn_tabdata (last input "y" is ignored) */
01346     BN_GET_TABDATA(data, tabp);
01347     for (i = 0; i < len-1; i++)  {
01348         data->y[i] = array[(i<<1)+1];
01349     }
01350 
01351     if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_tabdata_from_array(%p) = %p\n", (void *)array, (void *)data);
01352     return data;
01353 }
01354 
01355 /*
01356  *                      B N _ T A B D A T A _ F R E Q _ S H I F T
01357  *@brief
01358  *  Shift the data by a constant offset in the independent variable
01359  *  (often frequency), interpolating new sample values.
01360  */
01361 void
01362 bn_tabdata_freq_shift(struct bn_tabdata *out, const struct bn_tabdata *in, double offset)
01363 {
01364     const struct bn_table       *tabp;
01365     register size_t i;
01366 
01367     if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_tabdata_freq_shift(%p, %p, offset=%g)\n", (void *)out, (void *)in, offset);
01368 
01369     BN_CK_TABDATA(out);
01370     BN_CK_TABDATA(in);
01371     tabp = in->table;
01372 
01373     if (tabp != out->table)
01374         bu_bomb("bn_tabdata_freq_shift(): samples drawn from different tables\n");
01375     if (in->ny != out->ny)
01376         bu_bomb("bn_tabdata_freq_shift(): different tabdata lengths?\n");
01377 
01378     for (i=0; i < out->ny; i++) {
01379         out->y[i] = bn_table_lin_interp(in, tabp->x[i]+offset);
01380     }
01381 }
01382 
01383 /*
01384  *                      B N _ T A B L E _ I N T E R V A L _ N U M _ S A M P L E S
01385  *@brief
01386  *  Returns number of sample points between 'low' and 'hi', inclusive.
01387  */
01388 int
01389 bn_table_interval_num_samples(const struct bn_table *tabp, double low, double hi)
01390 {
01391     register size_t i;
01392     register size_t count = 0;
01393 
01394     BN_CK_TABLE(tabp);
01395 
01396     for (i=0; i < tabp->nx-1; i++)  {
01397         if (tabp->x[i] >= low && tabp->x[i] <= hi)  count++;
01398     }
01399 
01400     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);
01401     return count;
01402 }
01403 
01404 /*
01405  *                      B N _ T A B L E _ D E L E T E _ S A M P L E _ P T S
01406  *@brief
01407  *  Remove all sampling points between subscripts i and j, inclusive.
01408  *  Don't bother freeing the tiny bit of storage at the end of the array.
01409  *  Returns number of points removed.
01410  */
01411 int
01412 bn_table_delete_sample_pts(struct bn_table *tabp, unsigned int i, unsigned int j)
01413 {
01414     size_t tokill;
01415     unsigned int k;
01416 
01417     if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_table_delete_samples(%p, %d, %d)\n", (void *)tabp, i, j);
01418 
01419     BN_CK_TABLE(tabp);
01420 
01421     if (i >= tabp->nx || j >= tabp->nx)
01422         bu_bomb("bn_table_delete_sample_pts() index out of range\n");
01423 
01424     tokill = j - i + 1;
01425     if (tokill < 1)  bu_bomb("bn_table_delete_sample_pts(): nothing to delete\n");
01426     if (tokill >= tabp->nx) bu_bomb("bn_table_delete_sample_pts(): you can't kill 'em all!\n");
01427 
01428     tabp->nx -= tokill;
01429 
01430     for (k = i; k < tabp->nx; k++)  {
01431         tabp->x[k] = tabp->x[k+tokill];
01432     }
01433     return tokill;
01434 }
01435 
01436 /*
01437  *                      B N _ T A B L E _ M E R G E 2
01438  *@brief
01439  *  A new table is returned which has sample points at places from
01440  *  each of the input tables.
01441  */
01442 struct bn_table *
01443 bn_table_merge2(const struct bn_table *a, const struct bn_table *b) {
01444     struct bn_table *newtab;
01445     register size_t i, j, k;
01446 
01447     BN_CK_TABLE(a);
01448     BN_CK_TABLE(b);
01449 
01450     BN_GET_TABLE(newtab, a->nx + b->nx + 2);
01451 
01452     i = j = 0;          /* input subscripts */
01453     k = 0;                      /* output subscript */
01454     while (i <= a->nx || j <= b->nx)  {
01455         if (i > a->nx)  {
01456             while (j <= b->nx)
01457                 newtab->x[k++] = b->x[j++];
01458             break;
01459         }
01460         if (j > b->nx)  {
01461             while (i <= a->nx)
01462                 newtab->x[k++] = a->x[i++];
01463             break;
01464         }
01465         /* Both have remaining elements, take lower one */
01466         if (ZERO(a->x[i] - b->x[j]))  {
01467             newtab->x[k++] = a->x[i++];
01468             j++;                /* compress out duplicate */
01469             continue;
01470         }
01471         if (a->x[i] <= b->x[j])  {
01472             newtab->x[k++] = a->x[i++];
01473         } else {
01474             newtab->x[k++] = b->x[j++];
01475         }
01476     }
01477     if (k > newtab->nx)  bu_bomb("bn_table_merge2() assertion failed, k>nx?\n");
01478     newtab->nx = k-1;
01479 
01480     if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_table_merge2(%p, %p) = %p\n", (void *)a, (void *)b, (void *)newtab);
01481     return newtab;
01482 }
01483 
01484 /*
01485  *              B N _ T A B D A T A _ M K _ L I N E A R _ F I L T E R
01486  *@brief
01487  *  Create a filter to accept power in a given band.
01488  *  The first and last filter values will be in the range 0..1,
01489  *  while all the internal filter values will be 1.0,
01490  *  and all samples outside the given band will be 0.0.
01491  *
01492  *
01493  *  @return     NULL    if given band does not overlap input spectrum
01494  *  @return     tabdata*
01495  */
01496 struct bn_tabdata *
01497 bn_tabdata_mk_linear_filter(const struct bn_table *spectrum, double lower_wavelen, double upper_wavelen) {
01498     struct bn_tabdata *filt;
01499     int first;
01500     int last;
01501     fastf_t filt_range;
01502     fastf_t     cell_range;
01503     fastf_t     frac;
01504     size_t i;
01505 
01506     BN_CK_TABLE(spectrum);
01507 
01508     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);
01509 
01510     if (lower_wavelen < spectrum->x[0])
01511         if (upper_wavelen > spectrum->x[spectrum->nx])
01512             bu_log("bn_tabdata_mk_linear_filter() warning, upper_wavelen %g > hightest sampled wavelen %g\n",
01513                    upper_wavelen, spectrum->x[spectrum->nx]);
01514 
01515     /* First, find first (possibly partial) sample */
01516     first = bn_table_find_x(spectrum, lower_wavelen);
01517     if (first == -1)  {
01518         first = 0;
01519         bu_log("bn_tabdata_mk_linear_filter() warning, lower_wavelen %g < lowest sampled wavelen %g\n",
01520                lower_wavelen, spectrum->x[0]);
01521     } else if (first <= -2)  {
01522         bu_log("bn_tabdata_mk_linear_filter() ERROR, lower_wavelen %g > highest sampled wavelen %g\n",
01523                lower_wavelen, spectrum->x[spectrum->nx]);
01524         return NULL;
01525     }
01526 
01527     /* Second, find last (possibly partial) sample */
01528     last = bn_table_find_x(spectrum, upper_wavelen);
01529     if (last == -1)  {
01530         bu_log("bn_tabdata_mk_linear_filter() ERROR, upper_wavelen %g < lowest sampled wavelen %g\n",
01531                upper_wavelen, spectrum->x[0]);
01532         return NULL;
01533     } else if (last <= -2)  {
01534         last = spectrum->nx-1;
01535         bu_log("bn_tabdata_mk_linear_filter() warning, upper_wavelen %g > highest sampled wavelen %g\n",
01536                upper_wavelen, spectrum->x[spectrum->nx]);
01537     }
01538 
01539     /* 'filt' is filled with zeros by default */
01540     BN_GET_TABDATA(filt, spectrum);
01541 
01542     /* Special case:  first and last are in same sample cell */
01543     if (first == last)  {
01544         filt_range = upper_wavelen - lower_wavelen;
01545         cell_range = spectrum->x[first+1] - spectrum->x[first];
01546         frac = filt_range / cell_range;
01547 
01548         /* Could use a BU_ASSERT_RANGE_FLOAT */
01549         BU_ASSERT((frac >= 0.0) && (frac <= 1.0));
01550 
01551         filt->y[first] = frac;
01552         return filt;
01553     }
01554 
01555     /* Calculate fraction 0..1.0 for first and last samples */
01556     filt_range = spectrum->x[first+1] - lower_wavelen;
01557     cell_range = spectrum->x[first+1] - spectrum->x[first];
01558     frac = filt_range / cell_range;
01559     if (frac > 1)  frac = 1;
01560     BU_ASSERT((frac >= 0.0) && (frac <= 1.0));
01561     filt->y[first] = frac;
01562 
01563     filt_range = upper_wavelen - spectrum->x[last];
01564     cell_range = spectrum->x[last+1] - spectrum->x[last];
01565     frac = filt_range / cell_range;
01566     if (frac > 1)  frac = 1;
01567     BU_ASSERT((frac >= 0.0) && (frac <= 1.0));
01568     filt->y[last] = frac;
01569 
01570     /* Fill in range between with 1.0 values */
01571     for (i = first+1; i < (size_t)last; i++)
01572         filt->y[i] = 1.0;
01573 
01574     return filt;
01575 }
01576 /** @} */
01577 /*
01578  * Local Variables:
01579  * mode: C
01580  * tab-width: 8
01581  * indent-tabs-mode: t
01582  * c-file-style: "stroustrup"
01583  * End:
01584  * ex: shiftwidth=4 tabstop=8
01585  */
Generated on Tue Dec 11 13:14:28 2012 for LIBBN by  doxygen 1.6.3