[Insight-users] Problems with VersorRigid3Transform

[Anders Almås]

Hi!

I am performing registration of 3D images and a I am using
VersorRigid3DTransform. My goal is to create a test set by using
VersorRigid3DTransform and then try to registrate it back to the
starting point.

My problem is to do the set up for this test... It might be a dumb
question but... 

My code is the following:

 
  typedef itk::ResampleImageFilter<ImageType,ImageType> FilterType;
//The FilterType is for checking the transform/rotation..

  FilterType::Pointer filter = FilterType::New();
  
  const int offsetx = par[0]; //Translation parameters in
  const int offsety = par[1]; // x-, y- and z direction
  const int offsetz = par[2];
  

  typedef itk::VersorRigid3DTransform<double> TransformType;	

  TransformType::Pointer  rotation = TransformType::New();
  TransformType::VersorType VersorType;

 typedef VersorType::VectorType VectorType; 
 
  TransformType::Pointer transform = TransformType::New();
  TransformType::OutputVectorType translation;
  translation[0] = offsetx;
  translation[1] = offsety;
  translation[2] = offsetz;
  
 transform->SetTranslation(translation);
  
  typedef itk::LinearInterpolateImageFunction<ImageType,double>
InterpolateType;

  InterpolateType::Pointer interpolator = InterpolateType::New();

  filter->SetInterpolator(interpolator);

  filter->SetDefaultPixelValue(0);

  double spacing[ImageDimension]; // 3D...

  spacing[0] = 1.453; // Spacings between the pixels..
  spacing[1] = 1.453;
  spacing[2] = 2.0;
  
  filter->SetOutputSpacing(spacing);

  double origin[ImageDimension];

  origin[0]=0;
  origin[1]=0;
  origin[2]=0;

  filter->SetOutputOrigin(origin);

  ImageType::SizeType size;

  size[0] =  input->GetBufferedRegion().GetSize()[0];
  size[1] =  input->GetBufferedRegion().GetSize()[1];
  size[2] =  input->GetBufferedRegion().GetSize()[2];

  filter->SetSize(size);

  
  filter->SetInput(input);
  TransformType::OutputVectorType translation1;

//Calculate a imagecenter.. 

  const double imageCenterX = origin[0] + spacing[0] * size[0] / 2.0;
  const double imageCenterY = origin[1] + spacing[1] * size[1] / 2.0;
  const double imageCenterZ = origin[2] + spacing[2] * size[2] / 2.0;

  translation1[0] =   -imageCenterX;
  translation1[1] =   -imageCenterY;
  translation1[2] =   -imageCenterZ;
  
  transform->SetTranslation( translation1 );

  double pi = 3.14159265358979323846;

  const double degreesToRadians = pi/180;
  const double angle1 =par[3] * degreesToRadians; //par[3] = 
                                              rotation in x direction

  VersorType rotation_x;
  VersorType rotation_y;
  VersorType rotation_z;
 

 VectorType axis_x;
  axis_x[0] = 1; // Rotate around x-axis...
  axis_x[1] = 0;// Rotate around y-axis...
  axis_x[2] = 0;// Rotate around z-axis...
  
  rotation_x.Set(axis_x,angle1);
  transform->SetRotation(rotation_x);

  const double angle2 =par[4] * degreesToRadians; //Rotation in
y-direction

 VectorType axis_y;
  axis_y[0] = 0; // Rotate around x-axis...
  axis_y[1] = 1;// Rotate around y-axis...
  axis_y[2] = 0;// Rotate around z-axis...
  
  rotation_y.Set(axis_y,angle2);
  transform->SetRotation(rotation_y);

  const double angle3 =par[5] * degreesToRadians; // Rotation in z
                                                    direction
  VectorType axis_z;
  axis_z[0] = 0; // Rotate around x-axis...
  axis_z[1] = 0;// Rotate around y-axis...
  axis_z[2] = 1;// Rotate around z-axis...
  
  rotation_z.Set(axis_z,angle3);
  transform->SetRotation(rotation_z);

  TransformType::OutputVectorType translation2;
  translation2[0] =   imageCenterX;
  translation2[1] =   imageCenterY;
  translation2[2] =   imageCenterZ;	
  transform->SetTranslation( translation2); //Translate back

  filter->SetTransform( transform );
  filter->Update();
 // TransformType::MatrixType matrix = transform->GetRotationMatrix();
  TransformType::ParametersType param = transform->GetParameters();
  std::cout<<"Transformen: "<<param<<std::endl;
 // return filter->GetOutput();


What is wrong.. when i execute the code, it seems it is only performing
rotation in the z - direction. Because when I execute the result with
parameters: The rotation in all directions is 30 degrees.. and the
translation is 10 mm in each direction...
the result is:
[0,0,0,0.258819,185.984,185.984,176]

if I have understand it correctly:

the three first zero's are versors?? and 0.258819 is the total
rotation?? or what?? and the final three are the image center. Is this
correct?
Do I now have an image which is rotated in 30 degrees in every direction
and translated 10 mm in the same directions?? My second question is:
how to apply these parameters on the resampling filter?