# Class tutorial on optics and image formation (Psych 221)

This is a very general overview of image formation and a small tutorial on interacting with ISET structures.

Date: 01.02.96 (First drafted) Duration: 45 minutes

Matlab 5: Checked 01.06.98, BW 01.07.99: Added IjSpeert calls, BW 01.12.01 Added new graphs for ijspeert curves and Williams et al. data, BW 01.18.03 Started integration with vCamera-2.0, minor changes. 01.02.08 Verified with local files for single download - BW

## Contents

- Visual Angle
- The viewing geometry
- Geometric angle calculation
- The Westheimer linespread function
- Another way to calculate and plot the Westheimer function
- Convolution of the image and linespread
- Defocus in the frequency domain
- Cosinusoids at different spatial frequencies.
- The Fourier Transform and the linespread function.
- The pointspread and linespread
- How the linespread varies with wavelength
- Chromatic aberration in the frequency domain
- More modern Linespreads, Pointspreads, and MTFs
- Comparison of the MTFs with the Williams data
- What does the pointspread function look like in 2D?
- END

ieInit

## Visual Angle

% To think about the effect of an image on the eye, we must specify the % image in terms of degrees of visual angle. As an example for how to % compute spacing in terms of degrees of visual angle, consider a printer % whose dots are spaced (dots per inch) dpi = 600 % Suppose we read the paper at a viewing distance (inches) viewingDistance = 12

dpi = 600 viewingDistance = 12

## The viewing geometry

vcNewGraphWin; line([0 viewingDistance 0 0],[0 0 1 0]); axis equal set(gca,'xlim',[-2 20]), grid on xlabel('Viewing Distance (inch)') ylabel('Position between spots on paper (inch)')

## Geometric angle calculation

% In radians, the viewing angle, phi, satisfies % $tan(\phi) = (opposite/adjacent)$ deg2rad = 2*pi/360; rad2deg = 360/(2*pi) phi = atan(1/viewingDistance)*rad2deg % There are 600 dots per inch, so that each dot occupies % DegPerDot = phi/dpi % There are 60 min of visual angle per deg, % MinPerDot = 60*DegPerDot % and 60 sec of visual angle per min, % SecPerDot = 60*MinPerDot % As you will see later in the course, experiments have shown that people % can localize the position of a line to a spatial position of roughly 6 % sec of visual angle. Hence, at this viewing distance and with this many % dots per inch, the dot spacing is wider than the spacing that can be just % discriminated by the human eye.

rad2deg = 57.2958 phi = 4.7636 DegPerDot = 0.0079 MinPerDot = 0.4764 SecPerDot = 28.5819

## The Westheimer linespread function

% Westheimer calculated that the linespread function of the human % eye, specified in terms of minutes of arc and using a 3mm % pupil, should be approximated using the following formula % % LineSpread = 0.47*exp(-3.3 *(x.^2)) + 0.53*exp(-0.93*abs(x)); % % Suppose we wish to plot the function by defining the spatial % variable, x, in terms of seconds of arc, xSec = -300:1:300; xMin = xSec/60; ls = 0.47*exp(-3.3 *(xMin.^2)) + 0.53*exp(-0.93*abs(xMin)); ls = ls / sum(ls); vcNewGraphWin; plot(xSec,ls) set(gca,'xlim',[-240 240],'xtick',(-240:60:240)), grid on xlabel('Arc sec'), ylabel('Responsivity'), title('Westheimer Linespread') % From our previous calculation, we observed that the dots in a 600 dpi % printer, viewed at 12 inches, are spaced 28.5819 sec of visual angle % apart. At this distance, the linespread has fallen to about one-half of % its peak value. % Hence, if we could control the intensity and color of the printed dots -- % which we cannot do on conventional laser printers -- then at this viewing % distance we would be able to produce images that were very realistic in % their appearance.

## Another way to calculate and plot the Westheimer function

