从视觉响应到色品图

       很久很久之前，人们对颜色既熟悉又陌生，自从牛顿发现了赤橙黄绿青蓝紫，这个世界才开始丰富多彩起来。人类具有3种类型的视觉感官细胞，不同波段有不一样的亮度感觉，如下：

wave_l=csvread('./cie lms function/linCIE2008v10e_fine.csv',0,0,[0,0,4400,0]);
lumin =csvread('./cie lms function/linCIE2008v10e_fine.csv',0,1,[0,1,4400,1]);
plot(wave_l,lumin,'k');
xlabel('wave length /nm');
ylabel('Efficiency');
title('Luminous efficiency');

     科学家们为了确定人眼对光的敏感程度，他们找来一群人做测试，去观察感觉不同波段的光的亮度，那么在亮度平均一致的情况下，可以得到光的辐射功率。因为光是在实验室里面通过机器发出的，功率数值是显而易见的。为了达到同样的视觉亮度，很明显功率越大的波，越不容易被人眼感觉到，于是人眼的光谱视觉效应可以定义为功率谱的倒数。经过实验数据知道，波长为555纳米时，感觉最明显，因此取为最大值，其他波长的光归一化到此点上，得到了下图的数据。

wave_l=csvread('./cie lms function/linss10e_fine.csv',0,0,[0,0,4400,0]);
cone_l=csvread('./cie lms function/linss10e_fine.csv',0,1,[0,1,4400,1]);
cone_m=csvread('./cie lms function/linss10e_fine.csv',0,2,[0,2,4400,2]);
cone_s=csvread('./cie lms function/linss10e_fine.csv',0,3,[0,3,4400,3]);
figure
plot(wave_l,cone_s,'b');hold on
plot(wave_l,cone_m,'g');hold on
plot(wave_l,cone_l,'r');
legend('short wave','Medium wave','long wave','Location','east');
xlabel('wave length /nm');
ylabel('sensitivity');
title('Cone Fundaments');
 

上面的图就是三种视锥细胞的感光响应曲线，我一下想不到科学家是怎么做实验得到这种数据的。为了量化光谱的颜色，他们选择了三原色去拟合其他光的颜色，波长分别为700nm,546.1nm,435.8nm。得出了下面的参数曲线，出现负值的原因是无论怎么调节三基色的比例都无法适配待定光谱颜色，只有在待配色中加入一种颜色才可以等价，相当于三基色的一个分量基为负。

因为出现了负值，在应用上面造成了很大的不便。于是科学家们做了一种线性变换，近视的从rgb转换到xyz空间里面。

rgb2xyz=[2.7689    1.7522    1.1296;
               1.0005    4.5916    0.0602;
               0.0001    0.0554    5.5943];
tri_xyz = rgb2xyz*[color_r color_g color_b]';
figure
plot(wave_c,tri_xyz(1,:),'r');hold on
plot(wave_c,tri_xyz(2,:),'g');hold on
plot(wave_c,tri_xyz(3,:),'b');
legend('x component','y component','z component','Location','east');
xlabel('wave length /nm');
ylabel('amplitude');
title('color match functions');

随着科学技术的进步，科学测量越来越发达，科学家直接导出了lms空间到xyz的映射矩阵。这样以来得出的数据就更加精确，下面给出10度观察者的参数，并画出色品图。

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wave_l=csvread('./cie lms function/linss10e_fine.csv',0,0,[0,0,4400,0]);
cone_l=csvread('./cie lms function/linss10e_fine.csv',0,1,[0,1,4400,1]);
cone_m=csvread('./cie lms function/linss10e_fine.csv',0,2,[0,2,4400,2]);
cone_s=csvread('./cie lms function/linss10e_fine.csv',0,3,[0,3,4400,3]);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
lms2xyz=[1.93986443     -1.34664359 0.43044935 ;%10度观察
                0.69283932     0.34967567  0 ;
                0                      0                   2.14687945 ] ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c_xyz =lms2xyz*[cone_l cone_m cone_s]';
