diff --git a/.gitattributes b/.gitattributes new file mode 100644 index 0000000..67f1724 --- /dev/null +++ b/.gitattributes @@ -0,0 +1,22 @@ +# Custom attribute to mark sources as using our C++/C code style. +[attr]our-c-style whitespace=tab-in-indent,no-lf-at-eof hooks.style=KWStyle,clangformat +*.c our-c-style +*.h our-c-style +*.cxx our-c-style +*.hxx our-c-style +*.txx our-c-style + +# Custom attribute to mark sources as using our CMake code style. +[attr]our-cmake-style whitespace=tab-in-indent,no-lf-at-eof hooks.style=cmakeformat +*.txt our-cmake-style +*.cmake our-cmake-style +*.wrap our-cmake-style +CMakeLists.txt our-cmake-style + +# ExternalData content links must have LF newlines +*.md5 crlf=input +*.sha512 crlf=input +*.cid crlf=input + +# ghostflow-director GitHub automatic check for maximum repository file size +* hooks-max-size=100000 diff --git a/.github/workflows/build-test-package.yml b/.github/workflows/build-test-package.yml index 4b5f746..8b7d3fc 100644 --- a/.github/workflows/build-test-package.yml +++ b/.github/workflows/build-test-package.yml @@ -4,9 +4,9 @@ on: [push,pull_request] jobs: cxx-build-workflow: - uses: InsightSoftwareConsortium/ITKRemoteModuleBuildTestPackageAction/.github/workflows/build-test-cxx.yml@v5.4.0 + uses: InsightSoftwareConsortium/ITKRemoteModuleBuildTestPackageAction/.github/workflows/build-test-cxx.yml@v5.4.7 python-build-workflow: - uses: InsightSoftwareConsortium/ITKRemoteModuleBuildTestPackageAction/.github/workflows/build-test-package-python.yml@v5.4.0 + uses: InsightSoftwareConsortium/ITKRemoteModuleBuildTestPackageAction/.github/workflows/build-test-package-python.yml@v5.4.7 secrets: pypi_password: ${{ secrets.pypi_password }} diff --git a/CMakeLists.txt b/CMakeLists.txt index 75c578c..65575a3 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,10 +1,10 @@ -cmake_minimum_required(VERSION 3.10.2) +cmake_minimum_required(VERSION 3.16.3) project(ThinShellDemons) #set(ThinShellDemons_LIBRARIES ThinShellDemons) if(NOT ITK_SOURCE_DIR) - find_package(ITK REQUIRED) + find_package(ITK 5.2 REQUIRED) list(APPEND CMAKE_MODULE_PATH ${ITK_CMAKE_DIR}) include(ITKModuleExternal) else() diff --git a/README.md b/README.md index 8cddc1a..cf53782 100644 --- a/README.md +++ b/README.md @@ -11,11 +11,6 @@ This module implements the Thin Shell Demons regularization proposed in > MIUA 2015 - -> :warning: **This module requires to be compiled against an ITK version with** -> - [PointSetToPointSetMetricWithIndexv4](https://github.com/InsightSoftwareConsortium/ITK/pull/2385) -> -

