% Solve the cases in the J. Opt. Soc. Am. A 2005 paper: Wave-front reconstruction from its gradients.
% S2p package: functions for calculation of phase phi from shears Sx, Sy
% Inputs: 
% 	Nx, Ny: dimensions
% 	Sx,Sy = shears, differences in phase, of size (1:Nx,1:Ny)
% 	    Sx(1:Ny,1:Nx-1) = phi(1:Ny,2:Nx)-phi(1:Ny,1:Nx-1)
% 	    Sx(1:Ny,Nx) = 0  or some other method of extension
% 	    
% 	    Sy(1:Ny-1,1:Nx) = phi(2:Ny,1:Nx)-phi(1:Ny-1,1:Nx)
% 	    Sx(Ny,1:Nx) = 0  or some other method of extension
% Output:  
% 	phi(1:Ny,1:Nx) = phase
% 
% functions:
%   phi=cos_s2p (Sx,Sy) - current nonperiodic method
%   phi=fft_s2p (Sx,Sy) - previous periodic FFT
%   
%   [var,dsx,dsy]=err_s2p(phi,Sx,Sy) - calculate the errors (how good a solution)
%   var = variance = sum{ |dsx|^2 + |dsy|^2 }
%   dsx = difference = phi(y,x+1)-phi(y,x)-Sx(y,x)
%   dsy = difference = phi(y+1,x)-phi(y,x)-Sy(y,x)
%
% Amos Talmi and Erez Ribak 2005

% Hit any key to continue
help main_s2p, %pause
close all; clear all
Nx = 64; Ny = 128; 		% define array size and other useful arrays
X1 = [1:Nx]; Y1 = [1:Ny]; oneX = ones (1, Nx); oneY = ones (1, Ny);
S = zeros (Ny, Nx);  		% initialize slope arrays to zero
preparrays 			% create analytic and random arrays and calc slopes without and with noise 

sxc = ['oneY'' * oneX                         ' ; ...               % given x slopes
        'oneY'' * (X1+10)                      ' ; ...
        '(Y1+10)''*oneX                        ' ; ...
        '(Y1-16.5)'' * ((X1/64).* (1- (X1/64)))' ; ...
        'Px                                   ' ; ...
        'Ax                                   ' ; ...
        'Bx                                   '];
syc = ['S ' ; ...                                                  % given y slopes
        'S ' ; ...
        'S ' ; ...
        'S ' ; ...
        'Py' ; ...
        'Ay' ; ...
        'By'];
pc = ['S' ; ...                                                     % given y slopes
        'S' ; ...
        'S' ; ...
        'S' ; ...
        'P' ; ...
        'A' ; ...
        'A'];
cases = ['Sx = 1, Sy = 0                            ' ; ...    % list given slopes
        'Sx = n+10 ,  Sy = 0                       ' ; ...
        'Sx(n+½, m) = m+10, Sy = 0                 ' ; ...
        'Sx(n+½, m) = (m-16.5)*n/64(1-n/64), Sy = 0' ; ...
        'Analytic input                            ' ; ...
        'Random input                              ' ; ...
        'Random input with noise on slopes         '];

format short; format compact 
for c = 1:7
    Sx = eval (sxc (c, :));  		% use x slopes above
    Sy = eval (syc (c, :)); 		% use y slopes above
    phi = eval (pc (c, :));  		% use function above (analytic, random cases)
    
    phi1 = cos_s2p (Sx, Sy);  		% calc function from slopes, our method
    phi2 = fft_s2p (Sx, Sy);  		% calc function from slopes, periodic 
    
    [var1, dsx, dsy] = err_s2p (phi1, Sx, Sy);   % errors in our method
    max_dev_cos = max (max (max (abs (dsx))), max (max (abs (dsy))));
    [var2, dsx, dsy] = err_s2p (phi2, Sx, Sy);   % errors in FFT method
    max_dev_fft = max (max (max (abs (dsx))), max (max (abs (dsy))));
    
    std_cos = sqrt (abs (var1)/Nx/Ny); 		% normalised standard deviations
    std_fft = sqrt (abs (var2)/Nx/Ny); 		% normalised standard deviations
%   disp (['case ' int2str(c) ': ' cases(c,:)]) % list cases above  
%   disp ('STDcos MaxDevCos STDFFT MaxDevFFT' )
%   disp ([std_cos, max_dev_cos, std_fft, max_dev_fft])
    fprintf('Case %d: %s. STDcos: %8.4g; MaxDevCos: %8.4g; STDFFT: %8.4g; MaxDevFFT: %8.4g\n', c, cases(c,:),[max_dev_cos, max_dev_fft, std_cos, std_fft])
    
    if c > 4 % analytic and random functions, without and with slopes noise
        figure (c); 
        subplot (131); imagesc(phi ); axis image;axis xy; colorbar h; title (cases(c,:))
if c>6,   subplot (131); imagesc(phi1); axis image;axis xy; colorbar h; title (cases(c,:)), end
        phi1 = phi1 - phi; phi1 = (phi1 - min (phi1(:))) .* R2; 	% difference from input
        phi2 = phi2 - phi; phi2 = (phi2 - min (phi2(:))) .* R2; 	% difference from input
        subplot (132); imagesc(phi1); axis image;axis xy; colorbar h; title ('cos error')
        subplot (133); imagesc(phi2); axis image;axis xy; colorbar h; title ('fft error')
    end
    
end
format