可见光 人眼的可见光波谱范围约为370nm~730nm
蓝光：400nm,绿光：550nm,红光：650nm
光谱表示 SPD(spectral power distribution): 光谱功率分布
SPD是一个关于波长的函数，通常是做法是找到一组基函数，用于将无限维度的SPD映射到低纬度空间表示
1 2 3 4 5 6 7 8 template <int nSamples>class CoefficientSpecturm { public : <CoefficientSpecturm public method> protected : <CoefficientSpecturm protected data> };
1 2 <CoefficientSpecturm protected data> = float c[nSamples];
接下来介绍两种常用的表示方法
SampleSpectrum 这种方法是将光谱重新均匀采样来进行拟合
考虑波长为400nm~700nm的可见光，可以将区间分为N段，以30段为例，通常足以拟合比较复杂的光谱数据
1 2 3 4 5 6 7 8 9 10 11 12 13 <Spectrum Utility Declarations> = static const int sampledLambdaStart = 370 ; static const int sampledLambdaEnd = 700 ; static const int nSpectralSamples = 30 ; <Spectrum Declarations> = class SampledSpectrum :public CoefficientSpectrum<nSpectralSamples> { public : <SampledSpectrum public Methods> private : <SampledSpectrum Private Data> };
那么如何从一组光谱数据来构造光谱表示?
假设现在有一组光谱数据$(\lambda_i,v_i)$，表示第$i$个采样点的波长为$\lambda_i$，值为$v_i$,这些采样点往往分布不均匀
对于给定的一组光谱数据，可以构造出相应的折线图
然后将此折线图划分为30个区域，如上图：[400,410) , [410,420)….[690,700)
计算每个区域覆盖的面积和，求其平均值，得到新的采样值
1 2 3 4 5 6 7 8 9 10 static SampledSpectrum FromSampled (const float * lambda , const float * v , int n) { <Sort samples if unordered, use sorted for returned spectrum> SampledSpectrum r; for (int i = 0 ; i < nSpectralSamples; i++) { <Compute average value of given SPD over ith sample' s range> } }
1 2 3 4 5 <Compute average value of given SPD over ith sample' s range>= float lambda0 = Lerp(float (i) / float (nSpectralSamples) , sampledLambdaStart , sampledLambdaEnd); float lambda1 = Lerp(float (i+1 ) / float (nSpectralSamples) , sampledLambdaStart , sampledLambdaEnd); r.c[i] = AverageSpectrumSamples(lambda, v, n, lambda0, lambda1);
1 2 3 4 5 6 7 8 9 10 11 <Spectrum Method Definitions> = float AverageSpectrumSamples (const float * lambda , const float * vals , int n , float lambdaStart , float lambdaEnd) { <handles case with out-of bounds range or single sample only> float sum = 0.0f ; <Add contribute of constant segments before/after samples> <Advance to first relevant wavelength , segment> <Loop over wavelength sample segments and add contributes> return sum / (lambdaEnd - lanbdaStart); }
处理边界条件：数据的波长范围全部处于区间之外，或者原始数据只有一个采样点
1 2 3 4 <handles case with out-of bounds range or single sample only> = if (lambdaEnd <= lambda[0 ]) return lambda[0 ]; if (lambdaStart >= lambda[n - 1 ]) return lambda[n - 1 ]; if (n == 1 ) return lambda[0 ];
对于处于数据之外的区域，将其值视为常量,即区域为方形
1 2 3 4 5 <Add contribute of constant segments before/after samples> = if (lambdaStart < lambda[0 ]) sum += vals[0 ] * (lambda[0 ] - lambdaStart); if (lambdaEnd > lambda[n - 1 ]) sum += vals[n - 1 ] * (lambdaEnd - lambda[n - 1 ]);
找到区间起始波长所在的数据段
1 2 3 <Advance to first relevant wavelength , segment> = int i = 0 ; while (lambdaStart > lambeda[i+1 ]) ++i;
计算区域覆盖的梯形面积和
1 2 3 4 5 6 7 8 9 10 11 12 <Loop over wavelength sample segments and add contributes> #define INTERP(w,i)\ Lerp(((w) - lambda[(i)]) / (lambda[(i) + 1 ] - lambda[(i)]) , vals[(i)] , vals[(i) + 1 ]) #define SEG_AVG(wl0,wl1,i) (0.5f * (INTERP(wl0,i) + INTERP(wl1,i))) for (; i+1 < n && lambdaEnd >= lambda[i]; i++) { float segStart = max(lambdaStart , lambda[i]); float segEnd = min(lambdaEnd , lambda[i+1 ]); sum += SEG_AVG(segStart,segEnd,i) * (segEnd - segStart); } #undef INTERP #undef SEG_AVG
经过以上步骤，将光谱数据拟合为线性分段函数，然而，最终我们需要将光谱转化为RGB才能在显示器上显示出来
三色刺激理论认为所有人眼可视的SPD都可以精确的使用三个值$(x_\lambda,y_\lambda,z_\lambda)$表示,公式如下：
$$x_\lambda = \frac{1}{\int Y(\lambda)d\lambda}\int_\lambda S(\lambda)X(\lambda)d\lambda$$
$$y_\lambda = \frac{1}{\int Y(\lambda)d\lambda}\int_\lambda S(\lambda)Y(\lambda)d\lambda$$
$$z_\lambda = \frac{1}{\int Y(\lambda)d\lambda}\int_\lambda S(\lambda)Z(\lambda)d\lambda$$
$X(\lambda),Y(\lambda),Z(\lambda)$为光谱匹配曲线,是由国际照明委员会经过一系列实验制定的.
这三条曲线被认为与人类视网膜的三种感色视椎体细胞响应类似
两个由多种不同波长的光混合而成的光源可以表现出同样的颜色，这叫做“异谱同色”(metamerism)
尽管对于人类观察者而言，XYZ可以很好的表示SPD，然而对于光谱计算并不是一种好的方案
XYZ的乘积与SPD的乘积再计算得到的XYZ相差很大
在SampledSpectrum中，SPD将会被变换重采样进行降维，因此，光谱匹配曲线也需要在重采样到相同维度才能进行运算
1 2 3 4 5 6 <Spectral Data Declarations> = static const int nCIESamples = 471 ; extern const float CIE_X[nCIESamples]; extern const float CIE_Y[nCIESamples]; extern const float CIE_Z[nCIESamples]; extern const float CIE_lambda[nCIESamples];
CIE_lambda中记录的是每个采样点的波长，这里使用1nm为步长，CIE_X,Y,Z记录从360nm~830nm的光谱匹配曲线采样点
1 2 3 <SampledSpectrum Private Data> = static SampledSpectrum X , Y , Z; static float yInt;
由于光谱匹配曲线是固定的，因此使用static变量
XYZ匹配曲线计算放在初始化函数中，当程序启动时进行调用
1 2 3 4 5 6 <SampledSpectrum Public Methods> += static void Init () { <Compute XYZ matching functions for SampledSpectrum> <Compute RGB to spectrum functions for SampledSpectrum> }
1 2 3 4 5 6 7 8 9 10 <Compute XYZ matching functions for SampledSpectrum> = for (int i = 0 ; i < nSpectralSamples; i++) { float wl0 = Lerp(float (i) / float (nSpectrualSamples) , sampledLambdaStart , sampledLambdaEnd); float wl1 = Lerp(float (i+1 ) / float (nSpectrualSamples) , sampledLambdaStart , sampledLambdaEnd); X.c[i] = AverageSpectrumSamples(CIE_lambda, CIE_X, nCIESamples, wl0, wl1); Y.c[i] = AverageSpectrumSamples(CIE_lambda, CIE_Y, nCIESamples, wl0, wl1); Z.c[i] = AverageSpectrumSamples(CIE_lambda, CIE_Z, nCIESamples, wl0, wl1); yInt += Y.c[i]; }
现在我们有了XYZ匹配曲线的重采样数据，公式可以转化为用黎曼和表示:
$$x_\lambda \approx \frac{1} {\int Y(\lambda)d\lambda} \sum_i X_i c_i$$
$$y_\lambda \approx \frac{1} {\int Y(\lambda)d\lambda} \sum_i Y_i c_i$$
$$z_\lambda \approx \frac{1} {\int Y(\lambda)d\lambda} \sum_i Z_i c_i$$
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 <SampledSpectrum Public Methods> += void ToXYZ (float xyz[3 ]) const { xyz[0 ] = xyz[1 ] = xyz[2 ] = 0.0f ; for (int i = 0 ; i < nSpectralSamples; i++) { xyz[0 ] += X.c[i] * c[i]; xyz[1 ] += Y.c[i] * c[i]; xyz[2 ] += Z.c[i] * c[i]; } xyz[9 ] /= yInt; xyz[1 ] /= yInt; xyz[2 ] /= yInt; }
XYZ中的y与亮度密切相关.
1 2 3 4 5 6 7 8 9 10 <SampledSpectrum public Methods> += float y () const { float yy = 0.0f ; for (int i = 0 ; i < nSpectralSamples; i++) { yy += Y.c[i] * c[i]; } return yy / yInt; }
现在我们已经将SPD转换为了XYZ,那么如何将XYZ转化为RGB用于显示呢？
当我们将RGB显示在显示器上时，真实显示的光谱主要由三色光谱响应曲线的加权和确定.，由发光材料决定
LED,LCD,以及等离子体的三色光谱响应曲线不同，因此相同的SPD,在这三种情况下显示出来的RGB是不同的.
在LCD,LED显示RGB(0.6,0.3,0.2)的光谱分别如下图所示：
可见SPD相差很大.因此，RGB作为颜色只在知道显示器属性的情况下才是有意义的.
给定一个SPD的$x_{\lambda},y_\lambda,z_\lambda$
在给定显示器的情况下,也就是给定三色光谱响应曲线$$R(\lambda),G(\lambda),B(\lambda)$$
RGB按照如下公式计算所得：
$$
$$
对于高清电视而言：
计算得到的RGB是在linear space,需要Gamma矫正才能变换到sRGB空间
1 2 3 4 5 6 7 <Spectrum Utility Declarations> += inline void XYZToRGB (const float xyz[3 ] , float rgb[3 ]) { rgb[0 ] = 3.240479 * xyz[0 ] - 1.537150 * xyz[1 ] - 0.498535 * xyz[2 ]; rgb[1 ] =-0.969256 * xyz[0 ] + 1.875991 * xyz[1 ] + 0.041556 * xyz[2 ]; rgb[2 ] = 0.055648 * xyz[0 ] - 0.204043 * xyz[1 ] + 1.057311 * xyz[2 ]; }
逆矩阵可以将RGB转为XYZ
1 2 3 4 5 6 7 <Spectrum Utility Declarations> += inline void RGBToXYZ (const float rgb[3 ], float xyz[3 ]) { xyz[0 ] = 0.412453 * rgb[0 ] + 0.357580 * rgb[1 ] + 0.180423 * rgb[2 ]; xyz[1 ] = 0.212671 * rgb[0 ] + 0.715160 * rgb[1 ] + 0.072169 * rgb[2 ]; xyz[2 ] = 0.019334 * rgb[0 ] + 0.119193 * rgb[1 ] + 0.950227 * rgb[2 ]; }
因此将SPD转为RGB:
1 2 3 4 5 6 7 <SampledSpectrum Public Methods> += void ToRGB (flaot rgb[3 ]) const { float xyz[3 ]; ToXYZ(xyz); XYZToRGB(xyz,rgb); }
介绍完了如何使用SPD计算RGB，那么如何通过RGB得到SPD?
通过测量得到以下数据
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 <Spectrual Data Declaration> += static const int nRGB2SpectSamples = 32 ; extern const float RGB2SpectLambda[nRGB2SpectSamples]; extern const float RGBRefl2SpectWhite[nRGB2SpectSamples]; extern const float RGBRefl2SpectCyan[nRGB2SpectSamples]; extern const float RGBRefl2SpectMagenta[nRGB2SpectSamples]; extern const float RGBRefl2SpectYellow[nRGB2SpectSamples]; extern const float RGBRefl2SpectRed[nRGB2SpectSamples]; extern const float RGBRefl2SpecGreen[nRGB2SpectSamples]; extern const float RGBRefl2SpectBlue[nRGB2SpectSamples]; extern const float RGBIllum2SpectWhite[nRGB2SpectSamples]; extern const float RGBIllum2SpectCyan[nRGB2SpectSamples]; extern const float RGBIllum2SpectMagenta[nRGB2SpectSamples]; extern const float RGBIllum2SpectYellow[nRGB2SpectSamples]; extern const float RGBIllum2SpectRed[nRGB2SpectSamples]; extern const float RGBIllum2SpectGreen[nRGB2SpectSamples]; extern const float RGBIllum2SpectBlue[nRGB2SpectSamples];
如果RGB颜色用于描述光源的发射光，那么应该使用表示发射源的光谱能量分布的转换表进行计算.
RGBillum2Spect*使用的是D65光谱能量分布，为CIE用于表示正午日光的标准.
用于计算方便，可以将RGBRefl2Spect*,RGBillum2Spect*分布，在初始化阶段重采样为SampledSpectrum
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 <SampledSpectrum Private Data> static SampledSpectrum rgbRelf2SpectWhite; static SampledSpectrum rgbRelf2SpectCyan; static SampledSpectrum rgbRelf2SpectMagenta; static SampledSpectrum rgbRelf2SpectYellow; static SampledSpectrum rgbRelf2SpectRed; static SampledSpectrum rgbRelf2SpectGreen; static SampledSpectrum rgbRelf2SpectBlue; static SampledSpectrum rgbIllum2SpectWhite; static SampledSpectrum rgbIllum2SpectCyan; static SampledSpectrum rgbIllum2SpectMagenta; static SampledSpectrum rgbIllum2SpectYellow; static SampledSpectrum rgbIllum2SpectRed; static SampledSpectrum rgbIllum2SpectGreen; static SampledSpectrum rgbIllum2SpectBlue;
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 <Compute RGB to spectrum functions for SampledSpectrum> = float wl0 = Lerp(float (i) / float (nSpectralSamples),sampledLambdaStart, sampledLambdaEnd); float wl1 = Lerp(float (i+1 ) / float (nSpectralSamples),sampledLambdaStart, sampledLambdaEnd); rgbRefl2SpectWhite.c[i] = AverageSpectrumSamples(RGB2SpectLambda, RGBRefl2SpectWhite, nRGB2SpectSamples, wl0, wl1); rgbRefl2SpectCyan.c[i] = AverageSpectrumSamples(RGB2SpectLambda, RGBRefl2SpectCyan, nRGB2SpectSamples, wl0, wl1); rgbRefl2SpectMagenta.c[i] = AverageSpectrumSamples(RGB2SpectLambda, RGBRefl2SpectMagenta, nRGB2SpectSamples, wl0, wl1); rgbRefl2SpectYellow.c[i] = AverageSpectrumSamples(RGB2SpectLambda, RGBRefl2SpectYellow, nRGB2SpectSamples, wl0, wl1); rgbRefl2SpectRed.c[i] = AverageSpectrumSamples(RGB2SpectLambda, RGBRefl2SpectRed, nRGB2SpectSamples, wl0, wl1); rgbRefl2SpectGreen.c[i] = AverageSpectrumSamples(RGB2SpectLambda, RGBRefl2SpectGreen, nRGB2SpectSamples, wl0, wl1); rgbRefl2SpectBlue.c[i] = AverageSpectrumSamples(RGB2SpectLambda, RGBRefl2SpectBlue, nRGB2SpectSamples, wl0, wl1); rgbIllum2SpectWhite.c[i] = AverageSpectrumSamples(RGB2SpectLambda, RGBIllum2SpectWhite, nRGB2SpectSamples, wl0, wl1); rgbIllum2SpectCyan.c[i] = AverageSpectrumSamples(RGB2SpectLambda, RGBIllum2SpectCyan, nRGB2SpectSamples, wl0, wl1); rgbIllum2SpectMagenta.c[i] = AverageSpectrumSamples(RGB2SpectLambda, RGBIllum2SpectMagenta, nRGB2SpectSamples, wl0, wl1); rgbIllum2SpectYellow.c[i] = AverageSpectrumSamples(RGB2SpectLambda, RGBIllum2SpectYellow, nRGB2SpectSamples, wl0, wl1); rgbIllum2SpectRed.c[i] = AverageSpectrumSamples(RGB2SpectLambda, RGBIllum2SpectRed, nRGB2SpectSamples, wl0, wl1); rgbIllum2SpectGreen.c[i] = AverageSpectrumSamples(RGB2SpectLambda, RGBIllum2SpectGreen, nRGB2SpectSamples, wl0, wl1); rgbIllum2SpectBlue.c[i] = AverageSpectrumSamples(RGB2SpectLambda, RGBIllum2SpectBlue, nRGB2SpectSamples, wl0, wl1);
RGB转Spectrum
$$if \quad r < g < b$$
$$(r,g,b) = (r,r,r)+(0,g-r,g-r)+(0,0,b-g)=r(1,1,1)+(g-r)(0,1,1)+(b-g)(0,0,1)$$
$$if \quad r < b < g$$
$$(r,g,b)=(r,r,r)+(0,b-r,b-r)+(0,g-b,0)=r(1,1,1)+(b-r)(0,1,1)+(g-b)(0,1,0)$$
…
(1,1,1):White
(0,1,1):Cyan
(0,0,1):Blue
(0,1,0):Green
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 <Spectrum Utility Declaration> += enum SpectrumType {SPECTRUM_REFLECTANCE, SPECTRUM_ILLUMINANT}; <Spectrum Method Definitions> += SampledSpectrum SampledSpectrum::FromRGB (const float rgb[3 ] , SpectrumType type) { SampledSpectrum r; if (type == SPECTRUM_REFLECTANCE) { <Convert reflectance spectrum to RGB> } else { <Convert illuminant spectrum to RGB> } return r.Clamp(); } <Convert reflectance spectrum to RGB> = if (rgb[0 ] <= rgb[1 ] && rgb[0 ] <= rgb[2 ]) { <Compute reflectance SampledSpectrum with rgb[0 ] as minimum> } else if (rgb[1 ] <= rgb[0 ] && rgb[1 ] <= rgb[2 ]) { <Compute reflectance SampledSpectrum with rgb[1 ] as minimum> } else { <Compute reflectance SampledSpectrum with rgb[2 ] as minimum> } <Compute reflectance SampledSpectrum with rgb[0 ] as minimum> = r += rgb[0 ] * rgbRefl2SpectWhite; if (rgb[1 ] <= rgb[2 ]) { r += (rgb[1 ] - rgb[0 ]) * rgbRefl2SpectCyan; r += (rgb[2 ] - rgb[1 ]) * rgbRefl2SpectBlue; } else { r += (rgb[2 ] - rgb[0 ]) * rgbRefl2SpectCyan; r += (rgb[1 ] - rgb[2 ]) * rgbRefl2SpectGreen; }
RGBSpectrum RGBSpectrum使用R,G,B分量的加权和来表示SPDs
1 2 3 4 5 6 <Spectrum Declarations> += class RGBSpectrum :public CofficientSpectrum<3 > { public : <RGBSpectrum Public Methods> };
FromRGB这个函数就简单了，因此我们存储的就是rgb
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 <RGBSpectrum Public Methods> += static RGBSpectrum FromRGB (const float rgb[3 ], SpectrumType type = SPECTRUM_REFLECTION) { RGBSpectrum s; s.c[0 ] = rgb[9 ]; s.c[1 ] = rgb[1 ]; s.c[2 ] = rgb[2 ]; return s; } <RGBSpectrum Public Methods> += void TOrgb (float * rgb) const { rgb[0 ] = c[0 ]; rgb[1 ] = c[1 ]; rgb[2 ] = c[2 ]; }
根据SPD计算RGB,通过光谱匹配函数计算XYZ,然后XYZ转为RGB
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 <RGBSpectrum Public Methods> += static RGBSpectrum FromSampled (const float * lambda, const float * v, int n) { float xyz[3 ] = {0 ,0 ,0 }; float yint = 0.0f ; for (int i = 0 ; i < nCIESamples; i++) { yint += CIE_Y[i]; float val = InterpolateSpectrumSamples(lambda, v, n, CIE_lambda[i]); xyz[0 ] += val * CIE_X[i]; xyz[1 ] += val * CIE_Y[i]; xyz[2 ] += val * CIE_Z[i]; } xyz[0 ] /= yint; xyz[1 ] /= yint; xyz[2 ] /= yint; return FromXYZ(xyz); } <Spectrum Method Definitions> += float InterpolateSpectrumSamples (const float * lambda, const float * vals, int n, float l) { if (l <= lambda[0 ]) return vals[0 ]; if (l >= lambda[n-1 ]) return vals[n-1 ]; for (int i=0 ; i < n-1 ; i++) { if (l >= lambda[i] && l <= lambda[i+1 ]) { float t = (l - lambda[i]) / (lambda[i+1 ] - lambda[i]); return Lerp(t , vals[i] , vals[i+1 ]); } } }
LinearSpace RGB to sRGB $$
Reference
评论