Contents

% KisoHB_Check_PTB_cones_sp.m
%
% PTBのT_cones_spやCVRLのSP錐体感度の検討。
%

clear all; close all;

% PTBのSmith & Pokorny fundamentals
load T_cones_sp;
wls = SToWls(S_cones_sp);		% 380 - 780 nm [380 5 81]

% CVRLからダウンロードした Smith & Pokorny fundamentals
buf = csvread('sp.csv');		% 380 - 825 nm
wls_sp = buf(:,1)';
CVRL_sp = 10.^buf(:,2:4);	% 対数値なので変換
CVRL_sp = CVRL_sp ./ repmat(max(CVRL_sp), size(CVRL_sp,1), 1);	% max=1にnormalize

% DeMarco, Pokorny & Smith (1992) versions of the Smith-Pokorny fundamentals
buf = csvread('dpse.csv');
wls_dpse = buf(:,1)';
CVRL_dpse = 10.^buf(:,2:4);
CVRL_dpse = CVRL_dpse ./ repmat(max(CVRL_dpse), size(CVRL_dpse,1), 1);

CVRL_sp = CVRL_sp';
CVRL_dpse = CVRL_dpse';

PTBのT_cones_spの確認：Judd-Vos修正等色関数から計算して、一致度を確認する。

% PTBのT_xyzJuddVosが、CVRLのものと一致することは、確認済み。
load T_xyzJuddVos

% 下記はSP fundamentals用の変換行列としてSmith & Pokorny (1975)に示されているものだが、
% Judd修正等色関数用の変換行列であることに注意。
% m_Judd =	[0.15514, 0.54312, -0.03286; ...
%			-0.15514, 0.45684,  0.03286; ...
%			 0,       0,        0.00801];

% Judd-Voss修正等色関数用の変換行列は下記。
% Smith, V.C., and Pokorny, J. (1996)
% The design and use of a cone chromaticity space: A tutorial.
% Color Research and Application 21, 375-383.
%
% 中野靖久先生の解説も参照のこと
% Nakano, Y. (1996) Color vision mathematics: A tutorial. pp. 544-562
% Kaiser & Boynton (1996) Human Color Vision, Appendix Part III
% 変換行列は、p. 557, A.3.15式
m =		[0.15516, 0.54308, -0.03287; ...
		-0.15516, 0.45692,  0.03287; ...
		 0,       0,        0.00801];

% SP fundamentalsの導出
LMS = m * T_xyzJuddVos;
LMS = LMS ./ repmat(max(LMS,[],2), 1, size(LMS,2));	% max=1にnormalize
wls_LMS = SToWls(S_xyzJuddVos);

計算結果とPTB SP fundamentalsの比較

