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