%% KisoHB_ColorSpaceTransformation.m

% ある色空間から別の色空間への変換のデモ
% colorSpaceTutorial.m (06/18/98  dhb) に基づいている。

% Initialization
clear all; close all

%% Smith & Pokorny (1975) 錐体分光感度 と Judd-Vos修正等色関数を読み込む

% Smith & Pokorny (1975) 錐体分光感度を読み込む
load OurT_cones_sp
T_cones_sp = OurT_cones_sp;
S_cones_sp = OurS_cones_sp;

% Judd-Vos修正等色関数を読み込む
load T_xyzJuddVos

% 波長ステップの計算。
% PTBでは、波長データを[開始波長 ステップ 波長ステップ数]という独特の形式で表す。
% 例えば、S_cones_sp = [380, 5, 81]　これは、380nmから5nmステップで81波長という意味。
% このデータを波長ベクトルに変換する。
wls = SToWls(S_cones_sp);

%% 各データをプロットしてみる。

%% Fig. 1: Smith & Pokorny (1975) 錐体分光感度のプロット
figure(1);
plot(wls,T_cones_sp(1,:)','r');
hold on
plot(wls,T_cones_sp(2,:)','g');
plot(wls,T_cones_sp(3,:)','b');
xlabel('Wavelength (nm)', 'FontSize',14);
ylabel('Relative Sensitivity', 'FontSize',14);
title('Fig.1: Smith & Pokorny Cone Spectral Sensitivity Functions', 'FontSize',14);

%% Fig. 2: Judd-Vos修正等色関数のプロット
figure(2);
plot(wls,T_xyzJuddVos(1,:)','r');
hold on
plot(wls,T_xyzJuddVos(2,:)','g');
plot(wls,T_xyzJuddVos(3,:)','b');
xlabel('Wavelength (nm)', 'FontSize',14);
ylabel('Tristimulus value', 'FontSize',14);
title('Fig.2: Judd-Vos Modified Color Matching Functions', 'FontSize',14);

%% Judd-Vos修正等色関数からSmith & Pokorny (1975) 錐体分光感度への変換行列を求める
M_JuddVosToCones = ((T_xyzJuddVos')\(T_cones_sp'))';

% 求めた変換行列を使って、錐体分光感度を実際に計算してみる。
T_cones_fromJuddVos = M_JuddVosToCones*T_xyzJuddVos;

%% Fig. 3: 計算結果と元データを重ねて確認
figure(3);
plot(wls,T_cones_sp(1,:)','r');
hold on
plot(wls,T_cones_sp(2,:)','g');
plot(wls,T_cones_sp(3,:)','b');
xlabel('Wavelength (nm)', 'FontSize',14);
ylabel('Relative Sensitivity', 'FontSize',14);
title('Fig.3: Smith & Pokorny Cone Spectral Sensitivity Functions', 'FontSize',14);
plot(wls,T_cones_fromJuddVos(1,:)','rs');
plot(wls,T_cones_fromJuddVos(2,:)','gs');
plot(wls,T_cones_fromJuddVos(3,:)','bs');
axis([350 800 -0.05 1.1]);

%%

% 元々同じ標準観測者の錐体分光感度と等色関数なので、きれいに一致する。

%% Smith & Pokorny (1975) 錐体分光感度 と CIE 1931 等色関数

% 同じことを、Judd-Vosではなく、CIE1931でやってみる
% Smith & Pokorny (1975) 錐体分光感度を読み込む
load OurT_cones_sp

% CIE1931等色関数を読み込む
load T_xyz1931;

%% Fig4: CIE1931等色関数のプロット

% データをプロットしてみる。
figure(4);
plot(wls,T_xyz1931(1,:)','r');
hold on
plot(wls,T_xyz1931(2,:)','g');
plot(wls,T_xyz1931(3,:)','b');
xlabel('Wavelength (nm)', 'FontSize',14);
ylabel('Tristimulus value', 'FontSize',14);
title('Fig.4: CIE1931 Color Matching Functions', 'FontSize',14);

%% CIE1931等色関数からSmith & Pokorny (1975) 錐体分光感度への変換行列を無理矢理求める
M_1931ToCones = ((T_xyz1931')\(T_cones_sp'))';

% 求めた変換行列を使って、錐体分光感度を実際に計算してみる。
T_cones_from1931 = M_1931ToCones*T_xyz1931;

%% Fig.5: 計算結果を元データに重ねてみる
figure(5);
plot(wls,T_cones_sp(1,:)','r');
hold on
plot(wls,T_cones_sp(2,:)','g');
plot(wls,T_cones_sp(3,:)','b');
xlabel('Wavelength (nm)', 'FontSize',14);
ylabel('Relative Sensitivity', 'FontSize',14);
title('Fig.5: Smith&Pokorny vs. ConesFromCIE1931', 'FontSize',14);
plot(wls,T_cones_from1931(1,:)','ro');
plot(wls,T_cones_from1931(2,:)','go');
plot(wls,T_cones_from1931(3,:)','bo');
axis([350 800 -0.05 1.1]);

%%

% 異なる標準観測者のデータを無理矢理変換しているので、当然のことながら一致は良くない。
% 特にS錐体でのズレが大きいが、これは、Judd-Vosによる修正が主として短波長領域に関するものであるため。
%
% 色彩輝度計で測定されるxy色度座標と輝度値(Y値）は、CIE1931 XYZ等色関数を使用しているが、この等色関数を
% 線形変換してSmith & Pokornyの錐体分光感度を導出することはできない（三刺激値から錐体刺激値を正確には
% 計算できない）。
%
% 錐体刺激値を計算するには、分光放射輝度を計測し、それと錐体分光感度から計算する。
%

%% Stockman & Sharpe (2000) 錐体分光感度 と CIE 1931 等色関数

% 今度は同じことを、Stockman & Sharpe (2000)とCIE1931でやってみる
% CIE1931等色関数を読み込む
load T_xyz1931;

% Stockman & Sharpe (2000) 錐体分光感度を読み込む
load OurT_cones_ss2;
T_cones_ss2 = OurT_cones_ss2;
S_cones_ss2 = OurS_cones_ss2;

% 波長ステップがデータセット間で異なるので補間する
T_xyz1931 = SplineCmf(S_xyz1931,T_xyz1931,S_cones_ss2);
wls = SToWls(S_cones_ss2);

% CIE1931等色関数からStockman & Sharpe (2000)錐体分光感度への変換行列を無理矢理求める
M_1931ToCones_SS = ((T_xyz1931')\(T_cones_ss2'))';

% 求めた変換行列を使って、錐体分光感度を実際に計算してみる。
T_cones_SS2_from1931 = M_1931ToCones_SS*T_xyz1931;

%% Fig.6: 計算結果を元データに重ねてみる
figure(6);
ourStep = 1:5:S_cones_ss2(3);	% 見づらいので粗くプロットする
plot(wls,T_cones_ss2(1,:)','r');
hold on
plot(wls,T_cones_ss2(2,:)','g');
plot(wls,T_cones_ss2(3,:)','b');
xlabel('Wavelength (nm)', 'FontSize',14);
ylabel('Relative Sensitivity', 'FontSize',14);
title('Fig.6: Stockman&Sharpe vs. ConesFromCIE1931-SS', 'FontSize',14);
plot(wls(ourStep),T_cones_SS2_from1931(1,ourStep)','r*');
plot(wls(ourStep),T_cones_SS2_from1931(2,ourStep)','g*');
plot(wls(ourStep),T_cones_SS2_from1931(3,ourStep)','b*');
axis([350 800 -0.05 1.1]);

%%

% やはり、S錐体の分光感度でズレが大きい。

%% Stockman & Sharpe (2000) 錐体分光感度 と Judd-Vos修正等色関数

% 今度は同じことを、Stockman & Sharpe (2000)とJudd-Vosででやってみる
% Judd-Vos修正等色関数を読み込む
load T_xyzJuddVos

% Stockman & Sharpe (2000) 錐体分光感度を読み込む
load T_cones_ss2;

% 波長ステップがデータセット間で異なるので補間する
T_xyzJuddVos = SplineCmf(S_xyzJuddVos,T_xyzJuddVos,S_cones_ss2);
wls = SToWls(S_cones_ss2);

% Judd-Vos修正等色関数からStockman & Sharpe (2000) 錐体分光感度への変換行列を無理矢理求める
M_JuddVosToCones_SS = ((T_xyzJuddVos')\(T_cones_ss2'))';

% 求めた変換行列を使って、錐体分光感度を実際に計算してみる。
T_cones_SS_fromJuddVos = M_JuddVosToCones_SS*T_xyzJuddVos;

%% Fig.7: 計算結果を元データに重ねてみる
ourStep = 1:5:S_cones_ss2(3);
figure(7);
plot(wls,T_cones_ss2(1,:)','r');
hold on
plot(wls,T_cones_ss2(2,:)','g');
plot(wls,T_cones_ss2(3,:)','b');
xlabel('Wavelength (nm)', 'FontSize',14);
ylabel('Relative Sensitivity', 'FontSize',14);
title('Fig.7: Stockman&Sharpe vs. ConesFromJuddVos-SS', 'FontSize',14);
plot(wls(ourStep),T_cones_SS_fromJuddVos(1,ourStep)','rx');
plot(wls(ourStep),T_cones_SS_fromJuddVos(2,ourStep)','gx');
plot(wls(ourStep),T_cones_SS_fromJuddVos(3,ourStep)','bx');
axis([350 800 -0.05 1.1]);

%%

% やはり、S錐体の分光感度でズレが大きい。



