What is optimized FOV really?

[Aron H]

I'm a new member of the group, and I'm looking again at Hugin and
PanoTools after a long time.

I work on a product that combines panoramic photography with a range
scanner, our website is here if you're curious:
http://www.deltasphere.com/

We take color photos using a Nikon D300 and a Nikkor 14mm rectilinear
lens, covering 360 x 150 or so with 18 photos, portrait orientation, 9
aimed high and 9 low.

We calibrate the camera and lens combination in our lab, and use a
fixed focus and F-stop.
We use the Bouguet calib toolbox to do this:
http://www.vision.caltech.edu/bouguetj/calib_doc/index.html

What I'm trying to figure out is the correspondence between the
Bouguet camera calibration and the Hugin model. I can't figure out how
the FOV that is optimized and displayed in Hugin is used to get real-
world 3D rays.

Specifically, the a, b, c distortion numbers use the normalized
radius, which is set to be 1.0 at h/2 (if h is the smaller image
dimension). This doesn't refer to real-world size or angles at all.

Hugin figures out where each photo is pointing in the world. If I want
to find out where one of the pixels in a photo is pointing, I need to
use the FOV and the distortion numbers to make a ray that points to
that pixel. Do I use the FOV before the distortion, or after?

I'm having a hard time stating this clearly, but what I'd really like
is the math used to take a 3D point as seen from the nodal point of
the panorama, and figure out which pixel it hits in a photo.

We already do this with the Bouguet camera calibration, and I have
lots of notes and a spreadsheet comparing Hugin and Bouget, but I
can't get it to quite correspond...

P.S. Thanks to Carl von Einem for pointing me to the developers list,
I didn't know it existed from looking at the wiki, and I knew this
question wasn't good for the user's list...

Thanks,
Aron

[Klaus Foehl]

Hello Aron,

If I understand hugin and the undelying panotools correctly, then the
maths direction is already in the direction you require. In the output
image, equirectangular, you already have theta and phi angle in 3D for
each pixel. Distance to object is missing of course. With this (x,y)
which in equirectangular projection is (phi,theta) hugin, well
panotools needs to know which camera pixel it needs to look at.

So one way would be to look up the formulas in the source code.

I presume that output coordinates are first transformed into focal
plane coordinates using standard lense formula and for this the hfov
number is used. Then a new radius is computed from the a,b,c plus d
formula. Of course after that correction for a non-square input image
the hfov angle is not the real horizontal angle as captured with the
sensor. I would cling to points at h/2 radius as these are not
affected by the a+b+c+d correction. Given that your input images have
a know aspect ratio, for the maths to start I would convert the hfov
value into a vfov value which survives the correction unchanged.

Hope that helps

Klaus

[Aron H]

Hi Klaus,

On Jul 1, 1:31 pm, Klaus Foehl <k...@ph.ed.ac.uk> wrote:
> Hello Aron,
>
> If I understand hugin and the undelying panotools correctly, then the
> maths direction is already in the direction you require. In the output
> image, equirectangular, you already have theta and phi angle in 3D for
> each pixel. Distance to object is missing of course. With this (x,y)
> which in equirectangular projection is (phi,theta) hugin, well
> panotools needs to know which camera pixel it needs to look at.
>
> So one way would be to look up the formulas in the source code.

Can someone point me in the right direction? The doxygen pdf has 300+
pages :)

>
> I presume that output coordinates are first transformed into focal
> plane coordinates using standard lense formula and for this the hfov
> number is used. Then a new radius is computed from the a,b,c plus d
> formula. Of course after that correction for a non-square input image
> the hfov angle is not the real horizontal angle as captured with the
> sensor. I would cling to points at h/2 radius as these are not
> affected by the a+b+c+d correction. Given that your input images have
> a know aspect ratio, for the maths to start I would convert the hfov
> value into a vfov value which survives the correction unchanged.

Let me attempt some pseudo-code, which a developer can correct?

