Computing a homography between 2D scene and image

      #include <gandalf/vision/homog33_fit.h>

If a part of the viewed scene is planar, or the camera is undergoing a pure rotation (or both), the (part of the) scene can be reconstructed using 2D methods. Here we assume a point-cloud representation, so the scene is represented by $n$ points ${\bf X}_i$ in homogeneous coordinates, $i=1,\ldots,n$. The relationship between the ${\bf X}_i$ and points ${\bf x}_i$ in an image of the same (part of the) scene is a simple linear projective transformation or homography:

$${\bf x}_i = \lambda_i P {\bf X}_i$$ (5.3)

$P$ is a $3\times 3$ matrix representing the homography and $\lambda$ is a scale factor. This equation also assumes that the camera employed in projecting the points onto the image is linear, but if the camera is non-linear AND the camera parameters are known, the distortion can be removed first by applying the function gan_camera_remove_distortion_[qi]() to the image points ${\bf x}_i$ as described in Section 5.1. Given four or more point correspondences in the image (in general position), the homography matrix $P$ can be computed. This can be done by first eliminating $\lambda$ to obtain homogeneous linear equations for the elements of $P$. Given that ${\bf X}=(X\;Y\;Z)^{\top}$ and ${\bf x}=(x\;y\;z_h)^{\top}$, we can obtain the equations

$$x {\bf P}_3.{\bf X}- z_h {\bf P}_1.{\bf X}= 0,\;\;\; y {\bf P}_3.{\bf X}- z_h {\bf P}_2.{\bf X}= 0$$ (5.4)

where $P$ is separated into rows as

$$P = \left(\!\!\begin{array}{c} {\bf P}_1^{\top}\ {\bf P}_2^{\top}\ {\bf P}_3^{\top}\end{array}\!\!\right)$$

From four points we get eight such equations, which allows $P$ to be computed up to a scale factor using the same symmetric eigensystem routines as are used to solve for the fundamental and essential matrices above.

Note that this formulation differs from the normal formulation which considers the homographies between images. That is a special case of our formulation, because we can take an image as the projective ``scene'' representation ${\bf X}_i$. The scene/image formulation also allows us to represent the motion over a sequence of $k$ images in a compact way as the set of homographies $P{\scriptstyle (j)}$ for images $j=1,\ldots,k$ mapping the scene ${\bf X}_i$ to each set of image points ${\bf x}_i{\scriptstyle (j)}$, rather than as an arbitrary collection of pairwise homographies.

To start the calculation, define an accumulated symmetric matrix eigensystem structure and initialise it using the following routine:

      Gan_SymMatEigenStruct SymEigen;

      /* initialise eigensystem matrix */
      gan_homog33_init ( &SymEigen );

Then for each point correspondence, build the equations 5.4 and increment the accumulated symmetric eigensystem matrix by calling the following function:

      int iEqCount=0, iCount;
      Gan_Vector3 v3X, v3x; /* declare scene and image points X & x */

      for ( iCount = 0; iCount < 100; iCount++ )
      {
         /* ... build scene and image point coordinates into X and x ... */

         /* increment matrix using point correspondence */
         gan_homog33_increment_p ( &SymEigen, &v3X, &v3x, 1.0, &iEqCount );
      }

The fourth argument 1.0 is a weighting factor for the equations as described in Section 3.2.2.15. The last argument iEqCount is a running count of the total number of equations processed thus far, to be passed below to the function to solve for $P$.

Once the point correspondences have been processed in this way, you can solve the equations using

      Gan_Matrix33 m33P; /* homography matrix P */

      gan_homog33_solve ( &SymEigen, iEqCount, &m33P );

to compute the homography $P$. If you want to repeat the calculation of a homography with new data, you can start again by calling

      gan_homog33_reset ( &SymEigen );

At the end of the homography calculation(s) you can free the eigensystem structure using the function

      gan_homog33_free ( &SymEigen );

Given correspondences between lines, it is also possible to generate homogeneous linear equations for $P$ and either combine with points or compute $P$ purely from lines. To see how to derive the equations for lines, take the line equations

$${\bf L}.{\bf X}= 0,\;\;\;\;{\bf l}.{\bf x}= 0$$

define the homogeneous line parameters ${\bf L}$ in the scene and ${\bf l}$ in the image. We can derive the relationship between ${\bf L}$, ${\bf l}$ and $P$ using the point projection equation 5.3, yielding

$${\bf L}= \mu P^{\top}{\bf l}$$

for a scale factor $\mu$. Separating $P$ into columns as

$$P = \left(\!\!\begin{array}{ccc} {\bf P}'_1 & {\bf P}'_2 & {\bf P}'_3 \end{array}\!\!\right),$$

eliminating $\mu$ from the above equation, and writing ${\bf L}= (L_X\;L_Y\;L_Z)^{\top}$, we obtain the two homogeneous linear equations

$$L_X {\bf P}'_3.{\bf l}= L_Z {\bf P}'_1.{\bf l},\;\;\; L_Y {\bf P}'_3.{\bf l}= L_Z {\bf P}'_2.{\bf l}.$$

Given correspondence between a known scene line ${\bf L}$ and a known image line ${\bf l}$, the following routine generates these equations and accumulates them in the calculation of $P$:

         Gan_Vector3 v3L, v3l; /* declare scene line L and image line l */

         /* ... fill L and l with values for corresponding lines ... */

         /* increment matrix using line correspondence */
         gan_homog33_increment_l ( &SymEigen, &v3L, &v3l, 1.0, &iEqCount );

This is assuming that the endpoints of the scene line are unknown. In practice the scene line will normally be created from previous matching of image lines, which are line segments, so that the endpoints ${\bf X}_1$ and ${\bf X}_2$ of the line in scene coordinates will be approximately known. Note that we don't depend on locating the actual endpoints of the line accurately, which is a notoriously difficult problem. You should think of the two points ${\bf X}_1$ and ${\bf X}_2$ instead as representative points on the line. In this case there is an alternative way of incorporating the line information which seems to give better numerical performance. We note that the scene line endpoints ${\bf X}_1$ and ${\bf X}_2$ should project onto the image line ${\bf l}$, so we obtain

$${\bf l}.(P{\bf X}_1) = 0,\;\;\;\;{\bf l}.(P{\bf X}_2) = 0$$

These are homogeneous linear equations in the elements of $P$ which can be directly fed into the accumulated matrix calculation for $P$, using the routine

      Gan_Vector3 v3X1, v3X2; /* declare scene line endpoints X1 & X2 */
      Gan_Vector3 v3l; /* image line homogeneous coordinates l */

      /* ... set X1, X2 and l for corresponding scene line and image line ... */

      /* add equations for two endpoints */
      gan_homog33_increment_le ( &SymEigen, &v3X1, &v3l, 1.0, &iEqCount );
      gan_homog33_increment_le ( &SymEigen, &v3X2, &v3l, 1.0, &iEqCount );

Error detection: gan_homog33_init() returns a pointer to the initialised structure, and returns NULL on error. All the other routines except the void routine gan_homog33_free() return a boolean value, which is GAN_FALSE on error. The Gandalf error handler is invoked when an error occurs.