From 0af1a1991c8f179f55990fbde1f15d73f11ce6cc Mon Sep 17 00:00:00 2001 From: Bryn Lloyd Date: Sun, 6 Aug 2023 16:51:59 +0200 Subject: [PATCH] add FixTopologyCarveInside: grow from inside under topological constraints (#6) * refactor FixTopologyCarveOutside, so code can be shared * add carve from inside * add itkFixTopologyCarveInside * fix progress reporting --------- Co-authored-by: Bryn Lloyd --- .github/workflows/build-test-package.yml | 10 +- .github/workflows/clang-format-linter.yml | 1 + .../itkFixTopologyCarveOutsideTestOutput.mha | Bin 8777 -> 0 bytes include/TopologyInvariants.h | 334 +++++++++--------- include/itkFixTopologyBase.h | 157 ++++++++ include/itkFixTopologyBase.hxx | 133 +++++++ include/itkFixTopologyCarveInside.h | 101 ++++++ include/itkFixTopologyCarveInside.hxx | 175 +++++++++ include/itkFixTopologyCarveOutside.h | 89 +---- include/itkFixTopologyCarveOutside.hxx | 166 ++------- setup.py | 4 +- ...ixTopologyCarveInsideTestOutput.mha.sha512 | 1 + ...xTopologyCarveOutsideTestOutput.mha.sha512 | 1 + test/CMakeLists.txt | 14 +- test/itkFixTopologyCarveInsideTest.cxx | 147 ++++++++ test/itkFixTopologyCarveOutsideTest.cxx | 7 +- wrapping/itkFixTopologyBase.wrap | 6 + wrapping/itkFixTopologyCarveInside.wrap | 6 + 18 files changed, 962 insertions(+), 390 deletions(-) delete mode 100644 Data/Baseline/itkFixTopologyCarveOutsideTestOutput.mha create mode 100644 include/itkFixTopologyBase.h create mode 100644 include/itkFixTopologyBase.hxx create mode 100644 include/itkFixTopologyCarveInside.h create mode 100644 include/itkFixTopologyCarveInside.hxx create mode 100644 test/Baseline/itkFixTopologyCarveInsideTestOutput.mha.sha512 create mode 100644 test/Baseline/itkFixTopologyCarveOutsideTestOutput.mha.sha512 create mode 100644 test/itkFixTopologyCarveInsideTest.cxx create mode 100644 wrapping/itkFixTopologyBase.wrap create mode 100644 wrapping/itkFixTopologyCarveInside.wrap diff --git a/.github/workflows/build-test-package.yml b/.github/workflows/build-test-package.yml index 4bda130..48bc584 100644 --- a/.github/workflows/build-test-package.yml +++ b/.github/workflows/build-test-package.yml @@ -4,9 +4,13 @@ on: [push,pull_request] jobs: cxx-build-workflow: - uses: InsightSoftwareConsortium/ITKRemoteModuleBuildTestPackageAction/.github/workflows/build-test-cxx.yml@5083da2740617b78423ebf6083489e1e70ee8ca0 + uses: InsightSoftwareConsortium/ITKRemoteModuleBuildTestPackageAction/.github/workflows/build-test-cxx.yml@d4a5ce4f219b66b78269a15392e15c95f90e7e00 + with: + itk-git-tag: 'v5.3.0' python-build-workflow: - uses: InsightSoftwareConsortium/ITKRemoteModuleBuildTestPackageAction/.github/workflows/build-test-package-python.yml@5083da2740617b78423ebf6083489e1e70ee8ca0 + uses: InsightSoftwareConsortium/ITKRemoteModuleBuildTestPackageAction/.github/workflows/build-test-package-python.yml@d4a5ce4f219b66b78269a15392e15c95f90e7e00 + with: + itk-wheel-tag: 'v5.3.0' secrets: - pypi_password: ${{ secrets.PYPI_API_TOKEN }} \ No newline at end of file + pypi_password: ${{ secrets.PYPI_API_TOKEN }} diff --git a/.github/workflows/clang-format-linter.yml b/.github/workflows/clang-format-linter.yml index 899d579..59c661d 100644 --- a/.github/workflows/clang-format-linter.yml +++ b/.github/workflows/clang-format-linter.yml @@ -10,3 +10,4 @@ jobs: - uses: actions/checkout@v3 - uses: InsightSoftwareConsortium/ITKClangFormatLinterAction@master + \ No newline at end of file diff --git a/Data/Baseline/itkFixTopologyCarveOutsideTestOutput.mha b/Data/Baseline/itkFixTopologyCarveOutsideTestOutput.mha deleted file mode 100644 index 7bea79ea7ffaa7a870b6bd3a7278d93e9cbf85b5..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 8777 zcmeIwyGjE=6b9f#MPV>6APC#oSh*S{jhH(IldL3LD`Rptj_i%ui6m;_1K9WmHnw7A zt6*y(K1D>|K_?*zvA0V89hRLrXO`LjW37`jYACxKlCf))Jl-N!%?mHmt1XrlK9}7* zS6nY`*`|p~wyVgJO_G&*R(Ffsjfmv}FO(ET)SP_jZ?!Hi`e&t*+dIsb+>e|Q_Qo2tWV=5P$##AOL~s x5xBYx({nnG(fj#2Us#`>U`PT12tWV=5P$##AOHafKmY;|m|}sK{pZzm{1+7%)e8Us diff --git a/include/TopologyInvariants.h b/include/TopologyInvariants.h index a2fc01c..9564aa1 100644 --- a/include/TopologyInvariants.h +++ b/include/TopologyInvariants.h @@ -8,199 +8,187 @@ #include #include -namespace topology { +namespace topology +{ -template -inline int CountIf(const T& v) +template +inline int +CountIf(const T & v) { - static_assert(std::is_same::value, "CountIf argument must be a boolean."); - return (v != 0) ? 1 : 0; + static_assert(std::is_same::value, "CountIf argument must be a boolean."); + return (v != 0) ? 1 : 0; } /** \brief True if Euler characteristic of voxel set does not change -*/ -template -bool EulerInvariant(const TNeighborhood& neighbors, const TLabel label) + */ +template +bool +EulerInvariant(const TNeighborhood & neighbors, const TLabel label) { - // For a face-connected set, a vertex, edge, or facet is included - // only if all voxels it is incident to are in the set, and for a - // vertex-connected set it is included if any incident voxels are in the set - - // Here we count only the changes at the center voxel. - // ==> If we mark the center as BG, then P, F, E and V are '0'. - static const int c = 27 / 2; - - assert(neighbors[c] == label); - - // Voxels (parallelepipeds) - count center - const int P = 1; - - // Faces - count faces around center voxel - const int F = CountIf(neighbors[c - 1] == label) // -x - + CountIf(neighbors[c + 1] == label) // +x - + CountIf(neighbors[c - 3] == label) // -y - + CountIf(neighbors[c + 3] == label) // +y - + CountIf(neighbors[c - 9] == label) // -z - + CountIf(neighbors[c + 9] == label); // +z - - // Edges - count edges around center voxel - int E = 0; - { - static const short offsets[12][2] = { - {-1, -3}, {+1, -3}, {-1, +3}, {+1, +3}, {-1, -9}, {+1, -9}, {-1, +9}, {+1, +9}, {-3, -9}, {+3, -9}, {-3, +9}, {+3, +9}}; - for (int i = 0; i < 12; i++) - { - short o = offsets[i][0]; - short u = offsets[i][1]; - - E += CountIf(neighbors[c + o] == label && - neighbors[c + u] == label && - neighbors[c + o + u] == label); - } - } - - // Vertices - count verts around center voxel - int V = 0; - { - static const short offsets[8][3] = { - {-1, -3, -9}, {+1, -3, -9}, {-1, +3, -9}, {+1, +3, -9}, {-1, -3, +9}, {+1, -3, +9}, {-1, +3, +9}, {+1, +3, +9}}; - for (int i = 0; i < 8; i++) - { - short x = offsets[i][0]; - short y = offsets[i][1]; - short z = offsets[i][2]; - - V += CountIf(neighbors[c + x] == label && - neighbors[c + y] == label && - neighbors[c + z] == label && - neighbors[c + x + y] == label && - neighbors[c + x + z] == label && - neighbors[c + y + z] == label && - neighbors[c + x + y + z] == label); - } - } - - return (V - E + F - P) == 0; + // For a face-connected set, a vertex, edge, or facet is included + // only if all voxels it is incident to are in the set, and for a + // vertex-connected set it is included if any incident voxels are in the set + + // Here we count only the changes at the center voxel. + // ==> If we mark the center as BG, then P, F, E and V are '0'. + static constexpr int c = 27 / 2; + + assert(neighbors[c] == label); + + // Voxels (parallelepipeds) - count center + const int P = 1; + + // Faces - count faces around center voxel + const int F = CountIf(neighbors[c - 1] == label) // -x + + CountIf(neighbors[c + 1] == label) // +x + + CountIf(neighbors[c - 3] == label) // -y + + CountIf(neighbors[c + 3] == label) // +y + + CountIf(neighbors[c - 9] == label) // -z + + CountIf(neighbors[c + 9] == label); // +z + + // Edges - count edges around center voxel + int E = 0; + { + static const short offsets[12][2] = { { -1, -3 }, { +1, -3 }, { -1, +3 }, { +1, +3 }, { -1, -9 }, { +1, -9 }, + { -1, +9 }, { +1, +9 }, { -3, -9 }, { +3, -9 }, { -3, +9 }, { +3, +9 } }; + for (int i = 0; i < 12; i++) + { + short o = offsets[i][0]; + short u = offsets[i][1]; + + E += CountIf(neighbors[c + o] == label && neighbors[c + u] == label && neighbors[c + o + u] == label); + } + } + + // Vertices - count verts around center voxel + int V = 0; + { + static const short offsets[8][3] = { { -1, -3, -9 }, { +1, -3, -9 }, { -1, +3, -9 }, { +1, +3, -9 }, + { -1, -3, +9 }, { +1, -3, +9 }, { -1, +3, +9 }, { +1, +3, +9 } }; + for (int i = 0; i < 8; i++) + { + short x = offsets[i][0]; + short y = offsets[i][1]; + short z = offsets[i][2]; + + V += CountIf(neighbors[c + x] == label && neighbors[c + y] == label && neighbors[c + z] == label && + neighbors[c + x + y] == label && neighbors[c + x + z] == label && neighbors[c + y + z] == label && + neighbors[c + x + y + z] == label); + } + } + + return (V - E + F - P) == 0; } /** \brief Returns the number of connected components for selected label -*/ -template -unsigned ConnectedComponents(const TNeighborhood& neighbors, const TLabel label) + */ +template +unsigned +ConnectedComponents(const TNeighborhood & neighbors, const TLabel label) { - std::deque queue; - std::array visited; - visited.fill(false); - - static const short ijk_lut[27][3] = { - {0, 0, 0}, - {1, 0, 0}, - {2, 0, 0}, - {0, 1, 0}, - {1, 1, 0}, - {2, 1, 0}, - {0, 2, 0}, - {1, 2, 0}, - {2, 2, 0}, // - {0, 0, 1}, - {1, 0, 1}, - {2, 0, 1}, - {0, 1, 1}, - {1, 1, 1}, - {2, 1, 1}, - {0, 2, 1}, - {1, 2, 1}, - {2, 2, 1}, // - {0, 0, 2}, - {1, 0, 2}, - {2, 0, 2}, - {0, 1, 2}, - {1, 1, 2}, - {2, 1, 2}, - {0, 2, 2}, - {1, 2, 2}, - {2, 2, 2}, - }; - - unsigned new_label = 0; - for (unsigned n = 0; n < 27; ++n) - { - // skip other labels and visited voxels - if (neighbors[n] != label || visited[n]) - continue; - - queue.push_back(n); - visited[n] = true; - new_label++; - - while (!queue.empty()) - { - auto id = queue.front(); - queue.pop_front(); - auto ijk = ijk_lut[id]; - - if (ijk[0] < 2 && neighbors[id] == neighbors[id + 1] && !visited[id + 1]) // +x - { - queue.push_back(id + 1); - visited[id + 1] = true; - } - if (ijk[0] > 0 && neighbors[id] == neighbors[id - 1] && !visited[id - 1]) // -x - { - queue.push_back(id - 1); - visited[id - 1] = true; - } - if (ijk[1] < 2 && neighbors[id] == neighbors[id + 3] && !visited[id + 3]) // +y - { - queue.push_back(id + 3); - visited[id + 3] = true; - } - if (ijk[1] > 0 && neighbors[id] == neighbors[id - 3] && !visited[id - 3]) // -y - { - queue.push_back(id - 3); - visited[id - 3] = true; - } - if (ijk[2] < 2 && neighbors[id] == neighbors[id + 9] && !visited[id + 9]) // +z - { - queue.push_back(id + 9); - visited[id + 9] = true; - } - if (ijk[2] > 0 && neighbors[id] == neighbors[id - 9] && !visited[id - 9]) // -z - { - queue.push_back(id - 9); - visited[id - 9] = true; - } - } - } - - return new_label; + std::deque queue; + std::array visited; + visited.fill(false); + + static const short ijk_lut[27][3] = { + { 0, 0, 0 }, { 1, 0, 0 }, { 2, 0, 0 }, { 0, 1, 0 }, { 1, 1, 0 }, + { 2, 1, 0 }, { 0, 2, 0 }, { 1, 2, 0 }, { 2, 2, 0 }, // + { 0, 0, 1 }, { 1, 0, 1 }, { 2, 0, 1 }, { 0, 1, 1 }, { 1, 1, 1 }, + { 2, 1, 1 }, { 0, 2, 1 }, { 1, 2, 1 }, { 2, 2, 1 }, // + { 0, 0, 2 }, { 1, 0, 2 }, { 2, 0, 2 }, { 0, 1, 2 }, { 1, 1, 2 }, + { 2, 1, 2 }, { 0, 2, 2 }, { 1, 2, 2 }, { 2, 2, 2 }, + }; + + unsigned new_label = 0; + for (unsigned n = 0; n < 27; ++n) + { + // skip other labels and visited voxels + if (neighbors[n] != label || visited[n]) + continue; + + queue.push_back(n); + visited[n] = true; + new_label++; + + while (!queue.empty()) + { + auto id = queue.front(); + queue.pop_front(); + auto ijk = ijk_lut[id]; + + if (ijk[0] < 2 && neighbors[id] == neighbors[id + 1] && !visited[id + 1]) // +x + { + queue.push_back(id + 1); + visited[id + 1] = true; + } + if (ijk[0] > 0 && neighbors[id] == neighbors[id - 1] && !visited[id - 1]) // -x + { + queue.push_back(id - 1); + visited[id - 1] = true; + } + if (ijk[1] < 2 && neighbors[id] == neighbors[id + 3] && !visited[id + 3]) // +y + { + queue.push_back(id + 3); + visited[id + 3] = true; + } + if (ijk[1] > 0 && neighbors[id] == neighbors[id - 3] && !visited[id - 3]) // -y + { + queue.push_back(id - 3); + visited[id - 3] = true; + } + if (ijk[2] < 2 && neighbors[id] == neighbors[id + 9] && !visited[id + 9]) // +z + { + queue.push_back(id + 9); + visited[id + 9] = true; + } + if (ijk[2] > 0 && neighbors[id] == neighbors[id - 9] && !visited[id - 9]) // -z + { + queue.push_back(id - 9); + visited[id - 9] = true; + } + } + } + + return new_label; } /** \brief True if number of connected components of voxel set does not change -*/ -template -bool CCInvariant(TNeighborhood neighbors, const TLabel label) + */ +template +bool +CCInvariant(TNeighborhood neighbors, const TLabel label) { - neighbors[27 / 2] = 1; - unsigned cc_before = ConnectedComponents(neighbors, label); + neighbors[27 / 2] = 1; + unsigned cc_before = ConnectedComponents(neighbors, label); - neighbors[27 / 2] = 0; - unsigned cc_after = ConnectedComponents(neighbors, label); - return (cc_before == cc_after); + neighbors[27 / 2] = 0; + unsigned cc_after = ConnectedComponents(neighbors, label); + return (cc_before == cc_after); } /** \brief Will voxels be manifold after removing center voxel? */ -template -bool NonmanifoldRemove(const TNeighborhood& neighbors, const TLabel label) +template +bool +NonmanifoldRemove(const TNeighborhood & neighbors, const TLabel label) { - static const int c = 27 / 2; - - // test for non-manifold edges around center voxel - return (neighbors[c - 1] == label && neighbors[c - 3] == label && neighbors[c - 1 - 3] != label) || (neighbors[c + 1] == label && neighbors[c - 3] == label && neighbors[c + 1 - 3] != label) || (neighbors[c - 1] == label && neighbors[c + 3] == label && neighbors[c - 1 + 3] != label) || (neighbors[c + 1] == label && neighbors[c + 3] == label && neighbors[c + 1 + 3] != label) - // - || (neighbors[c - 1] == label && neighbors[c - 9] == label && neighbors[c - 1 - 9] != label) || (neighbors[c + 1] == label && neighbors[c - 9] == label && neighbors[c + 1 - 9] != label) || (neighbors[c - 1] == label && neighbors[c + 9] == label && neighbors[c - 1 + 9] != label) || (neighbors[c + 1] == label && neighbors[c + 9] == label && neighbors[c + 1 + 9] != label) - // - || (neighbors[c - 3] == label && neighbors[c - 9] == label && neighbors[c - 3 - 9] != label) || (neighbors[c + 3] == label && neighbors[c - 9] == label && neighbors[c + 3 - 9] != label) || (neighbors[c - 3] == label && neighbors[c + 9] == label && neighbors[c - 3 + 9] != label) || (neighbors[c + 3] == label && neighbors[c + 9] == label && neighbors[c + 3 + 9] != label); + static const int c = 27 / 2; + + // test for non-manifold edges around center voxel + return (neighbors[c - 1] == label && neighbors[c - 3] == label && neighbors[c - 1 - 3] != label) || + (neighbors[c + 1] == label && neighbors[c - 3] == label && neighbors[c + 1 - 3] != label) || + (neighbors[c - 1] == label && neighbors[c + 3] == label && neighbors[c - 1 + 3] != label) || + (neighbors[c + 1] == label && neighbors[c + 3] == label && neighbors[c + 1 + 3] != label) + // + || (neighbors[c - 1] == label && neighbors[c - 9] == label && neighbors[c - 1 - 9] != label) || + (neighbors[c + 1] == label && neighbors[c - 9] == label && neighbors[c + 1 - 9] != label) || + (neighbors[c - 1] == label && neighbors[c + 9] == label && neighbors[c - 1 + 9] != label) || + (neighbors[c + 1] == label && neighbors[c + 9] == label && neighbors[c + 1 + 9] != label) + // + || (neighbors[c - 3] == label && neighbors[c - 9] == label && neighbors[c - 3 - 9] != label) || + (neighbors[c + 3] == label && neighbors[c - 9] == label && neighbors[c + 3 - 9] != label) || + (neighbors[c - 3] == label && neighbors[c + 9] == label && neighbors[c - 3 + 9] != label) || + (neighbors[c + 3] == label && neighbors[c + 9] == label && neighbors[c + 3 + 9] != label); } } // namespace topology diff --git a/include/itkFixTopologyBase.h b/include/itkFixTopologyBase.h new file mode 100644 index 0000000..7c540da --- /dev/null +++ b/include/itkFixTopologyBase.h @@ -0,0 +1,157 @@ +/*========================================================================= + * + * 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 itkFixTopologyBase_h +#define itkFixTopologyBase_h + +#include "itkImageToImageFilter.h" +#include "itkProgressAccumulator.h" + +#include + +namespace itk +{ +/** \class FixTopologyBase + * + * \brief Base class for morphological closing/opening with topology constraints + * + * \author Bryn Lloyd + * \ingroup TopologyControl + */ +template +class ITK_TEMPLATE_EXPORT FixTopologyBase : public ImageToImageFilter +{ +public: + ITK_DISALLOW_COPY_AND_MOVE(FixTopologyBase); + + /** Extract dimension from input and output image. */ + itkStaticConstMacro(InputImageDimension, unsigned int, TInputImage::ImageDimension); + itkStaticConstMacro(OutputImageDimension, unsigned int, TOutputImage::ImageDimension); + + static const unsigned int ImageDimension = InputImageDimension; + + /** Standard class typedefs. */ + using Self = FixTopologyBase; + using Superclass = ImageToImageFilter; + using Pointer = SmartPointer; + using ConstPointer = SmartPointer; + + /** Run-time type information (and related methods). */ + itkTypeMacro(FixTopologyBase, ImageToImageFilter); + + /** Type for input image. */ + using InputImageType = TInputImage; + + /** Type for output image: Skeleton of the object. */ + using OutputImageType = TOutputImage; + + /** Type for mask image */ + using MaskImageType = TMaskImage; + + /** Type for the pixel type of the input image. */ + using InputImagePixelType = typename InputImageType::PixelType; + + /** Type for the pixel type of the input image. */ + using OutputImagePixelType = typename OutputImageType::PixelType; + + /** Pointer Type for input image. */ + using InputImagePointer = typename InputImageType::ConstPointer; + + /** Pointer Type for the output image. */ + using OutputImagePointer = typename OutputImageType::Pointer; + + /** Pointer Type for the mask image. */ + using MaskImageTypePointer = typename MaskImageType::Pointer; + + /** Optional mask (if none is provided the input mask is dilated by 'Radius') */ + void + SetMaskImage(const MaskImageType * mask); + const MaskImageType * + GetMaskImage() const; + + itkSetMacro(Radius, SizeValueType); + itkGetConstMacro(Radius, SizeValueType); + + itkSetMacro(InsideValue, InputImagePixelType); + itkGetConstMacro(InsideValue, InputImagePixelType); + +#ifdef ITK_USE_CONCEPT_CHECKING + /** Begin concept checking */ + itkConceptMacro(SameDimensionCheck, (Concept::SameDimension)); + itkConceptMacro(SameTypeCheck, (Concept::SameType)); + itkConceptMacro(InputAdditiveOperatorsCheck, (Concept::AdditiveOperators)); + itkConceptMacro(InputConvertibleToIntCheck, (Concept::Convertible)); + itkConceptMacro(IntConvertibleToInputCheck, (Concept::Convertible)); + itkConceptMacro(InputIntComparableCheck, (Concept::Comparable)); + /** End concept checking */ +#endif + +protected: + FixTopologyBase(); + ~FixTopologyBase() override = default; + void + PrintSelf(std::ostream & os, Indent indent) const override; + + void + GenerateData() override; + + void + PrepareData(ProgressAccumulator * progress); + + virtual MaskImageType * + CreateDefaultMask(ProgressAccumulator * progress) = 0; + + virtual void + ComputeThinImage(ProgressAccumulator * progress) = 0; + + std::vector + GetNeighborOffsets() + { + // 18-connectivity + using OffsetType = typename InputImageType::OffsetType; + return { OffsetType{ -1, 0, 0 }, OffsetType{ 1, 0, 0 }, OffsetType{ 0, -1, 0 }, OffsetType{ 0, 1, 0 }, + OffsetType{ 0, 0, -1 }, OffsetType{ 0, 0, 1 }, OffsetType{ -1, -1, 0 }, OffsetType{ 1, -1, 0 }, + OffsetType{ -1, 1, 0 }, OffsetType{ 1, 1, 0 }, OffsetType{ 0, -1, -1 }, OffsetType{ 0, 1, -1 }, + OffsetType{ 0, -1, 1 }, OffsetType{ 0, 1, 1 }, OffsetType{ -1, 0, -1 }, OffsetType{ 1, 0, -1 }, + OffsetType{ -1, 0, 1 }, OffsetType{ 1, 0, 1 } }; + } + + enum ePixelState : OutputImagePixelType + { + kBackground = 0, + kHardForeground, + kSoftForeground, + kQueued + }; + + MaskImageTypePointer m_PaddedOutput; + + using RealImageType = itk::Image; + RealImageType::Pointer m_DistanceMap; + + SizeValueType m_Radius = 1; + InputImagePixelType m_InsideValue = 1; +}; // end of FixTopologyBase class + +} // end namespace itk + +#ifndef ITK_MANUAL_INSTANTIATION +# include "itkFixTopologyBase.hxx" +#endif + +#endif // itkFixTopologyBase diff --git a/include/itkFixTopologyBase.hxx b/include/itkFixTopologyBase.hxx new file mode 100644 index 0000000..c7647a3 --- /dev/null +++ b/include/itkFixTopologyBase.hxx @@ -0,0 +1,133 @@ +/*========================================================================= + * + * 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 itkFixTopologyBase_hxx +#define itkFixTopologyBase_hxx + +#include "itkFixTopologyBase.h" + +#include "itkImageRegionConstIterator.h" +#include "itkImageRegionRange.h" +#include "itkSignedMaurerDistanceMapImageFilter.h" + +namespace itk +{ + +template +FixTopologyBase::FixTopologyBase() +{ + this->SetNumberOfRequiredOutputs(1); +} + +template +void +FixTopologyBase::SetMaskImage(const TMaskImage * mask) +{ + this->ProcessObject::SetNthInput(1, const_cast(mask)); +} + +template +const TMaskImage * +FixTopologyBase::GetMaskImage() const +{ + return itkDynamicCastInDebugMode(const_cast(this->ProcessObject::GetInput(1))); +} + +template +void +FixTopologyBase::PrepareData(ProgressAccumulator * progress) +{ + InputImagePointer input_image = dynamic_cast(ProcessObject::GetInput(0)); + OutputImagePointer thin_image = this->GetOutput(); + thin_image->SetBufferedRegion(thin_image->GetRequestedRegion()); + thin_image->Allocate(); + + // pad by 1 layer so we can get 1x1x1 neighborhood without checking if we are at boundary + auto region = thin_image->GetRequestedRegion(); + auto padded_region = region; + padded_region.PadByRadius(1); + + m_PaddedOutput = MaskImageType::New(); + m_PaddedOutput->SetRegions(padded_region); + m_PaddedOutput->SetSpacing(thin_image->GetSpacing()); + m_PaddedOutput->Allocate(); + m_PaddedOutput->FillBuffer(ePixelState::kBackground); + + ImageRegionConstIterator it(input_image, region); + ImageRegionIterator ot(m_PaddedOutput, region); + + // mark the input mask as kHardForeground + for (it.GoToBegin(), ot.GoToBegin(); !ot.IsAtEnd(); ++it, ++ot) + { + if (it.Get() == m_InsideValue) + { + ot.Set(ePixelState::kHardForeground); + } + } + + // compute distance map: used for priority queue + auto distance_filter = SignedMaurerDistanceMapImageFilter::New(); + progress->RegisterInternalFilter(distance_filter, 0.1); + distance_filter->SetInput(m_PaddedOutput); + distance_filter->SetUseImageSpacing(true); + distance_filter->SetInsideIsPositive(false); + distance_filter->SetSquaredDistance(false); + distance_filter->SetBackgroundValue(0); + distance_filter->Update(); + m_DistanceMap = distance_filter->GetOutput(); + + typename MaskImageType::ConstPointer mask_image = GetMaskImage(); + if (!mask_image) + { + // if no mask is provided we dilate the input mask + mask_image = this->CreateDefaultMask(progress); + } + + // mark dilated as '2', but don't copy padding layer + ImageRegionConstIterator mt(mask_image, region); + for (mt.GoToBegin(), ot.GoToBegin(); !ot.IsAtEnd(); ++mt, ++ot) + { + if (mt.Get() != ot.Get()) + { + ot.Set(ePixelState::kSoftForeground); + } + } +} + +template +void +FixTopologyBase::GenerateData() +{ + ProgressAccumulator::Pointer progress = ProgressAccumulator::New(); + progress->SetMiniPipelineFilter(this); + + this->PrepareData(progress); + + this->ComputeThinImage(progress); +} + +template +void +FixTopologyBase::PrintSelf(std::ostream & os, Indent indent) const +{ + Superclass::PrintSelf(os, indent); +} + +} // end namespace itk + +#endif // itkFixTopologyBase_hxx diff --git a/include/itkFixTopologyCarveInside.h b/include/itkFixTopologyCarveInside.h new file mode 100644 index 0000000..a9378dc --- /dev/null +++ b/include/itkFixTopologyCarveInside.h @@ -0,0 +1,101 @@ +/*========================================================================= + * + * 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 itkFixTopologyCarveInside_h +#define itkFixTopologyCarveInside_h + +#include "itkFixTopologyBase.h" + +namespace itk +{ +/** \class FixTopologyCarveInside + * + * \brief This filter does morphological opening with topology constraints + * + * It works by doing following steps: + * 1. ... + * 2. ... + * + * The first step closes holes and the second returns as close as possible to the input mask, while ensuring that the + * holes are not re-opened. + * + * This filter implements ideas from: Vanderhyde, James. "Topology control of volumetric data.", PhD dissertation, + * Georgia Institute of Technology, 2007.. + * + * \author Bryn Lloyd + * \ingroup TopologyControl + */ +template > +class ITK_TEMPLATE_EXPORT FixTopologyCarveInside : public FixTopologyBase +{ +public: + ITK_DISALLOW_COPY_AND_MOVE(FixTopologyCarveInside); + + /** Extract dimension from input and output image. */ + itkStaticConstMacro(InputImageDimension, unsigned int, TInputImage::ImageDimension); + itkStaticConstMacro(OutputImageDimension, unsigned int, TOutputImage::ImageDimension); + + static const unsigned int ImageDimension = InputImageDimension; + + /** Standard class typedefs. */ + using Self = FixTopologyCarveInside; + using Superclass = FixTopologyBase; + using Pointer = SmartPointer; + using ConstPointer = SmartPointer; + + /** Method for creation through the object factory */ + itkNewMacro(Self); + + /** Run-time type information (and related methods). */ + itkTypeMacro(FixTopologyCarveInside, FixTopologyBase); + + /** Type for input image. */ + using InputImageType = TInputImage; + + /** Type for output image: Skeleton of the object. */ + using OutputImageType = TOutputImage; + + /** Type for mask image */ + using MaskImageType = TMaskImage; + + /** Pointer Type for input image. */ + using InputImagePointer = typename InputImageType::ConstPointer; + + /** Pointer Type for the output image. */ + using OutputImagePointer = typename OutputImageType::Pointer; + +protected: + FixTopologyCarveInside() = default; + ~FixTopologyCarveInside() override = default; + + MaskImageType * + CreateDefaultMask(ProgressAccumulator * progress) override; + + void + ComputeThinImage(ProgressAccumulator * progress) override; + + using ePixelState = typename Superclass::ePixelState; +}; // end of FixTopologyCarveInside class + +} // end namespace itk + +#ifndef ITK_MANUAL_INSTANTIATION +# include "itkFixTopologyCarveInside.hxx" +#endif + +#endif // itkFixTopologyCarveInside diff --git a/include/itkFixTopologyCarveInside.hxx b/include/itkFixTopologyCarveInside.hxx new file mode 100644 index 0000000..4028323 --- /dev/null +++ b/include/itkFixTopologyCarveInside.hxx @@ -0,0 +1,175 @@ +/*========================================================================= + * + * 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 itkFixTopologyCarveInside_hxx +#define itkFixTopologyCarveInside_hxx + +#include "itkFixTopologyCarveInside.h" +#include "TopologyInvariants.h" + +#include "itkBinaryErodeImageFilter.h" +#include "itkFlatStructuringElement.h" +#include "itkNeighborhoodIterator.h" + +namespace itk +{ + +template +TMaskImage * +FixTopologyCarveInside::CreateDefaultMask(ProgressAccumulator * progress) +{ + // if no mask is provided we dilate the input mask + using kernel_type = itk::FlatStructuringElement<3>; + auto radius = this->GetRadius(); + auto ball = kernel_type::Ball({ radius, radius, radius }, false); + + auto erode = itk::BinaryErodeImageFilter::New(); + progress->RegisterInternalFilter(erode, 0.1); + erode->SetInput(this->m_PaddedOutput); + erode->SetKernel(ball); + erode->SetForegroundValue(ePixelState::kHardForeground); + erode->SetBackgroundValue(0); + erode->Update(); + return erode->GetOutput(); +} + +template +void +FixTopologyCarveInside::ComputeThinImage(ProgressAccumulator * accumulator) +{ + using IndexType = typename InputImageType::IndexType; + using NeighborhoodIteratorType = NeighborhoodIterator>; + + OutputImagePointer thin_image = this->GetOutput(); + + // note: all processing is done on padded_output to avoid + // checking if index is in region + auto padded_output = this->m_PaddedOutput; + auto region = thin_image->GetRequestedRegion(); + + // for progress reporting + size_t mask_size = 0; + itk::ImageRegionRange image_range(*padded_output, region); + for (auto && pixel : image_range) + { + mask_size += (pixel == ePixelState::kSoftForeground) ? 1 : 0; + } + + // process pixels further away from input background first + // - distance is negative inside + // - use min priority queue + using node = std::pair; + const auto cmp = [](const node & l, const node & r) { return l.first > r.first; }; + std::priority_queue, decltype(cmp)> queue(cmp); + + auto neighbors = this->GetNeighborOffsets(); + const size_t num_neighbors = neighbors.size(); + + auto add_neighbors = [&](const IndexType & idx) { + for (size_t k = 0; k < num_neighbors; ++k) + { + const IndexType n_id = idx + neighbors[k]; + const float n_dist = this->m_DistanceMap->GetPixel(n_id); + + if (padded_output->GetPixel(n_id) == ePixelState::kSoftForeground) + { + // mark as visited + padded_output->SetPixel(n_id, ePixelState::kQueued); + + // add to queue + queue.push(std::make_pair(n_dist, n_id)); + } + } + }; + + // initial seeds + itk::ImageRegionConstIterator it(padded_output, region); + for (it.GoToBegin(); !it.IsAtEnd(); ++it) + { + if (it.Get() == ePixelState::kHardForeground) + { + add_neighbors(it.GetIndex()); + } + } + + NeighborhoodIteratorType n_it({ 1, 1, 1 }, padded_output, region); + + const auto get_mask = [&n_it](const IndexType & idx) { + n_it.SetLocation(idx); + auto n = n_it.GetNeighborhood(); + for (auto & v : n.GetBufferReference()) + v = (v == ePixelState::kHardForeground) ? 1 : 0; + return n; + }; + + // dilate while topology does not change + ProgressReporter progress(this, 0, mask_size, 100); + while (true) + { + int num_changed = 0; + + while (!queue.empty()) + { + auto idx = queue.top().second; // node + queue.pop(); + + // skip if already processed + if (padded_output->GetPixel(idx) != ePixelState::kQueued) + continue; + + auto vals = get_mask(idx); + + // check if point is simple (deletion does not change connectivity in the 3x3x3 neighborhood) + if (topology::EulerInvariant(vals, 0) && topology::CCInvariant(vals, 0)) + { + padded_output->SetPixel(idx, ePixelState::kHardForeground); + progress.CompletedPixel(); + num_changed++; + } + + // add unvisited neighbors to queue + add_neighbors(idx); + } + + if (num_changed == 0) + break; + + for (it.GoToBegin(); !it.IsAtEnd(); ++it) + { + if (it.Get() == ePixelState::kQueued) + { + queue.push(std::make_pair(this->m_DistanceMap->GetPixel(it.GetIndex()), it.GetIndex())); + } + } + } + + // copy to output + InputImagePointer input_image = dynamic_cast(ProcessObject::GetInput(0)); + itk::ImageRegionConstIterator i_it(input_image, region); + itk::ImageRegionConstIterator s_it(padded_output, region); + itk::ImageRegionIterator o_it(thin_image, region); + + for (i_it.GoToBegin(), s_it.GoToBegin(), o_it.GoToBegin(); !s_it.IsAtEnd(); ++i_it, ++s_it, ++o_it) + { + o_it.Set(s_it.Get() == ePixelState::kHardForeground ? this->m_InsideValue : 0); + } +} + +} // end namespace itk + +#endif // itkFixTopologyCarveInside_hxx diff --git a/include/itkFixTopologyCarveOutside.h b/include/itkFixTopologyCarveOutside.h index a5b628b..de65b59 100644 --- a/include/itkFixTopologyCarveOutside.h +++ b/include/itkFixTopologyCarveOutside.h @@ -19,10 +19,7 @@ #ifndef itkFixTopologyCarveOutside_h #define itkFixTopologyCarveOutside_h -#include "itkImageToImageFilter.h" -#include "itkProgressAccumulator.h" - -#include +#include "itkFixTopologyBase.h" namespace itk { @@ -44,9 +41,11 @@ namespace itk * \ingroup TopologyControl */ template > -class ITK_TEMPLATE_EXPORT FixTopologyCarveOutside : public ImageToImageFilter +class ITK_TEMPLATE_EXPORT FixTopologyCarveOutside : public FixTopologyBase { public: + ITK_DISALLOW_COPY_AND_MOVE(FixTopologyCarveOutside); + /** Extract dimension from input and output image. */ itkStaticConstMacro(InputImageDimension, unsigned int, TInputImage::ImageDimension); itkStaticConstMacro(OutputImageDimension, unsigned int, TOutputImage::ImageDimension); @@ -55,7 +54,7 @@ class ITK_TEMPLATE_EXPORT FixTopologyCarveOutside : public ImageToImageFilter; + using Superclass = FixTopologyBase; using Pointer = SmartPointer; using ConstPointer = SmartPointer; @@ -63,7 +62,7 @@ class ITK_TEMPLATE_EXPORT FixTopologyCarveOutside : public ImageToImageFilter)); - itkConceptMacro(SameTypeCheck, (Concept::SameType)); - itkConceptMacro(InputAdditiveOperatorsCheck, (Concept::AdditiveOperators)); - itkConceptMacro(InputConvertibleToIntCheck, (Concept::Convertible)); - itkConceptMacro(IntConvertibleToInputCheck, (Concept::Convertible)); - itkConceptMacro(InputIntComparableCheck, (Concept::Comparable)); - /** End concept checking */ -#endif - protected: - FixTopologyCarveOutside(); + FixTopologyCarveOutside() = default; ~FixTopologyCarveOutside() override = default; - void - PrintSelf(std::ostream & os, Indent indent) const override; - void - GenerateData() override; + MaskImageType * + CreateDefaultMask(ProgressAccumulator * progress) override; void - PrepareData(ProgressAccumulator * progress); + ComputeThinImage(ProgressAccumulator * progress) override; - void - ComputeThinImage(ProgressAccumulator * progress); - - std::vector - GetNeighborOffsets() - { - // 18-connectivity - using OffsetType = typename InputImageType::OffsetType; - return { OffsetType{ -1, 0, 0 }, OffsetType{ 1, 0, 0 }, OffsetType{ 0, -1, 0 }, OffsetType{ 0, 1, 0 }, - OffsetType{ 0, 0, -1 }, OffsetType{ 0, 0, 1 }, OffsetType{ -1, -1, 0 }, OffsetType{ 1, -1, 0 }, - OffsetType{ -1, 1, 0 }, OffsetType{ 1, 1, 0 }, OffsetType{ 0, -1, -1 }, OffsetType{ 0, 1, -1 }, - OffsetType{ 0, -1, 1 }, OffsetType{ 0, 1, 1 }, OffsetType{ -1, 0, -1 }, OffsetType{ 1, 0, -1 }, - OffsetType{ -1, 0, 1 }, OffsetType{ 1, 0, 1 } }; - } - -private: - ITK_DISALLOW_COPY_AND_ASSIGN(FixTopologyCarveOutside); - - enum ePixelState : OutputImagePixelType - { - kBackground = 0, - kHardForeground, - kDilated, - kVisited - }; - - MaskImageTypePointer m_PaddedOutput; - - using RealImageType = itk::Image; - RealImageType::Pointer m_DistanceMap; - - SizeValueType m_Radius = 1; - InputImagePixelType m_InsideValue = 1; + using ePixelState = typename Superclass::ePixelState; }; // end of FixTopologyCarveOutside class } // end namespace itk diff --git a/include/itkFixTopologyCarveOutside.hxx b/include/itkFixTopologyCarveOutside.hxx index 64c84a2..87be310 100644 --- a/include/itkFixTopologyCarveOutside.hxx +++ b/include/itkFixTopologyCarveOutside.hxx @@ -23,105 +23,29 @@ #include "TopologyInvariants.h" #include "itkBinaryDilateImageFilter.h" -#include "itkConstantBoundaryCondition.h" #include "itkFlatStructuringElement.h" #include "itkImageLinearConstIteratorWithIndex.h" -#include "itkImageRegionConstIterator.h" -#include "itkImageRegionRange.h" #include "itkNeighborhoodIterator.h" -#include "itkSignedMaurerDistanceMapImageFilter.h" namespace itk { template -FixTopologyCarveOutside::FixTopologyCarveOutside() +TMaskImage * +FixTopologyCarveOutside::CreateDefaultMask(ProgressAccumulator * progress) { - this->SetNumberOfRequiredOutputs(1); -} - -template -void -FixTopologyCarveOutside::SetMaskImage(const TMaskImage * mask) -{ - this->ProcessObject::SetNthInput(1, const_cast(mask)); -} - -template -const TMaskImage * -FixTopologyCarveOutside::GetMaskImage() const -{ - return itkDynamicCastInDebugMode(const_cast(this->ProcessObject::GetInput(1))); -} - -template -void -FixTopologyCarveOutside::PrepareData(ProgressAccumulator * progress) -{ - InputImagePointer input_image = dynamic_cast(ProcessObject::GetInput(0)); - OutputImagePointer thin_image = this->GetOutput(); - thin_image->SetBufferedRegion(thin_image->GetRequestedRegion()); - thin_image->Allocate(); - - // pad by 1 layer so we can get 1x1x1 neighborhood without checking if we are at boundary - auto region = thin_image->GetRequestedRegion(); - auto padded_region = region; - padded_region.PadByRadius(1); - - m_PaddedOutput = MaskImageType::New(); - m_PaddedOutput->SetRegions(padded_region); - m_PaddedOutput->SetSpacing(thin_image->GetSpacing()); - m_PaddedOutput->Allocate(); - m_PaddedOutput->FillBuffer(ePixelState::kBackground); - - ImageRegionConstIterator it(input_image, region); - ImageRegionIterator ot(m_PaddedOutput, region); - - // mark the input mask as kHardForeground - for (it.GoToBegin(), ot.GoToBegin(); !ot.IsAtEnd(); ++it, ++ot) - { - if (it.Get() == m_InsideValue) - { - ot.Set(ePixelState::kHardForeground); - } - } - - // compute distance map: used for priority queue - auto distance_filter = SignedMaurerDistanceMapImageFilter::New(); - progress->RegisterInternalFilter(distance_filter, 0.1); - distance_filter->SetInput(m_PaddedOutput); - distance_filter->SetUseImageSpacing(true); - distance_filter->SetInsideIsPositive(false); - distance_filter->SetSquaredDistance(false); - distance_filter->SetBackgroundValue(0); - distance_filter->Update(); - m_DistanceMap = distance_filter->GetOutput(); - - typename MaskImageType::ConstPointer mask_image = GetMaskImage(); - if (!mask_image) - { - // if no mask is provided we dilate the input mask - using kernel_type = itk::FlatStructuringElement<3>; - auto ball = kernel_type::Ball({ m_Radius, m_Radius, m_Radius }, false); - - auto dilate = itk::BinaryDilateImageFilter::New(); - progress->RegisterInternalFilter(dilate, 0.1); - dilate->SetInput(m_PaddedOutput); - dilate->SetKernel(ball); - dilate->SetForegroundValue(ePixelState::kHardForeground); - dilate->Update(); - mask_image = dilate->GetOutput(); - } - - // mark dilated as '2', but don't copy padding layer - ImageRegionConstIterator mt(mask_image, region); - for (mt.GoToBegin(), ot.GoToBegin(); !ot.IsAtEnd(); ++mt, ++ot) - { - if (mt.Get() != ot.Get()) - { - ot.Set(ePixelState::kDilated); - } - } + // if no mask is provided we dilate the input mask + using kernel_type = itk::FlatStructuringElement<3>; + auto radius = this->GetRadius(); + auto ball = kernel_type::Ball({ radius, radius, radius }, false); + + auto dilate = itk::BinaryDilateImageFilter::New(); + progress->RegisterInternalFilter(dilate, 0.1); + dilate->SetInput(this->m_PaddedOutput); + dilate->SetKernel(ball); + dilate->SetForegroundValue(ePixelState::kHardForeground); + dilate->Update(); + return dilate->GetOutput(); } template @@ -132,13 +56,16 @@ FixTopologyCarveOutside::ComputeThinImage using NeighborhoodIteratorType = NeighborhoodIterator>; OutputImagePointer thin_image = this->GetOutput(); - // note output image not padded, all processing is done on m_PaddedOutput, i.e. padded image + + // note: all processing is done on padded_output to avoid + // checking if index is in region + auto padded_output = this->m_PaddedOutput; auto region = thin_image->GetRequestedRegion(); std::vector seeds; for (int direction = 0; direction < ImageDimension; ++direction) { - ImageLinearConstIteratorWithIndex it(m_PaddedOutput, region); + ImageLinearConstIteratorWithIndex it(padded_output, region); it.SetDirection(direction); it.GoToBegin(); for (; !it.IsAtEnd(); it.NextLine()) @@ -177,32 +104,34 @@ FixTopologyCarveOutside::ComputeThinImage // for progress reporting size_t mask_size = 0; - itk::ImageRegionRange image_range(*m_PaddedOutput, region); + itk::ImageRegionRange image_range(*padded_output, region); for (auto && pixel : image_range) { - mask_size += (pixel == kDilated) ? 1 : 0; + mask_size += (pixel == ePixelState::kSoftForeground) ? 1 : 0; } - // use max priority queue: process pixels further away from input mask first + // process pixels further away from input foreground first + // - distance is positive outside + // - use max priority queue using node = std::pair; const auto cmp = [](const node & l, const node & r) { return l.first < r.first; }; std::priority_queue, decltype(cmp)> queue(cmp); for (const auto & idx : seeds) { - m_PaddedOutput->SetPixel(idx, kVisited); - queue.push(std::make_pair(m_DistanceMap->GetPixel(idx), idx)); + padded_output->SetPixel(idx, ePixelState::kQueued); + queue.push(std::make_pair(this->m_DistanceMap->GetPixel(idx), idx)); } auto neighbors = this->GetNeighborOffsets(); const size_t num_neighbors = neighbors.size(); - NeighborhoodIteratorType n_it({ 1, 1, 1 }, m_PaddedOutput, region); + NeighborhoodIteratorType n_it({ 1, 1, 1 }, padded_output, region); - const auto get_mask = [&n_it](const IndexType & idx, const int FG) { + const auto get_mask = [&n_it](const IndexType & idx) { n_it.SetLocation(idx); auto n = n_it.GetNeighborhood(); for (auto & v : n.GetBufferReference()) - v = (v != 0) ? FG : 1 - FG; + v = (v != 0) ? 1 : 0; return n; }; @@ -213,66 +142,47 @@ FixTopologyCarveOutside::ComputeThinImage auto idx = queue.top().second; // node queue.pop(); - if (m_PaddedOutput->GetPixel(idx) != kVisited) + if (padded_output->GetPixel(idx) != ePixelState::kQueued) continue; - auto vals = get_mask(idx, 1); + auto vals = get_mask(idx); // check if point is simple (deletion does not change connectivity in the 3x3x3 neighborhood) if (topology::EulerInvariant(vals, 1) && topology::CCInvariant(vals, 1) && topology::CCInvariant(vals, 0)) { - m_PaddedOutput->SetPixel(idx, kBackground); + padded_output->SetPixel(idx, ePixelState::kBackground); + progress.CompletedPixel(); } // add unvisited neighbors to queue for (size_t k = 0; k < num_neighbors; ++k) { const IndexType n_id = idx + neighbors[k]; - const float n_dist = m_DistanceMap->GetPixel(n_id); + const float n_dist = this->m_DistanceMap->GetPixel(n_id); - if (m_PaddedOutput->GetPixel(n_id) == kDilated) + if (padded_output->GetPixel(n_id) == ePixelState::kSoftForeground) { // mark as visited - m_PaddedOutput->SetPixel(n_id, kVisited); + padded_output->SetPixel(n_id, ePixelState::kQueued); // add to queue queue.push(std::make_pair(n_dist, n_id)); } } - progress.CompletedPixel(); } // copy to output InputImagePointer input_image = dynamic_cast(ProcessObject::GetInput(0)); itk::ImageRegionConstIterator i_it(input_image, region); - itk::ImageRegionConstIterator s_it(m_PaddedOutput, region); + itk::ImageRegionConstIterator s_it(padded_output, region); itk::ImageRegionIterator o_it(thin_image, region); for (i_it.GoToBegin(), s_it.GoToBegin(), o_it.GoToBegin(); !s_it.IsAtEnd(); ++i_it, ++s_it, ++o_it) { - o_it.Set(s_it.Get() != ePixelState::kBackground ? m_InsideValue : i_it.Get()); + o_it.Set(s_it.Get() != ePixelState::kBackground ? this->m_InsideValue : i_it.Get()); } } -template -void -FixTopologyCarveOutside::GenerateData() -{ - ProgressAccumulator::Pointer progress = ProgressAccumulator::New(); - progress->SetMiniPipelineFilter(this); - - this->PrepareData(progress); - - this->ComputeThinImage(progress); -} - -template -void -FixTopologyCarveOutside::PrintSelf(std::ostream & os, Indent indent) const -{ - Superclass::PrintSelf(os, indent); -} - } // end namespace itk #endif // itkFixTopologyCarveOutside_hxx diff --git a/setup.py b/setup.py index a47b9ce..7b33c97 100644 --- a/setup.py +++ b/setup.py @@ -13,7 +13,7 @@ setup( name="itk-topologycontrol", - version="1.0.1", + version="1.1.0", author="Bryn Lloyd", author_email="lloyd@itis.swiss", packages=["itk"], @@ -43,5 +43,5 @@ license="Apache", keywords="ITK InsightToolkit", url=r"https://github.com/dyollb/ITKTopologyControl", - install_requires=[r"itk-filtering~=5.3.0"], + install_requires=[r"itk-filtering==5.3.0"], ) diff --git a/test/Baseline/itkFixTopologyCarveInsideTestOutput.mha.sha512 b/test/Baseline/itkFixTopologyCarveInsideTestOutput.mha.sha512 new file mode 100644 index 0000000..e0317f5 --- /dev/null +++ b/test/Baseline/itkFixTopologyCarveInsideTestOutput.mha.sha512 @@ -0,0 +1 @@ +ded0ee5ea129936d674658f7f8ee3e72b70b328dfc82cc2232d553c90b927c103ea0af0fe9b755e7dcaf5eff983cb871803c93893dd964880bdd91889f82be48 diff --git a/test/Baseline/itkFixTopologyCarveOutsideTestOutput.mha.sha512 b/test/Baseline/itkFixTopologyCarveOutsideTestOutput.mha.sha512 new file mode 100644 index 0000000..b350b04 --- /dev/null +++ b/test/Baseline/itkFixTopologyCarveOutsideTestOutput.mha.sha512 @@ -0,0 +1 @@ +b4f915c08a0453378b8d4bbd96fa2c91978efc38d78b603bc6d2b070bf0c426e7de02ce1eae32a50e919c518b3a5867fa867adfe6857f719193a47906e1904f0 diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 3e655de..1e33906 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -2,17 +2,27 @@ itk_module_test() set(TopologyControlTests itkFixTopologyCarveOutsideTest.cxx + itkFixTopologyCarveInsideTest.cxx ) CreateTestDriver(TopologyControl "${TopologyControl-Test_LIBRARIES}" "${TopologyControlTests}") -set(TEST_DATA_ROOT "${TopologyControl_SOURCE_DIR}/Data") +set(TEST_DATA_ROOT "${CMAKE_CURRENT_SOURCE_DIR}") itk_add_test(NAME itkFixTopologyCarveOutsideTest COMMAND TopologyControlTestDriver --compare - ${TEST_DATA_ROOT}/Baseline/itkFixTopologyCarveOutsideTestOutput.mha + DATA{Baseline/itkFixTopologyCarveOutsideTestOutput.mha} ${ITK_TEST_OUTPUT_DIR}/itkFixTopologyCarveOutsideTestOutput.mha itkFixTopologyCarveOutsideTest ${ITK_TEST_OUTPUT_DIR}/itkFixTopologyCarveOutsideTestOutput.mha ) + +itk_add_test(NAME itkFixTopologyCarveInsideTest + COMMAND TopologyControlTestDriver + --compare + DATA{Baseline/itkFixTopologyCarveInsideTestOutput.mha} + ${ITK_TEST_OUTPUT_DIR}/itkFixTopologyCarveInsideTestOutput.mha + itkFixTopologyCarveInsideTest + ${ITK_TEST_OUTPUT_DIR}/itkFixTopologyCarveInsideTestOutput.mha +) diff --git a/test/itkFixTopologyCarveInsideTest.cxx b/test/itkFixTopologyCarveInsideTest.cxx new file mode 100644 index 0000000..956d85c --- /dev/null +++ b/test/itkFixTopologyCarveInsideTest.cxx @@ -0,0 +1,147 @@ +/*========================================================================= + * + * 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. + * + *=========================================================================*/ + +#include "itkFixTopologyCarveInside.h" + +#include "itkCommand.h" +#include "itkImageFileReader.h" +#include "itkImageFileWriter.h" +#include "itkTestingMacros.h" + +namespace +{ +class ShowProgress : public itk::Command +{ +public: + itkNewMacro(ShowProgress); + + void + Execute(itk::Object * caller, const itk::EventObject & event) override + { + Execute((const itk::Object *)caller, event); + } + + void + Execute(const itk::Object * caller, const itk::EventObject & event) override + { + if (!itk::ProgressEvent().CheckEvent(&event)) + { + return; + } + const auto * processObject = dynamic_cast(caller); + if (!processObject) + { + return; + } + std::cout << " " << processObject->GetProgress(); + } +}; +} // namespace + +int +itkFixTopologyCarveInsideTest(int argc, char * argv[]) +{ + if (argc < 2) + { + std::cerr << "Missing parameters." << std::endl; + std::cerr << "Usage: " << itkNameOfTestExecutableMacro(argv); + std::cerr << " outputImage [inputImage]"; + std::cerr << std::endl; + return EXIT_FAILURE; + } + const char * outputImageFileName = (argc >= 2) ? argv[1] : "itkFixTopologyCarveInsideTestOutput.mha"; + + constexpr unsigned int Dimension = 3; + using PixelType = int; + using ImageType = itk::Image; + using RangeType = itk::ImageRegionRange; + using FilterType = itk::FixTopologyCarveInside; + + auto filter = FilterType::New(); + + ITK_EXERCISE_BASIC_OBJECT_METHODS(filter, FixTopologyCarveInside, FixTopologyBase); + + ImageType::Pointer image; + if (argc <= 2) + { + // Create input image to avoid test dependencies. + image = ImageType::New(); + image->SetRegions({ 128, 128, 128 }); + image->Allocate(); + image->FillBuffer(0); + + ImageType::RegionType region; + region.SetIndex({ 20, 20, 20 }); + region.SetSize({ 20, 20, 20 }); + for (auto & pixel : RangeType(*image, region)) + { + pixel = 1; + } + + region.SetIndex({ 42, 20, 20 }); + region.SetSize({ 20, 20, 20 }); + for (auto & pixel : RangeType(*image, region)) + { + pixel = 1; + } + + region.SetIndex({ 40, 34, 34 }); + region.SetSize({ 2, 5, 5 }); + for (auto & pixel : RangeType(*image, region)) + { + pixel = 1; + } + + image->SetPixel({ 40, 24, 24 }, 1); + image->SetPixel({ 41, 24, 24 }, 1); + } + else + { + const char * inputImageFileName = argv[2]; + auto reader = itk::ImageFileReader::New(); + reader->SetFileName(inputImageFileName); + reader->Update(); + image = reader->GetOutput(); + } + + auto mask = FilterType::MaskImageType::New(); + mask->SetRegions(image->GetBufferedRegion().GetSize()); + mask->Allocate(); + mask->CopyInformation(image); + mask->FillBuffer(0); + mask->SetPixel({ 30, 30, 30 }, 1); + mask->SetPixel({ 51, 30, 30 }, 1); + + auto showProgress = ShowProgress::New(); + filter->AddObserver(itk::ProgressEvent(), showProgress); + filter->SetInput(image); + if (argc <= 2) + filter->SetMaskImage(mask); + else + filter->SetRadius(5); + + auto writer = itk::ImageFileWriter::New(); + writer->SetFileName(outputImageFileName); + writer->SetInput(filter->GetOutput()); + writer->SetUseCompression(true); + + ITK_TRY_EXPECT_NO_EXCEPTION(writer->Update()); + + std::cout << "Test finished." << std::endl; + return EXIT_SUCCESS; +} diff --git a/test/itkFixTopologyCarveOutsideTest.cxx b/test/itkFixTopologyCarveOutsideTest.cxx index 39f4b51..b097bc3 100644 --- a/test/itkFixTopologyCarveOutsideTest.cxx +++ b/test/itkFixTopologyCarveOutsideTest.cxx @@ -73,7 +73,7 @@ itkFixTopologyCarveOutsideTest(int argc, char * argv[]) auto filter = itk::FixTopologyCarveOutside::New(); - ITK_EXERCISE_BASIC_OBJECT_METHODS(filter, FixTopologyCarveOutside, ImageToImageFilter); + ITK_EXERCISE_BASIC_OBJECT_METHODS(filter, FixTopologyCarveOutside, FixTopologyBase); ImageType::Pointer image; if (argc <= 2) @@ -108,13 +108,12 @@ itkFixTopologyCarveOutsideTest(int argc, char * argv[]) image = reader->GetOutput(); } - ShowProgress::Pointer showProgress = ShowProgress::New(); + auto showProgress = ShowProgress::New(); filter->AddObserver(itk::ProgressEvent(), showProgress); filter->SetInput(image); filter->SetRadius(3); - using WriterType = itk::ImageFileWriter; - WriterType::Pointer writer = WriterType::New(); + auto writer = itk::ImageFileWriter::New(); writer->SetFileName(outputImageFileName); writer->SetInput(filter->GetOutput()); writer->SetUseCompression(true); diff --git a/wrapping/itkFixTopologyBase.wrap b/wrapping/itkFixTopologyBase.wrap new file mode 100644 index 0000000..5f617ba --- /dev/null +++ b/wrapping/itkFixTopologyBase.wrap @@ -0,0 +1,6 @@ +itk_wrap_class("itk::FixTopologyBase" POINTER) + foreach(t ${WRAP_ITK_INT}) + itk_wrap_template("${ITKM_I${t}3}" "${ITKT_I${t}3},${ITKT_I${t}3},${ITKT_IUC3}") + endforeach() +itk_end_wrap_class() + diff --git a/wrapping/itkFixTopologyCarveInside.wrap b/wrapping/itkFixTopologyCarveInside.wrap new file mode 100644 index 0000000..57a3505 --- /dev/null +++ b/wrapping/itkFixTopologyCarveInside.wrap @@ -0,0 +1,6 @@ +itk_wrap_class("itk::FixTopologyCarveInside" POINTER) + foreach(t ${WRAP_ITK_INT}) + itk_wrap_template("${ITKM_I${t}3}" "${ITKT_I${t}3},${ITKT_I${t}3},${ITKT_IUC3}") + endforeach() +itk_end_wrap_class() +