~ubuntu-branches/ubuntu/trusty/psychtoolbox-3/trusty-proposed

« back to all changes in this revision

Viewing changes to Psychtoolbox/PsychColorimetricData/ComputeRawConeFundamentals.m

  • Committer: Package Import Robot
  • Author(s): Yaroslav Halchenko
  • Date: 2013-11-19 23:34:50 UTC
  • mfrom: (3.1.4 experimental)
  • Revision ID: package-import@ubuntu.com-20131119233450-f7nf92vb8qavjmk8
Tags: 3.0.11.20131017.dfsg1-3
Upload to unsable since fresh glew has arrived to sid!

Show diffs side-by-side

added added

removed removed

Lines of Context:
1
 
function T_quantal = ComputeRawConeFundamentals(params,staticParams)
 
 
b'% T_quantal = ComputeRawConeFundamentals(params,staticParams)'
 
 
b'%'
 
 
b'% Function to compute normalized cone quantal sensitivities'
 
 
b'% from underlying pieces and parameters.'
 
 
b'%'
 
 
b'% In the passed params structure, you can either pass'
 
 
b'% the lambdaMax values for the photopigment, in which'
 
 
b'% case the absorbance is computed from the specified'
 
 
b'% nomogram, or you can pass the absorbance values'
 
 
b'% directly in T_xxx format.  A typical choice in this'
 
 
b'% case would be 10.^T_lgo10coneabsorbance_ss for the'
 
 
b'% Stockman-Sharpe/CIE estimates.'
 
 
b'%'
 
 
b"% It's not clear that it is necessary to normalize; eventually"
 
 
b'% this should probably be changed to provide fundamentals that'
 
 
b'% give actual quantal efficiency.'
 
 
b'%'
 
 
b'% The typical use of this function is to be called by '
 
 
b'% ComputeCIEConeFundamentals, which sets up the'
 
 
b'% passed structures acording to the CIE standard. '
 
 
b'% This routine, however, could in principle be used'
 
 
b'% with a wide variety of choices of the component pieces.'
 
 
b'%'
 
 
b'% 8/12/11  dhb  Starting to make this actually work.'
 
 
b'% 8/14/11  dhb  Change name, expand comment.'
 
 
b'% Handle bad value'
 
 
b'index = find(params.axialDensity <= 0.0001);'
 
 
b'if (~isempty(index))'
2
1
   params.axialDensity(index) = 0.0001;
 
 
b'end'
 
 
b'% Handle optional values'
 
 
b"if (~isfield(params,'extraLens'))"
3
2
   params.extraLens = 0;
 
 
b'end'
 
 
b"if (~isfield(params,'extraMac'))"
4
3
   params.extraMac = 0;
 
 
b'end'
 
 
b'% Prereceptor transmittance.  Sometimes adjustments of peak density'
 
 
b'% recommended by various standards push the density less than'
 
 
b'% zero at some wavelengths, so we need to make sure we never have'
 
 
b'% transmittance greater than 1.'
 
 
b'lens = 10.^-(-log10(staticParams.lensTransmittance)+params.extraLens);'
 
 
b'lens(lens > 1) = 1;'
 
 
b'mac = 10.^-(-log10(staticParams.macularTransmittance)+params.extraMac);'
 
 
b'mac(mac > 1) = 1;'
 
 
b"% Compute nomogram if absorbance wasn't passed directly.  We detect"
 
 
b'% a direct pass by the existance of params.absorbance.'
 
 
b"if (isfield(params,'absorbance'))"
5
4
   absorbance = params.absorbance;
 
 
b'else'
6
5
   absorbance = PhotopigmentNomogram(staticParams.S,params.lambdaMax,staticParams.whichNomogram);
 
 
b'end'
 
 
b'% Compute absorbtance'
 
 
b'%'
 
 
b'% Handle special case where we deal with ser/ala polymorphism for L cone'
 
 
b'if (size(absorbance,1) == 4)'
7
6
   absorbtance = AbsorbanceToAbsorbtance(absorbance,staticParams.S,...
8
7
       [params.axialDensity(1) ; params.axialDensity(1) ; ...
9
8
       params.axialDensity(2) ; params.axialDensity(3)]);
 
 
b'elseif (size(absorbance,1) == 3)'
10
9
   absorbtance = AbsorbanceToAbsorbtance(absorbance,staticParams.S,...
11
10
       [params.axialDensity(1) ; params.axialDensity(2) ; ...
12
11
       params.axialDensity(3)]);
 
 
b'else'
13
12
   error('Unexpected number of photopigment lambda max values passed');
 
 
b'end'
 
 
b'%% Put together pre-receptor and receptor parts'
 
 
b'for i = 1:size(absorbtance,1)'
14
13
   absorbtance(i,:) = absorbtance(i,:) .* lens .* mac;
 
 
b'end'
 
 
b'% Put it into the right form'
 
 
b'T_quantal = zeros(3,staticParams.S(3));'
 
 
b'if (size(absorbtance,1) == 4)'
15
14
   T_quantal(1,:) = staticParams.LserWeight*absorbtance(1,:) + ...