xSec = -300:300; % 600 sec, total = 10 min westheimerOTF = abs(fft(westheimerLSF(xSec))); % One cycle spans 10 min of arc, so freq=1 is 6 c/deg freq = [0:11]*6; vcNewGraphWin; semilogy(freq,westheimerOTF([1:12])); grid on; xlabel('Freq (cpd)'); ylabel('Relative contrast'); set(gca,'ylim',[-.1 1.1])

## Convolution of the image and linespread

% We can estimate the visual image created by an image printed at 600 dpi % using the following simple convolution calculation. Let's create an image % that spans 0.2 deg and has a dot every 30 sec. (i.e., roughly 28.58) dotSpacing = 30; secPerDeg = 60*60; x = 1:0.2*secPerDeg; im = zeros(1,length(x)); im(1:dotSpacing:length(im)) = ones(size(im(1:dotSpacing:length(im)))); % Here is an image showing the sampled line positions % (The lines represent a printer dot) vcNewGraphWin; imshow(im(ones(1,128),1:512)) title('Image of line stimulus'); % Each line in the physical image adds a unit linespread to the retinal % image. We can compute the retinal image by forming the convolution of % the image with the Westheimer linespread function. Remember: we sampled % the linespread once every sec of arc. So, we can simply convolve the % image and the linespread function now as: retIm = conv2(ls,im,'full'); vcNewGraphWin; plot(retIm),grid on title('The one-dimensional retinal image') xlabel('Sec of arc'), ylabel('Image intensity') % (If the colors in the image look funny, make sure to place your % cursor inside the window. This may change the local color % map). % While the original image varies from black to white, after blurring by % the eye's optics, there is only a small amount of residual variation in % the retinal image. Because of the blurring, the retinal image is much % more likely the image of a bar than it is the image of a set of % individual lines. % In fact, the dots placed on the page are not perfect line samples. Each % ink line has some width. So, a more realistic input image might be % created by blurring the stimulus and then convolving with the linespread. gKernel = fspecial('gaussian',[1,30],2); % This produces a little Gaussian window for filtering blurIm = conv2(im,gKernel,'full'); blurIm = ieScale(blurIm,1,32); retIm = conv2(blurIm,ls,'full'); vcNewGraphWin; imshow(blurIm(ones(1,128),1:512),[]); title('Image of line stimulus blurred by ink width'); % Notice the very small ripples left in the curve after taking % into account the blurring by the physical display and by the % eye. vcNewGraphWin; plot(retIm), axis square, grid on title('Retinal image of blurred lines') xlabel('Sec of arc'), ylabel('Intensity') % The question you might ask yourself now is this: will those % small ripples be detectable by the observers? How can we tell? % You might also ask what will happen when we view the page at 6 % inches, or at 24 inches. What if we increase the printer % resolution to 1200 dpi? What if we introduce some ability to % modulate the density of the ink and hence the light scattered % back to the eye?

## Defocus in the frequency domain

% First, make a new linespread function that is smaller and % easier to compute with. Have it extend over 1 deg (60 min) so % the Fourier Transform is easier to interpret xMin = -30:1:29; ls = 0.47*exp(-3.3 *(xMin.^2)) + 0.53*exp(-0.93*abs(xMin)); ls = ls / sum(ls); vcNewGraphWin; plot(xMin,ls), grid on xlabel('Min of arc'),ylabel('Linespread value') % Now, consider the retinal image that is formed by some simple % harmonic functions. Here is a sinusoid that varies at 1 cycle % per degree of visual angle. nSamples = length(xMin); f = 1; harmonic = cos(2*pi*f*xMin/nSamples); vcNewGraphWin; plot(xMin,harmonic), grid on title('Sampled cosinusoid') xlabel('Arc sec'), ylabel('Intensity')

## Cosinusoids at different spatial frequencies.