color_x=c_xyz(1,:);
color_y=c_xyz(2,:);
color_z=c_xyz(3,:);
normal_x=color_x./(color_x + color_y + color_z);%归一化
normal_y=color_y./(color_x + color_y + color_z);
normal_z=color_z./(color_x + color_y + color_z);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
plot(normal_x,normal_y,'k');
xlabel('x coordinate');
ylabel('y coordinate');
title('chromaticity coordinates');

下面以sRGB为基准，将打印出光谱色彩，需要说明的是sRGB无法完全拟合上面的xyz空间。

sXYZ_R  = [0.64 0.33 0.03]';%R分量在色域中的坐标
sXYZ_G  = [0.30 0.60 0.10]';%G分量在色域中的坐标
sXYZ_B  = [0.15 0.06 0.79]';%B分量在色域中的坐标
sRGB_W  =[0.95047 1.00000 1.08883]';%白色混色比率

显然白色点的位置这样被标定：[sXYZ_R sXYZ_G sXYZ_B] * w=sRGB_W  
w=inv([sXYZ_R sXYZ_G sXYZ_B])*sRGB_W;

[sXYZ_R sXYZ_G sXYZ_B] =

   [0.6400    0.3000    0.1500
    0.3300    0.6000    0.0600
    0.0300    0.1000    0.7900]

w=[0.6445     1.1919    1.2029]'

于是任意颜色定义为C_xyz'=[sXYZ_R*w(1), sXYZ_G*w(2), sXYZ_B*w(3)] *C_rgb';

因此色彩空间rgb2xyz=[sXYZ_R*w(1), sXYZ_G*w(2), sXYZ_B*w(3)]；

rgb2xyz =0.4125    0.3576    0.1804
               0.2127    0.7152    0.0722
               0.0193    0.1192    0.9503

xyz2rgb=3.2405   -1.5371   -0.4985
               -0.9693    1.8760    0.0416
                0.0556   -0.2040    1.0572

将所得到的转换矩阵使用函数形式表达如下：

function [rgb] = csc_xyz2srgb(xyz)
    M = [ 3.2405   -1.5371   -0.4985;
         -0.9693    1.8760    0.0416;
         0.0556   -0.2040    1.0572];
     
    [rows cols ] = size(xyz);
    rgb = M*xyz;
    for c = 1:cols
        for ch = 1:rows
            %伽马校正
            if rgb(ch,c) <= 0.0031308
                rgb(ch,c) = 12.92*rgb(ch,c);
            else
                rgb(ch,c) = 1.055*(rgb(ch,c)^(1.0/2.4)) - 0.055;
            end
            %超边界截断
            if rgb(ch,c) < 0
                rgb(ch,c) = 0;
            elseif rgb(ch,c) > 1
                rgb(ch,c) = 1;
            end
        end
    end
end

代码如下:

%打印光谱的颜色
[rows cols] =size(normal_x);
face=[1 2 3 4];
vert=zeros(4,2,'double');