Given a pixel in the final panorama image, find what sub-pixel it maps
to
in a given input image:
-Given 0,0 in the image center, and FOV, find a yaw/pitch for the
pixel, and make a direction vector.
-Rotate that vector by the image's yaw-pitch-roll angles, so it is now
in camera coordinates
-Convert the hfov to a vfov using the aspect ration, as described
here: http://wiki.panotools.org/Field_of_View
-Use the condition that r=1 at the h/2 radius, find the depth of the
image plane, zf
tan(vfov/2) = r/zf
zf = r / tan(vfov/2)
- If our direction vector is (vx, vy, vz), mulitply it by zf/vz, to
get (ix, iy), a point on the image plane. Points behind the camera or
too far away can be thrown out.
- Calculate r = sqrt(ix^2 + iy^2), and distort by multplying (ix, iy)
* (a r^3 + b r^2 + c r + d) = (dx, dy)
- multiply (dx, dy) * h/2 to get a number in pixels
- Interpolate a color from the surrounding pixels in the source image.

Is that actually correct? Can someone point me to where this is done
in the source code?
BTW, I do OpenGL as part of my job, and I'd be happy to answer
questions/help on that subject....
>
> Hope that helps
>
> Klaus
Regards,
Aron

[Bruno Postle]

On Wed 01-Jul-2009 at 11:56 -0700, Aron H wrote:
>
>Given a pixel in the final panorama image, find what sub-pixel it maps
>to
>in a given input image:
>-Given 0,0 in the image center, and FOV, find a yaw/pitch for the
>pixel, and make a direction vector.
>-Rotate that vector by the image's yaw-pitch-roll angles, so it is now
>in camera coordinates

>-Convert the hfov to a vfov using the aspect ration, as described
>here: http://wiki.panotools.org/Field_of_View

>-Use the condition that r=1 at the h/2 radius, find the depth of the
>image plane, zf
> tan(vfov/2) = r/zf
> zf = r / tan(vfov/2)
>- If our direction vector is (vx, vy, vz), mulitply it by zf/vz, to
>get (ix, iy), a point on the image plane. Points behind the camera or
>too far away can be thrown out.
>- Calculate r = sqrt(ix^2 + iy^2), and distort by multplying (ix, iy)
>* (a r^3 + b r^2 + c r + d) = (dx, dy)
>- multiply (dx, dy) * h/2 to get a number in pixels
>- Interpolate a color from the surrounding pixels in the source image.

Correct, though your assumption is for landscape format images, it
is slightly simpler for portrait.

>Is that actually correct? Can someone point me to where this is done
>in the source code?

math.c in libpano13.

--
Bruno

[Aron H]

On Jul 1, 5:57 pm, Bruno Postle <br...@postle.net> wrote:

> Correct, though your assumption is for landscape format images, it
> is slightly simpler for portrait.

I think you mean that if the input images actually have a height
greater than their width, Hugin will still display the horizontal FOV,
which can then be used directly, because it's the FOV for the smaller
image dimension.

For our typical panorama, we take pictures with the camera in the
portrait orientation, but of course they are saved by the camera in
landscape. After optimizing, Hugin finds a roll angle near 90 degrees.

>
> >Is that actually correct? Can someone point me to where this is done
> >in the source code?
>
> math.c in libpano13.

Thanks! I'll check it out.
Is anyone interested in a writeup comparing Bouguet camera calibration
to Hugin calibration? Somewhere in the wiki?

>
> --
> Bruno

[Daniel Reetz]

> Thanks!  I'll check it out.
> Is anyone interested in a writeup comparing Bouguet camera calibration
> to Hugin calibration? Somewhere in the wiki?

I am very interested. I have a big research project (currently on
hold) that would benefit from such a discussion.

[Klaus Foehl]

On 2 July, 16:31, Aron H <aron.hel...@gmail.com> wrote:

> Is anyone interested in a writeup comparing Bouguet camera calibration
> to Hugin calibration? Somewhere in the wiki?

Yes indeed. I'd like to see how much is maths/physics/symmetry based
(like parameter b in hugin/panotools) and where phenomenology-only
(like a and c) is used in the Bouguet camera calibration (and yes, I
shall read the introductory page at Caltech...)

Cheers

Klaus

