Commit d308cfb8 authored by Floris Berendsen's avatar Floris Berendsen
Browse files

ENH: Split off ItkGaussianExponentialDiffeomorphicTransform

parent 57a3f660
/*=========================================================================
*
* Copyright Leiden University Medical Center, Erasmus University Medical
* Center and contributors
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0.txt
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*
*=========================================================================*/
#ifndef selxItkGaussianExponentialDiffeomorphicTransformComponent_h
#define selxItkGaussianExponentialDiffeomorphicTransformComponent_h
#include "ComponentBase.h"
#include "Interfaces.h"
#include "itkGaussianExponentialDiffeomorphicTransform.h"
#include <string.h>
#include "selxMacro.h"
namespace selx
{
template <class InternalComputationValueType, int Dimensionality>
class ItkGaussianExponentialDiffeomorphicTransformComponent :
public Implements<
Accepting< >,
Providing< itkTransformInterface<InternalComputationValueType,Dimensionality>>
>
{
public:
selxNewMacro(ItkGaussianExponentialDiffeomorphicTransformComponent, ComponentBase);
//itkStaticConstMacro(Dimensionality, unsigned int, Dimensionality);
ItkGaussianExponentialDiffeomorphicTransformComponent();
virtual ~ItkGaussianExponentialDiffeomorphicTransformComponent();
//typedef double InternalComputationValueType;
/** Type of the optimizer. */
using TransformType = typename itkTransformInterface<InternalComputationValueType,Dimensionality>::TransformType;
using TransformPointer = typename itkTransformInterface<InternalComputationValueType,Dimensionality>::TransformPointer;
typedef typename itk::GaussianExponentialDiffeomorphicTransform<InternalComputationValueType, Dimensionality> GaussianExponentialDiffeomorphicTransformType;
virtual TransformPointer GetItkTransform() override;
virtual bool MeetsCriterion(const ComponentBase::CriterionType &criterion) override;
//static const char * GetName() { return "ItkGaussianExponentialDiffeomorphicTransform"; } ;
static const char * GetDescription() { return "ItkGaussianExponentialDiffeomorphicTransform Component"; };
private:
typename GaussianExponentialDiffeomorphicTransformType::Pointer m_Transform;
protected:
/* The following struct returns the string name of computation type */
/* default implementation */
};
} //end namespace selx
#ifndef ITK_MANUAL_INSTANTIATION
#include "selxItkGaussianExponentialDiffeomorphicTransform.hxx"
#endif
#endif // #define selxItkGaussianExponentialDiffeomorphicTransformComponent_h
\ No newline at end of file
/*=========================================================================
*
* Copyright Leiden University Medical Center, Erasmus University Medical
* Center and contributors
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0.txt
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*
*=========================================================================*/
#include "selxItkGaussianExponentialDiffeomorphicTransform.h"
namespace selx
{
template<class InternalComputationValueType, int Dimensionality>
ItkGaussianExponentialDiffeomorphicTransformComponent< InternalComputationValueType,Dimensionality>::ItkGaussianExponentialDiffeomorphicTransformComponent()
{
m_Transform = GaussianExponentialDiffeomorphicTransformType::New();
typedef itk::Image<VectorType, Dimensionality> ConstantVelocityFieldType;
typename ConstantVelocityFieldType::Pointer displacementField = ConstantVelocityFieldType::New();
displacementField->CopyInformation(fixedImage);
displacementField->SetRegions(fixedImage->GetBufferedRegion());
displacementField->Allocate();
displacementField->FillBuffer(zeroVector);
typename ConstantVelocityFieldTransformType::Pointer fieldTransform = ConstantVelocityFieldTransformType::New();
fieldTransform->SetGaussianSmoothingVarianceForTheUpdateField(3.0);
fieldTransform->SetGaussianSmoothingVarianceForTheConstantVelocityField(6.0);
fieldTransform->SetConstantVelocityField(displacementField);
fieldTransform->SetCalculateNumberOfIntegrationStepsAutomatically(true);
fieldTransform->IntegrateVelocityField();
//TODO: instantiating the filter in the constructor might be heavy for the use in component selector factory, since all components of the database are created during the selection process.
// we could choose to keep the component light weighted (for checking criteria such as names and connections) until the settings are passed to the filter, but this requires an additional initialization step.
}
template<class InternalComputationValueType, int Dimensionality>
ItkGaussianExponentialDiffeomorphicTransformComponent< InternalComputationValueType, Dimensionality>::~ItkGaussianExponentialDiffeomorphicTransformComponent()
{
}
template<class InternalComputationValueType, int Dimensionality>
typename ItkGaussianExponentialDiffeomorphicTransformComponent< InternalComputationValueType, Dimensionality>::TransformPointer ItkGaussianExponentialDiffeomorphicTransformComponent< InternalComputationValueType, Dimensionality>::GetItkTransform()
{
return (TransformPointer) this->m_Transform;
}
template<class InternalComputationValueType, int Dimensionality>
bool
ItkGaussianExponentialDiffeomorphicTransformComponent< InternalComputationValueType, Dimensionality>
::MeetsCriterion(const ComponentBase::CriterionType &criterion)
{
bool hasUndefinedCriteria(false);
bool meetsCriteria(false);
if (criterion.first == "ComponentProperty")
{
meetsCriteria = true;
for (auto const & criterionValue : criterion.second) // auto&& preferred?
{
if (criterionValue != "SomeProperty") // e.g. "GradientDescent", "SupportsSparseSamples
{
meetsCriteria = false;
}
}
}
else if (criterion.first == "Dimensionality") //Supports this?
{
meetsCriteria = true;
for (auto const & criterionValue : criterion.second) // auto&& preferred?
{
if (std::stoi(criterionValue) != Dimensionality)
{
meetsCriteria = false;
}
}
}
return meetsCriteria;
}
} //end namespace selx
......@@ -42,6 +42,7 @@ namespace selx
public Implements<
Accepting< itkImageFixedInterface<Dimensionality, TPixel>,
itkImageMovingInterface<Dimensionality, TPixel>,
itkTransformInterface<double, Dimensionality>,
itkMetricv4Interface<Dimensionality, TPixel>,
itkOptimizerv4Interface<double>
>,
......@@ -63,14 +64,12 @@ namespace selx
// Get the type definitions from the interfaces
typedef typename itkImageFixedInterface<Dimensionality, TPixel>::ItkImageType FixedImageType;
typedef typename itkImageMovingInterface<Dimensionality, TPixel>::ItkImageType MovingImageType;
typedef typename itkTransformInterface<double, Dimensionality>::TransformType TransformType;
typedef typename itkTransformInterface<double, Dimensionality>::TransformPointer TransformPointer;
typedef typename itkOptimizerv4Interface<double>::InternalComputationValueType InternalComputationValueType;
// TODO for now we hard code the transform to be a stationary velocity field. See Set(*MetricInterface) for implementation
typedef double RealType;
typedef itk::GaussianExponentialDiffeomorphicTransform<RealType, Dimensionality> ConstantVelocityFieldTransformType;
//typedef itk::ImageRegistrationMethodv4<FixedImageType, MovingImageType, ConstantVelocityFieldTransformType> TheItkFilterType;
typedef typename itkOptimizerv4Interface<double>::InternalComputationValueType OptimizerInternalComputationValueType; //should be from class template
typedef typename itkTransformInterface<double, Dimensionality>::InternalComputationValueType TransformInternalComputationValueType; //should be from class template
typedef itk::ImageRegistrationMethodv4<FixedImageType, MovingImageType> TheItkFilterType;
typedef typename TheItkFilterType::ImageMetricType ImageMetricType;
typedef itk::RegistrationParameterScalesFromPhysicalShift<ImageMetricType> ScalesEstimatorType;
......@@ -78,8 +77,9 @@ namespace selx
//Accepting Interfaces:
virtual int Set(itkImageFixedInterface<Dimensionality, TPixel>*) override;
virtual int Set(itkImageMovingInterface<Dimensionality, TPixel>*) override;
virtual int Set(itkTransformInterface<TransformInternalComputationValueType, Dimensionality>*) override;
virtual int Set(itkMetricv4Interface<Dimensionality, TPixel>*) override;
virtual int Set(itkOptimizerv4Interface<InternalComputationValueType>*) override;
virtual int Set(itkOptimizerv4Interface<OptimizerInternalComputationValueType>*) override;
//Providing Interfaces:
virtual TransformPointer GetItkTransform() override;
......
......@@ -121,10 +121,6 @@ namespace selx
ItkImageRegistrationMethodv4Component< Dimensionality, TPixel>::ItkImageRegistrationMethodv4Component()
{
m_theItkFilter = TheItkFilterType::New();
typename ConstantVelocityFieldTransformType::Pointer transform = ConstantVelocityFieldTransformType::New();
m_theItkFilter->SetInitialTransform(transform);
m_theItkFilter->InPlaceOff();
//TODO: instantiating the filter in the constructor might be heavy for the use in component selector factory, since all components of the database are created during the selection process.
......@@ -157,6 +153,16 @@ int ItkImageRegistrationMethodv4Component< Dimensionality, TPixel>
return 0;
}
template<int Dimensionality, class TPixel>
int ItkImageRegistrationMethodv4Component< Dimensionality, TPixel>::Set(itkTransformInterface<TransformInternalComputationValueType, Dimensionality>* component)
{
this->m_theItkFilter->SetInitialTransform(component->GetItkTransform());
return 0;
}
template<int Dimensionality, class TPixel>
int ItkImageRegistrationMethodv4Component< Dimensionality, TPixel>::Set(itkMetricv4Interface<Dimensionality, TPixel>* component)
{
......@@ -166,7 +172,7 @@ int ItkImageRegistrationMethodv4Component< Dimensionality, TPixel>::Set(itkMetri
}
template<int Dimensionality, class TPixel>
int ItkImageRegistrationMethodv4Component< Dimensionality, TPixel>::Set(itkOptimizerv4Interface<InternalComputationValueType>* component)
int ItkImageRegistrationMethodv4Component< Dimensionality, TPixel>::Set(itkOptimizerv4Interface<OptimizerInternalComputationValueType>* component)
{
this->m_theItkFilter->SetOptimizer(component->GetItkOptimizerv4());
......@@ -239,39 +245,7 @@ void ItkImageRegistrationMethodv4Component< Dimensionality, TPixel>::RunRegistra
this->m_theItkFilter->SetOptimizer(optimizer);
// TODO: for now we hard code the transform to be a stationary velocity field. See template declaration.
//typedef itk::CompositeTransform<RealType, Dimensionality> CompositeTransformType;
//typename CompositeTransformType::Pointer compositeTransform = CompositeTransformType::New();
//typedef itk::IdentityTransform < RealType, Dimensionality> IdentityTransformType;
//typename IdentityTransformType::Pointer idTransform = IdentityTransformType::New();
//compositeTransform->AddTransform(idTransform);
typedef itk::Vector<RealType, Dimensionality> VectorType;
VectorType zeroVector(0.0);
typedef itk::Image<VectorType, Dimensionality> DisplacementFieldType;
typedef itk::Image<VectorType, Dimensionality> ConstantVelocityFieldType;
typename ConstantVelocityFieldType::Pointer displacementField = ConstantVelocityFieldType::New();
displacementField->CopyInformation(fixedImage);
displacementField->SetRegions(fixedImage->GetBufferedRegion());
displacementField->Allocate();
displacementField->FillBuffer(zeroVector);
typename ConstantVelocityFieldTransformType::Pointer fieldTransform = ConstantVelocityFieldTransformType::New();
//fieldTransform->SetGaussianSmoothingVarianceForTheUpdateField(0.75);
//fieldTransform->SetGaussianSmoothingVarianceForTheConstantVelocityField(1.5);
fieldTransform->SetGaussianSmoothingVarianceForTheUpdateField(3.0);
fieldTransform->SetGaussianSmoothingVarianceForTheConstantVelocityField(6.0);
fieldTransform->SetConstantVelocityField(displacementField);
fieldTransform->SetCalculateNumberOfIntegrationStepsAutomatically(true);
fieldTransform->IntegrateVelocityField();
//this->m_theItkFilter->SetMovingInitialTransform(compositeTransform);
//this->m_theItkFilter->SetMovingInitialTransform(idTransform);
this->m_theItkFilter->SetNumberOfLevels(3);
// Shrink the virtual domain by specified factors for each level. See documentation
......@@ -292,6 +266,11 @@ void ItkImageRegistrationMethodv4Component< Dimensionality, TPixel>::RunRegistra
smoothingSigmasPerLevel[2] = 1;
this->m_theItkFilter->SetSmoothingSigmasPerLevel(smoothingSigmasPerLevel);
// TODO for now we hard code the transform to be a stationary velocity field.
typedef double RealType;
typedef itk::GaussianExponentialDiffeomorphicTransform<RealType, Dimensionality> ConstantVelocityFieldTransformType;
typedef ConstantVelocityFieldTransformType::ConstantVelocityFieldType ConstantVelocityFieldType;
typedef itk::GaussianExponentialDiffeomorphicTransformParametersAdaptor<ConstantVelocityFieldTransformType> VelocityFieldTransformAdaptorType;
typename TheItkFilterType::TransformParametersAdaptorsContainerType adaptors;
......@@ -302,10 +281,10 @@ void ItkImageRegistrationMethodv4Component< Dimensionality, TPixel>::RunRegistra
// domain at each level. To speed up calculation and avoid unnecessary memory
// usage, we could calculate these fixed parameters directly.
typedef itk::ShrinkImageFilter<ConstantVelocityFieldType, ConstantVelocityFieldType> ShrinkFilterType;
typedef itk::ShrinkImageFilter<FixedImageType, FixedImageType> ShrinkFilterType;
typename ShrinkFilterType::Pointer shrinkFilter = ShrinkFilterType::New();
shrinkFilter->SetShrinkFactors(shrinkFactorsPerLevel[level]);
shrinkFilter->SetInput(fieldTransform->GetConstantVelocityField());
shrinkFilter->SetInput(fixedImage);
shrinkFilter->Update();
typename VelocityFieldTransformAdaptorType::Pointer fieldTransformAdaptor = VelocityFieldTransformAdaptorType::New();
......@@ -328,13 +307,10 @@ void ItkImageRegistrationMethodv4Component< Dimensionality, TPixel>::RunRegistra
*/
this->m_theItkFilter->SetTransformParametersAdaptorsPerLevel(adaptors);
this->m_theItkFilter->SetInitialTransform(fieldTransform);
this->m_theItkFilter->InPlaceOn();
typedef CommandIterationUpdate<TheItkFilterType> DisplacementFieldRegistrationCommandType;
typename DisplacementFieldRegistrationCommandType::Pointer displacementFieldObserver = DisplacementFieldRegistrationCommandType::New();
this->m_theItkFilter->AddObserver(itk::IterationEvent(), displacementFieldObserver);
typedef CommandIterationUpdate<TheItkFilterType> RegistrationCommandType;
typename RegistrationCommandType::Pointer registrationObserver = RegistrationCommandType::New();
this->m_theItkFilter->AddObserver(itk::IterationEvent(), registrationObserver);
// perform the actual registration
this->m_theItkFilter->Update();
......
/*=========================================================================
*
* Copyright Leiden University Medical Center, Erasmus University Medical
* Center and contributors
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0.txt
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*
*=========================================================================*/
#ifndef selxItkResampleFilter_h
#define selxItkResampleFilter_h
#include "ComponentBase.h"
#include "Interfaces.h"
#include "itkImageRegistrationMethodv4.h"
#include "itkGradientDescentOptimizerv4.h"
#include "itkImageSource.h"
#include <itkTransformToDisplacementFieldFilter.h>
#include <string.h>
#include "selxMacro.h"
#include "itkComposeDisplacementFieldsImageFilter.h"
#include "itkGaussianExponentialDiffeomorphicTransform.h"
#include "itkGaussianExponentialDiffeomorphicTransformParametersAdaptor.h"
namespace selx
{
template <int Dimensionality, class TPixel, class TInternalComputationValue>
class ItkResampleFilterComponent :
public Implements<
Accepting< itkTransformInterface<TInternalComputationValue, Dimensionality>,
itkImageFixedInterface<Dimensionality, TPixel>, //TODO should be FixedImageDomainInterface, we do not require intensities
itkImageMovingInterface<Dimensionality, TPixel>
>,
Providing< itkImageInterface<Dimensionality, TPixel>
>
>
{
public:
selxNewMacro(ItkResampleFilterComponent, ComponentBase);
//itkStaticConstMacro(Dimensionality, unsigned int, Dimensionality);
ItkResampleFilterComponent();
virtual ~ItkResampleFilterComponent();
typedef TPixel PixelType;
// Get the type definitions from the interfaces
typedef typename itkImageFixedInterface<Dimensionality, TPixel>::ItkImageType FixedImageType;
typedef typename itkImageMovingInterface<Dimensionality, TPixel>::ItkImageType MovingImageType;
typedef typename itkImageInterface<Dimensionality, TPixel>::ItkImageType ResultImageType;
typedef itk::ResampleImageFilter<MovingImageType, ResultImageType> ResampleFilterType;
//Accepting Interfaces:
virtual int Set(itkImageFixedInterface<Dimensionality, TPixel>*) override;
virtual int Set(itkImageMovingInterface<Dimensionality, TPixel>*) override;
virtual int Set(itkTransformInterface<TInternalComputationValue, Dimensionality>*) override;
//Providing Interfaces:
virtual typename ResultImageType::Pointer GetItkImage() override;
//BaseClass methods
virtual bool MeetsCriterion(const ComponentBase::CriterionType &criterion) override;
//static const char * GetName() { return "ItkResampleFilter"; } ;
static const char * GetDescription() { return "ItkResampleFilter Component"; };
private:
typename ResampleFilterType::Pointer m_ResampleFilter;
protected:
/* The following struct returns the string name of computation type */
/* default implementation */
static inline const std::string GetTypeNameString()
{
itkGenericExceptionMacro(<< "Unknown ScalarType" << typeid(TPixel).name());
// TODO: provide the user instructions how to enable the compilation of the component with the required template types (if desired)
// We might define an exception object that can communicate various error messages: for simple user, for developer user, etc
}
static inline const std::string GetPixelTypeNameString()
{
itkGenericExceptionMacro(<< "Unknown PixelType" << typeid(TPixel).name());
// TODO: provide the user instructions how to enable the compilation of the component with the required template types (if desired)
// We might define an exception object that can communicate various error messages: for simple user, for developer user, etc
}
};
// unfortunately partial specialization of member functions is not allowed, without partially specializing the entire class.
/*
template <int Dimensionality>
class ItkResampleFilterComponent < Dimensionality, double >
{
static inline const std::string GetPixelTypeNameString();
};
template <int Dimensionality>
inline const std::string
ItkResampleFilterComponent<Dimensionality, double>
::GetPixelTypeNameString()
{
return std::string("double");
}
*/
template <>
inline const std::string
ItkResampleFilterComponent<2, float, double>
::GetPixelTypeNameString()
{
return std::string("float");
}
template <>
inline const std::string
ItkResampleFilterComponent<2, double, double>
::GetPixelTypeNameString()
{
return std::string("double");
}
template <>
inline const std::string
ItkResampleFilterComponent<3, float, double>
::GetPixelTypeNameString()
{
return std::string("float");
}
template <>
inline const std::string
ItkResampleFilterComponent<3, double, double>
::GetPixelTypeNameString()
{
return std::string("double");
}
template <>
inline const std::string
ItkResampleFilterComponent<2, float, double>
::GetTypeNameString()
{
return std::string("2_float");
}
template <>
inline const std::string
ItkResampleFilterComponent<2, double, double>
::GetTypeNameString()
{
return std::string("2_double");
}
template <>
inline const std::string
ItkResampleFilterComponent<3, float, double>
::GetTypeNameString()
{
return std::string("3_float");
}
template <>
inline const std::string
ItkResampleFilterComponent<3, double, double>
::GetTypeNameString()
{
return std::string("3_double");
}
} //end namespace selx
#ifndef ITK_MANUAL_INSTANTIATION
#include "selxItkResampleFilter.hxx"
#endif
#endif // #define GDOptimizer3rdPartyComponent_h
\ No newline at end of file
/*=========================================================================
*
* Copyright Leiden University Medical Center, Erasmus University Medical
* Center and contributors
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0.txt
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*
*=========================================================================*/
#include "selxItkResampleFilter.h"
namespace selx
{
template <int Dimensionality, class TPixel, class TInternalComputationValue>
ItkResampleFilterComponent< Dimensionality, TPixel, TInternalComputationValue>::ItkResampleFilterComponent()
{
m_ResampleFilter = ResampleFilterType::New();
//TODO: instantiating the filter in the constructor might be heavy for the use in component selector factory, since all components of the database are created during the selection process.
// we could choose to keep the component light weighted (for checking criteria such as names and connections) until the settings are passed to the filter, but this requires an additional initialization step.
}
template <int Dimensionality, class TPixel, class TInternalComputationValue>
ItkResampleFilterComponent< Dimensionality, TPixel, TInternalComputationValue>::~ItkResampleFilterComponent()
{
}
template <int Dimensionality, class TPixel, class TInternalComputationValue>
int ItkResampleFilterComponent< Dimensionality, TPixel, TInternalComputationValue>
::Set(itkImageFixedInterface<Dimensionality, TPixel>* component)
{
auto fixedImage = component->GetItkImageFixed();
// connect the itk pipeline
//this->m_ResampleFilter->SetSize(fixedImage->GetBufferedRegion().GetSize()); //should be virtual image...
this->m_ResampleFilter->SetSize(fixedImage->GetLargestPossibleRegion().GetSize()); //should be virtual image...
this->m_ResampleFilter->SetOutputOrigin(fixedImage->GetOrigin());
this->m_ResampleFilter->SetOutputSpacing(fixedImage->GetSpacing());
this->m_ResampleFilter->SetOutputDirection(fixedImage->GetDirection());
this->m_ResampleFilter->SetDefaultPixelValue(0);
return 0;
}
template <int Dimensionality, class TPixel, class TInternalComputationValue>
int ItkResampleFilterComponent< Dimensionality, TPixel, TInternalComputationValue>
::Set(itkImageMovingInterface<Dimensionality, TPixel>* component)
{
auto movingImage = component->GetItkImageMoving();
// connect the itk pipeline
this->m_ResampleFilter->SetInput(movingImage);
return 0;