Image rectification algorithm in Matlab

Question

I have recently found an interesting article regarding image rectification for two stereo image pairs. I liked the algorithm because it was very compact and from what the article suggested it did the right thing. After I implemented the matlab version on two images, I didn't get a correct rectified image. I got an image that was pitch black apart from the left and down line which had pixels. In the image there also were some gray pixels from the original image but just a hand full. I posted below the matlab code, and the link to the article and also an example of the result I got for one image (for the other image it was the same)

This is the link to the article A compact algorithm for rectification of stereo pairs.

A screen shot with the initial images and the results is bellow:

The initial images are the following two(such that you do not have to search for another stereo pair) :

function [T1,T2,Pn1,Pn2] = rectify(Po1,Po2)

% RECTIFY: compute rectification matrices 

% factorize old PPMs
[A1,R1,t1] = art(Po1);
[A2,R2,t2] = art(Po2);

 % optical centers (unchanged)
c1 = - inv(Po1(:,1:3))*Po1(:,4);
c2 = - inv(Po2(:,1:3))*Po2(:,4);

% new x axis (= direction of the baseline)
 v1 = (c1-c2);
% new y axes (orthogonal to new x and old z)
v2 = cross(R1(3,:)',v1);
% new z axes (orthogonal to baseline and y)
v3 = cross(v1,v2);

% new extrinsic parameters 
R = [v1'/norm(v1)
   v2'/norm(v2)
   v3'/norm(v3)];
% translation is left unchanged

% new intrinsic parameters (arbitrary) 
A = (A1 + A2)./2;
A(1,2)=0; % no skew
A(1,3) = A(1,3) + 160;
% new projection matrices
Pn1 = A * [R -R*c1 ];
Pn2 = A * [R -R*c2 ];

% rectifying image transformation
T1 = Pn1(1:3,1:3)* inv(Po1(1:3,1:3));
T2 = Pn2(1:3,1:3)* inv(Po2(1:3,1:3));

function [A,R,t] = art(P)
% ART: factorize a PPM as  P=A*[R;t]
Q = inv(P(1:3, 1:3));
[U,B] = qr(Q);

R = inv(U);
t = B*P(1:3,4);
A = inv(B);
A = A ./A(3,3);

This is the "main" code from which I call my rectify function

img1 = imread('D:\imag1.png');
img2 = imread('D:\imag2.png');

im1 = rgb2gray(img1);
im2 = rgb2gray(img2);

im1 = im2double(im1);
im2 = im2double(im2);

figure; imshow(im1, 'border', 'tight')
figure; imshow(im2, 'border', 'tight')

%pair projection matrices obtained after the calibration P01,P02

a = double(9.765*(10^2))
b = double(5.790*(10^-1))
format bank;
Po1 = double([a 5.382*10 -2.398*(10^2) 3.875*(10^5); 
    9.849*10 9.333*(10^2) 1.574*(10^2) 2.428*(10^5);
    b 1.108*(10^(-1)) 8.077*(10^(-1)) 1.118*(10^3)]);
Po2 = [9.767*(10^2) 5.376*10 -2.400*(10^2) 4.003*(10^4);
    9.868*10 9.310*(10^2) 1.567*(10^2) 2.517*(10^5);
    5.766*(10^(-1)) 1.141*(10^(-1)) 8.089*(10^(-1)) 1.174*(10^3)];
[T1, T2, Pn1, Pn2] = rectify(Po1, Po2);

imnoua =  conv2(im1, T1);
imnoua2 = conv2(im2, T2);

fprintf('Imaginea noua e \n');

figure; imshow(imnoua, 'border', 'tight')
figure; imshow(imnoua2, 'border', 'tight')

Thank you for your time!

Answer 1

As Shai says, T1 and T2 are projective transformation matrices, not filter kernels. You should be using imwarp, rather than conv2:

imnoua = imwarp(im1, projective2d(T1));
imnoua2 = imwarp(im2, projective2d(T2));

Better yet, use rectifyStereoImages from the Computer Vision System Toolbox. Check out this example.