View source for Lens Correction in PanoTools

This article is a mathematical analysis of how the panotools library computes lens correction parameters, why those parameters are not portable, and how they could be made portable.  For a more general and use-oriented 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.

==Ideal Lens Models==

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: 	<math>\textstyle \frac R F  =  \tan(A)</math>
:	Equal-angle:	<math>\textstyle \frac R F  =  A</math>
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.  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  <math>\textstyle R = F A</math>  or  <math>\textstyle R = F \tan(A)</math>, 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:	<math>\textstyle \frac R F  = 2\sin\left(\frac A 2 \right)</math>
:	Stereographic:	<math>\textstyle \frac R F  =  2\tan\left(\frac A 2 \right)</math>.

The dimensionless quantity <math>\textstyle N = \frac R F</math> 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 should be the normalized ideal radius,
:       <math>\textstyle N = \frac R F</math>, 
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:
:	<math>\textstyle r = F n</math>.

Many calibration packages use a  polynomial with only even order terms beyond the first:
:	<math>\textstyle n =  N  +  a N^2  +  b N^4  +  c N^6</math>.
Equivalently
:	<math>\textstyle n  =  N ( 1 + a N  +  b N^3  +  c N^5 )</math>
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 to the ideal lens projection, and rescaling them according to the ratio of pixel sizes in the source and ideal images.  

Third, the correction polynomial is normalized to hold a certain radius, <math>\textstyle r_0</math>, constant.  It is a cubic polynomial, that computes the ratio of observed to ideal radius.  Its argument is <math>\textstyle \frac R {r_0}</math>, and its constant term is set so that the result is exactly 1 when the argument is 1, that is, when <math>\textstyle R = r_0</math>.  With
:	<math>\textstyle X  =  \frac R {r_0}</math>
The correction factor is
:	<math>\textstyle x  =  (1 - a - b - c)  +  a X  +  b X^2  +  c X^3</math>,
and the observed radius is given by
:	<math>\textstyle r   =  R x</math>.
The observed radius is thus formally a 4th order polynomial in R:
:	<math>\textstyle r  =  s R  +  t R^2  +  u R^3  +  v R^4</math>,
where  <math>\textstyle s = (1-a-b-c),\  t = \frac a {r_0},\  u = \frac b {{r_0}^2},\  v = \frac c {{r_0}^3}</math>.

The normalization makes the PanoTools polynomial equivalent to the generic one, but with different coefficients.  This can be seen as follows.  The ideal radius is
:	<math>\textstyle R  =  F N</math>
where F is the ideal focal length in pixels, so we can write the adjusted radius as
:	<math>\textstyle r  =  F N\ \operatorname{poly}\left(\frac {F N} {r_0} \right)</math>,
If <math> r_0 </math> is proportional to F, then the quotient is proportional to N, and the polynomial is equivalent to one whose argument is N. That is the case when <math> r_0 </math> is proportional to source image size, which is proportional to F by definition.  But the proportionality factor varies with source image format, so the PanoTools coefficients also depend on source format.
	
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  <math>\textstyle R = \sqrt{ h^2 + v^2 }</math>,  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.
:	<math>\textstyle d = \frac {half\ width\ of\ pano} {\operatorname{A2Npano}\ {(half\ hfov\ of\ pano)}}</math>,<br>
:	<math>\textstyle e  = \frac {half\ width\ of\ source} {\operatorname{A2Nsource}\ {(half\ hfov\ of\ source)}}</math>,
where A2Npano() and A2Nsource() are the ideal projection functions for panorama and lens.  There are many panorama projections but only two lens projections: 
:       <math>\textstyle {\operatorname{A2Nsource}( a )} = tan( a )</math>   for rectilinear lenses
:       <math>\textstyle {\operatorname{A2Nsource}( a )} =  a </math>   for fisheye lenses.
The scale factor from panorama to source coordinates is
:	<math>\textstyle s = \frac e d</math>.

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

The true source projection is unknown, so <math>e</math> is an estimate of <math>F_{source}</math> according to the fitted correction parameters.  Since hfov is one of those parameters, <math>e</math> will be proportional to the true <math>F_{source}</math>; the constant of proportionality will approach 1 as the fitted polynomial coefficients approach 0.  

In other words, <math>e</math> is a biased estimate of <math>F_{source}</math>.  However, the overall correction is equivalent to the generic one because the bias in the correction polynomial cancels the bias in the focal length.  The only real defect in the PanoTools scheme is that its parameters work for just one image format.

==Portable Correction Coefficients==

In the generic calibration scheme, dividing image coordinates by F makes it possible for the fitted correction parameters (apart from F) to be independent of both image format and physical pixel size, so that they apply to any image made with the given lens. As explained above, dividing image coordinates by any factor proportional to F is logically sufficient; however values other than F itself lead to non-portable parameter values.

In the PanoTools scheme, the "distance parameter" d, which is the focal length in panorama pixels, would be the appropriate divisor. That would make the argument of the radius scaling polynomial the ideal normalized radius,
:	<math>\textstyle N = \frac {R_{pano}} {F_{pano}}</math>
and the fitted coefficient values would be portable.

Alternatively the current non-portable coefficients can be converted using data available inside libpano.  With
:	<math>\textstyle k = \frac d {r_0}</math>,
:	<math>\textstyle w' =  w k </math>,
:	<math>\textstyle a' =  a k^2</math>,
:	<math>\textstyle b' =  b k^3</math>,
:	<math>\textstyle c' =  c k^4</math>
are the coefficients of a polynomial in <math>\textstyle N = \frac R d</math> 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
:       <math>\textstyle W  =  1</math>
:       <math>\textstyle A  =  a \frac k w</math>
:       <math>\textstyle B  =  b \frac {k^2} w</math>
:       <math>\textstyle C  =  c \frac {k^3} w</math>
So the portable radius mapping is
:       <math>\textstyle r = R ( 1 + A N + B N^2 + C N^3 )</math>
Along with the ideal function A2Nsource(), which gives N as a function of angle, this constitutes a portable lens correction function.

==Fully Portable Corrections==

To make a lens correction fully portable also requires expressing the fitted focal length in physical units rather than in pixels.  

The focal length in pixels must be known in order to compute, or to apply, any lens calibration, portable or not.  Physically, this quantity depends on lens focal length and camera properties.  Today, equipment manufacturers' specifications usually provide the needed data:
:	<math>\textstyle F_{pixels} = F_{mm} \frac {sensor\ width\ in\ pixels} {sensor\ width\ in\ mm}</math>.
Alternatively, the EXIF data from most modern high-end cameras includes the "focal plane resolution" field, which gives the physical pixel size directly.

In any practical calibration scheme F is actually an adjustable parameter.  However the fitted value is expected to be quite close to the one given by the physical specifications, which would of course be used as the initial value.  The main uncertainty is how accurately the nominal lens focal length reflects the true one, because for digital cameras the focal plane resolution is known exactly.  With <math> h </math> the width of a pixel in mm, the portable form of the fitted lens focal length is
:	<math>\textstyle F_{mm}  =  h F_{pixels} =  h e</math>, scale factor e defined above.
To adapt a portable correction to a given image it is only necessary to calculate <math>\textstyle F_{pixels} </math> from the calibrated <math>\textstyle F_{mm} </math> and the the pixel size associated with the image.
 
Unfortunately PanoTools does not now use any physical parameters, so fully portable corrections would have to be calculated, saved and restored by front-end software that has access to focal length and pixel size.  But if those were added to the PanoTools parameter set, libpano could handle fully portable corrections autonomously.


 -- 24 Jan 2010 T K Sharpless

[[Category:Community:Project]]