Blank Pixels When Correcting Wide-Angle (Barrel-Distortion) Images Using Polynomial Model

Hi everyone.

I have implemented a barrel-distotion correction algorithm (currently only experimenting with one radial-distortion coefficient, k1) to correct an image taken from a wide-angle (fish-eye) lens. However, the resulting corrected picture appears to have a pattern of blank pixels. I assume this is because the pixels get translated outwards rather than inwards. I can not seem to find any resource that faces this problem hence I am asking here. If anyone has faced this issue before, kindly explain where my algorithm is flawed. Comments in my code may help in understanding the steps I am taking. Original distorted photo and corrected photo are also shown to clarify what I mean. Thank you in advance.

Original Photo (3840x2160):

Resulting Corrected photo(7640x4480):

% imshow(a(:,:,1));
param =[0.0000002 0 0 0 0];
%param =[0.0000004 0 0 0 0];
%param = [0.0000002 -0.0000000000002 0 0 0];
%param = [0.0000002 -0.0000000000001 0 0 0];
k1 = param(1);
k2 = param(2);
k3 = param(3);
p1 = param(4);
p2 = param(5);

%get image and dimensions and rgb values
pic_in = imread('city.jpg');
[yres, xres, ch] = size(pic_in);

%xtrax and xtray are extra pixels that enlarge the corrected image 
%frame such that some of the corrected pixel locations that exceed 
%original picture dimensions are still visible
xtrax = 3800;
xtray = 2320;

corrected_pic = zeros( yres + xtray, xres + xtrax, ch);

%undistort
tic;
for yi = 1:1:yres
    for xi = 1:1:xres
        %get pixel value
        pixel = pic_in(yi,xi,:);

        %convert matrix index to coordinate
        x1 = cx(xi, xres);
        y1 = cy(yi, yres);

        %get coordinate of distortion
        x2 = floor(correct_x(x1,y1,k1,k2,k3,p1,p2));
        y2 = floor(correct_y(x1,y1,k1,k2,k3,p1,p2));

        %convert coordinate into matrix index
        x3 = mx(x2, xres, xtrax);
        y3 = my(y2, yres, xtray);

        %place pixel in new distorted coordinate
        if( x3 > 0 && x3 <= xres+xtrax && y3 > 0 && y3 <= yres+xtray )
            corrected_pic(y3,x3,:) = pixel;
        end
    end
end
toc;
figure();
corrected_pic = uint8(corrected_pic);
image(corrected_pic);
title('corrected');

figure();
imshow(pic_in);
title('distorted input');

%% calculate corrected coordinates based on lens parameters and actual coordinate
function x2 = correct_x(x1,y1,k1,k2,k3,p1,p2)
    x2 =  (x1*(1+(k1*((radius(x1,y1))^2))+(k2*((radius(x1,y1))^4))))+(2*p1*x1*y1)+(p2*(((radius(x1,y1))^2)+ 2*(x1^2)));
end

function y2 = correct_y(x1,y1,k1,k2,k3,p1,p2)
    y2 =  (y1*(1+(k1*((radius(x1,y1))^2))+(k2*((radius(x1,y1))^4))))+(2*p2*x1*y1)+(p1*(((radius(x1,y1))^2)+ 2*(y1^2)));
end

%% converts matrix index x to x coordinate(image center as origin)
function x = cx(xi, xres)
    x = xi - (xres/2);
end

%% converts matrix index y to y coordinate(image center as origin)
function y = cy(yi, yres)
    y = -1*yi + (yres/2);
end
%% gets radius
function r = radius(x,y)
    r = abs(sqrt((x^2)+(y^2)));
end

%% converts x coordinate to matrix index x
function x = mx(xi, xres, xtrax)
    x = xi + (xres/2) + xtrax/2;
end

%% converts y coordinate to matrix index y
function y = my(yi, yres, xtray)
    y = -1*yi + (yres/2) + xtray/2;
end

% %% finds extra frame space needed such that no pixel is outside frame after correction
% function xtrax = find_xtrax( xres, yres, k1, k2, k3, p1, p2)
%     xmax = xres/2;
%     ymax = yres/2;
%     xmax = correct_x(xmax,ymax,k1,k2,k3,p1,p2);
%     xtrax = ceil(xmax-(xres/2));
%     
% end
% 
% function xtray = find_xtray( xres, yres, k1, k2, k3, p1, p2)
%     xmax = xres/2;
%     ymax = yres/2;
%     ymax = correct_y(xmax,ymax,k1,k2,k3,p1,p2);
%     xtray = ceil(ymax-(yres/2));
% end

Answer 1

[Image Analyst]

You have holes in your image, because you're scanning the input image and deciding where the input pixel should go to in the output image. This is an intuitive, but wrong, approach that many beginners make. What happens is that you have some output pixels that never got an input pixel sent to them.

What you need to do is to scan the output image and decide where to pull the input pixel from. This will be a floating point number (between 4 input pixels) so you will have to use interpolation to figure out what the input pixel value would be at that location, since it most likely does not line up perfectly with an input pixel.

By the way, you might speed up the algorithm some by swapping the x and y loops. If y is the inside loop, you are scanning the memory in column major order (down rows in a particular column first, then move to to the next column) which is the way matrices are stored in memory in MATLAB.

