diff --git a/Code/ITK/include/mapNormalizedTotalGradientImageToImageMetric.h b/Code/ITK/include/mapNormalizedTotalGradientImageToImageMetric.h index 3008402..5b28e8e 100644 --- a/Code/ITK/include/mapNormalizedTotalGradientImageToImageMetric.h +++ b/Code/ITK/include/mapNormalizedTotalGradientImageToImageMetric.h @@ -1,119 +1,127 @@ /*========================================================================= * * Copyright NumFOCUS * * 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 * * https://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 mapNormalizedTotalGradientImageToImageMetric_h #define mapNormalizedTotalGradientImageToImageMetric_h #include "itkImageToImageMetric.h" #include "itkPoint.h" #include "itkIndex.h" #include // For unique_ptr. namespace itk { /** * \class NormalizedTotalGradientImageToImageMetric * \brief TODO * \ingroup ITKRegistrationCommon * * \sphinx * \sphinxexample{Registration/Common/ComputeMeanSquareBetweenTwoImages,Compute Mean Squares Metric Between Two Images} * \endsphinx */ template class ITK_TEMPLATE_EXPORT NormalizedTotalGradientImageToImageMetric : public ImageToImageMetric { public: ITK_DISALLOW_COPY_AND_MOVE(NormalizedTotalGradientImageToImageMetric); /** Standard class type aliases. */ using Self = NormalizedTotalGradientImageToImageMetric; using Superclass = ImageToImageMetric; using Pointer = SmartPointer; using ConstPointer = SmartPointer; /** Method for creation through the object factory. */ itkNewMacro(Self); /** Run-time type information (and related methods). */ itkTypeMacro(NormalizedTotalGradientImageToImageMetric, ImageToImageMetric); /** Types inherited from Superclass. */ using typename Superclass::TransformType; using typename Superclass::TransformPointer; using typename Superclass::TransformJacobianType; using typename Superclass::InterpolatorType; using typename Superclass::MeasureType; using typename Superclass::DerivativeType; using typename Superclass::ParametersType; using typename Superclass::FixedImageType; using typename Superclass::MovingImageType; using typename Superclass::MovingImagePointType; using typename Superclass::FixedImageConstPointer; using typename Superclass::MovingImageConstPointer; using typename Superclass::CoordinateRepresentationType; using typename Superclass::FixedImageSampleContainer; using typename Superclass::ImageDerivativesType; using typename Superclass::WeightsValueType; using typename Superclass::IndexValueType; // Needed for evaluation of Jacobian. using typename Superclass::FixedImagePointType; /** The moving image dimension. */ static constexpr unsigned int MovingImageDimension = MovingImageType::ImageDimension; /** * Initialize the Metric by * (1) making sure that all the components are present and plugged * together correctly, * (2) prepare everything "global" needed for computation*/ void Initialize() override; /** Get the value. */ MeasureType GetValue(const ParametersType& parameters) const override; /** Get the derivatives of the match measure. */ void GetDerivative(const ParametersType& parameters, DerivativeType& derivative) const override; /** Get the value and derivatives for single valued optimizers. */ void GetValueAndDerivative(const ParametersType& parameters, MeasureType& value, DerivativeType& derivative) const override; + /** Set/Get Delta value. This value is used as the differential in the + * computation of the metric derivative using the finite differences method. */ + itkGetConstMacro(Delta, double); + itkSetMacro(Delta, double); + protected: NormalizedTotalGradientImageToImageMetric(); ~NormalizedTotalGradientImageToImageMetric() override = default; void PrintSelf(std::ostream& os, Indent indent) const override; +private: + double m_Delta; + }; } // end namespace itk #ifndef ITK_MANUAL_INSTANTIATION # include "mapNormalizedTotalGradientImageToImageMetric.hxx" #endif #endif diff --git a/Code/ITK/include/mapNormalizedTotalGradientImageToImageMetric.hxx b/Code/ITK/include/mapNormalizedTotalGradientImageToImageMetric.hxx index 7793601..c33445b 100644 --- a/Code/ITK/include/mapNormalizedTotalGradientImageToImageMetric.hxx +++ b/Code/ITK/include/mapNormalizedTotalGradientImageToImageMetric.hxx @@ -1,133 +1,234 @@ /*========================================================================= * * Copyright NumFOCUS * * 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 * * https://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 itkNormalizedTotalGradientImageToImageMetric_hxx #define itkNormalizedTotalGradientImageToImageMetric_hxx #include "itkCovariantVector.h" #include "itkImageRegionIterator.h" #include "itkImageIterator.h" #include "itkMath.h" #include "itkMakeUniqueForOverwrite.h" +#include "itkShapedNeighborhoodIterator.h" +#include "itkSubtractImageFilter.h" namespace itk { /** * Constructor */ template NormalizedTotalGradientImageToImageMetric::NormalizedTotalGradientImageToImageMetric() { this->SetComputeGradient(true); this->m_WithinThreadPreProcess = false; this->m_WithinThreadPostProcess = false; // For backward compatibility, the default behavior is to use all the pixels // in the fixed image. // This should be fixed in ITKv4 so that this metric behaves as the others. this->SetUseAllPixels(true); + + m_Delta = 0.00011; } /** * Print out internal information about this class */ template void NormalizedTotalGradientImageToImageMetric::PrintSelf(std::ostream & os, Indent indent) const { Superclass::PrintSelf(os, indent); } /** * Initialize */ template void NormalizedTotalGradientImageToImageMetric::Initialize() { this->Superclass::Initialize(); } template auto NormalizedTotalGradientImageToImageMetric::GetValue(const ParametersType & parameters) const -> MeasureType { itkDebugMacro("GetValue( " << parameters << " ) "); if (!this->m_FixedImage) { itkExceptionMacro(<< "Fixed image has not been assigned"); } if (!this->m_MovingImage) { itkExceptionMacro(<< "Moving image has not been assigned"); } // Set up the parameters in the transform this->m_Transform->SetParameters(parameters); //do the computation based on Fixed and moving images //Mask support can be skipped for now I suppose, but should be statet if not implemented in the class docu. - //TODO + auto resampleFilter = itk::ResampleImageFilter::New(); + resampleFilter->SetReferenceImage(this->m_FixedImage); + resampleFilter->UseReferenceImageOn(); + resampleFilter->SetTransform(this->m_Transform); + resampleFilter->SetInput(this->m_MovingImage); + resampleFilter->Update(); + auto transformedImage = resampleFilter->GetOutput(); + + using subtractFilterType = itk::SubtractImageFilter; + auto subtractFilter = subtractFilterType::New(); + subtractFilter->SetInput1(this->m_FixedImage.GetPointer()); + subtractFilter->SetInput2(transformedImage); + subtractFilter->Update(); + auto diffImage = subtractFilter->GetOutput(); + + auto dirDerTarget = CalculateDirectionalDerivativeSum(this->m_FixedImage.GetPointer()); + auto dirDerMoving = CalculateDirectionalDerivativeSum(transformedImage); + auto dirDerDiff = CalculateDirectionalDerivativeSum(diffImage); + + return dirDerDiff / (dirDerTarget + dirDerMoving); } /** * Get the both Value and Derivative Measure */ template void NormalizedTotalGradientImageToImageMetric::GetValueAndDerivative(const ParametersType & parameters, MeasureType & value, DerivativeType & derivative) const { if (!this->m_FixedImage) { itkExceptionMacro(<< "Fixed image has not been assigned"); } - // Set up the parameters in the transform - this->m_Transform->SetParameters(parameters); - throw "not implemented yet"; + value = this->GetValue(parameters); + this->GetDerivative(parameters, derivative); } /** * Get the match measure derivative */ template void NormalizedTotalGradientImageToImageMetric::GetDerivative(const ParametersType & parameters, DerivativeType & derivative) const { if (!this->m_FixedImage) { itkExceptionMacro(<< "Fixed image has not been assigned"); } - MeasureType value; - // call the combined version - this->GetValueAndDerivative(parameters, value, derivative); + ParametersType testPoint; + + testPoint = parameters; + + const unsigned int numberOfParameters = this->GetNumberOfParameters(); + derivative = DerivativeType(numberOfParameters); + for (unsigned int i = 0; i < numberOfParameters; ++i) + { + testPoint[i] -= m_Delta; + const MeasureType valuep0 = this->GetValue(testPoint); + testPoint[i] += 2 * m_Delta; + const MeasureType valuep1 = this->GetValue(testPoint); + derivative[i] = (valuep1 - valuep0) / (2 * m_Delta); + testPoint[i] = parameters[i]; + } } } // end namespace itk +template +double CalculateDirectionalDerivativeSum(const TImage* image) +{ + using ImageType = TImage; + using IteratorType = itk::ConstShapedNeighborhoodIterator; + IteratorType::RadiusType radius; + radius.Fill(1); + IteratorType iterator(radius, image, image->GetLargestPossibleRegion()); + IteratorType::OffsetType top, bottom, left, right; + top = { { 0,-1} }; + bottom = { { 0, 1} }; + left = { {-1, 0} }; + right = { { 1, 0} }; + iterator.ActivateOffset(top); + iterator.ActivateOffset(bottom); + iterator.ActivateOffset(left); + iterator.ActivateOffset(right); + double sum = 0; + auto regionSize = image->GetLargestPossibleRegion().GetSize(); + + for (iterator.GoToBegin(); !iterator.IsAtEnd(); ++iterator) + { + auto centerIndex = iterator.GetIndex(); + if (centerIndex[0] > 0 && centerIndex[0] < regionSize[0] - 1 && centerIndex[1] > 0 && centerIndex[1] < regionSize[1] - 1) + { + auto centerVal = *(iterator.GetElement(0)); + IteratorType::ConstIterator localIterator = iterator.Begin(); + + while (!localIterator.IsAtEnd()) + { + sum += abs(centerVal - localIterator.Get()); + ++localIterator; + } + } + } + + // iterate over diagonal neighbors separately to include /SQRT2 factor + IteratorType::OffsetType topLeft, topRight, bottomLeft, bottomRight; + topLeft = { {-1,-1} }; + topRight = { { 1, 1} }; + bottomLeft = { {-1, 1} }; + bottomRight = { { 1, 1} }; + iterator.ActivateOffset(topLeft); + iterator.ActivateOffset(topRight); + iterator.ActivateOffset(bottomLeft); + iterator.ActivateOffset(bottomRight); + + const double SQRT2 = sqrt(2); + for (iterator.GoToBegin(); !iterator.IsAtEnd(); ++iterator) + { + auto centerIndex = iterator.GetIndex(); + if (centerIndex[0] > 0 && centerIndex[0] < regionSize[0] - 1 && centerIndex[1] > 0 && centerIndex[1] < regionSize[1] - 1) + { + auto centerVal = *(iterator.GetElement(0)); + IteratorType::ConstIterator localIterator = iterator.Begin(); + + while (!localIterator.IsAtEnd()) + { + sum += abs(centerVal - localIterator.Get()) / SQRT2; + ++localIterator; + } + } + } + + return sum; +} + #endif