[Insight-developers] ImageRegistration9

[Gunnar Farneback]

A few questions and notes about
Examples/Registration/ImageRegistration9.

1. The comments about testing it on two images say:

  //  The second image is the result of intentionally rotating the first
  //  image by $10$ degrees and then translating by $(-13,-17)$.  Both images
  //  have unit-spacing and are shown in Figure
  //  \ref{fig:FixedMovingImageRegistration9}. We execute the code using the
  //  following parameters: step length=1.0, translation scale= 0.0001 and
  //  maximum number of iterations = 300. With these images and parameters
  //  the registration takes $240$ iterations and produces
  //
  //  \begin{center}
  //  \begin{verbatim}
  //   239 44.3138   
  //  [0.984935, -0.172989, 0.172608, 0.984926, 
  //    123.249, 147.12, 9.58612, 17.9202]
  //  \end{verbatim}
  //  \end{center}

I can't quite reproduce this. With default parameter settings I end
with

237   50.7402   [0.986115, -0.170693, 0.168981, 0.985808, 111.202, 131.578, 12.4614, 16.0768]

which is close but the difference is larger than I would expect from
pure numerical variations. These numbers also seem to agree with a
random Dashboard I checked.

If I instead use the parameters indicated in the text I don't even get
close to the expected values.

Am I doing something wrong or has something changed since the
comments were written?

2. Later in the code there is a comment:

  //  The second component of the matrix values is usually associated with
  //  $\sin{\theta}$.

This is only valid if there's no significant scaling in the transform.
On the other hand it's not difficult to compute the rotation angle and
scaling parameters properly from the singular value decomposition of
the matrix, like this:

  vnl_matrix<double> p(2, 2);
  p[0][0] = (double) finalParameters[0];
  p[0][1] = (double) finalParameters[1];
  p[1][0] = (double) finalParameters[2];
  p[1][1] = (double) finalParameters[3];
  vnl_svd<double> svd(p);
  vnl_matrix<double> r(2, 2);
  r = svd.U() * vnl_transpose(svd.V());
  double angle = asin(r[1][0]);
  
  std::cout << " Scale 1         = " << svd.W(0)                 << std::endl;
  std::cout << " Scale 2         = " << svd.W(1)                 << std::endl;
  std::cout << " Angle (radians) = " << angle                    << std::endl;
  std::cout << " Angle (degrees) = " << angle * 45.0 / atan(1.0) << std::endl;

Unless someone thinks it's a bad idea to add this information to the
final output, I'll add a complete patch to the bug tracker.

3. Using the code above the affine parameters I get compute to a
   rotation of 9.77 degrees and scaling parameters 1.00135 and
   0.999613. The rotation angle doesn't look quite as impressive as
   the one in the code comments:

  // In this case the value $0.1729$ corresponds to a rotation 
  //  of $9.95$ degrees, which is approximately the intentional misalignment of 
  //  $10.0$ degrees.

/Gunnar