% Notice that the amplitude of the cosinusoid falls off as the spatial % frequency increases. We will store the amplitude of the cosinusoid in % the variable "peak". vcNewGraphWin; freq =[1 5 10 15]; peak = zeros(1,length(freq)); for i = 1:length(freq) harmonic = cos(2*pi*freq(i)*xMin/nSamples); retIm = convolvecirc(harmonic,ls); subplot(2,2,i) plot(retIm), grid on, set(gca,'ylim',[-1 1],'xlim',[0 64]); xlabel('Arc sec') peak(i) = max(retIm(:)) end % We can plot the amplitude of the retinal cosinusoid, and its % amplitude decreases with the input frequency. I also use the % fact that at f = 0 the amplitude = 1 (because the area under % the linespread is 1). vcNewGraphWin; plot([0 freq],[1 peak],'-') set(gca,'ylim',[0 1]) xlabel('Spatial freq (cpd)'), ylabel('Transfer') grid on

peak = 0.9916 0 0 0 peak = 0.9916 0.8389 0 0 peak = 0.9916 0.8389 0.6232 0 peak = 0.9916 0.8389 0.6232 0.4968

## The Fourier Transform and the linespread function.

% Remember, the linespread was built so that it spans 1 deg, hence % frequency is in cycles per degree. mtf = abs(fft(ls)); hold on, plot(freq,mtf(freq + 1),'ro'); hold off % The values we obtain from convolution are plotted as solid line, % whereas the amplitude of the Fourier Transform of the linespread % function is plotted as a red circles at each frequency. % The functions match, which should give you some intuition about % what the amplitude of the Fourier Transform represents.

## The pointspread and linespread

% When working with two-dimensional inputs, we must consider the % pointspread function, that is the response to an input that is % a point of light. A standard formula for the cross-section of % the pointspread function of the human eye for a 3mm pupil is % also provided by Westheimer. We can compare the linespread and % the cross-section of the pointspread in the following graphs. xSec = -300:300; xMin = xSec/60; ls = 0.47*exp(-3.3 *(xMin.^2)) + 0.53*exp(-0.93*abs(xMin)); ps = 0.952*exp(-2.59*abs(xMin).^1.36) + 0.048*exp(-2.43*abs(xMin).^1.74); vcNewGraphWin; p = plot(xSec,ps,'r-',xSec,ls,'b--'); grid on set(gca,'xlim',[-180 180]) xlabel('Arc sec'), ylabel('LS or PS amplitude') legend(p,'Pointspread','Linespread') % Next, we can create a graph of the pointspread. First, create % a matrix whose entries are the distance from the origin xSec = -240:10:240; xMin = xSec/60; X = xMin(ones(1,length(xMin)),:); Y = X'; D = X.^2 + Y.^2; D = D.^0.5; % If you want to see the distance from the origin, you might show % this image: colormap(gray(32)),imagesc(D), axis image % Then, compute the pointspread function and make a picture of it % ps = 0.952*exp(-2.59*abs(D).^1.36) + 0.048*exp(-2.43*abs(D).^1.74); vcNewGraphWin; colormap(cool(64)), mesh(ps) % To see the pointspread as an image, rather than as a mesh plot, % you might make this figure: colormap(gray(32)),imagesc(ps), axis image

## How the linespread varies with wavelength

