clear all;
close all;
im=imread('m640.tif','tiff');
%imshow(im);


[xs,ys]=size(im(:,:,1));
fovx=round(xs/2);
fovy=round(ys/2);
CT0 = 1/75;                     % constant from paper
%alpha = 0.106;                  % constant from paper
alpha = 0.106*1.5;                  % constant from paper
epsilon2 = 2.3/1.5;                 % constant from paper
dotpitch = 0.25*(10^-3);         % in meters
viewingdist = 0.25;              % in meters

% calculate ec, the eccentricity from the foveal center, for each
% point in the image.  ec is in degrees.
[ex, ey] = meshgrid(-fovx+1:xs-fovx, -fovy+1:ys-fovy);
eradius=dotpitch.*sqrt(ex.^2+ey.^2);
ec = 180*atan((dotpitch .* sqrt(ex.^2 + ey.^2)) ./ viewingdist)/pi+0.001;

% maximum spatial frequency which can be adequately represented onscreen:
maxfreq = pi./(atan(2.*dotpitch./sqrt(viewingdist^2+eradius.^2))*180); 
%maxfreq=pi./((atan((eradius+dotpitch)./viewingdist)-atan((eradius-dotpitch)./viewingdist)).*180);

for color=1:3

p=double(im(:,:,color));

[xs,ys]=size(p);

% Build the pyramid of the input image:

[pyr pind] = buildGpyr(p);

% find the maximum pyramid height:

[maxh junk]=size(pind);
levels=maxh;
% make a gaussian foveal blob, scale to height of 1.


% calculate the appropriate (fractional) level of the pyramid to use with
% each pixel in the foveated image.
%levelmat = maxfreq ./ (((epsilon2 ./ (alpha.*ec)).*log(1/CT0)) - epsilon2);
fmat = ((((epsilon2 ./ alpha).*log(1/CT0)))./(ec+epsilon2));
levelmat = maxfreq ./ fmat;
fovea=max(1,min(levels, max(levelmat,1)));

% make a 3-D matrix of the image size by the # of pyramid levels:

psc=zeros(xs,ys,maxh);

% fill the 3-D array with pyramid levels, upblurring levels to
% equalize heights:

for k=1:min(maxh,ceil(max(max(fovea))));
 psc(:,:,k)=imcrop(upBlur(pyrBand(pyr,pind,k),k-1),[1 1 xs-1 ys-1]);
end;

% Do a 3-D interpolation by foveation magnitude between
% the appropriate pyramid levels.

x=meshgrid(1:1:ys);
x=x(1:xs,:);
y=meshgrid(1:1:ys)';
y=y(1:xs,:);

b=interp3(psc,x,y,fovea,'linear');

if color==1
  r=b;
end;
if color==2;
  g=b;
end;
end;

rgb=zeros(xs,ys,3);
rgb(:,:,1)=r/255;
rgb(:,:,2)=g/255;
rgb(:,:,3)=b/255;
figure;
imshow(rgb);
imwrite(rgb,'out.tif','tiff');