Skip to content

Commit

Permalink
BUG: Fix Fiducial Registration module rigid transform changing scale
Browse files Browse the repository at this point in the history 
When Rigid transform type was selected in Fiducial Registration module, the resulting transform included scaling (scale factor != 1.0) as reported at https://discourse.slicer.org/t/rigid-transformation-is-scaling/31953

The problem occurred because ITK landmark initializer has been improved so that it can now compute scaling. Fixed this by using more specific transform types (rigid/similarity).

Since ITK now also supports "affine" transform initialization, this option is now exposed in the module as well.
  • Loading branch information
@lassoan
lassoan committed Oct 5, 2023
1 parent 9d50b77 commit 06e48e5
Show file tree
Hide file tree
Showing 2 changed files with 72 additions and 106 deletions.
177 changes: 71 additions & 106 deletions Modules/CLI/FiducialRegistration/FiducialRegistration.cxx
View file
Expand Up @@ -7,6 +7,7 @@
#include <itkLandmarkBasedTransformInitializer.h>
#include <itkSimilarity3DTransform.h>
#include <itkTransformFileWriter.h>
#include <itkVersorRigid3DTransform.h>

// STD includes
#include <numeric>
Expand Down Expand Up @@ -43,49 +44,14 @@ namespace

};

// Function to compute the scaling factor between two sets of points.
// This is the symmetric form given by
// Berthold K. P. Horn (1987),
// "Closed-form solution of absolute orientation using unit quaternions,"
// Journal of the Optical Society of America A, 4:629-642

double
computeSymmetricScale(const std::vector<itk::Point<double, 3> >& fixedPoints,
const std::vector<itk::Point<double, 3> >& movingPoints,
const itk::Point<double, 3>& fixedcenter,
const itk::Point<double, 3>& movingcenter)
{
std::vector<double> centeredFixedPoints(fixedPoints.size(), 0.0);
std::vector<double> centeredMovingPoints(movingPoints.size(), 0.0);

std::transform(fixedPoints.begin(), fixedPoints.end(),
centeredFixedPoints.begin(),
SquaredPointDistance(fixedcenter) );

std::transform(movingPoints.begin(), movingPoints.end(),
centeredMovingPoints.begin(),
SquaredPointDistance(movingcenter) );

double fixedmag = 0.0, movingmag = 0.0;
fixedmag = std::accumulate(centeredFixedPoints.begin(),
centeredFixedPoints.end(),
fixedmag);

movingmag = std::accumulate(centeredMovingPoints.begin(),
centeredMovingPoints.end(),
movingmag);

return sqrt(movingmag / fixedmag);
}

}

int main(int argc, char* argv[])
{
PARSE_ARGS;


double invalidRMS = -1;
const double INVALID_RMS = -1;

// Checking conditions.

Expand All @@ -95,7 +61,7 @@ int main(int argc, char* argv[])
std::cerr << "Fixed and moving landmark lists must be of the same size "
<< "and contain at least one point" << std::endl;
outputMessage = "Fixed and moving landmark lists must be of the same size and contain at least one point";
rms = invalidRMS;
rms = INVALID_RMS;
}

if( saveTransform == "" )
Expand All @@ -104,7 +70,7 @@ int main(int argc, char* argv[])
if (outputMessage == "")
{
outputMessage = "An output transform must be specified";
rms = invalidRMS;
rms = INVALID_RMS;
}
}

Expand All @@ -117,13 +83,13 @@ int main(int argc, char* argv[])
if (outputMessage == "")
{
outputMessage = "At least 3 fixed landmark fiducial points must be specified for Rigid or Similarity transforms";
rms = invalidRMS;
rms = INVALID_RMS;
}
}

// Return if conditions not met.
// Return if conditions not met.

if ( rms == invalidRMS )
if ( rms == INVALID_RMS )
{
// Write out the return parameters in "name = value" form
std::ofstream rts;
Expand All @@ -135,8 +101,7 @@ int main(int argc, char* argv[])
return EXIT_SUCCESS;
}


// only calculate if the above conditions hold
// Prerequisites are fulfilled.

typedef std::vector<itk::Point<double, 3> > PointList;

Expand All @@ -153,94 +118,94 @@ int main(int argc, char* argv[])
movingPoints.begin(),
convertStdVectorToITKPoint);

// Our input into landmark based initialize will be of this form
// The format for saving to slicer is defined later
typedef itk::Similarity3DTransform<double> SimilarityTransformType;
SimilarityTransformType::Pointer transform = SimilarityTransformType::New();
transform->SetIdentity();
// workaround a bug in older versions of ITK
transform->SetScale(1.0);


typedef itk::LandmarkBasedTransformInitializer<SimilarityTransformType,
itk::Image<short, 3>, itk::Image<short, 3> > InitializerType;
InitializerType::Pointer initializer = InitializerType::New();

// This expects a VersorRigid3D. The similarity transform works because
// it derives from that class
initializer->SetTransform(transform);

initializer->SetFixedLandmarks(fixedPoints);
initializer->SetMovingLandmarks(movingPoints);

initializer->InitializeTransform();


// Handle different transform types.
typedef itk::AffineTransform<double, 3> AffineTransform;
AffineTransform::Pointer fixedToMovingT = itk::AffineTransform<double, 3>::New();

