Re: Radial distortion model in PanoTools

[lireful]

I managed to understand the code from math.c (this is what I call unreadable code), but all the other questions are still open.

понедельник, 15 октября 2012 г., 19:15:37 UTC+6 пользователь lireful написал:

Hello,

I've read a lot on the PanoTools lens correction model in this user group and over the Internet, however I still have a lot of questions. All its disadvantages were discussed several times, but I could not a detailed description of the mechanism itself.

I have calibrated a fisheye lens using Kannala model, which is basically a polynomial R(theta), where R is the polar coordinates radius of current point. Now I have a need to work with this calibration in PTGui. We use Canon APS-C matrix crop-factor 1.62, 4272x2848 matrix in the landscape mode, 8 mm Sigma fisheyes.

I need to express PanoTools model as a function R(theta). To accomplish this, I can iterate all possible thetas in the [0, pi/2] range and find the radius R, corresponding to it. This is how I do it:

1) calculate undistorted equidistant radius (used both by Panotools and PTGui, not by Hugin, right?) R_undist = f*theta*m  as undistorted radius in pixels. I take m as a mm-to-pixel conversion factor from camera specification: matrix width is 22.2 mm = 4272 pixels.

2) convert it to PTGui units using the condition r=1 when r=width/2 (this is taken from PTGui documentation). r = R_undist  / (width/2).

3) solve quartic equation in PTGui units: a*r^4 + b*r^3 + c*r^2 + (1-a-b-c)*r - Rundist = 0 , which gives us r - the distorted radius

4) convert distorted radius to pixels: r_dist_px = r_dist*(width/2), so we have a r(theta) mapping.

For backprojection theta(R) we have to do the same it the opposite direction:

1) Convert R to PTGui units: r = R/(width/2)

2) Apply the undistorting polynomial: r_undist = a*r^4+b*r^3+c^r*2+(1-a-b-c)*r

3) Convert r_undist back to pixels, and then to mm : r_undist_mm = r_undist*(width/2)/m

4) Taking f*theta model into account, return r_undist_mm / f as a backprojection angle theta.

This looks logical, but doesn't work. The radiuses don't get distorted and undistorted as expected. I've also found this intriguing code in math.c of libpano13. Why are the coefficients used with r, and not polynomially? The same thing can be found in "inv_radial" function, where the polynomial is solved. Maybe I misunderstand the whole model?

int radial( double x_dest, double y_dest, double* x_src, double* y_src, void* params)
{
// params: double coefficients[4], scale, correction_radius
register double r, scale;
r = (sqrt( x_dest*x_dest + y_dest*y_dest )) / ((double*)params)[4]; // r = r_px/scale
if( r < ((double*)params)[5] )
{
scale = ((((double*)params)[3] * r + ((double*)params)[2]) * r + 
((double*)params)[1]) * r + ((double*)params)[0];                // scale = a*r+b*r+c*r+d
}
else
scale = 1000.0;

*x_src = x_dest * scale  ;                                            
*y_src = y_dest * scale  ;
    return 1;
}

I've also tried to understand the conversion from fov to f and back. This is what Joost, the author of PTGui, exactly told me:

Yes PTGui still uses the f*theta projection.
This is how PTGui calculates the horz fov:
        hfov=4*RAD2DEG*asin(cropwidth_mm/4.0/focallength);
cropwidth_mm is the diameter of the cropping circle, or the width of the image.
I realize that this is not consistent with the equidistant formula but it does match the fov for practical lenses.

However, when I use the diameter of crop circle as cropwidth_mm (4500 px in my case), it doesn't work correctly. I empirically found, that the mistake is linear, and the correct conversion formula is  f = 0.967555562*D/(4*m*sin(radians(fov)/4)) or hfov=4*asin(0.967555562*D/(4*m*focallength)), where D is the crop width in pixels, m the scale coefficient desribed earlier. The fact that the mistake is linear makes me think I use a different value of m? Or maybe it's more complicated.

I've contacted the author of PTGui himself, he has answered some of my questions, but was not willing to solve the problems, so I decided to try to get advice in this group. 

Thanks in advance for any help, 

Lire

[Erik Krause]

Am 15.10.2012 15:15, schrieb lireful:

> I have calibrated a fisheye lens using Kannala model, which is
> basically a polynomial R(theta), where R is the polar coordinates
> radius of current point. Now I have a need to work with this
> calibration in PTGui. We use Canon APS-C matrix crop-factor 1.62,
> 4272x2848 matrix in the landscape mode, 8 mm Sigma fisheyes.

You should have asked on the PTGui mailing list first of all...
Sorry, I don't understand the math you provided, only one point: r = 1
in panotools is not half the width, it's half the shorter side. And
perhaps its worth mentioning that the original panotools model (which is
used by PTGui) uses the abc parameters to handle the deviation from the
R=f*theta*m model, not the fisheye distortion itself.

However, I don't understand what you want to accomplish. If you simply
want to calibrate your lens for taking panoramas, do exactly this: take
a panorama (no zenith or nadir for that task), least-parallax-point set
up the best that is possible, no near detail, with a lot of overlap and
enough structure throughout the whole image, put it in hugin or PTGui,
create control points and optimize. That's all that is needed.

> Maybe I misunderstand the whole model?

I guess you read http://wiki.panotools.org/Lens_correction_model
To see the model in action while you play with the parameters you can
download abc.xls from http://4pi.org/downloads/

best regards
--
Erik Krause
http://www.erik-krause.de

[Bruno Postle]