[Taher Anjary]

Thank you kind sir/madam. I shall try this out and update the forum on its performance.

[Taher Anjary]

Hello everyone. By the advice given above by Image Analyst, I made necessary adjustments to my algorithm, code to which can be found below:

clear all;
clc;

 param =[-0.00000002 -0.000000000000001 0 0 0];
%param =[-0.00000002 -0.000000000000002 0 0 0];
%param =[-0.00000003 -0.000000000000001 0 0 0];
%param =[-0.00000002 -0.00000000000001 0 0 0];

k1 = param(1);
k2 = param(2);
k3 = param(3);
p1 = param(4);
p2 = param(5);

%get image and its dimensions and rgb values
pic_in = imread('city.jpg');
[yres, xres, ch] = size(pic_in);

%xtrax and xtray are extra pixels that enlarge the corrected image 
%frame such that some of the corrceted pixel locations that exceed 
%original picture dimensions are still visible
xtrax = 1800;
xtray = 2320;

corrected_pic = zeros( yres + xtray, xres + xtrax, ch);

%undistort
tic;
%scan through output photo frame-- backward
for xi = 1:1:xres+xtrax
    for yi = 1:1:yres+xtray
        
        %convert matrix index to coordinate
        [x1, y1] = cxy(xi, xres, xtrax, yi, yres, xtray);

        %get coordinate of pixel to extract based on distortion model
        [x2, y2] = forward_xy(x1,y1,k1,k2,k3,p1,p2);
        x2 = floor(x2);
        y2 = floor(y2);

        %convert above coordinate (x2,y2) into matrix index
        [x3,y3] = mxy(x2, xres, y2, yres);
        
        %get pixel value and place into corrected photo
        if( x3>0 && x3 <= xres && y3>0 && y3 <= yres)
            corrected_pic(yi,xi,:) = pic_in(y3,x3,:);
        end
    end
end
toc;

figure();
corrected_pic = uint8(corrected_pic);
image(corrected_pic);
title('corrected');
xlabel(['k1: ' num2str(k1) '  k2: ' num2str(k2)]);

figure();
imshow(pic_in);
title('distorted input');

%% FUNCTIONS

%% calculate corrected coordinates based on lens parameters and actual coordinate
function [x2, y2] = backward_xy(x1,y1,k1,k2,k3,p1,p2)
    r2 = (x1^2) + (y1^2);   %r2 means r squared
    x2 = x1 - x1*( k1*r2 + k2*(r2^2) + (k1^2)*(r2^2) + (k2^2)*(r2^4) + (2*k1*k2)*(r2^3) )/( 1 + 4*k1*r2 + 6*k2*(r2^2) );
    y2 = y1 - y1*( k1*r2 + k2*(r2^2) + (k1^2)*(r2^2) + (k2^2)*(r2^4) + (2*k1*k2)*(r2^3) )/( 1 + 4*k1*r2 + 6*k2*(r2^2) );
end

%% calculate corrected coordinates based on lens parameters and actual coordinate
function [x2, y2] = forward_xy(x1,y1,k1,k2,k3,p1,p2)
    r2 = (x1^2) + (y1^2);   %r2 means r squared
    x2 =  (x1*(1+(k1*r2)+(k2*(r2^2))))+(2*p1*x1*y1)+(p2*(r2+ 2*(x1^2)));
    y2 =  (y1*(1+(k1*r2)+(k2*(r2^2))))+(2*p2*x1*y1)+(p1*(r2+ 2*(y1^2)));
end

%% converts corrected output image's matrix index (xi,yi) to coordinate 
%  (x,y)(image center as origin)
function [x,y] = cxy(xi, xres, xtrax, yi, yres, xtray)
    x = xi - ((xres+xtrax)/2);
    y = -1*yi + ((yres+xtray)/2);
end

%% converts coordinate (xi,yi) to matrix index (x,y) of input image
function [x,y] = mxy(xi, xres, yi, yres)
    x = xi + (xres/2);
    y = -1*yi + (yres/2);
end

The approach works in that the holes have disappeared. This is because of a 'backward' approach in which you scan throuh an empty output image, and figure out which pixel from the input image to place there.

Note that the expression for forward and backward_xy() is the estimated inverse of forward_xy() and they can be used interchangeably depending on which research paper you use as reference.

I have however run into a different set of problems that I would like help with. I'll explain it slowly:

If I use :

param =[-0.00000002 -0.000000000000001 0 0 0];

... the resulting output image looks like this:

... which is great. I have all the original information intact. But the horizon and the buildings are still curved, therefore, I increase the magnitude of the parameters k1 and k2.

But then this starts to happen:

...when using...

param =[-0.00000002 -0.000000000000002 0 0 0];
param =[-0.00000003 -0.000000000000001 0 0 0];
param =[-0.00000002 -0.00000000000001 0 0 0];

... respectively.

As you may observe, I begin losing information at the corners, and the remaining portion of the corner begins mirroring, giving these weird shapes that I do not like.

So ultimately, my question is, how do I correct the distortion perfectly, without losing information at the corners and without getting these weird shapes? Is it about finding the right parameters k1 and k2? Is it impossible to get what I want? And how is this mirroring even possible?