[Klaus Foehl]

Hello,

Quick addition and summary from my side.

1) Barrel distortion correction: panotools uses r* = a*r^4 + b*r^3 +
c*r^2 + d*r,
with r being the normalised radius. If it were not for the condition a
+b+c+d=1
this parameter d would be fully collinear with the focal length.
focal_length_panotools * d = true_focal_length

2) http://www.vision.caltech.edu/bouguetj/calib_doc/htmls/parameters.html
x_d=x_n * (1+kc(1)*r^2+kc(2)*r^4+kc(5)*r^6) +dx, dx=[...]
Looking at the distortion formulas the Bouguet model looks
mathematically
sound to me, as the parameters x_d and dx only use odd and even
powers,
respectively, of (r,x,y) polynomial expressions.
That respects the lense symmetry. Parameters a and c in panotools do
not.

3) In the xd formula 1 is equivalent to d, kc(1)=b/d, kc(2) and kc(5)
are unaccounted for. (N.B. I suggested a few times to include
a kc(2) type parameter into hugin/panotools).

4) Sheer parameters are available in panotools but are not shown
in the hugin GUI. There are two parameters, and combined with
the roll parameter I suspect that this set of three only spans
a two-dimensional configuration space, save for possibly a
higher order correction that results in non-square-like pixels.

5) Tangential distortion vectors, kc(3) and kc(4) in Bouguet,
I have not seen such in hugin/panotools.

6) Principal point cc in Bouguet are parameters e and f
in panotools, in the GUI labelled x and y.

7) Rectangular non-square pixels, Bouguet using fc(1) and fc(2),
are not there in the hugin/panotools model.

Summary: models are reasonably similar, Bouguet is more powerful,
and parameters a and c of hugin/panotools, the only ones
I cannot find a match in Bougues, have a maths symmetry issue.

Cheers
Klaus

[Aron H]

On Jul 3, 4:08 am, Klaus Foehl <k...@ph.ed.ac.uk> wrote:
> Hello,
>
> Quick addition and summary from my side.
>
> 1) Barrel distortion correction: panotools uses r* = a*r^4 + b*r^3 +
> c*r^2 + d*r,
> with r being the normalised radius. If it were not for the condition a
> +b+c+d=1
> this parameter d would be fully collinear with the focal length.
> focal_length_panotools * d = true_focal_length
>

> 2)http://www.vision.caltech.edu/bouguetj/calib_doc/htmls/parameters.html

> x_d=x_n * (1+kc(1)*r^2+kc(2)*r^4+kc(5)*r^6) +dx, dx=[...]
> Looking at the distortion formulas the Bouguet model looks
> mathematically
> sound to me, as the parameters x_d and dx only use odd and even
> powers,
> respectively, of (r,x,y) polynomial expressions.
> That respects the lense symmetry. Parameters a and c in panotools do
> not.

I've looked through Bouguet's site for an explanation of the symmetry
issue, and I think he pushes it off to publications. I guess I don't
know what the argument is, except perhaps that r^2 and r^4 distortion
are expected from the lens design, and the r and r^3 are not. The
formulas only deal with positive r values, and apply distortion
totally in the radial direction, so they are 'symmetric' in that
sense. Is there a good explanation of this somewhere?

>
> 3) In the xd formula 1 is equivalent to d, kc(1)=b/d, kc(2) and kc(5)
> are unaccounted for. (N.B. I suggested a few times to include
> a kc(2) type parameter into hugin/panotools).

Note, however, that the 'r' parameter is scaled differently in Hugin
and Bouguet. I'll write up the derivation....

>
> 4) Sheer parameters are available in panotools but are not shown
> in the hugin GUI. There are two parameters, and combined with
> the roll parameter I suspect that this set of three only spans
> a two-dimensional configuration space, save for possibly a
> higher order correction that results in non-square-like pixels.
>
> 5) Tangential distortion vectors, kc(3) and kc(4) in Bouguet,
> I have not seen such in hugin/panotools.
>
> 6) Principal point cc in Bouguet are parameters e and f
> in panotools, in the GUI labelled x and y.
>
> 7) Rectangular non-square pixels, Bouguet using fc(1) and fc(2),
> are not there in the hugin/panotools model.
>
> Summary: models are reasonably similar, Bouguet is more powerful,
> and parameters a and c of hugin/panotools, the only ones
> I cannot find a match in Bougues, have a maths symmetry issue.
>
> Cheers
> Klaus

