Spectrum.cc

00001 /*
00002 File: Spectrum.cc
00003 
00004 Function: Provides classes for manipulating spectra
00005 
00006 Author: Andrew Willmott
00007 
00008 Notes: 
00009 
00010 */
00011 
00012 #include <stdlib.h>
00013 #include "Spectrum.h"
00014 #include "ColourSystem.h"
00015 #include "VecUtil.h"
00016 
00017 const Int Spectrum::kComponents = 100; // number of basis functions
00018 const Float Spectrum::kStart = 350.0; // start of spectrum
00019 const Float Spectrum::kEnd = 800.0; // end of spectrum
00020 
00021 // CIE colour matching functions, x_bar, y_bar and z_bar
00022 // 360-830nm in steps of 5nm
00023 // good source for this stuff is http://cvision.ucsd.edu/basicindex.htm
00024 
00025  const Float kCIE_Start = 360.0;
00026  const Float kCIE_End = 830.0;
00027  const Float kCIE_Delta = 5.0;
00028 
00029  static Float tCIE_xbar[95] =
00030 {
00031 #include "cie_xbar.spectrum" 
00032 };
00033 RegularSpectralCurve kCIE_xbar(tCIE_xbar, 360, 830, 95);
00034 
00035  static Float tCIE_ybar[95] =
00036 {
00037 #include "cie_ybar.spectrum" 
00038 };
00039 RegularSpectralCurve kCIE_ybar(tCIE_ybar, 360, 830, 95);
00040 
00041  static Float tCIE_zbar[95] =
00042 {
00043 #include "cie_zbar.spectrum" 
00044 };
00045 RegularSpectralCurve kCIE_zbar(tCIE_zbar, 360, 830, 95);
00046 
00047  Colour Spectrum::ToCIE_XYZ()
00048 {
00049 Int i;
00050 Colour c = vl_0;
00051 
00052 for (i = 0; i < kComponents; i++)
00053 {
00054 ClrReal x = (i + 0.5) / kComponents;
00055 ClrReal lambda = mix(kStart, kEnd, x);
00056 
00057 c[0] += kCIE_xbar.Sample(lambda) * SELF[i];
00058 c[1] += kCIE_ybar.Sample(lambda) * SELF[i];
00059 c[2] += kCIE_zbar.Sample(lambda) * SELF[i];
00060 }
00061 
00062 return(c);
00063 }
00064 
00065 // --- Spectral curve methods -------------------------------------------------
00066 
00067  SpectralCurve::~SpectralCurve()
00068 {
00069 }
00070 
00071  SpectralCurve::operator Spectrum()
00072 {
00073 Spectrum result;
00074 Float lambda, x;
00075 Int i;
00076 
00077 #ifdef USE_INT
00078 Float delta, integral, lastIntegral;
00079 cout << "INTEGRATING" << endl;
00080 delta = (Spectrum::kEnd - Spectrum::kStart) / Spectrum::kComponents;
00081 lambda = Spectrum::kStart;
00082 lastIntegral = Integral(lambda);
00083 
00084 for (i = 0; i < Spectrum::kComponents; i++)
00085 {
00086 lambda += delta;
00087 integral = Integral(lambda);
00088 // normalised box functions, so we divide by the support size
00089 result[i] = (integral - lastIntegral) / delta;
00090 lastIntegral = integral;
00091 }
00092 #else
00093 
00094 for (i = 0; i < Spectrum::kComponents; i++)
00095 {
00096 x = (i + 0.5) / Spectrum::kComponents;
00097 lambda = Spectrum::kStart + x * (Spectrum::kEnd - Spectrum::kStart);
00098 
00099 result[i] = Sample(lambda);
00100 }
00101 
00102 #endif
00103 
00104 return(result);
00105 }
00106 
00107  Float SpectralCurve::Integral(Float lambda) const
00108 {
00109 const Float kIntegrateDelta = 1.0;
00110 Float s, sum = 0.0;
00111 Int n = 0;
00112 
00113 for (s = Spectrum::kStart; s < lambda; s += kIntegrateDelta)
00114 {
00115 sum += Sample(s);
00116 n++;
00117 }
00118 
00119 if (n > 0)
00120 sum *= (lambda - Spectrum::kStart) / n;
00121 
00122 return(sum);
00123 }
00124 
00125  Float SpectralCurve::MaxCmpt() const
00126 {
00127 const Float kMaxDelta = 1.0;
00128 Float s, max = Sample(Spectrum::kStart);
00129 
00130 for (s = Spectrum::kStart + kMaxDelta; s <= Spectrum::kEnd; s += kMaxDelta)
00131 {
00132 max = Max(max, Sample(s));
00133 }
00134 
00135 return(max);
00136 }
00137 
00138  Void SpectralCurve::MakeImage(Image &image, ClrReal maxAmplitude)
00139 // Plot spectrum into image
00140 {
00141 Int i, j;
00142 Float *samples = new Float[image.Width()], ampMax = 0.0;
00143 Colour *waveColour = new Colour[image.Width()];
00144 
00145 image.Clear(cBlack);
00146 
00147 for (i = 0; i < image.Width(); i++)
00148 {
00149 Float x = (i + 0.5) / image.Width();
00150 Float lambda = mix(Spectrum::kStart, Spectrum::kEnd, x);
00151 
00152 samples[i] = Sample(lambda);
00153 waveColour[i] = WavelengthToXYZ(lambda);
00154 waveColour[i] = csDisplay.XYZToGamutClip(waveColour[i]);
00155 ampMax = Max(samples[i], ampMax);
00156 }
00157 
00158 if (maxAmplitude < 0.0)
00159 maxAmplitude = ampMax;
00160 
00161 for (j = 0; j < image.Height(); j++)
00162 { 
00163 Float y = maxAmplitude * (j + 0.5) / image.Height();
00164 
00165 for (i = 0; i < image.Width(); i++)
00166 if (samples[i] >= y)
00167 image.SetPixel(i, j, waveColour[i]);
00168 }
00169 
00170 delete[] samples;
00171 delete[] waveColour;
00172 }
00173 
00174  Colour SpectralCurve::ToCIE_XYZ()
00175 {
00176 Float lambda, sample;
00177 Colour c = vl_0;
00178 
00179 for (lambda = kCIE_Start; lambda <= kCIE_End; lambda += kCIE_Delta)
00180 {
00181 sample = Sample(lambda);
00182 c[0] += kCIE_xbar.Sample(lambda) * sample;
00183 c[1] += kCIE_ybar.Sample(lambda) * sample;
00184 c[2] += kCIE_zbar.Sample(lambda) * sample;
00185 }
00186 
00187 return(c);
00188 }
00189 
00190  static Float *MakeAccumData(Float *data, Int divisions)
00191 {
00192 Int i;
00193 Float *result, sum = 0;
00194 
00195 result = new Float[divisions];
00196 
00197 for (i = 0; i < divisions; i++)
00198 {
00199 sum += data[i];
00200 result[i] = sum;
00201 }
00202 
00203 return(result);
00204 }
00205 
00206  RegularSpectralCurve::RegularSpectralCurve(Float *sdata, Int start, Int end, Int divisions) : 
00207 data(sdata),
00208 waveDivs(divisions),
00209 waveStart(start),
00210 waveEnd(end)
00211 {
00212 waveDelta = (end - start) / (divisions - 1);
00213 accumData = MakeAccumData(data, divisions);
00214 }
00215 
00216  RegularSpectralCurve::~RegularSpectralCurve()
00217 {
00218 delete[] accumData;
00219 }
00220 
00221  Float RegularSpectralCurve::Sample(Float lambda) const
00222 // interpolate between two closest samples
00223 {
00224 Float result, x, dx;
00225 Int bin1, bin2;
00226 
00227 if (lambda < waveStart || lambda > waveEnd)
00228 return(0.0);
00229 
00230 x = (lambda - waveStart) / waveDelta;
00231 
00232 bin1 = (Int) x;
00233 dx = x - bin1;
00234 if (dx > 1e-8)
00235 bin2 = bin1 + 1;
00236 else
00237 bin2 = bin1;
00238 
00239 CheckRange(bin1, 0, waveDivs, "Bad bin number");
00240 CheckRange(bin2, 0, waveDivs, "Bad bin number");
00241 
00242 result = mix(data[bin1], data[bin2], dx);
00243 
00244 return(result);
00245 }
00246 
00247 #ifdef UNFINISHED
00248 Float RegularSpectralCurve::Integral(Float lambda) const
00249 // interpolate between two closest samples
00250 {
00251 Float result, x, dx;
00252 Int bin1, bin2;
00253 
00254 if (lambda < waveStart)
00255 return(0.0);
00256 else if (lambda > waveEnd)
00257 return(waveDelta * accumData[waveDivs - 1]);
00258 
00259 x = (lambda - waveStart) / waveDelta;
00260 
00261 bin1 = (Int) x;
00262 dx = x - bin1;
00263 if (dx > 1e-8)
00264 bin2 = bin1 + 1;
00265 else
00266 bin2 = bin1;
00267 
00268 CheckRange(bin1, 0, waveDivs, "Bad bin number");
00269 CheckRange(bin2, 0, waveDivs, "Bad bin number");
00270 
00271 result = waveDelta * mix(accumData[bin1], accumData[bin2], dx);
00272 
00273 return(result);
00274 }
00275 #endif
00276 
00277  IrregularSpectralCurve::IrregularSpectralCurve(Float *amps, Float *wavs, Int divs) :
00278 amplitudes(amps),
00279 wavelengths(wavs),
00280 divisions(divs)
00281 {
00282 accumAmplitudes = MakeAccumData(amps, divs);
00283 }
00284 
00285  IrregularSpectralCurve::IrregularSpectralCurve() :
00286 amplitudes(0),
00287 wavelengths(0),
00288 divisions(0),
00289 accumAmplitudes(0)
00290 {
00291 }
00292 
00293  IrregularSpectralCurve::~IrregularSpectralCurve()
00294 {
00295 delete[] accumAmplitudes;
00296 }
00297 
00298  IrregularSpectralCurve &IrregularSpectralCurve::ReadSDFile(FileName filename)
00299 {
00300 String line;
00301 StrConstArray tokens;
00302 ifstream s;
00303 Int i = 0;
00304 
00305 cout << "parsing " << filename.GetPath() << endl;
00306 
00307 s.open(filename.GetPath());
00308 
00309 if (!s)
00310 {
00311 cerr << "(ReadSDFile) Cannot access " << filename.GetPath() << endl;
00312 return(SELF);
00313 }
00314 
00315 while (s)
00316 {
00317 if (line.ReadLine(s))
00318 {
00319 Split(line, tokens);
00320 if (tokens.NumItems() == 0 || tokens[0][0] == '#')
00321 ;
00322 else if (tokens.NumItems() == 1)
00323 {
00324 divisions = atoi(tokens[0]);
00325 cout << "setting size to " << divisions << endl;
00326 delete[] wavelengths;
00327 delete[] amplitudes;
00328 wavelengths = new Float[divisions];
00329 amplitudes = new Float[divisions];
00330 i = 0;
00331 }
00332 else if (tokens.NumItems() >= 2)
00333 {
00334 if (!wavelengths)
00335 _Error("(ReadSDFile) didn't find number of divisions");
00336 wavelengths[i] = atof(tokens[0]);
00337 amplitudes[i] = atof(tokens[1]);
00338 i++;
00339 }
00340 }
00341 } 
00342 
00343 s.close();
00344 accumAmplitudes = MakeAccumData(amplitudes, divisions);
00345 
00346 return(SELF);
00347 }
00348 
00349  Float IrregularSpectralCurve::Sample(Float lambda) const
00350 // interpolate between two closest samples
00351 // need binary search or something to speed this up. (cache?)
00352 // currently we extend the endpoint samples out to infinity either
00353 // side. Perhaps should ramp down over X nm instead?
00354 {
00355 Float x, result;
00356 Int i;
00357 
00358 if (divisions == 0)
00359 return(0.0);
00360 else if (divisions == 1)
00361 return(amplitudes[0]);
00362 
00363 if (lambda < wavelengths[0])
00364 return(amplitudes[0]);
00365 
00366 i = 1;
00367 while (i < divisions)
00368 if (lambda < wavelengths[i])
00369 {
00370 x = (lambda - wavelengths[i - 1])
00371 / (wavelengths[i] - wavelengths[i - 1]);
00372 
00373 result = mix(amplitudes[i - 1], amplitudes[i], x);
00374 
00375 return(result);
00376 }
00377 else
00378 i++;
00379 
00380 return(amplitudes[divisions - 1]);
00381 }
00382 
00383 #ifdef UNFINISHED
00384 Float IrregularSpectralCurve::Integral(Float lambda) const
00385 // interpolate between two closest samples
00386 // need binary search or something to speed this up. (cache?)
00387 // currently we extend the endpoint samples out to infinity either
00388 // side. Perhaps should ramp down over X nm instead?
00389 {
00390 Float x, result;
00391 Int i;
00392 
00393 if (divisions == 0)
00394 return(0.0);
00395 else if (divisions == 1)
00396 return(amplitudes[0]);
00397 
00398 if (lambda < wavelengths[0])
00399 return(amplitudes[0]);
00400 
00401 i = 1;
00402 while (i < divisions)
00403 if (lambda < wavelengths[i])
00404 {
00405 x = (lambda - wavelengths[i - 1])
00406 / (wavelengths[i] - wavelengths[i - 1]);
00407 
00408 result = mix(amplitudes[i - 1], amplitudes[i], x);
00409 
00410 return(result);
00411 }
00412 else
00413 i++;
00414 
00415 return(amplitudes[divisions - 1]);
00416 }
00417 #endif
00418 
00419 // --- Utility routines -------------------------------------------------------
00420 
00421  Colour WavelengthToXYZ(ClrReal lambda)
00422 {
00423 Colour c;
00424 
00425 // calculate XYZ
00426 c[0] = kCIE_xbar.Sample(lambda);
00427 c[1] = kCIE_ybar.Sample(lambda);
00428 c[2] = kCIE_zbar.Sample(lambda);
00429 
00430 return(c);
00431 }
00432 
00433 // chromaticity stuff originally from utah's sunsky code
00434 
00435 // 300-830 10nm
00436  static Float tS0Amplitudes[54] =
00437 {
00438 0.04,6.0,29.6,55.3,57.3,
00439 61.8,61.5,68.8,63.4,65.8,
00440 94.8,104.8,105.9,96.8,113.9,
00441 125.6,125.5,121.3,121.3,113.5,
00442 113.1,110.8,106.5,108.8,105.3,
00443 104.4,100.0,96.0,95.1,89.1,
00444 90.5,90.3,88.4,84.0,85.1,
00445 81.9,82.6,84.9,81.3,71.9,
00446 74.3,76.4,63.3,71.7,77.0,
00447 65.2,47.7,68.6,65.0,66.0,
00448 61.0,53.3,58.9,61.9
00449 };
00450 
00451  static Float tS1Amplitudes[54] =
00452 {
00453 0.02,4.5,22.4,42.0,40.6,
00454 41.6,38.0,42.4,38.5,35.0,
00455 43.4,46.3,43.9,37.1,36.7,
00456 35.9,32.6,27.9,24.3,20.1,
00457 16.2,13.2,8.6,6.1,4.2,
00458 1.9,0.0,-1.6,-3.5,-3.5,
00459 -5.8,-7.2,-8.6,-9.5,-10.9,
00460 -10.7,-12.0,-14.0,-13.6,-12.0,
00461 -13.3,-12.9,-10.6,-11.6,-12.2,
00462 -10.2,-7.8,-11.2,-10.4,-10.6,
00463 -9.7,-8.3,-9.3,-9.8
00464 };
00465 
00466  static Float tS2Amplitudes[54] =
00467 {
00468 0.0,2.0,4.0,8.5,7.8,
00469 6.7,5.3,6.1,3.0,1.2,
00470 -1.1,-0.5,-0.7,-1.2,-2.6,
00471 -2.9,-2.8,-2.6,-2.6,-1.8,
00472 -1.5,-1.3,-1.2,-1.0,-0.5,
00473 -0.3,0.0,0.2,0.5,2.1,
00474 3.2,4.1,4.7,5.1,6.7,
00475 7.3,8.6,9.8,10.2,8.3,
00476 9.6,8.5,7.0,7.6,8.0,
00477 6.7,5.2,7.4,6.8,7.0,
00478 6.4,5.5,6.1,6.5
00479 };
00480 
00481 RegularSpectralCurve kS0Spectrum(tS0Amplitudes, 300, 830, 54);
00482 RegularSpectralCurve kS1Spectrum(tS1Amplitudes, 300, 830, 54);
00483 RegularSpectralCurve kS2Spectrum(tS2Amplitudes, 300, 830, 54);
00484 
00485  ChromaticitySpectrum::ChromaticitySpectrum(ClrReal x, ClrReal y)
00486 {
00487 M1 = (-1.3515 - 1.7703 * x + 5.9114 * y) / (0.0241 + 0.2562 * x - 0.7341 * y);
00488 M2 = ( 0.03 -31.4424 * x + 30.0717 * y) / (0.0241 + 0.2562 * x - 0.7341 * y);
00489 }
00490 
00491  Float ChromaticitySpectrum::Sample(Float lambda) const
00492 {
00493 return(kS0Spectrum.Sample(lambda) + M1 * kS1Spectrum.Sample(lambda)
00494 + M2 * kS2Spectrum.Sample(lambda));
00495 }
00496 
00497  static Coord WavelengthToChroma(Float lambda)
00498 {
00499 Colour c;
00500 
00501 // calculate XYZ
00502 c[0] = kCIE_xbar.Sample(lambda);
00503 c[1] = kCIE_ybar.Sample(lambda);
00504 c[2] = kCIE_zbar.Sample(lambda);
00505 
00506 c /= (c[0] + c[1] + c[2]);
00507 
00508 return(Coord(c[0], c[1]));
00509 }
00510 
00511  Void DrawChromaticityImage(Image &image)
00512 // Plot CIE colour diagram
00513 {
00514 Int i, j;
00515 Colour chromaClr;
00516 Coord a1, a2, b1, b2;
00517 Float aLambda, bLambda;
00518 Int startX, endX;
00519 Bool lastClip, clip;
00520 
00521 // the only tricky part to drawing the CIE
00522 // diagram is culling everything outside the (x,y)
00523 // curve defined by the pure wavelengths, and the
00524 // "line of purples" between the highest and lowest
00525 // wavelength in the spectrum.
00526 
00527 // to take care of this we use a1/a2/b1/b2 to step
00528 // through the spectral curve from either end, starting
00529 // from the beginning, and pick the correct start and
00530 // ending x values for each scanline.
00531 
00532 image.Clear(cGrey50);
00533 
00534 a1 = WavelengthToChroma(kCIE_Start);
00535 aLambda = kCIE_Start + kCIE_Delta;
00536 a2 = WavelengthToChroma(aLambda);
00537 b1 = a1;
00538 bLambda = kCIE_End;
00539 b2 = WavelengthToChroma(bLambda);
00540 
00541 for (j = 0; j < image.Height(); j++)
00542 { 
00543 Float y = (j + 0.5) / image.Height();
00544 
00545 if (y < a1[1])
00546 continue;
00547 // use 'while' here because some consecutive wavelength
00548 // samples will have the same (x,y) values
00549 while (y > a2[1])
00550 {
00551 a1 = a2;
00552 aLambda += kCIE_Delta;
00553 a2 = WavelengthToChroma(aLambda);
00554 }
00555 while (y > b2[1])
00556 {
00557 b1 = b2;
00558 bLambda -= kCIE_Delta;
00559 b2 = WavelengthToChroma(bLambda);
00560 }
00561 
00562 if (bLambda < aLambda)
00563 // all done
00564 break;
00565 
00566 startX = Int(image.Width() * mix(a1[0], a2[0], (y - a1[1]) / (a2[1] - a1[1])));
00567 endX = Int(image.Width() * mix(b1[0], b2[0], (y - b1[1]) / (b2[1] - b1[1])));
00568 if (endX >= image.Width())
00569 endX = image.Width() - 1;
00570 
00571 lastClip = true;
00572 
00573 for (i = startX; i < endX; i++)
00574 {
00575 Float x = (i + 0.5) / image.Width();
00576 
00577 chromaClr = Colour(x, y, (1.0 - x - y));
00578 // chromaClr = Colour(x / y, 1.0, (1.0 - x - y) / y); 
00579 clip = MinCmpt(csDisplay.fromXYZ * chromaClr) < 0.0;
00580 if (clip != lastClip)
00581 {
00582 image.SetPixel(i, j, cYellow);
00583 lastClip = clip;
00584 }
00585 else
00586 {
00587 chromaClr = csDisplay.XYZToGamutClip(chromaClr);
00588 image.SetPixel(i, j, chromaClr);
00589 }
00590 }
00591 }
00592 }
00593 
00594  Colour CIEXYZToLuv(const Colour &cXYZ, const Colour &cWhiteXYZ)
00595 {
00596 ClrReal lumRatio = cXYZ[1] / cWhiteXYZ[1];
00597 const Colour denomRat(1.0, 15.0, 3.0);
00598 ClrReal denom, u, v, lumStar;
00599 Colour result;
00600 
00601 if (lumRatio <= 0.008856)
00602 lumStar = 903.3 * lumRatio;
00603 else
00604 lumStar = 116.0 * pow(lumRatio, 1.0 / 3.0) - 16.0;
00605 
00606 denom = dot(denomRat, cXYZ);
00607 u = 4.0 * cXYZ[0] / denom;
00608 v = 9.0 * cXYZ[1] / denom;
00609 denom = dot(denomRat, cWhiteXYZ);
00610 u -= 4.0 * cWhiteXYZ[0] / denom;
00611 v -= 9.0 * cWhiteXYZ[1] / denom;
00612 
00613 result[0] = lumStar;
00614 result[1] = 13.0 * lumStar * u;
00615 result[2] = 13.0 * lumStar * v;
00616 
00617 return(result);
00618 }
00619 
00620  BlackBodySpectrum::BlackBodySpectrum(Float temperature)
00621 {
00622 temp = temperature;
00623 }
00624 
00625  Float BlackBodySpectrum::Sample(Float lambda) const
00626 {
00627 Float wavelength = lambda * 1e-9;
00628 Float result;
00629 
00630 result = 3.74183e-16 * pow(wavelength, -5.0);
00631 result /= exp(1.4388e-2 / (wavelength * temp)) - 1.0;
00632 
00633 return(result);
00634 }