psychtoolbox-3-common 3.0.14.20170103+git6-g605ff5c.dfsg1-1build1 / usr / share / psychtoolbox-3 / PsychDemos / IsomerizationsInEyeDemo.m

% IsomerizationsInEyeDemo
%
% Shows how to compute photoreceptor isomerizations using toolbox
% routines.  These calculations are for the human eye,
% starting with a spectrum as measured by the PR-650
% in watts/sr-m^2-wlinterval, or with a relative spectrum
% and a photopic troland value.
%
% NOTE, DHB, 7/19/13. This demo routine and its associated data routines 
% (DefaultPhotoreceptors, FillInPhotoreceptors, PrintPhotoreceptors)
% should be better integrated with the more recent code that
% implements the CIE physiological cone fundamentals, and the
% whole set of stuff should be better documented.  See also
%    IsomerizationsInDishDemo
%    CIEConeFundamentalsTest
%    ComputeCIEConeFundamentals
%    ComputeRawConeFundamentals
%    DefaultPhotoreceptors
%    FillInPhotoreceptors
%    PrintPhotoreceptors
%    RetIrradianceToIsoRecSec
% In particular, there should be some default for the 
% photoreceptors structure that gives one the CIE cone
% fundamentals in all their parametric glory, plus additional
% parameters that yield real energy/quantal sensitivites so
% that the resulting coordinates are isomerization rates in
% real units.  I think that we're close to having that, but
% better documentation and tidying is needed.
%
% 07/08/03 dhb  Wrote starting from IsomerizationsInDishDemo.
% 07/11/03 dhb  Grab data through subroutines.  Get rid of integration time.
% 07/15/03 dhb  Take eye size from function.
% 08/14/11 dhb  Comment out saving of T_dogrec at end.  Want to be careful when and where
%               this is done, but the template may be useful someday.
% 03/20/12 dhb  Update cal file for PTB 3.
% 04/09/12 dhb  Add test of irradiance to troland conversion.
% 04/27/13 dhb  More extensive comments.
% 7/19/13  dhb  Print out photoreceptors structure using PrintPhotoreceptors.
%          dhb  Add monochromatic light option to the section that starts with trolands.
% 8/11/13  dhb  Add test of AborbtanceToAbsorbance.
%          dhb  Protect against case when absorbance is provided directly.
% 05/26/14 dhb  Dusted off.
% 6/10/14  npc, dhb  Modifications for accessing calibration data using a @CalStruct object.
% 7/7/14   dhb  Make calStruct object code conditional on the support routines existing on the path.

%% Clear
clear; close all;

