Precise meaning of control point coordinates?

[johnfi...@gmail.com]

Control points coordinates allow fractions of a pixel (1/3 pixel at a time when moving via arrow keys.  Even finer detail for "fine tuning" results).

What is the precise meaning of the coordinate relative to the pixel?

Imagine each pixel as a square in a very magnified image.  I would expect best behavior if that pixel's integer coordinates represented the center of that pixel.  But maybe some code would be simpler if the pixel's exact coordinates represented the top left corner of that pixel (its coordinates plus 0.5,0.5 would be its center).

I'm pretty sure various parts of of the hugin code are inconsistent with each other regarding that decision (following neither of the above rules and differently from other sections of the hugin code).

I'm working on some enhancements to hugin++, that I can only code correctly if I know the precise answer to this question.  So I looked for it in the current code and experimentally in the current behavior and found contradictions.  I'd like to know which parts are following the intended behavior, vs which are subtle bugs.

[Florian Königstein]

This is in fact an important question. I havn't looked at the code good enough to be able to give an answer.

I suggest looking at first at the code of the geometrical optimizer (fastPTOptimizer or libpano13). Probably most of the code with the math concerning the CP coordinates is the original code of Prof. Helmut Dersch. Somewhere the pixel coordinates together with the hfov will be converted into angles on the unit sphere or maybe normed coordinates on some "image plane". This could be a good place to look. Maybe I find the time and will look myself.

If some code in Hugin++ is not consistent with it, it should be adjusted.

After all, it wouldn't make sense to estimate the CP position with sub-pixel accuracy (in the fine tuning) if at the same time there are systematic errors of 0.5 pixels due to inconsistencies in the code.

[Bruno Postle]

On Sun, 15 May 2022, 11:34 <johnfi...@gmail.com> wrote:

Imagine each pixel as a square in a very magnified image.  I would expect best behavior if that pixel's integer coordinates represented the center of that pixel.  But maybe some code would be simpler if the pixel's exact coordinates represented the top left corner of that pixel (its coordinates plus 0.5,0.5 would be its center).

I remember much discussion about this on the old panotools list, the consensus was that Helmut Dersch got this correct in the libpano library (and Hugin uses the same conventions).

The way to test is to construct some artificial images with hard edge patterns and try optimising, stitching, reverse functions etc..

-- 

Bruno

[John Fine]

Maybe I'm misunderstanding something important, but:

On Sun, May 15, 2022 at 1:29 PM Florian Königstein <niets...@gmail.com> wrote:

I suggest looking at first at the code of the geometrical optimizer (fastPTOptimizer or libpano13). Probably most of the code with the math concerning the CP coordinates is the original code of Prof. Helmut Dersch. Somewhere the pixel coordinates together with the hfov will be converted into angles on the unit sphere or maybe normed coordinates on some "image plane". This could be a good place to look. Maybe I find the time and will look myself.

I cannot see how any of that (the optimizer) needs to be consistent regarding this question with any other part of hugin.  I think the optimizer would only need to be consistent with itself.  The consistency question matters elsewhere.

