Distortion Correction without any camera information

I have a image which has pincushion distortion. I know nothing about the camera parameters, e.g. focal length, projection, etc.

I can, using SVD (singular value decomposition), successfully compute the coefficients of the polynomial that correct the image at the point centers, since I know the relative positions of the ideal points - point correspondences between points.

So, the polynomial is of the form

x = a0 + a1x + a2y + a3xy + a4x^2 + a5y^2;

y = b0 + b1x + b2y + b3xy + b4x^2 + b5y^2;

There are 6 different sets of of 12 coefficients, on set for each line of dots.

When I apply the coefficients with polynomials to the center of ground truth points, they map perfectly to the center of the distorted points.

So far so good.

My question is, what is the way to get the coefficients "between" the dots, to map the points between the dots to the ideal points?

My idea was to interpolate weights between each row:

weightA = 1 - n weightB = n

where n is from 0 - 1.0.

So for instance the first row would have a weight of 1, with n= 0. and row 2 coefficients would have a weight of 1, as weightA would be 0.

e.g. weightA *a0 + weightB b0, for all coeffients.

And the middle, would be half of each coefficient summed.

middle coefficients = 0.5*a0 + 0.5b0...0.5a5+0.5b5

So a linear weighting between the two.

This however, in my tests seems to give me the wrong results.

So, my question is, is the approach correct and my calculations incorrect, or is the approach for interpolating / weighting the coefficients wrong?

Thanks

Answer 1

[Bjorn Gustavsson]

To my understanding you are supposed to fit the polynomials to the full set of points, not row-by-row but all 36 points. That should give you 2 correction-polynomials to apply for each image coordinate.

HTH

[Richard]

thanks for answering. i am using the following function to make the measurement matrix

function [M] = makeMeasurementMatrix(app,points)
            [rows,~] = size(points);
            
            col1 = ones(rows,1);
            M=[col1(:,1),points(:,1),points(:,2),points(:,1).* points(:,2),points(:,1).^2,points(:,2).^2];
            end

and the following to get the coefficients

        function [aX,aY] = computeCoefficients(M,points)
            [U,S,V] = svd(M);
            xPattern = points(:,1);
            yPattern = points(:,2);
            % pseudo-inverse of diagonal matrix s
            sigma = eps^(1/2); % minimum value considered non-zero
            qqs = diag(S);
            qqs(abs(qqs)>=sigma) = 1./qqs(abs(qqs)>=sigma);
            qqs(abs(qqs)<sigma)=0;
            qqs = diag(qqs);
            aX = V*qqs'*U'*xPattern;
            aY = V*qqs'*U'*yPattern;
        end

i guess whats throwing me off is, if the measurement matrix has 36 measurements

1, xi, yi, xiyi, xi^2 yi2

...

1,xn,yn,xnyn,xn^2,yn^2

how does one solve for only 6 pairs of coefficients?

[Bjorn Gustavsson]

Perhaps it will be as simple as using polyfit. That way you will get a least-square fitting polynomial to your points. Typically it is desireable to have an overdetermined problem when doing fittings like this since that will (at least in general) reduce the impact of measurement noise of your results (the noise-contributions from multiple points will effect the results in different directions, while if you have an even-determined problem the noise will propagate through. And there is always noise in some general sense, your determination of the point centroids at least comes with some limited accuracy.)

If you use polyfit and polyval, those functions have some additional capabilities like scaling and centering of your data, which can be useful.

Finally, shouldn't the distortion-correction be done with a polynomial in the radial direction from the image-coordinate of the optical axis?

[Richard]

thank you for for your answer.

the data is being prescaled to be centered and mean distance of 2^1/2 from center to minimize the effects raising the pixel coordinate to the nth power, in this case n is 2.

the svd does the fit.

these are x ray calibration patterns and algorithms as described in this lecture @10 minutes in

Lecture

the algorithm is well explained, at least well enough, except how to get only 2 sets of coefficients as you mention, and he does too in the lecture, is at this time not understood by me.

thanks

[Richard]

i was just pointed to

Domingo Montgomery git hub, and his XVis package has similar code to mine, and i see i was setting up my matrix incorrectly - your first comment about getting only 1 set of pairs of cefficients was correct, of course, and i think i see where inset things up incorrectly.

thanks for your advice.