I agree with the rest, thanks for the nice summary,
Aron

[Tom Sharpless]

Thanks all for a useful discussion.

It is standard practice in lens calibration to compute a "corrected"
radius from the observed radius of an image point by means of a
polynomial. However I think it is becoming clear that PanoTools and
its clients such as Hugin are almost alone in including odd-order
terms in the correction polynomial. The even order terms correspond
to the assumed rotational symmetry of the lens, so are well behaved in
numerical optimization. Odd terms beyond the first do not correspond
to any realistic physical property, and tend to destabilize numerical
optimization.

So PT would be better off also using an even-orders-only polynomial
for radius correction. That may not just be a matter of defining b
to be always 0, however, because, as Klaus points out, in the present
scheme the focal length (equivalently fov) is mixed into the
correction polynomial rather than being a separate parameter (as it
should be).

As I see it, the right way to normalize the correction is to divide
the measured radius by the focal length before applying the
polynomial: let
f = current estimate of focal length;
u = r / f;
R = u + p2 * u ^2 + p4 * u^4;
then R is the proper argument to the trigonometric function that
defines the ideal view angle (in radians), either
Angle = atan( R ) for a rectilinear lens, or
Angle = k * asin( R / k ) for an equal-area type fisheye (design
parameter k), or even
Angle = R for a true equal-angle fisheye.

I believe that using this scheme would improve the stability (and so
the average quality) of the numerical results delivered by the PT
optimizer.

Regards, Tom

[Turk]

Look in Modern Optical Engineering, by Warren J. Smith. The aberration
polynomial is given on p. 58.

The terms called "distortion" are B2 (cubic term) and C12 (quintic
term), and correspond to the Bouguet's quadratic and quartic
correction.They are a function only of a power of h (distance from the
optical axis).

The quadratic term (linear PanoTools) is also a function of angle
(i.e. h^2 * cosθ^n), and has three parts: astigmatism (B3), Petzval
(B4), and oblique spherical distortion (C4, C5, C6).

The quartic term (cubic PanoTools) is also a function of angle:
Petzval (C10) and astigmatism (C11).

In other words, the terms corresponding to odd terms in PanoTools are
not simple functions of merely the distance from the optical axis, but
also are a function of the angle around the optical axis.

BTW, Jean-Yves Bouguet had updated the calibrator for fisheye lenses,
and promised to push the changes out into the open source. Did this
happen, or should I bug him?

-Ken

[Bruno Postle]

On Sun 05-Jul-2009 at 11:17 -0700, Tom Sharpless wrote:
>
>So PT would be better off also using an even-orders-only polynomial
>for radius correction. That may not just be a matter of defining b
>to be always 0, however, because, as Klaus points out, in the present
>scheme the focal length (equivalently fov) is mixed into the
>correction polynomial rather than being a separate parameter (as it
>should be).

The panotools polynomial is normalised to the angle of view, i.e. it
is normalised to the vertical or the horizontal angle of view,
whichever is smaller.

This _is_ confusing as it means the 'v' horizontal angle of view
parameter means something completely different for landscape and
portrait images.

e.g. for my Peleng fisheye images, the horizontal angle of view in
portrait orientation is 111 degrees and so is the 'v' parameter.

For landscape images the 'v' parameter is 168, but the lens
correction parameters mean that the actual angle of view covered is
173 degrees.

I think the portrait situation makes a lot more sense, it is sane to
normalise the lens distortion to the narrowest width of the image
since this is well behaved.

--
Bruno

[Tom Sharpless]

Hi all

Bruno is right that the way PT normalizes the lens correction is
confusing.