if( transformType == "Translation" )
if (transformType == "Translation" || transformType == "Rigid")
{
// Clear out the computed rotaitoin if we only requested translation
itk::Versor<double> v;
v.SetIdentity();
transform->SetRotation(v);
}
else if( transformType == "Rigid" )
{
// do nothing
typedef itk::VersorRigid3DTransform<double> TransformType;
TransformType::Pointer transform = TransformType::New();
transform->SetIdentity();

typedef itk::LandmarkBasedTransformInitializer<TransformType,
itk::Image<short, 3>, itk::Image<short, 3> > InitializerType;
InitializerType::Pointer initializer = InitializerType::New();
initializer->SetTransform(transform);
initializer->SetFixedLandmarks(fixedPoints);
initializer->SetMovingLandmarks(movingPoints);
initializer->InitializeTransform();

if (transformType == "Translation")
{
// Clear out the computed rotation if we only requested translation
itk::Versor<double> v;
v.SetIdentity();
transform->SetRotation(v);
}

// Convert into an affine transform for saving to Slicer.
fixedToMovingT->SetCenter(transform->GetCenter());
fixedToMovingT->SetMatrix(transform->GetMatrix());
fixedToMovingT->SetTranslation(transform->GetTranslation());
}
else if( transformType == "Similarity" )
{
// Compute the scaling factor and add that in
itk::Point<double, 3> fixedCenter(transform->GetCenter() );
itk::Point<double, 3> movingCenter(transform->GetCenter() + transform->GetTranslation() );
typedef itk::Similarity3DTransform<double> TransformType;
TransformType::Pointer transform = TransformType::New();
transform->SetIdentity();

typedef itk::LandmarkBasedTransformInitializer<TransformType,
itk::Image<short, 3>, itk::Image<short, 3> > InitializerType;
InitializerType::Pointer initializer = InitializerType::New();
initializer->SetTransform(transform);
initializer->SetFixedLandmarks(fixedPoints);
initializer->SetMovingLandmarks(movingPoints);
initializer->InitializeTransform();

double s = computeSymmetricScale(fixedPoints, movingPoints,
fixedCenter, movingCenter);
transform->SetScale(s);
// Convert into an affine transform for saving to Slicer.
fixedToMovingT->SetCenter(transform->GetCenter());
fixedToMovingT->SetMatrix(transform->GetMatrix());
fixedToMovingT->SetTranslation(transform->GetTranslation());
}
else if( transformType == "Affine" )
{
// itk::Matrix<double, 3> a = computeAffineTransform(fixedPoints, movingPoints, fixedCenter, movingCenter);
std::cerr << "Unsupported transform type: " << transformType << std::endl;
// return EXIT_FAILURE;
typedef itk::AffineTransform<double> TransformType;
TransformType::Pointer transform = TransformType::New();
transform->SetIdentity();

typedef itk::LandmarkBasedTransformInitializer<TransformType,
itk::Image<short, 3>, itk::Image<short, 3> > InitializerType;
InitializerType::Pointer initializer = InitializerType::New();
initializer->SetTransform(transform);
initializer->SetFixedLandmarks(fixedPoints);
initializer->SetMovingLandmarks(movingPoints);
initializer->InitializeTransform();

// Convert into an affine transform for saving to Slicer.
fixedToMovingT->SetCenter(transform->GetCenter());
fixedToMovingT->SetMatrix(transform->GetMatrix());
fixedToMovingT->SetTranslation(transform->GetTranslation());
}
else
{
std::cerr << "Unsupported transform type: " << transformType << std::endl;
return EXIT_FAILURE;
}


// Convert into an affine transform for saving to Slicer.

typedef itk::AffineTransform<double, 3> AffineTransform;
AffineTransform::Pointer fixedToMovingT = itk::AffineTransform<double, 3>::New();

fixedToMovingT->SetCenter( transform->GetCenter() );
fixedToMovingT->SetMatrix( transform->GetMatrix() );
fixedToMovingT->SetTranslation( transform->GetTranslation() );


// Compute RMS error in the target coordinate system.

AffineTransform::Pointer movingToFixedT = AffineTransform::New();
fixedToMovingT->GetInverse( movingToFixedT );

typedef InitializerType::LandmarkPointContainer LandmarkPointContainerType;

typedef PointList LandmarkPointContainerType;
typedef LandmarkPointContainerType::const_iterator PointsContainerConstIterator;
PointsContainerConstIterator mitr = movingPoints.begin();
PointsContainerConstIterator fitr = fixedPoints.begin();

SimilarityTransformType::OutputVectorType errortr;
SimilarityTransformType::OutputVectorType::RealValueType sum;
InitializerType::LandmarkPointType movingPointInFixed;
int counter;

sum = itk::NumericTraits< SimilarityTransformType::OutputVectorType::RealValueType >::ZeroValue();
counter = itk::NumericTraits< int >::ZeroValue();

AffineTransform::OutputVectorType::RealValueType sum = itk::NumericTraits< double >::ZeroValue();
AffineTransform::OutputVectorType errortr;
AffineTransform::OutputPointType movingPointInFixed;
int counter = itk::NumericTraits< int >::ZeroValue();
while( mitr != movingPoints.end() )
{
movingPointInFixed = movingToFixedT->TransformPoint( *mitr );
Expand Down
1 change: 1 addition & 0 deletions Modules/CLI/FiducialRegistration/FiducialRegistration.xml
View file
Expand Up @@ -42,6 +42,7 @@
<element>Translation</element>
<element>Rigid</element>
<element>Similarity</element>
<element>Affine</element>
<default>Rigid</default>
</string-enumeration>
<double>
Expand Down

0 comments on commit 06e48e5

Please sign in to comment.