function [img] = unwarp(warped_img,coord_frame, current_shape,orig_img)
img=orig_img;
base_shape = coord_frame.base_shape;
base_texture = coord_frame.base_texture;
triangles = coord_frame.triangles;
resolution = coord_frame.resolution;

num_of_triangles = size(triangles, 1);
A = zeros(1,6); % affine transformation for each triangle
O= warp_image(coord_frame, current_shape, orig_img);
m1=mean(O(O~=0));
m=mean(warped_img(warped_img~=0));
m=m1/m;
for t = 1:num_of_triangles
    %% Step 1:
    imagesc(img); colormap(gray);
    % Coordinates of the three vertices of each triangle in base shape
    X = base_shape(triangles(t,:),1) - 1;
    Y = base_shape(triangles(t,:),2) - 1;
    
    % Coordinates of the three vertices of the corresponding triangle in current shape
    U = current_shape(triangles(t,:),1) - 1;
    V = current_shape(triangles(t,:),2) - 1;
    
    % Compute A from the U,V and X,Y
    denominator = (U(2) - U(1)) * (V(3) - V(1)) - (V(2) - V(1)) * (U(3) - U(1));
    
    A(1) = X(1) + ((V(1) * (U(3) - U(1)) - U(1)*(V(3) - V(1))) * (X(2) - X(1)) + (U(1) * (V(2) - V(1)) - V(1)*(U(2) - U(1))) * (X(3) - X(1))) / denominator;
    A(2) = ((V(3) - V(1)) * (X(2) - X(1)) - (V(2) - V(1)) * (X(3) - X(1))) / denominator;
    A(3) = ((U(2) - U(1)) * (X(3) - X(1)) - (U(3) - U(1)) * (X(2) - X(1))) / denominator;
    
    A(4) = Y(1) + ((V(1) * (U(3) - U(1)) - U(1) * (V(3) - V(1))) * (Y(2) - Y(1)) + (U(1) * (V(2) - V(1)) - V(1)*(U(2) - U(1))) * (Y(3) - Y(1))) / denominator;
    A(5) = ((V(3) - V(1)) * (Y(2) - Y(1)) - (V(2) - V(1)) * (Y(3) - Y(1))) / denominator;
    A(6) = ((U(2) - U(1)) * (Y(3) - Y(1)) - (U(3) - U(1)) * (Y(2) - Y(1))) / denominator;
    
    
    %% Step 2
    % Get coordinates of all pixels within each triangle
%     [v, u] = find(base_texture == t);
    [u,v]= triag_list([U(1),V(1)],[U(2),V(2)],[U(3),V(3)]);
    if (~isempty(v) && ~isempty(u))
        
        ind_base = v + (u-1) *size(orig_img, 1); 
        
        %% Step 3
        v = v - 1; u = u - 1;
        warped_x = A(1) + A(2) .* u + A(3) .* v + 1;
        warped_y = A(4) + A(5) .* u + A(6) .* v + 1;
        
        % round the values
        ind_warped = round(warped_y) + (round(warped_x)-1) * resolution(1);
        
        %% Step 4
        % Find pixel coordinates in warped image
        max_ind_warped = resolution(1)*resolution(2);
        list = (ind_warped > 0 & ind_warped <= max_ind_warped);
        empty=find(warped_img(ind_warped(list))==0);
        list(empty)=0;
        img(ind_base(list))=warped_img(ind_warped(list))*m;
    end
end
% imshow(img,[],'Border','tight');
end

function y= intpoint_inside_trigon(s,a,b,c)
    as_x = s(1)-a(1);
    as_y = s(2)-a(2);
    s_ab = (b(1)-a(1))*as_y-(b(2)-a(2))*as_x > 0;
    y=0;
    if((c(1)-a(1))*as_y-(c(2)-a(2))*as_x > 0 == s_ab) 
        return;
    end;

    if((c(1)-b(1))*(s(2)-b(2))-(c(2)-b(2))*(s(1)-b(1)) > 0 ~= s_ab) 
        return;
    end;
    y=true;
    return;
end


function [u,v]= triag_list(a,b,c)
x1= round(min([a(1),b(1),c(1)]));
x2= round(max([a(1),b(1),c(1)]));
y1= round(min([a(2),b(2),c(2)]));
y2= round(max([a(2),b(2),c(2)]));
idx=1;
u=zeros((x2-x1)*(y2-y1),1);
v=zeros((x2-x1)*(y2-y1),1);
for i=x1:x2
    for j=y1:y2
        if (intpoint_inside_trigon([i,j],a,b,c))
            u(idx)=i;
            v(idx)=j;
            idx=idx+1;
        end
    end
end;
u(u==0)=[];
v(v==0)=[];
end