00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024 #if defined(_MSC_VER)
00025 #pragma once
00026 #endif
00027
00028 #ifndef PBRT_CORE_SPECTRUM_H
00029 #define PBRT_CORE_SPECTRUM_H
00030
00031
00032 #include "pbrt.h"
00033 #include "parallel.h"
00034
00035
00036 static const int sampledLambdaStart = 400;
00037 static const int sampledLambdaEnd = 700;
00038 static const int nSpectralSamples = 30;
00039 extern bool SpectrumSamplesSorted(const float *lambda, const float *vals, int n);
00040 extern void SortSpectrumSamples(float *lambda, float *vals, int n);
00041 extern float AverageSpectrumSamples(const float *lambda, const float *vals,
00042 int n, float lambdaStart, float lambdaEnd);
00043 inline void XYZToRGB(const float xyz[3], float rgb[3]) {
00044 rgb[0] = 3.240479f*xyz[0] - 1.537150f*xyz[1] - 0.498535f*xyz[2];
00045 rgb[1] = -0.969256f*xyz[0] + 1.875991f*xyz[1] + 0.041556f*xyz[2];
00046 rgb[2] = 0.055648f*xyz[0] - 0.204043f*xyz[1] + 1.057311f*xyz[2];
00047 }
00048
00049
00050 inline void RGBToXYZ(const float rgb[3], float xyz[3]) {
00051 xyz[0] = 0.412453f*rgb[0] + 0.357580f*rgb[1] + 0.180423f*rgb[2];
00052 xyz[1] = 0.212671f*rgb[0] + 0.715160f*rgb[1] + 0.072169f*rgb[2];
00053 xyz[2] = 0.019334f*rgb[0] + 0.119193f*rgb[1] + 0.950227f*rgb[2];
00054 }
00055
00056
00057 enum SpectrumType { SPECTRUM_REFLECTANCE, SPECTRUM_ILLUMINANT };
00058 extern void Blackbody(const float *wl, int n, float temp, float *vals);
00059 extern float InterpolateSpectrumSamples(const float *lambda, const float *vals,
00060 int n, float l);
00061
00062
00063 static const int nCIESamples = 471;
00064 extern const float CIE_X[nCIESamples];
00065 extern const float CIE_Y[nCIESamples];
00066 extern const float CIE_Z[nCIESamples];
00067 extern const float CIE_lambda[nCIESamples];
00068 static const int nRGB2SpectSamples = 32;
00069 extern const float RGB2SpectLambda[nRGB2SpectSamples];
00070 extern const float RGBRefl2SpectWhite[nRGB2SpectSamples];
00071 extern const float RGBRefl2SpectCyan[nRGB2SpectSamples];
00072 extern const float RGBRefl2SpectMagenta[nRGB2SpectSamples];
00073 extern const float RGBRefl2SpectYellow[nRGB2SpectSamples];
00074 extern const float RGBRefl2SpectRed[nRGB2SpectSamples];
00075 extern const float RGBRefl2SpectGreen[nRGB2SpectSamples];
00076 extern const float RGBRefl2SpectBlue[nRGB2SpectSamples];
00077 extern const float RGBIllum2SpectWhite[nRGB2SpectSamples];
00078 extern const float RGBIllum2SpectCyan[nRGB2SpectSamples];
00079 extern const float RGBIllum2SpectMagenta[nRGB2SpectSamples];
00080 extern const float RGBIllum2SpectYellow[nRGB2SpectSamples];
00081 extern const float RGBIllum2SpectRed[nRGB2SpectSamples];
00082 extern const float RGBIllum2SpectGreen[nRGB2SpectSamples];
00083 extern const float RGBIllum2SpectBlue[nRGB2SpectSamples];
00084
00085
00086 template <int nSamples> class CoefficientSpectrum {
00087 public:
00088
00089 CoefficientSpectrum(float v = 0.f) {
00090 for (int i = 0; i < nSamples; ++i)
00091 c[i] = v;
00092 Assert(!HasNaNs());
00093 }
00094 #ifdef DEBUG
00095 CoefficientSpectrum(const CoefficientSpectrum &s) {
00096 Assert(!s.HasNaNs());
00097 for (int i = 0; i < nSamples; ++i)
00098 c[i] = s.c[i];
00099 }
00100
00101 CoefficientSpectrum &operator=(const CoefficientSpectrum &s) {
00102 Assert(!s.HasNaNs());
00103 for (int i = 0; i < nSamples; ++i)
00104 c[i] = s.c[i];
00105 return *this;
00106 }
00107 #endif // DEBUG
00108 void Print(FILE *f) const {
00109 fprintf(f, "[ ");
00110 for (int i = 0; i < nSamples; ++i) {
00111 fprintf(f, "%f", c[i]);
00112 if (i != nSamples-1) fprintf(f, ", ");
00113 }
00114 fprintf(f, "]");
00115 }
00116 CoefficientSpectrum &operator+=(const CoefficientSpectrum &s2) {
00117 Assert(!s2.HasNaNs());
00118 for (int i = 0; i < nSamples; ++i)
00119 c[i] += s2.c[i];
00120 return *this;
00121 }
00122 CoefficientSpectrum operator+(const CoefficientSpectrum &s2) const {
00123 Assert(!s2.HasNaNs());
00124 CoefficientSpectrum ret = *this;
00125 for (int i = 0; i < nSamples; ++i)
00126 ret.c[i] += s2.c[i];
00127 return ret;
00128 }
00129 CoefficientSpectrum operator-(const CoefficientSpectrum &s2) const {
00130 Assert(!s2.HasNaNs());
00131 CoefficientSpectrum ret = *this;
00132 for (int i = 0; i < nSamples; ++i)
00133 ret.c[i] -= s2.c[i];
00134 return ret;
00135 }
00136 CoefficientSpectrum operator/(const CoefficientSpectrum &s2) const {
00137 Assert(!s2.HasNaNs());
00138 CoefficientSpectrum ret = *this;
00139 for (int i = 0; i < nSamples; ++i)
00140 ret.c[i] /= s2.c[i];
00141 return ret;
00142 }
00143 CoefficientSpectrum operator*(const CoefficientSpectrum &sp) const {
00144 Assert(!sp.HasNaNs());
00145 CoefficientSpectrum ret = *this;
00146 for (int i = 0; i < nSamples; ++i)
00147 ret.c[i] *= sp.c[i];
00148 return ret;
00149 }
00150 CoefficientSpectrum &operator*=(const CoefficientSpectrum &sp) {
00151 Assert(!sp.HasNaNs());
00152 for (int i = 0; i < nSamples; ++i)
00153 c[i] *= sp.c[i];
00154 return *this;
00155 }
00156 CoefficientSpectrum operator*(float a) const {
00157 CoefficientSpectrum ret = *this;
00158 for (int i = 0; i < nSamples; ++i)
00159 ret.c[i] *= a;
00160 Assert(!ret.HasNaNs());
00161 return ret;
00162 }
00163 CoefficientSpectrum &operator*=(float a) {
00164 for (int i = 0; i < nSamples; ++i)
00165 c[i] *= a;
00166 Assert(!HasNaNs());
00167 return *this;
00168 }
00169 friend inline
00170 CoefficientSpectrum operator*(float a, const CoefficientSpectrum &s) {
00171 Assert(!isnan(a) && !s.HasNaNs());
00172 return s * a;
00173 }
00174 CoefficientSpectrum operator/(float a) const {
00175 Assert(!isnan(a));
00176 CoefficientSpectrum ret = *this;
00177 for (int i = 0; i < nSamples; ++i)
00178 ret.c[i] /= a;
00179 Assert(!ret.HasNaNs());
00180 return ret;
00181 }
00182 CoefficientSpectrum &operator/=(float a) {
00183 Assert(!isnan(a));
00184 for (int i = 0; i < nSamples; ++i)
00185 c[i] /= a;
00186 return *this;
00187 }
00188 bool operator==(const CoefficientSpectrum &sp) const {
00189 for (int i = 0; i < nSamples; ++i)
00190 if (c[i] != sp.c[i]) return false;
00191 return true;
00192 }
00193 bool operator!=(const CoefficientSpectrum &sp) const {
00194 return !(*this == sp);
00195 }
00196 bool IsBlack() const {
00197 for (int i = 0; i < nSamples; ++i)
00198 if (c[i] != 0.) return false;
00199 return true;
00200 }
00201 friend CoefficientSpectrum Sqrt(const CoefficientSpectrum &s) {
00202 CoefficientSpectrum ret;
00203 for (int i = 0; i < nSamples; ++i)
00204 ret.c[i] = sqrtf(s.c[i]);
00205 Assert(!ret.HasNaNs());
00206 return ret;
00207 }
00208 template <int n> friend inline CoefficientSpectrum<n> Pow(const CoefficientSpectrum<n> &s, float e);
00209 CoefficientSpectrum operator-() const {
00210 CoefficientSpectrum ret;
00211 for (int i = 0; i < nSamples; ++i)
00212 ret.c[i] = -c[i];
00213 return ret;
00214 }
00215 friend CoefficientSpectrum Exp(const CoefficientSpectrum &s) {
00216 CoefficientSpectrum ret;
00217 for (int i = 0; i < nSamples; ++i)
00218 ret.c[i] = expf(s.c[i]);
00219 Assert(!ret.HasNaNs());
00220 return ret;
00221 }
00222 CoefficientSpectrum Clamp(float low = 0, float high = INFINITY) const {
00223 CoefficientSpectrum ret;
00224 for (int i = 0; i < nSamples; ++i)
00225 ret.c[i] = ::Clamp(c[i], low, high);
00226 Assert(!ret.HasNaNs());
00227 return ret;
00228 }
00229 bool HasNaNs() const {
00230 for (int i = 0; i < nSamples; ++i)
00231 if (isnan(c[i])) return true;
00232 return false;
00233 }
00234 bool Write(FILE *f) const {
00235 for (int i = 0; i < nSamples; ++i)
00236 if (fprintf(f, "%f ", c[i]) < 0) return false;
00237 return true;
00238 }
00239 bool Read(FILE *f) {
00240 for (int i = 0; i < nSamples; ++i)
00241 if (fscanf(f, "%f ", &c[i]) != 1) return false;
00242 return true;
00243 }
00244 protected:
00245
00246 float c[nSamples];
00247 };
00248
00249
00250 class SampledSpectrum : public CoefficientSpectrum<nSpectralSamples> {
00251 public:
00252
00253 SampledSpectrum(float v = 0.f) {
00254 for (int i = 0; i < nSpectralSamples; ++i) c[i] = v;
00255 }
00256 SampledSpectrum(const CoefficientSpectrum<nSpectralSamples> &v)
00257 : CoefficientSpectrum<nSpectralSamples>(v) { }
00258 static SampledSpectrum FromSampled(const float *lambda,
00259 const float *v, int n) {
00260
00261 if (!SpectrumSamplesSorted(lambda, v, n)) {
00262 vector<float> slambda(&lambda[0], &lambda[n]);
00263 vector<float> sv(&v[0], &v[n]);
00264 SortSpectrumSamples(&slambda[0], &sv[0], n);
00265 return FromSampled(&slambda[0], &sv[0], n);
00266 }
00267 SampledSpectrum r;
00268 for (int i = 0; i < nSpectralSamples; ++i) {
00269
00270 float lambda0 = Lerp(float(i) / float(nSpectralSamples),
00271 sampledLambdaStart, sampledLambdaEnd);
00272 float lambda1 = Lerp(float(i+1) / float(nSpectralSamples),
00273 sampledLambdaStart, sampledLambdaEnd);
00274 r.c[i] = AverageSpectrumSamples(lambda, v, n, lambda0, lambda1);
00275 }
00276 return r;
00277 }
00278 static void Init() {
00279
00280 for (int i = 0; i < nSpectralSamples; ++i) {
00281 float wl0 = Lerp(float(i) / float(nSpectralSamples),
00282 sampledLambdaStart, sampledLambdaEnd);
00283 float wl1 = Lerp(float(i+1) / float(nSpectralSamples),
00284 sampledLambdaStart, sampledLambdaEnd);
00285 X.c[i] = AverageSpectrumSamples(CIE_lambda, CIE_X, nCIESamples,
00286 wl0, wl1);
00287 Y.c[i] = AverageSpectrumSamples(CIE_lambda, CIE_Y, nCIESamples,
00288 wl0, wl1);
00289 Z.c[i] = AverageSpectrumSamples(CIE_lambda, CIE_Z, nCIESamples,
00290 wl0, wl1);
00291 yint += Y.c[i];
00292 }
00293
00294
00295 for (int i = 0; i < nSpectralSamples; ++i) {
00296 float wl0 = Lerp(float(i) / float(nSpectralSamples),
00297 sampledLambdaStart, sampledLambdaEnd);
00298 float wl1 = Lerp(float(i+1) / float(nSpectralSamples),
00299 sampledLambdaStart, sampledLambdaEnd);
00300 rgbRefl2SpectWhite.c[i] = AverageSpectrumSamples(RGB2SpectLambda, RGBRefl2SpectWhite,
00301 nRGB2SpectSamples, wl0, wl1);
00302 rgbRefl2SpectCyan.c[i] = AverageSpectrumSamples(RGB2SpectLambda, RGBRefl2SpectCyan,
00303 nRGB2SpectSamples, wl0, wl1);
00304 rgbRefl2SpectMagenta.c[i] = AverageSpectrumSamples(RGB2SpectLambda, RGBRefl2SpectMagenta,
00305 nRGB2SpectSamples, wl0, wl1);
00306 rgbRefl2SpectYellow.c[i] = AverageSpectrumSamples(RGB2SpectLambda, RGBRefl2SpectYellow,
00307 nRGB2SpectSamples, wl0, wl1);
00308 rgbRefl2SpectRed.c[i] = AverageSpectrumSamples(RGB2SpectLambda, RGBRefl2SpectRed,
00309 nRGB2SpectSamples, wl0, wl1);
00310 rgbRefl2SpectGreen.c[i] = AverageSpectrumSamples(RGB2SpectLambda, RGBRefl2SpectGreen,
00311 nRGB2SpectSamples, wl0, wl1);
00312 rgbRefl2SpectBlue.c[i] = AverageSpectrumSamples(RGB2SpectLambda, RGBRefl2SpectBlue,
00313 nRGB2SpectSamples, wl0, wl1);
00314
00315 rgbIllum2SpectWhite.c[i] = AverageSpectrumSamples(RGB2SpectLambda, RGBIllum2SpectWhite,
00316 nRGB2SpectSamples, wl0, wl1);
00317 rgbIllum2SpectCyan.c[i] = AverageSpectrumSamples(RGB2SpectLambda, RGBIllum2SpectCyan,
00318 nRGB2SpectSamples, wl0, wl1);
00319 rgbIllum2SpectMagenta.c[i] = AverageSpectrumSamples(RGB2SpectLambda, RGBIllum2SpectMagenta,
00320 nRGB2SpectSamples, wl0, wl1);
00321 rgbIllum2SpectYellow.c[i] = AverageSpectrumSamples(RGB2SpectLambda, RGBIllum2SpectYellow,
00322 nRGB2SpectSamples, wl0, wl1);
00323 rgbIllum2SpectRed.c[i] = AverageSpectrumSamples(RGB2SpectLambda, RGBIllum2SpectRed,
00324 nRGB2SpectSamples, wl0, wl1);
00325 rgbIllum2SpectGreen.c[i] = AverageSpectrumSamples(RGB2SpectLambda, RGBIllum2SpectGreen,
00326 nRGB2SpectSamples, wl0, wl1);
00327 rgbIllum2SpectBlue.c[i] = AverageSpectrumSamples(RGB2SpectLambda, RGBIllum2SpectBlue,
00328 nRGB2SpectSamples, wl0, wl1);
00329 }
00330 }
00331 void ToXYZ(float xyz[3]) const {
00332 xyz[0] = xyz[1] = xyz[2] = 0.f;
00333 for (int i = 0; i < nSpectralSamples; ++i) {
00334 xyz[0] += X.c[i] * c[i];
00335 xyz[1] += Y.c[i] * c[i];
00336 xyz[2] += Z.c[i] * c[i];
00337 }
00338 xyz[0] /= yint;
00339 xyz[1] /= yint;
00340 xyz[2] /= yint;
00341 }
00342 float y() const {
00343 float yy = 0.f;
00344 for (int i = 0; i < nSpectralSamples; ++i)
00345 yy += Y.c[i] * c[i];
00346 return yy / yint;
00347 }
00348 void ToRGB(float rgb[3]) const {
00349 float xyz[3];
00350 ToXYZ(xyz);
00351 XYZToRGB(xyz, rgb);
00352 }
00353 RGBSpectrum ToRGBSpectrum() const;
00354 static SampledSpectrum FromRGB(const float rgb[3],
00355 SpectrumType type = SPECTRUM_REFLECTANCE);
00356 static SampledSpectrum FromXYZ(const float xyz[3],
00357 SpectrumType type = SPECTRUM_REFLECTANCE) {
00358 float rgb[3];
00359 XYZToRGB(xyz, rgb);
00360 return FromRGB(rgb, type);
00361 }
00362 SampledSpectrum(const RGBSpectrum &r, SpectrumType type = SPECTRUM_REFLECTANCE);
00363 private:
00364
00365 static SampledSpectrum X, Y, Z;
00366 static float yint;
00367 static SampledSpectrum rgbRefl2SpectWhite, rgbRefl2SpectCyan;
00368 static SampledSpectrum rgbRefl2SpectMagenta, rgbRefl2SpectYellow;
00369 static SampledSpectrum rgbRefl2SpectRed, rgbRefl2SpectGreen;
00370 static SampledSpectrum rgbRefl2SpectBlue;
00371 static SampledSpectrum rgbIllum2SpectWhite, rgbIllum2SpectCyan;
00372 static SampledSpectrum rgbIllum2SpectMagenta, rgbIllum2SpectYellow;
00373 static SampledSpectrum rgbIllum2SpectRed, rgbIllum2SpectGreen;
00374 static SampledSpectrum rgbIllum2SpectBlue;
00375 };
00376
00377
00378 class RGBSpectrum : public CoefficientSpectrum<3> {
00379 using CoefficientSpectrum<3>::c;
00380 public:
00381
00382 RGBSpectrum(float v = 0.f) : CoefficientSpectrum<3>(v) { }
00383 RGBSpectrum(const CoefficientSpectrum<3> &v)
00384 : CoefficientSpectrum<3>(v) { }
00385 RGBSpectrum(const RGBSpectrum &s, SpectrumType type = SPECTRUM_REFLECTANCE) {
00386 *this = s;
00387 }
00388 static RGBSpectrum FromRGB(const float rgb[3],
00389 SpectrumType type = SPECTRUM_REFLECTANCE) {
00390 RGBSpectrum s;
00391 s.c[0] = rgb[0];
00392 s.c[1] = rgb[1];
00393 s.c[2] = rgb[2];
00394 Assert(!s.HasNaNs());
00395 return s;
00396 }
00397 void ToRGB(float *rgb) const {
00398 rgb[0] = c[0];
00399 rgb[1] = c[1];
00400 rgb[2] = c[2];
00401 }
00402 const RGBSpectrum &ToRGBSpectrum() const {
00403 return *this;
00404 }
00405 void ToXYZ(float xyz[3]) const {
00406 RGBToXYZ(c, xyz);
00407 }
00408 static RGBSpectrum FromXYZ(const float xyz[3],
00409 SpectrumType type = SPECTRUM_REFLECTANCE) {
00410 RGBSpectrum r;
00411 XYZToRGB(xyz, r.c);
00412 return r;
00413 }
00414 float y() const {
00415 const float YWeight[3] = { 0.212671f, 0.715160f, 0.072169f };
00416 return YWeight[0] * c[0] + YWeight[1] * c[1] + YWeight[2] * c[2];
00417 }
00418 static RGBSpectrum FromSampled(const float *lambda, const float *v,
00419 int n) {
00420
00421 if (!SpectrumSamplesSorted(lambda, v, n)) {
00422 vector<float> slambda(&lambda[0], &lambda[n]);
00423 vector<float> sv(&v[0], &v[n]);
00424 SortSpectrumSamples(&slambda[0], &sv[0], n);
00425 return FromSampled(&slambda[0], &sv[0], n);
00426 }
00427 float xyz[3] = { 0, 0, 0 };
00428 float yint = 0.f;
00429 for (int i = 0; i < nCIESamples; ++i) {
00430 yint += CIE_Y[i];
00431 float val = InterpolateSpectrumSamples(lambda, v, n,
00432 CIE_lambda[i]);
00433 xyz[0] += val * CIE_X[i];
00434 xyz[1] += val * CIE_Y[i];
00435 xyz[2] += val * CIE_Z[i];
00436 }
00437 xyz[0] /= yint;
00438 xyz[1] /= yint;
00439 xyz[2] /= yint;
00440 return FromXYZ(xyz);
00441 }
00442 };
00443
00444
00445
00446
00447 template <int nSamples> inline CoefficientSpectrum<nSamples>
00448 Pow(const CoefficientSpectrum<nSamples> &s, float e) {
00449 CoefficientSpectrum<nSamples> ret;
00450 for (int i = 0; i < nSamples; ++i)
00451 ret.c[i] = powf(s.c[i], e);
00452 Assert(!ret.HasNaNs());
00453 return ret;
00454 }
00455
00456
00457 inline Spectrum Lerp(float t, const Spectrum &s1, const Spectrum &s2) {
00458 return (1.f - t) * s1 + t * s2;
00459 }
00460
00461
00462
00463 #endif // PBRT_CORE_SPECTRUM_H