[ITK] [ITK-users] QuaternionRigidTransformGradientDescentOptimizer and QuaternionRigidTransform

Matias Montroull matimontg at gmail.com
Sat May 2 11:32:13 EDT 2015


Hi,

I'm trying to do 3D registration using the Quaternion Classes and for some
reason I can't..

I get an error and I think it is because of the OptimizerScales.. I'm
unable to find an example that shows the use of Quaternions transforms so
I'm a little lost.

The error I get is at this line of code within
itkQuaternionRigidTransformGradientDescentOptimizer.cxx:

const unsigned int spaceDimension =
 m_CostFunction->GetNumberOfParameters();

#include "itkAffineTransform.h"
#include "itkQuaternionRigidTransform.h"
#include "itkQuaternionRigidTransformGradientDescentOptimizer.h"
#include "itkCenteredTransformInitializer.h"
#include "itkMultiResolutionImageRegistrationMethod.h"
#include "itkMattesMutualInformationImageToImageMetric.h"
#include "itkRegularStepGradientDescentOptimizer.h"
#include "itkRecursiveMultiResolutionPyramidImageFilter.h"
#include "itkImage.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkResampleImageFilter.h"
#include "itkCastImageFilter.h"
#include "itkNormalizeImageFilter.h"
#include "itkCheckerBoardImageFilter.h"

#include "itkCommand.h"
class CommandIterationUpdate : public itk::Command
{
public:
typedef  CommandIterationUpdate   Self;
typedef  itk::Command             Superclass;
typedef  itk::SmartPointer<Self>  Pointer;
itkNewMacro(Self);
protected:
CommandIterationUpdate() : m_CumulativeIterationIndex(0) {};
public:
typedef   itk::RegularStepGradientDescentOptimizer  OptimizerType;
typedef   const OptimizerType *                     OptimizerPointer;
void Execute(itk::Object *caller, const itk::EventObject & event)
ITK_OVERRIDE
{
Execute((const itk::Object *)caller, event);
}
void Execute(const itk::Object * object, const itk::EventObject & event)
ITK_OVERRIDE
{
OptimizerPointer optimizer = static_cast<OptimizerPointer>(object);
if (!(itk::IterationEvent().CheckEvent(&event)))
{
return;
}
std::cout << optimizer->GetCurrentIteration() << "   ";
std::cout << optimizer->GetValue() << "   ";
std::cout << optimizer->GetCurrentPosition() << "  " <<
m_CumulativeIterationIndex++ << std::endl;
}
private:
unsigned int m_CumulativeIterationIndex;
};

