function Ufinal = unfold(Uoriginal,jumpvalue,dim)
% Unwrapping algorithm using the unfolding procedure 
%
% Starts with function UORIGINAL, delivers function UFINAL as result of the
% unwrapping process with jumps larger than PI to the 2PI complement.
% UORIGINAL can be a vector, matrix or N-D array
%
% unfold(Uoriginal,jumpvalue) returns the unwrapped face with JUMPVALUE
% tolerance; default value is jumpvalue = 2*pi
%
% unfold(Uoriginal,[],dim) unwrapps along the DIM dimension for 
% discontinuities of jumpvalue magnitude. unfold(Uoriginal,jumpvalue,dim) uses 
% the tolerance jumpvalue
%
% Ref: "Unfolding wrapped phase" 
%       Opt. Eng. 54(11), 110503 (Nov 20, 2015). doi:10.1117/1.OE.54.11.110503
%
% Carlos G Treviño Palacios
% INAOE abril,2015

ni = nargin;

% Treat row vector as a column vector (unless DIM is specified)
rflag = 0;
if ni<3 && (ndims(Uoriginal)==2) && (size(Uoriginal,1)==1), 
   rflag = 1; 
   Uoriginal = Uoriginal.';
end

% Initialize parameters.
nshifts = 0;
perm = 1:ndims(Uoriginal);
switch ni
case 1
   [Uoriginal,nshifts] = shiftdim(Uoriginal);
   jumpvalue = 2*pi;     
case 2
   [Uoriginal,nshifts] = shiftdim(Uoriginal);
otherwise    % nargin == 3
   perm = [dim:max(ndims(Uoriginal),dim) 1:dim-1];
   Uoriginal = permute(Uoriginal,perm);
   if isempty(jumpvalue),
      jumpvalue = 2*pi; 
   end
end

% Reshape UORIGINAL to a matrix.
siz = size(Uoriginal);
Uoriginal = reshape(Uoriginal, [siz(1) prod(siz(2:end))]);

%Unwrap each column of vector UORIGINAL

for otro=1:size(Uoriginal,2)

jump = jumpvalue; %Reset magnitude of the discontinuity; default value 2*pi
jump2=jump/2;
jump4=jump/4;

%ORIGINAL function column by column

UM=Uoriginal(:,otro);

for rep=1:100 % Iterate for functions with amplitude up to 2^100 
    
    MN=UM-mod(UM+jump2,jump);

    if (max(MN)-min(MN))<jump4
     break
    end
    
    shift=0;
    flag=0;
    for jj=2:size(Uoriginal,1)
        if MN(jj)-MN(jj-1)~=0
%Identify the first discontinuity
            if flag==0
                signo=sign(UM(jj)-UM(jj-1));
                flag=1;
            end
%Unfolding condition
            if abs(UM(jj)-UM(jj-1)+shift)>jump2
                if shift==0
                    shift=signo*jump;
                else
                    shift=0;
                end
            end
        end
        UM(jj)=UM(jj)+shift;
    end
 
    if signo~=1 %Displacement for minus sign
        UM=UM+jump;
    end
    
    jump2= jump; 
    jump = jump*2; %Scales the discontinuity magnitude
    jump4= jump/4;
    
end % Unfold iterations
Uout(:,otro)=UM;
end %All columns

%Reshapes the matrix to the result UFINAL

Uout = reshape(Uout,siz);
Uout = ipermute(Uout,perm);
Uout = shiftdim(Uout,-nshifts);
if rflag, 
   Uout = Uout.'; 
end

Ufinal=Uout;