16
15
       (1-staticParams.LserWeight)*absorbtance(2,:);
17
16
   T_quantal(2,:) = absorbtance(3,:);
18
17
   T_quantal(3,:) = absorbtance(4,:);
 
 
b'elseif (size(absorbtance,1) == 3)'
19
18
   T_quantal(1,:) = absorbtance(1,:);
20
19
   T_quantal(2,:) = absorbtance(2,:);
21
20
   T_quantal(3,:) = absorbtance(3,:);
 
 
b'else'
22
21
   error('Unexpected number of photopigment lambda max values passed');
 
 
b'end'
 
 
b'% Normalize to max of one for each receptor'
 
 
b'for i = 1:size(T_quantal,1)'
23
22
   T_quantal(i,:) = T_quantal(i,:)/max(T_quantal(i,:));
 
 
b'end'
 
 
b'\\ No newline at end of file'
 
23
function [T_quantalAbsorptionsNormalized,T_quantalAbsorptions,T_quantalIsomerizations] = ComputeRawConeFundamentals(params,staticParams)
 
24
% [T_quantalAbsorptionsNormalized,T_quantalAbsorptions,T_quantalIsomerizations] = ComputeRawConeFundamentals(params,staticParams)
 
25
%
 
26
% Function to compute normalized cone quantal sensitivities
 
27
% from underlying pieces and parameters.
 
28
%
 
29
% Note that this routine returns quantal sensitivities.  You
 
30
% may want energy sensitivities.  In that case, use EnergyToQuanta to convert
 
