Contents

% KisoHB_DKLDemo
%
% DKLDemo
% https://github.com/Psychtoolbox-3/Psychtoolbox-3/blob/master/Psychtoolbox/PsychDemos/DKLDemo.m
% を実験心理学ハンドブック用に変更
%
% 04/08/15	ek		DKLDemo -->  KisoHB_DKLDemo

% Clear
clear; close all

錐体分光感度の設定

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

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

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

Brainard (1996)に載っている例題で検算してみる

bgLMS = [2 4 3]';
test_LMSInc = [2 -2.5 1]';
dkl = OurConeIncToDKL(test_LMSInc, bgLMS, T_cones, T_Y);
fprintf('L+M = %g, L-M = %g, S = %g\n', dkl(1), dkl(2), dkl(3));

% 正しければ、以下のように表示されるはず。
% L+M = -0.144338, L-M = 1.2112, S = 0.416667
%
% 分光視感効率へのスケーリングを行わないと、異なる値となることを確認。

L+M = -0.144338, L-M = 1.2112, S = 0.416667

DKL_Demoの計算：キャリブレーションデータを使って、実際に計算してみる

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

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

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

% incremental cone coordinatesをDKL座標に変換する行列を求める
% M_ConeIncToDKLが、Brainard (1996)の式(A.4.11)の変換行列
M_ConeIncToDKL = ComputeDKL_M(bgLMS,T_cones,T_Y);
% M_DKLToConeInc、Brainard (1996)の式(A.4.12)の変換行列
% M_DKLToConeInc = inv(M_ConeIncToDKL);

% 上記の変換行列は、pooled cone contrastに基づくスケーリングがなされた結果を表す。
% これでよければ、この変換行列をそのまま利用する（OurConeIncToDKL.mを参照）。
%
% ディスプレイで表示できる最大値に基づきスケーリングを行う例がDKLDemo.mにあるので、
% 参照のこと。
% それを関数の形に直したのがOurComputeDKL_M_Device.m
%

ディスプレイで表示できる最大値に基づきスケーリングした変換行列で計算

% incremental cone coordinatesをDKL座標に変換する行列を求める
% 今度は、ディスプレイで表示できる最大値に基づきスケーリングしたもの
M_ConeIncToDKL_Device = OurComputeDKL_M_Device(bgLMS,T_cones,T_Y,calLMS);

% この座標はディスプレイに依存するので、L+M, L-M, S-(L+M)のそれぞれが
% 1.0の時のxyYを求めておく
coneInc = M_ConeIncToDKL_Device\eye(3);
DKL_unit_LMS = coneInc + repmat(bgLMS,1,3);

% LMSからディスプレイのRGB primary coordinatesを求める
DKL_unit_Primary = SensorToPrimary(calLMS, DKL_unit_LMS);

% xyY計算用にデータを読み込み、計算
load T_xyz1931;
T_xyz = SplineCmf(S_xyz1931,683*T_xyz1931,S);
cal_XYZ = SetSensorColorSpace(cal,T_xyz,S);
DKL_unit_xyY = XYZToxyY(PrimaryToSensor(cal_XYZ, DKL_unit_Primary));
bg_xyY = XYZToxyY(PrimaryToSensor(cal_XYZ,bgPrimary));

% 座標値の表示
fprintf('\nxyY coordinates\n');
fprintf('1.0 L+M = (%.4f, %.4f, %.4f)\n', DKL_unit_xyY(1,1),DKL_unit_xyY(2,1),DKL_unit_xyY(3,1));
fprintf('1.0 L-M = (%.4f, %.4f, %.4f)\n', DKL_unit_xyY(1,2),DKL_unit_xyY(2,2),DKL_unit_xyY(3,2));
fprintf('1.0 S   = (%.4f, %.4f, %.4f)\n', DKL_unit_xyY(1,3),DKL_unit_xyY(2,3),DKL_unit_xyY(3,3));
fprintf('Bkg     = (%.4f, %.4f, %.4f)\n', bg_xyY(1),bg_xyY(2),bg_xyY(3));

xyY coordinates
1.0 L+M = (0.2825, 0.2990, 144.4568)
1.0 L-M = (0.3285, 0.2797, 77.0038)
1.0 S   = (0.2502, 0.2125, 76.5901)
Bkg     = (0.2825, 0.2990, 77.0074)

