%% interp_x2_3D.m
%
% E. Koenig, 08/04/2010
%
% Interpolation by a factor 2 of a 3D field.
%
% The interpolation can be done using 3 methods :
%
% *  *nearest*: the new point values are copies of their left nearest 
% neighboor
% *  *linear*:  the new point values are linearly interpolated from the
% neighboors
% *  *splines*: the new point values are interpolated using splines. This
% method is based on the works of Unser et al [1]. 
%
%
% Needs the function |mexinterp|, which has to be compiled, and is placed in
% the folder interpol_2D.
%
%
%% References
%
% |[1] Unser, M., Aldroubi, A. Eden, M.: Enlargement or reduction of digital 
% images with minimum loss of information. IEEE Transactions on Image 
% Processing 4(3), 247-258 (1995)|
%
%
%% Use
%
% I_interp = interp_x2_3D(I, method, conservation);
%
%
%% Inputs
%
% * *I*:  input 3D field
% * *method*: method used for the interpolation : 'nearest', 'linear' or
% 'splines'
% * *conservation*: 1 if the interpolation is conservative, 0 if not
%
% By default, the interpolation is made by the linear method with no
% conservation.
%
%
%% Outputs
%
% * *I_interp*: interpolated 3D field. The size of this resulting field is
% twice the size of I
%
%
%% Example
%
% Interpolation of a 4x4x4 random field using the splines method :
%
%   I = randn(4,4,4);
%   I_interp = interp_x2_3D(I, 'splines', 0);
%




function I_interp = interp_x2_3D(I, method, conservation)

% Default
if nargin<3
    conservation = 0 ;
    if nargin<2
        method = '*linear';
    end
end

SS = size(I);

% The intermediate field is the same field as the initial one with the
% first and last row, columns and depth doubled. This leads to a field of
% size SS+2
I_tmp = zeros(SS(1)+2,SS(2)+2,SS(3)+2);
I_tmp(2:end-1, 2:end-1, 2:end-1) = I;
I_tmp(1,:,:) = I_tmp(2,:,:);
I_tmp(end,:,:) = I_tmp(end-1,:,:);
I_tmp(:,1,:) = I_tmp(:,2,:);
I_tmp(:,end,:) = I_tmp(:,end-1,:);
I_tmp(:,:,1) = I_tmp(:,:,2);
I_tmp(:,:,end) = I_tmp(:,:,end-1);

S = size(I_tmp);


% Interpolation step
if(length(method) == length('*splines') & method == '*splines')
    % In the case of the splines method
    
    xi = 3/4:1/2:(S(1)-2+1/4);
    yi = 3/4:1/2:(S(2)-2+1/4);
    zi = 3/4:1/2:(S(3)-2+1/4);

    I_interp = zeros(SS.*2);

    % a first interpolation is made in the planes of the field taken along
    % the z-axis
    for kp = 1:S(3)
        I_int(:,:,kp) = mexinterp(squeeze(I_tmp(:,:,kp)), xi, yi, 3);
    end

    % a second interpolation is made for each depth previously interpolated
    for ip = 1:size(I_interp,1)
        for jp = 1:size(I_interp,2)
            I_interp(ip, jp, :) = mexinterp(squeeze(I_int(ip, jp, :)), zi, 3);
        end
    end
    
    clear xi yi zi I_int I_tmp S
    
else
    % In the case of the others methods
    
    x = 0:S(1)-1;
    y = 0:S(2)-1;
    z = 0:S(3)-1;
    xi = 3/4:1/2:(S(1)-2+1/4);
    yi = 3/4:1/2:(S(2)-2+1/4);
    zi = 3/4:1/2:(S(3)-2+1/4);
    [XI,YI,ZI] = meshgrid(xi,yi,zi);
    
    I_interp = interpn(x,y,z,I_tmp,XI,YI,ZI,method);

    clear x y z XI YI ZI xi yi zi S I_tmp

end

% Conservation step
if conservation
    % zoom out to the initial size
    I_moy_red = zoom_arriere(I_interp,2,0) ;
    I_norm = I(2:end-1,2:end-1) ./ I_moy_red ;
    clear I_moy_red
    % back to the interpolated resolution using the nearest method
    I_norm = interpol_image_fois2_3D_tmp(I_norm,'*nearest') ;
    I_interp = I_interp .* I_norm ;
    clear I_norm
end

clear method conservation I