template <typename TRegistration>
class RegistrationInterfaceCommand : public itk::Command
{
public:
typedef  RegistrationInterfaceCommand   Self;
typedef  itk::Command                   Superclass;
typedef  itk::SmartPointer<Self>        Pointer;
itkNewMacro(Self);
protected:
RegistrationInterfaceCommand() {};
public:
typedef   TRegistration                              RegistrationType;
typedef   RegistrationType *                         RegistrationPointer;
typedef   itk::RegularStepGradientDescentOptimizer   OptimizerType;
typedef   OptimizerType *                            OptimizerPointer;
void Execute(itk::Object * object, const itk::EventObject & event)
ITK_OVERRIDE
{
if (!(itk::IterationEvent().CheckEvent(&event)))
{
return;
}
RegistrationPointer registration = static_cast<RegistrationPointer>(object);
OptimizerPointer optimizer =
static_cast<OptimizerPointer>(registration->GetModifiableOptimizer());
std::cout << "-------------------------------------" << std::endl;
std::cout << "MultiResolution Level : "
<< registration->GetCurrentLevel() << std::endl;
std::cout << std::endl;
if (registration->GetCurrentLevel() == 0)
{
optimizer->SetMaximumStepLength(1.00);
optimizer->SetMinimumStepLength(0.01);
}
else
{
optimizer->SetMaximumStepLength(optimizer->GetMaximumStepLength() / 4.0);
optimizer->SetMinimumStepLength(optimizer->GetMinimumStepLength() / 10.0);
}
}
void Execute(const itk::Object *, const itk::EventObject &) ITK_OVERRIDE
{ return; }
};
int main(int argc, char *argv[])
{
if (argc < 4)
{
std::cerr << "Missing Parameters " << std::endl;
std::cerr << "Usage: " << argv[0];
std::cerr << " imagenfija  imagenflotante ";
std::cerr << " salida";
return EXIT_FAILURE;
}
const    unsigned int    Dimension = 3;
typedef  signed short  PixelType;
typedef itk::Image< PixelType, Dimension >  FixedImageType;
typedef itk::Image< PixelType, Dimension >  MovingImageType;
typedef   float                                    InternalPixelType;
typedef itk::Image< InternalPixelType, Dimension > InternalImageType;

typedef itk::QuaternionRigidTransform< double> TransformType;

typedef itk::QuaternionRigidTransformGradientDescentOptimizer
OptimizerType;
typedef itk::LinearInterpolateImageFunction<
InternalImageType,
double             > InterpolatorType;
typedef itk::MattesMutualInformationImageToImageMetric<
InternalImageType,
InternalImageType >    MetricType;
typedef OptimizerType::ScalesType       OptimizerScalesType;
typedef itk::MultiResolutionImageRegistrationMethod<
InternalImageType,
InternalImageType    > RegistrationType;
typedef itk::ImageFileReader< FixedImageType  > FixedImageReaderType;
typedef itk::ImageFileReader< MovingImageType > MovingImageReaderType;
typedef itk::NormalizeImageFilter<FixedImageType, InternalImageType>
NormalizeFilterType;
typedef itk::DiscreteGaussianImageFilter<InternalImageType,
InternalImageType> GaussianFilterType;
typedef itk::CenteredTransformInitializer<TransformType, FixedImageType,
MovingImageType >  TransformInitializerType;

OptimizerType::Pointer      optimizer = OptimizerType::New();
InterpolatorType::Pointer   interpolator = InterpolatorType::New();
RegistrationType::Pointer   registration = RegistrationType::New();
MetricType::Pointer         metric = MetricType::New();
TransformType::Pointer   transform = TransformType::New();
FixedImageReaderType::Pointer  fixedImageReader =
FixedImageReaderType::New();
MovingImageReaderType::Pointer movingImageReader =
MovingImageReaderType::New();
NormalizeFilterType::Pointer fixedNormalizer = NormalizeFilterType::New();
NormalizeFilterType::Pointer movingNormalizer = NormalizeFilterType::New();
GaussianFilterType::Pointer fixedSmoother = GaussianFilterType::New();
GaussianFilterType::Pointer movingSmoother = GaussianFilterType::New();
TransformInitializerType::Pointer initializer =
TransformInitializerType::New();

fixedImageReader->SetFileName(argv[1]);
movingImageReader->SetFileName(argv[2]);

registration->SetOptimizer(optimizer);
registration->SetInterpolator(interpolator);
registration->SetMetric(metric);
registration->SetTransform(transform);

fixedNormalizer->SetInput(fixedImageReader->GetOutput());
movingNormalizer->SetInput(movingImageReader->GetOutput());

fixedSmoother->SetVariance(2.0);
movingSmoother->SetVariance(2.0);

fixedSmoother->SetInput(fixedNormalizer->GetOutput());
movingSmoother->SetInput(movingNormalizer->GetOutput());

registration->SetFixedImage(fixedSmoother->GetOutput());
registration->SetMovingImage(movingSmoother->GetOutput());

fixedNormalizer->Update();
registration->SetFixedImageRegion(fixedNormalizer->GetOutput()->GetBufferedRegion());

initializer->SetTransform(transform);
initializer->SetFixedImage(fixedImageReader->GetOutput());
initializer->SetMovingImage(movingImageReader->GetOutput());
initializer->MomentsOn();
initializer->InitializeTransform();
//transform->SetIdentity(); //Esto generó Too Many Samples outside mapping..

registration->SetInitialTransformParameters(transform->GetParameters());


OptimizerScalesType optimizerScales(transform->GetNumberOfParameters());

if (Dimension == 3)
{
optimizerScales[0] = 1.0;
optimizerScales[1] = 1.0;
optimizerScales[2] = 1.0;
optimizerScales[3] = 1.0;
optimizerScales[4] = 0.0001;
optimizerScales[5] = 0.0001;
optimizerScales[6] = 0.0001;
optimizerScales[7] = 0.0001;
}

optimizer->SetScales(optimizerScales);
 optimizer->SetNumberOfIterations(20);
optimizer->MaximizeOn();

metric->SetNumberOfSpatialSamples(50);

CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
optimizer->AddObserver(itk::IterationEvent(), observer);


typedef RegistrationInterfaceCommand<RegistrationType> CommandType;
CommandType::Pointer command = CommandType::New();
registration->AddObserver(itk::IterationEvent(), command);

try
{
registration->Update();
std::cout << "Optimizer stop condition: "
<< registration->GetOptimizer()->GetStopConditionDescription()
<< std::endl;
}
catch (itk::ExceptionObject & err)
{
std::cout << "ExceptionObject caught !" << std::endl;
std::cout << err << std::endl;
return EXIT_FAILURE;
}
std::cout << "Optimizer Stopping Condition = "
<< optimizer->GetStopCondition() << std::endl;
typedef RegistrationType::ParametersType ParametersType;
ParametersType finalParameters = registration->GetLastTransformParameters();
double TranslationAlongX = finalParameters[4];
double TranslationAlongY = finalParameters[5];
unsigned int numberOfIterations = optimizer->GetCurrentIteration();
double bestValue = optimizer->GetValue();

std::cout << "Result = " << std::endl;
std::cout << " Translation X = " << TranslationAlongX << std::endl;
std::cout << " Translation Y = " << TranslationAlongY << std::endl;
std::cout << " Iterations    = " << numberOfIterations << std::endl;
std::cout << " Metric value  = " << bestValue << std::endl;


typedef itk::ResampleImageFilter<
MovingImageType,
FixedImageType >    ResampleFilterType;
TransformType::Pointer finalTransform = TransformType::New();
finalTransform->SetParameters(finalParameters);
finalTransform->SetFixedParameters(transform->GetFixedParameters());
ResampleFilterType::Pointer resample = ResampleFilterType::New();
resample->SetTransform(finalTransform);
resample->SetInput(movingImageReader->GetOutput());
FixedImageType::Pointer fixedImage = fixedImageReader->GetOutput();
PixelType backgroundGrayLevel = 0;

resample->SetSize(fixedImage->GetLargestPossibleRegion().GetSize());
resample->SetOutputOrigin(fixedImage->GetOrigin());
resample->SetOutputSpacing(fixedImage->GetSpacing());
resample->SetOutputDirection(fixedImage->GetDirection());
resample->SetDefaultPixelValue(backgroundGrayLevel);
typedef  signed short                           OutputPixelType;
typedef itk::Image< OutputPixelType, Dimension > OutputImageType;
typedef itk::CastImageFilter<
FixedImageType,
OutputImageType >          CastFilterType;
typedef itk::ImageFileWriter< OutputImageType >  WriterType;
WriterType::Pointer      writer = WriterType::New();
CastFilterType::Pointer  caster = CastFilterType::New();
writer->SetFileName(argv[3]);
caster->SetInput(resample->GetOutput());
writer->SetInput(caster->GetOutput());
writer->Update();

typedef itk::CheckerBoardImageFilter< FixedImageType >
CheckerBoardFilterType;
CheckerBoardFilterType::Pointer checker = CheckerBoardFilterType::New();
checker->SetInput1(fixedImage);
checker->SetInput2(resample->GetOutput());
caster->SetInput(checker->GetOutput());
writer->SetInput(caster->GetOutput());
resample->SetDefaultPixelValue(0);

resample->SetTransform(finalTransform);
if (argc > 4)
{
writer->SetFileName(argv[11]);
writer->Update();
}
return EXIT_SUCCESS;
}
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://public.kitware.com/pipermail/community/attachments/20150502/fd0d0fd6/attachment-0001.html>
-------------- next part --------------
_____________________________________
Powered by www.kitware.com

Visit other Kitware open-source projects at
http://www.kitware.com/opensource/opensource.html

Kitware offers ITK Training Courses, for more information visit:
http://www.kitware.com/products/protraining.php

Please keep messages on-topic and check the ITK FAQ at:
http://www.itk.org/Wiki/ITK_FAQ

Follow this link to subscribe/unsubscribe:
http://public.kitware.com/mailman/listinfo/insight-users


More information about the Community mailing list