% 背景とL+Mでは、色度座標が一致し、輝度だけが異なる。
% L+M軸と言っているが、L+Mのみが変化する軸ではなく、L+M過程を単離する軸
% （つまり、L+M以外のL-MとS-(L+M)が変化しない方向、deltaL : deltaM : deltaS が一定）
% であることに注意。
%
% L+MとSでは、色度が異なり、輝度が一致...しない。それはなぜかというと...

各錐体感度と対応した分光視感効率で輝度を計算する

DKL_unit_Lum = PrimaryToSensor(calLum, DKL_unit_Primary);
bg_Lum = PrimaryToSensor(calLum, bgPrimary);

fprintf('\nLuminances calculated with the modified luminous efficiency\n')
fprintf('1.0 L+M = %.4f\n', DKL_unit_Lum(1));
fprintf('1.0 L-M = %.4f\n', DKL_unit_Lum(2));
fprintf('1.0 S   = %.4f\n', DKL_unit_Lum(3));
fprintf('Bkg     = %.4f\n', bg_Lum);

Luminances calculated with the modified luminous efficiency
1.0 L+M = 145.3638
1.0 L-M = 77.4907
1.0 S   = 77.4907
Bkg     = 77.4907

% 今度は、背景とL+MとSで輝度が一致する。
% DKLの等輝度面は、CIE 1931XYZ表色系でいうところの等輝度面ではない。

極座標系での表示

% radius, azimuth, and elevation
dkl_polar = LinearDKLToPolarDeg(eye(3));

fprintf('\nPolar coordinates (Radius, Azimuth, Elevation)\n')
fprintf('1.0 L+M = (%.2f, %4.1f, %4.1f)\n', dkl_polar(1,1),dkl_polar(2,1),dkl_polar(3,1));
fprintf('1.0 L-M = (%.2f, %4.1f, %4.1f)\n', dkl_polar(1,2),dkl_polar(2,2),dkl_polar(3,2));
fprintf('1.0 S   = (%.2f, %4.1f, %4.1f)\n', dkl_polar(1,3),dkl_polar(2,3),dkl_polar(3,3));

Polar coordinates (Radius, Azimuth, Elevation)
1.0 L+M = (1.00,  0.0, 90.0)
1.0 L-M = (1.00,  0.0,  0.0)
1.0 S   = (1.00, 90.0,  0.0)

等輝度面の色表示

% 最後に、DKL空間の等輝度面を表示
% 以下の計算のほとんどは、DKLDemoより抜粋。丸く表示するようにしただけ。
imageSize = 513;
[X,Y] = meshgrid(0:imageSize-1,0:imageSize-1);
X = X-(imageSize-1)/2; Y = Y-(imageSize-1)/2;
X = X/max(abs(X(:))); Y = Y/max(abs(Y(:)));
XVec = reshape(X,1,imageSize^2); YVec = -reshape(Y,1,imageSize^2);
imageLMS = bgLMS*ones(size(XVec))+coneInc(:,2)*XVec+coneInc(:,3)*YVec;

% RGB settings（RGBの設定値）を求める
[imageRGB,badIndex] = SensorToSettings(calLMS,imageLMS);
bgRGB = SensorToSettings(calLMS,bgLMS);

vecLength = hypot(XVec, YVec);	% どうせなら、丸く表示する
imageRGB(:,find(vecLength > 1)) = bgRGB(:,ones(size(find(vecLength > 1))));

rPlane = reshape(imageRGB(1,:),imageSize,imageSize);
gPlane = reshape(imageRGB(2,:),imageSize,imageSize);
bPlane = reshape(imageRGB(3,:),imageSize,imageSize);
theImage = zeros(imageSize,imageSize,3);
theImage(:,:,1) = rPlane;
theImage(:,:,2) = gPlane;
theImage(:,:,3) = bPlane;

% 等輝度面を色付きで表示し、説明を追加
figure; clf;
image(theImage);
hold on;
plot([imageSize/2, imageSize/2], [0, imageSize], 'k:');
plot([0, imageSize], [imageSize/2, imageSize/2], 'k:');
title('DKL equiluminant plane', 'FontSize',16);
axis off; axis equal;
text(imageSize/2,imageSize+30,'L-M', 'FontSize', 18, 'HorizontalAlignment','center');
text(-50,imageSize/2,'S-(L+M)', 'FontSize', 18, 'HorizontalAlignment','center', 'Rotation',90);
axis([1 imageSize 1 imageSize]);

Published with MATLAB® R2014b