function outputIm = backward_geometry(inputIm, A, type)
% inputIm = 输入的图像
%       A = 仿射变换的系数，一个2x3的矩阵

% 获取输入图像的大小
inputSize = size(inputIm);
if(size(inputIm, 3) == 1)
   inputSize(3) = 1; 
end
h=inputSize(1);
w=inputSize(2);

% 计算输出图像的画布大小
[outputSize, deltaShift] = calcOutputSize(inputSize, A, type);
outputIm=zeros(outputSize(2),outputSize(1));

% expand the previous image
inputImEx=zeros(h+2,w+2,1);
inputImEx(2:h+1,2:w+1,1)=inputIm(:,:,1);
inputImEx(1,2:w+1,1)=inputIm(1,:,1);
inputImEx(2:h+1,1,1)=inputIm(:,1,1);
inputImEx(h+2,2:w+1,1)=inputIm(h,:,1);
inputImEx(2:h+1,w+2,1)=inputIm(:,w,1);
inputImEx(1,1,1)=inputIm(1,1,1);
inputImEx(1,w+2,1)=inputIm(1,w,1);
inputImEx(h+2,1,1)=inputIm(h,1,1);
inputImEx(h+2,w+2,1)=inputIm(h,w,1);

% 根据确定的输出画布大小来进行遍历
for i = 1 : outputSize(1)
    for j = 1 : outputSize(2)
        y = j;
        x = i;
        % 进行逆向变换，计算当前点(x,y)在输入图像中的坐标
        A(3,:)=[0,0,1];
        vec=[x-deltaShift(1);y-deltaShift(2);1];
        vec0 = A\vec;
        x0=vec0(1); y0=vec0(2);
        % 进行双线性插值获取像素点的值
        if x0>0 && x0<=w && y0>0 && y0<=h
            xf=floor(x0)+1; xc=xf+1;
            yf=floor(y0)+1; yc=yf+1;
            u=x0+1-xf; v=y0+1-yf;
            res=u*v*inputImEx(yc,xc)+u*(1-v)*inputImEx(yf,xc)+...
                    (1-u)*v*inputImEx(yc,xf)+(1-u)*(1-v)*inputImEx(yf,xf);
            outputIm(y,x,1)=round(res);
        end
    end
end
outputIm=uint8(outputIm);
end


function [outputSize, deltaShift] = calcOutputSize(inputSize, A, type)
% type 有两种，一种是 loose， 一种是crop，参考imrotate命令的帮助文件
% 需要实现这两种
% 'crop'
% Make output image B the same size as the input image A, cropping the rotated image to fit
% {'loose'}
% Make output image B large enough to contain the entire rotated image. B is larger than A

% 获取图像的行和列的总数，其中行方向对应着y方向，列方向对应着x方向    
ny = inputSize(1);
nx = inputSize(2);

% 计算四个顶点的齐次坐标
inputBoundingBox = [ 1  1 1;...
                    nx  1 1;...
                    nx ny 1;...
                     1 ny 1];
inputBoundingBox = inputBoundingBox';

% 获取输入图像经过仿射变换后在输出图像中的框
outputBoundingBox = A * inputBoundingBox;

% 找到输出图像的紧致的框
xlo = floor(min(outputBoundingBox(1,:)));
xhi =  ceil(max(outputBoundingBox(1,:)));
ylo = floor(min(outputBoundingBox(2,:)));
yhi =  ceil(max(outputBoundingBox(2,:)));

if strcmp(type,'loose')==1
    outputSize(1) = xhi-xlo+1;
    outputSize(2) = yhi-ylo+1;
    deltaShift(1) = -xlo+1;
    deltaShift(2) = -ylo+1;
elseif strcmp(type,'crop')==1
    outputSize(1) = inputSize(2);
    outputSize(2) = inputSize(1);
    deltaShift(1) = -(xlo+xhi-outputSize(1)-1)/2;
    deltaShift(2) = -(ylo+yhi-outputSize(2)-1)/2;
end
end

type为'crop'或者'loose'，最终输出灰度图。