getrectifiedimage.asv

这是我用Delphi和Matlab写的一个程序

function [imH0 imH1] = GetRectifiedImage(im0, H0, im1, H1);
%  GETRECTIFIEDIMAGE  Return the rectified image
%  im0、im1: 图像数据
%  H0、H1: the projective transformation matrices
%  imH0、imH1: 校正后的图像
%
%  Ran Liu: [e-mail]
%  College of Computer Science, Chongqing University
%  Panovasic Technology Co.,Ltd

%%  预处理图像im0，得到Rect0(是一个MBR)
Rect0 = Preprocess(im0, H0);
Rect1 = Preprocess(im1, H1);

%%  求解平移矩阵
[Trans0, Trans1] = GetTransMats(Rect0, Rect1);

%%  生成图像
imH0 = ProduceImg(im0, H0, Rect0, Trans0);
imH1 = ProduceImg(im1, H1, Rect1, Trans1);

function Rect = Preprocess(Ii, Hi)
%  PREPROCESS  Return the Preprocessed image
%  计算图像HiIi
[Hight Width Clr] = size(Ii);  %  求Ii各维的大小
Img = zeros(3, Hight * Width);  %  构造一个3×(Hight*Width)的零矩阵
for i = 1 : Hight  %  Img的每一列是点经过变换后的齐次坐标，按行扫描
    for j = 1 : Width  % j: 列
      Img(:, i * j)  = Hi * [j; i; 1];  %  [j; i; 1]: 第i行第j列的像素在像素图像坐标系中的齐次坐标（列向量）
                                       %  理论上不必让让齐次坐标的h等于1，参见论文《立体图像对的极线校正》
    end
end;

%  求解HiIi的最小边界矩形
TopLeftX = round(min(Img(1, :)));
TopLeftY = round(min(Img(2, :)));
BottomRightX = round(max(Img(1, :)));
BottomRightY = round(max(Img(2, :)));
clear Img;
Rect = [TopLeftX, TopLeftY, BottomRightX, BottomRightY];

function [Trans0, Trans1] = GetTransMats(Rect0, Rect1)
%  GETTRANSMATS  Return the translation  matrixes
if (abs(Rect0(2)) > abs(Rect1(2)))
    TransY = Rect0(2);
else
    TransY = Rect1(2);
end
Trans0 = [1   0  (1 - Rect0(1))
          0   1  (1 - TransY)
          0   0   1   ];
      
Trans1 = [1   0  (1 - Rect1(1))
          0   1  (1 - TransY)
          0   0   1   ];     

function imH = ProduceImg(im, H, Rect, Trans, Hight, Width, Clr);
% 以从左到右、从上到下的顺序，计算Recti中的每一个点在原始图像平面中的对应点
% Rect = [TopLeftX, TopLeftY, BottomRightX, BottomRightY];
[Hight Width Clr] = size(Ii); 
TopLeftX = Rect(1);
TopLeftY = Rect(2);
BottomRightX = Rect(3);
BottomRightY = Rect(4);

imH = zeros(BottomRightY - TopLeftY + 1, BottomRightX - TopLeftX + 1, Clr); %  申请空间 
TransMat = inv(H) * inv(Trans);
for i = TopLeftY : BottomRightY  %  按行扫描
    for j = TopLeftX : BottomRightX
        ui = TransMat * [j; i; 1];  % ui为原始图像平面中的点
        ui = round(ui);  %  取整
        if (ui(2) >= 1) && (ui(2) <= Hight) && (ui(1) >= 1) && (ui(1) <= Width)  % ui在原始图像平面范围内
            try
                imH(i, j, :) = im(ui(2), ui(1), :);   %  将 ui对应的像素拷贝到点(j, i)的位置
            catch
                disp([i, j]);
            end;   
        end;
    end;
end;
imH = uint8(imH);  %  将图像变为24位位图