% laplacian.m
%
% Bill Overall
% Mar 12, 1999
%
% Generates foveated and non-foveated output files from an imput
% image.  The compression method used in both cases is a
% Laplacian multiresolution pyramid as described in Burt and
% Adelson (1983).  Uses formula from Geisler and Perry (1998) to
% determine spatial dropoff around a point of gaze.  Saves image
% as integer data, then invokes the unix gzip command in order to
% compress the images.

clear all;
close all;

file = 'out-frnt.jpg' 			% file to foveate
outfile = 1; 				% generate output files from movie?
use_ycbcr = 1; 				% convert RGB to YCbCr?
levels = 8; 				% number of pyramid levels
fovx = 325; 				% x-pixel value of fovea center
fovy = 200; 				% y-pixel value of fovea center
CT0 = 1/75; 				% constant from Geisler&Perry
alpha = 0.106; 				% constant from Geisler&Perry
epsilon2 = 2.3; 			% constant from Geisler&Perry
dotpitch = .25*(10^-3); 		% monitor dot pitch in meters
viewingdist = .3; 			% viewing distance in meters

% read file and store statistics from it
color_img = double(imread(file));

% normalize the values to some maximum value; for most MATLAB
% functions, this value must be one.
max_cval = 1.0;
range = [min(color_img(:)) max(color_img(:))];
img_size = size(color_img(:,:,1));
color_img = max_cval.*(color_img-range(1)) ./ (range(2)-range(1));

% if we're using YCbCr values for foveation, then convert over to them now.
if use_ycbcr == 1
  color_img = rgb2ycbcr(color_img);
end

% ex and ey are the x- and y- offsets of each pixel compared to
% the point of focus (fovx,fovy) in pixels.
[ex, ey] = meshgrid(-fovx+1:img_size(2)-fovx, -fovy+1:img_size(1)-fovy);

% eradius is the radial distance between each point and the point
% of gaze.  This is in meters.
eradius = dotpitch .* sqrt(ex.^2+ey.^2); % in meters
clear ex ey;

% calculate ec, the eccentricity from the foveal center, for each
% point in the image.  ec is in degrees.
ec = 180*atan(eradius ./ viewingdist)/pi;

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

clear eradius;
	      
% calculate the appropriate (fractional) level of the pyramid to use with
% each pixel in the foveated image.

% eyefreq is a matrix of the maximum spatial resolutions (in cpd)
% which can be resolved by the eye
eyefreq = ((epsilon2 ./(alpha*(ec+epsilon2))).*log(1/CT0));

% pyrlevel is the fractional level of the pyramid which must be
% used at each pixel in order to match the foveal resolution
% function defined above.
pyrlevel = maxfreq ./ eyefreq;

% constrain pyrlevel in order to conform to the levels of the
% pyramid which have been computed.
pyrlevel = floor(max(1,min(levels,pyrlevel)));

clear ec maxfreq;

savedata = [];
pyr_file = [];

for color_idx = 1:3
  img = color_img(:,:,color_idx);
  
  % build Gaussian pyramid
  [pyr,indices] = buildLpyr(img,levels);

  % output only those points which are necessary for foveated display
  % in the savedata structure.  For comparison, save all points
  % in the pyramid into the pyr_file structure
  point = 1;
  for n=1:levels
    nextpoint = point + prod(indices(n,:)) - 1;
    show = reshape(pyr(point:nextpoint),indices(n,:))./2^(n-1);
    point = nextpoint + 1;
    this_pyrlevel = pyrlevel(1:2^(n-1):img_size(1),1:2^(n-1):img_size(2));
    requiredpoints = show(find(this_pyrlevel <= n));
    savedata = [savedata requiredpoints(:)'];
    pyr_file = [pyr_file show(:)'];
  end
  
  clear this_pyrlevel requiredpoints;
  clear pyr indices;
  clear show;
  
end

% show uncompressed data savings
pyr_size = size(pyr_file)
savedata_size=size(savedata)

qlevels = 256;

% convert floating-point data range to that of 8-bit integer
% data; i.e., scaled from 0-255.  This quantizes the data and
% therefore introduces quality degradation.  Until this point,
% all manipulations had been percceptually lossless.

pyr_file = (pyr_file + 1.0).*(qlevels/2.0);
savedata = (savedata + 1.0).*(qlevels/2.0);

clear pyrlevel;
clear color_img img;

% write data out to a binary file:
if(outfile == 1)
  u_str = sprintf('U%s',file);
  c_str = sprintf('C%s',file);
  ufid = fopen(u_str,'w');
  cfid = fopen(c_str,'w');
  fwrite(ufid, pyr_file, 'int8');
  fwrite(cfid, savedata, 'int8');
  fclose(ufid);
  fclose(cfid);
end

% Compress them both

tempstr = sprintf('gzip -c %s > %s.gz',u_str,u_str);
unix(tempstr);
tempstr = sprintf('gzip -c %s > %s.gz',c_str,c_str);
unix(tempstr);

% Look at their compressed sizes.
tempstr = sprintf('ls -l %s* %s*',u_str,c_str);
unix(tempstr);
