% KisoHB_DKL_vs_ConeContrast.m
%
% DKL空間と錐体コントラスト空間の関係
%
% 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,0);
cal_Lum = SetSensorColorSpace(cal,T_Y,S_Y);

% 背景をRGB=(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);

%% DKL空間の L-M vs. L+M 平面で円を描き、それを錐体コントラスト空間にプロットしてみる

angle = [0:10:90, 80:-10:0, -10:-10:-90, -80:10:-10];
nAngle = length(angle);
radius = 0.7;
circle_DKL_polar = [radius*ones(1,nAngle); [zeros(1,10), 180*ones(1,18), zeros(1,8)]; angle];
circle_DKL = PolarDegToLinearDKL(circle_DKL_polar);

circle_coneInc = M_ConeIncToDKL_Device\circle_DKL;
circle_cones = circle_coneInc + repmat(bg_cones,1,nAngle);
circle_settings = SensorToSettings(cal_cones, circle_cones);

% bg_Primary_all=repmat(bg_Primary,1,nAngle);
% circlePrimaryInc = SensorToPrimary(cal_cones, circle_cones) - bg_Primary_all;
% circleScale = OurMaximizeGamutContrast(circlePrimaryInc, bg_Primary);
% maxRadius = min(circleScale);

%%  DKL空間の L-M vs. L+M 平面上の円
figure;
ourMarkerSize = 8;
plot([0, 0], [-10,10],'k:');
hold on;
plot([-10, 10],[0, 0],'k:');

for co = 1:nAngle,
	plot(circle_DKL(2,co), circle_DKL(1,co), 'ko', 'MarkerSize',ourMarkerSize, ...
		'MarkerEdgeColor', circle_settings(:,co), 'MarkerFaceColor', circle_settings(:,co));
end;
axis([-1 1 -1 1]);
axis square;
title('DKL L-M vs. L+M plane', 'FontSize',16);
xlabel('L-M', 'FontSize', 14);
ylabel('L+M', 'FontSize', 14);

index1 = 2;
c1 = polyfit([0,circle_DKL(2,index1)],[0,circle_DKL(1,index1)],1);
plot([0,1],[0,c1(1)],'-','Color',circle_settings(:,index1));
index2 = 18;
c2 = polyfit([0,circle_DKL(2,index2)],[0,circle_DKL(1,index2)],1);
plot([0,-1],[0,-c2(1)],'-','Color',circle_settings(:,index2));
index3 = 20;
c3 = polyfit([0,circle_DKL(2,index3)],[0,circle_DKL(1,index3)],1);
plot([0,-1],[0,-c3(1)],'-','Color',circle_settings(:,index3));
index4 = 36;
c4 = polyfit([0,circle_DKL(2,index4)],[0,circle_DKL(1,index4)],1);
plot([0,1],[0,c4(1)],'-','Color',circle_settings(:,index4));

fill([0, 1, 1, 0], [0, c1(1), c4(1), 0],[0.9,0.9,0.9]);
fill([0, -1, -1, 0], [0, -c2(1), -c3(1), 0],[0.9,0.9,0.9]);

% プロットを再度描いて、見栄えを良くする
plot([0, 0], [-10,10],'k:');
plot([-10, 10],[0, 0],'k:');
for co = 1:nAngle,
	plot(circle_DKL(2,co), circle_DKL(1,co), 'ko', 'MarkerSize',ourMarkerSize, ...
		'MarkerEdgeColor', circle_settings(:,co), 'MarkerFaceColor', circle_settings(:,co));
end;
plot([0,1],[0,c1(1)],'-','Color',circle_settings(:,index1));
plot([0,-1],[0,-c2(1)],'-','Color',circle_settings(:,index2));
plot([0,-1],[0,-c3(1)],'-','Color',circle_settings(:,index3));
plot([0,1],[0,c4(1)],'-','Color',circle_settings(:,index4));

%%  錐体コントラスト空間にプロットしてみる
circle_coneContrast = circle_coneInc ./ repmat(bg_cones,1,nAngle);

figure;
plot([0, 0], [-1,1],'k:');
hold on;
plot([-1, 1],[0, 0],'k:');
for co = 1:nAngle,
	plot(circle_coneContrast(1,co), circle_coneContrast(2,co), 'ko', 'MarkerSize',ourMarkerSize, ...
		'MarkerEdgeColor', circle_settings(:,co), 'MarkerFaceColor', circle_settings(:,co));
end;
axis([-1 1 -1 1]);
axis square;
title('Cone Contrast (\DeltaL/L vs. \DeltaM/M)', 'FontSize',16);
xlabel('\DeltaL/L', 'FontSize', 14);
ylabel('\DeltaM/M', 'FontSize', 14);

index1 = 2;
c1 = polyfit([0,circle_coneContrast(1,index1)],[0,circle_coneContrast(2,index1)],1);
index2 = 18;
c2 = polyfit([0,circle_coneContrast(1,index2)],[0,circle_coneContrast(2,index2)],1);
index3 = 20;
c3 = polyfit([0,circle_coneContrast(1,index3)],[0,circle_coneContrast(2,index3)],1);
index4 = 36;
c4 = polyfit([0,circle_coneContrast(1,index4)],[0,circle_coneContrast(2,index4)],1);

fill([0, 1, 1, -1, 0], [0, c1(1), -1, -c4(1), 0],[0.9,0.9,0.9]);
fill([0, 1, -1, -1, 0], [0, c2(1), 1, -c3(1), 0],[0.9,0.9,0.9]);

% プロットを再度描いて、見栄えを良くする
plot([0, 0], [-1,1],'k:');
plot([-1, 1],[0, 0],'k:');
for co = 1:nAngle,
	plot(circle_coneContrast(1,co), circle_coneContrast(2,co), 'ko', 'MarkerSize',ourMarkerSize, ...
		'MarkerEdgeColor', circle_settings(:,co), 'MarkerFaceColor', circle_settings(:,co));
end;
plot([0,1],[0,c1(1)],'-','Color',circle_settings(:,index1));
plot([0,1],[0,c2(1)],'-','Color',circle_settings(:,index2));
plot([0,-1],[0,-c3(1)],'-','Color',circle_settings(:,index3));
plot([0,-1],[0,-c4(1)],'-','Color',circle_settings(:,index4));

%%

% DKL空間でのL-M軸付近のごく狭い範囲（グレーで塗った領域）が、
% 錐体コントラスト空間では広がって表現されていることがわかる
