function ReconImg = polar2cart( PolarImg, Row, Col, radRatio, InterpMethod)
% transfrom image from Polar Coordinates to Cartesian Coordinates

% circle padding for processing the boundary  between 0 degee and 360 degree
ThetaExtSize  = 5;
PolarImg = [PolarImg(:,end-ThetaExtSize+1:end), PolarImg,  PolarImg(:,1:ThetaExtSize)];
% pad for computing center
RadExtSize  = 1;
PolarImg = [PolarImg(RadExtSize:-1:1,:); PolarImg];   

%% ******************************
[numRad,numAngle] = size(PolarImg);
CenterRow = (Row + 1)/2;
CenterCol = (Col + 1)/2;
ReconImg = zeros( Row, Col );
%% nearest interpolation
if strcmp(InterpMethod, 'nearest')
    for j = 1 : Col
        for i  = 1 : Row
            x = j - CenterCol;
            y = i - CenterRow;
            radi = round( sqrt( x^2 + y^2 )*radRatio );
            
            if i >= CenterRow && j > CenterCol        % first quadrant
                theta = atan(y/x);
            elseif i > CenterRow && j < CenterCol     % second quadrant
                theta = pi - abs( atan(y/x) );
            elseif i <= CenterRow && j < CenterCol     % third quadrant
                theta = pi + atan(y/x);
            elseif i < CenterRow && j > CenterCol     % fourth quadrant
                theta = 2*pi - abs( atan(y/x) );
            elseif i >= CenterRow && j == CenterCol
                theta = 0.5*pi;
            elseif i < CenterRow  && j == CenterCol
                theta = 1.5*pi;
            else
                fprintf('Transform polar image to coordinate image error!\n');
                return;
            end
            thetai = round( theta/(2*pi) * (numAngle - ThetaExtSize*2) ) ;
            
            if radi>=0 && thetai>=0 && radi<=numRad && thetai<=numAngle
                ReconImg(i,j) = PolarImg(radi+RadExtSize, thetai+ThetaExtSize);
            end
            
        end
    end
end

%% bilinear interpolation
if strcmp(InterpMethod,  'bilinear')  
    for j = 1 : Col
        for i = 1 : Row
            x = j - CenterCol;
            y = i - CenterRow;
            radi = sqrt( x^2 + y^2 )*radRatio;
            
            if i >= CenterRow && j > CenterCol        % first quadrant
                theta = atan(y/x);
            elseif i > CenterRow && j <CenterCol     % second quadrant
                theta = pi - abs( atan(y/x) );
            elseif i <= CenterRow && j < CenterCol     % third quadrant
                theta = pi + atan(y/x);
            elseif i < CenterRow && j > CenterCol     % fourth quadrant
                theta = 2*pi - abs( atan(y/x) );
            elseif i >= CenterRow && j == CenterCol
                theta = 0.5*pi;
            elseif i < CenterRow  && j == CenterCol
                theta = 1.5*pi;
            else
                fprintf('Transform polar image to coordinate image error!\n');
                return;
            end
            thetai =  theta/(2*pi) * (numAngle - ThetaExtSize*2) ;
            
            rr = floor(radi);
            cc = floor(thetai);
            x = radi - rr;
            y = thetai - cc;
            cc = cc + ThetaExtSize;
            rr = rr + RadExtSize;
            if rr>0 && rr<=numRad && cc>=0  && cc<=numAngle
                Ia = PolarImg(rr,cc);
                Ib = PolarImg(rr+1,cc);
                Ic = PolarImg(rr+1,cc+1);
                Id = PolarImg(rr,cc+1);
                ReconImg(i,j) = Ia + (Ib - Ia)*x + (Id - Ia)*y + (Ia - Ib + Ic - Id)*x*y;
            end
            
        end
    end
end

%% DDT interpolation
if strcmp(InterpMethod,  'DDT' )   
    Mesh = triangulation_mesh(PolarImg);
    Mesh = [Mesh,Mesh(:,1)];
    for j = 1 : Col
        for i = 1 : Row
            x = j - CenterCol;
            y = i - CenterRow;
            radi = sqrt( x^2 + y^2 )*radRatio;
            
            if i >= CenterRow && j > CenterCol        % first quadrant
                theta = atan(y/x);
            elseif i > CenterRow && j <CenterCol     % second quadrant
                theta = pi - abs( atan(y/x) );
            elseif i <= CenterRow && j < CenterCol     % third quadrant
                theta = pi + atan(y/x);
            elseif i < CenterRow && j > CenterCol     % fourth quadrant
                theta = 2*pi - abs( atan(y/x) );
            elseif i >= CenterRow && j == CenterCol
                theta = 0.5*pi;
            elseif i < CenterRow  && j == CenterCol
                theta = 1.5*pi;
            else
                fprintf('Transform polar image to coordinate image error!\n');
                return;
            end
            thetai =  theta/(2*pi) * (numAngle - ThetaExtSize*2);
            
            rr = floor(radi);
            cc = floor(thetai);
            x = radi - rr;
            y = thetai - cc;
            cc = cc + ThetaExtSize;
            rr = rr + RadExtSize;
            if rr>0 && cc>=0 && rr<=numRad && cc<=numAngle
                Ia = PolarImg(rr,cc);
                Ib = PolarImg(rr+1,cc);
                Ic = PolarImg(rr+1,cc+1);
                Id = PolarImg(rr,cc+1);
                if Mesh(rr,cc)==1 && x>y
                    ReconImg(i,j) = Ia + ( Ib - Ia + ( Ic - Ib) *y) * x;
                elseif Mesh(rr,cc) == 1 && x<=y
                    ReconImg(i,j) = Ia + (Ic - Ia)*x + (Id - Ia)*y + (Ia - Id)*x*y;
                elseif Mesh(rr,cc) == 0 && x+y<1
                    ReconImg(i,j) = Ia + (Ib - Ia)*x + (Id - Ia)*y + (Ia - Id)*x*y;
                elseif Mesh(rr,cc) == 0 && x+y>=1
                    ReconImg(i,j) = Id + ( Ib - Id + (Ic - Ib) *y) * x;
                end
            end
            
        end
    end
end

%% function triangulation_mesh for PLDDT
function ex_Mesh = triangulation_mesh(Img)
% Reference: Dan Su and Philip Willis. Image Interpolation by Pixel-Level Data-Dependent Triangulation[J]. Computer graphics forum, 2004,23(2): 189-201.
% Input    : Img  - image to be interpolated
% Output   : Mesh - binary image which represent triangulation mesh direction in the image
% Notation: We define NW--SE = 1; SW--NE = 0

[ROW,COL] = size(Img);
Mesh = zeros(ROW-1,COL-1,'uint8');

% calculate square grid direction of four pixels
for j = 1 : COL-1
    for i = 1 : ROW-1
        if abs( Img(i,j) - Img(i+1,j+1) ) < abs( Img(i+1,j) - Img(i,j+1) )
            Mesh(i,j) = 1;
        end
    end
end

% extend method, neighborhood of 3x3 squares
ex_Mesh = Mesh;
for j = 2 : size(Mesh,2)-1
    for i = 2 : size(Mesh,1)-1
        if sum( sum( Mesh(i-1:i+1, j-1:j+1) )) >= 7
            ex_Mesh(i,j) = 1;
        elseif sum( sum( Mesh(i-1:i+1, j-1:j+1) )) <= 2
            ex_Mesh(i,j) = 0;
        end
    end
end