diff --git a/include/itkThinShellDemonsMetricv4.h b/include/itkThinShellDemonsMetricv4.h index ab0e29a..25a9912 100644 --- a/include/itkThinShellDemonsMetricv4.h +++ b/include/itkThinShellDemonsMetricv4.h @@ -55,20 +55,19 @@ namespace itk * * \ingroup ThinShellDemons */ -template< typename TFixedMesh, typename TMovingMesh = TFixedMesh, - class TInternalComputationValueType = double > -class ITK_TEMPLATE_EXPORT ThinShellDemonsMetricv4: - public PointSetToPointSetMetricWithIndexv4< TFixedMesh, TMovingMesh, TInternalComputationValueType> +template +class ITK_TEMPLATE_EXPORT ThinShellDemonsMetricv4 + : public PointSetToPointSetMetricWithIndexv4 { public: ITK_DISALLOW_COPY_AND_MOVE(ThinShellDemonsMetricv4); /** Standard class typedefs. */ - typedef ThinShellDemonsMetricv4 Self; - typedef PointSetToPointSetMetricWithIndexv4< TFixedMesh, TMovingMesh > Superclass; + typedef ThinShellDemonsMetricv4 Self; + typedef PointSetToPointSetMetricWithIndexv4 Superclass; - typedef SmartPointer< Self > Pointer; - typedef SmartPointer< const Self > ConstPointer; + typedef SmartPointer Pointer; + typedef SmartPointer ConstPointer; /** Method for creation through the object factory. */ itkNewMacro(Self); @@ -77,8 +76,8 @@ class ITK_TEMPLATE_EXPORT ThinShellDemonsMetricv4: itkTypeMacro(ThinShellDemonsMetricv4, PointSetToPointSetMetricWithIndexv4); /** Types transferred from the base class. */ - typedef typename Superclass::FixedPointSetType FixedPointSetType; - typedef typename Superclass::MovingPointSetType MovingPointSetType; + typedef typename Superclass::FixedPointSetType FixedPointSetType; + typedef typename Superclass::MovingPointSetType MovingPointSetType; /** Types transferred from the base class */ using MeasureType = typename Superclass::MeasureType; @@ -108,16 +107,20 @@ class ITK_TEMPLATE_EXPORT ThinShellDemonsMetricv4: using PointSetPointer = typename Superclass::FixedPointSetType::ConstPointer; - void Initialize(void) override; + void + Initialize(void) override; MeasureType - GetLocalNeighborhoodValueWithIndex(const PointIdentifier &, const PointType &, - const PixelType & pixel = 0) const override; + GetLocalNeighborhoodValueWithIndex(const PointIdentifier &, + const PointType &, + const PixelType & pixel = 0) const override; void - GetLocalNeighborhoodValueAndDerivativeWithIndex(const PointIdentifier &, const PointType &, - MeasureType &, LocalDerivativeType &, - const PixelType & pixel = 0) const override; + GetLocalNeighborhoodValueAndDerivativeWithIndex(const PointIdentifier &, + const PointType &, + MeasureType &, + LocalDerivativeType &, + const PixelType & pixel = 0) const override; /** * Stretching penalty weight @@ -139,7 +142,7 @@ class ITK_TEMPLATE_EXPORT ThinShellDemonsMetricv4: itkSetMacro(GeometricFeatureWeight, double); itkGetConstReferenceMacro(GeometricFeatureWeight, double); - /** + /** * Update feature match at each iteration. * * When used in conjunction with UseConfidenceWeighting and @@ -188,8 +191,8 @@ class ITK_TEMPLATE_EXPORT ThinShellDemonsMetricv4: ThinShellDemonsMetricv4(); virtual ~ThinShellDemonsMetricv4() override = default; - //Create a points locator for feature matching - using FeaturePointSetType = PointSet< double, FixedPointDimension+1>; + // Create a points locator for feature matching + using FeaturePointSetType = PointSet; using FeaturePointSetPointer = typename FeaturePointSetType::Pointer; using FeaturePointType = typename FeaturePointSetType::PointType; using FeaturePointsContainer = typename FeaturePointSetType::PointsContainer; @@ -204,31 +207,36 @@ class ITK_TEMPLATE_EXPORT ThinShellDemonsMetricv4: * * Override to use geometric features */ - virtual void InitializePointSets() const override; - void InitializeFeaturePointsLocators() const; + virtual void + InitializePointSets() const override; + void + InitializeFeaturePointsLocators() const; /** * This class uses it's own Points locators to * accomodate feature matching */ - bool RequiresMovingPointsLocator() const override + bool + RequiresMovingPointsLocator() const override { return false; }; - bool RequiresFixedPointsLocator() const override + bool + RequiresFixedPointsLocator() const override { return false; }; - void PrintSelf(std::ostream & os, Indent indent) const override; + void + PrintSelf(std::ostream & os, Indent indent) const override; private: typedef std::vector> NeighborhoodMap; - NeighborhoodMap neighborMap; + NeighborhoodMap neighborMap; - typedef std::vector< std::vector > EdgeLengthMap; - EdgeLengthMap edgeLengthMap; + typedef std::vector> EdgeLengthMap; + EdgeLengthMap edgeLengthMap; mutable MeshTypePointer fixedITKMesh; mutable MeshTypePointer movingITKMesh; @@ -236,35 +244,42 @@ class ITK_TEMPLATE_EXPORT ThinShellDemonsMetricv4: CurvatureFilterTypePointer curvature_filter; - double m_StretchWeight; - double m_BendWeight; - double m_GeometricFeatureWeight; + double m_StretchWeight; + double m_BendWeight; + double m_GeometricFeatureWeight; mutable double m_ConfidenceSigma; - bool m_UseConfidenceWeighting; - bool m_UpdateFeatureMatchingAtEachIteration; - bool m_UseMaximalDistanceConfidenceSigma; - - void FillPointAndCell(PointSetPointer &pointset, MeshTypePointer ¤tITKMesh); - double ComputeConfidenceValueAndDerivative(const VectorType &v, - VectorType &derivative) const; - void ComputeStretchAndBend(const PointIdentifier &index, - double &stretchEnergy, - double &bendEnergy, - VectorType &stretch, - VectorType &bend) const; - void ComputeNeighbors(); - void ComputeMaximalDistanceSigma() const; - FeaturePointType GetFeaturePoint(const double *v, const double &c) const; - FeaturePointType GetFeaturePoint(const PointType &v, const double &c) const; - VectorType GetMovingDirection(const PointIdentifier &identifier) const; - FeaturePointSetPointer GenerateFeaturePointSets(bool fixed) const; + bool m_UseConfidenceWeighting; + bool m_UpdateFeatureMatchingAtEachIteration; + bool m_UseMaximalDistanceConfidenceSigma; + void + FillPointAndCell(PointSetPointer & pointset, MeshTypePointer & currentITKMesh); + double + ComputeConfidenceValueAndDerivative(const VectorType & v, VectorType & derivative) const; + void + ComputeStretchAndBend(const PointIdentifier & index, + double & stretchEnergy, + double & bendEnergy, + VectorType & stretch, + VectorType & bend) const; + void + ComputeNeighbors(); + void + ComputeMaximalDistanceSigma() const; + FeaturePointType + GetFeaturePoint(const double * v, const double & c) const; + FeaturePointType + GetFeaturePoint(const PointType & v, const double & c) const; + VectorType + GetMovingDirection(const PointIdentifier & identifier) const; + FeaturePointSetPointer + GenerateFeaturePointSets(bool fixed) const; }; } // end namespace itk #ifndef ITK_MANUAL_INSTANTIATION -#include "itkThinShellDemonsMetricv4.hxx" +# include "itkThinShellDemonsMetricv4.hxx" #endif #endif diff --git a/include/itkThinShellDemonsMetricv4.hxx b/include/itkThinShellDemonsMetricv4.hxx index 61f8f07..0302e77 100644 --- a/include/itkThinShellDemonsMetricv4.hxx +++ b/include/itkThinShellDemonsMetricv4.hxx @@ -26,9 +26,8 @@ namespace itk { -template< typename TFixedMesh, typename TMovingMesh, typename TInternalComputationValueType > -ThinShellDemonsMetricv4< TFixedMesh, TMovingMesh, TInternalComputationValueType > -::ThinShellDemonsMetricv4() +template +ThinShellDemonsMetricv4::ThinShellDemonsMetricv4() { m_BendWeight = 1; m_StretchWeight = 1; @@ -39,7 +38,7 @@ ThinShellDemonsMetricv4< TFixedMesh, TMovingMesh, TInternalComputationValueType m_UseConfidenceWeighting = true; m_UpdateFeatureMatchingAtEachIteration = false; m_MovingTransformedFeaturePointsLocator = nullptr; - + fixedITKMesh = nullptr; movingITKMesh = nullptr; fixedCurvature = nullptr; @@ -48,8 +47,9 @@ ThinShellDemonsMetricv4< TFixedMesh, TMovingMesh, TInternalComputationValueType /* Set the points and cells for the mesh */ template void -ThinShellDemonsMetricv4< TFixedMesh, TMovingMesh, TInternalComputationValueType > -::FillPointAndCell(PointSetPointer &pointset, MeshTypePointer ¤tITKMesh) +ThinShellDemonsMetricv4::FillPointAndCell( + PointSetPointer & pointset, + MeshTypePointer & currentITKMesh) { /* Insert points and cells in the currentITKMesh */ for (unsigned int n = 0; n < pointset->GetNumberOfPoints(); n++) @@ -57,15 +57,15 @@ ThinShellDemonsMetricv4< TFixedMesh, TMovingMesh, TInternalComputationValueType PointType point = pointset->GetPoint(n); currentITKMesh->SetPoint(n, point); } - + for (unsigned int n = 0; n < pointset->GetNumberOfCells(); n++) { MeshCellAutoPointer tri_cell; pointset->GetCell(n, tri_cell); - // Creating a Cell from the Triangle Cell and inserting it into the Mesh + // Creating a Cell from the Triangle Cell and inserting it into the Mesh auto * triangleCell = new MeshTriangleCellType; - + itk::Array point_ids = tri_cell->GetPointIdsContainer(); for (unsigned int k = 0; k < 3; ++k) { @@ -76,15 +76,13 @@ ThinShellDemonsMetricv4< TFixedMesh, TMovingMesh, TInternalComputationValueType t_cell.TakeOwnership(triangleCell); currentITKMesh->SetCell(n, t_cell); } - } /** Initialize the metric */ -template< typename TFixedMesh, typename TMovingMesh, typename TInternalComputationValueType > +template void -ThinShellDemonsMetricv4< TFixedMesh, TMovingMesh, TInternalComputationValueType > -::Initialize() +ThinShellDemonsMetricv4::Initialize() { if (!this->m_FixedPointSet) @@ -116,11 +114,11 @@ ThinShellDemonsMetricv4< TFixedMesh, TMovingMesh, TInternalComputationValueType FillPointAndCell(this->m_FixedPointSet, this->fixedITKMesh); /* fill points and cells in moving mesh */ FillPointAndCell(this->m_MovingPointSet, this->movingITKMesh); - + /* Build the Cell Links for the ITK Mesh for calculating the neighbours*/ this->fixedITKMesh->BuildCellLinks(); this->movingITKMesh->BuildCellLinks(); - + this->curvature_filter = CurvatureFilterType::New(); /* Compute Neighbors which will be used to calculate the stretch and bend energy*/ @@ -128,98 +126,100 @@ ThinShellDemonsMetricv4< TFixedMesh, TMovingMesh, TInternalComputationValueType Superclass::Initialize(); - //Compute confidence sigma - if( this->m_UseMaximalDistanceConfidenceSigma ) - { + // Compute confidence sigma + if (this->m_UseMaximalDistanceConfidenceSigma) + { this->ComputeMaximalDistanceSigma(); - } + } } /* Iterate over all the cells in which a point belongs and get the points present in those cells*/ template void -ThinShellDemonsMetricv4< TFixedMesh, TMovingMesh, TInternalComputationValueType > -::ComputeNeighbors() +ThinShellDemonsMetricv4::ComputeNeighbors() { this->neighborMap.resize(fixedITKMesh->GetNumberOfPoints()); this->edgeLengthMap.resize(fixedITKMesh->GetNumberOfPoints()); - + for (PointIdentifier id = 0; id < fixedITKMesh->GetNumberOfPoints(); id++) { /* For iterating over the cells for a given point */ const std::set link_set = this->fixedITKMesh->GetCellLinks()->ElementAt(id); - std::set pointIdSet; + std::set pointIdSet; /* Iterate over the cells and get the neighbouring points */ - for (auto elem : link_set){ - MeshCellAutoPointer tri_cell; - this->fixedITKMesh->GetCell(elem, tri_cell); - MeshCellPointIdConstIterator point_ids = tri_cell->GetPointIds(); - for (int ik = 0; ik < 3; ++ik){ - if (point_ids[ik] != id){ - pointIdSet.insert(point_ids[ik]); - } + for (auto elem : link_set) + { + MeshCellAutoPointer tri_cell; + this->fixedITKMesh->GetCell(elem, tri_cell); + MeshCellPointIdConstIterator point_ids = tri_cell->GetPointIds(); + for (int ik = 0; ik < 3; ++ik) + { + if (point_ids[ik] != id) + { + pointIdSet.insert(point_ids[ik]); } - } + } + } // Convert Set to Vector for later use - std::vector pointIdList( pointIdSet.begin(), pointIdSet.end() ); - - //Store edge lengths + std::vector pointIdList(pointIdSet.begin(), pointIdSet.end()); + + // Store edge lengths edgeLengthMap[id].resize(pointIdList.size()); - + const PointType & p = this->m_FixedPointSet->GetPoint(id); - for (unsigned long int j=0; j < pointIdList.size(); ++j) + for (unsigned long int j = 0; j < pointIdList.size(); ++j) { - PointIdentifier nid = pointIdList[j]; - const PointType &pn = this->m_FixedPointSet->GetPoint(nid); - edgeLengthMap[id][j] = p.EuclideanDistanceTo(pn); - //Avoid division by zero - if( edgeLengthMap[id][j] < itk::NumericTraits::epsilon()) - { + PointIdentifier nid = pointIdList[j]; + const PointType & pn = this->m_FixedPointSet->GetPoint(nid); + edgeLengthMap[id][j] = p.EuclideanDistanceTo(pn); + // Avoid division by zero + if (edgeLengthMap[id][j] < itk::NumericTraits::epsilon()) + { edgeLengthMap[id][j] = itk::NumericTraits::epsilon(); - } } + } neighborMap[id] = pointIdList; - } + } } -template< typename TFixedMesh, typename TMovingMesh, typename TInternalComputationValueType > +template double -ThinShellDemonsMetricv4< TFixedMesh, TMovingMesh, TInternalComputationValueType > -::ComputeConfidenceValueAndDerivative(const VectorType &v, VectorType &derivative) const +ThinShellDemonsMetricv4::ComputeConfidenceValueAndDerivative( + const VectorType & v, + VectorType & derivative) const { double variance = m_ConfidenceSigma * m_ConfidenceSigma; double dist = v.GetSquaredNorm(); - double confidence = exp( -dist / (2*variance) ); - if( m_UpdateFeatureMatchingAtEachIteration ) - { - derivative = (-confidence/variance) * v; - } + double confidence = exp(-dist / (2 * variance)); + if (m_UpdateFeatureMatchingAtEachIteration) + { + derivative = (-confidence / variance) * v; + } return confidence; } -template< typename TFixedMesh, typename TMovingMesh, typename TInternalComputationValueType > -typename ThinShellDemonsMetricv4< TFixedMesh, TMovingMesh, TInternalComputationValueType > -::VectorType -ThinShellDemonsMetricv4< TFixedMesh, TMovingMesh, TInternalComputationValueType > -::GetMovingDirection(const PointIdentifier &identifier) const +template +typename ThinShellDemonsMetricv4::VectorType +ThinShellDemonsMetricv4::GetMovingDirection( + const PointIdentifier & identifier) const { PointType p1 = this->m_FixedPointSet->GetPoint(identifier); PointType p2 = this->m_FixedTransformedPointSet->GetPoint(identifier); return p2 - p1; } -template< typename TFixedMesh, typename TMovingMesh, typename TInternalComputationValueType > +template void -ThinShellDemonsMetricv4< TFixedMesh, TMovingMesh, TInternalComputationValueType > -::ComputeStretchAndBend( const PointIdentifier &identifier, - double &stretchEnergy, - double &bendEnergy, - VectorType &stretch, - VectorType &bend) const +ThinShellDemonsMetricv4::ComputeStretchAndBend( + const PointIdentifier & identifier, + double & stretchEnergy, + double & bendEnergy, + VectorType & stretch, + VectorType & bend) const { stretchEnergy = 0; bendEnergy = 0; @@ -228,15 +228,15 @@ ThinShellDemonsMetricv4< TFixedMesh, TMovingMesh, TInternalComputationValueType // Collect all neighbors std::vector pointIdList = this->neighborMap[identifier]; - int degree = pointIdList.size(); - VectorType v = this->GetMovingDirection(identifier); - VectorType bEnergy; + int degree = pointIdList.size(); + VectorType v = this->GetMovingDirection(identifier); + VectorType bEnergy; bEnergy.Fill(0); - for (long unsigned int i=0; i < pointIdList.size(); ++i) + for (long unsigned int i = 0; i < pointIdList.size(); ++i) { PointIdentifier neighborIdx = pointIdList[i]; - int nDegree = this->neighborMap[neighborIdx].size(); + int nDegree = this->neighborMap[neighborIdx].size(); VectorType vn = this->GetMovingDirection(neighborIdx); VectorType dx = (v - vn); @@ -244,14 +244,14 @@ ThinShellDemonsMetricv4< TFixedMesh, TMovingMesh, TInternalComputationValueType // times 4 because edge appears two times in the energy function // and the derivative has another factor of 2 from the squared norm // divided by the vertex degrees of current and nieghbor vertex - stretch += dx * (4 / (degree+nDegree)); + stretch += dx * (4 / (degree + nDegree)); stretchEnergy += dx.GetSquaredNorm(); - //Normalize bending by edge length + // Normalize bending by edge length dx /= edgeLengthMap[identifier][i]; bEnergy += dx; - bend += dx * (degree * 4 / (degree+nDegree)); - } + bend += dx * (degree * 4 / (degree + nDegree)); + } if (degree > 0) { @@ -263,14 +263,13 @@ ThinShellDemonsMetricv4< TFixedMesh, TMovingMesh, TInternalComputationValueType /* Function definition of the original method definition in itkPointSetToPointSetMetric*/ /* Performs the computation in a multi-threaded manner */ template -typename ThinShellDemonsMetricv4< TFixedMesh, TMovingMesh, TInternalComputationValueType > -::MeasureType -ThinShellDemonsMetricv4< TFixedMesh, TMovingMesh, TInternalComputationValueType > -::GetLocalNeighborhoodValueWithIndex(const PointIdentifier &identifier, - const PointType &point, - const PixelType & pixel) const +typename ThinShellDemonsMetricv4::MeasureType +ThinShellDemonsMetricv4::GetLocalNeighborhoodValueWithIndex( + const PointIdentifier & identifier, + const PointType & point, + const PixelType & pixel) const { - MeasureType value = 0; + MeasureType value = 0; LocalDerivativeType derivative; this->GetLocalNeighborhoodValueAndDerivativeWithIndex(identifier, point, value, derivative, pixel); return value; @@ -281,47 +280,46 @@ ThinShellDemonsMetricv4< TFixedMesh, TMovingMesh, TInternalComputationValueType /* This method is called inside the CalculateValueAndDerivative in itkPointSetToPointSetMetricWithIndexv4.hxx */ template void -ThinShellDemonsMetricv4< TFixedMesh, TMovingMesh, TInternalComputationValueType > -::GetLocalNeighborhoodValueAndDerivativeWithIndex(const PointIdentifier &identifier, - const PointType &point, - MeasureType &value, - LocalDerivativeType &derivative, - const PixelType & pixel) const +ThinShellDemonsMetricv4:: + GetLocalNeighborhoodValueAndDerivativeWithIndex(const PointIdentifier & identifier, + const PointType & point, + MeasureType & value, + LocalDerivativeType & derivative, + const PixelType & pixel) const { - + FeaturePointType fpoint = this->GetFeaturePoint(point, fixedCurvature->GetPointData()->ElementAt(identifier)); - + PointIdentifier mPointId = this->m_MovingTransformedFeaturePointsLocator->FindClosestPoint(fpoint); - PointType closestPoint = this->m_MovingTransformedPointSet->GetPoint(mPointId); + PointType closestPoint = this->m_MovingTransformedPointSet->GetPoint(mPointId); VectorType direction = closestPoint - point; - double dist = direction.GetSquaredNorm(); - double confidence = 1; + double dist = direction.GetSquaredNorm(); + double confidence = 1; VectorType confidenceDerivative{}; - if(this->m_UseConfidenceWeighting) - { + if (this->m_UseConfidenceWeighting) + { confidence = this->ComputeConfidenceValueAndDerivative(direction, confidenceDerivative); - } - double sE = 0; - double bE = 0; + } + double sE = 0; + double bE = 0; VectorType sD; VectorType bD; this->ComputeStretchAndBend(identifier, sE, bE, sD, bD); - VectorType dx = direction * confidence * 2 - m_StretchWeight*sD - bD * m_BendWeight; + VectorType dx = direction * confidence * 2 - m_StretchWeight * sD - bD * m_BendWeight; /* Refer to Equation 2 in the MIUA2015 paper */ value = confidence * dist + m_StretchWeight * sE + m_BendWeight * bE; - if(this->m_UseConfidenceWeighting && this->m_UpdateFeatureMatchingAtEachIteration) - { + if (this->m_UseConfidenceWeighting && this->m_UpdateFeatureMatchingAtEachIteration) + { dx += dist * confidenceDerivative; - } + } derivative = dx; } -template< typename TFixedMesh, typename TMovingMesh, typename TInternalComputationValueType > +template void -ThinShellDemonsMetricv4< TFixedMesh, TMovingMesh, TInternalComputationValueType > -::InitializePointSets() const +ThinShellDemonsMetricv4::InitializePointSets() const { /* The call to Superclass initializes the m_MovingTransformedPointSet */ Superclass::InitializePointSets(); @@ -329,18 +327,18 @@ ThinShellDemonsMetricv4< TFixedMesh, TMovingMesh, TInternalComputationValueType } -template< typename TFixedMesh, typename TMovingMesh, typename TInternalComputationValueType > -typename ThinShellDemonsMetricv4< TFixedMesh, TMovingMesh, TInternalComputationValueType >::FeaturePointSetPointer -ThinShellDemonsMetricv4< TFixedMesh, TMovingMesh, TInternalComputationValueType > -::GenerateFeaturePointSets(bool fixed) const +template +typename ThinShellDemonsMetricv4::FeaturePointSetPointer +ThinShellDemonsMetricv4::GenerateFeaturePointSets( + bool fixed) const { MeshTypePointer currentMesh; - //Update meshes according to current transforms - if(fixed) + // Update meshes according to current transforms + if (fixed) + { + for (PointIdentifier i = 0; i < this->m_FixedTransformedPointSet->GetNumberOfPoints(); i++) { - for(PointIdentifier i=0; im_FixedTransformedPointSet->GetNumberOfPoints(); i++ ) - { PointType data1 = this->m_FixedTransformedPointSet->GetPoint(i); fixedITKMesh->SetPoint(i, data1); } @@ -362,26 +360,27 @@ ThinShellDemonsMetricv4< TFixedMesh, TMovingMesh, TInternalComputationValueType curvature_filter->SetCurvatureTypeToGaussian(); curvature_filter->Compute(); auto curvature_output = curvature_filter->GetGaussCurvatureData(); - - FeaturePointSetPointer features = FeaturePointSetType::New(); - if( fixed ) + FeaturePointSetPointer features = FeaturePointSetType::New(); + + if (fixed) + { + /* Instantiate first time and re-use it for later iterations */ + if (!this->fixedCurvature) { - /* Instantiate first time and re-use it for later iterations */ - if(!this->fixedCurvature){ - this->fixedCurvature = MeshType::New(); - PointDataContainerPointer pointData = PointDataContainer::New(); - pointData->Reserve(currentMesh->GetNumberOfPoints()); - this->fixedCurvature->SetPointData(pointData); - } + this->fixedCurvature = MeshType::New(); + PointDataContainerPointer pointData = PointDataContainer::New(); + pointData->Reserve(currentMesh->GetNumberOfPoints()); + this->fixedCurvature->SetPointData(pointData); + } - for (PointIdentifier i = 0; i < currentMesh->GetNumberOfPoints(); i++) - { - this->fixedCurvature->SetPointData(i, curvature_output->GetElement(i)); - } + for (PointIdentifier i = 0; i < currentMesh->GetNumberOfPoints(); i++) + { + this->fixedCurvature->SetPointData(i, curvature_output->GetElement(i)); } + } else - { + { auto fPoints = features->GetPoints(); for (PointIdentifier i = 0; i < currentMesh->GetNumberOfPoints(); i++) { @@ -393,43 +392,39 @@ ThinShellDemonsMetricv4< TFixedMesh, TMovingMesh, TInternalComputationValueType return features; } -template< typename TFixedMesh, typename TMovingMesh, typename TInternalComputationValueType > +template void -ThinShellDemonsMetricv4< TFixedMesh, TMovingMesh, TInternalComputationValueType > -::ComputeMaximalDistanceSigma() - const +ThinShellDemonsMetricv4::ComputeMaximalDistanceSigma() const { - FeaturePointsContainerPointer mpoints = - this->m_MovingTransformedFeaturePointsLocator->GetPoints(); - double maximalDistance = 0; + FeaturePointsContainerPointer mpoints = this->m_MovingTransformedFeaturePointsLocator->GetPoints(); + double maximalDistance = 0; for (PointIdentifier i = 0; i < fixedITKMesh->GetNumberOfPoints(); i++) { - FeaturePointType fpoint = this->GetFeaturePoint(fixedITKMesh->GetPoint(i), fixedCurvature->GetPointData()->ElementAt(i)); - PointIdentifier id = this->m_MovingTransformedFeaturePointsLocator->FindClosestPoint(fpoint); + FeaturePointType fpoint = + this->GetFeaturePoint(fixedITKMesh->GetPoint(i), fixedCurvature->GetPointData()->ElementAt(i)); + PointIdentifier id = this->m_MovingTransformedFeaturePointsLocator->FindClosestPoint(fpoint); FeaturePointType cpoint = mpoints->GetElement(id); - double dist = cpoint.SquaredEuclideanDistanceTo(fpoint); - if( dist > maximalDistance ) + double dist = cpoint.SquaredEuclideanDistanceTo(fpoint); + if (dist > maximalDistance) { maximalDistance = dist; - } } - this->m_ConfidenceSigma = sqrt(maximalDistance)/3; + } + this->m_ConfidenceSigma = sqrt(maximalDistance) / 3; } -template< typename TFixedMesh, typename TMovingMesh, typename TInternalComputationValueType > +template void -ThinShellDemonsMetricv4< TFixedMesh, TMovingMesh, TInternalComputationValueType > -::InitializeFeaturePointsLocators() - const +ThinShellDemonsMetricv4::InitializeFeaturePointsLocators() const { - //Update fixed curvature - if(!fixedCurvature || this->m_UpdateFeatureMatchingAtEachIteration){ + // Update fixed curvature + if (!fixedCurvature || this->m_UpdateFeatureMatchingAtEachIteration) + { this->GenerateFeaturePointSets(true); } - //Update moving curvature feature locator - if( !this->m_MovingTransformedFeaturePointsLocator - || this->m_UpdateFeatureMatchingAtEachIteration ) + // Update moving curvature feature locator + if (!this->m_MovingTransformedFeaturePointsLocator || this->m_UpdateFeatureMatchingAtEachIteration) { if (!this->m_MovingTransformedPointSet) { @@ -442,12 +437,11 @@ ThinShellDemonsMetricv4< TFixedMesh, TMovingMesh, TInternalComputationValueType // Only for the moving mesh, pass false to the GenerateFeaturePointSets FeaturePointSetPointer features = this->GenerateFeaturePointSets(false); - this->m_MovingTransformedFeaturePointsLocator->SetPoints( - features->GetPoints()); + this->m_MovingTransformedFeaturePointsLocator->SetPoints(features->GetPoints()); this->m_MovingTransformedFeaturePointsLocator->Initialize(); } - //Compute confidence sigma + // Compute confidence sigma /* if( this->m_UpdateFeatureMatchingAtEachIteration && this->m_UseMaximalDistanceConfidenceSigma ) @@ -458,41 +452,39 @@ ThinShellDemonsMetricv4< TFixedMesh, TMovingMesh, TInternalComputationValueType } /* returns point with values [x, y, z, feature] */ -template< typename TFixedMesh, typename TMovingMesh, typename TInternalComputationValueType > -typename ThinShellDemonsMetricv4< TFixedMesh, TMovingMesh, TInternalComputationValueType > -::FeaturePointType -ThinShellDemonsMetricv4< TFixedMesh, TMovingMesh, TInternalComputationValueType > -::GetFeaturePoint(const double *v, const double &c) const +template +typename ThinShellDemonsMetricv4::FeaturePointType +ThinShellDemonsMetricv4::GetFeaturePoint(const double * v, + const double & c) const { FeaturePointType fpoint; - for(unsigned int i=0; i -typename ThinShellDemonsMetricv4< TFixedMesh, TMovingMesh, TInternalComputationValueType > -::FeaturePointType -ThinShellDemonsMetricv4< TFixedMesh, TMovingMesh, TInternalComputationValueType > -::GetFeaturePoint(const PointType &v, const double &c) const +template +typename ThinShellDemonsMetricv4::FeaturePointType +ThinShellDemonsMetricv4::GetFeaturePoint(const PointType & v, + const double & c) const { FeaturePointType fpoint; - for(unsigned int i=0; i +template void -ThinShellDemonsMetricv4< TFixedMesh, TMovingMesh, TInternalComputationValueType > -::PrintSelf(std::ostream & os, Indent indent) const +ThinShellDemonsMetricv4::PrintSelf(std::ostream & os, + Indent indent) const { Superclass::PrintSelf(os, indent); } diff --git a/test/itkThinShellDemonsTestv4_Affine.cxx b/test/itkThinShellDemonsTestv4_Affine.cxx index 5e36adc..aab9647 100644 --- a/test/itkThinShellDemonsTestv4_Affine.cxx +++ b/test/itkThinShellDemonsTestv4_Affine.cxx @@ -27,47 +27,52 @@ #include "itkMeshFileReader.h" #include "itkMeshFileWriter.h" -template +template class CommandIterationUpdate : public itk::Command { public: - typedef CommandIterationUpdate Self; - typedef itk::Command Superclass; - typedef itk::SmartPointer Pointer; - itkNewMacro( Self ); + typedef CommandIterationUpdate Self; + typedef itk::Command Superclass; + typedef itk::SmartPointer Pointer; + itkNewMacro(Self); + protected: CommandIterationUpdate() {}; + public: - void Execute(itk::Object *caller, const itk::EventObject & event) override - { - Execute( (const itk::Object *) caller, event); - } + void + Execute(itk::Object * caller, const itk::EventObject & event) override + { + Execute((const itk::Object *)caller, event); + } - void Execute(const itk::Object * object, const itk::EventObject & event) override + void + Execute(const itk::Object * object, const itk::EventObject & event) override + { + if (typeid(event) != typeid(itk::IterationEvent)) { - if( typeid( event ) != typeid( itk::IterationEvent ) ) - { return; - } - const auto * optimizer = dynamic_cast< const TFilter * >( object ); + } + const auto * optimizer = dynamic_cast(object); - if( !optimizer ) - { - itkGenericExceptionMacro( "Error dynamic_cast failed" ); - } + if (!optimizer) + { + itkGenericExceptionMacro("Error dynamic_cast failed"); + } std::cout << "It: " << optimizer->GetCurrentIteration(); std::cout << " metric value: " << optimizer->GetCurrentMetricValue(); std::cout << std::endl; - } + } }; -int itkThinShellDemonsTestv4_Affine( int args, char *argv []) +int +itkThinShellDemonsTestv4_Affine(int args, char * argv[]) { const unsigned int Dimension = 3; - + using MeshType = itk::Mesh; using PointsContainerPointer = MeshType::PointsContainerPointer; - + using ReaderType = itk::MeshFileReader; using WriterType = itk::MeshFileWriter; @@ -82,7 +87,7 @@ int itkThinShellDemonsTestv4_Affine( int args, char *argv []) { fixedPolyDataReader->Update(); } - catch( itk::ExceptionObject & excp ) + catch (itk::ExceptionObject & excp) { std::cerr << "Error during Fixed Mesh Update() " << std::endl; std::cerr << excp << std::endl; @@ -93,13 +98,13 @@ int itkThinShellDemonsTestv4_Affine( int args, char *argv []) /* Initialize moving mesh polydata reader */ - ReaderType::Pointer movingPolyDataReader = ReaderType::New(); + ReaderType::Pointer movingPolyDataReader = ReaderType::New(); movingPolyDataReader->SetFileName(argv[2]); try { movingPolyDataReader->Update(); } - catch( itk::ExceptionObject & excp ) + catch (itk::ExceptionObject & excp) { std::cerr << "Error during Moving Mesh Update() " << std::endl; std::cerr << excp << std::endl; @@ -116,43 +121,43 @@ int itkThinShellDemonsTestv4_Affine( int args, char *argv []) using MovingImageType = itk::Image; - FixedImageType::SizeType fixedImageSize; - FixedImageType::PointType fixedImageOrigin; + FixedImageType::SizeType fixedImageSize; + FixedImageType::PointType fixedImageOrigin; FixedImageType::DirectionType fixedImageDirection; - FixedImageType::SpacingType fixedImageSpacing; + FixedImageType::SpacingType fixedImageSpacing; using PointIdentifier = MeshType::PointIdentifier; using BoundingBoxType = itk::BoundingBox; BoundingBoxType::Pointer boundingBox = BoundingBoxType::New(); - PointsContainerPointer points = movingMesh->GetPoints(); + PointsContainerPointer points = movingMesh->GetPoints(); boundingBox->SetPoints(points); boundingBox->ComputeBoundingBox(); typename BoundingBoxType::PointType minBounds = boundingBox->GetMinimum(); typename BoundingBoxType::PointType maxBounds = boundingBox->GetMaximum(); - int imageDiagonal = 5; + int imageDiagonal = 5; double spacing = sqrt(boundingBox->GetDiagonalLength2()) / imageDiagonal; - auto diff = maxBounds - minBounds; - fixedImageSize[0] = ceil( 1.2 * diff[0] / spacing ); - fixedImageSize[1] = ceil( 1.2 * diff[1] / spacing ); - fixedImageSize[2] = ceil( 1.2 * diff[2] / spacing ); + auto diff = maxBounds - minBounds; + fixedImageSize[0] = ceil(1.2 * diff[0] / spacing); + fixedImageSize[1] = ceil(1.2 * diff[1] / spacing); + fixedImageSize[2] = ceil(1.2 * diff[2] / spacing); fixedImageOrigin[0] = minBounds[0] - 0.1 * diff[0]; fixedImageOrigin[1] = minBounds[1] - 0.1 * diff[1]; fixedImageOrigin[2] = minBounds[2] - 0.1 * diff[2]; fixedImageDirection.SetIdentity(); - fixedImageSpacing.Fill( spacing ); + fixedImageSpacing.Fill(spacing); FixedImageType::Pointer fixedImage = FixedImageType::New(); - fixedImage->SetRegions( fixedImageSize ); - fixedImage->SetOrigin( fixedImageOrigin ); - fixedImage->SetDirection( fixedImageDirection ); - fixedImage->SetSpacing( fixedImageSpacing ); + fixedImage->SetRegions(fixedImageSize); + fixedImage->SetOrigin(fixedImageOrigin); + fixedImage->SetDirection(fixedImageDirection); + fixedImage->SetSpacing(fixedImageSpacing); fixedImage->Allocate(); using TransformType = itk::AffineTransform; TransformType::Pointer transform = TransformType::New(); transform->SetIdentity(); - transform->SetCenter(minBounds + (maxBounds - minBounds)/2); + transform->SetCenter(minBounds + (maxBounds - minBounds) / 2); using PointSetMetricType = itk::ThinShellDemonsMetricv4; @@ -164,21 +169,21 @@ int itkThinShellDemonsTestv4_Affine( int args, char *argv []) metric->UseMaximalDistanceConfidenceSigmaOn(); metric->UpdateFeatureMatchingAtEachIterationOff(); metric->SetMovingTransform(transform); - //Reversed due to using points instead of an image - //to keep semantics the same as in itkThinShellDemonsTest.cxx - //For the ThinShellDemonsMetricv4 the fixed mesh is - //regularized + // Reversed due to using points instead of an image + // to keep semantics the same as in itkThinShellDemonsTest.cxx + // For the ThinShellDemonsMetricv4 the fixed mesh is + // regularized metric->SetFixedPointSet(movingMesh); metric->SetMovingPointSet(fixedMesh); metric->SetVirtualDomainFromImage(fixedImage); metric->Initialize(); // Scales estimator - using ScalesType = itk::RegistrationParameterScalesFromPhysicalShift< PointSetMetricType >; + using ScalesType = itk::RegistrationParameterScalesFromPhysicalShift; ScalesType::Pointer shiftScaleEstimator = ScalesType::New(); shiftScaleEstimator->SetMetric(metric); // Needed with pointset metrics - shiftScaleEstimator->SetVirtualDomainPointSet( metric->GetVirtualTransformedPointSet() ); + shiftScaleEstimator->SetVirtualDomainPointSet(metric->GetVirtualTransformedPointSet()); // optimizer @@ -191,19 +196,19 @@ int itkThinShellDemonsTestv4_Affine( int args, char *argv []) optimizer->SetScalesEstimator( shiftScaleEstimator ); */ typedef itk::ConjugateGradientLineSearchOptimizerv4 OptimizerType; - OptimizerType::Pointer optimizer = OptimizerType::New(); - optimizer->SetNumberOfIterations( numberOfIterations ); - optimizer->SetScalesEstimator( shiftScaleEstimator ); - optimizer->SetMaximumStepSizeInPhysicalUnits( 0.5 ); - optimizer->SetMinimumConvergenceValue( 0.0 ); - optimizer->SetConvergenceWindowSize( 10 ); + OptimizerType::Pointer optimizer = OptimizerType::New(); + optimizer->SetNumberOfIterations(numberOfIterations); + optimizer->SetScalesEstimator(shiftScaleEstimator); + optimizer->SetMaximumStepSizeInPhysicalUnits(0.5); + optimizer->SetMinimumConvergenceValue(0.0); + optimizer->SetConvergenceWindowSize(10); using CommandType = CommandIterationUpdate; CommandType::Pointer observer = CommandType::New(); - optimizer->AddObserver( itk::IterationEvent(), observer ); + optimizer->AddObserver(itk::IterationEvent(), observer); - using AffineRegistrationType = itk::ImageRegistrationMethodv4; + using AffineRegistrationType = + itk::ImageRegistrationMethodv4; AffineRegistrationType::Pointer registration = AffineRegistrationType::New(); registration->SetNumberOfLevels(1); registration->SetObjectName("registration"); @@ -215,14 +220,14 @@ int itkThinShellDemonsTestv4_Affine( int args, char *argv []) std::cout << "Start Value= " << metric->GetValue() << std::endl; try - { + { registration->Update(); - } - catch( itk::ExceptionObject &e ) - { + } + catch (itk::ExceptionObject & e) + { std::cerr << "Exception caught: " << e << std::endl; return EXIT_FAILURE; - } + } TransformType::Pointer tx = registration->GetModifiableTransform(); metric->SetTransform(tx); diff --git a/test/itkThinShellDemonsTestv4_Displacement.cxx b/test/itkThinShellDemonsTestv4_Displacement.cxx index 2ebac44..5fd06b2 100644 --- a/test/itkThinShellDemonsTestv4_Displacement.cxx +++ b/test/itkThinShellDemonsTestv4_Displacement.cxx @@ -21,53 +21,58 @@ #include "itkThinShellDemonsMetricv4.h" #include "itkImageRegistrationMethodv4.h" #include -//#include "itkMeshDisplacementTransform.h" +// #include "itkMeshDisplacementTransform.h" #include "itkConjugateGradientLineSearchOptimizerv4.h" #include "itkLBFGS2Optimizerv4.h" #include "itkMesh.h" #include "itkMeshFileReader.h" #include "itkMeshFileWriter.h" -template +template class CommandIterationUpdate : public itk::Command { public: - typedef CommandIterationUpdate Self; - typedef itk::Command Superclass; - typedef itk::SmartPointer Pointer; - itkNewMacro( Self ); + typedef CommandIterationUpdate Self; + typedef itk::Command Superclass; + typedef itk::SmartPointer Pointer; + itkNewMacro(Self); + protected: CommandIterationUpdate() {}; + public: - void Execute(itk::Object *caller, const itk::EventObject & event) override - { - Execute( (const itk::Object *) caller, event); - } + void + Execute(itk::Object * caller, const itk::EventObject & event) override + { + Execute((const itk::Object *)caller, event); + } - void Execute(const itk::Object * object, const itk::EventObject & event) override + void + Execute(const itk::Object * object, const itk::EventObject & event) override + { + if (typeid(event) != typeid(itk::IterationEvent)) { - if( typeid( event ) != typeid( itk::IterationEvent ) ) - { return; - } - const auto * optimizer = dynamic_cast< const TFilter * >( object ); + } + const auto * optimizer = dynamic_cast(object); - if( !optimizer ) - { - itkGenericExceptionMacro( "Error dynamic_cast failed" ); - } + if (!optimizer) + { + itkGenericExceptionMacro("Error dynamic_cast failed"); + } std::cout << "It: " << optimizer->GetCurrentIteration(); std::cout << " metric value: " << optimizer->GetCurrentMetricValue(); std::cout << std::endl; - } + } }; -int itkThinShellDemonsTestv4_Displacement( int args, char *argv []) +int +itkThinShellDemonsTestv4_Displacement(int args, char * argv[]) { const unsigned int Dimension = 3; using MeshType = itk::Mesh; using PointsContainerPointer = MeshType::PointsContainerPointer; - + using ReaderType = itk::MeshFileReader; using WriterType = itk::MeshFileWriter; @@ -80,7 +85,7 @@ int itkThinShellDemonsTestv4_Displacement( int args, char *argv []) { fixedPolyDataReader->Update(); } - catch( itk::ExceptionObject & excp ) + catch (itk::ExceptionObject & excp) { std::cerr << "Error during Fixed Mesh Update() " << std::endl; std::cerr << excp << std::endl; @@ -91,13 +96,13 @@ int itkThinShellDemonsTestv4_Displacement( int args, char *argv []) /* Initialize moving mesh polydata reader */ - ReaderType::Pointer movingPolyDataReader = ReaderType::New(); + ReaderType::Pointer movingPolyDataReader = ReaderType::New(); movingPolyDataReader->SetFileName(argv[2]); try { movingPolyDataReader->Update(); } - catch( itk::ExceptionObject & excp ) + catch (itk::ExceptionObject & excp) { std::cerr << "Error during Moving Mesh Update() " << std::endl; std::cerr << excp << std::endl; @@ -114,52 +119,52 @@ int itkThinShellDemonsTestv4_Displacement( int args, char *argv []) using MovingImageType = itk::Image; - FixedImageType::SizeType fixedImageSize; - FixedImageType::PointType fixedImageOrigin; + FixedImageType::SizeType fixedImageSize; + FixedImageType::PointType fixedImageOrigin; FixedImageType::DirectionType fixedImageDirection; - FixedImageType::SpacingType fixedImageSpacing; + FixedImageType::SpacingType fixedImageSpacing; using PointIdentifier = MeshType::PointIdentifier; using BoundingBoxType = itk::BoundingBox; BoundingBoxType::Pointer boundingBox = BoundingBoxType::New(); - PointsContainerPointer points = movingMesh->GetPoints(); + PointsContainerPointer points = movingMesh->GetPoints(); boundingBox->SetPoints(points); boundingBox->ComputeBoundingBox(); typename BoundingBoxType::PointType minBounds = boundingBox->GetMinimum(); typename BoundingBoxType::PointType maxBounds = boundingBox->GetMaximum(); - int imageDiagonal = 200; + int imageDiagonal = 200; double spacing = sqrt(boundingBox->GetDiagonalLength2()) / imageDiagonal; - auto diff = maxBounds - minBounds; - fixedImageSize[0] = ceil( 1.2 * diff[0] / spacing ); - fixedImageSize[1] = ceil( 1.2 * diff[1] / spacing ); - fixedImageSize[2] = ceil( 1.2 * diff[2] / spacing ); - fixedImageOrigin[0] = minBounds[0] - 0.1*diff[0]; - fixedImageOrigin[1] = minBounds[1] - 0.1*diff[1]; - fixedImageOrigin[2] = minBounds[2] - 0.1*diff[2]; + auto diff = maxBounds - minBounds; + fixedImageSize[0] = ceil(1.2 * diff[0] / spacing); + fixedImageSize[1] = ceil(1.2 * diff[1] / spacing); + fixedImageSize[2] = ceil(1.2 * diff[2] / spacing); + fixedImageOrigin[0] = minBounds[0] - 0.1 * diff[0]; + fixedImageOrigin[1] = minBounds[1] - 0.1 * diff[1]; + fixedImageOrigin[2] = minBounds[2] - 0.1 * diff[2]; fixedImageDirection.SetIdentity(); - fixedImageSpacing.Fill( spacing ); + fixedImageSpacing.Fill(spacing); FixedImageType::Pointer fixedImage = FixedImageType::New(); - fixedImage->SetRegions( fixedImageSize ); - fixedImage->SetOrigin( fixedImageOrigin ); - fixedImage->SetDirection( fixedImageDirection ); - fixedImage->SetSpacing( fixedImageSpacing ); + fixedImage->SetRegions(fixedImageSize); + fixedImage->SetOrigin(fixedImageOrigin); + fixedImage->SetDirection(fixedImageDirection); + fixedImage->SetSpacing(fixedImageSpacing); fixedImage->Allocate(); using TransformType = itk::DisplacementFieldTransform; auto transform = TransformType::New(); - using DisplacementFieldType = TransformType::DisplacementFieldType; + using DisplacementFieldType = TransformType::DisplacementFieldType; DisplacementFieldType::Pointer field = DisplacementFieldType::New(); - field->SetRegions( fixedImageSize ); - field->SetOrigin( fixedImageOrigin ); - field->SetDirection( fixedImageDirection ); - field->SetSpacing( fixedImageSpacing ); + field->SetRegions(fixedImageSize); + field->SetOrigin(fixedImageOrigin); + field->SetDirection(fixedImageDirection); + field->SetSpacing(fixedImageSpacing); field->Allocate(); transform->SetDisplacementField(field); - using PointSetMetricType = itk::ThinShellDemonsMetricv4 ; + using PointSetMetricType = itk::ThinShellDemonsMetricv4; PointSetMetricType::Pointer metric = PointSetMetricType::New(); metric->SetStretchWeight(1); metric->SetBendWeight(5); @@ -167,22 +172,22 @@ int itkThinShellDemonsTestv4_Displacement( int args, char *argv []) metric->UseConfidenceWeightingOn(); metric->UseMaximalDistanceConfidenceSigmaOn(); metric->UpdateFeatureMatchingAtEachIterationOn(); - metric->SetMovingTransform( transform ); - //Reversed due to using points instead of an image - //to keep semantics the same as in itkThinShellDemonsTest.cxx - //For the ThinShellDemonsMetricv4 the fixed mesh is - //regularized - metric->SetFixedPointSet( movingMesh ); - metric->SetMovingPointSet( fixedMesh ); - metric->SetVirtualDomainFromImage( fixedImage ); + metric->SetMovingTransform(transform); + // Reversed due to using points instead of an image + // to keep semantics the same as in itkThinShellDemonsTest.cxx + // For the ThinShellDemonsMetricv4 the fixed mesh is + // regularized + metric->SetFixedPointSet(movingMesh); + metric->SetMovingPointSet(fixedMesh); + metric->SetVirtualDomainFromImage(fixedImage); metric->Initialize(); // Scales estimator - using ScalesType = itk::RegistrationParameterScalesFromPhysicalShift< PointSetMetricType >; + using ScalesType = itk::RegistrationParameterScalesFromPhysicalShift; ScalesType::Pointer shiftScaleEstimator = ScalesType::New(); - shiftScaleEstimator->SetMetric( metric ); + shiftScaleEstimator->SetMetric(metric); // Needed with pointset metrics - shiftScaleEstimator->SetVirtualDomainPointSet( metric->GetVirtualTransformedPointSet() ); + shiftScaleEstimator->SetVirtualDomainPointSet(metric->GetVirtualTransformedPointSet()); // optimizer @@ -190,26 +195,26 @@ int itkThinShellDemonsTestv4_Displacement( int args, char *argv []) // Does currently not support local transform // but change requested in: // https://github.com/InsightSoftwareConsortium/ITK/pull/2372 -/* - typedef itk::LBFGS2Optimizerv4 OptimizerType; - OptimizerType::Pointer optimizer = OptimizerType::New(); - optimizer->SetScalesEstimator( shiftScaleEstimator ); -*/ + /* + typedef itk::LBFGS2Optimizerv4 OptimizerType; + OptimizerType::Pointer optimizer = OptimizerType::New(); + optimizer->SetScalesEstimator( shiftScaleEstimator ); + */ typedef itk::ConjugateGradientLineSearchOptimizerv4 OptimizerType; - OptimizerType::Pointer optimizer = OptimizerType::New(); - optimizer->SetNumberOfIterations( 50 ); - optimizer->SetScalesEstimator( shiftScaleEstimator ); - optimizer->SetMaximumStepSizeInPhysicalUnits( 0.5 ); - optimizer->SetMinimumConvergenceValue( 0.0 ); - optimizer->SetConvergenceWindowSize( 10 ); + OptimizerType::Pointer optimizer = OptimizerType::New(); + optimizer->SetNumberOfIterations(50); + optimizer->SetScalesEstimator(shiftScaleEstimator); + optimizer->SetMaximumStepSizeInPhysicalUnits(0.5); + optimizer->SetMinimumConvergenceValue(0.0); + optimizer->SetConvergenceWindowSize(10); using CommandType = CommandIterationUpdate; CommandType::Pointer observer = CommandType::New(); - optimizer->AddObserver( itk::IterationEvent(), observer ); + optimizer->AddObserver(itk::IterationEvent(), observer); - using RegistrationType = itk::ImageRegistrationMethodv4; + using RegistrationType = + itk::ImageRegistrationMethodv4; RegistrationType::Pointer registration = RegistrationType::New(); registration->SetNumberOfLevels(1); registration->SetObjectName("registration"); @@ -221,14 +226,14 @@ int itkThinShellDemonsTestv4_Displacement( int args, char *argv []) std::cout << "Start Value= " << metric->GetValue() << std::endl; try - { + { registration->Update(); - } - catch( itk::ExceptionObject &e ) - { + } + catch (itk::ExceptionObject & e) + { std::cerr << "Exception caught: " << e << std::endl; return EXIT_FAILURE; - } + } TransformType::Pointer tx = registration->GetModifiableTransform(); metric->SetTransform(tx); diff --git a/test/itkThinShellDemonsTestv4_SyN.cxx b/test/itkThinShellDemonsTestv4_SyN.cxx index 637bda1..6b50ea7 100644 --- a/test/itkThinShellDemonsTestv4_SyN.cxx +++ b/test/itkThinShellDemonsTestv4_SyN.cxx @@ -32,12 +32,13 @@ * to comnine the thin shell regularization with for * example the SyN diffeomoprhic transformations. */ -int itkThinShellDemonsTestv4_SyN( int args, char *argv []) +int +itkThinShellDemonsTestv4_SyN(int args, char * argv[]) { const unsigned int Dimension = 3; using MeshType = itk::Mesh; using PointsContainerPointer = MeshType::PointsContainerPointer; - + using ReaderType = itk::MeshFileReader; using WriterType = itk::MeshFileWriter; @@ -50,7 +51,7 @@ int itkThinShellDemonsTestv4_SyN( int args, char *argv []) { fixedPolyDataReader->Update(); } - catch( itk::ExceptionObject & excp ) + catch (itk::ExceptionObject & excp) { std::cerr << "Error during Fixed Mesh Update() " << std::endl; std::cerr << excp << std::endl; @@ -61,13 +62,13 @@ int itkThinShellDemonsTestv4_SyN( int args, char *argv []) /* Initialize moving mesh polydata reader */ - ReaderType::Pointer movingPolyDataReader = ReaderType::New(); + ReaderType::Pointer movingPolyDataReader = ReaderType::New(); movingPolyDataReader->SetFileName(argv[2]); try { movingPolyDataReader->Update(); } - catch( itk::ExceptionObject & excp ) + catch (itk::ExceptionObject & excp) { std::cerr << "Error during Moving Mesh Update() " << std::endl; std::cerr << excp << std::endl; @@ -84,37 +85,37 @@ int itkThinShellDemonsTestv4_SyN( int args, char *argv []) using MovingImageType = itk::Image; - FixedImageType::SizeType fixedImageSize; - FixedImageType::PointType fixedImageOrigin; + FixedImageType::SizeType fixedImageSize; + FixedImageType::PointType fixedImageOrigin; FixedImageType::DirectionType fixedImageDirection; - FixedImageType::SpacingType fixedImageSpacing; + FixedImageType::SpacingType fixedImageSpacing; using PointIdentifier = MeshType::PointIdentifier; using BoundingBoxType = itk::BoundingBox; BoundingBoxType::Pointer boundingBox = BoundingBoxType::New(); - PointsContainerPointer points = movingMesh->GetPoints(); + PointsContainerPointer points = movingMesh->GetPoints(); boundingBox->SetPoints(points); boundingBox->ComputeBoundingBox(); typename BoundingBoxType::PointType minBounds = boundingBox->GetMinimum(); typename BoundingBoxType::PointType maxBounds = boundingBox->GetMaximum(); - int imageDiagonal = 100; + int imageDiagonal = 100; double spacing = sqrt(boundingBox->GetDiagonalLength2()) / imageDiagonal; - auto diff = maxBounds - minBounds; - fixedImageSize[0] = ceil( 1.2 * diff[0] / spacing ); - fixedImageSize[1] = ceil( 1.2 * diff[1] / spacing ); - fixedImageSize[2] = ceil( 1.2 * diff[2] / spacing ); - fixedImageOrigin[0] = minBounds[0] - 0.1*diff[0]; - fixedImageOrigin[1] = minBounds[1] - 0.1*diff[1]; - fixedImageOrigin[2] = minBounds[2] - 0.1*diff[2]; + auto diff = maxBounds - minBounds; + fixedImageSize[0] = ceil(1.2 * diff[0] / spacing); + fixedImageSize[1] = ceil(1.2 * diff[1] / spacing); + fixedImageSize[2] = ceil(1.2 * diff[2] / spacing); + fixedImageOrigin[0] = minBounds[0] - 0.1 * diff[0]; + fixedImageOrigin[1] = minBounds[1] - 0.1 * diff[1]; + fixedImageOrigin[2] = minBounds[2] - 0.1 * diff[2]; fixedImageDirection.SetIdentity(); - fixedImageSpacing.Fill( spacing ); + fixedImageSpacing.Fill(spacing); FixedImageType::Pointer fixedImage = FixedImageType::New(); - fixedImage->SetRegions( fixedImageSize ); - fixedImage->SetOrigin( fixedImageOrigin ); - fixedImage->SetDirection( fixedImageDirection ); - fixedImage->SetSpacing( fixedImageSpacing ); + fixedImage->SetRegions(fixedImageSize); + fixedImage->SetOrigin(fixedImageOrigin); + fixedImage->SetDirection(fixedImageDirection); + fixedImage->SetSpacing(fixedImageSpacing); fixedImage->Allocate(); using VectorType = itk::Vector; @@ -136,10 +137,8 @@ int itkThinShellDemonsTestv4_SyN( int args, char *argv []) using TransformType = itk::DisplacementFieldTransform; using DisplacementFieldRegistrationType = - itk::SyNImageRegistrationMethod; - DisplacementFieldRegistrationType::Pointer registration = - DisplacementFieldRegistrationType::New(); + itk::SyNImageRegistrationMethod; + DisplacementFieldRegistrationType::Pointer registration = DisplacementFieldRegistrationType::New(); using OutputTransformType = DisplacementFieldRegistrationType::OutputTransformType; OutputTransformType::Pointer outputTransform = OutputTransformType::New(); @@ -151,11 +150,11 @@ int itkThinShellDemonsTestv4_SyN( int args, char *argv []) using AffineTransformType = itk::AffineTransform; AffineTransformType::Pointer transform = AffineTransformType::New(); transform->SetIdentity(); -/* - using PointSetMetricType = itk::EuclideanDistancePointSetToPointSetMetricv4; - PointSetMetricType::Pointer metric = PointSetMetricType::New(); -*/ - using PointSetMetricType = itk::ThinShellDemonsMetricv4 ; + /* + using PointSetMetricType = itk::EuclideanDistancePointSetToPointSetMetricv4; + PointSetMetricType::Pointer metric = PointSetMetricType::New(); + */ + using PointSetMetricType = itk::ThinShellDemonsMetricv4; PointSetMetricType::Pointer metric = PointSetMetricType::New(); metric->SetStretchWeight(1); metric->SetBendWeight(5); @@ -163,22 +162,20 @@ int itkThinShellDemonsTestv4_SyN( int args, char *argv []) metric->UseConfidenceWeightingOn(); metric->UseMaximalDistanceConfidenceSigmaOn(); metric->UpdateFeatureMatchingAtEachIterationOff(); - metric->SetMovingTransform( transform ); - //Reversed due to using points instead of an image - //to keep semantics the same as in itkThinShellDemonsTest.cxx - //For the ThinShellDemonsMetricv4 the fixed mesh is - //regularized - metric->SetFixedPointSet( movingMesh ); - metric->SetMovingPointSet( fixedMesh ); - metric->SetVirtualDomainFromImage( fixedImage ); + metric->SetMovingTransform(transform); + // Reversed due to using points instead of an image + // to keep semantics the same as in itkThinShellDemonsTest.cxx + // For the ThinShellDemonsMetricv4 the fixed mesh is + // regularized + metric->SetFixedPointSet(movingMesh); + metric->SetMovingPointSet(fixedMesh); + metric->SetVirtualDomainFromImage(fixedImage); metric->Initialize(); - double varianceForUpdateField = spacing*spacing*25; + double varianceForUpdateField = spacing * spacing * 25; double varianceForTotalField = 0.0; - registration->SetGaussianSmoothingVarianceForTheUpdateField( - varianceForUpdateField); - registration->SetGaussianSmoothingVarianceForTheTotalField( - varianceForTotalField); + registration->SetGaussianSmoothingVarianceForTheUpdateField(varianceForUpdateField); + registration->SetGaussianSmoothingVarianceForTheTotalField(varianceForTotalField); registration->SetFixedPointSet(movingMesh); registration->SetMovingPointSet(fixedMesh); registration->SetMovingInitialTransform(transform); @@ -230,14 +227,14 @@ int itkThinShellDemonsTestv4_SyN( int args, char *argv []) std::cout << "Start Value= " << metric->GetValue() << std::endl; try - { + { registration->Update(); - } - catch( itk::ExceptionObject &e ) - { + } + catch (itk::ExceptionObject & e) + { std::cerr << "Exception caught: " << e << std::endl; return EXIT_FAILURE; - } + } std::cout << "Solution Value= " << metric->GetValue() << std::endl; OutputTransformType::Pointer tx = registration->GetModifiableTransform(); @@ -250,6 +247,6 @@ int itkThinShellDemonsTestv4_SyN( int args, char *argv []) writer->SetInput(movingMesh); writer->SetFileName("synMovingMesh.vtk"); writer->Update(); - + return EXIT_SUCCESS; }