Lens Correction in PanoTools

This article is about a possible new way to correct for lens distortions in panotools library. If you look for a general description of the current way panotools deals with lens distortion see Lens correction model

Lens Correction in PanoTools

The PanoTools library implements an effective, but rather idiosyncratic method for correcting lens projections, that causes a good deal of puzzlement. Lens parameters optimized for one image format generally do not work for a different format; even rotating a set of images 90 degrees before aligning them produces different and incompatible lens parameters. One would expect that there must be a way to convert either of those parameter sets to a common form, that would apply equally well to both formats, or indeed to any image taken with the same lens. To see how that might be done, I have made a detailed analysis of PanoTools lens correction computations, based on the code in historic as well as current versions of libpano and helpful discussions with Helmut Dersch.

Why Lens Correction?

To make a panoramic image from photographs, it is essential to be able to calculate the direction in space corresponding to any given position in a given photo. Specifically, we need to know the angles between the view directions of the photos (the alignment of the images), and a radial projection function that relates the distance of a point from image center to the true angle of view, measured from the optical axis of the lens. Given a set of control points linking the images, PanoTools estimates both the alignment and the lens projection by a nonlinear least squares fitting procedure -- optimization. Using the fitted lens parameters, the stitcher can correct each image to match the ideal geometry of the scene, according to whatever projection is chosen for the panorama. Done right, that makes all the images fit together perfectly; moreover, it yields a panoramic image that seems to have been made with a perfect lens.

Mapping View Angle <=> Radius

The radial projection curve of a real lens may approximate some known mathematical function, but in practice it must be determined experimentally, a process known as calibrating the lens. A calibration is a parametrized mathematical model, fitted to experimental data. The typical model consists of an ideal angle-to-radius function, and a polynomial that converts the ideal radius to the actual radius measured on the image.

Like many lens calibration programs, libpano uses just two ideal functions to model lenses: rectilinear, for 'normal' lenses, and 'fisheye', for all others. The rectilinear projection has radius proportional to the tangent of the view angle. PT's 'fisheye', better known as the equal-angle spherical projection, has radius proportional to the angle itself. The constant of proportionality is the lens focal length, F. With angle A in radians, and R the ideal radius, the formulas are

Rectilinear: ${\frac {R}{F}}=\tan(A)$

Equal-angle: ${\frac {R}{F}}=A$

Of course R and F have to be measured in the same units. If we have F in mm, then R is in mm also. If we want to measure R in pixels, then we need F in pixels; for example

$F\ (pixels)=FL\ in\ mm*{\frac {image\ width\ in\ pixels}{sensor\ width\ in\ mm}}$.

In any case, F is the constant of proportionality between the actual radius and the value of a trigonometric function that defines the basic shape of the projection.

In physical optics, focal length is defined as the first derivative of R by A, at A = 0. That is easy to see if we write $R=FA$ or $R=F\tan(A)$, because the slopes of A and tan(A) are both 1 at A = 0. This is also true of other trigonometric functions commonly used as ideal lens projections:

Equal-Area: ${\frac {R}{F}}=2\sin \left({\frac {A}{2}}\right)$

Stereographic: ${\frac {R}{F}}=2\tan \left({\frac {A}{2}}\right)$.

The dimensionless quantity $N={\frac {R}{F}}$ is the normalized ideal radius. Multiplying N by the focal length, in any units, gives the ideal image radius in the same units.

Generic Correction Scheme

The difference between the real lens projection and the ideal one is modeled by an adjustable correction function that gives the observed radius as a function of the ideal radius. The adjustable part is almost always a polynomial, because it it easy to fit polynomials to experimental data. The argument to the polynomial is usually the normalized ideal radius,

$N={\frac {R}{F}}$,

because that makes the polynomial coefficients independent of how image size is measured. The constant term is 0 because both radii are zero at the same point. If the coefficient of the linear term is 1, so that the first derivative at 0 is 1, then the value of the polynomial will be the normalized observed radius, n = r / F. Multiplying n by the focal length, in any units, gives the observed image radius in the same units:

$r=Fn$.

Many calibration packages use a polynomial with only even order terms beyond the first:

$n=N+aN^{2}+bN^{4}+cN^{6}$.

Equivalently

$n=N(1+aN+bN^{3}+cN^{5})$

The expression in parentheses is the ratio of observed to ideal radius, which is expected to be close to 1 everywhere if the ideal model function is well chosen.

PanoTools Correction Scheme

Lens correction in PanoTools is unusual in several respects. First, it ignores the physical parameters of the lens (focal length) and camera (pixel size). Instead, it computes angle-to-radius scale factors from image dimensions and fields of view, as described below. All correction computations are in terms of image radii, measured in pixels, rather than the normalized radii described above. However, normalized radii are evaluated implicitly.

Second, the correction is computed in equal-angle spherical coordinates, rather than camera coordinates. Observed image points are found by remapping those coordinates according to the ideal lens projection, and rescaling them according to the ratio of pixel sizes in the source and ideal images.

Third, the correction is normalized to hold a certain radius, $r_{0}$, constant. It essentially consists of a cubic polynomial that computes the ratio of observed to ideal radius. The argument to this polynomial is ${\frac {R}{r_{0}}}$, and its constant term is set so that the result is exactly 1 when the argument is 1, that is, when $R=r_{0}$. With

$X={\frac {R}{r_{0}}}$

The correction factor is

$x=(1-a-b-c)+aX+bX^{2}+cX^{3}$,

and the observed radius is given by

$r=Rx$.

The observed radius is thus formally a 4th order polynomial in R:

$r=sR+tR^{2}+uR^{3}+vR^{4}$,