% The linespread varies quite strongly with wavelength. When the % eye is in good focus at 580 nm (yellow-part of the spectrum) % the light in the short-wavelength (400-450nm) is blurred quite % strongly and light in the long-wavelength part of the spectrum % is blurred, too, though somewhat less. We can calculate the % linespread as a function of wavelength (Marimont and Wandell, % 1993) from basic principles. The linespreads a various % wavelengths are contained in the data file: load lineSpread; % linespread 370-730nm, in 1nm steps. % This was generated by the script ChromAb.m by Wandell and % Marimont. xDim is in degrees of visual angle. % We select three wavelengths and plot their linespread functions % together. Notice that for the shorter wavelength, the % linespread function is much more spread-out than for the middle % and long wavelengths. vcNewGraphWin; plot(xDim, lineSpread(80, :), 'b-', xDim, lineSpread(200,:), ... 'g:', xDim, lineSpread(361, :), 'r--' ); legend('wavelength 450', 'wavelength 570','wavelength 730'); xlabel('Degrees'); ylabel('Image Intensity'); title('Linespread functions for three wavelengths'); % Look at the line spread functions for all wavelengths lw = 1:10:length(wave); vcNewGraphWin; colormap(hot(32)); mesh(xDim, wave(lw), lineSpread(lw,:)); set(gca,'xlim',[-1 1],'ylim',[350 730]) ylabel('wavelength (nm)'); xlabel('degrees'); zlabel('intensity'); % Different wavelength components of an image are blurred to % different extents by the eye. We use a set of lines again as % an example. For this computation, we will assume that the % input begins with equal energy at all wavelengths from 370 to % 730 nm. % Here, we create and display the image of sample lines. im = reshape([0 0 0 1 0 0 0 0]' * ones(1, 16), 1, 128); imshow(im(ones(100,1), 1:128)); % To calculate the retinal image for this pattern, we % convolve each wavelength component of the image with the % appropriate linespread function. The routine conv2 takes the % input image (im, size(im) = 1 128) and convolves it with each % of the linespread functions in lineSpread size(lineSpread = % 361,65). This results in 361 images (one for each % wavelength). retIm = conv2(im, lineSpread, 'full'); % We must remember the size (in deg) of each sample point this way. % X = (-size(retIm,2)/2 : ((size(retIm, 2)/2) - 1)) / 64; % We can plot the retinal image corresponding to two wavelengths % this way. vcNewGraphWin; subplot(2,1,1) plot(X,retIm(201,:),'g-') set(gca,'ylim',[0 0.5]) title(sprintf('Retinal Image for %d nm', wave(201))); grid on subplot(2,1,2) plot(X,retIm(51,:),'b-') set(gca,'ylim',[0 0.5]) title(sprintf('Retinal Image for %d nm', wave(51))); grid on % Notice that the two images have the same mean, they only differ % in terms of the contrast: the green image has a lot more % contrast, but the same mean. mean(retIm(50,:),2) mean(retIm(200,:),2) % The short wavelength (420 nm) component is blurred much more than % the longer wavelength component of the image. Hence, the image % has very low amplitude ripples, and appears almost like a % single, uniform bar. The 570nm component, however, has high % amplitude ripples that are quite distinct. Hence, the % short-wavelength variation would be very hard to detect in the % retinal image, while the 570 nm component would be quite easy % to detect.

ans = 0.0833 ans = 0.0833

## Chromatic aberration in the frequency domain

% Finally, let's make a few graphs of the modulation transfer % function of the eye's optical system for individual % wavelengths. For short wavelength lights, high spatial % frequency contrast is attenuated a lot by the optical path of % the eye. % Load the MTFs for wavelengths from 370-730nm. These were % calculated using the methods in Marimont and Wandell that are % in a script in the /local/class/psych221/tutorials/chromAb sub-directory. load combinedOtf; % Here is a graph of a few of the MTFs % vcNewGraphWin; plot(sampleSf, combinedOtf(81, :), 'b-', ... sampleSf, combinedOtf(201,:), ... 'g:', sampleSf, combinedOtf(361, :), 'r--' ); legend('wavelength 450', 'wavelength 570','wavelength 730'); xlabel('Frequency (CPD)'); ylabel('Scale factor'); grid on title('Modulation transfer functions for 3 wavelengths'); % Notice that the amplitude of the short-wavelength becomes % negative. This occurs because the blurring is so severe that % the harmonic function is reproduced in the opposite phase % compared to the input harmonic. Hence, the amplitude is % represented by a negative number. This was illustrated in % class using the slide projector, and the phenomenon is called % "spurious resolution."

## More modern Linespreads, Pointspreads, and MTFs

% In recent years, Ijspeert and others in the Netherlands % developed a more extensive set of functions to predict the % basic image formation variables in the average human eye. The % curves they derived were based on empirical inspection of data % sets, and do not have any strong theoretical foundation. % Still, they are probably more accurate than the Westheimer % function and they are parameterized for the subject's age, % pigmentation, and pupil diameter. Hence, for practical % computation, these are probably more useful. % A student in Psych 221 developed the code to compute these % values in the function "ijspeert." Here is an example % age = 20; % Subject's age pupil = 1.5; % diameter in mm pigmentation = 0.142; % Caucasian freqIndexRange = 1:50; % The spatial frequency range % Set the span to be 1 deg, so we know that 1 cycle in the MTF % corresponds to 1 cycle per deg % angleInDeg = (-.25:.005:.25); angleInSec = angleInDeg*3600; angleInRad = angleInDeg*deg2rad; [iMTF, iPSF, iLSF] = ijspeert(age, pupil, pigmentation, ... freqIndexRange, angleInRad); % The functions should be normalized so that the area under the % linespread and the first value of the MTF are one. % iMTF = iMTF/iMTF(1); iLSF = iLSF/sum(iLSF); % These are the modulation transfer function and linespread % for the ijspeert data vcNewGraphWin; %clf subplot(1,2,1), plot(iMTF); grid on set(gca,'xtick',(0:10:80),'xlim',[0 80]); xlabel('Spatial frequency (cpd)'); ylabel('Amplitude'); title('MTF') % Here is the linespread function subplot(1,2,2), plot(angleInSec,iLSF); grid on xlabel('Position (sec)'); ylabel('Intensity'); set(gca,'xtick',(-500:250:500),'xlim',[-500 500]); title('Line spread') % We can compare the ijspeert and westheimer linespread functions % xMin = angleInSec/60; ls = 0.47*exp(-3.3 *(xMin.^2)) + 0.53*exp(-0.93*abs(xMin)); ls = ls / sum(ls); westMTF = abs(fft(ls)); clf subplot(1,2,1) plot(angleInSec/60,iLSF,'b-',angleInSec/60,ls,'r-') set(gca,'xtick',(-8:2:8),'xlim',[-6 6]); xlabel('Position (min)'), ylabel('Intensity'), title('Linespread') grid on

## Comparison of the MTFs with the Williams data

load williams subplot(1,2,2) n = length(iMTF); freq = 0:(n-1); plot(freq,iMTF(1:n),'b-',freq,westMTF(1:n),'r-') hold on plot(dataF,dhb,'ro',dataF,drw,'go',dataF,rnb,'bo') xlabel('Spatial frequency (cpd)'), ylabel('Amplitude'), title('MTF') set(gca,'ylim',[0 1]) grid on hold off legend('Ijspeert','Westheimer','Data')

## What does the pointspread function look like in 2D?

% The PSF is circularly symmetric. So, we can accumulate % the 1D values into a 2D surface. % angleInRad2D = linspace(min(angleInRad)/4,max(angleInRad)/4,length(angleInRad)/4); nSamples = length(angleInRad2D); [X,Y] = meshgrid(angleInRad2D,angleInRad2D); % Because the pointspread function is symmetric, we can calculate the value by % building a matrix D that measures only the distance from the center of matrix. % Try plotting this distance matrix using the command: % % imagesc(D); colorbar; axis image % D = sqrt(X.^2 + Y.^2); % Now, loop through the rows of D to calculate the pointspread values % as a function of angle. % iPSF2D = zeros(size(D)); for ii=1:nSamples a = D(ii,:); [iMTF, iPSF2D(ii,:)] = ijspeert(age, pupil, pigmentation, ... freqIndexRange, a); end vcNewGraphWin; colormap(cool(64)); surf(angleInRad2D,angleInRad2D,iPSF2D);