Skip to content

Commit 90040cb

Browse files
michaeabmichaeab
committed
2 parents 848b5c0 + 53cbe7b commit 90040cb

File tree

3 files changed

+262
-129
lines changed

3 files changed

+262
-129
lines changed

‎code/analysis/monitorGamutInLSPlane/LSDstimulusGeneration.m‎

Lines changed: 0 additions & 128 deletions
This file was deleted.

‎code/analysis/monitorGamutInLSPlane/makeMonitorGamutFigure.m‎

Lines changed: 261 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,261 @@
1+
% Make a figure showing our monitor gamut in the LS contrast plane
2+
3+
%% Initialize
4+
clear; close all;
5+
6+
%% Where to write figure
7+
figureDir = getpref('ColorTracking','figureSavePath');
8+
9+
%% Load typical calibration file from the experiment
10+
% whichCalFile = 'ViewSonicG220fb_670.mat';
11+
whichCalFile = 'ViewSonicG220fb.mat';
12+
whichCalNumber = 1;
13+
resourcesDir = getpref('ColorTracking','CalDataFolder');
14+
load(fullfile(resourcesDir,whichCalFile),'cals');
15+
cal = cals{whichCalNumber};
16+
17+
%% Cone cal object
18+
calObjCones = ObjectToHandleCalOrCalStruct(cal);
19+
20+
% Make this a 12-bit device as far as the calibration file goes.
21+
%
22+
% HARDWARE FOR DETECTION WAS 14-BIT, BUT FOR TRACKING WAS LOWER.
23+
nDeviceBits = 12;
24+
nDeviceLevels = 2^nDeviceBits;
25+
CalibrateFitGamma(calObjCones, nDeviceLevels);
26+
nPrimaries = calObjCones.get('nDevices');
27+
28+
% Can change gammaMethod to 1 for unquantized analysis
29+
gammaMethod = 2;
30+
SetGammaMethod(calObjCones,gammaMethod);
31+
32+
% Set wavelength support.
33+
Scolor = calObjCones.get('S');
34+
35+
% Zero out ambient?
36+
% [HOW WAS THIS SET FOR THE EXPERIMENT???]
37+
NOAMBIENT = false;
38+
if (NOAMBIENT)
39+
calObjCones.set('P_ambient',zeros(size(calObjCones.get('P_ambient'))));
40+
end
41+
42+
% Cone fundamentals. Methods say cones were SS 2-deg. That is
43+
% what was loaded in the tracking code. In the detection code,
44+
% it was ComputeObserverFundamentals with parameters below, which
45+
% are not quite T_cones_ss2.
46+
%
47+
% psiParamsStruct.coneParams = DefaultConeParams('cie_asano');
48+
% psiParamsStruct.coneParams.fieldSizeDegrees = 2;
49+
% psiParamsStruct.coneParams.ageYears = 30;
50+
%T_cones = ComputeObserverFundamentals(psiParamsStruct.coneParams,Scolor);
51+
% [WHICH CONES WERE REALLY USED? SUBJECT AGES WERE 28, 29, and 33]
52+
load T_cones_ss2
53+
T_cones = SplineCmf(S_cones_ss2,T_cones_ss2,Scolor);
54+
SetSensorColorSpace(calObjCones,T_cones,Scolor);
55+
56+
%% XYZ cal object
57+
calObjXYZ = ObjectToHandleCalOrCalStruct(cal);
58+
59+
% Get gamma correct
60+
CalibrateFitGamma(calObjXYZ, nDeviceLevels);
61+
SetGammaMethod(calObjXYZ,gammaMethod);
62+
if (NOAMBIENT)
63+
calObjXYZ.set('P_ambient',zeros(size(calObjCones.get('P_ambient'))));
64+
end
65+
66+
%% Image scale factor
67+
%
68+
% We report a nominal max contrast which is the contrast that would have been
69+
% shown had the 100% contrast Gabor image been in cosine phase. But ours
70+
% were in sin phase. The maximum value of the Gabor was 0.9221. We divide
71+
% our maximum gamut contrasts by this to get the real maximum contrast we
72+
% report.
73+
imageScaleFactor = 0.9221;
74+
75+
% XYZ
76+
USE1931XYZ = true;
77+
if (USE1931XYZ)
78+
load T_xyz1931.mat
79+
T_xyz = 683*SplineCmf(S_xyz1931,T_xyz1931,Scolor);
80+
else
81+
load T_xyzCIEPhys2.mat
82+
T_xyz = 683*SplineCmf(S_xyzCIEPhys2,T_xyzCIEPhys2,Scolor);
83+
end
84+
SetSensorColorSpace(calObjXYZ,T_xyz,Scolor);
85+
86+
%% Compute ambient
87+
ambientCones = SettingsToSensor(calObjCones,[0 0 0]');
88+
ambientXYZ = SettingsToSensor(calObjXYZ,[0 0 0']');
89+
90+
%% Compute the background, taking quantization into account
91+
%
92+
% The paper says the background was 30.75 cd/m2, x = 0.326, y = 0.372;
93+
SPECIFIEDBG = false;
94+
if (SPECIFIEDBG)
95+
bgxyYTarget = [0.326, 0.372 30.75]';
96+
bgXYZTarget = xyYToXYZ(bgxyYTarget);
97+
bgPrimary = SettingsToPrimary(calObjXYZ,SensorToSettings(calObjXYZ,bgXYZTarget));
98+
bgXYZ = PrimaryToSensor(calObjXYZ,bgPrimary);
99+
else
100+
bgPrimary = SettingsToPrimary(calObjCones,PrimaryToSettings(calObjCones,[0.5 0.5 0.5]'));
101+
bgXYZ = PrimaryToSensor(calObjXYZ,bgPrimary);
102+
end
103+
bgxyY = XYZToxyY(bgXYZ);
104+
bgCones = PrimaryToSensor(calObjCones,bgPrimary);
105+
fprintf('\nBackground x,y = %0.4f, %0.4f\n',bgxyY(1),bgxyY(2));
106+
fprintf('Background Y = %0.2f cd/m2, ambient %0.3f cd/m2\n\n',bgXYZ(2),ambientXYZ(2));
107+
108+
%% Max contrast
109+
%
110+
% Find maximum in gamut contrast for a set of color directions. We
111+
% are not going to worry about device quantization since we used
112+
% high-bit depth hardware
113+
nAngles = 1000;
114+
theAngles = linspace(0,2*pi,nAngles);
115+
for aa = 1:nAngles
116+
% Get a unit contrast vector at the specified angle
117+
unitContrastDir = [cos(theAngles(aa)) 0 sin(theAngles(aa))]';
118+
119+
% Convert from cone contrast to cone excitation direction.
120+
% Don't care about length here as that is handled by the contrast
121+
% maximization code below.
122+
unitConesDir = unitContrastDir .* bgCones;
123+
124+
% Convert the direction to the desired direction in primary space.
125+
% Since this is desired, we do not go into settings here. Adding
126+
% and subtracting the background handles the ambient correctly.
127+
unitPrimaryDir = SensorToPrimary(calObjCones,unitConesDir + bgCones) - SensorToPrimary(calObjCones,bgCones);
128+
129+
% Find out how far we can go in the desired direction and scale the
130+
% unitPrimaryDir by that amount
131+
[s,sPos,sNeg] = MaximizeGamutContrast(unitPrimaryDir,bgPrimary);
132+
gamutPrimaryDir = s*unitPrimaryDir;
133+
if (any(gamutPrimaryDir+bgPrimary < -1e-3) | any(gamutPrimaryDir+bgPrimary > 1+1e-3))
134+
error('Somehow primaries got too far out of gamut\n');
135+
end
136+
if (any(-gamutPrimaryDir+bgPrimary < -1e-3) | any(-gamutPrimaryDir+bgPrimary > 1+1e-3))
137+
error('Somehow primaries got too far out of gamut\n');
138+
end
139+
gamutDevPos1 = abs(gamutPrimaryDir+bgPrimary - 1);
140+
gamutDevNeg1 = abs(gamutPrimaryDir+bgPrimary);
141+
gamutDevPos2 = abs(-gamutPrimaryDir+bgPrimary - 1);
142+
gamutDevNeg2 = abs(-gamutPrimaryDir+bgPrimary);
143+
gamutDev = min([gamutDevPos1 gamutDevNeg1 gamutDevPos2 gamutDevNeg2]);
144+
if (gamutDev > 1e-3)
145+
error('Did not get primaries close enough to gamut edge');
146+
end
147+
148+
% Get the settings that as closely as possible approximate what we
149+
% want. One of these should be very close to 1 or 0, and none should
150+
% be less than 0 or more than 1.
151+
gamutSettings = PrimaryToSettings(calObjCones,gamutPrimaryDir + bgPrimary);
152+
if (any(gamutSettings < 0) | any(gamutSettings > 1))
153+
error('Somehow settings got out of gamut\n');
154+
end
155+
156+
% Figure out the cone excitations for the settings we computed, and
157+
% then convert to contrast as our maximum contrast in this direction.
158+
%
159+
% Dividing by imageScaleFactor handles the sine phase of the Gabor
160+
gamutCones = SettingsToSensor(calObjCones,gamutSettings);
161+
gamutContrast(:,aa) = ((gamutCones-bgCones) ./ bgCones)/imageScaleFactor;
162+
vectorLengthContrast(aa) = norm(gamutContrast(:,aa));
163+
end
164+
165+
%gamutContrast
166+
167+
% Make a plot of the gamut in the LS contrast plane
168+
figure; clf; hold on;
169+
%plot(100*gamutContrast(1,:),100*gamutContrast(3,:),'ko','MarkerSize',8);
170+
plot([-100 100],[0 0],'k:','LineWidth',0.5);
171+
plot([0 0],[-100 100],'k:','LineWidth',0.5);
172+
plot(100*gamutContrast(1,:),100*gamutContrast(3,:),'k','LineWidth',2);
173+
xlim([-100 100])
174+
ylim([-100 100]);
175+
axis('square');
176+
xlabel('L Cone Contrast (%)')
177+
ylabel('S Cone Contrast (%)');
178+
saveas(gcf,fullfile(figureDir,'MonitorGamutFigure.pdf'),'pdf')
179+
180+
% This does the same computation for specific stimulus angles
181+
theSpecificAngles = [0 90.0000 75.0000 -75.0000 45.0000 -45.0000 78.7500 82.5000 86.2000 -78.7500 -82.5000 -86.2000 89.6000 88.6000 87.6000 22.5000 -1.4000 -22.5000];
182+
for aa = 1:length(theSpecificAngles)
183+
184+
% Convert from cone contrast to cone excitation direction.
185+
% Don't care about length here as that is handled by the contrast
186+
% maximization code below.
187+
unitContrastDir = [cosd(theSpecificAngles(aa)) 0 sind(theSpecificAngles(aa))]';
188+
189+
% Convert from cone contrast to cone excitation direction.
190+
% Don't care about length here as that is handled by the contrast
191+
% maximization code below.
192+
unitConesDir = unitContrastDir .* bgCones;
193+
194+
% Convert the direction to the desired direction in primary space.
195+
% Since this is desired, we do not go into settings here. Adding
196+
% and subtracting the background handles the ambient correctly.
197+
unitPrimaryDir = SensorToPrimary(calObjCones,unitConesDir + bgCones) - SensorToPrimary(calObjCones,bgCones);
198+
199+
% Find out how far we can go in the desired direction and scale the
200+
% unitPrimaryDir by that amount
201+
[s,sPos(aa),sNeg(aa)] = MaximizeGamutContrast(unitPrimaryDir,bgPrimary);
202+
if (sPos(aa) < sNeg(aa))
203+
gamutLimitSign(aa) = 1;
204+
else
205+
gamutLimitSign(aa) = -1;
206+
end
207+
gamutPrimaryDir = s*unitPrimaryDir;
208+
if (any(gamutPrimaryDir+bgPrimary < -1e-3) | any(gamutPrimaryDir+bgPrimary > 1+1e-3))
209+
error('Somehow primaries got too far out of gamut\n');
210+
end
211+
if (any(-gamutPrimaryDir+bgPrimary < -1e-3) | any(-gamutPrimaryDir+bgPrimary > 1+1e-3))
212+
error('Somehow primaries got too far out of gamut\n');
213+
end
214+
gamutDevPos1 = abs(gamutPrimaryDir+bgPrimary - 1);
215+
gamutDevNeg1 = abs(gamutPrimaryDir+bgPrimary);
216+
gamutDevPos2 = abs(-gamutPrimaryDir+bgPrimary - 1);
217+
gamutDevNeg2 = abs(-gamutPrimaryDir+bgPrimary);
218+
gamutDev = min([gamutDevPos1 gamutDevNeg1 gamutDevPos2 gamutDevNeg2]);
219+
if (gamutDev > 1e-3)
220+
error('Did not get primaries close enough to gamut edge');
221+
end
222+
223+
% Get the settings that as closely as possible approximate what we
224+
% want. One of these should be very close to 1 or 0, and none should
225+
% be less than 0 or more than 1.
226+
gamutSettings = PrimaryToSettings(calObjCones,gamutPrimaryDir + bgPrimary);
227+
if (any(gamutSettings < 0) | any(gamutSettings > 1))
228+
error('Somehow settings got out of gamut\n');
229+
end
230+
231+
% Figure out the cone excitations for the settings we computed, and
232+
% then convert to contrast as our maximum contrast in this direction.
233+
%
234+
% Dividing by imageScaleFactor handles the sine phase of the Gabor
235+
gamutCones = SettingsToSensor(calObjCones,gamutSettings);
236+
specificGamutContrast(:,aa) = ((gamutCones-bgCones) ./ bgCones)/imageScaleFactor;
237+
specificVectorLengthContrast(aa) = norm(specificGamutContrast(:,aa));
238+
fprintf('Angle %0.1f, L cone contrast %0.3f%%, M, %0.3f%%, S %0.3f%%, vector length %0.1f%%\n', ...
239+
theSpecificAngles(aa),100*specificGamutContrast(1,aa),100*specificGamutContrast(2,aa),100*specificGamutContrast(3,aa), ...
240+
100*specificVectorLengthContrast(aa));
241+
end
242+
specificVectorLengthContrast
243+
244+
% CF MAB data from Tracking
245+
%
246+
% uniqueColorDirs(:)'
247+
%
248+
% ans =
249+
%
250+
% 0 90.0000 75.0000 -75.0000 45.0000 -45.0000 78.7500 82.5000 86.2000 -78.7500 -82.5000 -86.2000 89.6000 88.6000 87.6000 22.5000 -1.4000 -22.5000
251+
%
252+
% matrixContrasts(1,:)
253+
%
254+
% ans =
255+
%
256+
% 0.1800 0.8500 0.6500 0.7800 0.2500 0.2600 0.8300 0.8500 0.8500 0.8400 0.8400 0.8400 0.8500 0.8500 0.8500 0.1900 0.1800 0.1900
257+
%
258+
% cals{1} - August 31 cal, used for tracking experiment
259+
% specificVectorLengthContrast =
260+
%
261+
% 0.1863 0.8513 0.6556 0.7990 0.2569 0.2710 0.8438 0.8734 0.8605 0.8462 0.8442 0.8458 0.8521 0.8542 0.8567 0.1999 0.1865 0.2039

0 commit comments

Comments
 (0)