% KisoHB_MB_diagram_demo.m
% MB色度図内にスペクトル軌跡を描き、等エネルギー白色点をプロットする

clear all; close all;

%% MB色度図内にスペクトル軌跡を描く
whichCones = 'StockmanSharpe';
[S_cones, T_cones, S_Y, T_Y] = SetFundamentals2(whichCones);

% sMBのピークが1.0になるようにスケーリング（from CVRL for SS fundamentals）
% このようにスケーリングすると、lMBもsMBも、可視スペクトルに対して0.0～1.0の値を
% とることになる。
[val,index] = max(T_cones(3,:)./T_Y);
s_scaler = 1/(T_cones(3,index)/T_Y(index));		% s_scaler = 0.0371597 * 683;

% 各波長の単色光に対するMB色度座標を計算する
ls = OurLMSToMacBoyn2(T_cones,T_cones,T_Y,s_scaler);

% プロットしてみる
ourStep = 1:4:S_cones(3);
figure;
plot(ls(1,:),ls(2,:),'k-'); hold on;
plot(ls(1,ourStep),ls(2,ourStep),'ko');
axis([-0.02 1.02 -0.02 1.02]);
axis equal;
xlabel('l_M_B', 'FontSize', 16');
ylabel('s_M_B', 'FontSize', 16');

% 等エネルギー白色点をプロット
spd_EEW = ones(S_cones(3),1);
cones_EEW = T_cones*spd_EEW;
ls_EEW = OurLMSToMacBoyn2(cones_EEW,T_cones,T_Y,s_scaler);
plot(ls_EEW(1),ls_EEW(2),'ks','MarkerEdgeColor','k','MarkerFaceColor',[0,0,0]);

% 標準照明光D65についても計算
% load spd_D65
% spd_D65 = SplineCmf(S_D65,spd_D65',S_cones);
% spd_D65 = spd_D65';
% cones_D65 = T_cones*spd_D65;
% ls_D65 = OurLMSToMacBoyn2(cones_D65,T_cones,T_Y,s_scaler);
% plot(ls_D65(1),ls_D65(2),'ro','MarkerEdgeColor','r','MarkerFaceColor',[1.0,0,0]);

%%

% sMBのピークが1.0になるようにスケーリングすると、lMBもsMBも、
% 可視スペクトルに対して0.0～1.0の値をとることになる。
%
% 範囲が一緒という意味では公平なスケーリングとも言えそうだが、
% これだと、等エネルギー白色（黒い四角）が一番下の方にプロットされ、白色点より下の黄色領域が
% ほとんど表示されていないことになる。


%% 白色点付近を拡大

ourStep = 1:4:S_cones(3);
figure;
plot(ls(1,:),ls(2,:),'k-'); hold on;
plot(ls(1,ourStep),ls(2,ourStep),'ko');
axis([0.4 1.02 -0.02 0.2]);
axis square;
xlabel('l_M_B', 'FontSize', 16');
ylabel('s_M_B', 'FontSize', 16');
plot(ls_EEW(1),ls_EEW(2),'ks','MarkerEdgeColor','k','MarkerFaceColor',[0,0,0]);

%%

% 白色点付近を拡大してみると、白色点に対するsMBの値が非常に小さいことがわかる。


%% 等エネルギー白色に対してsMB=1.0とする。monitor primaryもプロット。

% Boynton & Kambe (1980) 方式（オリジナルでは網膜照度を使用）
% Boynton, R.M., and Kambe, N. (1980)
% Chromatic difference steps of moderate size measured along theoretically critical axes.
% Color Research and Application, 5, 13-23.
% 
% 等エネルギー白色を背景野として使用する場合には、背景野でsMB=1.0となるので都合が良い。

% s_scalerを1として、等エネルギー白色のsMB値を計算し、係数を決定
s_scaler = 1 * 683;
spd_EEW = ones(S_cones(3),1);
cones_EEW = T_cones*spd_EEW;
ls_EEW = OurLMSToMacBoyn2(cones_EEW,T_cones,T_Y,s_scaler);
s_scaler = 683/ls_EEW(2);

% 新たなスケーリング係数で再計算
ls = OurLMSToMacBoyn2(T_cones,T_cones,T_Y,s_scaler);
ls_EEW = OurLMSToMacBoyn2(cones_EEW,T_cones,T_Y,s_scaler);

% プロットしてみる
ourStep = 1:4:S_cones(3);
figure;
plot(ls(1,:),ls(2,:),'k-'); hold on;
plot(ls(1,ourStep),ls(2,ourStep),'ko');
xlabel('l_M_B', 'FontSize', 16');
ylabel('s_M_B', 'FontSize', 16');
plot(ls_EEW(1),ls_EEW(2),'ks','MarkerEdgeColor','k','MarkerFaceColor',[0,0,0]);

% PTB3のキャリブレーションファイルから、モニターのRGB蛍光体の色度座標をプロットしてみる。
% キャリブレーションデータの読み込みと必要な設定を行う
cal = LoadCalFile('PTB3TestCal');
% S = cal.S_device;
% cal.T_device = WlsToT(S);
% cal.T_ambient = WlsToT(S);
cal = SetGammaMethod(cal,0);
cal_cones = SetSensorColorSpace(cal,T_cones,S_cones);

% RGBがそれぞれmax(1.0)のときのMB色度座標を計算してプロット
settings_RGB = eye(3);
cones_RGB = SettingsToSensor(cal_cones, settings_RGB);
ls_RGB = OurLMSToMacBoyn2(cones_RGB,T_cones,T_Y,s_scaler);
plot([ls_RGB(1,:) ls_RGB(1,1)],[ls_RGB(2,:) ls_RGB(2,1)],'k');
plot(ls_RGB(1,1),ls_RGB(2,1),'ro','MarkerEdgeColor','r','MarkerFaceColor',[1,0,0]);
plot(ls_RGB(1,2),ls_RGB(2,2),'go','MarkerEdgeColor','g','MarkerFaceColor',[0,1,0]);
plot(ls_RGB(1,3),ls_RGB(2,3),'bo','MarkerEdgeColor','b','MarkerFaceColor',[0,0,1]);
axis([0.5 0.9 -0.02 8]);
axis square;

%%

% やはり、白から青方向が相対的に広い。
%
% MB色度図はこういう特徴をもつのは、
% S錐体の分光感度の位置が、M, L錐体から短波長方向にかなり離れているため。
%

%% 比較のため、xy色度図で同様の計算を行う
load T_xyz1931;

% 等エネルギー白色の分光エネルギー分布を設定
spd_EEW = ones(S_xyz1931(3),1);

% 各波長の単色光に対するxy色度座標を計算する
CIE_x = T_xyz1931(1,:) ./ sum(T_xyz1931);
CIE_y = T_xyz1931(2,:) ./ sum(T_xyz1931);

% プロットしてみる
ourStep = 1:4:S_xyz1931(3);
figure;
plot(CIE_x,CIE_y,'k-'); hold on;
plot(CIE_x(ourStep),CIE_y(ourStep),'ko');
axis equal;
axis([-0.02 1.02 -0.02 1.02]);
xlabel('x', 'FontSize', 16');
ylabel('y', 'FontSize', 16');

% モニターのRGB蛍光体の色度座標を計算し、プロット
cal_XYZ = SetSensorColorSpace(cal,T_xyz1931,S_xyz1931);
XYZ_RGB = SettingsToSensor(cal_XYZ, settings_RGB);
xyY_RGB = XYZToxyY(XYZ_RGB);
plot([xyY_RGB(1,:) xyY_RGB(1,1)],[xyY_RGB(2,:) xyY_RGB(2,1)],'k');
plot(xyY_RGB(1,1),xyY_RGB(2,1),'ro','MarkerEdgeColor','r','MarkerFaceColor',[1,0,0]);
plot(xyY_RGB(1,2),xyY_RGB(2,2),'go','MarkerEdgeColor','g','MarkerFaceColor',[0,1,0]);
plot(xyY_RGB(1,3),xyY_RGB(2,3),'bo','MarkerEdgeColor','b','MarkerFaceColor',[0,0,1]);

% 等エネルギー白色の色度座標を計算し、プロット
XYZ_EEW = T_xyz1931*spd_EEW;
xyY_EEW = XYZToxyY(XYZ_EEW);
plot(xyY_EEW(1),xyY_EEW(2),'ks','MarkerEdgeColor','k','MarkerFaceColor',[0,0,0]);

% 標準照明光D65についても計算
% load spd_D65
% xyY_D65 = XYZToxyY(T_xyz1931*spd_D65);
% plot(xyY_D65(1),xyY_D65(2),'ro','MarkerEdgeColor','r','MarkerFaceColor',[1.0,0,0]);

%%

% RGBの位置がMB色度図と異なることに注意
% 
% 等エネルギー白色点が、モニターで表示できる色範囲の真ん中付近にある。
%

%% おまけで、uv色度図にもプロットしてみる
figure;

% XYZ値をuvY値に変換する（スペクトル軌跡用）
uvY_CIE = XYZTouvY(T_xyz1931);

% プロットしてみる
ourStep = 1:4:S_xyz1931(3);
plot(uvY_CIE(1,:),uvY_CIE(2,:),'k-'); hold on;
plot(uvY_CIE(1,ourStep),uvY_CIE(2,ourStep),'ko');
axis equal;
axis([-0.02 0.9 -0.02 0.7]);
xlabel('u', 'FontSize', 16');
ylabel('v', 'FontSize', 16');

% モニターのRGB蛍光体のuv色度座標を計算し、プロット
uvY_RGB = XYZTouvY(XYZ_RGB);
plot([uvY_RGB(1,:) uvY_RGB(1,1)],[uvY_RGB(2,:) uvY_RGB(2,1)],'k');
plot(uvY_RGB(1,1),uvY_RGB(2,1),'ro','MarkerEdgeColor','r','MarkerFaceColor',[1,0,0]);
plot(uvY_RGB(1,2),uvY_RGB(2,2),'go','MarkerEdgeColor','g','MarkerFaceColor',[0,1,0]);
plot(uvY_RGB(1,3),uvY_RGB(2,3),'bo','MarkerEdgeColor','b','MarkerFaceColor',[0,0,1]);

% 等エネルギー白色の色度座標を計算し、プロット
uvY_EEW = XYZTouvY(XYZ_EEW);
plot(uvY_EEW(1),uvY_EEW(2),'ko','MarkerEdgeColor','k','MarkerFaceColor',[0,0,0]);

%%

% xy色度図と比較すると、白から青方向が拡大されたことがわかる。