31
%   T_energy = EnergyToQuanta(S,T_quantal')'
 
32
% and then renormalize.  (You call EnergyToQuanta because you're converting
 
33
% sensitivities, which go the opposite directoin from spectra.)
 
34
%
 
35
% The routine also returns two quantal sensitivity functions.  The first gives
 
36
% the probability that a photon will be absorbed.  The second is the probability
 
37
% that the photon will cause a photopigment isomerization.  It is the latter
 
38
% that is what you want to compute isomerization rates from retinal illuminance.
 
39
% See note at the end of function FillInPhotoreceptors for some information about
 
40
% convention.  In particular, this routine takes pre-retinal absorption into
 
41
% account in its computation of probability of absorptions and isomerizations,
 
42
% so that the relevant retinal illuminant is one computed without accounting for
 
43
% those factors.  This routine does not account for light attenuation due to
 
44
% the pupil, however.  The only use of pupil size here is becuase of its
 
45
% slight effect on lens density as accounted for in the CIE standard.
 
46
%
 
47
% In the passed params structure, you can either pass
 
48
% the lambdaMax values for the photopigment, in which
 
49
% case the absorbance is computed from the specified
 
50
% nomogram, or you can pass the absorbance values
 
51
% directly in T_xxx format.  A typical choice in this
 
52
% case would be 10.^T_lgo10coneabsorbance_ss for the
 
53
% Stockman-Sharpe/CIE estimates.
 
54
%
 
55
% The typical use of this function is to be called by 
 
56
% ComputeCIEConeFundamentals, which sets up the
 
57
% passed structures acording to the CIE standard. 
 
58
% This routine, however, could in principle be used
 
59
% with a wide variety of choices of the component pieces.
 
60
%
 
61
% The other place this function is used is in attempts to
 
62
% fit a set of cone fundamentals by doing parameter search
 
63
% over the pieces.  It is this second use that led to the
 
64
% parameters being split over two separate structures, one
 
65
% that is held static during the fit and the other which
 
66
% contains the parameters that could be searched over by a calling
 
67
% routine.  For examples, see:
 
68
%   FitConeFundamentalsWithNomogram, FitConeFundamentalsTest.
 
69
% Looking around today (8/10/13), I (DHB) don't see any examples where
 
70
% this routine is called directly.  Rather, it is a subfunction
 
71
% called by ComputeCIEConeFundamentals.  The search routines above
 
72
% use ComputeCIEConeFundamentals, and only search over lambdaMax
 
73
% values.  I think I wrote this with the thought that I might one
 
74
% day search over more parameters, but lost motivation to carry it
 
75
% throught.
 
76
%
 
77
% The computations done here are very similar to those done in
 
78
% routine FillInPhotoreceptors.  I (DHB) think that I forgot about
 
79
% what FillInPhotoreceptors could do when I wrote this, which has
 
80
% led to some redundancy. FillInPhotoreceptors returns a field
 
81
% called effectiveAbsorbtance, which are the actual quantal efficiencies
 
82
% (not normalized) referred to the cornea.  FillInPhotoceptors also
 
83
% computes a field isomerizationAbsorbtance, which takes the quantal
 
84
% efficiency of isomerizations (probability of an isomerization given
 
85
% an absorption into acount.  This routine does not do that.
 
86
%
 
87
% It would probably be clever to unify the two sets of routines a
 
88
% little more, but we may never get to it.  The routine ComputeCIEConeFundamentals
 
89
% does contain a check that this routine and what is returned by FillInPhotoreceptors
 
90
% agree with each other, for cases where the parameters match.
 
91
%
 
92
% See also: ComputeCIEConeFundamentals, CIEConeFundamentalsTest, FitConeFundamentalsWithNomogram,
 
93
%           FitConeFundamentalsTest, DefaultPhotoreceptors, FillInPhotoreceptors.
 
94
%
 
95
% 8/12/11  dhb  Starting to make this actually work.
 
96
% 8/14/11  dhb  Change name, expand comments.
 
97
% 8/10/13  dhb  Expand comments.  Return unscaled quantal efficiencies too.
 
98
 
 
99
% Handle bad value
 
100
index = find(params.axialDensity <= 0.0001);
 
101
if (~isempty(index))
 
102
    params.axialDensity(index) = 0.0001;
 
103
end
 
104
 
 
105
% Handle optional values
 
106
if (~isfield(params,'extraLens'))
 
107
    params.extraLens = 0;
 
108
end
 
109
if (~isfield(params,'extraMac'))
 
110
    params.extraMac = 0;
 
111
end
 
112
 
 
113
% Prereceptor transmittance.  Sometimes adjustments of peak density
 
114
% recommended by various standards push the density less than
 
115
% zero at some wavelengths, so we need to make sure we never have
 
116
% transmittance greater than 1.
 
117
lens = 10.^-(-log10(staticParams.lensTransmittance)+params.extraLens);
 
118
lens(lens > 1) = 1;
 
119
mac = 10.^-(-log10(staticParams.macularTransmittance)+params.extraMac);
 
120
mac(mac > 1) = 1;
 
121
 
 
122
% Compute nomogram if absorbance wasn't passed directly.  We detect
 
123
% a direct pass by the existance of params.absorbance.
 
124
if (isfield(params,'absorbance'))
 
125
    absorbance = params.absorbance;
 
126
else
 
127
    absorbance = PhotopigmentNomogram(staticParams.S,params.lambdaMax,staticParams.whichNomogram);
 
128
end
 
129
 
 
130
% Compute absorbtance
 
131
%
 
132
% Handle special case where we deal with ser/ala polymorphism for L cone
 
133
if (size(absorbance,1) == 4)
 
134
    absorbtance = AbsorbanceToAbsorbtance(absorbance,staticParams.S,...
 
135
        [params.axialDensity(1) ; params.axialDensity(1) ; ...
 
136
        params.axialDensity(2) ; params.axialDensity(3)]);
 
137
elseif (size(absorbance,1) == 3)
 
138
    absorbtance = AbsorbanceToAbsorbtance(absorbance,staticParams.S,...
 
139
        [params.axialDensity(1) ; params.axialDensity(2) ; ...
 
140
        params.axialDensity(3)]);
 
141
elseif (size(absorbance,1) == 1 && params.DORODS)
 
142
    absorbtance = AbsorbanceToAbsorbtance(absorbance,staticParams.S,...
 
143
        params.axialDensity(1));
 
144
else
 
145
    error('Unexpected number of photopigment lambda max values passed');
 
146
end
 
147
 
 
148
%% Put together pre-receptor and receptor parts
 
149
for i = 1:size(absorbtance,1)
 
150
    absorbtance(i,:) = absorbtance(i,:) .* lens .* mac;
 
151
end
 
152
 
 
153
%% Put it into the right form
 
154
if (size(absorbtance,1) == 4)
 
155
    T_quantalAbsorptions = zeros(3,staticParams.S(3));
 
156
    T_quantalAbsorptions(1,:) = staticParams.LserWeight*absorbtance(1,:) + ...
 
157
        (1-staticParams.LserWeight)*absorbtance(2,:);
 
158
    T_quantalAbsorptions(2,:) = absorbtance(3,:);
 
159
    T_quantalAbsorptions(3,:) = absorbtance(4,:);
 
160
elseif (size(absorbtance,1) == 3)
 
161
    T_quantalAbsorptions = zeros(3,staticParams.S(3));
 
162
    T_quantalAbsorptions(1,:) = absorbtance(1,:);
 
163
    T_quantalAbsorptions(2,:) = absorbtance(2,:);
 
164
    T_quantalAbsorptions(3,:) = absorbtance(3,:);
 
165
elseif (size(absorbtance,1) == 1 && params.DORODS)
 
166
    T_quantalAbsorptions = zeros(1,staticParams.S(3));
 
167
    T_quantalAbsorptions(1,:) = absorbtance(1,:);
 
168
else
 
169
    error('Unexpected number of photopigment lambda max values passed');
 
170
end
 
171
 
 
172
%% Normalize to max of one for each receptor, and also compute isomerization quantal efficiency.
 
173
for i = 1:size(T_quantalAbsorptions,1)
 
174
    T_quantalIsomerizations = T_quantalAbsorptions*staticParams.quantalEfficiency(i);
 
175
    T_quantalAbsorptionsNormalized(i,:) = T_quantalAbsorptions(i,:)/max(T_quantalAbsorptions(i,:));
 
176
end
 
177
 
 
178