%svsim - OMNIVIEWS: SVAVISCA simulation
%
%       function oim = svsim(im,par[,f0,spar])
%
%       im           = input image or N x 2 matrix of points to be transformed
%       oim          = output image
%       par          = [rx dr dt cx cy]
%
%			rx  ... the maximum radius in the image (the boundaries of the mirror)
%			dr  ... the resolution on the r axis in the high resolution image
%			dt  ... the resolution on the theta axis in the high resolution image
%			cx,
%                       cy  ... offset that should be added to image
%		     	        to set the coordinates the way that the center of the 
%			        mirror is also the center of the image.
%
%       spar         = SVAVISCA parameters  = [px py mpx mpy gf nf nr m]
%                    
%                       px  ... pixel size in x axis
%                       py  ... pixel size in y axis
%                       mpx ... minimum pixel size in x axis
%                       mpy ... minimum pixel size in y axis
%                       gf  ... grow factor 
%                       nf  ... number of the rings in fovea
%                       nr  ... number of the rings in perifery  
%                       m   ... number of pixels in each perrifery ring
% 
%                    if not defined, implicit values 
%
%                     spar = [6.87 6.52 6.8 6.45 1.02337 42 110 252] 
%
%                    are used
%
%	f0	     = angle at which starts the unwarping,
%	               imiplicitly 0
%
%       See also:  . 

%       Author       : Hadas Roth, shadasr@techst02.technion.ac.il
%                      12/27/2000 Center for Machine Perception, 
%                      Czech Technical University in Prague
%       Documentation:                            
%       Language     : Matlab 5.2, (c) MathWorks                         
%       Last change  : $Id$
%       Status       : Ready ($Source$)
%
function [oim] = svsim(im,par,f0,spar) 

max_r          = par(1);
res_dr         = par(2);
res_dt         = par(3);
center_image_x = par(4);
center_image_y = par(5);

% general specs of svavisca camera
if nargin < 4
  pixel_f     = [6.87,6.52];    
  min_pixel_r = [6.8,6.45];
  grow_factor = 1.02337;
  Nf          = 42;   		
  Nr          = 110;  		
  M           = 252;  		
else
  pixel_f     = spar(1:2);
  min_pixel_r = spar(3:4);
  grow_factor = spar(5);
  Nf          = spar(6);
  Nr          = spar(7);
  M           = spar(8);
end

if nargin < 3
  f0 = 0;
end

