OTF and pupil size

[npcottaris]

The following panel assemblies depict the human OTF as a function of pupil diameter from 1.0 mm to 6.0 mm. The code used to generate these panels is attached at the end.

(1) Perhaps the code is not setup correctly, but if it is, I think that the built-in OTF does not do a very good job for low frequencies. It needs to be computed at a much finer grid. Which brings me to my earlier issue (#26), of how to specify a fine support grid for the OTF.

(2) Also, assuming I did not do something stupid in the code, it appears that there is significant noise in the OTF for certain pupil sizes (e.g., diam=3.0 mm). Is this expected? I think these are also due to the coarse support grid.

(3) Finally, should't the amplitude of the human OTF be much higher than what is shown here ? I am comparing these curves to those shown in Fig. 2 of Watson's JOV(2013) 13(6) paper (http://www.journalofvision.org/content/13/6/18.short?related-urls=yes&legid=jov;13/6/18) It this attenuation because of another optical element, like a diffuser? The diffuser method is shown as 'skip'.

Nicolas

Here is the code:

function testOTF
    s_initISET;

    h = figure(1);
    set(h, 'Position', [100 100 650 760]);
    clf;

    % Pupil diameters to test
    pupilDiametersInMillimeters = (1:0.5:3.0);
    pupilDiametersInMillimeters = (3.5:0.5:6);

    for pupilSizeIndex = 1:numel(pupilDiametersInMillimeters)

        %% Retrieve examined pupil radius 
        pupilRadiusInMeters = pupilDiametersInMillimeters(pupilSizeIndex)/2.0/1000.0;

        %% Create human optics with given pupil radius
        optics = opticsCreate('human', pupilRadiusInMeters);

        %% Initialize optical image with above optics
        oi = oiCreate('human');
        oi = oiSet(oi, 'optics', optics);

        %% Compute optical image for given scene
        scene = sceneCreate('macbethd65');
        %% Make the scene angular size = 1 deg and place it at a distance = 1.0 m
        sceneAngularSizeInDeg = 0.5;
        sceneDistanceInMeters = 1.0;
        scene = sceneSet(scene,'wangular', sceneAngularSizeInDeg);
        scene = sceneSet(scene,'distance', sceneDistanceInMeters);

        %% Compute optical image
        oi = oiCompute(scene,oi); 

        %% Compute RGB rendition of optical image
        opticalRGBImage = oiGet(oi, 'rgb image');

        %% Retrieve the full OTF
        optics = oiGet(oi, 'optics');
        OTF    = fftshift(abs(opticsGet(optics,'otf data')));

        %% Retrieve the wavelength axis
        OTFwavelengths = opticsGet(optics,'otf wave');

        %% Retrieve the spatial frequency support. This is in cycles/micron
        OTFsupport = opticsGet(optics,'otf support', 'um');

        otf_sfXInCyclesPerMicron = OTFsupport{1};
        otf_sfYInCyclesPerMicron = OTFsupport{2};

        %% Convert to cycles/deg.
        % In human retina, 1 deg of visual angle is about 288 microns
        micronsPerDegee = 288;
        otf_sfX = otf_sfXInCyclesPerMicron * micronsPerDegee;
        otf_sfY = otf_sfYInCyclesPerMicron * micronsPerDegee;

        %% Get the 2D slice at 550 nm
        [~,waveIndex]   = min(abs(OTFwavelengths - 550));
        OTF550 = squeeze(OTF(:,:,waveIndex));

        %% Get a 2D slice through origin
        [~, sfIndex] = min(abs(otf_sfY - 0));
        OTFslice = squeeze(OTF550(sfIndex,:));

        %% Generate plots
        plotWidth = 0.3;
        plotHeight = 0.85/numel(pupilDiametersInMillimeters);
        margin     = 0.021;

        %% Plot the 2D OTF
        %subplot(numel(pupilDiametersInMillimeters),3, (pupilSizeIndex-1)*3+1);
        subplot('Position', [0.03 1+margin/2-pupilSizeIndex*(plotHeight+margin) plotWidth plotHeight]);
        imagesc(otf_sfX, otf_sfY, OTF550);
        axis 'image'
        axis 'xy'
        if (pupilSizeIndex == numel(pupilDiametersInMillimeters))
           xlabel('cycles/deg');
        else
           xlabel(''); 
        end
        set(gca, 'XLim', [-50 50], 'YLim', [-50 50]);
        colormap(gray(256));

        %% Plot the 1D OTF slice
        %subplot(numel(pupilDiametersInMillimeters),3, (pupilSizeIndex-1)*3+2);
        subplot('Position', [0.06+plotWidth 1+margin/2-pupilSizeIndex*(plotHeight+margin) plotWidth plotHeight]);
        indices = find(otf_sfX >= 0);
        OTFsfX = otf_sfX(indices)+0.1;
        plot(OTFsfX, OTFslice(indices), 'rs-');
        if (pupilSizeIndex == numel(pupilDiametersInMillimeters))
           xlabel('cycles/deg');
        else
           xlabel(''); 
        end
        text(0.12, 0.005, sprintf('Pupil:%2.1f mm', pupilDiametersInMillimeters(pupilSizeIndex)), 'FontSize', 12, 'FontWeight', 'bold');
        set(gca, 'XLim', [0.1 100], 'YLim', [0.001 1]);
        set(gca, 'XScale', 'log', 'YScale', 'log', 'XTick', [0.1 1 2 5 10 20 50 100], 'YTick', [0.001 0.002 0.005 0.01 0.02 0.05 0.10 0.20 0.50 1.0]);
        set(gca, 'XTickLabel', [0.1 1 2 5 10 20 50 100], 'YTickLabel',  [0.001 0.002 0.005 0.01 0.02 0.05 0.10 0.20 0.50 1.0]);
        box on;
        grid on;
        colormap(gray(256));

        %% Plot the RGB rendered optical image
        %subplot(numel(pupilDiametersInMillimeters),3, (pupilSizeIndex-1)*3+3);
        subplot('Position', [0.09+2*plotWidth 1+margin/2-pupilSizeIndex*(plotHeight+margin) plotWidth plotHeight]);
        imshow(opticalRGBImage);
        axis 'image'

        drawnow; 
    end
end

[wandell]

Hi Nicolas,

Great you are doing this. The code I have that analyzes the human OTF is in the script s_humanOptics.m

The values computed there don't fall off like the code you wrote. It might be best for you to compare what's in there with your script (which is much longer and has lots of moves in it that would be hard for me to track down). In the s_humanOptics.m script, at 550 nm with a pupil radius of 3 mm the OTF is down to about 0.15. Clearly, this is much higher than the values you are calculating.

That script could be improved, too, of course. It relies only on humanCore(), and the value of your script is that it goes all the way to the OI and Optics. We will try to stick with you on this, but don't give up hope. I think the core calculation is about right. There is a lot of complexity in tracking the variables and maintaining units.

[wandell]

This code also produces a reasonable OTF, using reasonably high level calls.

Brian

scene = sceneCreate('line d65'); scene = sceneSet(scene,'h fov',1); oi = oiCreate('human'); oi = oiCompute(oi,scene); vcAddObject(oi); oiWindow; vcNewGraphWin; plotOI(oi,'otf wavelength')

[npcottaris]

Hi Brian,

The problem was on my part. I was calling fftshift on the 3D OTF instead on the 2D slice at the desired wavelength. Stupid.

Now the OTF slices look much more reasonable. I will validate these against those appearing in Watson's paper.

Thanks, Nicolas

[wandell]

But what you are doing is really important. When you are done, we will all be feeling much more confident about the numbers. So far, it has been our lab and a few colleagues from industry who checked (mainly in cyc/mm space). When you are done, we will have a good validation procedure. Very cool.