Activity 13 - Correcting Geometric Distortion

In this activity, we correct an image with geometric distortion. We used an image that has grid patterns in the image. This is to determine the kind of distortion present in the image and correct the image by making curve lines make straight. There are two kinds of distortions: 1. pincushion and 2. barrel distortion. Barrel distortion introduces bulging effect in the images while pincushion distortion introduces inward bending of lines. To correct the image, a regularly occurring pattern is required to be part of the image to be able to visualize the distortion effect.

In the correcting process, a transformation of the pixel coordinates and interpolation of gray level values are required to complete the correction. Pixel coordinate transformation allows one to straighten the coordinates and gray level interpolation to visualize the transformation. It can be expressed in matrix notation,

and solving for the coefficients,

.

After which, we can now compute for the undistorted pixel coordinates using:

.

Then, we project the grayscale values from the original image to the undistorted using:

where the coefficients can be solved similar to the matrix notation above.

By generating an ideal grid of known vertex coordinates, we can proceed to pixel coordinates transformation.

Shown to the left is the original image with barrel distortion. At the right is the ideal rectangle vertices and below them is the corrected image. Note that there are some pixels that are not interpolated in the corrected image.

I give myself 9 pts for understanding the activity but blogging it late.

Reference:
A13 - Correcting Geometric distortion. M. Soriano. 2008.

Acknowledgment:
I acknowledge the help of Jaya in understanding the notations in the procedure.

Code:
chdir('E:\Documents and Settings\vergara\My Documents\acads\1st sem 09-10\186\act13');
//165, 92 //182, 92 //165, 109 //182, 109
stacksize(50000000)
img = gray_imread('barrel1.png');

// for visualization of ideal rect
sx = size(img,1); sy = size(img,2);
z = zeros(sx, sy);
X = 1:110:sx; Y = 1:110:sy; //17
//z(X, Y) = 1;
//imwrite(z, 'ideal_rect.jpg');
//scf(), imshow(z, [])


// locate real vertices
pts = 25; //25, 49, 81
scf(), imshow(img, [])
l = locate(pts); l = l';
l = [l(pts:-1:1,2) l(pts:-1:1,1)];

fprintfMat('barrel.txt', l, '%0.2f');
l = fscanfMat('barrel.txt');

// for visualization of real rect
z1 = zeros(sx, sy);
for i = 1:pts
z1(l(i,1), l(i,2)) = 1;
end
scf(), imshow(z1, [])
//imwrite(z1, 'vertices_rect.jpg');

// reshaping of data
ideal = [X' mtlb_repmat(Y(size(Y,2))',size(Y,2),1)];
rreal = l;
for i = size(X,2) - 1:-1:1
ideal = cat(1,ideal,[X' mtlb_repmat(Y(i)',size(X,2),1)])
end

// test
for j = 1:4// row
img(ix(j,1), iy(j,1)) = 1;
end

//reshape again
rx = matrix(rreal(:,1),size(X,2),size(Y,2)); rx = rx(:,size(Y,2):-1:1);
ry = matrix(rreal(:,2),size(X,2),size(Y,2)); ry = ry(:,size(Y,2):-1:1);

ix = matrix(ideal(:,1),size(X,2),size(Y,2)); ix = ix(:,size(Y,2):-1:1);
iy = matrix(ideal(:,2),size(X,2),size(Y,2)); iy = iy(:,size(Y,2):-1:1);

//
undist = zeros(sx, sy);
for i = 1:size(X,2) - 1
for j = 1:size(Y,2) - 1
T = [ix(i,j) iy(i,j) ix(i,j)*iy(i,j) 1;...
ix(i,j+1) iy(i,j+1) ix(i,j+1)*iy(i,j+1) 1;...
ix(i+1,j) iy(i+1,j) ix(i+1,j)*iy(i+1,j) 1;...
ix(i+1,j+1) iy(i+1,j+1) ix(i+1,j+1)*iy(i+1,j+1) 1];
Xx = [rx(i,j); rx(i,j+1); rx(i+1,j); rx(i+1,j+1)];
Yy = [ry(i,j); ry(i,j+1); ry(i+1,j); ry(i+1,j+1)];
C14 = inv(T)*Xx;
C58 = inv(T)*Yy;

for kx = ix(i,j):ix(i+1,j) - 1
for ky = iy(i,j):iy(i,j+1) - 1
xh = C14(1)*kx + C14(2)*ky + C14(3)*kx*ky + C14(4);
yh = C58(1)*kx + C58(2)*ky + C58(3)*kx*ky + C58(4);
// not this easy, solve for the 4 unknowns and 4 equations
if modulo(xh, round(xh)) == 0 & modulo(yh, round(yh)) == 0 // integer coordinates
undist(round(xh), round(yh)) = img(round(xh), round(yh));
else // non-integer coordinates
if xh == 1 | xh == sx | yh == 1 | yh == sy // at boundaries
undist(xh,yh) = 0;
else // not at boundaries
U = [xh-1 yh-1 (xh-1)*(yh-1) 1;...
xh-1 yh+1 (xh-1)*(yh+1) 1;...
xh+1 yh-1 (xh+1)*(yh-1) 1;...
xh+1 yh+1 (xh+1)*(yh+1) 1];
V = [img(xh-1, yh-1) img(xh-1, yh+1) img(xh+1, yh-1) img(xh+1, yh+1)]';
A = inv(U)*V;
undist(xh,yh) = A(1)*xh + A(2)*yh + A(3)*xh*yh + A(4);
end
end
end
end
end
end

scf(), imshow(undist, []);