% 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 primary = (0.5,0.5,0.5)として、錐体刺激値を計算する
bg_Primary = [0.5 0.5 0.5]';
bg_cones = PrimaryToSensor(cal_cones,bg_Primary);

% MB色度図計算用の準備
% 上で設定した背景に対してsMB=1.0とする
s_scaler = 1 * 683;
bg_ls = OurLMSToMacBoyn2(bg_cones,T_cones,T_Y,s_scaler);
s_scaler = 683/bg_ls(2);

% 新しい係数で計算し直し。
bg_ls = OurLMSToMacBoyn2(bg_cones,T_cones,T_Y,s_scaler);

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

%% DKL空間の L-M vs. S 平面で円を描き、それをMacLeod-Boynton色度図にプロットしてみる

angle = 0:10:350;
nAngle = length(angle);
radius = 0.7;
circle_DKL_polar = [radius*ones(1,nAngle); angle; zeros(1,nAngle)];
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);

circle_ls = OurLMSToMacBoyn2(circle_cones,T_cones,T_Y,s_scaler);

%%  DKL空間の L-M vs. S 平面上の円
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(3,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. S plane', 'FontSize',16);
xlabel('L-M', 'FontSize', 14);
ylabel('S', 'FontSize', 14);

%%  MacLeod-Boynton色度図上にプロット(縦横等スケール)
figure;
hold on;
for co = 1:nAngle,
	plot(circle_ls(1,co), circle_ls(2,co), 'ko', 'MarkerSize',ourMarkerSize, ...
		'MarkerEdgeColor', circle_settings(:,co), 'MarkerFaceColor', circle_settings(:,co));
end;
scaler = 1.2;
halfRange = scaler*(max(circle_ls,[],2)-bg_ls);
axis([bg_ls(1)-halfRange(1) bg_ls(1)+halfRange(1) bg_ls(2)-halfRange(2) bg_ls(2)+halfRange(2)]);
% axis square;
axis equal;
set(gca,'xtick', -0.4:0.2:2);
set(gca,'ytick', 0.4:0.2:2.4);
lims = axis;
plot([bg_ls(1), bg_ls(1)], [lims(3),lims(4)],'k:');
plot([lims(1), lims(2)],[bg_ls(2), bg_ls(2)],'k:');

% 見栄えを良くするため、もう一度プロット
for co = 1:nAngle,
	plot(circle_ls(1,co), circle_ls(2,co), 'ko', 'MarkerSize',ourMarkerSize, ...
		'MarkerEdgeColor', circle_settings(:,co), 'MarkerFaceColor', circle_settings(:,co));
end;
title('MacLeod-Boynton diagram (axis equal)', 'FontSize',16);
xlabel('l_M_B', 'FontSize', 16');
ylabel('s_M_B', 'FontSize', 16');

%% 

% 縦軸と横軸が等スケールだと、プロットは縦長となる


%%  MacLeod-Boynton色度図上にプロット(円になるようにプロット)
figure;
hold on;
for co = 1:nAngle,
	plot(circle_ls(1,co), circle_ls(2,co), 'ko', 'MarkerSize',ourMarkerSize, ...
		'MarkerEdgeColor', circle_settings(:,co), 'MarkerFaceColor', circle_settings(:,co));
end;
scaler = 1.2;
halfRange = scaler*(max(circle_ls,[],2)-bg_ls);
axis([bg_ls(1)-halfRange(1) bg_ls(1)+halfRange(1) bg_ls(2)-halfRange(2) bg_ls(2)+halfRange(2)]);
axis square;
lims = axis;
plot([bg_ls(1), bg_ls(1)], [lims(3),lims(4)],'k:');
plot([lims(1), lims(2)],[bg_ls(2), bg_ls(2)],'k:');

% 見栄えを良くするため、もう一度プロット
for co = 1:nAngle,
	plot(circle_ls(1,co), circle_ls(2,co), 'ko', 'MarkerSize',ourMarkerSize, ...
		'MarkerEdgeColor', circle_settings(:,co), 'MarkerFaceColor', circle_settings(:,co));
end;
title('MacLeod-Boynton diagram (axis square)', 'FontSize',16);
xlabel('l_M_B', 'FontSize', 16');
ylabel('s_M_B', 'FontSize', 16');

%%

% スケールを調整すると、きれいな円になる。
% DKLの結果と比較してみてください。