%% Set photoreceptor properties.
%
% The photoreceptors structure gets filled with
% key parameters values (pupil size, eye length,
% pre-retinal absorbance, etc.)
%
% The routine DefaultPhotoreceptors is a high level
% call.  It fills in the 'source' fields and some
% values according to high-level descriptor (e.g.,
% ('CIE2Deg').  See help for that routine
% for available options.
%
% The routine FillInPhotoreceptors fetches the actual
% values for various fields, depending on the source.
%
% To get a feel for this, check what is in the photoreceptors
% structure after the first call, and then after the second.
whatCalc = 'CIE2Deg';
photoreceptors = DefaultPhotoreceptors(whatCalc);
photoreceptors.eyeLengthMM.source = 'LeGrand';
photoreceptors = FillInPhotoreceptors(photoreceptors);

%% Check AbsorptancetoAbsorbance
%
% Simple check that this routine does what we expect, since
% we never use it anywhere else and this seems like as good
% a place to test it as any.
%
% We omit the normalization, because sometimes the wavelength
% sampling we use leads to a maximimum initial absorbance that
% is not unity, and letting AbsorptanceToAbsorbance normalize
% causes disagreement.
testAbsorbance = photoreceptors.absorbance;
testAbsorptance = photoreceptors.absorptance;
checkAbsorbance = AbsorptanceToAbsorbance(testAbsorptance, photoreceptors.nomogram.S, photoreceptors.axialDensity.bleachedValue,false);
diffs = testAbsorbance-checkAbsorbance;
if (max(abs(diffs(:))) > 1e-7)
    error('Cannot properly invert absorbance/absorptance computations');
end

%% Define common wavelength sampling for this script.
% 
% S is [start delta nsamples] for the wavelengths in nm.
% This is standard PTB convention.
S = photoreceptors.nomogram.S;

%% XYZ color matching functions
load T_xyz1931
T_xyz = SplineCmf(S_xyz1931,683*T_xyz1931,S);
T_Y = T_xyz(2,:);
        
%% Get light spectrum.  You can choose various illustrative examples.
%
% Available options:
%  'fromTrolands'
%  'fromMonitorRadiance'
%  'fromUniformQuantalSpd'
whichInputType = 'fromTrolands';
switch (whichInputType)
    
    % Start with troland value and a relative spectrum
    case 'fromTrolands'
        fprintf('Computing from troland value and relative spectrum\n');
        % Give troland value and type
        %
        % Type may be 'Photopic', 'Scotopic', or 'JuddVos'
        trolands = 1;
        trolandType = 'Photopic';
        switch (trolandType)
            case 'Photopic'
                fprintf('Using photopic trolands\n');
            case 'Scotopic'
                fprintf('Using scotopic trolands\n');
            case 'JuddVos'
                fprintf('Using Judd-Vos luminosity function in troland calculations\n');
            otherwise
                fprintf('Unknown troland type specified');
        end
        
        %% Pupil.
        %
        % We do these computations for a fixed pupil size, ignoring the
        % pupilDiamter field of the photoreceptors structure.  That field
        % is set up to specify a source formula that estimates pupil diameter
        % from luminance.
        %
        % Since we are starting in trolands, the pupil size shouldn't actually
        % effect the calculations, except for finding the radiance that is
        % equivalent to the specified troland value.  
        % 
        % We remove the pupilDiameter.source field to make sure we aren't sending
        % mixed messages about how we want to handle pupil diameter.
        photoreceptors.pupilDiameter.value = 2;
        if (isfield(photoreceptors.pupilDiameter,'source'))
            photoreceptors.pupilDiameter = rmfield(photoreceptors.pupilDiameter,'source');
        end
        pupilAreaMm2 = pi*(photoreceptors.pupilDiameter.value/2)^2;
        
        % Specify relative spectrum to be used in
        % conversion to a full spectrum.
        %
        % Choices are:
        %  'Monochromatic'
        %  'XenonArc'
        %
        % If type is 'Monochromatic', must specify
        % wavelengthNm.
        spectrumType = 'Monochromatic';
        switch (spectrumType)
            case 'Monochromatic'
                monoWavelengthNm = 550;  
                wls = SToWls(S);
                monoWavelengthIndex = find(wls == monoWavelengthNm);
                if (isempty(monoWavelengthIndex))
                    error('No sample wavelength matches desired wavelength');
                end
                spd_fromTrolands = zeros(size(wls));
                spd_fromTrolands(monoWavelengthIndex) = 1;
                fprintf('Using monochromatic %d nm light as relative spectrum\n',monoWavelengthNm);
            case 'XenonArc'
                load spd_xenonArc;
                spd_fromTrolands = SplineSpd(S_xenonArc,spd_xenonArc,S);
                clear S_xenonArc spd_xenonArc
                fprintf('Using the spectrum of a xenon arc lamp as relative spectrum\n');
            otherwise
                error('Unknown spectrum type specified');
        end
        
        % Convert trolands back to spectral retinal irradiance.  This
        % depends on the pupil size and eye length specified.
		irradianceWattsPerUm2 = TrolandsToRetIrradiance(spd_fromTrolands,S,trolands, ...
			trolandType,photoreceptors.species,photoreceptors.eyeLengthMM.value);
        irradianceTrolandsCheck = RetIrradianceToTrolands(irradianceWattsPerUm2,S,trolandType, ...
            photoreceptors.species,photoreceptors.eyeLengthMM.value);
        trolandsCheck = sum(irradianceTrolandsCheck);
        fprintf('Input troland value is %0.1f, checked value is %0.1f\n',trolands,trolandsCheck);
		irradianceWattsPerUm2 = SplineSpd(S,irradianceWattsPerUm2,S);
	
		% Another way to do this calculation.  Pupil size should cancel out.  Should get
	    % same answer as above.  This has as a byproduct computing a stimulus radiance,
        % which is useful for some of the common printout below.
        luminanceCdM2FromTrolands = TrolandsToLum(trolands,pupilAreaMm2);
		radianceWattsPerM2Sr = LumToRadiance(spd_fromTrolands,S,luminanceCdM2FromTrolands,trolandType); 
        photopicLuminanceCdM2 = T_Y*radianceWattsPerM2Sr;
		irradianceWattsCheck = RadianceToRetIrradiance(radianceWattsPerM2Sr,S,pupilAreaMm2,photoreceptors.eyeLengthMM.value);
		figure(1); clf; hold on
		set(plot(SToWls(S),irradianceWattsPerUm2,'r'),'LineWidth',2);
		set(plot(SToWls(S),irradianceWattsCheck,'k'),'LineWidth',2);
		set(title('Check of trolands to irradiance calculation'),'FontSize',14);
		set(xlabel('Wavelength (mm)'),'FontSize',14);
		set(ylabel('Irradiance'),'FontSize',14);
        
        % For case if monochromatic light, can check retinal irradiance against
        % the direct formulae provided in W&S, 2cd edition, p. 105, eqs 2.4.4 
        % Note that these equations are missing a factor of
        % the wavelength in the numerator for the quantal conversions, as
        % pointed out by Makous in his 1997 JOSA paper.
        %
        % This only works for 'Photopic' and 'Scotopic' trolands.
        %
        % There must be an eye length implicit in this calculation. Using
        % the LeGrand model eye length gives good agreement between these
        % and our more general calculations done in the main section below.
        %
        % Makous (1997), JOSA A, 14, p. 2331 gives retinal illuminances of
        % 5.44 quanta/[um2-sec] for 1 scotopic troland, and 14.65 quanta/[um2-sec]
        % for 1 photopic troland, with the calulations specified for 510 nm.
        % The calculations here, done in several different ways, yield
        % 5.442 (scotopic, agrees) and 26.85 (photopic, does not agree).  But at 
        % 550 nm, the code here yields 14.64 quanta/[um2-sec], which seems close
        % enough to the provided 14.65 to make me think that Makous' value is
        % actually for a wavelength close to 550 nm.  That would be a more typical
        % wavelength at which to do photopic calculations, despite what the text
        % in the paper says.
        %
        % Note that there are also errors in Tables 3 and 4 of the Makous
        % paper, and that corrected values appear in Tables 3 and 2 of
        % Makous (2004), Scotopic vision, In The Visual Neurosciences,
        % Werner and Chalupa (eds).  What does not seem to be specified
        % in either place is the wavelengths used in the calculations of
        % the two tables.
        if (strcmp(spectrumType,'Monochromatic'))
            switch (trolandType)
                case 'Photopic'
                    load T_xyz1931;
                    T_vLambda = SplineCmf(S_xyz1931,T_xyz1931(2,:),S);
                    clear T_xyz1931 S_xyz1931
                    magicFactorE = 5.261e-12;
                    magicFactorQ = 2.649e13;                  
                case 'Scotopic'
                    load T_rods;
                    T_vLambda = SplineCmf(S_rods,T_rods,S);
                    magicFactorE = 2.114e-12;
                    magicFactorQ = 1.064e13;
            end
        end
        if (strcmp(spectrumType,'Monochromatic'))
            switch (trolandType)
                case {'Photopic','Scotopic'}
                    irradianceDirectWattsPerMm2Check = magicFactorE*trolands/T_vLambda(monoWavelengthIndex);
                    irradianceDirectQuantaPerMm2SecCheck = 1e-9*monoWavelengthNm*magicFactorQ*trolands/T_vLambda(monoWavelengthIndex);
                    fprintf('Retinal irradiance in area units computed from from trolands via Wyszecki and Stiles formulae\n');
                    fprintf('\t%0.4g Watts/mm2\n\t%0.4g quanta/[mm2-sec]\n',sum(irradianceDirectWattsPerMm2Check),sum(irradianceDirectQuantaPerMm2SecCheck));
            end
        end
        
	% Start with radiance measurements, which we just
	% pull out of the Toolbox's default calibration file.
	case 'fromMonitorRadiance'
		% Load light radiance.  We'll use a monitor white.
		% The original units are watts/sr-m^2-wlinterval.
		cal = LoadCalFile('PTB3TestCal');
        
        % If we're running using BrainardLabToolbox and support CalStruct objects, use those.
        % But if we don't know about those, use old style direct access of fields in cal
        % structure.  The test is the existance of the routine ObjectToHandleCalOrCalStruct().
        if (exist('ObjectToHandleCalOrCalStruct','file'))
            [calStructOBJ, inputArgIsACalStructOBJ] = ObjectToHandleCalOrCalStruct(cal); clear 'cal';
            P_device = calStructOBJ.get('P_device');
            S_device = calStructOBJ.get('S');
            radianceWattsPerM2Sr = SplineSpd(S_device,sum(P_device,2),S);
        else
            P_device = cal.P_device;
            S_device = cal.S_device;
        end
        radianceWattsPerM2Sr = SplineSpd(S_device,sum(P_device,2),S);
        clear S_device P_device

		% Find pupil area, needed to get retinal irradiance.  We compute
		% pupil area based on the luminance of stimulus according to the
        % algorithm specified in the photoreceptors structure.
		theXYZ = T_xyz*radianceWattsPerM2Sr; theLuminance = theXYZ(2);
        if (isfield(photoreceptors.pupilDiameter,'value'))
            pupilAreaMm2 = pi*(photoreceptors.pupilDiameter.value/2).^2;
        else
            [nil,pupilAreaMm2] = PupilDiameterFromLum(theLuminance,photoreceptors.pupilDiameter.source);
        end
        photopicLuminanceCdM2 = T_Y*radianceWattsPerM2Sr;
		
		% Convert radiance of source to retinal irradiance and convert to quantal units.
		irradianceWattsPerUm2 = RadianceToRetIrradiance(radianceWattsPerM2Sr,S, ...
			pupilAreaMm2,photoreceptors.eyeLengthMM.value);

	% This light as well as some parameter tweaking are here to match a parameterization that Brian Wandell supplied
	% to match what his code to do these computations produces.  Note also
	% the mucking with the photoreceptors structure.  Wandell estimates
	% L, M, S isomerizations/cone-sec of 16.5, 12.68, 2.27.  These are very close to the numbers
    % we get here.
	case 'fromUniformQuantalSpd',
		% Load corneal cone sensitivities in energy units, convert to quantal sensitivities
		% and set specified peak absorptance.
		%
		% Note that overwriting the isomerizationAbsorptance in the photoreceptors structure
		% makes the isomerization computation work, but not the absorptions calculation, which
		% will be done with what was produced by FillInPhotoreceptors called above.  This is
		% not a recommended compute path for the toolbox code, but is done here to match Wandell's
		% parameterization.
		load T_cones_ss2; T_cones = T_cones_ss2; S_cones = S_cones_ss2;
		% load T_cones_ss10; T_cones = T_cones_ss10; S_cones = S_cones_ss10;
		% load T_cones_smj; T_cones = T_cones_smj; S_cones = S_cones_smj;
		% load T_cones_sp; T_cones = T_cones_sp; S_cones = S_cones_sp;
        
		peakIsomerizationEfficiency = [0.27 0.23 0.07]';
		T_cones = SplineCmf(S_cones,QuantaToEnergy(S_cones,T_cones')',S);
		T_cones(1,:) = T_cones(1,:)/max(T_cones(1,:));
		T_cones(2,:) = T_cones(2,:)/max(T_cones(2,:));
		T_cones(3,:) = T_cones(3,:)/max(T_cones(3,:));
		T_cones = diag(peakIsomerizationEfficiency)*T_cones;
		photoreceptors.isomerizationAbsorptance = T_cones;

		% Create a spectrally uniform spd (in quantal units), and convert
		% to energy units.
		uniformSpd = QuantaToEnergy(S,ones(S(3),1));

		% Normalize to radiance corresponding to 1 cd/m2.
		normConst = T_Y*uniformSpd;
		radianceWattsPerM2Sr = uniformSpd/normConst;
        photopicLuminanceCdM2 = T_Y*radianceWattsPerM2Sr;

		% Set pupil diameter for 1 mm2 pupil area, photoreceptor diameter for 4 mm2 collecting
		% area.  Set eye length to 17 mm.
		photoreceptors.pupilDiameter.value = 2*sqrt(1/pi);
		pupilAreaMm2 = pi*(photoreceptors.pupilDiameter.value/2)^2;
		photoreceptors.ISdiameter.value = [2*sqrt(4/pi) 2*sqrt(4/pi) 2*sqrt(4/pi)]';
		photoreceptors.eyeLengthMM.value = 17;
		irradianceWattsPerUm2 = RadianceToRetIrradiance(radianceWattsPerM2Sr,S,pupilAreaMm2,photoreceptors.eyeLengthMM.value );
end

%% Print out a whole bunch of quantities that are equivalent to the radiance, given
% other eye parameters.
radianceWattsPerCm2Sr = (10.^-4)*radianceWattsPerM2Sr;
radianceQuantaPerCm2SrSec = EnergyToQuanta(S,radianceWattsPerCm2Sr);
degPerMm = RetinalMMToDegrees(1,photoreceptors.eyeLengthMM.value);
irradianceWattsPerUm2Check = RadianceToRetIrradiance(radianceWattsPerM2Sr,S,pupilAreaMm2,photoreceptors.eyeLengthMM.value);
if (any(abs(irradianceWattsPerUm2 - irradianceWattsPerUm2Check) > 1e-10))
    error('Back computation of retinal irradiance from radiance does not check');
end
irradianceScotTrolands = RetIrradianceToTrolands(irradianceWattsPerUm2, S, 'Scotopic', [], num2str(photoreceptors.eyeLengthMM.value));
irradiancePhotTrolands = RetIrradianceToTrolands(irradianceWattsPerUm2, S, 'Photopic', [], num2str(photoreceptors.eyeLengthMM.value));
irradianceQuantaPerUm2Sec = EnergyToQuanta(S,irradianceWattsPerUm2);
irradianceWattsPerCm2 = (10.^8)*irradianceWattsPerUm2;
irradianceQuantaPerCm2Sec = (10.^8)*irradianceQuantaPerUm2Sec;
irradianceQuantaPerMm2Sec = (10.^-2)*irradianceQuantaPerCm2Sec;
irradianceQuantaPerUm2Sec = (10.^-6)*irradianceQuantaPerMm2Sec;
irradianceQuantaPerDeg2Sec = (degPerMm^2)*irradianceQuantaPerMm2Sec;

% Print out photoreceptor stucture information
fprintf('\n');
PrintPhotoreceptors(photoreceptors);
fprintf('\n');

% Radiometric iformation
fprintf('Luminance %0.3f cd/m2\n',photopicLuminanceCdM2);
fprintf('Stimulus retinal irradiance %0.4g (%0.1f log10) watts/cm2\n',sum(irradianceWattsPerCm2),log10(sum(irradianceWattsPerCm2)));
fprintf('Stimulus retinal irradiance %0.4g (%0.1f log10) watts/mm2\n',1e-2*sum(irradianceWattsPerCm2),log10(sum(irradianceWattsPerCm2)));
fprintf('Stimulus retinal irradiance %0.4g (%0.1f log10) quanta/[cm2-sec]\n',sum(irradianceQuantaPerCm2Sec),log10(sum(irradianceQuantaPerCm2Sec)));
fprintf('Stimulus retinal irradiance %0.4g (%0.1f log10) quanta/[mm2-sec]\n',sum(irradianceQuantaPerMm2Sec),log10(sum(irradianceQuantaPerMm2Sec)));
fprintf('Stimulus retinal irradiance %0.4g (%0.1f log10) quanta/[um2-sec]\n',sum(irradianceQuantaPerUm2Sec),log10(sum(irradianceQuantaPerUm2Sec)));
fprintf('Stimulus retinal irradiance %0.4g (%0.1f log10) quanta/[deg2-sec]\n',sum(irradianceQuantaPerDeg2Sec),log10(sum(irradianceQuantaPerDeg2Sec)));
fprintf('\n');
        
%% Get retinal irradiance in quanta/[sec-um^2-wlinterval]
irradianceQuanta = EnergyToQuanta(S,irradianceWattsPerUm2);
figure(2); clf; set(gcf,'Position',[100 400 700 300]);
subplot(1,2,1); hold on
set(plot(SToWls(S),irradianceQuanta,'r'),'LineWidth',2);
set(title('Light Spectrum'),'FontSize',14);
set(xlabel('Wavelength (nm)'),'FontSize',12);
set(ylabel('Quanta/sec-um^2-wlinterval'),'FontSize',12);

%% Do the work in toolbox function
[isoPerConeSec,absPerConeSec,photoreceptors] = ...
	RetIrradianceToIsoRecSec(irradianceWattsPerUm2,S,photoreceptors);

% Make a plot showing the effective photoreceptor sensitivities in quantal
% units, expressed as probability of isomerization.
subplot(1,2,2); hold on
set(plot(SToWls(S),photoreceptors.isomerizationAbsorptance(1,:),'r'),'LineWidth',2);
set(plot(SToWls(S),photoreceptors.isomerizationAbsorptance(2,:),'g'),'LineWidth',2);
set(plot(SToWls(S),photoreceptors.isomerizationAbsorptance(3,:),'b'),'LineWidth',2);
set(title('Isomerization Absorptance'),'FontSize',14);
set(xlabel('Wavelength (nm)'),'FontSize',12);
set(ylabel('Probability'),'FontSize',12);
axis([300 800 0 1]);

% Print out a table summarizing the calculation.
fprintf('***********************************************\n');
fprintf('Isomerization calculations for living human retina\n');
fprintf('\n');
fprintf('Photoreceptor Type             |\t       L\t       M\t     S\n');
fprintf('______________________________________________________________________________________\n');
fprintf('\n');
if (isfield(photoreceptors.nomogram,'lambdaMax'))
    fprintf('Lambda max                     |\t%8.1f\t%8.1f\t%8.1f\t nm\n',photoreceptors.nomogram.lambdaMax);
end
if (isfield(photoreceptors,'OSlength') && ~isempty(photoreceptors.OSlength.value))
    fprintf('Outer Segment Length           |\t%8.1f\t%8.1f\t%8.1f\t um\n',photoreceptors.OSlength.value);
end
if (isfield(photoreceptors,'OSdiameter') && ~isempty(photoreceptors.OSdiameter.value))
    fprintf('Outer Segment Diameter         |\t%8.1f\t%8.1f\t%8.1f\t um\n',photoreceptors.OSdiameter.value);
end
fprintf('Inner Segment Diameter         |\t%8.1f\t%8.1f\t%8.1f\t um\n',photoreceptors.ISdiameter.value);
fprintf('\n');
if (isfield(photoreceptors,'specificDensity') && ~isempty(photoreceptors.specificDensity.value))
    fprintf('Axial Specific Density         |\t%8.3f\t%8.3f\t%8.3f\t /um\n',photoreceptors.specificDensity.value);
end
fprintf('Axial Optical Density          |\t%8.3f\t%8.3f\t%8.3f\n',photoreceptors.axialDensity.value);
fprintf('Bleached Axial Optical Density |\t%8.3f\t%8.3f\t%8.3f\n',photoreceptors.axialDensity.bleachedValue);
fprintf('Peak isomerization prob.       |\t%8.3f\t%8.3f\t%8.3f\n',max(photoreceptors.isomerizationAbsorptance,[],2));
fprintf('______________________________________________________________________________________\n');
fprintf('\n');
fprintf('Absorption Rate                |\t%4.2e\t%4.2e\t%4.2e\t quanta/photoreceptor-sec\n',...
	absPerConeSec);
fprintf('Isomerization Efficiency       |\t%8.3f\t%8.3f\t%8.3f\n',...
	photoreceptors.quantalEfficiency.value);
fprintf('Isomerization Rate             |\t%4.2e\t%4.2e\t%4.2e\t iso/photoreceptor-sec\n',...
	 isoPerConeSec);
fprintf('In log10 units                 |\t%8.2f\t%8.2f\t%8.2f\t log10(iso)/photoreceptor-sec\n',...
	 log10(isoPerConeSec));
fprintf('______________________________________________________________________________________\n');

% Allow dumping out of photoreceptor sensitivities into a file for use elsewhere.  We want energy sensitivities
% for this purpose
% switch (whatCalc)
%     % Dog receptors (L, S, rod) in energy units, normalized to max of 1.
%     case 'LivingDog'
%         T_dogrec = EnergyToQuanta(S,photoreceptors.isomerizationAbsorptance')';
%         for i = 1:3
%             T_dogrec(i,:) = T_dogrec(i,:)/max(T_dogrec(i,:));
%         end
%         S_dogrec = S;
%         save T_dogrec T_dogrec S_dogrec
%     otherwise
% end