%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Class:     Psych 221/EE 362
% Tutorial:  Image Formation
% Author:    Wandell
% Purpose:   Introduce visual angle, linespread, pointspread,
% 	     and chromatic aberration
% Date:      01.02.96	
% Duration:  60 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.14.02   Adjusted convolution, adjusted directory location.
%            More fixes shortly, but this is usable. 
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Change the matlab process so that it is executing from inside
% the tutorials directory.  This will make my life easier.

% cd /usr/class/psych221/tutorials
cd ImageFormation

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% 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 occupy (dots per inch)
dpi = 600

% Suppose we read the paper at a viewing distance (inches)
viewingDistance = 12

% The viewing geometry is illustrated in this figure:
%
f = figure
set(f,'Position',[700 600 384 256])
line([0 viewingDistance 0 0],[0 0 1 0]); 
axis equal, axis square
set(gca,'xlim',[-2 20]), grid on
xlabel('Viewing Distance (inch)')
ylabel('Spacing on Paper (inch)')

% In radians, the viewing angle, phi, satisfies 
% tan(phi) = (opposite/adjacent).
deg2rad = 2*pi/360;
rad2deg = 360/(2*pi);
phi = atan(1/viewingDistance)*rad2deg  % phi in degrees

% 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.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% 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);

figure
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.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% 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.

secPerDeg = 60*60;
x = 1:0.2*secPerDeg;
im = zeros(1,length(x));
im(1:30:length(im)) = ones(size(im(1:30:length(im))));

% Here is an image showing the sampled line positions
%
imshow(im(ones(1,128),1:512),32)

% 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');
figure
plot(retIm),grid on
title('The one-dimensional retinal image')
xlabel('Sec of arc'), ylabel('Image intensity')

figure; 
imagesc(retIm(ones(1,128),1:length(retIm))); colormap(gray(256));


% (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 like 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 = mkGaussKernel([1,30],[1 8]);
blurIm = conv2(im,gKernel,'full');
blurIm = scale(blurIm,1,32);
retIm = conv2(blurIm,ls,'full');

% Notice the very small ripples left in the curve after taking
% into account the blurring by the physical display and by the
% eye. 
figure
plot(retIm), axis square, grid on
title('Retinal image of blurred lines')
xlabel('Sec of arc'), ylabel('Intensity')

figure; 
imagesc(retIm(ones(1,128),1:length(retIm))); colormap(gray(256));

% 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?

close all

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Looking at 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);

figure
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);

figure
plot(xMin,harmonic), grid on
title('Sampled cosinusoid')
xlabel('Arc sec'), ylabel('Intensity')

% Now, we will show what happens to 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".

figure
peak = []; freq =[1 5 10 15]; 
for i = 1:length(freq)
 harmonic = cos(2*pi*freq(i)*xMin/nSamples);
% retIm = convolvecirc(harmonic,ls);
 retIm = conv2([harmonic,harmonic],ls,'same');
 subplot(2,2,i)
 plot(retIm), grid on, set(gca,'ylim',[-1 1],'xlim',[0 64]);
 xlabel('Arc sec')
 peak = [peak 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 linesspread is 1).

figure
plot([0 freq],[1 peak],'-')
set(gca,'ylim',[0 1])
xlabel('Spatial freq (cpd)'), ylabel('Transfer')

% Let's compare the values we obtain from the convolution with
% the values we obtain by calculating the amplitude of the
% Fourier Transform of 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.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
% Comparison of 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);

figure
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);
figure
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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
% Chromatic aberration:  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 chromAb/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.

figure
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);
figure
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), 32);

% 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.

figure

subplot(2,1,1)
plot(X,retIm(200,:),'g-')
set(gca,'ylim',[0 0.5])
grid on

subplot(2,1,2)
plot(X,retIm(50,:),'b-')
set(gca,'ylim',[0 0.5])
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,:)')
mean(retIm(200,:)')

% 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.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% 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 chromAb/combinedOtf;

% Here is a graph of a few of the MTFs
% 
figure
plot(sampleSf, combinedOtf(80, :), 'b-', ...
    sampleSf, combinedOtf(200,:), ...
    '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  (01.07.99, BW)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 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 call
% 
age = 20; 				% Subject's age
pupil = 3.0; 				% 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] = ijspeertMTF(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
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

% And here is a comparison of the MTFs with the Williams data superimposed
% 
load data/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,:)] = ijspeertMTF(age, pupil, pigmentation, ...
      freqIndexRange, a);
end

clf;
colormap(cool(64));
surf(angleInRad2D,angleInRad2D,iPSF2D);

% END OF TUTORIAL