From 7a793af5d663a7a2bf2f18d42ca50f84a3e88a8f Mon Sep 17 00:00:00 2001 From: gokhangg Date: Sat, 4 Apr 2020 12:28:48 +0200 Subject: [PATCH] ADD: Addition of point to surface metric to Elastix. COMP: Remove cmake* from .gitignore Having cmake* in there stopped Gokhan Gunay from adding his CMakeLists file to pull request #236 (Addition of PointToSurface metric). Addition of PointToSurface metric ADD: Addition of point to surface metric to Elastix. In this commit addition of point to surface metric is done. This metric essentially considers how surface of an masked organ is aligned with the points provided from the surface of the organ. Therefore, it is used to improve alignment of the organ under interest. --- .gitignore | 1 - .../PointToSurfaceDistance/CMakeLists.txt | 8 + .../elxPointToSurfaceDistanceMetric.cxx | 21 ++ .../elxPointToSurfaceDistanceMetric.h | 169 +++++++++ .../elxPointToSurfaceDistanceMetric.hxx | 217 +++++++++++ .../itkPointToSurfaceDistanceMetric.h | 130 +++++++ .../itkPointToSurfaceDistanceMetric.hxx | 341 ++++++++++++++++++ 7 files changed, 886 insertions(+), 1 deletion(-) create mode 100644 Components/Metrics/PointToSurfaceDistance/CMakeLists.txt create mode 100644 Components/Metrics/PointToSurfaceDistance/elxPointToSurfaceDistanceMetric.cxx create mode 100644 Components/Metrics/PointToSurfaceDistance/elxPointToSurfaceDistanceMetric.h create mode 100644 Components/Metrics/PointToSurfaceDistance/elxPointToSurfaceDistanceMetric.hxx create mode 100644 Components/Metrics/PointToSurfaceDistance/itkPointToSurfaceDistanceMetric.h create mode 100644 Components/Metrics/PointToSurfaceDistance/itkPointToSurfaceDistanceMetric.hxx diff --git a/.gitignore b/.gitignore index 262393bf0..1114109ba 100644 --- a/.gitignore +++ b/.gitignore @@ -33,7 +33,6 @@ # Generated files and help/ -cmake* Makefile CMakeFiles/ CTestTestfile.cmake diff --git a/Components/Metrics/PointToSurfaceDistance/CMakeLists.txt b/Components/Metrics/PointToSurfaceDistance/CMakeLists.txt new file mode 100644 index 000000000..5fc5bef01 --- /dev/null +++ b/Components/Metrics/PointToSurfaceDistance/CMakeLists.txt @@ -0,0 +1,8 @@ + +ADD_ELXCOMPONENT( PointToSurfaceDistanceMetric ON + elxPointToSurfaceDistanceMetric.h + elxPointToSurfaceDistanceMetric.hxx + elxPointToSurfaceDistanceMetric.cxx + itkPointToSurfaceDistanceMetric.h + itkPointToSurfaceDistanceMetric.hxx ) + diff --git a/Components/Metrics/PointToSurfaceDistance/elxPointToSurfaceDistanceMetric.cxx b/Components/Metrics/PointToSurfaceDistance/elxPointToSurfaceDistanceMetric.cxx new file mode 100644 index 000000000..d6be492a1 --- /dev/null +++ b/Components/Metrics/PointToSurfaceDistance/elxPointToSurfaceDistanceMetric.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 "elxPointToSurfaceDistanceMetric.h" + +elxInstallMacro( PointToSurfaceDistanceMetric ); diff --git a/Components/Metrics/PointToSurfaceDistance/elxPointToSurfaceDistanceMetric.h b/Components/Metrics/PointToSurfaceDistance/elxPointToSurfaceDistanceMetric.h new file mode 100644 index 000000000..51d0bdbee --- /dev/null +++ b/Components/Metrics/PointToSurfaceDistance/elxPointToSurfaceDistanceMetric.h @@ -0,0 +1,169 @@ +/*========================================================================= + * + * 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 __elxPointToSurfaceDistanceMetric_H__ +#define __elxPointToSurfaceDistanceMetric_H__ + +#include "elxIncludes.h" // include first to avoid MSVS warning +#include "itkPointToSurfaceDistanceMetric.h" + +namespace elastix +{ + +/** + * + * For more information check the paper:\n + * Gunay, G. , Luu, M. H., Moelker, A. , Walsum, T. and Klein, S. (2017), + * "Semiautomated registration of pre‐ and intraoperative CT for image‐guided percutaneous + * liver tumor ablation interventions," Med. Phys., 44: 3718-3725. + * + * The parameters used in this class are:\n + * \parameter Metric: Select this metric as follows: + * (Metric "PointToSurfaceDistance") + * + * \parameter PointToSurfaceDistanceAverage: A parameter to get average of the points by dividing the weight of the penalty with the number of the points included. + * example: (PointToSurfaceDistanceAverage "true") \n + * Default is "true". + * + * \Command line parameters + * \parameter -fp "file" : Elastix pointset file. + * example: + + * \parameter -dt "file" : Distance transform (file) of segmentation of the organ under interest in the moving image. + * \parameter -seg "file" : Segmentation (file) of the organ under interest in the moving image. With this parameter the algorithm internally computes the distance transform. + * \parameter -dtout "file" : If the user wants to get the distance transform from organ segmentation, this parameter is called with the name of the distance transform image file to write . + * + * + * \ingroup Metrics + * + */ + +template< class TElastix > +class PointToSurfaceDistanceMetric : + public + itk::PointToSurfaceDistanceMetric < + typename MetricBase< TElastix >::FixedPointSetType, + typename MetricBase< TElastix >::MovingPointSetType >, + public MetricBase< TElastix > +{ +public: + + using Self = PointToSurfaceDistanceMetric; + + /** Standard ITK-stuff. */ + using Superclass1 = itk::PointToSurfaceDistanceMetric::FixedPointSetType, typename MetricBase< TElastix >::MovingPointSetType >; + using Superclass2 = MetricBase< TElastix >; + using Pointer = itk::SmartPointer< Self >; + using ConstPointer = itk::SmartPointer< const Self >; + + ITK_DISALLOW_COPY_AND_ASSIGN(PointToSurfaceDistanceMetric); + /** Method for creation through the object factory. */ + itkNewMacro( Self ); + + /** Run-time type information (and related methods). */ + itkTypeMacro( PointToSurfaceDistanceMetric, itk::PointToSurfaceDistanceMetric ); + + /** Name of this class. + * Use this name in the parameter file to select this specific metric. \n + * example: (Metric "PointToSurfaceDistance")\n + */ + elxClassNameMacro( "PointToSurfaceDistance" ); + + /** Typedefs from the superclass. */ + using CoordinateRepresentationType = typename Superclass1::CoordinateRepresentationType; + using FixedPointSetType = typename Superclass1::FixedPointSetType; + using FixedPointSetConstPointer = typename Superclass1::FixedPointSetConstPointer; + using MovingPointSetType = typename Superclass1::MovingPointSetType; + using MovingPointSetConstPointer = typename Superclass1::MovingPointSetConstPointer; + using TransformType = typename Superclass1::TransformType; + using TransformPointer = typename Superclass1::TransformPointer; + using InputPointType = typename Superclass1::InputPointType; + using OutputPointType = typename Superclass1::OutputPointType; + using TransformParametersType = typename Superclass1::TransformParametersType; + using TransformJacobianType = typename Superclass1::TransformJacobianType; + using FixedImageMaskType = typename Superclass1::FixedImageMaskType; + using FixedImageMaskPointer = typename Superclass1::FixedImageMaskPointer; + using MovingImageMaskType = typename Superclass1::MovingImageMaskType; + using MovingImageMaskPointer = typename Superclass1::MovingImageMaskPointer; + using MeasureType = typename Superclass1::MeasureType; + using DerivativeType = typename Superclass1::DerivativeType; + using ParametersType = typename Superclass1::ParametersType; + + /** Typedefs inherited from elastix. */ + using ElastixType = typename Superclass2::ElastixType; + using ElastixPointer = typename Superclass2::ElastixPointer; + using ConfigurationType = typename Superclass2::ConfigurationType; + using ConfigurationPointer = typename Superclass2::ConfigurationPointer; + using RegistrationType = typename Superclass2::RegistrationType; + using RegistrationPointer = typename Superclass2::RegistrationPointer; + using ITKBaseType = typename Superclass2::ITKBaseType; + using FixedImageType = typename Superclass2::FixedImageType; + using MovingImageType = typename Superclass2::MovingImageType; + + /** The fixed image dimension. */ + itkStaticConstMacro( FixedImageDimension, unsigned int, FixedImageType::ImageDimension ); + + /** The moving image dimension. */ + itkStaticConstMacro( MovingImageDimension, unsigned int, MovingImageType::ImageDimension ); + + using PointSetType = FixedPointSetType; + using ImageType = typename Superclass1::ImageType; + + + /** Sets up a timer to measure the initialisation time and + * calls the Superclass' implementation. + */ + virtual void Initialize(); + + /** + * Do some things before all: + * \li Check and print the command line arguments fp and mp. + * This should be done in BeforeAllBase and not BeforeAll. + */ + virtual int BeforeAllBase(); + + /** + * Do some things before registration: + * \li Load and set the pointsets. + */ + virtual void BeforeRegistration(); + + using ImageConstPointer = typename Superclass1::ImageType::ConstPointer; + + unsigned int ReadLandmarks( const std::string & landmarkFileName, typename PointSetType::Pointer & pointSet, typename FixedImageType::ConstPointer image ); + + /** Overwrite to silence warning. */ + virtual void SelectNewSamples( void ){} + +protected: + + /** The constructor. */ + PointToSurfaceDistanceMetric() = default; + /** The destructor. */ + virtual ~PointToSurfaceDistanceMetric() = default; + +}; + + + +} // end namespace elastix + +#ifndef ITK_MANUAL_INSTANTIATION +#include "elxPointToSurfaceDistanceMetric.hxx" +#endif + +#endif diff --git a/Components/Metrics/PointToSurfaceDistance/elxPointToSurfaceDistanceMetric.hxx b/Components/Metrics/PointToSurfaceDistance/elxPointToSurfaceDistanceMetric.hxx new file mode 100644 index 000000000..13da14198 --- /dev/null +++ b/Components/Metrics/PointToSurfaceDistance/elxPointToSurfaceDistanceMetric.hxx @@ -0,0 +1,217 @@ +/*========================================================================= + * + * 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 __elxPointToSurfaceDistanceMetric_HXX__ +#define __elxPointToSurfaceDistanceMetric_HXX__ + +#include "elxPointToSurfaceDistanceMetric.h" +#include "itkTransformixInputPointFileReader.h" + +namespace elastix +{ + +/** + * ******************* Initialize *********************** + */ + +template< class TElastix > +void +PointToSurfaceDistanceMetric< TElastix > +::Initialize( void ) +{ +} + +/** + * ***************** BeforeAllBase *********************** + */ + +template< class TElastix > +int +PointToSurfaceDistanceMetric< TElastix > +::BeforeAllBase() +{ + this->Superclass2::BeforeAllBase(); + /** Check if the current configuration uses this metric. */ + unsigned int count = 0; + for( unsigned int i = 0; i < this->m_Configuration->CountNumberOfParameterEntries( "Metric" ); ++i ) + { + std::string metricName = ""; + this->m_Configuration->ReadParameter( metricName, "Metric", i ); + if( metricName == "PointToSurfaceDistance" ) { count++; } + } + if( count == 0 ) { return 0; } + + /** Check Command line options and print them to the log file. */ + elxout << "Command line options from PointToSurfaceDistanceMetric:" << std::endl; + + /** Check for appearance of parameter "PointToSurfaceDistanceAverage". */ + std::string PointToSurfaceDistanceAverageStr = "true"; + this->m_AvPointWeigh = true; + if (this->m_Configuration->CountNumberOfParameterEntries( "PointToSurfaceDistanceAverage" ) == 1) + { + this->m_Configuration->ReadParameter( PointToSurfaceDistanceAverageStr, "PointToSurfaceDistanceAverage", 0 ); + if (PointToSurfaceDistanceAverageStr == "false") this->m_AvPointWeigh = false; + } + + elxout << "\nAverage of points in annotation set : " + << PointToSurfaceDistanceAverageStr <<"\n" + << std::endl; + + /** Check for appearance of "-fp". */ + auto _fp = this->m_Configuration->GetCommandLineArgument( "-fp" ); + if( _fp.empty() ) + { + elxout << "-fp unspecified" << std::endl; + } + else + { + elxout << "-fp " << _fp << std::endl; + } + /** Check for appearance of "-dt". */ + auto _dt = this->m_Configuration->GetCommandLineArgument( "-dt" ); + if( _dt.empty() ) + { + elxout << "-dt unspecified" << std::endl; + } + else + { + elxout << "-dt " << _dt << std::endl; + this->Superclass1::SetDTImageIn(_dt); + } + + /** Check for appearance of "-seg". */ + auto _seg = this->m_Configuration->GetCommandLineArgument( "-seg" ); + if( _seg.empty() ) + { + elxout << "-seg unspecified" << std::endl; + } + else + { + elxout << "-seg " << _seg << std::endl; + this->Superclass1::SetSegImageIn(_seg); + } + + /** Check for appearance of "-dtout". */ + auto _dtout = this->m_Configuration->GetCommandLineArgument( "-dtout" ); + if( _dtout.empty() ) + { + elxout << "-dtout unspecified" << std::endl; + } + else + { + elxout << "-dtout " << _dtout << std::endl; + this->Superclass1::SetDTImageOut(_dtout); + } + + this->Superclass1::Initialize(); + return 0; +} + +/** + * ***************** BeforeRegistration *********************** + */ + +template< class TElastix > +void +PointToSurfaceDistanceMetric< TElastix > +::BeforeRegistration() +{ + /** Read and set the fixed pointset. */ + auto fixedPointsetFileName = this->GetConfiguration()->GetCommandLineArgument( "-fp" ); + typename PointSetType::Pointer fixedPointSet; + const typename FixedImageType::ConstPointer Image = this->GetElastix()->GetFixedImage(); + ReadLandmarks( fixedPointsetFileName, fixedPointSet, Image ); + + this->SetFixedPointSet( fixedPointSet );////this is pointset interface for the layer a code +} + +/** + * ***************** ReadLandmarks *********************** + */ + +template< class TElastix > +unsigned int +PointToSurfaceDistanceMetric< TElastix > +::ReadLandmarks( const std::string & landmarkFileName, typename PointSetType::Pointer & pointSet, const typename FixedImageType::ConstPointer image ) +{ + using PointSetReaderType = itk::TransformixInputPointFileReader; + + elxout << "Loading landmarks for " << this->GetComponentLabel() + << ":" << this->elxGetClassName() << "." << std::endl; + + /** Read the landmarks. */ + auto reader = PointSetReaderType::New(); + + reader->SetFileName( landmarkFileName.c_str() ); + elxout << " Reading landmark file: " << landmarkFileName << std::endl; + try + { + reader->Update(); + } + catch( itk::ExceptionObject & err ) + { + xl::xout[ "error" ] << " Error while opening " << landmarkFileName << std::endl; + xl::xout[ "error" ] << err << std::endl; + itkExceptionMacro( << "ERROR: unable to configure " << this->GetComponentLabel() ); + } + + /** Some user-feedback. */ + const auto nrofpoints = reader->GetNumberOfPoints(); + if( reader->GetPointsAreIndices() ) + { + elxout << " Landmarks are specified as image indices." << std::endl; + } + else + { + elxout << " Landmarks are specified in world coordinates." << std::endl; + } + elxout << " Number of specified points: " << nrofpoints << std::endl; + + /** Get the pointset. */ + pointSet = reader->GetOutput(); + /** Convert from index to point if necessary */ + pointSet->DisconnectPipeline(); + + if( reader->GetPointsAreIndices() ) + { + /** Convert to world coordinates */ + for( auto j = 0u; j < nrofpoints; ++j ) + { + /** The landmarks from the pointSet are indices. We first cast to the + * proper type, and then convert it to world coordinates. + */ + typename ImageType::PointType point; + typename ImageType::IndexType index; + pointSet->GetPoint( j, &point ); + for( auto d = 0u; d < FixedImageDimension; ++d ) + { + index[ d ] = static_cast( itk::Math::Round< double >( point[ d ] ) ); + } + + /** Compute the input point in physical coordinates. */ + image->TransformIndexToPhysicalPoint( index, point ); + pointSet->SetPoint( j, point ); + } // end for all points + } // end for points are indices + + return nrofpoints; +} + +} // end namespace elastix + +#endif // end #ifndef __elxPointToSurfaceDistanceMetric_HXX__ diff --git a/Components/Metrics/PointToSurfaceDistance/itkPointToSurfaceDistanceMetric.h b/Components/Metrics/PointToSurfaceDistance/itkPointToSurfaceDistanceMetric.h new file mode 100644 index 000000000..4cf61117e --- /dev/null +++ b/Components/Metrics/PointToSurfaceDistance/itkPointToSurfaceDistanceMetric.h @@ -0,0 +1,130 @@ +/*========================================================================= + * + * 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 __itkPointToSurfaceDistanceMetric_h +#define __itkPointToSurfaceDistanceMetric_h + +#include "itkSingleValuedPointSetToPointSetMetric.h" + +/*These are included just for the metric. Presence of these stuffs also signifies ad-hoc definition or being less generic of the metric*/ +#include "itkDanielssonDistanceMapImageFilter.h" +#include "itkImage.h" +#include "itkImageFileReader.h" +#include "itkImageFileWriter.h" +#include "itkRescaleIntensityImageFilter.h" +#include "itkInvertIntensityImageFilter.h" +#include "itkAddImageFilter.h" +#include "itkBSplineInterpolateImageFunction.h" +/*These are included just for the metric. Presence of these stuffs also signifies ad-hoc definition of the metric*/ + +namespace itk +{ + +template< class TFixedPointSet, class TMovingPointSet > +class PointToSurfaceDistanceMetric : + public SingleValuedPointSetToPointSetMetric< TFixedPointSet, TMovingPointSet > +{ +public: + using Self = PointToSurfaceDistanceMetric; + using Superclass = SingleValuedPointSetToPointSetMetric< TFixedPointSet, TMovingPointSet >; + using Pointer = SmartPointer< Self >; + using ConstPointer = SmartPointer< const Self >; + + /** Types transferred from the base class */ + using TransformPointer = typename Superclass::TransformPointer; + using TransformParametersType = typename Superclass::TransformParametersType; + using TransformJacobianType = typename Superclass::TransformJacobianType; + using MeasureType = typename Superclass::MeasureType; + using DerivativeType = typename Superclass::DerivativeType; + using DerivativeValueType = typename Superclass::DerivativeValueType; + using FixedPointSetType = typename Superclass::FixedPointSetType; + using MovingPointSetType = typename Superclass::MovingPointSetType; + using FixedPointSetConstPointer = typename Superclass::FixedPointSetConstPointer; + using MovingPointSetConstPointer = typename Superclass::MovingPointSetConstPointer; + using PointIterator = typename Superclass::PointIterator; + using PointDataIterator = typename Superclass::PointDataIterator; + using InputPointType = typename Superclass::InputPointType; + using OutputPointType = typename Superclass::OutputPointType; + using CoordRepType = typename OutputPointType::CoordRepType; + using VnlVectorType = vnl_vector< CoordRepType >; + using NonZeroJacobianIndicesType = typename Superclass::NonZeroJacobianIndicesType; + + /** Constants for the pointset dimensions. */ + itkStaticConstMacro( FixedPointSetDimension, unsigned int, TFixedPointSet::PointDimension ); + + using ImageType = itk::Image< float, itkGetStaticConstMacro( FixedPointSetDimension )>; + using InputImageType = itk::Image< unsigned char, itkGetStaticConstMacro( FixedPointSetDimension ) >; + using RescalerType = itk::RescaleIntensityImageFilter< ImageType, ImageType >; + using DtFilterType = itk::DanielssonDistanceMapImageFilter< InputImageType, ImageType, ImageType >; + using NegateFilterType = itk::InvertIntensityImageFilter< InputImageType, InputImageType >; + using AddImageFilt = itk::AddImageFilter< ImageType, ImageType, ImageType >; + using SegReaderType =itk::ImageFileReader< InputImageType >; + using DTReaderType = itk::ImageFileReader< ImageType >; + using WriterType = itk::ImageFileWriter< ImageType >; + using InterpolatorType = itk::BSplineInterpolateImageFunction< ImageType, float, float >; + using ConstDTimage = typename ImageType::ConstPointer; + /** Method for creation through the object factory. */ + itkNewMacro( Self );//OK + + /** Run-time type information (and related methods). */ + itkTypeMacro( PointToSurfaceDistanceMetric, SingleValuedPointSetToPointSetMetric );//OK + /** Initialize. */ + virtual void Initialize(); + + /** Get the value for single valued optimizers. */ + MeasureType GetValue( const TransformParametersType & parameters ) const; + + /** Get the derivatives of the match measure. */ + void GetDerivative( const TransformParametersType & parameters, DerivativeType & Derivative ) const; + + /** Get value and derivatives for multiple valued optimizers. */ + void GetValueAndDerivative( const TransformParametersType & parameters, MeasureType & Value, DerivativeType & Derivative ) const; + + /** Set input Segmented Image File **/ + void SetSegImageIn( const std::string str ) ; + + /** Set input Segmented Image File **/ + void SetDTImageIn( const std::string str ) ; + + /** Set output Distance Transform File **/ + void SetDTImageOut( const std::string str ); + + /** Get Distance Transform Image*/ + ConstDTimage GetDTImage() const; + +protected: + + PointToSurfaceDistanceMetric(); + virtual ~PointToSurfaceDistanceMetric() = default; + bool m_AvPointWeigh{true}; + +private: + + PointToSurfaceDistanceMetric( const Self & ); + std::string m_segmentationFileIn, m_distanceTransformFileIn, m_distanceTransformFileOut; + typename ImageType::Pointer m_internalDistanceTransformImage; + typename InterpolatorType::Pointer m_interpolator; + +}; + +} // end namespace itk + +#ifndef ITK_MANUAL_INSTANTIATION +#include "itkPointToSurfaceDistanceMetric.hxx" +#endif + +#endif diff --git a/Components/Metrics/PointToSurfaceDistance/itkPointToSurfaceDistanceMetric.hxx b/Components/Metrics/PointToSurfaceDistance/itkPointToSurfaceDistanceMetric.hxx new file mode 100644 index 000000000..abce21dec --- /dev/null +++ b/Components/Metrics/PointToSurfaceDistance/itkPointToSurfaceDistanceMetric.hxx @@ -0,0 +1,341 @@ +/*========================================================================= + * + * 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 __itkPointToSurfaceDistanceMetric_hxx +#define __itkPointToSurfaceDistanceMetric_hxx + +#include "itkPointToSurfaceDistanceMetric.h" +#include + +namespace itk +{ + +/** + * ******************* Constructor ******************* + */ + +template< class TFixedPointSet, class TMovingPointSet > +PointToSurfaceDistanceMetric< TFixedPointSet, TMovingPointSet > +::PointToSurfaceDistanceMetric() + : m_internalDistanceTransformImage{ nullptr } + , m_interpolator{ nullptr } + +{ +} + +template< class TFixedPointSet, class TMovingPointSet > +void PointToSurfaceDistanceMetric< TFixedPointSet, TMovingPointSet > +::Initialize() +{ + /** Initialize transform, interpolator, etc. + Superclass::Initialize(); + + /***********************ITK Stuffs******************************/ + /** Sanity checks. */ + + if ((m_segmentationFileIn.length() <= 0) && (m_distanceTransformFileIn.length() <= 0)) + { + itkExceptionMacro( << "Neither Distance Transform nor Segmentation file is not found" ); + } + + if ((m_segmentationFileIn.length() > 0) && (m_distanceTransformFileIn.length() > 0)) + { + itkExceptionMacro( << "Distance Transform and Segmentation files can not be loaded at the same time" ); + } + /***********************ITK Stuffs******************************/ + try + { + if ((m_segmentationFileIn.length() > 0)) + { + /***********************Absolute DT computation*****************/ + m_internalDistanceTransformImage = ImageType::New(); + auto DTfilter = DtFilterType::New(); + auto DTfilter2 = DtFilterType::New(); + auto NegateFilter = NegateFilterType::New(); + auto scaler = RescalerType::New(); + auto AddImage = AddImageFilt::New(); + auto segReader = SegReaderType::New(); + auto segReader2 = SegReaderType::New(); + segReader->SetFileName( m_segmentationFileIn ); + segReader2->SetFileName( m_segmentationFileIn ); + + /******This values for rescaling were put trivially***********/ + scaler->SetOutputMaximum( 10L );//These values should be changed after comprehensive experiments + scaler->SetOutputMinimum( 0L ); + /******This values for rescaling were put trivially***********/ + + segReader->Update(); + segReader2->Update(); + + DTfilter->SetInput( segReader->GetOutput() ); + NegateFilter->SetInput( segReader2->GetOutput() ); + DTfilter2->SetInput( NegateFilter->GetOutput() ); + + AddImage->SetInput1(DTfilter->GetOutput()); + AddImage->SetInput2(DTfilter2->GetOutput()); + AddImage->Update(); + + scaler->SetInput( AddImage->GetOutput() ); + scaler->Update(); + + m_internalDistanceTransformImage = scaler->GetOutput(); + m_internalDistanceTransformImage->DisconnectPipeline(); + /***********************Absolute DT computation*****************/ + } + + // Read an ADT image outside. + if ((m_distanceTransformFileIn.length()>0)) + { + m_internalDistanceTransformImage = ImageType::New(); + auto dtReader = DTReaderType::New(); + dtReader->SetFileName( m_distanceTransformFileIn ); + dtReader->Update(); + m_internalDistanceTransformImage = dtReader->GetOutput(); + } + + // Write an ADT image outside, just for checking. + if (m_distanceTransformFileOut.length()>0) + { + auto writer = WriterType::New(); + writer->SetFileName( m_distanceTransformFileOut); + writer->SetInput( m_internalDistanceTransformImage ); + writer->Update(); + } + } + catch(...) + { + itkExceptionMacro(<< "There is a problem in the creation of distance transform image."); + } + + if (m_internalDistanceTransformImage) + { + m_interpolator = InterpolatorType::New(); + m_interpolator->SetSplineOrder(3); + m_interpolator->SetInputImage(m_internalDistanceTransformImage); + } + else + { + itkExceptionMacro(<< "There is no distance transform image."); + } +} + +/** + * ******************* GetValue ******************* + */ + +template< class TFixedPointSet, class TMovingPointSet > +auto +PointToSurfaceDistanceMetric< TFixedPointSet, TMovingPointSet > +::GetValue( const TransformParametersType & parameters ) const -> MeasureType +{ + /***********************ITK Stuffs******************************/ + /** Sanity checks. */ + auto fixedPointSet = this->GetFixedPointSet(); + if( !fixedPointSet ) + { + itkExceptionMacro( << "Fixed point set has not been assigned" ); + } + + /***********************ITK Stuffs******************************/ + auto measure = NumericTraits< MeasureType >::Zero; + + /** Make sure the transform parameters are up to date. */ + this->SetTransformParameters( parameters ); + + /*********************ADT cost computation**********************/ + /** Loop over the corresponding points. */ + { + this->m_NumberOfPointsCounted = 0; + OutputPointType fixedPoint, mappedPoint; + const auto pointItFixedEnd = fixedPointSet->GetPoints()->End(); + for (auto pointItFixed = fixedPointSet->GetPoints()->Begin(); pointItFixed != pointItFixedEnd; ++pointItFixed) + { + /** Get the current corresponding points. */ + fixedPoint = pointItFixed.Value(); + + /** Transform point and check if it is inside the B-spline support region. */ + mappedPoint = this->m_Transform->TransformPoint(fixedPoint); + + double meas = m_interpolator->Evaluate(mappedPoint); + ++this->m_NumberOfPointsCounted; + measure += meas * meas; + } // end loop over all corresponding points + } + /*********************ADT cost computation**********************/ + + if ((this->m_NumberOfPointsCounted > 0) && (this->m_AvPointWeigh)) + { + measure /= this->m_NumberOfPointsCounted; + } + + return (measure); +} + +/** + * ******************* GetDerivative ******************* + */ + +template< class TFixedPointSet, class TMovingPointSet > +void +PointToSurfaceDistanceMetric< TFixedPointSet, TMovingPointSet > +::GetDerivative( const TransformParametersType & parameters, DerivativeType & derivative ) const +{ + MeasureType dummyvalue = NumericTraits< MeasureType >::Zero; + this->GetValueAndDerivative( parameters, dummyvalue, derivative ); +} + +/** + * ******************* GetValueAndDerivative ******************* + */ + +template< class TFixedPointSet, class TMovingPointSet > +void +PointToSurfaceDistanceMetric< TFixedPointSet, TMovingPointSet > +::GetValueAndDerivative( const TransformParametersType & parameters, MeasureType & value, DerivativeType & derivative ) const +{ + /***********************ITK Stuffs******************************/ + /** Sanity checks. */ + + FixedPointSetConstPointer fixedPointSet = this->GetFixedPointSet(); + if( !fixedPointSet ) + { + itkExceptionMacro( << "Fixed point set has not been assigned" ); + } + + /***********************ITK Stuffs******************************/ + + /***************Initialization of some variables****************/ + /** Initialize some variables. */ + this->m_NumberOfPointsCounted = 0; + MeasureType measure = NumericTraits< MeasureType >::Zero; + + /** Make sure the transform parameters are up to date. */ + this->SetTransformParameters( parameters ); + + this->BeforeThreadedGetValueAndDerivative( parameters ); + + derivative = DerivativeType( this->GetNumberOfParameters() ); + derivative.Fill( NumericTraits< DerivativeValueType >::Zero ); + /***************Initialization of some variables****************/ + TransformJacobianType jacobian; + typename InterpolatorType::CovariantVectorType covVector; + + { + OutputPointType fixedPoint, mappedPoint; + /** Loop over the corresponding points. */ + const auto pointEnd = fixedPointSet->GetPoints()->End(); + for (auto pointItFixed = fixedPointSet->GetPoints()->Begin(); pointItFixed != pointEnd; ++pointItFixed) + { + + NonZeroJacobianIndicesType nzji(this->m_Transform->GetNumberOfNonZeroJacobianIndices()); + /** Get the current corresponding points. */ + fixedPoint = pointItFixed.Value(); + /** Transform point and check if it is inside the B-spline support region. */ + mappedPoint = this->m_Transform->TransformPoint(fixedPoint); + + const auto meas = m_interpolator->Evaluate(mappedPoint); + const auto Coeff = 2 * meas; + measure += meas * meas;//*interpolator->Evaluate(mappedPoint); + + /** Get the TransformJacobian dT/dmu. */ + this->m_Transform->GetJacobian(fixedPoint, jacobian, nzji); + + VnlVectorType diff = fixedPoint.GetVnlVector(); + + /********Finding the image derivative***********/ + covVector = m_interpolator->EvaluateDerivative(mappedPoint); + /********Finding the image derivative***********/ + + /***********Converting to derivative Type*******/ + for (auto i = 0u; i < covVector.GetCovariantVectorDimension(); i++) + { + diff.put(0, covVector.GetElement(0)); + diff.put(1, covVector.GetElement(1)); + diff.put(2, covVector.GetElement(2)); + } + + const auto sizeJI = nzji.size(); + /***********Converting to derivative Type*******/ + if (sizeJI == this->GetNumberOfParameters()) + { + /** Loop over all Jacobians. */ + derivative += Coeff * diff * jacobian; + } + else + { + /** Only pick the nonzero Jacobians. */ + for (auto i = 0u; i < sizeJI; ++i) + { + const unsigned int index = nzji[i]; + VnlVectorType column = jacobian.get_column(i); + derivative[index] += Coeff * dot_product(diff, column); + } + } + + ++this->m_NumberOfPointsCounted; + } // end loop over all corresponding points + } + + value = measure; + /** Taking average of grad and measure*/ + if( (this->m_NumberOfPointsCounted > 0) && (this->m_AvPointWeigh)) + { + derivative /= this->m_NumberOfPointsCounted; + value /= this->m_NumberOfPointsCounted; + } + +} + +/** Set input Segmented Image File **/ +template< class TFixedPointSet, class TMovingPointSet > +void +PointToSurfaceDistanceMetric< TFixedPointSet, TMovingPointSet > +::SetSegImageIn( const std::string str ) +{ + this->m_segmentationFileIn = str; +} + +/** Set input Distance Transform Image File **/ +template< class TFixedPointSet, class TMovingPointSet > +void +PointToSurfaceDistanceMetric< TFixedPointSet, TMovingPointSet > +::SetDTImageIn( const std::string str ) +{ + this->m_distanceTransformFileIn = str; +} + +/** Set output Distance Transform File **/ +template< class TFixedPointSet, class TMovingPointSet > +void +PointToSurfaceDistanceMetric< TFixedPointSet, TMovingPointSet > +::SetDTImageOut( const std::string str ) +{ + this->m_distanceTransformFileOut = str; +} + +/** Get Distance Transform Image*/ +template< class TFixedPointSet, class TMovingPointSet > +auto +PointToSurfaceDistanceMetric< TFixedPointSet, TMovingPointSet > +::GetDTImage() const -> ConstDTimage +{ + return static_cast::ImageType::ConstPointer>(m_internalDistanceTransformImage); +} + +} // end namespace itk + +#endif // end #ifndef __itkPointToSurfaceDistanceMetric_hxx