function PolarImg = cart2polar(img, numAngle, radRatio, InterpMethod )
%-------------------------------------------------------------------
% transfrom image from Cartesian coordinates to Polar coordinates
% Input:
%    img       -- image in Cartesian coordinates
%    numAngle  -- rows in polar image
%    numRad    -- columns in polar image
% Output:
%    PolarImg  -- image in polar coordinates
%-------------------------------------------------------------------

[Row,Col] = size(img);

Radius = ceil( sqrt( (Row/2)^2 + (Col/2)^2 ) );
if radRatio ~= 1
    Radius = Radius*radRatio;
end
CenterRow = (Row + 1)/2;
CenterCol = (Col + 1)/2;
PolarImg = zeros( Radius, numAngle );

costheta = cos ( (1:numAngle)/numAngle * (2*pi) );
sintheta = sin ( (1:numAngle)/numAngle * (2*pi) );

% nearest interpolation
if strcmp(InterpMethod, 'nearest')         
    for radi = 1 : Radius
        for thetai = 1 : numAngle
            x = round( radi/radRatio*costheta(thetai) + CenterCol );
            y = round( radi/radRatio*sintheta(thetai) + CenterRow );
            if x>0 && y>0 && x<=Col && y<=Row
                PolarImg(radi,thetai) = img(y,x);
            end
        end
    end
end

% bilinear interpolation
if strcmp(InterpMethod, 'bilinear')     
    for radi = 1 : Radius
        for thetai = 1 : numAngle
            rad = radi/radRatio*sintheta(thetai) + CenterRow;
            c   = radi/radRatio*costheta(thetai) + CenterCol;
            rr = floor(rad);
            cc = floor(c);
            x = rad - rr;
            y = c - cc;
            if rr>0 && cc>0 && cc<Col && rr<Row
                Ia = img(rr,cc);
                Ib = img(rr+1,cc);
                Ic = img(rr+1,cc+1);
                Id = img(rr,cc+1);
                PolarImg(radi,thetai) = 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(img);
    for radi = 1 : Radius
        for thetai = 1 : numAngle
            rad =  radi/radRatio*sintheta( thetai ) + CenterRow;
            c    = radi/radRatio*costheta( thetai ) + CenterCol;
            rr = floor(rad);
            cc = floor(c);
            x = rad - rr;
            y = c - cc;
            if rr>0 && cc>0 && cc<Col && rr<Row
                Ia = img(rr,cc);
                Ib = img(rr+1,cc);
                Ic = img(rr+1,cc+1);
                Id = img(rr,cc+1);
                if Mesh(rr,cc)==1 && x>y
                    PolarImg(radi,thetai) = Ia + ( Ib - Ia + ( Ic - Ib) *y) * x;
                elseif Mesh(rr,cc)==1 && x<=y
                    PolarImg(radi,thetai) = Ia + (Ic - Ia)*x + (Id - Ia)*y + (Ia - Id)*x*y;
                elseif Mesh(rr,cc)==0 && x+y<1
                    PolarImg(radi,thetai) = Ia + (Ib - Ia)*x + (Id - Ia)*y + (Ia - Id)*x*y;
                elseif Mesh(rr,cc)==0 && x+y>=1
                    PolarImg(radi,thetai) = 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