for i=1:cols
    rgb=csc_xyz2srgb([normal_x(i);normal_y(i);normal_z(i)]);
    vert(1,1)=(i-1)/cols*800;
    vert(1,2)=0;
    vert(2,1)=(i-1)/cols*800;
    vert(2,2)=800;    
    vert(4,1)=i/cols*800;
    vert(4,2)=0;
    vert(3,1)=i/cols*800;
    vert(3,2)=800;
   patch('Vertices',vert, 'Faces',face, 'EdgeColor','none', 'FaceVertexCData',[rgb,rgb,rgb,rgb]','FaceColor','interp');
end

代码如下：

load locus.mat
plot(locus(:,1),locus(:,2),'k','LineWidth',2)
grid on
hold on
axis([0.0 0.85 0.0 0.85])
xlabel('x')
ylabel('y')
% plot the non-spectral locus
plot([locus(1,1) locus(end,1)], [locus(1,2) locus(end,2)],'k','LineWidth',2)
% chromaticity coordinates of spectrum locus
x = [ 0.175596 0.172787 0.170806 0.170085 0.160343 0.146958 0.139149 ...
      0.133536 0.126688 0.115830 0.109616 0.099146 0.091310 0.078130 ...
      0.068717 0.054675 0.040763 0.027497 0.016270 0.008169 0.004876 ...
      0.003983 0.003859 0.004646 0.007988 0.013870 0.022244 0.027273 ...
      0.032820 0.038851 0.045327 0.052175 0.059323 0.066713 0.074299 ...
      0.089937 0.114155 0.138695 0.154714 0.192865 0.229607 0.265760 ...
      0.301588 0.337346 0.373083 0.408717 0.444043 0.478755 0.512467 ...
      0.544767 0.575132 0.602914 0.627018 0.648215 0.665746 0.680061 ...
      0.691487 0.700589 0.707901 0.714015 0.719017 0.723016 0.734674 ]';
y = [ 0.005295 0.004800 0.005472 0.005976 0.014496 0.026643 0.035211 ...
      0.042704 0.053441 0.073601 0.086866 0.112037 0.132737 0.170464 ...
     0.200773 0.254155 0.317049 0.387997 0.463035 0.538504 0.587196 ...
      0.610526 0.654897 0.675970 0.715407 0.750246 0.779682 0.792153 ...
      0.802971 0.812059 0.819430 0.825200 0.829460 0.832306 0.833833 ...
      0.833316 0.826231 0.814796 0.805884 0.781648 0.754347 0.724342 ...
      0.692326 0.658867 0.624470 0.589626 0.554734 0.520222 0.486611 ...
      0.454454 0.424252 0.396516 0.372510 0.351413 0.334028 0.319765 ...
      0.308359 0.299317 0.292044 0.285945 0.280951 0.276964 0.265326 ]';
N = length(x);
i = 1;
e = 1/3;
steps = 25;
xy4rgb = zeros(N*steps*4,5,'double');
for w = 1:N                                     % wavelength
    w2 = mod(w,N)+1;
    a1 = atan2(y(w)  -e,x(w)  -e);              % start angle
    a2 = atan2(y(w2) -e,x(w2) -e);              % end angle
    r1 = ((x(w)  - e)^2 + (y(w)  - e)^2)^0.5;   % start radius
    r2 = ((x(w2) - e)^2 + (y(w2) - e)^2)^0.5;   % end radius
    for c = 1:steps                               % colourfulness
        % patch polygon
        xyz(1,1) = e+r1*cos(a1)*c/steps;
        xyz(1,2) = e+r1*sin(a1)*c/steps;
        xyz(1,3) = 1 - xyz(1,1) - xyz(1,2);
        xyz(2,1) = e+r1*cos(a1)*(c-1)/steps;
        xyz(2,2) = e+r1*sin(a1)*(c-1)/steps;
        xyz(2,3) = 1 - xyz(2,1) - xyz(2,2);
        xyz(3,1) = e+r2*cos(a2)*(c-1)/steps;
        xyz(3,2) = e+r2*sin(a2)*(c-1)/steps;
        xyz(3,3) = 1 - xyz(3,1) - xyz(3,2);
        xyz(4,1) = e+r2*cos(a2)*c/steps;
        xyz(4,2) = e+r2*sin(a2)*c/steps;
        xyz(4,3) = 1 - xyz(4,1) - xyz(4,2);
        % compute sRGB for vertices
        rgb = csc_xyz2srgb(xyz');
        % store the results
        xy4rgb(i:i+3,1:2) = xyz(:,1:2);
        xy4rgb(i:i+3,3:5) = rgb';
        i = i + 4;
    end
end
[rows cols] = size(xy4rgb);
f = [1 2 3 4];
v = zeros(4,3,'double');
for i = 1:4:rows
    v(:,1:2) = xy4rgb(i:i+3,1:2);
    patch('Vertices',v, 'Faces',f, 'EdgeColor','none', ...
    'FaceVertexCData',xy4rgb(i:i+3,3:5),'FaceColor','interp')
end