% KisoHB_DKLDemo_CompSpace
%
% 04/08/15	ek		wrote it.

% Clear
clear; close all

%% 同じ刺激の座標を、異なる色空間で表示して比較する

%% 錐体分光感度の設定

% calibrationデータの読み込み（分光測色のデータ例）
cal = LoadCalFile('PTB3TestCal');
S = cal.S_device;

whichCones = 'SmithPokorny';
% whichCones = 'StockmanSharpe';

[S_cones, T_cones, S_Y, T_Y] = SetFundamentals2(whichCones,S);

% 各錐体分光感度はピークが1.0となっている。
% L錐体分光感度とM錐体分光感度の和が分光視感効率に一致するようにスケーリング
factors = (T_cones(1:2,:)'\T_Y');
T_cones = diag([factors ; 1])*T_cones;

%% 計算に必要な初期設定

% 座標計算に必要な設定を行う
T_cones = SplineCmf(S_cones,T_cones,S);
cal_cones = SetSensorColorSpace(cal,T_cones,S);
cal_cones = SetGammaMethod(cal_cones,1);
cal_Lum = SetSensorColorSpace(cal,T_Y,S_Y);

% 背景をRGB primary = (0.5,0.5,0.5)として、錐体刺激値を計算する
bg_Primary = [0.5 0.5 0.5]';
bg_cones = PrimaryToSensor(cal_cones,bg_Primary);

% DKL座標計算用の変換行列を求める
% ディスプレイで表示できる最大値に基づきスケーリングしたもの
M_ConeIncToDKL_Device = OurComputeDKL_M_Device(bg_cones,T_cones,T_Y,cal_cones);

% MB色度図計算用の準備
% 等エネルギー白色に対してsMB=1.0とする
s_scaler = 1 * 683;
EEW_spd = ones(S_cones(3),1);
EEW_cones = T_cones*EEW_spd;
EEW_ls = OurLMSToMacBoyn2(EEW_cones,T_cones,T_Y,s_scaler);
s_scaler = 683/EEW_ls(2);
bg_ls = OurLMSToMacBoyn2(bg_cones,T_cones,T_Y,s_scaler);

% xy色度座標計算用にデータを読み込む
load T_xyz1931;
T_xyz = SplineCmf(S_xyz1931,683*T_xyz1931,S);
cal_XYZ = SetSensorColorSpace(cal,T_xyz,S);
bg_xyY = XYZToxyY(PrimaryToSensor(cal_XYZ,bg_Primary));

% load T_xyzJuddVos;
% T_xyzJuddVos = SplineCmf(S_xyzJuddVos,683*T_xyzJuddVos,S);
% cal_XYZJuddVos = SetSensorColorSpace(cal,T_xyzJuddVos,S);

%% xy色度座標が等しく、輝度の異なる刺激をプロットしてみる

% 計算の詳細はそれぞれのデモを参照のこと。

% 適当に値を設定する
test1_xyY = [0.35, 0.35, 85.0]';
test2_xyY = [0.35, 0.35, 70.0]';

test1_Primary = SensorToPrimary(cal_XYZ, xyYToXYZ(test1_xyY));
test2_Primary = SensorToPrimary(cal_XYZ, xyYToXYZ(test2_xyY));

test1_cones = PrimaryToSensor(cal_cones, test1_Primary);
test2_cones = PrimaryToSensor(cal_cones, test2_Primary);

% MB coordinates
test1_ls = OurLMSToMacBoyn2(test1_cones,T_cones,T_Y,s_scaler);
test2_ls = OurLMSToMacBoyn2(test2_cones,T_cones,T_Y,s_scaler);

% DKL coordinates
test1_coneInc = test1_cones - bg_cones;
test2_coneInc = test2_cones - bg_cones;

test1_dkl = M_ConeIncToDKL_Device*test1_coneInc;
test2_dkl = M_ConeIncToDKL_Device*test2_coneInc;

% cone contrast
test1_coneContrast = test1_coneInc ./ bg_cones;
test2_coneContrast = test2_coneInc ./ bg_cones;

% 変換結果をプロットする
% xy 色度図
ourMarkerSize = 8;
figure;
subplot(2,2,1);
plot(bg_xyY(1), bg_xyY(2),'kx', 'MarkerSize',ourMarkerSize);
hold on;
plot([bg_xyY(1), bg_xyY(1)], [0,1],'k:');
plot([0, 1],[bg_xyY(2), bg_xyY(2)],'k:');
plot(test1_xyY(1), test1_xyY(2), 'bo', 'MarkerSize',ourMarkerSize);
plot(test2_xyY(1), test2_xyY(2), 'r*', 'MarkerSize',ourMarkerSize);
axis equal;
axis([0 0.6 0 0.6]);	% 見やすいようにスケーリング。適当です。
title('xy chromaticity diagram', 'FontSize',16);
xlabel('x', 'FontSize', 14);
ylabel('y', 'FontSize', 14);

% MB色度図
subplot(2,2,2);
plot(bg_ls(1), bg_ls(2),'kx', 'MarkerSize',ourMarkerSize);
hold on;
plot([bg_ls(1), bg_ls(1)], [0,100],'k:');
plot([0, 1],[bg_ls(2), bg_ls(2)],'k:');
plot(test1_ls(1), test1_ls(2), 'bo', 'MarkerSize',ourMarkerSize);
plot(test2_ls(1), test2_ls(2), 'r*', 'MarkerSize',ourMarkerSize);
axis([0.55 0.75 0 4]);	% 見やすいようにスケーリング。適当です。
axis square;
title('MacLeod-Boynton diagram', 'FontSize',16);
xlabel('l_M_B', 'FontSize', 14);
ylabel('s_M_B', 'FontSize', 14);

% DKL
subplot(2,2,3)
plot(0,0, 'kx', 'MarkerSize',ourMarkerSize);
hold on;
plot([0, 0], [-10,10],'k:');
plot([-10, 10],[0, 0],'k:');
plot(test1_dkl(1), test1_dkl(2), 'bo', 'MarkerSize',ourMarkerSize);
plot(test2_dkl(1), test2_dkl(2), 'r*', 'MarkerSize',ourMarkerSize);
axis([-1 1 -1 1]);
axis square;
title('DKL equiluminant plane', 'FontSize',16);
xlabel('L-M', 'FontSize', 14);
ylabel('S-(L+M)', 'FontSize', 14);

% cone contrast
subplot(2,2,4)
plot(0,0, 'kx', 'MarkerSize',ourMarkerSize);
hold on;
plot([0, 0], [-1,1],'k:');
plot([-1, 1],[0, 0],'k:');
plot(test1_coneContrast(1), test1_coneContrast(2), 'bo', 'MarkerSize',ourMarkerSize);
plot(test2_coneContrast(1), test2_coneContrast(2), 'r*', 'MarkerSize',ourMarkerSize);
axis([-0.2 0.2 -0.2 0.2]);
axis square;
title('Cone Contrast (\DeltaL/L vs. \DeltaM/M)', 'FontSize',16);
xlabel('\DeltaL/L', 'FontSize', 14);
ylabel('\DeltaM/M', 'FontSize', 14);

