diff --git a/Components/Transforms/AdvancedLimitedEulerTransform/CMakeLists.txt b/Components/Transforms/AdvancedLimitedEulerTransform/CMakeLists.txt new file mode 100644 index 000000000..9f71cc264 --- /dev/null +++ b/Components/Transforms/AdvancedLimitedEulerTransform/CMakeLists.txt @@ -0,0 +1,9 @@ + +ADD_ELXCOMPONENT( AdvancedLimitedEulerTransformElastix + itkLimitedEulerTransform.h + itkAdvancedLimitedEuler3DTransform.h + itkAdvancedLimitedEuler3DTransform.hxx + elxAdvancedLimitedEulerTransform.h + elxAdvancedLimitedEulerTransform.hxx + elxAdvancedLimitedEulerTransform.cxx ) + diff --git a/Components/Transforms/AdvancedLimitedEulerTransform/elxAdvancedLimitedEulerTransform.cxx b/Components/Transforms/AdvancedLimitedEulerTransform/elxAdvancedLimitedEulerTransform.cxx new file mode 100644 index 000000000..128406e7c --- /dev/null +++ b/Components/Transforms/AdvancedLimitedEulerTransform/elxAdvancedLimitedEulerTransform.cxx @@ -0,0 +1,21 @@ +/*========================================================================= + * + * Copyright UMC Utrecht 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 "elxAdvancedLimitedEulerTransform.h" + +elxInstallMacro( AdvancedLimitedEulerTransformElastix ); diff --git a/Components/Transforms/AdvancedLimitedEulerTransform/elxAdvancedLimitedEulerTransform.h b/Components/Transforms/AdvancedLimitedEulerTransform/elxAdvancedLimitedEulerTransform.h new file mode 100644 index 000000000..126a4ba67 --- /dev/null +++ b/Components/Transforms/AdvancedLimitedEulerTransform/elxAdvancedLimitedEulerTransform.h @@ -0,0 +1,267 @@ +/*========================================================================= + * + * Copyright UMC Utrecht 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 __elxAdvancedLimitedEulerTransform_H__ +#define __elxAdvancedLimitedEulerTransform_H__ + +#include "elxIncludes.h" // include first to avoid MSVS warning +#include "itkAdvancedCombinationTransform.h" +#include "itkLimitedEulerTransform.h" +#include "itkCenteredTransformInitializer.h" + +namespace elastix +{ + +/** + * \class AdvancedLimitedEulerTransformElastix + * \brief A transform based on the itk LimitedEulerTransforms. + * + * This transform is a rigid body transformation. + * + * The parameters used in this class are: + * \parameter Transform: Select this transform as follows:\n + * (%Transform "LimitedEulerTransform") + * \parameter UpperLimits: the upper limit value for the rotations and translations, + * used in the optimizer. \n + * example: (UpperLimits 3.14 3.14 3.14 100.00 100.00 100.00) \n + * The number of arguments should be equal to the number of parameters: + * for each parameter its upper limit value. Arguments are in the order as for + * the Euler transform parameters: . Set the limit + * value to extreme positive values to (e.g. 1e+3) to effectively disable the upper limits. + * \parameter LowerLimits: the lower limit values for the rotations and translations, + * used in the optimizer. \n + * example: (LowerLimits -3.14 -3.14 -3.14 -100.00 -100.00 -100.00) \n + * The number of arguments should be equal to the number of parameters: + * for each parameter its lower limit value. Arguments are in the order as for + * the Euler transform parameters: . Set the limit + * value to extreme negative values to (e.g. -1e+3) to effectively disable the upper limits. + * \parameter SharpnessOfLimits: the scale factor which determines the limit cutoff sharpness. \n + * example: (LimitSharpness 25.0) \n + * High values impose a sharp limit (e.g. 10.0--100.0), while lower values (1.0-10.0). + * Minimum value is 1.0, lower values will be clamped to 1.0. If this parameter option + * is not used, by default for the sharpness limit value is 25.0.* + * \parameter Scales: the scale factor between the rotations and translations, + * used in the optimizer. \n + * example: (Scales 200000.0) \n + * example: (Scales 100000.0 60000.0 ... 80000.0) \n + * If only one argument is given, that factor is used for the rotations. + * If more than one argument is given, then the number of arguments should be + * equal to the number of parameters: for each parameter its scale factor. + * If this parameter option is not used, by default the rotations are scaled + * by a factor of 100000.0. See also the AutomaticScalesEstimation parameter. + * \parameter AutomaticScalesEstimation: if this parameter is set to "true" the Scales + * parameter is ignored and the scales are determined automatically. \n + * example: ( AutomaticScalesEstimation "true" ) \n + * Default: "false" (for backwards compatibility). Recommended: "true". + * \parameter CenterOfRotation: an index around which the image is rotated. \n + * example: (CenterOfRotation 128 128 90) \n + * By default the CenterOfRotation is set to the geometric center of the image. + * \parameter AutomaticTransformInitialization: whether or not the initial translation + * between images should be estimated as the distance between their centers.\n + * example: (AutomaticTransformInitialization "true") \n + * By default "false" is assumed. So, no initial translation. + * \parameter AutomaticTransformInitializationMethod: how to initialize this + * transform. Should be one of {GeometricalCenter, CenterOfGravity}.\n + * example: (AutomaticTransformInitializationMethod "CenterOfGravity") \n + * By default "GeometricalCenter" is assumed.\n + * + * The transform parameters necessary for transformix, additionally defined by this class, are: + * \transformparameter CenterOfRotation: stores the center of rotation as an index. \n + * example: (CenterOfRotation 128 128 90)\n + * depecrated! From elastix version 3.402 this is changed to CenterOfRotationPoint! + * \transformparameter CenterOfRotationPoint: stores the center of rotation, expressed in world coordinates. \n + * example: (CenterOfRotationPoint 10.555 6.666 12.345) + * + * \ingroup Transforms + */ + +template< class TElastix > +class AdvancedLimitedEulerTransformElastix : + public itk::AdvancedCombinationTransform< + typename elx::TransformBase< TElastix >::CoordRepType, + elx::TransformBase< TElastix >::FixedImageDimension >, + public elx::TransformBase< TElastix > +{ +public: + + /** Standard ITK-stuff.*/ + typedef AdvancedLimitedEulerTransformElastix Self; + + typedef itk::AdvancedCombinationTransform< + typename elx::TransformBase< TElastix >::CoordRepType, + elx::TransformBase< TElastix >::FixedImageDimension > Superclass1; + + typedef elx::TransformBase< TElastix > Superclass2; + + /** The ITK-class that provides most of the functionality, and + * that is set as the "CurrentTransform" in the CombinationTransform */ + typedef itk::LimitedEulerTransform< + typename elx::TransformBase< TElastix >::CoordRepType, + elx::TransformBase< TElastix >::FixedImageDimension > LimitedEulerTransformType; + + // typedef itk::LimitedEuler3DTransform< + // typename elx::TransformBase< TElastix >::CoordRepType > LimitedEulerTransformType; + + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + /** Method for creation through the object factory. */ + itkNewMacro( Self ); + + /** Run-time type information (and related methods). */ + //itkTypeMacro( AdvancedLimitedEulerTransformElastix, LimitedEulerTransform ); + itkTypeMacro( AdvancedLimitedEulerTransformElastix, itk::AdvancedCombinationTransform ); + + /** Name of this class. + * Use this name in the parameter file to select this specific transform. \n + * example: (Transform "AdvancedLimitedEulerTransform")\n + */ + elxClassNameMacro( "AdvancedLimitedEulerTransform" ); + + /** Dimension of the fixed image. */ + itkStaticConstMacro( SpaceDimension, unsigned int, Superclass2::FixedImageDimension ); + + /** Typedefs inherited from the superclass. */ + + /** These are both in Euler2D and Euler3D. */ + typedef typename Superclass1::ScalarType ScalarType; + typedef typename Superclass1::ParametersType ParametersType; + typedef typename Superclass1::NumberOfParametersType NumberOfParametersType; + typedef typename Superclass1::JacobianType JacobianType; + + typedef typename Superclass1::InputPointType InputPointType; + typedef typename Superclass1::OutputPointType OutputPointType; + typedef typename Superclass1::InputVectorType InputVectorType; + typedef typename Superclass1::OutputVectorType OutputVectorType; + typedef typename Superclass1::InputCovariantVectorType InputCovariantVectorType; + typedef typename Superclass1::OutputCovariantVectorType OutputCovariantVectorType; + typedef typename Superclass1::InputVnlVectorType InputVnlVectorType; + typedef typename Superclass1::OutputVnlVectorType OutputVnlVectorType; + + typedef typename LimitedEulerTransformType::Pointer LimitedEulerTransformPointer; + typedef typename LimitedEulerTransformType::OffsetType OffsetType; + + /** Typedef's inherited from TransformBase. */ + typedef typename Superclass2::ElastixType ElastixType; + typedef typename Superclass2::ElastixPointer ElastixPointer; + typedef typename Superclass2::ParameterMapType ParameterMapType; + typedef typename Superclass2::ConfigurationType ConfigurationType; + typedef typename Superclass2::ConfigurationPointer ConfigurationPointer; + typedef typename Superclass2::RegistrationType RegistrationType; + typedef typename Superclass2::RegistrationPointer RegistrationPointer; + typedef typename Superclass2::CoordRepType CoordRepType; + typedef typename Superclass2::FixedImageType FixedImageType; + typedef typename Superclass2::MovingImageType MovingImageType; + typedef typename Superclass2::ITKBaseType ITKBaseType; + typedef typename Superclass2::CombinationTransformType CombinationTransformType; + + /** Other typedef's. */ + typedef typename FixedImageType::IndexType IndexType; + typedef typename IndexType::IndexValueType IndexValueType; + typedef typename FixedImageType::SizeType SizeType; + typedef typename FixedImageType::PointType PointType; + typedef typename FixedImageType::SpacingType SpacingType; + typedef typename FixedImageType::RegionType RegionType; + typedef typename FixedImageType::DirectionType DirectionType; + + typedef itk::CenteredTransformInitializer< + LimitedEulerTransformType, FixedImageType, MovingImageType > TransformInitializerType; + typedef typename TransformInitializerType::Pointer TransformInitializerPointer; + + /** For scales setting in the optimizer */ + typedef typename Superclass2::ScalesType ScalesType; + + /** Execute stuff before the actual registration: + * \li Call InitializeTransform + * \li Set the scales. + */ + void BeforeRegistration( void ) override; + + /** Initialize Transform. + * \li Set all parameters to zero. + * \li Set center of rotation: + * automatically initialized to the geometric center of the image, or + * assigned a user entered voxel index, given by the parameter + * (CenterOfRotation ...); + * If an initial transform is present and HowToCombineTransforms is + * set to "Compose", the initial transform is taken into account + * while setting the center of rotation. + * \li Set initial translation: + * the initial translation between fixed and moving image is guessed, + * if the user has set (AutomaticTransformInitialization "true"). + * + * It is not yet possible to enter an initial rotation angle. + */ + virtual void InitializeTransform( void ); + + /** Set the scales + * \li If AutomaticScalesEstimation is "true" estimate scales + * \li If scales are provided by the user use those, + * \li Otherwise use some default value + * This function is called by BeforeRegistration, after + * the InitializeTransform function is called + */ + virtual void SetScales( void ); + + /** Function to read transform-parameters from a file. + * + * It reads the center of rotation and calls the superclass' implementation. + */ + void ReadFromFile( void ) override; + + /** Function to create transform-parameters map. + * Creates the TransformParametersmap + */ + ParameterMapType + CreateDerivedTransformParametersMap(void) const override; + + void + AfterRegistration(void) override; + +protected: + + /** The constructor. */ + AdvancedLimitedEulerTransformElastix(); + /** The destructor. */ + ~AdvancedLimitedEulerTransformElastix() override {} + + /** Try to read the CenterOfRotationPoint from the transform parameter file + * The CenterOfRotationPoint is already in world coordinates. + * Transform parameter files generated by elastix version > 3.402 + * save the center of rotation in this way. + */ + virtual bool ReadCenterOfRotationPoint( InputPointType & rotationPoint ) const; + +private: + elxOverrideGetSelfMacro; + + /** The private constructor. */ + AdvancedLimitedEulerTransformElastix( const Self & ); // purposely not implemented + /** The private copy constructor. */ + void operator=( const Self & ); // purposely not implemented + + LimitedEulerTransformPointer m_LimitedEulerTransform; + +}; + +} // end namespace elastix + +#ifndef ITK_MANUAL_INSTANTIATION +#include "elxAdvancedLimitedEulerTransform.hxx" +#endif + +#endif // end #ifndef __elxAdvancedLimitedEulerTransform_H__ diff --git a/Components/Transforms/AdvancedLimitedEulerTransform/elxAdvancedLimitedEulerTransform.hxx b/Components/Transforms/AdvancedLimitedEulerTransform/elxAdvancedLimitedEulerTransform.hxx new file mode 100644 index 000000000..d9c87c9ac --- /dev/null +++ b/Components/Transforms/AdvancedLimitedEulerTransform/elxAdvancedLimitedEulerTransform.hxx @@ -0,0 +1,588 @@ +/*========================================================================= + * + * Copyright UMC Utrecht 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 __elxAdvancedLimitedEulerTransform_HXX_ +#define __elxAdvancedLimitedEulerTransform_HXX_ + +#include "elxAdvancedLimitedEulerTransform.h" +#include "itkImageGridSampler.h" +#include "itkContinuousIndex.h" + +namespace elastix +{ + +/** + * ********************* Constructor **************************** + */ + +template< class TElastix > +AdvancedLimitedEulerTransformElastix< TElastix > +::AdvancedLimitedEulerTransformElastix() +{ + this->m_LimitedEulerTransform = LimitedEulerTransformType::New(); + this->SetCurrentTransform( this->m_LimitedEulerTransform ); + +} // end Constructor + + +/** + * ******************* BeforeRegistration *********************** + */ + +template< class TElastix > +void +AdvancedLimitedEulerTransformElastix< TElastix > +::BeforeRegistration( void ) +{ + /** Set center of rotation and initial translation. */ + this->InitializeTransform(); + + /** Set the scales. */ + this->m_LimitedEulerTransform->SetScalesEstimation(true); + this->SetScales(); + this->m_LimitedEulerTransform->SetScalesEstimation(false); + +} // end BeforeRegistration() + + +/** + * ************************* ReadFromFile ************************ + */ + +template< class TElastix > +void +AdvancedLimitedEulerTransformElastix< TElastix > +::ReadFromFile( void ) +{ + /** Variables. */ + InputPointType centerOfRotationPoint; + centerOfRotationPoint.Fill( 0.0 ); + + /** Try first to read the CenterOfRotationPoint from the + * transform parameter file, this is the new, and preferred + * way, since elastix 3.402. + */ + bool pointRead = this->ReadCenterOfRotationPoint( centerOfRotationPoint ); + + if( !pointRead ) + { + xl::xout[ "error" ] << "ERROR: No center of rotation is specified in " + << "the transform parameter file" << std::endl; + itkExceptionMacro( << "Transform parameter file is corrupt." ) + } + + /** Set the center in this Transform. */ + this->m_LimitedEulerTransform->SetCenter( centerOfRotationPoint ); + + /** Read the ComputeZYX. */ + if( SpaceDimension == 3 ) + { + std::string computeZYX = "false"; + this->m_Configuration->ReadParameter( computeZYX, "ComputeZYX", 0 ); + if( computeZYX == "true" ) + { + this->m_LimitedEulerTransform->SetComputeZYX( true ); + } + } + + /** Read SharpnessOfLimits parameter */ + ScalarType sharpnessOfLimits = 50.0f; + this->m_Configuration->ReadParameter( sharpnessOfLimits, "SharpnessOfLimits", 0 ); + this->m_LimitedEulerTransform->SetSharpnessOfLimits( sharpnessOfLimits ); + + /** Read and set UpperLimits parameter */ + ParametersType upperLimits; + upperLimits.SetSize(6); + for( unsigned int i = 0; i < 6; i++ ) + { + upperLimits[ i ] = i < 3 ? 3.14f : 200.0f; + + /** Returns zero when parameter was in the parameter file */ + bool found = this->m_Configuration->ReadParameter( + upperLimits[ i ], "UpperLimits", i, false ); + if( !found ) + { + for( unsigned int i = 0; i < 6; i++ ) + { + upperLimits[ i ] = i < 3 ? 3.14f : 200.0f; + } + xl::xout[ "warning" ] << "WARNING: Upper limits should contain six values! " + << "Using default values instead!" << std::endl; + break; + } + } + this->m_LimitedEulerTransform->SetUpperLimits( upperLimits ); + + /** Read and set LowerLimits parameter */ + ParametersType lowerLimits; + lowerLimits.SetSize(6); + for( unsigned int i = 0; i < 6; i++ ) + { + lowerLimits[ i ] = i < 3 ? -3.14f : -200.0f; + + /** Returns zero when parameter was in the parameter file */ + bool found = this->m_Configuration->ReadParameter( + lowerLimits[ i ], "LowerLimits", i, false ); + if( !found ) + { + for( unsigned int i = 0; i < 6; i++ ) + { + lowerLimits[ i ] = i < 3 ? -3.14f : -200.0f; + } + xl::xout[ "warning" ] << "WARNING: Lower limits should contain six values! " + << "Using default values instead!" << std::endl; + break; + } + } + this->m_LimitedEulerTransform->SetLowerLimits( lowerLimits ); + + /** Call the ReadFromFile from the TransformBase. + * BE AWARE: Only call Superclass2::ReadFromFile() after CenterOfRotation + * is set, because it is used in the SetParameters()-function of this transform. + */ + this->Superclass2::ReadFromFile(); + +} // end ReadFromFile() + + +/** + * ************************* CreateDerivedTransformParametersMap ************************ + */ +template +auto +AdvancedLimitedEulerTransformElastix::CreateDerivedTransformParametersMap(void) const -> ParameterMapType +{ + ParameterMapType parameterMap{ + { "CenterOfRotationPoint", Conversion::ToVectorOfStrings(m_LimitedEulerTransform->GetCenter()) }, + { "SharpnessOfLimits", { Conversion::ToString(m_LimitedEulerTransform->GetSharpnessOfLimits()) } }, + { "UpperLimits", Conversion::ToVectorOfStrings(m_LimitedEulerTransform->GetUpperLimits()) }, + { "LowerLimits", Conversion::ToVectorOfStrings(m_LimitedEulerTransform->GetLowerLimits()) }, + { "UpperLimitsReached", Conversion::ToVectorOfStrings(m_LimitedEulerTransform->GetUpperLimitsReached()) }, + { "LowerLimitsReached", Conversion::ToVectorOfStrings(m_LimitedEulerTransform->GetLowerLimitsReached()) } + }; + + if (SpaceDimension == 3) + { + parameterMap["ComputeZYX"] = { Conversion::ToString(m_LimitedEulerTransform->GetComputeZYX()) }; + } + return parameterMap; +} // end CreateDerivedTransformParametersMap() + +/** + * ************************* CreateDerivedTransformParametersMap ************************ + */ +template +auto +AdvancedLimitedEulerTransformElastix::AfterRegistration(void) -> void +{ +} // end AfterRegistration() + + +/** + * ************************* InitializeTransform ********************* + */ + +template< class TElastix > +void +AdvancedLimitedEulerTransformElastix< TElastix > +::InitializeTransform( void ) +{ + /** Set all parameters to zero (no rotations, no translation). */ + this->m_LimitedEulerTransform->SetIdentity(); + + /** Read and set UpperLimits parameter */ + ParametersType upperLimits; + upperLimits.SetSize(6); + for (unsigned int i = 0; i < 6; i++) + { + /** Returns zero when parameter was in the parameter file */ + bool found = this->m_Configuration->ReadParameter( + upperLimits[i], "UpperLimits", i, false); + if (!found) + { + for (unsigned int j = 0; j < 6; j++) + { + upperLimits[j] = j < 3 ? 3.14f : 200.0f; + } + xl::xout["warning"] << "WARNING: Upper limits should contain six values! " + << "Using default values instead!" << std::endl; + break; + } + } + this->m_LimitedEulerTransform->SetUpperLimits(upperLimits); + + /** Read and set LowerLimits parameter */ + ParametersType lowerLimits; + lowerLimits.SetSize(6); + for (unsigned int i = 0; i < 6; i++) + { + /** Returns zero when parameter was in the parameter file */ + bool found = this->m_Configuration->ReadParameter( + lowerLimits[i], "LowerLimits", i, false); + if (!found) + { + for (unsigned int j = 0; j < 6; j++) + { + lowerLimits[j] = j < 3 ? -3.14f : -200.0f; + } + xl::xout["warning"] << "WARNING: Lower limits should contain six values! " + << "Using default values instead!" << std::endl; + break; + } + } + this->m_LimitedEulerTransform->SetLowerLimits(lowerLimits); + + + /** Verift that upper limits are larger than lower limits, otherwise throw error */ + for (unsigned int j = 0; j < 6; j++) + { + if (upperLimits[j] < lowerLimits[j]) + { + itkExceptionMacro(<< "ERROR: Entry in the UpperLimits is smaller than " + << " its corresponding entry in the LowerLimits."); + } + } + + /** Read SharpnessOfLimits parameter */ + ScalarType sharpnessOfLimits = 50.0f; + this->m_Configuration->ReadParameter(sharpnessOfLimits, "SharpnessOfLimits", 0); + this->m_LimitedEulerTransform->SetSharpnessOfLimits(sharpnessOfLimits); + + /** Try to read CenterOfRotationIndex from parameter file, + * which is the rotationPoint, expressed in index-values. + */ + IndexType centerOfRotationIndex; + InputPointType centerOfRotationPoint; + bool centerGivenAsIndex = true; + bool centerGivenAsPoint = true; + for( unsigned int i = 0; i < SpaceDimension; i++ ) + { + /** Initialize. */ + centerOfRotationIndex[ i ] = 0; + centerOfRotationPoint[ i ] = 0.0; + + /** Check COR index: Returns zero when parameter was in the parameter file. */ + bool foundI = this->m_Configuration->ReadParameter( + centerOfRotationIndex[ i ], "CenterOfRotation", i, false ); + if( !foundI ) + { + centerGivenAsIndex &= false; + } + + /** Check COR point: Returns zero when parameter was in the parameter file. */ + bool foundP = this->m_Configuration->ReadParameter( + centerOfRotationPoint[ i ], "CenterOfRotationPoint", i, false ); + if( !foundP ) + { + centerGivenAsPoint &= false; + } + + } // end loop over SpaceDimension + + /** Check if CenterOfRotation has index-values within image. */ + bool CORIndexInImage = true; + bool CORPointInImage = true; + if( centerGivenAsIndex ) + { + CORIndexInImage = this->m_Registration->GetAsITKBaseType() + ->GetFixedImage()->GetLargestPossibleRegion().IsInside( + centerOfRotationIndex ); + } + + if( centerGivenAsPoint ) + { + typedef itk::ContinuousIndex< double, SpaceDimension > ContinuousIndexType; + ContinuousIndexType cindex; + CORPointInImage = this->m_Registration->GetAsITKBaseType() + ->GetFixedImage()->TransformPhysicalPointToContinuousIndex( + centerOfRotationPoint, cindex ); + } + + /** Give a warning if necessary. */ + if( !CORIndexInImage && centerGivenAsIndex ) + { + xl::xout[ "warning" ] << "WARNING: Center of Rotation (index) is not " + << "within image boundaries!" << std::endl; + } + + /** Give a warning if necessary. */ + if( !CORPointInImage && centerGivenAsPoint && !centerGivenAsIndex ) + { + xl::xout[ "warning" ] << "WARNING: Center of Rotation (point) is not " + << "within image boundaries!" << std::endl; + } + + /** Check if user wants automatic transform initialization; false by default. + * If an initial transform is given, automatic transform initialization is + * not possible. + */ + bool automaticTransformInitialization = false; + bool tmpBool = false; + this->m_Configuration->ReadParameter( tmpBool, + "AutomaticTransformInitialization", 0 ); + if( tmpBool && this->Superclass1::GetInitialTransform() == 0 ) + { + automaticTransformInitialization = true; + } + + /** + * Run the itkTransformInitializer if: + * - No center of rotation was given, or + * - The user asked for AutomaticTransformInitialization + */ + bool centerGiven = centerGivenAsIndex || centerGivenAsPoint; + if( !centerGiven || automaticTransformInitialization ) + { + + /** Use the TransformInitializer to determine a center of + * of rotation and an initial translation. + */ + TransformInitializerPointer transformInitializer + = TransformInitializerType::New(); + transformInitializer->SetFixedImage( + this->m_Registration->GetAsITKBaseType()->GetFixedImage() ); + transformInitializer->SetMovingImage( + this->m_Registration->GetAsITKBaseType()->GetMovingImage() ); + transformInitializer->SetTransform( this->m_LimitedEulerTransform ); + + /** Select the method of initialization. Default: "GeometricalCenter". */ + transformInitializer->GeometryOn(); + std::string method = "GeometricalCenter"; + this->m_Configuration->ReadParameter( method, + "AutomaticTransformInitializationMethod", 0 ); + if( method == "CenterOfGravity" ) + { + transformInitializer->MomentsOn(); + } + + transformInitializer->InitializeTransform(); + } + + /** Set the translation to zero, if no AutomaticTransformInitialization + * was desired. + */ + if( !automaticTransformInitialization ) + { + OutputVectorType noTranslation; + noTranslation.Fill( 0.0 ); + this->m_LimitedEulerTransform->SetTranslation( noTranslation ); + } + + /** Set the center of rotation if it was entered by the user. */ + if( centerGiven ) + { + if( centerGivenAsIndex ) + { + /** Convert from index-value to physical-point-value. */ + this->m_Registration->GetAsITKBaseType()->GetFixedImage() + ->TransformIndexToPhysicalPoint( + centerOfRotationIndex, centerOfRotationPoint ); + } + this->m_LimitedEulerTransform->SetCenter( centerOfRotationPoint ); + } + + /** Apply the initial transform to the center of rotation, if + * composition is used to combine the initial transform with the + * the current (euler) transform. + */ + if( this->GetUseComposition() + && this->Superclass1::GetInitialTransform() != 0 ) + { + InputPointType transformedCenterOfRotationPoint + = this->Superclass1::GetInitialTransform()->TransformPoint( + this->m_LimitedEulerTransform->GetCenter() ); + this->m_LimitedEulerTransform->SetCenter( transformedCenterOfRotationPoint ); + } + + /** Set the initial parameters in this->m_Registration. */ + this->m_Registration->GetAsITKBaseType() + ->SetInitialTransformParameters( this->GetParameters() ); + + /** Give feedback. */ + // \todo: should perhaps also print fixed parameters + elxout << "Transform parameters are initialized as: " + << this->GetParameters() << std::endl; + +} // end InitializeTransform() + + +/** + * ************************* SetScales ********************* + */ + +template< class TElastix > +void +AdvancedLimitedEulerTransformElastix< TElastix > +::SetScales( void ) +{ + /** Create the new scales. */ + const NumberOfParametersType N = this->GetNumberOfParameters(); + ScalesType newscales( N ); + newscales.Fill( 1.0 ); + + /** Check if automatic scales estimation is desired. */ + bool automaticScalesEstimation = false; + this->m_Configuration->ReadParameter( automaticScalesEstimation, + "AutomaticScalesEstimation", 0 ); + + if( automaticScalesEstimation ) + { + elxout << "Scales are estimated automatically." << std::endl; + this->AutomaticScalesEstimation( newscales ); + } + else + { + + /** Here is an heuristic rule for estimating good values for + * the rotation/translation scales. + * + * 1) Estimate the bounding box of your points (in physical units). + * 2) Take the 3D Diagonal of that bounding box + * 3) Multiply that by 10.0. + * 4) use 1.0 /[ value from (3) ] as the translation scaling value. + * 5) use 1.0 as the rotation scaling value. + * + * With this operation you bring the translation units + * to the range of rotations (e.g. around -1 to 1). + * After that, all your registration parameters are + * in the relaxed range of -1:1. At that point you + * can start setting your optimizer with step lengths + * in the ranges of 0.001 if you are conservative, or + * in the range of 0.1 if you want to live dangerously. + * (0.1 radians is about 5.7 degrees). + * + * This heuristic rule is based on the naive assumption + * that your registration may require translations as + * large as 1/10 of the diagonal of the bounding box. + */ + + const double defaultScalingvalue = 100000.0; + + /** If the Dimension is 3, the first 3 parameters represent rotations. + * If the Dimension is 2, only the first parameter represent a rotation. + */ + unsigned int RotationPart = 3; + if( SpaceDimension == 2 ) { RotationPart = 1; } + + /** this->m_Configuration->ReadParameter() returns 0 if there is a value given + * in the parameter-file, and returns 1 if there is no value given in the + * parameter-file. + * Check which option is used: + * - Nothing given in the parameter-file: rotations are scaled by the default + * value 100000.0 + * - Only one scale given in the parameter-file: rotations are scaled by this + * value. + * - All scales are given in the parameter-file: each parameter is assigned its + * own scale. + */ + + std::size_t count + = this->m_Configuration->CountNumberOfParameterEntries( "Scales" ); + + /** Check which of the above options is used. */ + if( count == 0 ) + { + /** In this case the first option is used. */ + for( unsigned int i = 0; i < RotationPart; i++ ) + { + newscales[ i ] = defaultScalingvalue; + } + } + else if( count == 1 ) + { + /** In this case the second option is used. */ + double scale = defaultScalingvalue; + this->m_Configuration->ReadParameter( scale, "Scales", 0 ); + for( unsigned int i = 0; i < RotationPart; i++ ) + { + newscales[ i ] = scale; + } + } + else if( count == this->GetNumberOfParameters() ) + { + /** In this case the third option is used. */ + for( unsigned int i = 0; i < this->GetNumberOfParameters(); i++ ) + { + this->m_Configuration->ReadParameter( newscales[ i ], "Scales", i ); + } + } + else + { + /** In this case an error is made in the parameter-file. + * An error is thrown, because using erroneous scales in the optimizer + * can give unpredictable results. + */ + itkExceptionMacro( << "ERROR: The Scales-option in the parameter-file" + << " has not been set properly." ); + } + + } // end else: no automaticScalesEstimation + + elxout << "Scales for transform parameters are: " << newscales << std::endl; + + /** Set the scales into the optimizer. */ + this->m_Registration->GetAsITKBaseType()->GetModifiableOptimizer()->SetScales( newscales ); + +} // end SetScales() + + +/** + * ******************** ReadCenterOfRotationPoint ********************* + */ + +template< class TElastix > +bool +AdvancedLimitedEulerTransformElastix< TElastix > +::ReadCenterOfRotationPoint( InputPointType & rotationPoint ) const +{ + /** Try to read CenterOfRotationPoint from the transform parameter + * file, which is the rotationPoint, expressed in world coordinates. + */ + InputPointType centerOfRotationPoint; + bool centerGivenAsPoint = true; + for( unsigned int i = 0; i < SpaceDimension; i++ ) + { + centerOfRotationPoint[ i ] = 0; + + /** Returns zero when parameter was in the parameter file */ + bool found = this->m_Configuration->ReadParameter( + centerOfRotationPoint[ i ], "CenterOfRotationPoint", i, false ); + if( !found ) + { + centerGivenAsPoint &= false; + } + } + + if( !centerGivenAsPoint ) + { + return false; + } + + /** Copy the temporary variable into the output of this function, + * if everything went ok. + */ + rotationPoint = centerOfRotationPoint; + + /** Successfully read centerOfRotation as Point. */ + return true; + +} // end ReadCenterOfRotationPoint() + +} // end namespace elastix + +#endif // end #ifndef __elxAdvancedLimitedEulerTransform_HXX_ diff --git a/Components/Transforms/AdvancedLimitedEulerTransform/itkAdvancedLimitedEuler3DTransform.h b/Components/Transforms/AdvancedLimitedEulerTransform/itkAdvancedLimitedEuler3DTransform.h new file mode 100644 index 000000000..6a4cbd25e --- /dev/null +++ b/Components/Transforms/AdvancedLimitedEulerTransform/itkAdvancedLimitedEuler3DTransform.h @@ -0,0 +1,221 @@ +/*========================================================================= + * + * Copyright UMC Utrecht 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. + * + *=========================================================================*/ +/*========================================================================= + + Program: Insight Segmentation & Registration Toolkit + Module: $RCSfile: itkAdvancedLimitedEuler3DTransform.h,v $ + Language: C++ + Date: $Date: 2008-10-13 15:36:31 $ + Version: $Revision: 1.14 $ + + Copyright (c) Insight Software Consortium. All rights reserved. + See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details. + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notices for more information. + +=========================================================================*/ +#ifndef __itkAdvancedLimitedEuler3DTransform_h +#define __itkAdvancedLimitedEuler3DTransform_h + +#include +#include "itkAdvancedRigid3DTransform.h" + +namespace itk +{ + +/** \class AdvancedLimitedEuler3DTransform + * + * \brief AdvancedLimitedEuler3DTransform of a vector space (e.g. space coordinates) + * + * This transform applies a rotation and translation to the space given 3 euler + * angles and a 3D translation. Rotation is about a user specified center. + * + * The parameters for this transform can be set either using individual Set + * methods or in serialized form using SetParameters() and SetFixedParameters(). + * + * The serialization of the optimizable parameters is an array of 6 elements. + * The first 3 represents three euler angle of rotation respectively about + * the X, Y and Z axis. The last 3 parameters defines the translation in each + * dimension. + * + * The serialization of the fixed parameters is an array of 3 elements defining + * the center of rotation. + * + * \ingroup Transforms + */ +template< class TScalarType = double > +// Data type for scalars (float or double) +class ITK_EXPORT AdvancedLimitedEuler3DTransform : + public AdvancedRigid3DTransform< TScalarType > +{ +public: + + /** Standard class typedefs. */ + typedef AdvancedLimitedEuler3DTransform Self; + typedef AdvancedRigid3DTransform< TScalarType > Superclass; + typedef SmartPointer< Self > Pointer; + typedef SmartPointer< const Self > ConstPointer; + + /** New macro for creation of through a Smart Pointer. */ + itkNewMacro( Self ); + + /** Run-time type information (and related methods). */ + itkTypeMacro( AdvancedLimitedEuler3DTransform, AdvancedRigid3DTransform ); + + /** Dimension of the space. */ + itkStaticConstMacro( SpaceDimension, unsigned int, 3 ); + itkStaticConstMacro( InputSpaceDimension, unsigned int, 3 ); + itkStaticConstMacro( OutputSpaceDimension, unsigned int, 3 ); + itkStaticConstMacro( ParametersDimension, unsigned int, 6 ); + + typedef typename Superclass::ParametersType ParametersType; + typedef typename Superclass::NumberOfParametersType NumberOfParametersType; + typedef typename Superclass::JacobianType JacobianType; + typedef typename Superclass::ScalarType ScalarType; + typedef typename Superclass::InputVectorType InputVectorType; + typedef typename Superclass::OutputVectorType OutputVectorType; + typedef typename Superclass::InputCovariantVectorType InputCovariantVectorType; + typedef typename Superclass::OutputCovariantVectorType OutputCovariantVectorType; + typedef typename Superclass::InputVnlVectorType InputVnlVectorType; + typedef typename Superclass::OutputVnlVectorType OutputVnlVectorType; + typedef typename Superclass::InputPointType InputPointType; + typedef typename Superclass::OutputPointType OutputPointType; + typedef typename Superclass::MatrixType MatrixType; + typedef typename Superclass::InverseMatrixType InverseMatrixType; + typedef typename Superclass::CenterType CenterType; + typedef typename Superclass::TranslationType TranslationType; + typedef typename Superclass::OffsetType OffsetType; + typedef typename Superclass::ScalarType AngleType; + + typedef typename Superclass + ::NonZeroJacobianIndicesType NonZeroJacobianIndicesType; + typedef typename Superclass::SpatialJacobianType SpatialJacobianType; + typedef typename Superclass + ::JacobianOfSpatialJacobianType JacobianOfSpatialJacobianType; + typedef typename Superclass::SpatialHessianType SpatialHessianType; + typedef typename Superclass + ::JacobianOfSpatialHessianType JacobianOfSpatialHessianType; + typedef typename Superclass::InternalMatrixType InternalMatrixType; + + /** Set/Get the transformation from a container of parameters + * This is typically used by optimizers. There are 6 parameters. The first + * three represent the angles to rotate around the coordinate axis, and the + * last three represents the offset. */ + void SetParameters( const ParametersType & parameters ) override; + + const ParametersType & GetParameters( void ) const override; + + /** Set the rotational part of the transform. */ + void SetRotation( ScalarType angleX, ScalarType angleY, ScalarType angleZ ); + + itkGetConstMacro( AngleX, ScalarType ); + itkGetConstMacro( AngleY, ScalarType ); + itkGetConstMacro( AngleZ, ScalarType ); + + /** Compute the Jacobian of the transformation. */ + void GetJacobian( + const InputPointType &, + JacobianType &, + NonZeroJacobianIndicesType & ) const override; + + /** Set/Get the order of the computation. Default ZXY */ + itkSetMacro( ComputeZYX, bool ); + itkGetConstMacro( ComputeZYX, bool ); + + void SetIdentity( void ) override; + + /** Set/Get the whether the scales estimation is in progress */ + itkSetMacro(ScalesEstimation, bool); + + /** Set/Get the sharpness aka smoothness of parameter limits */ + itkSetMacro(SharpnessOfLimits, double); + itkGetConstMacro(SharpnessOfLimits, double); + + /** Update sharpness of limits parameter for each axis based on the inverval width */ + virtual void UpdateSharpnessOfLimitsVector(); + + /** Setters/getters for the upper/lower limits */ + virtual void SetUpperLimits(const ParametersType & upperLimits); + virtual const ParametersType& GetUpperLimits(); + + virtual void SetLowerLimits(const ParametersType & lowerLimits); + virtual const ParametersType& GetLowerLimits(); + + /** Upper/Lower limits reached indicator. Values close to 1 indicate that limit was reached. */ + virtual const ParametersType& GetUpperLimitsReached(); + virtual const ParametersType& GetLowerLimitsReached(); + +protected: + + AdvancedLimitedEuler3DTransform(); + AdvancedLimitedEuler3DTransform( const MatrixType & matrix, + const OutputPointType & offset ); + AdvancedLimitedEuler3DTransform( unsigned int paramsSpaceDims ); + + ~AdvancedLimitedEuler3DTransform() override {} + + void PrintSelf( std::ostream & os, Indent indent ) const override; + + /** Set values of angles directly without recomputing other parameters. */ + void SetVarRotation( ScalarType angleX, ScalarType angleY, ScalarType angleZ ); + + /** Compute the components of the rotation matrix in the superclass. */ + void ComputeMatrix( void ) override; + + void ComputeMatrixParameters( void ) override; + + /** Update the m_JacobianOfSpatialJacobian. */ + virtual void PrecomputeJacobianOfSpatialJacobian( void ); + + /** Initialize limits parameters to default values */ + void InitLimitParameters(); + + /** Apply Softplus limiting function to input value and based on range parameters */ + ScalarType SoftplusLimit(const ScalarType & inputValue, const ScalarType & limitValue, const ScalarType & sharpnessOfLimits, bool setHighLimit) const; + + /** Compute derivative of Softplus function at input value and based on range parameters */ + ScalarType DerivativeSoftplusLimit(const ScalarType & inputValue, const ScalarType & limitValue, const ScalarType & sharpnessOfLimits, bool setHighLimit) const; + +private: + + AdvancedLimitedEuler3DTransform( const Self & ); //purposely not implemented + void operator=( const Self & ); //purposely not implemented + + ScalarType m_AngleX; + ScalarType m_AngleY; + ScalarType m_AngleZ; + bool m_ComputeZYX; + + ScalarType m_SharpnessOfLimits; + ParametersType m_SharpnessOfLimitsVector; + ParametersType m_UpperLimits; + ParametersType m_LowerLimits; + ParametersType m_UpperLimitsReached; + ParametersType m_LowerLimitsReached; + bool m_ScalesEstimation; +}; //class AdvancedLimitedEuler3DTransform + +} // namespace itk + +#ifndef ITK_MANUAL_INSTANTIATION +#include "itkAdvancedLimitedEuler3DTransform.hxx" +#endif + +#endif /* __itkAdvancedLimitedEuler3DTransform_h */ diff --git a/Components/Transforms/AdvancedLimitedEulerTransform/itkAdvancedLimitedEuler3DTransform.hxx b/Components/Transforms/AdvancedLimitedEulerTransform/itkAdvancedLimitedEuler3DTransform.hxx new file mode 100644 index 000000000..90b525645 --- /dev/null +++ b/Components/Transforms/AdvancedLimitedEulerTransform/itkAdvancedLimitedEuler3DTransform.hxx @@ -0,0 +1,622 @@ +/*========================================================================= + * + * Copyright UMC Utrecht 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. + * + *=========================================================================*/ +/*========================================================================= + + Program: Insight Segmentation & Registration Toolkit + Module: $RCSfile: itkAdvancedLimitedEuler3DTransform.txx,v $ + Language: C++ + Date: $Date: 2008-10-13 15:36:31 $ + Version: $Revision: 1.24 $ + + Copyright (c) Insight Software Consortium. All rights reserved. + See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details. + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notices for more information. + +=========================================================================*/ +#ifndef __itkAdvancedLimitedEuler3DTransform_hxx +#define __itkAdvancedLimitedEuler3DTransform_hxx + +#include "itkAdvancedLimitedEuler3DTransform.h" + +namespace itk +{ + // Constructor with default arguments + template< class TScalarType > + AdvancedLimitedEuler3DTransform< TScalarType > + ::AdvancedLimitedEuler3DTransform() : + Superclass(ParametersDimension) + { + m_ComputeZYX = false; + m_AngleX = m_AngleY = m_AngleZ = NumericTraits< ScalarType >::Zero; + this->PrecomputeJacobianOfSpatialJacobian(); + + this->InitLimitParameters(); + } + + + // Constructor with default arguments + template< class TScalarType > + AdvancedLimitedEuler3DTransform< TScalarType > + ::AdvancedLimitedEuler3DTransform(const MatrixType& matrix, + const OutputPointType& offset) + { + m_ComputeZYX = false; + this->SetMatrix(matrix); + + OffsetType off; + off[0] = offset[0]; + off[1] = offset[1]; + off[2] = offset[2]; + this->SetOffset(off); + + this->PrecomputeJacobianOfSpatialJacobian(); + + this->InitLimitParameters(); + } + + + // Constructor with arguments + template< class TScalarType > + AdvancedLimitedEuler3DTransform< TScalarType > + ::AdvancedLimitedEuler3DTransform(unsigned int parametersDimension) : + Superclass(parametersDimension) + { + m_ComputeZYX = false; + m_AngleX = m_AngleY = m_AngleZ = NumericTraits< ScalarType >::Zero; + this->PrecomputeJacobianOfSpatialJacobian(); + + this->InitLimitParameters(); + } + + + // Set Angles + template< class TScalarType > + void + AdvancedLimitedEuler3DTransform< TScalarType > + ::SetVarRotation(ScalarType angleX, ScalarType angleY, ScalarType angleZ) + { + this->m_AngleX = angleX; + this->m_AngleY = angleY; + this->m_AngleZ = angleZ; + } + + + // Set Parameters + template< class TScalarType > + void + AdvancedLimitedEuler3DTransform< TScalarType > + ::SetParameters(const ParametersType& parameters) + { + itkDebugMacro(<< "Setting parameters " << parameters); + + // Set angles with parameters + m_AngleX = parameters[0]; + m_AngleY = parameters[1]; + m_AngleZ = parameters[2]; + this->ComputeMatrix(); + + // Transfer the translation part + OutputVectorType newTranslation; + newTranslation[0] = parameters[3]; + newTranslation[1] = parameters[4]; + newTranslation[2] = parameters[5]; + + this->SetVarTranslation(newTranslation); + this->ComputeOffset(); + + // Modified is always called since we just have a pointer to the + // parameters and cannot know if the parameters have changed. + this->Modified(); + + itkDebugMacro(<< "After setting parameters "); + } + + + // Get Parameters + template< class TScalarType > + const typename AdvancedLimitedEuler3DTransform< TScalarType >::ParametersType + & AdvancedLimitedEuler3DTransform< TScalarType > + ::GetParameters(void) const + { + this->m_Parameters[0] = m_AngleX; + this->m_Parameters[1] = m_AngleY; + this->m_Parameters[2] = m_AngleZ; + this->m_Parameters[3] = this->GetTranslation()[0]; + this->m_Parameters[4] = this->GetTranslation()[1]; + this->m_Parameters[5] = this->GetTranslation()[2]; + return this->m_Parameters; + } + + // Set Rotational Part + template< class TScalarType > + void + AdvancedLimitedEuler3DTransform< TScalarType > + ::SetRotation(ScalarType angleX, ScalarType angleY, ScalarType angleZ) + { + m_AngleX = angleX; + m_AngleY = angleY; + m_AngleZ = angleZ; + this->ComputeMatrix(); + this->ComputeOffset(); + } + + // Compose + template< class TScalarType > + void + AdvancedLimitedEuler3DTransform< TScalarType > + ::SetIdentity(void) + { + Superclass::SetIdentity(); + m_AngleX = 0; + m_AngleY = 0; + m_AngleZ = 0; + this->PrecomputeJacobianOfSpatialJacobian(); + } + + // Compute angles from the rotation matrix + template< class TScalarType > + void + AdvancedLimitedEuler3DTransform< TScalarType > + ::ComputeMatrixParameters(void) + { + if (m_ComputeZYX) + { + m_AngleY = -asin(this->GetMatrix()[2][0]); + double C = std::cos(m_AngleY); + if (std::fabs(C) > 0.00005) + { + double x = this->GetMatrix()[2][2] / C; + double y = this->GetMatrix()[2][1] / C; + m_AngleX = std::atan2(y, x); + x = this->GetMatrix()[0][0] / C; + y = this->GetMatrix()[1][0] / C; + m_AngleZ = std::atan2(y, x); + } + else + { + m_AngleX = 0; + double x = this->GetMatrix()[1][1]; + double y = -this->GetMatrix()[0][1]; + m_AngleZ = std::atan2(y, x); + } + } + else + { + m_AngleX = std::asin(this->GetMatrix()[2][1]); + double A = std::cos(m_AngleX); + if (std::fabs(A) > 0.00005) + { + double x = this->GetMatrix()[2][2] / A; + double y = -this->GetMatrix()[2][0] / A; + m_AngleY = std::atan2(y, x); + + x = this->GetMatrix()[1][1] / A; + y = -this->GetMatrix()[0][1] / A; + m_AngleZ = std::atan2(y, x); + } + else + { + m_AngleZ = 0; + double x = this->GetMatrix()[0][0]; + double y = this->GetMatrix()[1][0]; + m_AngleY = std::atan2(y, x); + } + } + this->ComputeMatrix(); + } + + + // Compute the matrix + template< class TScalarType > + void + AdvancedLimitedEuler3DTransform< TScalarType > + ::ComputeMatrix(void) + { + // need to check if angles are in the right order + const double cx = std::cos(m_AngleX); + const double sx = std::sin(m_AngleX); + const double cy = std::cos(m_AngleY); + const double sy = std::sin(m_AngleY); + const double cz = std::cos(m_AngleZ); + const double sz = std::sin(m_AngleZ); + + Matrix< TScalarType, 3, 3 > RotationX; + RotationX[0][0] = 1; RotationX[0][1] = 0; RotationX[0][2] = 0; + RotationX[1][0] = 0; RotationX[1][1] = cx; RotationX[1][2] = -sx; + RotationX[2][0] = 0; RotationX[2][1] = sx; RotationX[2][2] = cx; + + Matrix< TScalarType, 3, 3 > RotationY; + RotationY[0][0] = cy; RotationY[0][1] = 0; RotationY[0][2] = sy; + RotationY[1][0] = 0; RotationY[1][1] = 1; RotationY[1][2] = 0; + RotationY[2][0] = -sy; RotationY[2][1] = 0; RotationY[2][2] = cy; + + Matrix< TScalarType, 3, 3 > RotationZ; + RotationZ[0][0] = cz; RotationZ[0][1] = -sz; RotationZ[0][2] = 0; + RotationZ[1][0] = sz; RotationZ[1][1] = cz; RotationZ[1][2] = 0; + RotationZ[2][0] = 0; RotationZ[2][1] = 0; RotationZ[2][2] = 1; + + // Aply the rotation first around Y then X then Z + if (m_ComputeZYX) + { + this->SetVarMatrix(RotationZ * RotationY * RotationX); + } + else + { + // Like VTK transformation order + this->SetVarMatrix(RotationZ * RotationX * RotationY); + } + + this->PrecomputeJacobianOfSpatialJacobian(); + } + + + // Get Jacobian + template< class TScalarType > + void + AdvancedLimitedEuler3DTransform< TScalarType >::GetJacobian(const InputPointType& p, + JacobianType& j, + NonZeroJacobianIndicesType& nzji) const + { + // Compute limit derivatives based on current parameter values + ParametersType parameterValues, dUdL; + parameterValues.SetSize(ParametersDimension); + parameterValues[0] = m_AngleX; + parameterValues[1] = m_AngleY; + parameterValues[2] = m_AngleZ; + parameterValues[3] = this->GetTranslation()[0]; + parameterValues[4] = this->GetTranslation()[1]; + parameterValues[5] = this->GetTranslation()[2]; + + dUdL.SetSize(ParametersDimension); + for (unsigned int i = 0; i < ParametersDimension; i++) + { + dUdL[i] = m_ScalesEstimation ? 1.0f : DerivativeSoftplusLimit(parameterValues[i], m_UpperLimits[i], m_SharpnessOfLimitsVector[i], true) + * DerivativeSoftplusLimit(parameterValues[i], m_LowerLimits[i], m_SharpnessOfLimitsVector[i], false); + } + + // Initialize the Jacobian. Resizing is only performed when needed. + // Filling with zeros is needed because the lower loops only visit + // the nonzero positions. + j.SetSize(OutputSpaceDimension, ParametersDimension); + j.Fill(0.0); + + // Compute dR/dmu * (p-c) + const InputVectorType pp = p - this->GetCenter(); + const JacobianOfSpatialJacobianType& jsj = this->m_JacobianOfSpatialJacobian; + for (unsigned int dim = 0; dim < SpaceDimension; ++dim) + { + const InputVectorType column = jsj[dim] * pp; + for (unsigned int i = 0; i < SpaceDimension; ++i) + { + j(i, dim) = column[i]; + } + } + + // compute derivatives for the translation part + const unsigned int blockOffset = 3; + for (unsigned int dim = 0; dim < SpaceDimension; ++dim) + { + j[dim][0] *= dUdL[0]; + j[dim][1] *= dUdL[1]; + j[dim][2] *= dUdL[2]; + j[dim][blockOffset + dim] = dUdL[SpaceDimension + dim]; + } + + // Copy the constant nonZeroJacobianIndices + nzji = this->m_NonZeroJacobianIndices; + } + + + // Precompute Jacobian of Spatial Jacobian + template< class TScalarType > + void + AdvancedLimitedEuler3DTransform< TScalarType > + ::PrecomputeJacobianOfSpatialJacobian(void) + { + if (ParametersDimension < 6) + { + // Some subclass has a different number of parameters + return; + } + + // The Jacobian of spatial Jacobian is constant over inputspace, so is precomputed + JacobianOfSpatialJacobianType& jsj = this->m_JacobianOfSpatialJacobian; + jsj.resize(ParametersDimension); + const double cx = std::cos(m_AngleX); + const double sx = std::sin(m_AngleX); + const double cy = std::cos(m_AngleY); + const double sy = std::sin(m_AngleY); + const double cz = std::cos(m_AngleZ); + const double sz = std::sin(m_AngleZ); + + // derivatives: + const double cxd = -std::sin(m_AngleX); + const double sxd = std::cos(m_AngleX); + const double cyd = -std::sin(m_AngleY); + const double syd = std::cos(m_AngleY); + const double czd = -std::sin(m_AngleZ); + const double szd = std::cos(m_AngleZ); + + // rotation matrices + Matrix< TScalarType, 3, 3 > RotationX; + RotationX[0][0] = 1; RotationX[0][1] = 0; RotationX[0][2] = 0; + RotationX[1][0] = 0; RotationX[1][1] = cx; RotationX[1][2] = -sx; + RotationX[2][0] = 0; RotationX[2][1] = sx; RotationX[2][2] = cx; + + Matrix< TScalarType, 3, 3 > RotationY; + RotationY[0][0] = cy; RotationY[0][1] = 0; RotationY[0][2] = sy; + RotationY[1][0] = 0; RotationY[1][1] = 1; RotationY[1][2] = 0; + RotationY[2][0] = -sy; RotationY[2][1] = 0; RotationY[2][2] = cy; + + Matrix< TScalarType, 3, 3 > RotationZ; + RotationZ[0][0] = cz; RotationZ[0][1] = -sz; RotationZ[0][2] = 0; + RotationZ[1][0] = sz; RotationZ[1][1] = cz; RotationZ[1][2] = 0; + RotationZ[2][0] = 0; RotationZ[2][1] = 0; RotationZ[2][2] = 1; + + // derivative matrices + Matrix< TScalarType, 3, 3 > RotationXd; + RotationXd[0][0] = 0; RotationXd[0][1] = 0; RotationXd[0][2] = 0; + RotationXd[1][0] = 0; RotationXd[1][1] = cxd; RotationXd[1][2] = -sxd; + RotationXd[2][0] = 0; RotationXd[2][1] = sxd; RotationXd[2][2] = cxd; + + Matrix< TScalarType, 3, 3 > RotationYd; + RotationYd[0][0] = cyd; RotationYd[0][1] = 0; RotationYd[0][2] = syd; + RotationYd[1][0] = 0; RotationYd[1][1] = 0; RotationYd[1][2] = 0; + RotationYd[2][0] = -syd; RotationYd[2][1] = 0; RotationYd[2][2] = cyd; + + Matrix< TScalarType, 3, 3 > RotationZd; + RotationZd[0][0] = czd; RotationZd[0][1] = -szd; RotationZd[0][2] = 0; + RotationZd[1][0] = szd; RotationZd[1][1] = czd; RotationZd[1][2] = 0; + RotationZd[2][0] = 0; RotationZd[2][1] = 0; RotationZd[2][2] = 0; + + // Aply the rotation first around Y then X then Z + if (m_ComputeZYX) + { + jsj[0] = RotationZ * RotationY * RotationXd; + jsj[1] = RotationZ * RotationYd * RotationX; + jsj[2] = RotationZd * RotationY * RotationX; + + } + else + { + // Like VTK transformation order + jsj[0] = RotationZ * RotationXd * RotationY; + jsj[1] = RotationZ * RotationX * RotationYd; + jsj[2] = RotationZd * RotationX * RotationY; + } + + // Translation parameters: + for (unsigned int par = 3; par < ParametersDimension; ++par) + { + jsj[par].Fill(0.0); + } + } + + + // Print self + template< class TScalarType > + void + AdvancedLimitedEuler3DTransform< TScalarType >::PrintSelf(std::ostream& os, Indent indent) const + { + Superclass::PrintSelf(os, indent); + + os << indent << "Euler's angles: AngleX=" << m_AngleX + << " AngleY=" << m_AngleY + << " AngleZ=" << m_AngleZ + << std::endl; + os << indent << "m_ComputeZYX = " << m_ComputeZYX << std::endl; + } + + + template + void + AdvancedLimitedEuler3DTransform::InitLimitParameters() + { + m_ScalesEstimation = false; + m_SharpnessOfLimits = 50.0f; + + m_UpperLimits.SetSize(ParametersDimension); + m_LowerLimits.SetSize(ParametersDimension); + for (unsigned int i = 0; i < ParametersDimension / 2; i++) + { + m_UpperLimits[i] = 3.14f; + m_LowerLimits[i] = -3.14f; + m_UpperLimits[i + ParametersDimension / 2] = 200.0f; + m_LowerLimits[i + ParametersDimension / 2] = -200.0f; + } + this->UpdateSharpnessOfLimitsVector(); + } + + + template + void + AdvancedLimitedEuler3DTransform::UpdateSharpnessOfLimitsVector() + { + m_SharpnessOfLimitsVector.SetSize(ParametersDimension); + for (unsigned int i = 0; i < ParametersDimension / 2; i++) + { + m_SharpnessOfLimitsVector[i] = m_SharpnessOfLimits / ((m_UpperLimits[i] - m_LowerLimits[i]) + 1e-3); + m_SharpnessOfLimitsVector[i + ParametersDimension / 2] = + m_SharpnessOfLimits / ((m_UpperLimits[i + ParametersDimension / 2] - m_LowerLimits[i + ParametersDimension / 2]) + 1e-3); + } + } + + template + typename AdvancedLimitedEuler3DTransform::ScalarType + AdvancedLimitedEuler3DTransform::SoftplusLimit(const ScalarType& inputValue, const ScalarType& limitValue, const ScalarType& sharpnessOfLimits, bool setHighLimit) const + { + ScalarType axisDirection = setHighLimit ? -1.0f : 1.0f; + ScalarType nom = std::exp(sharpnessOfLimits * axisDirection * (inputValue - limitValue)); + if (std::isinf(nom)) + { + if (setHighLimit) + { + if (inputValue > limitValue) + { + return limitValue; + } + else + { + return inputValue; + } + } + else + { + if (inputValue < limitValue) + { + return limitValue; + } + else + { + return inputValue; + } + } + } + else + { + return axisDirection * (std::log(1.0f + nom) / sharpnessOfLimits + axisDirection * limitValue); + } + + } + + template + typename AdvancedLimitedEuler3DTransform::ScalarType + AdvancedLimitedEuler3DTransform::DerivativeSoftplusLimit(const ScalarType& inputValue, const ScalarType& limitValue, const ScalarType& sharpnessOfLimits, bool setHighLimit) const + { + ScalarType axisDirection = setHighLimit ? -1.0f : 1.0f; + ScalarType nom = std::exp(sharpnessOfLimits * axisDirection * (inputValue - limitValue)); + if (std::isinf(nom)) + { + if (setHighLimit) + { + if (inputValue > limitValue) + { + return 0.0f; + } + else + { + return 1.0f; + } + } + else + { + if (inputValue < limitValue) + { + return 0.0f; + } + else + { + return 1.0f; + } + } + } + else + { + return (nom / (1.0f + nom)); + } + + } + + template + void + AdvancedLimitedEuler3DTransform::SetUpperLimits(const ParametersType& upperLimits) + { + for (unsigned int i = 0; i < ParametersDimension; i++) + { + m_UpperLimits[i] = upperLimits[i]; + } + } + + template + void + AdvancedLimitedEuler3DTransform::SetLowerLimits(const ParametersType& lowerLimits) + { + for (unsigned int i = 0; i < ParametersDimension; i++) + { + m_LowerLimits[i] = lowerLimits[i]; + } + } + + template + const typename AdvancedLimitedEuler3DTransform::ParametersType& + AdvancedLimitedEuler3DTransform::GetUpperLimits() + { + return this->m_UpperLimits; + } + + template + const typename AdvancedLimitedEuler3DTransform::ParametersType& + AdvancedLimitedEuler3DTransform::GetLowerLimits() + { + return this->m_LowerLimits; + } + + template + const typename AdvancedLimitedEuler3DTransform::ParametersType& + AdvancedLimitedEuler3DTransform::GetUpperLimitsReached() + { + // Compute limit derivatives based on current parameter values + ParametersType parameterValues; + parameterValues.SetSize(ParametersDimension); + parameterValues[0] = m_AngleX; + parameterValues[1] = m_AngleY; + parameterValues[2] = m_AngleZ; + parameterValues[3] = this->GetTranslation()[0]; + parameterValues[4] = this->GetTranslation()[1]; + parameterValues[5] = this->GetTranslation()[2]; + + m_UpperLimitsReached.SetSize(ParametersDimension); + for (unsigned int i = 0; i < ParametersDimension; i++) + { + m_UpperLimitsReached[i] = 1.0f - DerivativeSoftplusLimit(parameterValues[i], m_UpperLimits[i], m_SharpnessOfLimitsVector[i], true); + } + + return m_UpperLimitsReached; + } + + template + const typename AdvancedLimitedEuler3DTransform::ParametersType& + AdvancedLimitedEuler3DTransform::GetLowerLimitsReached() + { + // Compute limit derivatives based on current parameter values + ParametersType parameterValues; + parameterValues.SetSize(ParametersDimension); + parameterValues[0] = m_AngleX; + parameterValues[1] = m_AngleY; + parameterValues[2] = m_AngleZ; + parameterValues[3] = this->GetTranslation()[0]; + parameterValues[4] = this->GetTranslation()[1]; + parameterValues[5] = this->GetTranslation()[2]; + + m_LowerLimitsReached.SetSize(ParametersDimension); + for (unsigned int i = 0; i < ParametersDimension; i++) + { + m_LowerLimitsReached[i] = 1.0f - DerivativeSoftplusLimit(parameterValues[i], m_LowerLimits[i], m_SharpnessOfLimitsVector[i], false); + } + + return m_LowerLimitsReached; + } + + +} // namespace + + +#endif \ No newline at end of file diff --git a/Components/Transforms/AdvancedLimitedEulerTransform/itkLimitedEulerTransform.h b/Components/Transforms/AdvancedLimitedEulerTransform/itkLimitedEulerTransform.h new file mode 100644 index 000000000..7a31ec308 --- /dev/null +++ b/Components/Transforms/AdvancedLimitedEulerTransform/itkLimitedEulerTransform.h @@ -0,0 +1,486 @@ +/*========================================================================= + * + * Copyright UMC Utrecht 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 __itkLimitedEulerTransform_H__ +#define __itkLimitedEulerTransform_H__ + +#include "itkAdvancedRigid2DTransform.h" +#include "itkAdvancedLimitedEuler3DTransform.h" +#include "itkAdvancedMatrixOffsetTransformBase.h" + +namespace itk +{ + + /** + * \class LimitedEulerGroup + * \brief This class only contains a dummy class. + * + */ + + template< unsigned int Dimension > + class LimitedEulerGroup + { + public: + + template< class TScalarType > + class Dummy + { + public: + + /** Typedef's.*/ + typedef AdvancedMatrixOffsetTransformBase< TScalarType, Dimension, Dimension > LimitedEulerTransform_tmp; + + }; + + }; + + /** + * \class LimitedEulerGroup<2> + * \brief This class only contains a dummy class for the 2D case. + * + */ + + template< > + class LimitedEulerGroup< 2 > + { + public: + + template< class TScalarType > + class Dummy + { + public: + + /** Typedef's.*/ + typedef AdvancedRigid2DTransform< TScalarType > LimitedEulerTransform_tmp; + + }; + + }; + + /** + * \class LimitedEulerGroup<3> + * \brief This class only contains a dummy class for the 3D case. + * + */ + + template< > + class LimitedEulerGroup< 3 > + { + public: + + template< class TScalarType > + class Dummy + { + public: + + /** Typedef's.*/ + typedef AdvancedLimitedEuler3DTransform< TScalarType > LimitedEulerTransform_tmp; + }; + + }; + + /** + * \class LimitedEulerGroupTemplate + * \brief This class templates the EulerGroup over its dimension. + * + */ + + template< class TScalarType, unsigned int Dimension > + class LimitedEulerGroupTemplate + { + public: + + typedef LimitedEulerGroupTemplate Self; + typedef TScalarType ScalarType; + itkStaticConstMacro(SpaceDimension, unsigned int, Dimension); + + // This declaration of 'LimitedEuler' does not work with the GCC compiler + // typedef LimitedEulerGroup< itkGetStaticConstMacro( SpaceDimension ) > LimitedEuler; + // The following trick works though: + template< unsigned int D > + class LimitedEulerGroupWrap + { + public: + + typedef LimitedEulerGroup< D > Euler; + }; + + typedef LimitedEulerGroupWrap< Dimension > LimitedEulerGroupWrapInstance; + typedef typename LimitedEulerGroupWrapInstance::Euler LimitedEuler; + + typedef typename LimitedEuler::template Dummy< ScalarType > LimitedEulerDummy; + typedef typename LimitedEulerDummy::LimitedEulerTransform_tmp LimitedEulerTransform_tmp; + + }; + + /** + * \class LimitedEulerTransform + * \brief This class combines the Euler2DTransform with the Euler3DTransform. + * + * This transform is a rigid body transformation. + * + * \ingroup Transforms + */ + + template< class TScalarType, unsigned int Dimension > + class LimitedEulerTransform : + public LimitedEulerGroupTemplate< + TScalarType, Dimension >::LimitedEulerTransform_tmp + { + public: + + /** Standard ITK-stuff. */ + typedef LimitedEulerTransform Self; + typedef typename LimitedEulerGroupTemplate< + TScalarType, Dimension > + ::LimitedEulerTransform_tmp Superclass; + typedef SmartPointer< Self > Pointer; + typedef SmartPointer< const Self > ConstPointer; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Run-time type information (and related methods). */ + itkTypeMacro(LimitedEulerTransform, LimitedEulerGroupTemplate); + + /** Dimension of the domain space. */ + itkStaticConstMacro(SpaceDimension, unsigned int, Dimension); + + /** Typedefs inherited from the superclass. */ + + /** These are both in Rigid2D and Euler3D. */ + typedef typename Superclass::ScalarType ScalarType; + typedef typename Superclass::ParametersType ParametersType; + typedef typename Superclass::NumberOfParametersType NumberOfParametersType; + typedef typename Superclass::JacobianType JacobianType; + typedef typename Superclass::OffsetType OffsetType; + typedef typename Superclass::InputPointType InputPointType; + typedef typename Superclass::OutputPointType OutputPointType; + typedef typename Superclass::InputVectorType InputVectorType; + typedef typename Superclass::OutputVectorType OutputVectorType; + typedef typename Superclass::InputCovariantVectorType InputCovariantVectorType; + typedef typename Superclass::OutputCovariantVectorType OutputCovariantVectorType; + typedef typename Superclass::InputVnlVectorType InputVnlVectorType; + typedef typename Superclass::OutputVnlVectorType OutputVnlVectorType; + + typedef typename Superclass + ::NonZeroJacobianIndicesType NonZeroJacobianIndicesType; + typedef typename Superclass::SpatialJacobianType SpatialJacobianType; + typedef typename Superclass + ::JacobianOfSpatialJacobianType JacobianOfSpatialJacobianType; + typedef typename Superclass::SpatialHessianType SpatialHessianType; + typedef typename Superclass + ::JacobianOfSpatialHessianType JacobianOfSpatialHessianType; + typedef typename Superclass::InternalMatrixType InternalMatrixType; + + /** Make sure SetComputeZYX() is available, also in 2D, + * in which case, its just a dummy function. + */ + void SetComputeZYX(const bool) // No override. + { + static_assert(SpaceDimension != 3, "This is not the specialization is 3D!"); + } + + /** Make sure GetComputeZYX() is available, also in 2D, + * in which case, it just returns false. + */ + bool GetComputeZYX(void) const // No override. + { + static_assert(SpaceDimension != 3, "This is not the specialization is 3D!"); + return false; + } + + /** Make sure setters/getters are available, also in 2D, + * in which case, it just issues an assert error. + */ + void SetScalesEstimation(const bool arg) // No override. + { + static_assert(SpaceDimension != 3, "Not implemented in 2D!"); + } + + void SetSharpnessOfLimits(const ScalarType sharpnessOfLimits) // No override. + { + static_assert(SpaceDimension != 3, "Not implemented in 2D!"); + } + + const ScalarType& GetSharpnessOfLimits() // No override. + { + static_assert(SpaceDimension != 3, "Not implemented in 2D!"); + return ScalarType(); + } + + void SetUpperLimits(const ParametersType & upperLimits) // No override. + { + static_assert(SpaceDimension != 3, "Not implemented in 2D!"); + } + + const ParametersType& GetUpperLimits() + { + static_assert(SpaceDimension != 3, "Not implemented in 2D!"); + return ParametersType(); + } + + void SetLowerLimits(const ParametersType & lowerLimits) + { + static_assert(SpaceDimension != 3, "Not implemented in 2D!"); + } + + const ParametersType& GetLowerLimits() + { + static_assert(SpaceDimension != 3, "Not implemented in 2D!"); + return ParametersType(); + } + + const ParametersType& GetUpperLimitsReached() + { + static_assert(SpaceDimension != 3, "Not implemented in 2D!"); + return ParametersType(); + } + + const ParametersType& GetLowerLimitsReached() + { + static_assert(SpaceDimension != 3, "Not implemented in 2D!"); + return ParametersType(); + } + + protected: + + LimitedEulerTransform() {} + ~LimitedEulerTransform() override {} + + private: + + LimitedEulerTransform(const Self&); // purposely not implemented + void operator=(const Self&); // purposely not implemented + + }; + + template< class TScalarType > + class LimitedEulerTransform : + public LimitedEulerGroupTemplate< + TScalarType, 3 >::LimitedEulerTransform_tmp + { + public: + + /** Standard ITK-stuff. */ + typedef LimitedEulerTransform Self; + typedef typename LimitedEulerGroupTemplate< + TScalarType, 3 > + ::LimitedEulerTransform_tmp Superclass; + typedef SmartPointer< Self > Pointer; + typedef SmartPointer< const Self > ConstPointer; + + typedef typename Superclass::ScalarType ScalarType; + typedef typename Superclass::ParametersType ParametersType; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Run-time type information (and related methods). */ + itkTypeMacro(LimitedEulerTransform, LimitedEulerGroupTemplate); + + /** Dimension of the domain space. */ + itkStaticConstMacro(SpaceDimension, unsigned int, 3); + + + /** Make sure SetComputeZYX() is available, also in 2D, + * in which case, its just a dummy function. + * \note This member function is only an `override` in 3D. + */ + void SetComputeZYX(const bool arg) override + { + static_assert(SpaceDimension == 3, "This specialization is for 3D only!"); + + typedef AdvancedLimitedEuler3DTransform< TScalarType > Euler3DTransformType; + typename Euler3DTransformType::Pointer transform + = dynamic_cast(this); + if (transform) + { + transform->Euler3DTransformType::SetComputeZYX(arg); + } + } + + + /** Make sure GetComputeZYX() is available, also in 2D, + * in which case, it just returns false. + * \note This member function is only an `override` in 3D. + */ + bool GetComputeZYX(void) const override + { + static_assert(SpaceDimension == 3, "This specialization is for 3D only!"); + + typedef AdvancedLimitedEuler3DTransform< TScalarType > Euler3DTransformType; + typename Euler3DTransformType::ConstPointer transform + = dynamic_cast(this); + + if (transform) + { + return transform->Euler3DTransformType::GetComputeZYX(); + } + return false; + } + + + /** Make sure setters/getters are available, also in 3D, + * in which case, it sets/returns values accordingly. + */ + void SetScalesEstimation(const bool arg) override + { + static_assert(SpaceDimension == 3, "This specialization is for 3D only!"); + + typedef AdvancedLimitedEuler3DTransform< TScalarType > Euler3DTransformType; + typename Euler3DTransformType::Pointer transform + = dynamic_cast(this); + + if (transform) + { + return transform->Euler3DTransformType::SetScalesEstimation(arg); + } + } + + void SetSharpnessOfLimits(const ScalarType sharpnessOfLimits) + { + static_assert(SpaceDimension == 3, "This specialization is for 3D only!"); + + typedef AdvancedLimitedEuler3DTransform< TScalarType > Euler3DTransformType; + typename Euler3DTransformType::Pointer transform + = dynamic_cast< Euler3DTransformType * >( this ); + + if( transform ) + { + transform->Euler3DTransformType::SetSharpnessOfLimits(sharpnessOfLimits); + transform->Euler3DTransformType::UpdateSharpnessOfLimitsVector(); + } + } + + const ScalarType& GetSharpnessOfLimits() // const override + { + static_assert(SpaceDimension == 3, "This specialization is for 3D only!"); + + typedef AdvancedLimitedEuler3DTransform< TScalarType > Euler3DTransformType; + typename Euler3DTransformType::Pointer transform + = dynamic_cast< Euler3DTransformType * >( this ); + + if( transform ) + { + return transform->Euler3DTransformType::GetSharpnessOfLimits(); + } + return ScalarType(); + } + + void SetUpperLimits(const ParametersType & upperLimits) override + { + static_assert(SpaceDimension == 3, "This specialization is for 3D only!"); + + typedef AdvancedLimitedEuler3DTransform< TScalarType > Euler3DTransformType; + typename Euler3DTransformType::Pointer transform + = dynamic_cast< Euler3DTransformType * >( this ); + + if( transform ) + { + transform->Euler3DTransformType::SetUpperLimits( upperLimits ); + } + } + + const ParametersType& GetUpperLimits() //const override + { + static_assert(SpaceDimension == 3, "This specialization is for 3D only!"); + + typedef AdvancedLimitedEuler3DTransform< TScalarType > Euler3DTransformType; + typename Euler3DTransformType::Pointer transform + = dynamic_cast< Euler3DTransformType * >( this ); + + if( transform ) + { + return transform->Euler3DTransformType::GetUpperLimits(); + } + return ParametersType(); + } + + void SetLowerLimits(const ParametersType & lowerLimits) override + { + static_assert(SpaceDimension == 3, "This specialization is for 3D only!"); + + typedef AdvancedLimitedEuler3DTransform< TScalarType > Euler3DTransformType; + typename Euler3DTransformType::Pointer transform + = dynamic_cast< Euler3DTransformType * >( this ); + + if( transform ) + { + transform->Euler3DTransformType::SetLowerLimits( lowerLimits ); + } + } + + const ParametersType& GetLowerLimits() //const override + { + static_assert(SpaceDimension == 3, "This specialization is for 3D only!"); + + typedef AdvancedLimitedEuler3DTransform< TScalarType > Euler3DTransformType; + typename Euler3DTransformType::Pointer transform + = dynamic_cast< Euler3DTransformType * >( this ); + + if( transform ) + { + return transform->Euler3DTransformType::GetLowerLimits(); + } + return ParametersType(); + } + + const ParametersType& GetUpperLimitsReached() //const override + { + static_assert(SpaceDimension == 3, "This specialization is for 3D only!"); + + typedef AdvancedLimitedEuler3DTransform< TScalarType > Euler3DTransformType; + typename Euler3DTransformType::Pointer transform + = dynamic_cast(this); + + if (transform) + { + return transform->Euler3DTransformType::GetUpperLimitsReached(); + } + return ParametersType(); + } + + const ParametersType& GetLowerLimitsReached() //const override + { + static_assert(SpaceDimension == 3, "This specialization is for 3D only!"); + + typedef AdvancedLimitedEuler3DTransform< TScalarType > Euler3DTransformType; + typename Euler3DTransformType::Pointer transform + = dynamic_cast(this); + + if (transform) + { + return transform->Euler3DTransformType::GetLowerLimitsReached(); + } + return ParametersType(); + } + + protected: + + LimitedEulerTransform() {} + ~LimitedEulerTransform() override {} + + private: + + LimitedEulerTransform(const Self&); // purposely not implemented + void operator=(const Self&); // purposely not implemented + }; + +} // end namespace itk + +#endif // end #ifndef __itkLimitedEulerTransform_H__