function [T_quantalAbsorptionsNormalized,T_quantalAbsorptions,T_quantalIsomerizations,adjIndDiffParams] = ComputeRawConeFundamentals(params,staticParams)
% [T_quantalAbsorptionsNormalized,T_quantalAbsorptions,T_quantalIsomerizations,adjIndDiffParams] = ComputeRawConeFundamentals(params,staticParams)
%
% Function to compute normalized cone quantal sensitivities from underlying
% pieces and parameters.
%
% Note that 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 types of quantal sensitivity functions. The
% first gives the probability that a photon will be absorbed. These are
% returned in variable T_quantalAbsorptionsNormalized and
% T_quantalAbsorptions, with the first being normalized. The second is the
% probability that the photon will cause a photopigment isomerization. This
% is returned in T_quantalIsomerizations.
%
% It is T_quantalIsomerizations that you want to use to compute
% isomerization rates from retinal illuminance. See note at the end of
% function FillInPhotoreceptors for some information about conventions. In
% particular, this routine takes pre-retinal absorption into account in its
% computation of probability of absorptions and isomerizations, so that the
% relevant retinal illuminance 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. Nor does it
% account for the collecting area of a photoreceptor, for cones the inner
% segment diameter.
%
% In the passed params structure, you can either pass the lambdaMax values
% for the photopigment, in which case the absorbance is computed from the
% specified nomogram, or you can pass the absorbance values directly in
% T_xxx format. A typical choice in this case would be
% 10.^T_log10coneabsorbance_ss for the Stockman-Sharpe/CIE estimates.
%
% The typical use of this function is to be called by
% ComputeCIEConeFundamentals, which sets up the passed structures acording
% to the CIE standard. This routine, however, could in principle be used
% with a wide variety of choices of the component pieces.
%
% The other place this function is used is in attempts to fit a set of cone
% fundamentals by doing parameter search over the pieces. It is this
% second use that led to the parameters being split over two separate
% structures, one that is held static during the fit and the other which
% contains the parameters that could be searched over by a calling routine.
% For examples, see:
% FitConeFundamentalsWithNomogram, FitConeFundamentalsTest.
% Looking around today (8/10/13), I (DHB) don't see any examples where this
% routine is called directly. Rather, it is a subfunction called by
% ComputeCIEConeFundamentals. The search routines above use
% ComputeCIEConeFundamentals, and only search over lambdaMax values. I
% think I wrote this with the thought that I might one day search over more
% parameters, but lost motivation to carry it throught.
%
% The computations done here are very similar to those done in routine
% FillInPhotoreceptors. I (DHB) think that I forgot about what
% FillInPhotoreceptors could do when I wrote this, which has led to some
% redundancy. FillInPhotoreceptors returns a field called
% effectiveAbsorptance, which are the actual quantal efficiencies (not
% normalized) referred to the cornea. FillInPhotoceptors also computes a
% field isomerizationAbsorptance, which takes the quantal efficiency of
% isomerizations (probability of an isomerization given an absorption into
% acount.
%
% It would probably be clever to unify the two sets of routines a little
% more, but we may never get to it. The routine ComputeCIEConeFundamentals
% does contain a check that this routine and what is returned by
% FillInPhotoreceptors agree with each other, for cases where the
% parameters match.
%
% See ComputeCIEConeFundamentals for the breakdown of how the Asano et al.
% (2016) individual differences model is specified in params.indDiffParams.
%
% See ComputeCIEConeFundamentals for documentation of the adjIndDiffParams
% output argument.
%
% See also: ComputeCIEConeFundamentals, CIEConeFundamentalsTest,
% FitConeFundamentalsWithNomogram,
% FitConeFundamentalsTest, DefaultPhotoreceptors,
% FillInPhotoreceptors.
%
% 8/12/11 dhb Starting to make this actually work.
% 8/14/11 dhb Change name, expand comments.
% 8/10/13 dhb Expand comments. Return unscaled quantal efficiencies too.
% 2/26/16 dhb, ms Add in Asano et al. (2016) individual observer adjustments
% 3/30/17 ms Added output argument returning adjusted ind differences
% 6/4/18 ms Included absorbance spectrum into adjusted ind diff output arg
% Handle bad value
index = find(params.axialDensity <= 0.0001);
if (~isempty(index))
params.axialDensity(index) = 0.0001;
end
% Figure out how many receptor classes we're handling
if (isfield(params,'absorbance'))
nReceptorTypes = size(params.absorbance,1);
else
nReceptorTypes = length(params.lambdaMax);
end
% Fill in null individual differences parameters if they are not passed
if isempty(params.indDiffParams)
params.indDiffParams.lambdaMaxShift = zeros(nReceptorTypes,1);
params.indDiffParams.shiftType = 'linear';
params.indDiffParams.dlens = 0;
params.indDiffParams.dmac = 0;
params.indDiffParams.dphotopigment = zeros(nReceptorTypes,1);
end
% Handle optional values for lens and macular pigment density
% in the parameters structure. There are two mutually exclusive
% ways to do this. For historical reasons, we can pass an additive
% density adjustment. More recently we implemented the Asano et al. (2016)
% parameterization. We don't allow both to happen at once.
%
% The logic here is a little hairy, because the way that we used to
% adjust lens and mac density was additive, but Asano et al. (2016) do
% it in a multiplicative fashion, so we need a flag to keep track of what
% we're going to do with the numbers down below.
if (~isfield(params,'extraLens'))
params.extraLens = 0;
end
if (~isfield(params,'extraMac'))
params.extraMac = 0;
end
if (params.extraLens ~= 0 && params.indDiffParams.dlens ~= 0)
error('Cannot specify lens density adjustment two ways');
end
if (params.extraMac ~= 0 && params.indDiffParams.dmac ~= 0)
error('Cannot specify macular pigment density adjustment two ways');
end
OLDLENSWAY = true;
if (params.extraLens == 0)
OLDLENSWAY = false;
params.extraLens = params.indDiffParams.dlens/100;
end
OLDMACWAY = true;
if (params.extraMac == 0)
OLDMACWAY = false;
params.extraMac = params.indDiffParams.dmac/100;
end
% Prereceptor transmittance. Check that passed parameters are not so weird
% as to lead to transmittances greater than 1, and throw error if so.
if (OLDLENSWAY)
fprintf('Using old way of adjusting lens density. Consider switching to newer implementation via the params.indDiffParams field\n');
lens = 10.^-(-log10(staticParams.lensTransmittance)+params.extraLens);
else
lens = 10.^-(-log10(staticParams.lensTransmittance) * (1 + params.extraLens));
end
if (any(lens > 1))
error('You have passed parameters that make lens transmittance greater than 1');
end
%lens(lens > 1) = 1;
adjIndDiffParams.lens = lens;
if (OLDMACWAY)
fprintf('Using old way of adjusting macular pigment density. Consider switching to newer implementation via the params.indDiffParams field\n');
mac = 10.^-(-log10(staticParams.macularTransmittance)+params.extraMac);
else
mac = 10.^-(-log10(staticParams.macularTransmittance) * (1 + params.extraMac));
end
if (any(mac > 1))
error('You have passed parameters that make macular pigment transmittance greater than 1');
end
adjIndDiffParams.mac = mac;
%mac(mac > 1) = 1;
% Compute nomogram if absorbance wasn't passed directly. We detect
% a direct pass by the existance of params.absorbance.
if (isfield(params,'absorbance'))
absorbance = params.absorbance;
else
absorbance = PhotopigmentNomogram(staticParams.S,params.lambdaMax,staticParams.whichNomogram);
end
% Shift absorbance, if desired
if (~isempty(params.indDiffParams.lambdaMaxShift))
if (length(params.indDiffParams.lambdaMaxShift) ~= size(absorbance,1))
error('Length of passed lambdaMaxShift does not match number of absorbances available to shift');
end
absorbance = ShiftPhotopigmentAbsorbance(staticParams.S,absorbance,params.indDiffParams.lambdaMaxShift,params.indDiffParams.shiftType);
end
adjIndDiffParams.absorbance = absorbance;
% Compute absorptance
%
% Handle special case where we deal with ser/ala polymorphism for L cone
%
% We've put in the Asano et al. (2016) multiplicative adjustment. Since we
% weren't adjusting photopigment density previously (except for
% self-screening), we don't have to deal with two threads here.
%
% Note that density can also get adjusted according to light level to
% account for photopigment bleaching. That happens in
% FillInPhotoreceptors, which would typically be called before we get here.
% We think this is OK, because both the Asano et al. and the fraction
% bleaching adjustment are multiplicative adjustments of axial density, and
% multiplication commutes so it doesn't matter what order we do things in.
if (size(absorbance,1) == 4)
if (any(params.indDiffParams.dphotopigment ~= 0))
error('Cannot use Asano et al. individual cone model with our weird 4 cone calling mode');
end
absorptance = AbsorbanceToAbsorptance(absorbance,staticParams.S,...
[params.axialDensity(1) ; params.axialDensity(1) ; ...
params.axialDensity(2) ; params.axialDensity(3)]);
elseif (size(absorbance,1) == 3)
if (length(params.indDiffParams.dphotopigment) ~= 3)
error('Density adjustment parameter length not right for cones');
end
LDensity = params.axialDensity(1) * (1 + params.indDiffParams.dphotopigment(1)/100);
MDensity = params.axialDensity(2) * (1 + params.indDiffParams.dphotopigment(2)/100);
SDensity = params.axialDensity(3) * (1 + params.indDiffParams.dphotopigment(3)/100);
absorptance = AbsorbanceToAbsorptance(absorbance,staticParams.S,[LDensity ; MDensity ; SDensity]);
adjIndDiffParams.dphotopigment = [LDensity MDensity SDensity];
elseif (size(absorbance,1) == 1 && params.DORODS)
if (length(params.indDiffParams.dphotopigment) ~= 1)
error('Density adjustment parameter length not right for rods');
end
RodDensity = params.axialDensity(1) + params.indDiffParams.dphotopigment(1)/100;
absorptance = AbsorbanceToAbsorptance(absorbance,staticParams.S,RodDensity);
adjIndDiffParams.dphotopigment = RodDensity;
else
error('Unexpected number of photopigment lambda max values passed');
end
adjIndDiffParams.absorptance = absorptance;
%% Put together pre-receptor and receptor parts
for i = 1:size(absorptance,1)
absorptance(i,:) = absorptance(i,:) .* lens .* mac;
end
%% Put it into the right form
if (size(absorptance,1) == 4)
T_quantalAbsorptions = zeros(3,staticParams.S(3));
T_quantalAbsorptions(1,:) = staticParams.LserWeight*absorptance(1,:) + ...
(1-staticParams.LserWeight)*absorptance(2,:);
T_quantalAbsorptions(2,:) = absorptance(3,:);
T_quantalAbsorptions(3,:) = absorptance(4,:);
elseif (size(absorptance,1) == 3)
T_quantalAbsorptions = zeros(3,staticParams.S(3));
T_quantalAbsorptions(1,:) = absorptance(1,:);
T_quantalAbsorptions(2,:) = absorptance(2,:);
T_quantalAbsorptions(3,:) = absorptance(3,:);
elseif (size(absorptance,1) == 1 && params.DORODS)
T_quantalAbsorptions = zeros(1,staticParams.S(3));
T_quantalAbsorptions(1,:) = absorptance(1,:);
else
error('Unexpected number of photopigment lambda max values passed');
end
%% Normalize to max of one for each receptor, and also compute isomerization quantal efficiency.
for i = 1:size(T_quantalAbsorptions,1)
T_quantalIsomerizations = T_quantalAbsorptions*staticParams.quantalEfficiency(i);
T_quantalAbsorptionsNormalized(i,:) = T_quantalAbsorptions(i,:)/max(T_quantalAbsorptions(i,:));
end