On Mon 15-Oct-2012 at 21:01 +0200, Erik Krause wrote:
>Am 15.10.2012 15:15, schrieb lireful:
>
>>I have calibrated a fisheye lens using Kannala model, which is
>>basically a polynomial R(theta), where R is the polar coordinates
>>radius of current point. Now I have a need to work with this
>>calibration in PTGui. We use Canon APS-C matrix crop-factor 1.62,
>>4272x2848 matrix in the landscape mode, 8 mm Sigma fisheyes.
>
>You should have asked on the PTGui mailing list first of all...
>Sorry, I don't understand the math you provided, only one point: r =
>1 in panotools is not half the width, it's half the shorter side. And
>perhaps its worth mentioning that the original panotools model (which
>is used by PTGui) uses the abc parameters to handle the deviation
>from the R=f*theta*m model, not the fisheye distortion itself.

Also PTGui handles photo cropping differently: Hugin fits the radial
distortion and angle of view to the edge of the image frame,
cropping is simply a mask that excludes pixels from the stitch.
PTGui fits the radial distortion to the boundary of the crop, the
crop boundary actually expands or contracts the 'photo' - This means
that for circular fisheyes (which are usually cropped) it isn't
possible to transfer lens parameters between the tools.

--
Bruno

[Carlos Eduardo G. Carvalho (Cartola)]

I made myself the same question: what do you want to accomplish with all that math? Maybe we can help you to achieve your goal in a simpler way (I also didn't understand the math).

I once have calibrated a fish eye lens with the hugin tutorial. It took me some hours to make the pictures, understand the concepts and achieve a good result. Later I have discovered that I didn't need to do that to make a full sphere panorama. Hugin - and probably PTGui - can do the necessary distortions based on control points and in a defined angle - like the 360 of the full sphere.

Here is the hugin tutorial I mentioned: http://hugin.sourceforge.net/tutorials/calibration/en.shtml

Bests,

Carlos E G Carvalho (Cartola)
http://cartola.org/360
http://www.panoforum.com.br/

2012/10/15 Erik Krause <erik....@gmx.de>

--
You received this message because you are subscribed to the Google Groups "Hugin and other free panoramic software" group.
A list of frequently asked questions is available at: http://wiki.panotools.org/Hugin_FAQ
To post to this group, send email to hugi...@googlegroups.com
To unsubscribe from this group, send email to hugin-ptx+unsubscribe@googlegroups.com
For more options, visit this group at http://groups.google.com/group/hugin-ptx

[lireful]

Thanks for the answers.

I do understand the differences of PTGui vs Hugin and PanoTools. I deliberately decided to ask the question here, because I've seen that here there are the people who "understand the maths".

What I'm asking for is a part of my much larger calibration project, where using the calibration in PTGui is only one of the steps. I do know about other calibration options, and this is not about that.

To state it in short, I really need to understand how the projection and backprojection of the rays work in Panotools/Hugin/PTGui fish/lens correction model. Not the calibration, not the exact piece of software. The model.

Hope one of the developers would come across the message and share the knowledge :)

[Bruno Postle]

On Mon 15-Oct-2012 at 14:54 -0700, lireful wrote:
>
>To state it in short, I really need to understand how the projection and
>backprojection of the rays work in Panotools/Hugin/PTGui fish/lens
>correction model. Not the calibration, not the exact piece of software. The
>model.

I'm not sure what you want.

In the fisheye model the pixel distance from the centre of the image
is directly proportional to the angle from the view axis.

The lens correction model is purely a radial distortion of the
bitmap image, unrelated to the fisheye/rectilinear/etc... projection
and unrelated to any view angles.

This lens correction is normalised to the width or height of the
image whichever is narrowest.

--
Bruno

[Erik Krause]

Am 15.10.2012 23:54, schrieb lireful:
> To state it in short, I really need to understand how the projection and
> backprojection of the rays work in Panotools/Hugin/PTGui fish/lens
> correction model. Not the calibration, not the exact piece of software.
> The model.
>
> Hope one of the developers would come across the message and share the
> knowledge :)

The original model was implemented by Helmut Dersch about 13 years ago.
Unfortunately he is not active on this list. Occasionally he writes (and
presumably reads) the PanotoolsNG list.

The model has not changed since then, it only has been extended. You
should find all information you need on the panotools wiki, on the
mirror of the old Dersch page at http://www.panotools.org/dersch/ and
eventually in the panotools mailing list archive which covers the very
first proj-imim list as well: http://www.panotools.org/mailarchive.php

[lireful]

Hello,

I would like to report that I've localized and solved all my problems. Problem with fov-focal length correction was due to the use of individual crop parameters for each image. In this case PTGui seems to use a default crop diameter Dcrop for this conversion. As far as I understood, it's the value which binds theoretical fov to theoretical focal length.

Moreover, I found that my problems with radial correction were connected to the incorrect fov estimation, as it is only the PTGui interface which calculates f from fov as stated above, but for radial distortion correction one should use the value calculated from equidistant projection.

Also, I found out that the polynomial in the forward direction should be used in projection, and solving the polynomial equation - in backprojection procedure.

In all the rest my statements were OK.

I'll just leave it here for next generation of developers :)

[Erik Krause]

Am 19.10.2012 17:05, schrieb lireful:
> Moreover, I found that my problems with radial correction were connected
> to the incorrect fov estimation, as it is only the PTGui interface which
> calculates f from fov as stated above, but for radial distortion
> correction one should use the value calculated from equidistant projection.

You should indeed use the same crop circle diameter for all images in
PTGui unless you specify individual lens parameters for the differing
diameters. Not so in hugin, where all correction parameters always refer
to the image width.

Both ways of cropping are found in the underlying panotools model
though, as described on
http://wiki.panotools.org/Panotools_internals#Cropping