I am trying to undo the radial distortion from three known straight lines (in the image projected to arcs due to the radial distortion) as described eg. in http://cv.cs.nthu.edu.tw/upload/activities/19/28.pdf. However, the approach does not work for me: the radial distortion center that I get lies outside of the input image & the remapping of the coordinates does not produce the expected result.
I am using the following image of the contours (the long horizontal arcs should be made straight):
I use Fitzgibbon Division model for the distortion. I fit the circles into of three arcs that I was to use & verify that the radius and centers of the circles are indeed fitted correctly. I use Taubin fit for that (from https://github.com/PmagPy/PmagPy/blob/master/SPD/lib/lib_curvature.py).
Then, each circle should correspond to the following equation:
$$x_d^2 + y_d^2 + D x_d + E y_d + F = 0,$$
where ($x_d,y_d$) are distorted coordinates and: $$D=\frac{A}{C\lambda} - 2x_0, E = \frac{B}{C\lambda}-2y_0, F = x_0^2 + y_0^2 - \frac{A}{C\lambda}x_0 - \frac{B}{C\lambda}y_0 + \frac{1}{\lambda},$$ where ($x_0, y_0$) is the distortion center, $\lambda$ is the distortion parameter from the Division model, and $A, B, C$ come from the equation of the line (projected to an arc due to the distortion): $A \frac{x_d}{1 + \lambda r_d^2} + B \frac{y_d}{1 + \lambda r_d^2} + C = 0$ ($r_d = \sqrt{(x_d - x_0)^2 + (y_d - y_0)^2}$).
So, Taubin fit gives me center of each circle ($x_c, y_c$) and its radius $r_c$ instead of $D$, $E$ and $F$. I derive those from the center and radius:
$$D = -2x_c, E = -2y_c, F = x_c^2 + y_c^2 - r_c^2.$$
I do these for all three arcs/circles. Then, according to the articles, I can use them to compute the center of radial distortion by solving a system of linear equations:
$$(D_1 - D_2)x_0 + (E_1 - E_2)y_0 + (F_1 - F_2) = 0,$$ $$(D_1 - D_3)x_0 + (E_1 - E_3)y_0 + (F_1 - F_3) = 0.$$
From this system, I get ($x_0, y_0$) and then also $\lambda$ from: $$\lambda = \frac{1}{x_0^2 + y_0^2 + Dx_0 + Ey_0 + F}.$$
However, this approach gives me the distortion center outside of my image ($x_0$ seems to make sense and is around the image center, $y_0$ is however 20 times higher than the height of the image). Because of that, the resulting $\lambda$ is very small and the final transformation basically just enlarges the image twice (not sure why) but does not really undo any distortion.
What is wrong in my approach? What am I missing? Thank you.