It is also inappropriate, because it makes the correction coefficients
depend on the focal length. Those should be independent, not just to
reduce confusion, but to help stabilize optimization. As Ken has
pointed out, the proper definition of focal length is 1 /
(magnification at center), or in calculus terms, f = limit as r->0 of
dr/da, a being view angle. This is independent of the field of view,
which of course also depends on the projection function. So you have
to know the projection function in order to calculate f from fov,
which is a bit odd when what you are trying to do is determine an
unknown projection function.

You need f to calculate view angles from measured radii, no matter
what ideal lens function and correction scheme you choose to use,
because it is the basic scale factor relating angles to distances in
the image plane. In

Dersch's decision to ask for the fov instead of the f.l. is perfectly
reasonable, given that it can be hard to discover the pixel dimensions
for many digital cameras (we almost always know the focal length in
mm, but we need f in pixels). But that information is much more
likely to be available nowadays, and should be used when it is. And
we should be very clear that an f calculated from a stated fov and
image dimension is just a first guess. The true fov can only be
computed from the true f and the true projection function, including
correction terms.

Even with f properly removed from the correction function, it is still
necessary to normalize that function so that it does not arbitrarily
change the apparent fov. An even order polynomial can't change the
central magnification, of course, but it can change it everywhere
else. This normalization amounts to choosing a radius that does not
change when the correction coefficients are adjusted. I see no strong
reason not to keep doing this Derch's way, that is, fixing the radius
corresponding to the smaller dimension of the image (though on second
thought there might be problems with very long narrow images -- the
[harmonic?] average of the two dimensions might be a better choice).
But it really doesn't matter much what fixed radius one chooses, the
important thing is that choosing one makes the correction coefficients
independent of f. Any particular choice just adds a small constant
"error" to the fitted f -- which doesn't matter a bit, as ultimately
we don't care what f is, but what angles our image really covers.

Regards, Tom

[Klaus Foehl]

Hello all,

One further thought about the cubic and quintic terms only in
distortion correction.

Let us assume that the mapping from the solid angle space onto the
focal plane sensor is mathematically well behaved. One standard way is
to parametrise this into a polynomial, a standard two-dimensional
polynomial which uses cartesian coordinates to avoid the coordinate
singularity in polar coordinates. In a nutshell we are talking a
smooth function (x',y')=f(x,y). In general the coordinate systems need
not even be centered along the optical axis.

If we now define a special point in this remapping and further ask
that our remapping conserves the polar angle with reference to this
point, it means the function f has to be composed of f_i=(x,y)*r^n
elements. With the ansatz that we want to write these as polynomials
as well, holomorphic functions, n=2m has to be even and r^n=(x^2+y^2)
^m.

r=sqrt(x^2+y^2) cannot be written as a polynomial, and is not
holomorphic for r=0.

Wording this aspect as a symmetry invariant condition, the mapping
formula when given in polar coordinates should map (+radius,angle) and
(-radius,angle+180deg) onto the same point.

Cheers
Klaus

[Aron H]

On Jul 10, 4:38 am, Klaus Foehl <k...@ph.ed.ac.uk> wrote:
> Hello all,
>
> One further thought about the cubic and quintic terms only in
> distortion correction.
>
> Let us assume that the mapping from the solid angle space onto the
> focal plane sensor is mathematically well behaved. One standard way is
> to parametrise this into a polynomial, a standard two-dimensional
> polynomial which uses cartesian coordinates to avoid the coordinate
> singularity in polar coordinates. In a nutshell we are talking a
> smooth function (x',y')=f(x,y). In general the coordinate systems need
> not even be centered along the optical axis.
>
> If we now define a special point in this remapping and further ask
> that our remapping conserves the polar angle with reference to this
> point, it means the function f has to be composed of f_i=(x,y)*r^n
> elements. With the ansatz that we want to write these as polynomials
> as well, holomorphic functions, n=2m has to be even and r^n=(x^2+y^2)
> ^m.
>
> r=sqrt(x^2+y^2) cannot be written as a polynomial, and is not
> holomorphic for r=0.