where $s=(1-a-b-c),\ t={\frac {a}{r_{0}}},\ u={\frac {b}{{r_{0}}^{2}}},\ v={\frac {c}{{r_{0}}^{3}}}$.

The normalization makes this relation stable under optimization of a, b, and c. It is also essential to the correctness of the result, which can be seen as follows. The ideal radius is

$R=FN$

where F is the ideal focal length in pixels. We can thus write the adjusted radius as

$r=FA\ \operatorname {poly} \left({\frac {FA}{r_{0}}}\right)$,

The normalized observed radius $n={\frac {r}{F}}$, so

$n=A\ \operatorname {poly} \left({\frac {FA}{r_{0}}}\right)$.

Since n is a dimensionless number, this can only be true if the argument to poly() is also a dimensionless number. F has the dimension of pixels, therefore it must be divided by another parameter dimensioned in pixels.

The overall computation proceeds as follows. PanoTools computes the ideal radius R by mapping a point in the panorama (which plays the role of the ideal image) to equal angle spherical projection. Then $R={\sqrt {h^{2}+v^{2}}}$, where h and v are the pixel coordinates relative to the center of the equal-angle projection. Then PT's radius() function computes x as described, and returns scaled coordinates ( h x, v x ). If the lens is rectilinear, PT next remaps those coordinates to rectilinear; if it is a fisheye, no remapping is needed. In either case the coordinates are finally rescaled to account for any difference in resolution between the panorama and the source image. The scale factor is computed from the dimensions and angular fields of view of the panorama and the source image, as follows.

$d={\frac {half\ width\ of\ pano}{\operatorname {A2Npano} \ {(half\ hfov\ of\ pano)}}}$,

$e={\frac {half\ width\ of\ source}{\operatorname {A2Nsource} \ {(half\ hfov\ of\ source)}}}$,

where A2Npano and A2Nsource are the ideal functions for panorama and lens. Then

$R_{source}=R_{pano}{\frac {e}{d}}$.

The scale factors d and e are focal lengths in pixels, because A2N() yields the normalized radius, equal to ${\frac {R}{F}}$. For the panorama, which follows an ideal projection, d is identical to F. In fact d, under the name “distance factor”, is used in many of libpano's coordinate transformation functions to convert radius in pixels to the ideal normalized radius in trigonometric units.

For the source image, whose true projection is only approximately known, e is an estimate of F according to the fitted correction parameters. Since hfov is one of those parameters, the fitted value of e will be proportional to the true F; the constant of proportionality will approach 1 as the fitted polynomial coefficients approach 0.

Portable Lens Parameters

The focal length in pixels must be known in order to calibrate a lens. This quantity clearly depends on both lens and camera properties. In most cases today, equipment manufacturers' specifications can provide the needed data:

$F={(FL\ in\ mm)}{\frac {image\ width\ in\ pixels}{sensor\ width\ in\ mm}}$.

Computing F this way makes it possible for the fitted correction coefficients to be independent of camera format, as explained above. In any practical calibration scheme F is actually an adjustable parameter. However the fitted value is expected to be close to the one implied by these physical specifications, the main uncertainty being how accurately the nominal lens focal length reflects the true one.

Focal length and projection function are separable lens properties. In fact many schemes determine F from different data than those used to fit the lens curve. PanoTools is somewhat unique in fitting all lens parameters to one set of experimental values.

The normalization of the radial polynomial makes it possible to equate the correction computed by PanoTools to the generic one. As discussed above, that is true for any choice of r0 that is proportional to image size. The best choice would be

$r_{0}=F_{pano}=d$, the distance parameter defined above.

That would make the argument of the radius scaling polynomial

${\frac {R_{pano}}{F_{pano}}}$

which is equal to the ideal normalized radius, N. Then the polynomial coefficients fitted by PanoTools would be independent of image dimensions, and could be used directly in a generic function, to compute the radius correction in camera coordinates. The value actually used is merely proportional to $F_{pano}$; the proportionality factor is a function of image dimensions, so the correction coefficients depend on the image format and are not portable.

To convert these PT lens parameters to a fully portable form requires two things: computing portable radial coefficients, and defining F in a way that is not tied to an image format. The latter requires access to at least one physical parameter, not in the PT parameter set. Ideally that would be the physical size of a pixel. With w the pixel width in mm,

$F(mm)=wF(pixels)=we$, scale factor e defined above.

This could be calculated by front-end software using the value of e from libpano.

The coefficients can be converted using data available inside libpano. With

$k={\frac {d}{r_{0}}}$,

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle w' = w k = (1 – a – b – c) * k} ,

$a'=ak^{2}$,

$b'=bk^{3}$,

$c'=ck^{4}$

are the coefficients of a polynomial in $N={\frac {R}{d}}$ that computes the same radius correction factor as the PT polynomial. The constant term w' is no longer a simple function of the other three, however it can be reduced to 1 by dividing all coefficients by w'. The reduced coefficients are

$W=1$

$A=a{\tfrac {k}{w}}$

$B=b{\tfrac {k^{2}}{w}}$

$C=c{\tfrac {k^{3}}{w}}$

So the portable radius mapping is

$r=R(1+AN+BN^{2}+CN^{3})$

Along with the ideal function A2Nsource(), which gives N as a function of angle, this constitutes a portable lens projection function.

This function could be very useful even without a portable version of PT's focal length estimate. Ultimately, focal length is just a linear scale factor, that can be fitted or estimated in many different ways according to the problem at hand, while the radial function represents the essential 'lens curve' that is independent of image scale.

-- 19 Jan 2010 T K Sharpless