[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/insight-users/attachments/20150502/fd0d0fd6/attachment.html>
More information about the Insight-users
mailing list