[Insight-users] Affine transform based registration in 3D

[marco giordano]

Hi,

I am running an Affine transform registration in 3D (basically
ImageRegistration9.cxx).

In order to support 3D mode, appart from some obvious changes, I did the
following:

Extend the optimizer scale:

....
 optimizerScales[0] =  1.0;
  optimizerScales[1] =  1.0;
  optimizerScales[2] =  1.0;
  optimizerScales[3] =  1.0;
  optimizerScales[4] =  1.0;
  optimizerScales[5] =  1.0;
  optimizerScales[6] =  1.0;
  optimizerScales[7] =  1.0;
  optimizerScales[8] =  1.0;
  optimizerScales[9] =  translationScale;
  optimizerScales[10] =  translationScale;
  optimizerScales[11] =  translationScale;
...


Get the full set of final parameters to print out

....
const double finalRotationCenterX = transform->GetCenter()[0];
  const double finalRotationCenterY = transform->GetCenter()[1];
  const double finalRotationCenterZ = transform->GetCenter()[2];
  const double finalTranslationX    = finalParameters[9];
  const double finalTranslationY    = finalParameters[10];
  const double finalTranslationZ    = finalParameters[11];
....

Then I guess I must retrieve the rotation angle by SVD
which I guess it is run on the 3x3 matrix but I do not know in which order I
must pick up the value. I am not sure that I got
this right but my guess would be:

(Guess)
 vnl_matrix<double> p(3, 3);
  p[0][0] = (double) finalParameters[0];
  p[0][1] = (double) finalParameters[1];
  p[0][2] = (double) finalParameters[2];

  p[1][0] = (double) finalParameters[3];
  p[1][1] = (double) finalParameters[4];
  p[1][2] = (double) finalParameters[5];

  p[2][0] = (double) finalParameters[6];
  p[2][1] = (double) finalParameters[7];
  p[2][2] = (double) finalParameters[8];


  vnl_svd<double> svd(p);
  vnl_matrix<double> r(3, 3);
  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 << " Scale 3         = " << svd.W(3)                 <<
std::endl;
  std::cout << " Angle (degrees) = " << angle * 45.0 / atan(1.0) <<
std::endl;


Does that make sense ?

I would appreciate some help on that.

Thank you

-- 
Marco G.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20091211/493bb896/attachment.htm>