PBRT
/home/felix/UBC/projects/AdaptiveLightfieldSampling/pbrt_v2/src/core/spectrum.h
00001 
00002 /*
00003     pbrt source code Copyright(c) 1998-2012 Matt Pharr and Greg Humphreys.
00004 
00005     This file is part of pbrt.
00006 
00007     Redistribution and use in source and binary forms, with or without
00008     modification, are permitted provided that the following conditions are
00009     met:
00010 
00011     - Redistributions of source code must retain the above copyright
00012       notice, this list of conditions and the following disclaimer.
00013 
00014     - Redistributions in binary form must reproduce the above copyright
00015       notice, this list of conditions and the following disclaimer in the
00016       documentation and/or other materials provided with the distribution.
00017 
00018     THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
00019     IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
00020     TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
00021     PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
00022     HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
00023     SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
00024     LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
00025     DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
00026     THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
00027     (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
00028     OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00029 
00030  */
00031 
00032 #if defined(_MSC_VER)
00033 #pragma once
00034 #endif
00035 
00036 #ifndef PBRT_CORE_SPECTRUM_H
00037 #define PBRT_CORE_SPECTRUM_H
00038 
00039 // core/spectrum.h*
00040 #include "pbrt.h"
00041 #include "parallel.h"
00042 
00043 // Spectrum Utility Declarations
00044 static const int sampledLambdaStart = 400;
00045 static const int sampledLambdaEnd = 700;
00046 static const int nSpectralSamples = 30;
00047 extern bool SpectrumSamplesSorted(const float *lambda, const float *vals, int n);
00048 extern void SortSpectrumSamples(float *lambda, float *vals, int n);
00049 extern float AverageSpectrumSamples(const float *lambda, const float *vals,
00050     int n, float lambdaStart, float lambdaEnd);
00051 inline void XYZToRGB(const float xyz[3], float rgb[3]) {
00052     rgb[0] =  3.240479f*xyz[0] - 1.537150f*xyz[1] - 0.498535f*xyz[2];
00053     rgb[1] = -0.969256f*xyz[0] + 1.875991f*xyz[1] + 0.041556f*xyz[2];
00054     rgb[2] =  0.055648f*xyz[0] - 0.204043f*xyz[1] + 1.057311f*xyz[2];
00055 }
00056 
00057 
00058 inline void RGBToXYZ(const float rgb[3], float xyz[3]) {
00059     xyz[0] = 0.412453f*rgb[0] + 0.357580f*rgb[1] + 0.180423f*rgb[2];
00060     xyz[1] = 0.212671f*rgb[0] + 0.715160f*rgb[1] + 0.072169f*rgb[2];
00061     xyz[2] = 0.019334f*rgb[0] + 0.119193f*rgb[1] + 0.950227f*rgb[2];
00062 }
00063 
00064 
00065 enum SpectrumType { SPECTRUM_REFLECTANCE, SPECTRUM_ILLUMINANT };
00066 extern void Blackbody(const float *wl, int n, float temp, float *vals);
00067 extern float InterpolateSpectrumSamples(const float *lambda, const float *vals,
00068                                         int n, float l);
00069 
00070 // Spectral Data Declarations
00071 static const int nCIESamples = 471;
00072 extern const float CIE_X[nCIESamples];
00073 extern const float CIE_Y[nCIESamples];
00074 extern const float CIE_Z[nCIESamples];
00075 extern const float CIE_lambda[nCIESamples];
00076 static const int nRGB2SpectSamples = 32;
00077 extern const float RGB2SpectLambda[nRGB2SpectSamples];
00078 extern const float RGBRefl2SpectWhite[nRGB2SpectSamples];
00079 extern const float RGBRefl2SpectCyan[nRGB2SpectSamples];
00080 extern const float RGBRefl2SpectMagenta[nRGB2SpectSamples];
00081 extern const float RGBRefl2SpectYellow[nRGB2SpectSamples];
00082 extern const float RGBRefl2SpectRed[nRGB2SpectSamples];
00083 extern const float RGBRefl2SpectGreen[nRGB2SpectSamples];
00084 extern const float RGBRefl2SpectBlue[nRGB2SpectSamples];
00085 extern const float RGBIllum2SpectWhite[nRGB2SpectSamples];
00086 extern const float RGBIllum2SpectCyan[nRGB2SpectSamples];
00087 extern const float RGBIllum2SpectMagenta[nRGB2SpectSamples];
00088 extern const float RGBIllum2SpectYellow[nRGB2SpectSamples];
00089 extern const float RGBIllum2SpectRed[nRGB2SpectSamples];
00090 extern const float RGBIllum2SpectGreen[nRGB2SpectSamples];
00091 extern const float RGBIllum2SpectBlue[nRGB2SpectSamples];
00092 
00093 // Spectrum Declarations
00094 template <int nSamples> class CoefficientSpectrum {
00095 public:
00096     // CoefficientSpectrum Public Methods
00097     CoefficientSpectrum(float v = 0.f) {
00098         for (int i = 0; i < nSamples; ++i)
00099             c[i] = v;
00100         Assert(!HasNaNs());
00101     }
00102 #ifdef DEBUG
00103     CoefficientSpectrum(const CoefficientSpectrum &s) {
00104         Assert(!s.HasNaNs());
00105         for (int i = 0; i < nSamples; ++i)
00106             c[i] = s.c[i];
00107     }
00108     
00109     CoefficientSpectrum &operator=(const CoefficientSpectrum &s) {
00110         Assert(!s.HasNaNs());
00111         for (int i = 0; i < nSamples; ++i)
00112             c[i] = s.c[i];
00113         return *this;
00114     }
00115 #endif // DEBUG
00116     void Print(FILE *f) const {
00117         fprintf(f, "[ ");
00118         for (int i = 0; i < nSamples; ++i) {
00119             fprintf(f, "%f", c[i]);
00120             if (i != nSamples-1) fprintf(f, ", ");
00121         }
00122         fprintf(f, "]");
00123     }
00124     CoefficientSpectrum &operator+=(const CoefficientSpectrum &s2) {
00125         Assert(!s2.HasNaNs());
00126         for (int i = 0; i < nSamples; ++i)
00127             c[i] += s2.c[i];
00128         return *this;
00129     }
00130     CoefficientSpectrum operator+(const CoefficientSpectrum &s2) const {
00131         Assert(!s2.HasNaNs());
00132         CoefficientSpectrum ret = *this;
00133         for (int i = 0; i < nSamples; ++i)
00134             ret.c[i] += s2.c[i];
00135         return ret;
00136     }
00137     CoefficientSpectrum operator-(const CoefficientSpectrum &s2) const {
00138         Assert(!s2.HasNaNs());
00139         CoefficientSpectrum ret = *this;
00140         for (int i = 0; i < nSamples; ++i)
00141             ret.c[i] -= s2.c[i];
00142         return ret;
00143     }
00144     CoefficientSpectrum operator/(const CoefficientSpectrum &s2) const {
00145         Assert(!s2.HasNaNs());
00146         CoefficientSpectrum ret = *this;
00147         for (int i = 0; i < nSamples; ++i)
00148             ret.c[i] /= s2.c[i];
00149         return ret;
00150     }
00151     CoefficientSpectrum operator*(const CoefficientSpectrum &sp) const {
00152         Assert(!sp.HasNaNs());
00153         CoefficientSpectrum ret = *this;
00154         for (int i = 0; i < nSamples; ++i)
00155             ret.c[i] *= sp.c[i];
00156         return ret;
00157     }
00158     CoefficientSpectrum &operator*=(const CoefficientSpectrum &sp) {
00159         Assert(!sp.HasNaNs());
00160         for (int i = 0; i < nSamples; ++i)
00161             c[i] *= sp.c[i];
00162         return *this;
00163     }
00164     CoefficientSpectrum operator*(float a) const {
00165         CoefficientSpectrum ret = *this;
00166         for (int i = 0; i < nSamples; ++i)
00167             ret.c[i] *= a;
00168         Assert(!ret.HasNaNs());
00169         return ret;
00170     }
00171     CoefficientSpectrum &operator*=(float a) {
00172         for (int i = 0; i < nSamples; ++i)
00173             c[i] *= a;
00174         Assert(!HasNaNs());
00175         return *this;
00176     }
00177     friend inline
00178     CoefficientSpectrum operator*(float a, const CoefficientSpectrum &s) {
00179         Assert(!isnan(a) && !s.HasNaNs());
00180         return s * a;
00181     }
00182     CoefficientSpectrum operator/(float a) const {
00183         Assert(!isnan(a));
00184         CoefficientSpectrum ret = *this;
00185         for (int i = 0; i < nSamples; ++i)
00186             ret.c[i] /= a;
00187         Assert(!ret.HasNaNs());
00188         return ret;
00189     }
00190     CoefficientSpectrum &operator/=(float a) {
00191         Assert(!isnan(a));
00192         for (int i = 0; i < nSamples; ++i)
00193             c[i] /= a;
00194         return *this;
00195     }
00196     bool operator==(const CoefficientSpectrum &sp) const {
00197         for (int i = 0; i < nSamples; ++i)
00198             if (c[i] != sp.c[i]) return false;
00199         return true;
00200     }
00201     bool operator!=(const CoefficientSpectrum &sp) const {
00202         return !(*this == sp);
00203     }
00204     bool IsBlack() const {
00205         for (int i = 0; i < nSamples; ++i)
00206             if (c[i] != 0.) return false;
00207         return true;
00208     }
00209     friend CoefficientSpectrum Sqrt(const CoefficientSpectrum &s) {
00210         CoefficientSpectrum ret;
00211         for (int i = 0; i < nSamples; ++i)
00212             ret.c[i] = sqrtf(s.c[i]);
00213         Assert(!ret.HasNaNs());
00214         return ret;
00215     }
00216     template <int n> friend inline CoefficientSpectrum<n> Pow(const CoefficientSpectrum<n> &s, float e);
00217     CoefficientSpectrum operator-() const {
00218         CoefficientSpectrum ret;
00219         for (int i = 0; i < nSamples; ++i)
00220             ret.c[i] = -c[i];
00221         return ret;
00222     }
00223     friend CoefficientSpectrum Exp(const CoefficientSpectrum &s) {
00224         CoefficientSpectrum ret;
00225         for (int i = 0; i < nSamples; ++i)
00226             ret.c[i] = expf(s.c[i]);
00227         Assert(!ret.HasNaNs());
00228         return ret;
00229     }
00230     CoefficientSpectrum Clamp(float low = 0, float high = INFINITY) const {
00231         CoefficientSpectrum ret;
00232         for (int i = 0; i < nSamples; ++i)
00233             ret.c[i] = ::Clamp(c[i], low, high);
00234         Assert(!ret.HasNaNs());
00235         return ret;
00236     }
00237     bool HasNaNs() const {
00238         for (int i = 0; i < nSamples; ++i)
00239             if (isnan(c[i])) return true;
00240         return false;
00241     }
00242     bool Write(FILE *f) const {
00243         for (int i = 0; i < nSamples; ++i)
00244             if (fprintf(f, "%f ", c[i]) < 0) return false;
00245         return true;
00246     }
00247     bool Read(FILE *f) {
00248         for (int i = 0; i < nSamples; ++i)
00249             if (fscanf(f, "%f ", &c[i]) != 1) return false;
00250         return true;
00251     }
00252 protected:
00253     // CoefficientSpectrum Protected Data
00254     float c[nSamples];
00255 };
00256 
00257 
00258 class SampledSpectrum : public CoefficientSpectrum<nSpectralSamples> {
00259 public:
00260     // SampledSpectrum Public Methods
00261     SampledSpectrum(float v = 0.f) {
00262         for (int i = 0; i < nSpectralSamples; ++i) c[i] = v;
00263     }
00264     SampledSpectrum(const CoefficientSpectrum<nSpectralSamples> &v)
00265         : CoefficientSpectrum<nSpectralSamples>(v) { }
00266     static SampledSpectrum FromSampled(const float *lambda,
00267                                        const float *v, int n) {
00268         // Sort samples if unordered, use sorted for returned spectrum
00269         if (!SpectrumSamplesSorted(lambda, v, n)) {
00270             vector<float> slambda(&lambda[0], &lambda[n]);
00271             vector<float> sv(&v[0], &v[n]);
00272             SortSpectrumSamples(&slambda[0], &sv[0], n);
00273             return FromSampled(&slambda[0], &sv[0], n);
00274         }
00275         SampledSpectrum r;
00276         for (int i = 0; i < nSpectralSamples; ++i) {
00277             // Compute average value of given SPD over $i$th sample's range
00278             float lambda0 = Lerp(float(i) / float(nSpectralSamples),
00279                                  sampledLambdaStart, sampledLambdaEnd);
00280             float lambda1 = Lerp(float(i+1) / float(nSpectralSamples),
00281                                  sampledLambdaStart, sampledLambdaEnd);
00282             r.c[i] = AverageSpectrumSamples(lambda, v, n, lambda0, lambda1);
00283         }
00284         return r;
00285     }
00286     static void Init() {
00287         // Compute XYZ matching functions for _SampledSpectrum_
00288         for (int i = 0; i < nSpectralSamples; ++i) {
00289             float wl0 = Lerp(float(i) / float(nSpectralSamples),
00290                              sampledLambdaStart, sampledLambdaEnd);
00291             float wl1 = Lerp(float(i+1) / float(nSpectralSamples),
00292                              sampledLambdaStart, sampledLambdaEnd);
00293             X.c[i] = AverageSpectrumSamples(CIE_lambda, CIE_X, nCIESamples,
00294                                             wl0, wl1);
00295             Y.c[i] = AverageSpectrumSamples(CIE_lambda, CIE_Y, nCIESamples,
00296                                             wl0, wl1);
00297             Z.c[i] = AverageSpectrumSamples(CIE_lambda, CIE_Z, nCIESamples,
00298                                             wl0, wl1);
00299             yint += Y.c[i];
00300         }
00301 
00302         // Compute RGB to spectrum functions for _SampledSpectrum_
00303         for (int i = 0; i < nSpectralSamples; ++i) {
00304             float wl0 = Lerp(float(i) / float(nSpectralSamples),
00305                              sampledLambdaStart, sampledLambdaEnd);
00306             float wl1 = Lerp(float(i+1) / float(nSpectralSamples),
00307                              sampledLambdaStart, sampledLambdaEnd);
00308             rgbRefl2SpectWhite.c[i] = AverageSpectrumSamples(RGB2SpectLambda, RGBRefl2SpectWhite,
00309                 nRGB2SpectSamples, wl0, wl1);
00310             rgbRefl2SpectCyan.c[i] = AverageSpectrumSamples(RGB2SpectLambda, RGBRefl2SpectCyan,
00311                 nRGB2SpectSamples, wl0, wl1);
00312             rgbRefl2SpectMagenta.c[i] = AverageSpectrumSamples(RGB2SpectLambda, RGBRefl2SpectMagenta,
00313                 nRGB2SpectSamples, wl0, wl1);
00314             rgbRefl2SpectYellow.c[i] = AverageSpectrumSamples(RGB2SpectLambda, RGBRefl2SpectYellow,
00315                 nRGB2SpectSamples, wl0, wl1);
00316             rgbRefl2SpectRed.c[i] = AverageSpectrumSamples(RGB2SpectLambda, RGBRefl2SpectRed,
00317                 nRGB2SpectSamples, wl0, wl1);
00318             rgbRefl2SpectGreen.c[i] = AverageSpectrumSamples(RGB2SpectLambda, RGBRefl2SpectGreen,
00319                 nRGB2SpectSamples, wl0, wl1);
00320             rgbRefl2SpectBlue.c[i] = AverageSpectrumSamples(RGB2SpectLambda, RGBRefl2SpectBlue,
00321                 nRGB2SpectSamples, wl0, wl1);
00322         
00323             rgbIllum2SpectWhite.c[i] = AverageSpectrumSamples(RGB2SpectLambda, RGBIllum2SpectWhite,
00324                 nRGB2SpectSamples, wl0, wl1);
00325             rgbIllum2SpectCyan.c[i] = AverageSpectrumSamples(RGB2SpectLambda, RGBIllum2SpectCyan,
00326                 nRGB2SpectSamples, wl0, wl1);
00327             rgbIllum2SpectMagenta.c[i] = AverageSpectrumSamples(RGB2SpectLambda, RGBIllum2SpectMagenta,
00328                 nRGB2SpectSamples, wl0, wl1);
00329             rgbIllum2SpectYellow.c[i] = AverageSpectrumSamples(RGB2SpectLambda, RGBIllum2SpectYellow,
00330                 nRGB2SpectSamples, wl0, wl1);
00331             rgbIllum2SpectRed.c[i] = AverageSpectrumSamples(RGB2SpectLambda, RGBIllum2SpectRed,
00332                 nRGB2SpectSamples, wl0, wl1);
00333             rgbIllum2SpectGreen.c[i] = AverageSpectrumSamples(RGB2SpectLambda, RGBIllum2SpectGreen,
00334                 nRGB2SpectSamples, wl0, wl1);
00335             rgbIllum2SpectBlue.c[i] = AverageSpectrumSamples(RGB2SpectLambda, RGBIllum2SpectBlue,
00336                 nRGB2SpectSamples, wl0, wl1);
00337         }
00338     }
00339     void ToXYZ(float xyz[3]) const {
00340         xyz[0] = xyz[1] = xyz[2] = 0.f;
00341         for (int i = 0; i < nSpectralSamples; ++i) {
00342             xyz[0] += X.c[i] * c[i];
00343             xyz[1] += Y.c[i] * c[i];
00344             xyz[2] += Z.c[i] * c[i];
00345         }
00346         xyz[0] /= yint;
00347         xyz[1] /= yint;
00348         xyz[2] /= yint;
00349     }
00350     float y() const {
00351         float yy = 0.f;
00352         for (int i = 0; i < nSpectralSamples; ++i)
00353             yy += Y.c[i] * c[i];
00354         return yy / yint;
00355     }
00356     void ToRGB(float rgb[3]) const {
00357         float xyz[3];
00358         ToXYZ(xyz);
00359         XYZToRGB(xyz, rgb);
00360     }
00361     RGBSpectrum ToRGBSpectrum() const;
00362     static SampledSpectrum FromRGB(const float rgb[3],
00363         SpectrumType type = SPECTRUM_REFLECTANCE);
00364     static SampledSpectrum FromXYZ(const float xyz[3],
00365             SpectrumType type = SPECTRUM_REFLECTANCE) {
00366         float rgb[3];
00367         XYZToRGB(xyz, rgb);
00368         return FromRGB(rgb, type);
00369     }
00370     SampledSpectrum(const RGBSpectrum &r, SpectrumType type = SPECTRUM_REFLECTANCE);
00371 private:
00372     // SampledSpectrum Private Data
00373     static SampledSpectrum X, Y, Z;
00374     static float yint;
00375     static SampledSpectrum rgbRefl2SpectWhite, rgbRefl2SpectCyan;
00376     static SampledSpectrum rgbRefl2SpectMagenta, rgbRefl2SpectYellow;
00377     static SampledSpectrum rgbRefl2SpectRed, rgbRefl2SpectGreen;
00378     static SampledSpectrum rgbRefl2SpectBlue;
00379     static SampledSpectrum rgbIllum2SpectWhite, rgbIllum2SpectCyan;
00380     static SampledSpectrum rgbIllum2SpectMagenta, rgbIllum2SpectYellow;
00381     static SampledSpectrum rgbIllum2SpectRed, rgbIllum2SpectGreen;
00382     static SampledSpectrum rgbIllum2SpectBlue;
00383 };
00384 
00385 
00386 class RGBSpectrum : public CoefficientSpectrum<3> {
00387     using CoefficientSpectrum<3>::c;
00388 public:
00389     // RGBSpectrum Public Methods
00390     RGBSpectrum(float v = 0.f) : CoefficientSpectrum<3>(v) { }
00391     RGBSpectrum(const CoefficientSpectrum<3> &v)
00392         : CoefficientSpectrum<3>(v) { }
00393     RGBSpectrum(const RGBSpectrum &s, SpectrumType type = SPECTRUM_REFLECTANCE) {
00394         *this = s;
00395     }
00396     static RGBSpectrum FromRGB(const float rgb[3],
00397             SpectrumType type = SPECTRUM_REFLECTANCE) {
00398         RGBSpectrum s;
00399         s.c[0] = rgb[0];
00400         s.c[1] = rgb[1];
00401         s.c[2] = rgb[2];
00402         Assert(!s.HasNaNs());
00403         return s;
00404     }
00405     void ToRGB(float *rgb) const {
00406         rgb[0] = c[0];
00407         rgb[1] = c[1];
00408         rgb[2] = c[2];
00409     }
00410     const RGBSpectrum &ToRGBSpectrum() const {
00411         return *this;
00412     }
00413     void ToXYZ(float xyz[3]) const {
00414         RGBToXYZ(c, xyz);
00415     }
00416     static RGBSpectrum FromXYZ(const float xyz[3],
00417             SpectrumType type = SPECTRUM_REFLECTANCE) {
00418         RGBSpectrum r;
00419         XYZToRGB(xyz, r.c);
00420         return r;
00421     }
00422     float y() const {
00423         const float YWeight[3] = { 0.212671f, 0.715160f, 0.072169f };
00424         return YWeight[0] * c[0] + YWeight[1] * c[1] + YWeight[2] * c[2];
00425     }
00426     static RGBSpectrum FromSampled(const float *lambda, const float *v,
00427                                    int n) {
00428         // Sort samples if unordered, use sorted for returned spectrum
00429         if (!SpectrumSamplesSorted(lambda, v, n)) {
00430             vector<float> slambda(&lambda[0], &lambda[n]);
00431             vector<float> sv(&v[0], &v[n]);
00432             SortSpectrumSamples(&slambda[0], &sv[0], n);
00433             return FromSampled(&slambda[0], &sv[0], n);
00434         }
00435         float xyz[3] = { 0, 0, 0 };
00436         float yint = 0.f;
00437         for (int i = 0; i < nCIESamples; ++i) {
00438             yint += CIE_Y[i];
00439             float val = InterpolateSpectrumSamples(lambda, v, n,
00440                                                    CIE_lambda[i]);
00441             xyz[0] += val * CIE_X[i];
00442             xyz[1] += val * CIE_Y[i];
00443             xyz[2] += val * CIE_Z[i];
00444         }
00445         xyz[0] /= yint;
00446         xyz[1] /= yint;
00447         xyz[2] /= yint;
00448         return FromXYZ(xyz);
00449     }
00450 };
00451 
00452 
00453 
00454 // Spectrum Inline Functions
00455 template <int nSamples> inline CoefficientSpectrum<nSamples>
00456 Pow(const CoefficientSpectrum<nSamples> &s, float e) {
00457     CoefficientSpectrum<nSamples> ret;
00458     for (int i = 0; i < nSamples; ++i)
00459         ret.c[i] = powf(s.c[i], e);
00460     Assert(!ret.HasNaNs());
00461     return ret;
00462 }
00463 
00464 
00465 inline Spectrum Lerp(float t, const Spectrum &s1, const Spectrum &s2) {
00466     return (1.f - t) * s1 + t * s2;
00467 }
00468 
00469 
00470 
00471 #endif // PBRT_CORE_SPECTRUM_H