if size(im,2)==2 % transform coordinates from cartesian SVAVISCA polar image
  
  % polar coordinates in original image units
  u     = im(:,1) - center_image_x;
  v     = im(:,2) - center_image_y;
  [t,r] = cart2pol(u,v);
  t     = t+f0;
  
  % Get SAVAVISCA theta and r distributions
  %
  % generating theta and r, should generate seperately for the retina
  % and fovea
  T = cell(1,Nf + Nr);
  R = (0:Nf)*pixel_f(2);
  for n = 0:Nf   % n runs from zero for a pixel in radius zero
    T{n+1} = linspace(0,2*pi,6*n+1);
  end  
  % generating retina
  for n = Nf+1:Nr+Nf
    T{n+1} = linspace(0,2*pi,M+1);
  end
  m     = 1:Nr;
  R     = [R,max(R) + cumsum(min_pixel_r(2)*grow_factor.^m)];
  % normalizing R to fit the image
  max_R = max(R);
  R     = R*max_r/max_R;  
  
  % the number of the ring
  ro    = interp1(R,1:length(R),r,'nearest');
  % theta index
  ti    = T(ro);
  for i=1:length(ti)
    tt     = ti{i};
    if length(ti{i})<2
      to(i)  = 1;
    else
      to(i)  = interp1(tt,1:length(ti{i}),t(i),'nearest'); 
    end
  end
  oim = [ro to'];  
else
% defining the high resolutions image for averaging 
[t_res,r_res] = meshgrid(f0:res_dt:(2*pi+f0),max_r:-res_dr:0);
[x1,y1]       = pol2cart(t_res,r_res);

x1  = x1 + center_image_x; 
y1  = y1 + center_image_y;

% the linear index for resampling to get the high resolution image
ix  = round(y1(:))+(floor(x1(:)))*size(im,1);

% the image in the high resolution coordinates
imx = im(ix);
clear ix;
imo = flipud(reshape(imx,size(t_res,1),size(t_res,2)));

% generating theta and r, should generate seperately for the retina
% and fovea

T = cell(1,Nf + Nr);
R = (0:Nf)*pixel_f(2);

for n = 0:Nf   % n runs from zero for a pixel in radius zero
  T{n+1} = linspace(0,2*pi,6*n+1);
end  

% generating retina
for n = Nf+1:Nr+Nf
  T{n+1} = linspace(0,2*pi,M+1);
end
m     = 1:Nr;
R     = [R,max(R) + cumsum(min_pixel_r(2)*grow_factor.^m)];

% normalizing R to fit the image
max_R = max(R);
R     = R*max_r/max_R;  

% generating  a seperate meshgrid for each line.
t = cell(1,Nf+Nr+1);
r = cell(1,Nf+Nr+1);
x = cell(1,Nf+Nr+1);
y = cell(1,Nf+Nr+1);

for n = 0:Nf+Nr
  [t_temp,r_temp] = meshgrid(T{n+1},R(n+1));
  [x_temp,y_temp] = pol2cart(t_temp,r_temp);
  x_temp          = x_temp + center_image_x;
  y_temp          = y_temp + center_image_y;
  [t{n+1},r{n+1},x{n+1},y{n+1}] = deal(t_temp,r_temp,x_temp,y_temp);
end


% Generating neighbors
%
% defining an array which holds for each pixel the four
% coordinations of the closest neighbors.
% for the central part theta is changed from pixel to pixel, 

R_dif  = (R + [0,R(1:length(R)-1)])/2;
R_dif  = R_dif(2:length(R_dif));

for m = 2:Nf+Nr+1
  for n = 1:(length(t{m})-1)
    if m == Nf+Nr+1
      if n == 1
	neighbor{m,n} = [(t{m}(length(t{m}))+t{m}(length(t{m})-1))/2,(t{m}(n)+t{m}(n+1))/2,R_dif(m-1),max_r];
      else
	neighbor{m,n} = [(t{m}(n)+t{m}(n-1))/2,(t{m}(n)+t{m}(n+1))/2,R_dif(m-1),max_r];
      end
    else
      if n == 1
	neighbor{m,n} = [(t{m}(length(t{m}))+t{m}(length(t{m})-1))/2,(t{m}(n)+t{m}(n+1))/2,R_dif(m-1),R_dif(m)];
      else
	neighbor{m,n} = [(t{m}(n)+t{m}(n-1))/2,(t{m}(n)+t{m}(n+1))/2,R_dif(m-1),R_dif(m)];
      end
    end    
  end
end

% should take care both edges regarding to m (radius)
neighbor{1,1} = R_dif(1);   % the radius of circle that must be averaged.


% Generating the output
for k = 1:Nf+Nr+1
  for l = 1:(length(t{k})-1)
    if l == 1
      if k==1
	output(k,l) = mean(mean(imo(round((0:(size(imo,1)-1)*neighbor{1,1}/max_r)+1),:)));
      else
	output(k,l) = mean(mean(...
		imo(round(((size(imo,1)-1)*neighbor{k,l}(3)/max_r+1):((size(imo,1)-1)*neighbor{k,l}(4)/max_r+1)),...
		    round([((size(imo,2)-1)*neighbor{k,l}(1)/2/pi+1):(size(imo,2)-1),(0:size(imo,2)-1)*neighbor{k,l}(2)/2/pi+1]))...
		));
      end
    else
      output(k,l) = ...
	  mean(mean(...
	      imo(round(((size(imo,1)-1)*neighbor{k,l}(3)/max_r+1):((size(imo,1)-1)*neighbor{k,l}(4)/max_r+1)),...
		  round(((size(imo,2)-1)*neighbor{k,l}(1)/2/pi+1) :((size(imo,2)-1)*neighbor{k,l}(2)/2/pi+1)))...
	      ));      
    end
  end
end
oim = output;
end