I am trying to fit a simple model over 2D data points, in the frame of an image formation model with perspective and optical distortion.
My model is the ratio of a second degree polynomial over a first degree one (with unit constant term). The numerator is meant to represent distortion, and the denominator, perspective.
$$X=\frac{P(x,y)}{ax+by+1}\\Y=\frac{Q(x,y)}{ax+by+1}.$$
As the equations are easily linearized by rewriting $$P(x,y)-axX-byX=X\\Q(x,y)-axY-byY=Y,$$ it is tempting to solve for the 6+6+2 unknown coefficients by linear least-squares. I have on the order of $100$ data points, with little noise.
Unfortunately, when working on real-world cases, I am facing the problem of the vertical asymptote, when the denominator vanishes, $ax+by+1=0$. Even in the case of weak perspective (i.e. true $a$ and $b$ small), the fitting is such that this line passes well through the data points, making the model unusable.
I am conscious that linearizing the equations does change the weighting of the fitting errors, but I expected this effect to be neglectible when $ax+by+1\approx1$.
Can you suggest a way to cope with this ? I would prefer to keep the simplicity of the linear fit.
UPDATE:
In the meantime, I understood that a quadratic distortion model is inappropriate, as the deformation is quasi-symmetric with respect to the image center and behaves like an odd function of the coordinates. So the least squares machine does its best to fit a bad model on the data.
A cubic model could probably do, but it would introduce (too) many new DOFs.
For those willing to reproduce, here is a typical sample set $(x,y,X,Y)$:
147.986786,36.568508,20.000000,80.000000
92.672272,90.325333,30.000000,70.000000
200.918198,35.073734,10.000000,60.000000
144.500000,89.492027,20.000000,70.000000
89.631821,143.853409,30.000000,80.000000
86.550903,198.709274,40.000000,80.000000
255.325500,34.015659,10.000000,50.000000
83.629425,253.992264,50.000000,80.000000
197.982376,88.702644,20.000000,60.000000
81.311066,309.768890,60.000000,80.000000
79.301971,365.847931,70.000000,80.000000
141.760300,143.966385,30.000000,70.000000
139.109436,199.038620,40.000000,70.000000
310.748657,33.539936,10.000000,40.000000
367.114899,33.271278,10.000000,30.000000
195.724686,144.002106,30.000000,60.000000
136.774155,255.275208,50.000000,70.000000
130.931183,425.586029,80.000000,70.000000
193.464584,200.022919,40.000000,60.000000
134.386063,312.168579,60.000000,70.000000
253.010391,88.072762,20.000000,50.000000
132.415466,368.963928,70.000000,70.000000
424.283508,33.481441,10.000000,20.000000
251.211655,143.932510,30.000000,50.000000
309.129211,87.976601,20.000000,40.000000
191.345078,256.566406,50.000000,60.000000
366.304596,88.049095,20.000000,30.000000
185.931000,429.153992,80.000000,60.000000
189.223557,314.327332,60.000000,60.000000
423.759918,88.164680,20.000000,20.000000
307.969360,144.222336,30.000000,40.000000
249.376968,200.434052,40.000000,50.000000
187.411423,371.838593,70.000000,60.000000
481.585938,88.737305,20.000000,10.000000
365.532166,144.476608,30.000000,30.000000
242.521271,432.500977,80.000000,50.000000
247.423187,258.299530,50.000000,50.000000
245.589081,316.289276,60.000000,50.000000
243.970413,374.488556,70.000000,50.000000
306.321594,201.368225,40.000000,40.000000
364.377716,202.303116,40.000000,30.000000
482.124298,145.624298,30.000000,10.000000
424.092957,144.699524,30.000000,20.000000
304.931702,259.518250,50.000000,40.000000
300.205048,435.520050,80.000000,40.000000
301.970215,376.738373,70.000000,40.000000
363.443298,260.768585,50.000000,30.000000
303.207611,318.124176,60.000000,40.000000
482.015625,203.984375,40.000000,10.000000
423.061401,202.981674,40.000000,20.000000
481.908661,263.041107,50.000000,10.000000
361.925385,319.847137,60.000000,30.000000
481.071808,322.264557,60.000000,10.000000
359.055298,437.726196,80.000000,30.000000
360.749084,379.024353,70.000000,30.000000
479.758972,382.034943,70.000000,10.000000
418.567078,440.028625,80.000000,20.000000
422.423218,261.834808,50.000000,20.000000
420.091003,380.665466,70.000000,20.000000
421.500000,321.085266,60.000000,20.000000
This is not a complete answer to your question but it is too long for a comment.
The problem would be simple if parameters $a,b$ were known. So, we could define a function $$\Phi(a,b)=SSQ_X+SSQ_Y$$ and try to minimize it with respect to $a$ and $b$ (even computing a regular grid would not be very time consuming); the major part of the work involves linear regression.
What you could notice is that $(a=b=0)$ does not lead to too bad results. The quick and dirty calculations I did only reduced the objective function from $318$ to $237$ which is far to be impressive for two extra parameters.
I shall continue working this problem and let you know. May be, meanwhile, could you have a look to a grid or a contour plot of the function $\Phi(a,b)$ ?
Added later to this answer
I started new calculations "linearizing" the model. The solution is immediate, leading in particular to $a=-5.606922 10^{-3}$, $b=-3.436547 10^{-3}$. Injected in the full model, the initial value of $\Phi(a,b)$ is huge $(87463)$ but iterations reduce it down to $150$ and the corresponding crucial parameters are $a=-6.152448 10^{-3}$, $b=-2.538128 10^{-3}$. As I said earlier, $(a=b=0)$ lead to a sum of squares equal to $318$. The fitted model seems to be quite good.