figure(1);
plot(wls_LMS, LMS','k-'); hold on;
plot(wls, T_cones_sp(1,:)','rx');
plot(wls, T_cones_sp(2,:)','gx');
plot(wls, T_cones_sp(3,:)','bx');
xlim([350 750]);
xlabel('Wavelength (nm)', 'FontSize',14);
ylabel('Relative Sensitivity', 'FontSize',14);
title('PTB SP fundamentals', 'FontSize',14);

% 380-395 nm はゼロになっている（一番左端の方の青シンボルに注目）

計算結果とCVRL SP fundamentalsの比較

figure(2);
plot(wls_LMS, LMS','k-'); hold on;
plot(wls_sp, CVRL_sp(1,:)','rx');
plot(wls_sp, CVRL_sp(2,:)','gx');
plot(wls_sp, CVRL_sp(3,:)','bx');
xlim([350 750]);
xlabel('Wavelength (nm)', 'FontSize',14);
ylabel('Relative Sensitivity', 'FontSize',14);
title('CVRL SP fundamentals', 'FontSize',14);

計算結果とCVRL DPSE fundamentalsの比較

figure(3);
plot(wls_LMS, LMS','k-'); hold on;
plot(wls_dpse, CVRL_dpse(1,:)','rx');
plot(wls_dpse, CVRL_dpse(2,:)','gx');
plot(wls_dpse, CVRL_dpse(3,:)','bx');
xlim([350 750]);
xlabel('Wavelength (nm)', 'FontSize',14);
ylabel('Relative Sensitivity', 'FontSize',14);
title('CVRL DPSE fundamentals', 'FontSize',14);

% 400-700 nm しか値がない

計算結果と他のデータの差を計算してみる

計算結果とCVRL SP fundamentals

diff_CVRL_sp = sum(sum(abs(LMS - CVRL_sp(:,1:81))));
fprintf('Judd-Vosに基づく計算結果とCVRL SP fundamentalsの差 = %f\n', diff_CVRL_sp);

Judd-Vosに基づく計算結果とCVRL SP fundamentalsの差 = 0.001593

% ほとんど差がない、と書こうと思ったら、それほどでもない。
% これに関する検討結果は、後ろのセクションに。

計算結果とPTB SP fundamentals

diff_T_cones_sp = sum(sum(abs(LMS - T_cones_sp)));
fprintf('Judd-Vosに基づく計算結果とPTB SP fundamentalsの差 = %f\n', diff_T_cones_sp);

Judd-Vosに基づく計算結果とPTB SP fundamentalsの差 = 0.194557

% 差が大きい

計算したL, M錐体分光感度の和と、Judd-Vos修正分光視感効率一致するか？確認してみる。

factors = (LMS(1:2,:)'\T_xyzJuddVos(2,:)');
LMS2 = diag([factors ; 1])*LMS;
% LMS2 = diag([0.6373 0.3924 1]')*LMS;	% PTB LMSToMacBoyn.m内に書かれている係数
V = LMS2(1,:) + LMS2(2,:);

diff_V = sum(sum(abs(V-T_xyzJuddVos(2,:))));
fprintf('計算したL, M錐体分光感度の和とJudd-Vos修正分光視感効率の差 = %e\n', diff_V);

計算したL, M錐体分光感度の和とJudd-Vos修正分光視感効率の差 = 5.041294e-15

% 差はほぼゼロ。まぁそうですよね。

計算結果とCVRL SP fundamentalsの比較：詳細その１

% 予想よりも差が大きかったので、どこで差が大きいか調べてみると、
diff_LMS = LMS - CVRL_sp(:,1:81);
figure;
plot(wls_LMS,diff_LMS(1,:),'r-'); hold on;
plot(wls_LMS,diff_LMS(2,:),'g-');
plot(wls_LMS,diff_LMS(3,:),'b-');
xlabel('Wavelength (nm)', 'FontSize',14);
ylabel('Difference', 'FontSize',14);
title('CVRL SP fundamentals vs. calculated LMS', 'FontSize',14);

[diff_values, index] = max(abs(((LMS - CVRL_sp(:,1:81)))),[],2);
cone_str = 'LMS';
for co = 1:3,
	fprintf('Max difference for %s cone: %e at %d nm\n', cone_str(co), diff_values(co), wls_LMS(index(co)));
end

Max difference for L cone: 3.027765e-05 at 610 nm
Max difference for M cone: 2.261073e-05 at 440 nm
Max difference for S cone: 1.249499e-04 at 490 nm

% S錐体の490nmで相対的に大きな差があることが判明。他とは１桁違う。
% S錐体の分光感度は、変換行列mを見ればわかる通り、等色関数のZと同じはず。
% ここにズレがあるということは、数値が一貫していないということ。

% S錐体でのズレはJudd－Vos修正等色関数データの問題かと思い、確認。
% T_xyzJuddVos(:,23)
% ans =
%    0.031731000000000
%    0.208020000000000
%    0.446690000000000
%
% Vos(1978)と照らし合わせ、問題ないことを確認。

計算結果とCVRL SP fundamentalsの比較：詳細その２

% ちなみに、CVRLのホームページにも、上記mではなくm_Juddの変換式が載っているので、
% 試しに計算してみる：SP fundamentalsの導出
m_Judd =	[0.15514, 0.54312, -0.03286; ...
			-0.15514, 0.45684,  0.03286; ...
			 0,       0,        0.00801];

LMS_Judd = m_Judd * T_xyzJuddVos;
LMS_Judd = LMS_Judd ./ repmat(max(LMS_Judd,[],2), 1, size(LMS_Judd,2));
diff_LMS_Judd_vs_CVRL_sp = sum(sum(abs(LMS_Judd - CVRL_sp(:,1:81))));
fprintf('Judd用係数に基づく計算結果とCVRL SP fundamentalsの差 = %f\n', diff_LMS_Judd_vs_CVRL_sp);

diff_LMS = LMS_Judd - CVRL_sp(:,1:81);
figure;
plot(wls_LMS,diff_LMS(1,:),'r-'); hold on;
plot(wls_LMS,diff_LMS(2,:),'g-');
plot(wls_LMS,diff_LMS(3,:),'b-');
xlabel('Wavelength (nm)', 'FontSize',14);
ylabel('Difference', 'FontSize',14);
title('CVRL SP fundamentals vs. LMS derived with m-Judd', 'FontSize',14);

[diff_values, index] = max(abs(((LMS_Judd - CVRL_sp(:,1:81)))),[],2);
for co = 1:3,
	fprintf('Max difference for %s cone: %e at %d nm\n', cone_str(co), diff_values(co), wls_LMS(index(co)));
end

Judd用係数に基づく計算結果とCVRL SP fundamentalsの差 = 0.000532
Max difference for L cone: 1.256124e-05 at 595 nm
Max difference for M cone: 1.988493e-05 at 490 nm
Max difference for S cone: 1.249499e-04 at 490 nm

% 今度は、M錐体でも490nmでズレが大きいことになった。
% ただし、全体的には、L、M錐体に関する差は小さくなる。
% ここで計算結果とWebのデータとの間の差が小さくなるということは、
% Webデータの導出にmではなくm_Juddを使っているということではないか？

Judd(1951)での計算とCVRL SP fundamentalsの比較

% Judd(1951)による修正等色関数が計算に使われたのか、確認。
% CVRLからダウンロードした Judd修正等色関数を読み込む
buf = csvread('ciexyzj.csv');
wls_Judd = buf(:,1)';
T_Judd = buf(:,2:4)';

% Judd修正等色関数とm_Juddを用いて、錐体分光感度の導出。
LMS_Judd2 = m_Judd * T_Judd;
LMS_Judd2 = LMS_Judd2 ./ repmat(max(LMS_Judd2,[],2), 1, size(LMS_Judd2,2));

% プロットしてみる
figure;
plot(wls_Judd(2:end), LMS_Judd2(:,2:end)','k-');
hold on;
plot(wls_Judd(2:end), CVRL_sp(1,1:2:79)','rx');
plot(wls_Judd(2:end), CVRL_sp(2,1:2:79)','gx');
plot(wls_Judd(2:end), CVRL_sp(3,1:2:79)','bx');
ylim([0 1]);
xlabel('Wavelength (nm)', 'FontSize',14);
ylabel('Relative Sensitivity', 'FontSize',14);
title('CVRL SP fundamentals vs. LMS derived with Judd(1951)', 'FontSize',14);

diff_LMS = LMS_Judd2(:,2:end) - CVRL_sp(:,1:2:79);
figure;
plot(wls_Judd(2:end),diff_LMS(1,:),'r-'); hold on;
plot(wls_Judd(2:end),diff_LMS(2,:),'g-');
plot(wls_Judd(2:end),diff_LMS(3,:),'b-');
xlabel('Wavelength (nm)', 'FontSize',14);
ylabel('Difference', 'FontSize',14);
title('CVRL SP fundamentals vs. LMS derived with Judd(1951)', 'FontSize',14);

[diff_values, index] = max(abs(((LMS_Judd2(:,2:end) - CVRL_sp(:,1:2:79)))),[],2);
for co = 1:3,
	fprintf('Max difference for %s cone: %e at %d nm\n', cone_str(co), diff_values(co), wls_Judd(1+index(co)));
end

Max difference for L cone: 2.552606e-03 at 570 nm
Max difference for M cone: 2.253989e-03 at 540 nm
Max difference for S cone: 7.836241e-02 at 410 nm

% S錐体の短波長側でズレが大きい。
% これで、等色関数としてはJudd-Vos修正等色関数を使ったことがわかる。
% (JuddとJudd-Vosの比較は、下の方の図を参照のこと。）

% 結局、ズレの原因は不明。元のSPデータが対数値であることと関係があるのかもしれない。

% 結局のところ、自分で計算した結果が一番良いのではないかと考え、保存し活用することに。
% 以下、保存用のスクリプト。

おまけ

% 計算した SP data を保存する
% OurT_cones_sp = LMS;
% OurS_cones_sp = S_cones_sp;
% save('OurT_cones_sp', 'OurS_cones_sp', 'OurT_cones_sp');

Golz & MacLeod (2003): Supplementary data

% Golz, J., and MacLeod, D.I. (2003). Colorimetry for CRT displays.
% Journal of the Optical Society of America A 20, 769-781.
% http://www.uni-kiel.de/psychologie/golz/publications/2003a/supplement.html
%
% この論文のsupplementary filesとしてデータがあることを発見したので、さらに検討。
%
% ホームページよりデータを入手し、ファイル名を変更後（先頭に'GM_'を追加）、
% カレントディレクトリーに置いてあるものとする。

% Golz & MacLeodによるJudd-Vos修正等色関数の読み込み。
% これが、5nmステップであれば、PTB & CVRLのものと一致することは確認済み。
Golz_JV = dlmread('GM_XYZ_Vos.dat');
wls_Golz_JV = Golz_JV(:,1)';
Golz_JV = Golz_JV(:,2:4);
Golz_JV = Golz_JV';
% disp(sum(sum(T_xyzJuddVos-Golz_JV(:,1:5:end))))  % これはゼロになる

% % 1nmステップで補間した結果は少し違うが、違いの合計は3e-06より小さいので無視
% T_xyzJuddVos1nm = SplineCmf(S_xyzJuddVos, T_xyzJuddVos, [380 1 401]);
% diff_JuddVos_vs_Golz_JuddVos = sum(sum(abs(Golz_JV - T_xyzJuddVos1nm)));
% fprintf('Golz&MacLeod JuddVosと1nmステップに補間したT_xyzJuddVosの差 = %e\n', diff_JuddVos_vs_Golz_JuddVos);
% [diff_values, index] = max(abs(((Golz_JV - T_xyzJuddVos1nm))),[],2);
% for co = 1:3,
% 	fprintf('Max difference for %s cone: %e at %d nm\n', cone_str(co), diff_values(co), wls_Golz_JV(index(co)));
% end

% Golz & MacLeodのJudd修正等色関数
Golz_J = dlmread('GM_XYZ_Judd.dat');
wls_Golz_J = Golz_J(:,1)';
Golz_J = Golz_J(:,2:4);
Golz_J = Golz_J';

Golz & MacLeod (2003): まずは、Juddによる修正等色関数とVosによるさらなる修正とを比較してみる

figure;
plot(wls_Golz_J, Golz_J','k-');
hold on;
plot(wls_Golz_JV, Golz_JV(1,:)','rx');
plot(wls_Golz_JV, Golz_JV(2,:)','gx');
plot(wls_Golz_JV, Golz_JV(3,:)','bx');
% ylim([0 1]);
xlabel('Wavelength (nm)', 'FontSize',14);
ylabel('Tristimulus value', 'FontSize',14);
title('Golz: Judd(1951) vs. Vos(1978)', 'FontSize',14);
legend('Judd-X','Judd-Y','Judd-Z','Vos-X','Vos-Y','Vos-Z');

% 400nm付近、あるいは、それより短波長側で差あり

Golz & MacLeod (2003): Juddに基づく計算とJudd-Vosに基づく計算を比較してみる

Golz_LMS_J = m_Judd * Golz_J;
Golz_LMS_J = Golz_LMS_J ./ repmat(max(Golz_LMS_J,[],2), 1, size(Golz_LMS_J,2));

Golz_LMS_JV = m * Golz_JV;
Golz_LMS_JV = Golz_LMS_JV ./ repmat(max(Golz_LMS_JV,[],2), 1, size(Golz_LMS_JV,2));

figure;
plot(wls_Golz_J, Golz_LMS_J','k-');
hold on;
plot(wls_Golz_JV, Golz_LMS_JV(1,:)','rx');
plot(wls_Golz_JV, Golz_LMS_JV(2,:)','gx');
plot(wls_Golz_JV, Golz_LMS_JV(3,:)','bx');
ylim([0 1]);
xlabel('Wavelength (nm)', 'FontSize',14);
ylabel('Relative Sensitivity', 'FontSize',14);
title('Golz: LMS from Judd(1951) vs. LMS from Vos(1978)', 'FontSize',14);

% 目に見える差があるのは、S錐体だけ

Golz & MacLeod (2003): SP錐体感度に関する確認

Golz_sp = dlmread('GM_LMS_SP.dat');
wls_Golz_sp = Golz_sp(:,1)';
Golz_sp = Golz_sp(:,2:4);
Golz_sp = Golz_sp ./ repmat(max(Golz_sp), size(Golz_sp,1), 1);
Golz_sp = Golz_sp';

figure;
plot(wls_Golz_J, Golz_LMS_JV','k-');
hold on;
plot(wls_Golz_sp, Golz_sp(1,:)','rx');
plot(wls_Golz_sp, Golz_sp(2,:)','gx');
plot(wls_Golz_sp, Golz_sp(3,:)','bx');
ylim([0 1]);
xlabel('Wavelength (nm)', 'FontSize',14);
ylabel('Relative Sensitivity', 'FontSize',14);
title('Golz: LMS from Judd-Vos vs. SP fundamentals', 'FontSize',14);

diff_Golz_LMS_JV_vs_Golz_sp = sum(sum(abs(Golz_LMS_JV - Golz_sp)));
fprintf('Golz&MacLeodのJudd-Vos修正等色関数に基づく計算結果とGolz_spの差 = %f\n', diff_Golz_LMS_JV_vs_Golz_sp);

diff_LMS = Golz_LMS_JV - Golz_sp;
figure;
plot(wls_Golz_sp,diff_LMS(1,:),'r-'); hold on;
plot(wls_Golz_sp,diff_LMS(2,:),'g-');
plot(wls_Golz_sp,diff_LMS(3,:),'b-');
xlabel('Wavelength (nm)', 'FontSize',14);
ylabel('Difference', 'FontSize',14);
title('Golz: LMS from Judd-Vos vs. SP fundamentals', 'FontSize',14);

Golz&MacLeodのJudd-Vos修正等色関数に基づく計算結果とGolz_spの差 = 0.018581

% うーむ、何か不思議な結果に。
% これは一致すると思ったのですが、S錐体でズレが相対的に大きい。
% 計算してから補間したとか、データの補間方法の問題でしょうか？
% ここでも、Judd-Vos修正等色関数のZとSP錐体分光感度のSの値が一貫していない。

Golz & MacLeod (2003): SP錐体感度と自分でJudd-Vosから計算した結果の比較

figure;
plot(wls_Golz_sp, Golz_sp','k-');
hold on;
plot(wls_LMS, LMS(1,:)','rx');
plot(wls_LMS, LMS(2,:)','gx');
plot(wls_LMS, LMS(3,:)','bx');
ylim([0 1]);
xlabel('Wavelength (nm)', 'FontSize',14);
ylabel('Relative Sensitivity', 'FontSize',14);
title('Golz SP fundamentals vs. Our LMS derived with Judd-Vos', 'FontSize',14);

diff_Golz_sp_vs_LMS = sum(sum(abs(Golz_sp(:,1:5:end) - LMS)));
fprintf('Golz&MacLeodのSP錐体感度とJudd-Vos修正等色関数に基づく自前の計算結果との差 = %f\n', diff_Golz_sp_vs_LMS);

diff_LMS = Golz_sp(:,1:5:end) - LMS;
figure;
plot(wls_LMS,diff_LMS(1,:),'r-'); hold on;
plot(wls_LMS,diff_LMS(2,:),'g-');
plot(wls_LMS,diff_LMS(3,:),'b-');
xlabel('Wavelength (nm)', 'FontSize',14);
ylabel('Difference', 'FontSize',14);
title('Golz SP fundamentals vs. Our LMS derived with Judd-Vos', 'FontSize',14);

Golz&MacLeodのSP錐体感度とJudd-Vos修正等色関数に基づく自前の計算結果との差 = 0.023057

% ここでも微妙な差が。

% 細かいところでいろいろとあるもので、
% 結局のところ、自分で計算した結果が一番良いのではないかという結論は変わらず。

Published with MATLAB® R2014b