One significant need for consistency of this definition should be between the fine tuning code (which I haven't looked at myself yet) and the magnifier display.  On a low res monitor, the consistency between the magnifier and the 200% zoom mode should matter.  On my higher res displays, with the hugin++ options of 400% and 800% zoom, consistency matters across the 3-way: fine tuning vs. zoom vs magnifier.

If there is any significant difference in effective resolution between images (a Z translation value or an actual resolution difference), then the blending code has a very serious need for consistent definition vs fine tuning.  I'm hoping to make some changes in blending creating a need for that consistency even without any resolution difference.

[Gunter Königsmann]

If control points would link a pixel in one image to a pixel in another they could have integer coordinates.
But if they link a 4.5 pixel wide and 2.5 pixel high leave of a tree to the same leaf in a different photo the leaf might be 300.45 pixels more to the left and 2.4 pixels higher in one picture than in the other.

If you drop the digits after the comma the world might get a playstation-1-look.

[Florian Königstein]

... Sorry for repeated deletion and re-posting, but there seems to be no possibility the edit the message ...

One question is e.g.: What is the exact meaning of the hfov (horizontal field of view): Is it the angle measured from the left of the leftmost pixel area to the right of the rightmost pixel or from the mid of the left to the mid of the right pixel ?

The formulas in the optimizer are faily complicated. But assume that for a pixel the angle between the optical axis and the ray going through the "reference point" of the pixel area was calcualted, call it alpha. Assume the y coordinate is so that the y position is in the middle   [[  y = (H-1)/2 if the reference point of each pixel is in the middle of the pixel area ]].

W and H are the width and height of the image in pixels.

Assume in the optimizer the angle alpha is calculated as:

alpha = arctan( (x- (W-1)/2) / ((W-1)/2) * tan(hfov / 2) )

This would indicate that the reference point is in the middle of the pixels.

If the calculation would be:

alpha = arctan( (x - W/2) / (W/2) * tan(hfov / 2) ) 

This would indicate that the reference point is in the left (and upper) pixel area.

The actual formulas are more complicated, but it should be possible to get this information out of them.

[Florian Königstein]

Ok, another correction:

The second formula would indicate that the hfov is measured from the left of the left pixel area to the right of the right pixel area. The first formula would indicate that the hfov is calculated from the mid of the left pixel to the mid of the right pixel.

[Florian Königstein]

As already said, for the hfov the following possibilities exist:

1. The hfov is calculated from the mid of the left pixel to the mid of the right pixel

2. The hfov is measured from the left of the left pixel area to the right of the right pixel area

Further more, for the reference point of the pixels is possible:

A   The reference point is in the mid of the pixels

B   The reference point is at the top, left point of the pixels

Combining this, 4 possibilities would exist ( 1A,  1B,  2A,  2B ).

In my example formulas I assumed (but didn't know surely) that only the combinations 1A and 2B are possible / reasonable.

Of course, I didn't derive these formulas and it may be difficult to see what I wanted to say.

But you can simply insert the values x = 0 and x = W-1 (first formula) in order to get alpha = - hfov / 2 and alpha = + hfov / 2 .

Or x = 0 and x = W (second formula) in order to also get alpha = - hfov / 2 and alpha = + hfov / 2 .

The difference is that for alpha = + hfov you must insert x = W and not x = W-1 in the second formula.

This corresponds to the combination 2 B : The x coordinate of the right point in the rightmost pixel area is W is this case.

In the case 1A the x coordinate of the middle of the right pixel is W-1 and the x coordinate of the right of the right pixel area would be W-1 + 1/2  = W - 1/2

If I find the formulas of Helmut Dersch too complicated, I could try to find out by numerical experiments which case (1A, 2B) is used there.

[Florian Königstein]

But also the cases 1B and 2A could be identified in this way:

If the combination 1B was realized in the optimizer, the x-coordinate 0,5 must yield alpha = -hfov / 2, the coordinate x=W-0,5 must yield alpha =+hfov / 2 and x = W/2 must yield alpha = 0.

In this case the example formula would be:

alpha = arctan( (x- W/2) / ((W-1)/2) * tan(hfov / 2) )

If the combination 2A was realized in the optimizer, the x-coordinate -0,5 must yield alpha = -hfov / 2, the coordinate x=W-0,5 must yield alpha =+hfov / 2 and x = (W-1)/2 must yield alpha = 0.

Then the formula would be:

alpha = arctan( (x - (W-1)/2) / (W/2) * tan(hfov / 2) ) 

[Florian Königstein]

Without having looked at all the complicated code of Helmut Dersch, I found two lines of code from which I assume that the optimizer realizes  the combination 2A, i.e.

The hfov is measured from the left of the left pixel area to the right of the right pixel area .  The reference point is in the mid of the pixels.

These lines are in the file adjust.c:

First (in function distSphere() ):

                h2      = (double)optInfo->im[ n[j] ].height / 2.0 - 0.5;
                w2      = (double)optInfo->im[ n[j] ].width  / 2.0 - 0.5;
                execute_stack_new(      (double)optInfo->cpt[num].x[j] - w2,            // cartesian x-coordinate src
                                                (double)optInfo->cpt[num].y[j] - h2,            // cartesian y-coordinate src
                                                &x, &y, stack);

This corresponds to my terms  x - (W-1)/2 that occur in the cases 2A and 1A.

Second (in function SetInvMakeParams() ):

        case _rectilinear:
            mp->distance        = (double) pn->width / (2.0 * tan(b/2.0));
            break;

Probably this mp->distance is some kind of distance D that is used when modelling the camera with a simple pinhole model.

In the pinhole model the angle alpha would be calculated as:   alpha = arctan( x / D )  where D is the distance of the optical center (the pinhole) to a plane parallel to the image plane and going through the point in 3D space.

Here  D = mp->distance  .

If you insert it, you get:

     alpha = arctan (  ( x - (W-1)/2 ) / mp->distance ) = arctan (  ( x - (W-1)/2 ) / (W / (2*tan(b / 2))) )  .

Here b is hfov in radians, so I call it hfov_rad  ( in the C-Code hfov is used with the unit degrees ).

So:

    alpha = arctan (  ( x - (W-1)/2 ) / (W / (2*tan(hfov_rad / 2))) )  = arctan( (x - (W-1)/2) / (W/2) * tan(hfov_rad / 2) )

This is exactly the same as in my case 2A, so I assume that the case 2A is realized in the optimizer.

Maybe I will also test it with numerical values.

[Florian Königstein]

I also tested it numerically and could confirm that the optimizer realizes combination 2A.

I defined a fictive image with size 3000 x 2000 and placed one (part of a) CP onto it.

The y coordinate is in the y-middle, the hfov is 50 degrees and I tested it for the x-coordinates of the CP x = -0.5  ,   x = 2999.5  and   x = 2999.0 / 2.

According to case 2A they should yield the angles alpha = -25, alpha = +25 and alpha = 0.

The function execute_stack_new() yields two parameters x and y which are angles and not coordinates.

In this special case the x angle is equal to my angle alpha.

I used the following program to test it:

#include "filter.h"
#include <stdio.h>
#include <time.h>

void TestDistSphere(double x_cpt)
{
    double          x, y;  // Coordinates of control point in panorama
    double          w2, h2;
    int j;
    Image sph;
    struct  MakeParams      mp;
    struct  fDesc           stack[15];
    CoordInfo b[2];
    CoordInfo cp;
    double lat[2], lon[2];  // latitude & longitude
    double dlon;
    double dangle;
    double dist;
    double radiansToPixelsFactor;

    int num = 0;

    SetImageDefaults(&sph);

    sph.width = 360;
    sph.height = 180;
    sph.format = _equirectangular;
    sph.hfov = 360.0;

    Image img;
    optVars opt;
    controlPoint cpt;
    AlignInfo ai;
    ai.im = &img;
    ai.opt = &opt;
    ai.cpt = &cpt;
    ai.numIm = 1;
    ai.numPts = 1;
    ai.numParam = 3;
    ai.fcn = fcnPano;
    AlignInfo* optInfo = &ai;

    double distanceComponents[2];

    SetImageDefaults(&img);

    img.width = 3000;  // just a fictive value
    img.height = 2000;  // just a fictive value
    img.hfov = 50;     // just a fictive value
    img.yaw = 0;
    img.pitch = 0;
    img.roll = 0;
    cpt.num[0] = cpt.num[1] = 0;
    cpt.type = 0;
    cpt.x[0] = cpt.x[1] = x_cpt;
    cpt.y[0] = cpt.y[1] = 0.5 * (img.height - 1.0);

    j = 0;
    {
        SetInvMakeParams(stack, &mp, &optInfo->im[0], &sph, 0);

        h2 = (double)optInfo->im[0].height / 2.0 - 0.5;
        w2 = (double)optInfo->im[0].width / 2.0 - 0.5;


        execute_stack_new((double)optInfo->cpt[num].x[j] - w2,            // cartesian x-coordinate src

            (double)optInfo->cpt[num].y[j] - h2,            // cartesian y-coordinate src
            &x, &y, stack);

        printf("x = %.16g    y = %.16g\n", x, y);
    }
}

void main()
{
    TestDistSphere(-0.5);
    TestDistSphere(-0.5 + 3000.0);
    TestDistSphere(2999.0 / 2.0);
}

Output:

x = -25    y = 0
x = 25    y = 0
x = 0    y = 0