BRL-CAD
spectrum.c
Go to the documentation of this file.
1 /* S P E C T R U M . 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 /** @addtogroup librt */
21 /** @{ */
22 /** @file librt/spectrum.c
23  *
24  * An application of the 'tabdata' package to spectral data.
25  *
26  * (Note that there is also a rttherm/spectrum.c)
27  *
28  * Inspired by -
29  * Roy Hall and his book "Illumination and Color in Computer
30  * Generated Imagery", Springer Verlag, New York, 1989.
31  * ISBN 0-387-96774-5
32  *
33  * With thanks to Russ Moulton Jr, EOSoft Inc. for his "rad.c" module.
34  */
35 /** @} */
36 
37 #include "common.h"
38 
39 #include <math.h>
40 #include "bio.h"
41 
42 #include "bu/debug.h"
43 #include "vmath.h"
44 #include "raytrace.h"
45 #include "spectrum.h"
46 
47 
48 /**
49  * This is the data for the CIE_XYZ curves take from Judd and Wyszecki
50  * (1975), table 2.6, these are for the 1931 standard observer with a
51  * 2-degree visual field. From Roy Hall, pg 228.
52  */
53 static const double rt_CIE_XYZ[81][4] = {
54  {380, 0.0014, 0.0000, 0.0065}, {385, 0.0022, 0.0001, 0.0105},
55  {390, 0.0042, 0.0001, 0.0201}, {395, 0.0076, 0.0002, 0.0362},
56  {400, 0.0143, 0.0004, 0.0679}, {405, 0.0232, 0.0006, 0.1102},
57  {410, 0.0435, 0.0012, 0.2074}, {415, 0.0776, 0.0022, 0.3713},
58  {420, 0.1344, 0.0040, 0.6456}, {425, 0.2148, 0.0073, 1.0391},
59  {430, 0.2839, 0.0116, 1.3856}, {435, 0.3285, 0.0168, 1.6230},
60  {440, 0.3483, 0.0230, 1.7471}, {445, 0.3481, 0.0298, 1.7826},
61  {450, 0.3362, 0.0380, 1.7721}, {455, 0.3187, 0.0480, 1.7441},
62  {460, 0.2908, 0.0600, 1.6692}, {465, 0.2511, 0.0739, 1.5281},
63  {470, 0.1954, 0.0910, 1.2876}, {475, 0.1421, 0.1126, 1.0419},
64  {480, 0.0956, 0.1390, 0.8130}, {485, 0.0580, 0.1693, 0.6162},
65  {490, 0.0320, 0.2080, 0.4652}, {495, 0.0147, 0.2586, 0.3533},
66  {500, 0.0049, 0.3230, 0.2720}, {505, 0.0024, 0.4073, 0.2123},
67  {510, 0.0093, 0.5030, 0.1582}, {515, 0.0291, 0.6082, 0.1117},
68  {520, 0.0633, 0.7100, 0.0782}, {525, 0.1096, 0.7932, 0.0573},
69  {530, 0.1655, 0.8620, 0.0422}, {535, 0.2257, 0.9149, 0.0298},
70  {540, 0.2904, 0.9540, 0.0203}, {545, 0.3597, 0.9803, 0.0134},
71  {550, 0.4334, 0.9950, 0.0087}, {555, 0.5121, 1.0000, 0.0057},
72  {560, 0.5945, 0.9950, 0.0039}, {565, 0.6784, 0.9786, 0.0027},
73  {570, 0.7621, 0.9520, 0.0021}, {575, 0.8425, 0.9154, 0.0018},
74  {580, 0.9163, 0.8700, 0.0017}, {585, 0.9786, 0.8163, 0.0014},
75  {590, 1.0263, 0.7570, 0.0011}, {595, 1.0567, 0.6949, 0.0010},
76  {600, 1.0622, 0.6310, 0.0008}, {605, 1.0456, 0.5668, 0.0006},
77  {610, 1.0026, 0.5030, 0.0003}, {615, 0.9384, 0.4412, 0.0002},
78  {620, 0.8544, 0.3810, 0.0002}, {625, 0.7514, 0.3210, 0.0001},
79  {630, 0.6424, 0.2650, 0.0000}, {635, 0.5419, 0.2170, 0.0000},
80  {640, 0.4479, 0.1750, 0.0000}, {645, 0.3608, 0.1382, 0.0000},
81  {650, 0.2835, 0.1070, 0.0000}, {655, 0.2187, 0.0816, 0.0000},
82  {660, 0.1649, 0.0610, 0.0000}, {665, 0.1212, 0.0446, 0.0000},
83  {670, 0.0874, 0.0320, 0.0000}, {675, 0.0636, 0.0232, 0.0000},
84  {680, 0.0468, 0.0170, 0.0000}, {685, 0.0329, 0.0119, 0.0000},
85  {690, 0.0227, 0.0082, 0.0000}, {695, 0.0158, 0.0057, 0.0000},
86  {700, 0.0114, 0.0041, 0.0000}, {705, 0.0081, 0.0029, 0.0000},
87  {710, 0.0058, 0.0021, 0.0000}, {715, 0.0041, 0.0015, 0.0000},
88  {720, 0.0029, 0.0010, 0.0000}, {725, 0.0020, 0.0007, 0.0000},
89  {730, 0.0014, 0.0005, 0.0000}, {735, 0.0010, 0.0004, 0.0000},
90  {740, 0.0007, 0.0002, 0.0000}, {745, 0.0005, 0.0002, 0.0000},
91  {750, 0.0003, 0.0001, 0.0000}, {755, 0.0002, 0.0001, 0.0000},
92  {760, 0.0002, 0.0001, 0.0000}, {765, 0.0001, 0.0000, 0.0000},
93  {770, 0.0001, 0.0000, 0.0000}, {775, 0.0001, 0.0000, 0.0000},
94  {780, 0.0000, 0.0000, 0.0000}
95 };
96 
97 
98 /**
99  * Given as input a spectral sampling distribution, generate the 3
100  * curves to match the human eye's response in CIE color parameters X,
101  * Y, and Z. XYZ space can be readily converted to RGB with a 3x3
102  * matrix.
103  *
104  * The tabulated data is linearly interpolated.
105  *
106  * Pointers to the three spectral weighting functions are "returned",
107  * storage for the X, Y, and Z curves is allocated by this routine and
108  * must be freed by the caller.
109  */
110 void
111 rt_spect_make_CIE_XYZ(struct bn_tabdata **x, struct bn_tabdata **y, struct bn_tabdata **z, const struct bn_table *tabp)
112 {
113  struct bn_tabdata *a, *b, *c;
114  fastf_t xyz_scale;
115  int i;
116  size_t j;
117 
118  BN_CK_TABLE(tabp);
119 
120  i = bn_table_interval_num_samples(tabp, 430., 650.);
121  if (i <= 4) bu_log("rt_spect_make_CIE_XYZ: insufficient samples (%d) in visible band\n", i);
122 
123  BN_GET_TABDATA(a, tabp);
124  BN_GET_TABDATA(b, tabp);
125  BN_GET_TABDATA(c, tabp);
126  *x = a;
127  *y = b;
128  *z = c;
129 
130  /* No CIE data below 380 nm */
131  for (j=0; tabp->x[j] < 380 && j < tabp->nx; j++) {
132  a->y[j] = b->y[j] = c->y[j] = 0;
133  }
134 
135  /* Traverse the CIE table. Produce as many output values as
136  * possible before advancing to next CIE table entry.
137  */
138  for (i = 0; i < 81-1; i++) {
139  fastf_t fract; /* fraction from [i] to [i+1] */
140 
141  again:
142  if (j >= tabp->nx) break;
143  if (tabp->x[j] < rt_CIE_XYZ[i][0]) bu_bomb("rt_spect_make_CIE_XYZ assertion1 failed\n");
144  if (tabp->x[j] >= rt_CIE_XYZ[i+1][0]) continue;
145  /* The CIE table has 5nm spacing */
146  fract = (tabp->x[j] - rt_CIE_XYZ[i][0]) / 5;
147  if (fract < 0 || fract > 1) bu_bomb("rt_spect_make_CIE_XYZ assertion2 failed\n");
148  a->y[j] = (1-fract) * rt_CIE_XYZ[i][1] + fract * rt_CIE_XYZ[i+1][1];
149  b->y[j] = (1-fract) * rt_CIE_XYZ[i][2] + fract * rt_CIE_XYZ[i+1][2];
150  c->y[j] = (1-fract) * rt_CIE_XYZ[i][3] + fract * rt_CIE_XYZ[i+1][3];
151  j++;
152  goto again;
153  }
154 
155  /* No CIE data above 780 nm */
156  for (; j < tabp->nx; j++) {
157  a->y[j] = b->y[j] = c->y[j] = 0;
158  }
159 
160  /* Normalize the curves so that area under Y curve is 1.0 */
161  xyz_scale = bn_tabdata_area2(b);
162  if (fabs(xyz_scale) < VDIVIDE_TOL) {
163  bu_log("rt_spect_make_CIE_XYZ(): Area = 0 (no luminance) in this part of the spectrum, skipping normalization step\n");
164  return;
165  }
166  xyz_scale = 1 / xyz_scale;
167  bn_tabdata_scale(a, a, xyz_scale);
168  bn_tabdata_scale(b, b, xyz_scale);
169  bn_tabdata_scale(c, c, xyz_scale);
170 }
171 
172 
173 /**
174  * Given reflectance data (in range 0..1) in terms of RGB color,
175  * convert that to a spectral reflectance curve.
176  *
177  * The assumption here is that the spectrum is made up of exactly
178  * three non-overlapping bands, and the reflectance is constant over
179  * each:
180  *
181  * red ===> 572nm to 1, 000, 000nm (includes the full IR band)
182  * green => 492nm to 572nm (just green)
183  * blue ==> 1nm to 492nm (includes Ultraviolet)
184  *
185  * As the caller may be doing a lot of this, the caller is expected to
186  * provide a pointer to a valid bn_tabdata structure which is to be
187  * filled in. Allowing caller to re-cycle them rather than doing
188  * constant malloc/free cycle.
189  */
190 void
191 rt_spect_reflectance_rgb(struct bn_tabdata *curve, const float *rgb)
192 {
193  register size_t i;
194  register const struct bn_table *tabp;
195 
196  BN_CK_TABDATA(curve);
197  tabp = curve->table;
198  BN_CK_TABLE(tabp);
199 
200  /* Fill in blue values, everything up to but not including 492nm */
201  for (i=0; i < tabp->nx; i++) {
202  if (tabp->x[i] >= 492) break;
203  curve->y[i] = rgb[2];
204  }
205 
206  /* Fill in green values, everything up to but not including 572nm */
207  for (; i < tabp->nx; i++) {
208  if (tabp->x[i] >= 572) break;
209  curve->y[i] = rgb[1];
210  }
211 
212  /* Fill in red values, everything from here up to end of table */
213  for (; i < tabp->nx; i++) {
214  curve->y[i] = rgb[0];
215  }
216 }
217 
218 
219 #define PLANCK_C1 37415 /* watts um^4 cm^-2 */
220 #define PLANCK_C2 14387.86 /* um K */
221 /* Russ gives these values at 37, 415 and 14, 388 */
222 /* Handbook of Physics and Chem gives these values as 37, 403 and 14, 384 */
223 /* Aircraft Combat Surv gives these values as 37, 483.2 and 14, 387.86 */
224 
225 /* Requires wavelength _w in um, not nm, returns units: W / cm**2 / um */
226 #define PLANCK(_w, _tempK) \
227  (PLANCK_C1/(_w*_w*_w*_w*_w*(exp(PLANCK_C2/(_w*_tempK))-1)))
228 
229 /**
230  * Integrate Planck's Radiation Formula for a black body radiator
231  * across the given spectrum. Returns radiant emittance in W/cm**2
232  * for each wavelength interval.
233  *
234  * Based upon code kindly provided by Russ Moulton, Jr., EOSoft Inc.
235  * Compute at 'n-1' wavelengths evenly spaced between ax and bx.
236  */
237 void
238 rt_spect_black_body(struct bn_tabdata *data, double temp, unsigned int n)
239 
240 /* Degrees Kelvin */
241 /* # wavelengths to eval at */
242 {
243  const struct bn_table *tabp;
244  size_t j;
245 
246  BN_CK_TABDATA(data);
247  tabp = data->table;
248  BN_CK_TABLE(tabp);
249 
251  bu_log("rt_spect_black_body(%p, %g degK) %g um to %g um\n",
252  (void *)data, temp,
253  tabp->x[0] * 0.001, /* nm to um */
254  tabp->x[tabp->nx] * 0.001 /* nm to um */
255  );
256  }
257 
258  if (n < 3) n = 3;
259 
260  for (j = 0; j < tabp->nx; j++) {
261  double ax; /* starting wavelength, um */
262  double bx; /* ending wavelength, um */
263  double dx; /* wavelength interval, um */
264  double w_sum; /* sum over wavelengths */
265  double wavlen; /* current wavelength */
266  unsigned long i;
267 
268  ax = tabp->x[j] * 0.001; /* nm to um */
269  bx = tabp->x[j+1] * 0.001; /* nm to um */
270  dx = (bx - ax) / (double)n;
271 
272  w_sum = 0;
273  wavlen = ax;
274  for (i=0; i<n; i++) {
275  w_sum += PLANCK(wavlen, temp);
276  wavlen += dx;
277  }
278  w_sum *= dx;
279 
280  data->y[j] = w_sum;
281  }
282 }
283 
284 
285 /**
286  * Returns radiant emittance for each spectral interval in the given
287  * spectrum in units of watts/cm**2. Integrate each wavelength
288  * interval of spectral radiant emittance, by fitting with a rectangle
289  * (approximating curve with a horizontal line). For narrow spacing
290  * in wavelength this is OK, but with large spacing this tends to
291  * over-predict the power by 20%, due to the sharp (exponential) slope
292  * of the curve. With coarse spacing, or when unsure, use
293  * rt_spect_black_body().
294  */
295 void
297 
298 /* Degrees Kelvin */
299 {
300  const struct bn_table *tabp;
301  size_t j;
302 
303  BN_CK_TABDATA(data);
304  tabp = data->table;
305  BN_CK_TABLE(tabp);
306 
308  bu_log("rt_spect_black_body_fast(%p, %g degK)\n",
309  (void *)data, temp);
310  }
311 
312  for (j = 0; j < tabp->nx; j++) {
313  data->y[j] = PLANCK((tabp->x[j]*0.001), temp) *
314  (tabp->x[j+1] - tabp->x[j]) * 0.001;
315  }
316 }
317 
318 
319 /**
320  * Returns point-sampled values of spectral radiant emittance, in
321  * units of watts/cm**2/um, straight from Planck's black-body
322  * radiation formula.
323  */
324 void
326 
327 /* Degrees Kelvin */
328 {
329  const struct bn_table *tabp;
330  size_t j;
331 
332  BN_CK_TABDATA(data);
333  tabp = data->table;
334  BN_CK_TABLE(tabp);
335 
337  bu_log("rt_spect_black_body_points(%p, %g degK)\n",
338  (void *)data, temp);
339  }
340 
341  for (j = 0; j < tabp->nx; j++) {
342  data->y[j] = PLANCK((tabp->x[j]*0.001), temp);
343  }
344 }
345 
346 
347 /*
348  * Local Variables:
349  * mode: C
350  * tab-width: 8
351  * indent-tabs-mode: t
352  * c-file-style: "stroustrup"
353  * End:
354  * ex: shiftwidth=4 tabstop=8
355  */
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
#define BN_GET_TABDATA(_data, _table)
Definition: tabdata.h:114
void rt_spect_black_body(struct bn_tabdata *data, double temp, unsigned int n)
Definition: spectrum.c:238
Header file for the BRL-CAD common definitions.
double bn_tabdata_area2(const struct bn_tabdata *in)
Definition: tabdata.c:433
#define BN_CK_TABDATA(_p)
Definition: tabdata.h:106
COMPLEX data[64]
Definition: fftest.c:34
size_t nx
Definition: tabdata.h:75
#define BN_CK_TABLE(_p)
Definition: tabdata.h:79
void bn_tabdata_scale(struct bn_tabdata *out, const struct bn_tabdata *in1, double scale)
fastf_t y[1]
array of ny samples, dynamically sized
Definition: tabdata.h:104
void rt_spect_black_body_points(struct bn_tabdata *data, double temp)
Definition: spectrum.c:325
void rt_spect_reflectance_rgb(struct bn_tabdata *curve, const float *rgb)
Definition: spectrum.c:191
#define PLANCK(_w, _tempK)
Definition: spectrum.c:226
int bu_debug
Definition: globals.c:87
fastf_t x[1]
array of nx+1 wavelengths, dynamically sized
Definition: tabdata.h:76
void rt_spect_black_body_fast(struct bn_tabdata *data, double temp)
Definition: spectrum.c:296
#define BU_DEBUG_TABDATA
Definition: debug.h:73
void bu_bomb(const char *str) _BU_ATTR_NORETURN
Definition: bomb.c:91
double fastf_t
Definition: defines.h:300
const struct bn_table * table
Up pointer to definition of X axis.
Definition: tabdata.h:103
void rt_spect_make_CIE_XYZ(struct bn_tabdata **x, struct bn_tabdata **y, struct bn_tabdata **z, const struct bn_table *tabp)
Definition: spectrum.c:111