I think I get it. You are requiring that we can write the remapping
function as a polynomial in x and y, not just in r. That rules out
using odd powers of r, because " r=sqrt(x^2+y^2) cannot be written as
a polynomial". Even powers of r can be written as a pure polynomial in
x and y.

> Wording this aspect as a symmetry invariant condition, the mapping
> formula when given in polar coordinates should map (+radius,angle) and
> (-radius,angle+180deg) onto the same point.
>
> Cheers
> Klaus

I appreciate the theory presented here, but I think we have drifted a
bit from some practical considerations. I suspect it would be a large
amount of work to change or extract the distortion equations from the
Hugin / PanoTools toolchain, across the multiple tools that use them.
So any developer or Google SOC student :) which might undertake this
would need to be convinced of the benefits. In theory, the benefit is
numerical stability, and matching best-practice academic research in
this area. However, are there practical examples of PanoTools failing
or acting badly because of its distortion model? If this change can
solve real-world problems, it might be worth it....
Regards,
Aron

[Aron H]

I'm happy to report that I am able to duplicate the Hugin transform
and camera distortion model in my software, using the basic outline I
described earlier.

I also wanted to report that I compared the results of camera
calibration in Hugin and Bouguet, and they agree very well, in a
limited case.
Allow Hugin to optimize FOV, only the 'b' distortion parameter, and
the 'd' and 'e' center shift.
Allow Bouguet to optimize FOV and center shift, and only the k0
distortion parameter.
k0 and b are the coefficient on the r^2 distortion in both camera
models.

I found that the center shift and r^2 distortion (after scaling)
agreed almost exactly. Here are the numbers:
Center shift:
Bouguet 17.32 27.47
Hugin 17.3 27.5
r^2 distortion
Bouguet -0.01138
Hugin -0.01138
Amazing, considering the input to Bouguet is 21 pictures of a
calibration grid, and the input to Hugin is an 18 picture 360 deg
panorama! I haven't figured out the formula to compare FOV/focal
length directly yet, although
if I follow Klaus' suggestion that

> focal_length_panotools * d = true_focal_length

then I get this for focal length, in pixels:
Bouguet 1950.0
Hugin 1954.91

I will get that promised comparison up on the wiki shortly, explaining
the scaling I did between the two systems.

Regards,
Aron

[Klaus Foehl]

Hello,

On 13 July, 20:31, Aron H <aron.hel...@gmail.com> wrote:
> [...]

> In theory, the benefit is
> numerical stability, and matching best-practice academic research in
> this area. However, are there practical examples of PanoTools failing
> or acting badly because of its distortion model? If this change can
> solve real-world problems, it might be worth it....

Badly failing it may not be, but on the pixel level it does.
I have always wondered, why with my Canon A610
my alignment is good at typically 0.3-1.0 pixel width
with y,p,r,b,x,y optimisation. Hence yesterday I took
a test image pair, 50% overlap, and first placed
control points in the upper half of the overlap area.

Then in a second step I placed control points
all over the overlap area. Both times I optimised
y,p,r,v,a,b,c,x,y. And I took care that the
control points were good ones.

Mean deviation and average deviation increases
by 10% and 20%. Parameters a and c change
by ballpark 10%. Most dramatically adding
the extra control points decreases the y value by half.

It means that determining the lense parameters these values a,b,c are
not independent on where one places the control points.

This means that with my camera lense there is some distortion which is
not accounted for in the optical model, and the currently available
parameters are not fully orthogonal to that missing degree of freedom.

While not explicitly investigated this time, I'd expect largest
alignment problems outside the area covered with control points.
Extrapolation is known in maths to have such problems.

Trying to understand this, from the outset three parameters is not
that small a number. But non-zero parameters a and c introduce a
mathematical singularity at r=0. A usual lense may have rotational
symmetry at r=0, but a singularity there is not. Hence a and c cover
something else than the expected lense behaviour.

So I would suggest to introduce at least a quintic term and name it aa
so the sequence looks like
aa,a,b,c and for the time keep a and c as quite a few tools still rely
on these.

If anything, I have learned that my imperfect alignment is most
probably not due to a shortcoming in quality of the "Finetune"
algorithm results.

Cheers

Klaus