function [T_quantalAbsorptionsNormalized,T_quantalAbsorptions,T_quantalIsomerizations,adjIndDiffParams,params,staticParams] = ...
ComputeCIEConeFundamentals(S,fieldSizeDegrees,ageInYears,pupilDiameterMM,lambdaMax,whichNomogram,LserWeight, ...
DORODS,rodAxialDensity,fractionPigmentBleached,indDiffParams)
% [T_quantalAbsorptionsNormalized,T_quantalAbsorptions,T_quantalIsomerizations,adjIndDiffParams,params,staticParams] = ...
% ComputeCIEConeFundamentals(S,fieldSizeDegrees,ageInYears,pupilDiameterMM,[lambdaMax],[whichNomogram],[LserWeight], ...
% [DORODS],[rodAxialDensity],[fractionPigmentBleached],[indDiffParams])
%
% Function to compute normalized cone quantal sensitivities
% from underlying pieces, as specified in CIE 170-1:2006.
%
% IMPORTANT: This routine returns quantal sensitivities. You
% may want energy sensitivities. In that case, use EnergyToQuanta to convert
% T_energy = EnergyToQuanta(S,T_quantal')'
% and then renormalize. (You call EnergyToQuanta because you're converting
% sensitivities, which go the opposite direction from spectra.)
% The routine also returns two quantal sensitivity functions. The first gives
% the probability that a photon will be absorbed. The second is the probability
% that the photon will cause a photopigment isomerization. It is the latter
% that is what you want to compute isomerization rates from retinal illuminance.
% See note at the end of function FillInPhotoreceptors for some information about
% convention. In particular, this routine takes pre-retinal absorption into
% account in its computation of probability of absorptions and isomerizations,
% so that the relevant retinal illuminant is one computed without accounting for
% those factors. This routine does not account for light attenuation due to
% the pupil, however. The only use of pupil size here is becuase of its
% slight effect on lens density as accounted for in the CIE standard.
%
% This standard allows customizing the fundamentals for
% field size, observer age, and pupil size in mm.
%
% To get the Stockman-Sharpe/CIE 2-deg fundamentals, use
% fieldSizeDegrees = 2;
% ageInYears = 32;
% pupilDiameterMM = 3;
% and don't pass the rest of the arguments.
%
% To get the Stockman-Sharpe/CIE 10-deg fundamentals, use
% fieldSizeDegrees = 10;
% ageInYears = 32;
% pupilDiameterMM = 3;
% and don't pass the rest of the arguments.
%
% Although this routine will compute something over any wavelength
% range, I'd (DHB) recommend not going lower than 390 or above about 780 without
% thinking hard about how various pieces were extrapolated out of the range
% that they are specified in the standard. Indeed, the lens optical
% density measurements only go down to 400 nm and these are extropolated
% to go below 400.
%
% This routine will compute from tabulated absorbance or absorbance based on a nomogram, where
% whichNomogram can be any source understood by the routine PhotopigmentNomogram. To obtain
% the nomogram behavior, pass a lambdaMax vector. You can then also optionally pass a nomogram
% source (default: StockmanSharpe). This option (using shifted nomograms) is not part of the
% CIE standard. See NOTE below for another way to handle individual differences
%
% The nominal values of lambdaMax to fit the CIE 2-degree fundamentals with the
% Stockman-Sharpe nomogram are 558.9, 530.3, and 420.7 nm for the LMS cones respectively.
% These in fact do a reasonable job of reconstructing the CIE 2-degree fundamentals, although
% there are small deviations from what you get if you simply read in the tabulated cone
% absorbances. Thus starting with these as nominal values and shifting is one way to
% produce fundamentals tailored to observers with different known photopigments.
%
% If you pass lambaMax and its length is 4, then first two values are treated as
% the peak wavelengths of the ser/ala variants of the L cone pigment, and these
% are then weighted according to LserWeight and (1-LserWeight). The default
% for LserWeight is 0.56. After travelling it for a distance to try to get better
% agreement between the nomogram based fundamentals and the tabulated fundamentals
% I (DHB) gave up and decided that using a single lambdaMax is as good as anything
% else I could come up with. If you are interested, see FitConeFundamentalsTest.
%
% NOTE 1: When we first implemented the CIE standard, adding this shifting feature
% seemed like a good idea to allow exploration of individual differences in photopigments.
% But, with 0 shift, none of the nomograms exactly reproduce the tabulated photopigment absorbance
% spectral sensitivities, and this is not so good. We are phasing out our
% use of this feature in favor of simply shifting the tabulated
% photopigment absorbances, and indeed in favor of adopting the method
% published by Asano, Fairchild, & Blonde (2016), PLOS One, doi: 10.1371/journal.pone.0145671
% to tailor the CIE fundamentals to individual observers. This is done by
% passing the argument indDiffParams, which is a structure as follows.
% indDiffParams.dlens - deviation in % from CIE computed peak lens density
% indDiffParams.dmac - deviation in % from CIE peak macular pigment density
% indDiffParams.dphotopigment - vector of deviations in % from CIE photopigment peak density.
% indDiffParams.lambdaMaxShift - vector of values (in nm) to shift lambda max of each photopigment absorbance by.
% indDiffParams.shiftType - 'linear' (default) or 'log'. 'linear' gets the Asano et al. behavior
%
% You also can shift the absorbances along a wavenumber axis after you have
% obtained them. To do this, pass argument lambdaMaxShift with the same
% number of entries as the number of absorbances that are used.
%
% The adjIndDiffParams output is a struct which is populated by ComputeRawConeFundamentals.
% It contains the actual parameter values for the parameters adjusted using the indDiffParams
% input. It contains the following fields:
% adjIndDiffParams.mac - the adjusted macular pigment transmittance as a function of wavelength
% as calculated in line 151 of ComputeRawConeFundamentals.
% adjIndDiffParams.lens - the adjusted lens transmittance as a function of wavelength as calculated
% in line 41 of ComputeRawConeFundamentals.
% adjIndDiffParams.dphotopigment - 3-vector of the adjusted photopigment axial density for
% L, M and S cones (in that order), as calculated in lines
% 200-202 of ComputeRawConeFundamentals; or rods, as calculated
% in line 216 of ComputeRawConeFundamentals if params.DORODS is true.
% adjIndDiffParams.absorbance - Photopigment absorbance as given in line 188 of ComputeRawConeFundamentals
% adjIndDiffParams.absorptance - Photopigment absorptance as given in line 230 of ComputeRawConeFundamentals
%
% For both adjIndDiffParams.mac and adjIndDiffParams.lens, the wavelength
% spacing is the same as in the S input variable of this function.
%
% The params and staticParams outputs are the argument strucutures that
% were passed to ComputeRawConeFundamentals by this routine to do the work.
% These can be useful if you'd like, say, to susequently use
% ComputeRawConeFundamentals to produce estimates for (e.g.) melanopsin or
% the rods, where you keep everything else as consistent as possible to
% what this routine does. Note that this is all a bit klugy for historical
% reasons, as there is redundancy between what you can/might do with
% adjIndDiffParams and with these two return outputs. In particular, these
% two return outputs would let you call ComputeRawConeFundamentals and get
% adjIndDiffParams directly from there.
%
% This function also has an option to compute rod spectral sensitivities, using
% the pre-retinal values that come from the CIE standard. Set DORODS to true on
% call. You then need to explicitly pass a single lambdaMax value. You can
% also pass an optional rodAxialDensity value. If you don't pass that, the
% routine uses the 'Alpern' estimate for 'Human'/'Rod' embodied in routine
% PhotopigmentAxialDensity. The default nomogram for the rod spectral
% absorbance is 'StockmanSharpe', but you can override with any of the
% others available in routine PhotopigmentNomogram. Use of this requires
% good choices for lambdaMax, rodAxialDensity, and the nomogram. We are
% working on identifying those values more precisely.
%
% Finally, you can adjust the returned spectral sensitivities to account for
% the possibility that some of the pigment in the cones is bleached. Pass
% a column vector with same length as number of spectral sensitivities beingt
% computed. You need to estimate the fraction elsewhere.
%
% Relevant to individual differences, S & S (2000) estimate the wavelength difference
% between the ser/ala variants to be be 2.7 nm (ser longer).
%
% NOTE 2. The CIE standard is specified for field sizes between 1 and 10
% degrees. Our code will extrapolate using the given formulae to larger
% field sizes without complaining. We think this is reasonable; see
% CIEConeFundamentalsFieldSizeTest and its header comments, but be aware
% that you have sailed into little charted territory if you do this.
%
% See also: ComputeRawConeFundamentals, CIEConeFundamentalsTest, CIEConeFundamentalsFieldSizeTest,
% FitConeFundamentalsTest, FitConeFundamentalsWithNomogram, StockmanSharpeNomogram,
% ComputePhotopigmentBleaching.
%
% 8/13/11 dhb Wrote it.
% 8/14/11 dhb Clean up a little.
% 12/16/12 dhb, ms Add rod option.
% 08/10/13 dhb Test for consistency between what's returned by FillInPhotoreceptors and
% what's returned by ComputeRawConeFundamentals.
% 05/24/14 dhb Add fractionPigmentBleached optional arg.
% 05/26/14 dhb Comment improvements.
% 02/08/16 dhb, ms Add lambdaMaxShift argument.
% ms Don't do two way check when lambdaMax is shifted.
% 02/24/16 dhb, ms Started to implement Asano et al. individual difference model
% 3/30/17 ms Added output argument returning adjusted ind differences
% 8/1/17 dhb, ms Added return of params and staticParams.
%% Are we doing rods rather than cones?
if (nargin < 8 || isempty(DORODS))
DORODS = 0;
end
%% Check whether we'll adjust axial density for bleaching
if (nargin < 10 || isempty(fractionPigmentBleached))
DOBLEACHING = 0;
else
DOBLEACHING = 1;
end
%% Check for passed lambdaMaxShift
if (nargin < 11 || isempty(indDiffParams))
params.indDiffParams = [];
else
params.indDiffParams = indDiffParams;
end
%% Get some basic parameters.
%
%
% We start with default CIE parameters in
% the photoreceptors structure, and then override
% as necessary.
% then override to match the CIE standard.
if (fieldSizeDegrees <= 4)
whatCalc = 'CIE2Deg';
else
whatCalc = 'CIE10Deg';
end
photoreceptors = DefaultPhotoreceptors(whatCalc);
%% Override default values so that FillInPhotoreceptors does
% our work for us. The CIE standard uses field size,
% age, and pupil diameter to computer other values.
% to compute other quantities.
photoreceptors.nomogram.S = S;
photoreceptors.fieldSizeDegrees = fieldSizeDegrees;
photoreceptors.pupilDiameter.value = pupilDiameterMM;
photoreceptors.ageInYears = ageInYears;
% Absorbance. Use tabulated CIE values (which are in the
% default CIE photoreceptors structure) unless a nomogram and
% lambdaMax values are passed.
SET_ABSORBANCE = false;
if (nargin > 4 && ~isempty(lambdaMax))
if (nargin < 6 || isempty(whichNomogram))
whichNomogram = 'StockmanSharpe';
end
photoreceptors = rmfield(photoreceptors,'absorbance');
photoreceptors.nomogram.source = whichNomogram;
photoreceptors.nomogram.lambdaMax = lambdaMax;
params.lambdaMax = lambdaMax;
staticParams.whichNomogram = whichNomogram;
else
% Absorbance is going to be specified directly. We get
% it after the call to FillInPhotoreceptors below,
% which will convert a file containing the absorbance into
% the needed data at the needed wavelength spacing.
SET_ABSORBANCE = true;
end
%% Are we doing the rods? In that case, a little more
% mucking is necessary.
if (DORODS)
if (isempty(lambdaMax) || length(lambdaMax) ~= 1)
error('When computing for rods, must specify exactly one lambda max');
end
photoreceptors.types = {'Rod'};
photoreceptors.nomogram.lambdaMax = lambdaMax;
photoreceptors.OSlength.source = 'None';
photoreceptors.specificDensity.source = 'None';
photoreceptors.axialDensity.source = 'Alpern';
params.DORODS = true;
end
%% Pigment bleaching
%
% Hope for the best with respect to dimensionality of what is passed.
% FillInPhotoreceptors will throw an error if the dimension isn't
% matched to that of the axialDensity value field.
if (DOBLEACHING)
photoreceptors.fractionPigmentBleached.value = fractionPigmentBleached;
end
%% Do the work. Note that to modify this code, you'll want a good
% understanding of the order of precedence enforced by FillInPhotoreceptors.
% This is non-trivial, although the concept is that if a quantity that
% can be computed is specified directly in the passed structure is
% actually specified, the speciefied value overrides what could be computed.
photoreceptors = FillInPhotoreceptors(photoreceptors);
if (SET_ABSORBANCE)
params.absorbance = photoreceptors.absorbance;
end
%% Set up for call into the low level routine that computes the CIE fundamentals.
staticParams.S = photoreceptors.nomogram.S;
staticParams.fieldSizeDegrees = photoreceptors.fieldSizeDegrees;
staticParams.ageInYears = photoreceptors.ageInYears;
staticParams.pupilDiameterMM = photoreceptors.pupilDiameter.value;
staticParams.lensTransmittance = photoreceptors.lensDensity.transmittance;
staticParams.macularTransmittance = photoreceptors.macularPigmentDensity.transmittance;
staticParams.quantalEfficiency = photoreceptors.quantalEfficiency.value;
CHECK_FOR_AGREEMENT = true;
if (nargin < 7 || isempty(LserWeight))
staticParams.LserWeight = 0.56;
else
staticParams.LserWeight = LserWeight;
end
if (DORODS && nargin >= 9 && ~isempty(rodAxialDensity))
params.axialDensity = rodAxialDensity;
CHECK_FOR_AGREEMENT = false;
else
params.axialDensity = photoreceptors.axialDensity.bleachedValue;
end
if (~isfield(params,'absorbance'))
if (length(params.lambdaMax) ~= 3 && length(params.lambdaMax) ~= 1)
CHECK_FOR_AGREEMENT = false;
end
end
% Shift in lambda max bookkeeping.
if isfield(params,'indDiffParams')
CHECK_FOR_AGREEMENT = false;
end
%% Drop into more general routine to compute
%
% See comment in ComputeRawConeFundamentals about the fact that
% we ought to unify this routine and what FillInPhotoreceptors does.
[T_quantalAbsorptionsNormalized,T_quantalAbsorptions,T_quantalIsomerizations,adjIndDiffParams] = ComputeRawConeFundamentals(params,staticParams);
%% A little reality check.
%
% The call to FillInPhotoreceptors also computes what here is called
% T_quantal. It is in the field effectiveAbsorptance. For cases where
% we aren't playing games with the parameters after the call to
% FillInPhotoreceptors, we can check for agreement.
if (CHECK_FOR_AGREEMENT)
diffs = abs(T_quantalAbsorptions(:)-photoreceptors.effectiveAbsorptance(:));
if (max(diffs(:)) > 1e-7)
error('Two ways of computing absorption quantal efficiency referred to the cornea DO NOT AGREE');
end
diffs = abs(T_quantalIsomerizations(:)-photoreceptors.isomerizationAbsorptance(:));
if (max(diffs(:)) > 1e-7)
error('Two ways of computing isomerization quantal efficiency referred to the cornea DO NOT AGREE');
end
end
end