diff --git a/Modules/Filtering/Cuberille/CMakeLists.txt b/Modules/Filtering/Cuberille/CMakeLists.txt new file mode 100644 index 000000000000..e17b3840ee69 --- /dev/null +++ b/Modules/Filtering/Cuberille/CMakeLists.txt @@ -0,0 +1,10 @@ +cmake_minimum_required(VERSION 3.16.3) +project(Cuberille) + +if(NOT ITK_SOURCE_DIR) + find_package(ITK 4.10 REQUIRED) + list(APPEND CMAKE_MODULE_PATH ${ITK_CMAKE_DIR}) + include(ITKModuleExternal) +else() + itk_module_impl() +endif() diff --git a/Modules/Filtering/Cuberille/LICENSE b/Modules/Filtering/Cuberille/LICENSE new file mode 100644 index 000000000000..62589edd12a3 --- /dev/null +++ b/Modules/Filtering/Cuberille/LICENSE @@ -0,0 +1,202 @@ + + Apache License + Version 2.0, January 2004 + https://www.apache.org/licenses/ + + TERMS AND CONDITIONS FOR USE, REPRODUCTION, AND DISTRIBUTION + + 1. Definitions. + + "License" shall mean the terms and conditions for use, reproduction, + and distribution as defined by Sections 1 through 9 of this document. + + "Licensor" shall mean the copyright owner or entity authorized by + the copyright owner that is granting the License. + + "Legal Entity" shall mean the union of the acting entity and all + other entities that control, are controlled by, or are under common + control with that entity. For the purposes of this definition, + "control" means (i) the power, direct or indirect, to cause the + direction or management of such entity, whether by contract or + otherwise, or (ii) ownership of fifty percent (50%) or more of the + outstanding shares, or (iii) beneficial ownership of such entity. + + "You" (or "Your") shall mean an individual or Legal Entity + exercising permissions granted by this License. + + "Source" form shall mean the preferred form for making modifications, + including but not limited to software source code, documentation + source, and configuration files. + + "Object" form shall mean any form resulting from mechanical + transformation or translation of a Source form, including but + not limited to compiled object code, generated documentation, + and conversions to other media types. + + "Work" shall mean the work of authorship, whether in Source or + Object form, made available under the License, as indicated by a + copyright notice that is included in or attached to the work + (an example is provided in the Appendix below). + + "Derivative Works" shall mean any work, whether in Source or Object + form, that is based on (or derived from) the Work and for which the + editorial revisions, annotations, elaborations, or other modifications + represent, as a whole, an original work of authorship. For the purposes + of this License, Derivative Works shall not include works that remain + separable from, or merely link (or bind by name) to the interfaces of, + the Work and Derivative Works thereof. + + "Contribution" shall mean any work of authorship, including + the original version of the Work and any modifications or additions + to that Work or Derivative Works thereof, that is intentionally + submitted to Licensor for inclusion in the Work by the copyright owner + or by an individual or Legal Entity authorized to submit on behalf of + the copyright owner. For the purposes of this definition, "submitted" + means any form of electronic, verbal, or written communication sent + to the Licensor or its representatives, including but not limited to + communication on electronic mailing lists, source code control systems, + and issue tracking systems that are managed by, or on behalf of, the + Licensor for the purpose of discussing and improving the Work, but + excluding communication that is conspicuously marked or otherwise + designated in writing by the copyright owner as "Not a Contribution." + + "Contributor" shall mean Licensor and any individual or Legal Entity + on behalf of whom a Contribution has been received by Licensor and + subsequently incorporated within the Work. + + 2. Grant of Copyright License. Subject to the terms and conditions of + this License, each Contributor hereby grants to You a perpetual, + worldwide, non-exclusive, no-charge, royalty-free, irrevocable + copyright license to reproduce, prepare Derivative Works of, + publicly display, publicly perform, sublicense, and distribute the + Work and such Derivative Works in Source or Object form. + + 3. Grant of Patent License. Subject to the terms and conditions of + this License, each Contributor hereby grants to You a perpetual, + worldwide, non-exclusive, no-charge, royalty-free, irrevocable + (except as stated in this section) patent license to make, have made, + use, offer to sell, sell, import, and otherwise transfer the Work, + where such license applies only to those patent claims licensable + by such Contributor that are necessarily infringed by their + Contribution(s) alone or by combination of their Contribution(s) + with the Work to which such Contribution(s) was submitted. If You + institute patent litigation against any entity (including a + cross-claim or counterclaim in a lawsuit) alleging that the Work + or a Contribution incorporated within the Work constitutes direct + or contributory patent infringement, then any patent licenses + granted to You under this License for that Work shall terminate + as of the date such litigation is filed. + + 4. Redistribution. You may reproduce and distribute copies of the + Work or Derivative Works thereof in any medium, with or without + modifications, and in Source or Object form, provided that You + meet the following conditions: + + (a) You must give any other recipients of the Work or + Derivative Works a copy of this License; and + + (b) You must cause any modified files to carry prominent notices + stating that You changed the files; and + + (c) You must retain, in the Source form of any Derivative Works + that You distribute, all copyright, patent, trademark, and + attribution notices from the Source form of the Work, + excluding those notices that do not pertain to any part of + the Derivative Works; and + + (d) If the Work includes a "NOTICE" text file as part of its + distribution, then any Derivative Works that You distribute must + include a readable copy of the attribution notices contained + within such NOTICE file, excluding those notices that do not + pertain to any part of the Derivative Works, in at least one + of the following places: within a NOTICE text file distributed + as part of the Derivative Works; within the Source form or + documentation, if provided along with the Derivative Works; or, + within a display generated by the Derivative Works, if and + wherever such third-party notices normally appear. The contents + of the NOTICE file are for informational purposes only and + do not modify the License. You may add Your own attribution + notices within Derivative Works that You distribute, alongside + or as an addendum to the NOTICE text from the Work, provided + that such additional attribution notices cannot be construed + as modifying the License. + + You may add Your own copyright statement to Your modifications and + may provide additional or different license terms and conditions + for use, reproduction, or distribution of Your modifications, or + for any such Derivative Works as a whole, provided Your use, + reproduction, and distribution of the Work otherwise complies with + the conditions stated in this License. + + 5. Submission of Contributions. Unless You explicitly state otherwise, + any Contribution intentionally submitted for inclusion in the Work + by You to the Licensor shall be under the terms and conditions of + this License, without any additional terms or conditions. + Notwithstanding the above, nothing herein shall supersede or modify + the terms of any separate license agreement you may have executed + with Licensor regarding such Contributions. + + 6. Trademarks. This License does not grant permission to use the trade + names, trademarks, service marks, or product names of the Licensor, + except as required for reasonable and customary use in describing the + origin of the Work and reproducing the content of the NOTICE file. + + 7. Disclaimer of Warranty. Unless required by applicable law or + agreed to in writing, Licensor provides the Work (and each + Contributor provides its Contributions) on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or + implied, including, without limitation, any warranties or conditions + of TITLE, NON-INFRINGEMENT, MERCHANTABILITY, or FITNESS FOR A + PARTICULAR PURPOSE. You are solely responsible for determining the + appropriateness of using or redistributing the Work and assume any + risks associated with Your exercise of permissions under this License. + + 8. Limitation of Liability. In no event and under no legal theory, + whether in tort (including negligence), contract, or otherwise, + unless required by applicable law (such as deliberate and grossly + negligent acts) or agreed to in writing, shall any Contributor be + liable to You for damages, including any direct, indirect, special, + incidental, or consequential damages of any character arising as a + result of this License or out of the use or inability to use the + Work (including but not limited to damages for loss of goodwill, + work stoppage, computer failure or malfunction, or any and all + other commercial damages or losses), even if such Contributor + has been advised of the possibility of such damages. + + 9. Accepting Warranty or Additional Liability. While redistributing + the Work or Derivative Works thereof, You may choose to offer, + and charge a fee for, acceptance of support, warranty, indemnity, + or other liability obligations and/or rights consistent with this + License. However, in accepting such obligations, You may act only + on Your own behalf and on Your sole responsibility, not on behalf + of any other Contributor, and only if You agree to indemnify, + defend, and hold each Contributor harmless for any liability + incurred by, or claims asserted against, such Contributor by reason + of your accepting any such warranty or additional liability. + + END OF TERMS AND CONDITIONS + + APPENDIX: How to apply the Apache License to your work. + + To apply the Apache License to your work, attach the following + boilerplate notice, with the fields enclosed by brackets "[]" + replaced with your own identifying information. (Don't include + the brackets!) The text should be enclosed in the appropriate + comment syntax for the file format. We also recommend that a + file or class name and description of purpose be included on the + same "printed page" as the copyright notice for easier + identification within third-party archives. + + Copyright [yyyy] [name of copyright owner] + + 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 + + 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. diff --git a/Modules/Filtering/Cuberille/README.md b/Modules/Filtering/Cuberille/README.md new file mode 100644 index 000000000000..2a208fcc6324 --- /dev/null +++ b/Modules/Filtering/Cuberille/README.md @@ -0,0 +1,22 @@ +# Cuberille + +In-tree ITK module providing the cuberille implicit-surface +extraction algorithm: given a binary or scalar volume and an +iso-value, produce a quad mesh approximating the iso-surface, with +optional vertex normal interpolation and project-to-iso-surface +refinement. + +The flagship class is `itk::CuberilleImageToMeshFilter`. + +## Origin + +Ingested from the standalone remote module +[**InsightSoftwareConsortium/ITKCuberille**](https://github.com/InsightSoftwareConsortium/ITKCuberille) +on 2026-05-04 via the v4 ingestion pipeline. The upstream repository +will be archived read-only after this PR merges; it remains +reachable at the URL above for historical reference. + +## References + +- Herman G.T., Liu H.K. *Three-dimensional display of human organs from computed tomograms.* Computer Graphics and Image Processing. 1979. +- Mueller D. *Cuberille implicit surface polygonization for ITK.* The Insight Journal. 2010. diff --git a/Modules/Filtering/Cuberille/include/itkCuberilleImageToMeshFilter.h b/Modules/Filtering/Cuberille/include/itkCuberilleImageToMeshFilter.h new file mode 100644 index 000000000000..57a5ed4923ca --- /dev/null +++ b/Modules/Filtering/Cuberille/include/itkCuberilleImageToMeshFilter.h @@ -0,0 +1,436 @@ +/*========================================================================= + * + * 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 itkCuberilleImageToMeshFilter_h +#define itkCuberilleImageToMeshFilter_h + +#include "itkMacro.h" +#include "itkMesh.h" +#include "itkImageToMeshFilter.h" +#include "itkCellInterface.h" +#include "itkTriangleCell.h" +#include "itkQuadrilateralCell.h" +#include "itkDefaultStaticMeshTraits.h" +#include "itkConstShapedNeighborhoodIterator.h" +#include "itkLinearInterpolateImageFunction.h" +#include "itkGradientImageFilter.h" +#include "itkGradientRecursiveGaussianImageFilter.h" +#include "itkVectorLinearInterpolateImageFunction.h" + +#include + +namespace itk +{ + +/** \class CuberilleImageToMeshFilter + * + * This filter uses the 'cuberille' method to convert an implicit surface + * (image) to a mesh. + * + * The 'cuberille' model was proposed over 30 years ago [1,2]. + * A basic summary of the algorithm is as follows: + * step over all pixels, for each pixel determine if it lies on the surface, + * center a cube on the surface pixel, create quadrilateral faces aligned with + * the cube (taking care of neighbouring pixels also on the surface), use a + * gradient descent based method to project cube vertices to iso-surface [3]. + * + * \par Parmeters + * (Required) Input: specifies the input image containing the implicit surface + * for polygonization. Currently, the input must NOT have iso-surface pixels + * on the edge of the image. If this is the case, you MUST use + * itkConstantPadImageFilter to pad the edges by at least 1 pixel. + * + * (Required) IsoSurfaceValue: specifies the value of the iso-surface for which + * to generate the mesh. Pixels equal to or greater than this value are considered + * to lie on the surface or inside the resultant mesh. + * + * (Optional) GenerateTriangleFaces: specifies whether triangle or + * quadrilateral faces should be generated. The default is to generate + * triangle faces. + * + * (Optional) ProjectVerticesToIsoSurface: specifies whether the vertices + * should be projected onto the iso-surface. If projection is disabled, the + * resultant mesh exhibits the traditional blocky features. Projection takes + * roughly half of the algorithm time. The default is to project the vertices. + * + * (Optional) ProjectVertexSurfaceDistanceThreshold: specifies the threshold + * for the 'distance' from iso-surface during vertex projection. Note that the + * distance is actually measured in pixel value units (not space). + * The smaller this value, the closer the vertices will be to the iso-surface. + * Small values result in longer convergence time (i.e. slower). + * Values are clamped to the range [0.0, max pixel value]. + * The default value is 0.5. + * + * (Optional) ProjectVertexStepLength: specifies the threshold for the step + * length during vertex projection. + * The smaller this value, the more likely the vertices will end up closer to + * the surface. Small values cause the projection to take longer to converge. + * Values are clamped to the range [0.0, large]. + * The default value is max spacing * 0.25 (expressed in physical space). + * + * (Optional) ProjectVertexStepLengthRelaxationFactor: specifies the step + * length relaxation factor during vertex projection. The step length is + * multiplied by this factor each iteration to allow convergence. + * Values are clamped to the range [0.0, 1.0].The default value is 0.95. + * + * (Optional) ProjectVertexMaximumNumberOfSteps: specifies the maximum number + * of steps used during vertex projection. The default value is 50. + * + * \par References + * [1] G. Herman and H. Liu, "Three-dimensional Display of Human organs + * from Computed Tomograms", Computer Graphics and Images Processing, + * Volume 9, Issue 1, Pages 1-21, 1979. + * [2] D. Gordon and J.K. Udupa, "Fast surface tracking in Three-dimensional + Binary Images", Computer Vision, Graphics and Image Processing, + Volume 45, Pages 196-214, 1989. + * [3] https://www2.imm.dtu.dk/~jab/gallery/polygonization.html + * + * This implementation was taken from the Insight Journal: + * https://hdl.handle.net/10380/3186 + * + * \author Dan Mueller, Philips Healthcare, dan dot muel at gmail dot com + * + * \ingroup Cuberille + * + */ +template > +class CuberilleImageToMeshFilter : public ImageToMeshFilter +{ +public: + ITK_DISALLOW_COPY_AND_MOVE(CuberilleImageToMeshFilter); + + /** Standard "Self" type alias. */ + using Self = CuberilleImageToMeshFilter; + using Superclass = ImageToMeshFilter; + using Pointer = SmartPointer; + using ConstPointer = SmartPointer; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Run-time type information (and related methods). */ + itkTypeMacro(CuberilleImageToMeshFilter, ImageToMeshFilter); + + /** Compile-time configuration flags (formerly preprocessor macros). + * Flip a flag to true and rebuild to enable the corresponding alternative + * algorithm path; the dead branch is dropped by the C++ compiler via + * `if constexpr`. */ + /** @ITKStartGrouping */ + static constexpr bool UseGradientRecursiveGaussian = false; + static constexpr bool UseAdvancedProjection = false; + static constexpr bool UseLineSearchProjection = false; + /** @ITKEndGrouping */ + + /** Some convenient type alias. */ + using OutputMeshType = TOutputMesh; + using OutputMeshPointer = typename OutputMeshType::Pointer; + using OutputMeshTraits = typename OutputMeshType::MeshTraits; + using OutputPointType = typename OutputMeshType::PointType; + using OutputPixelType = typename OutputMeshTraits::PixelType; + using CellTraits = typename OutputMeshType::CellTraits; + using PointsContainerPointer = typename OutputMeshType::PointsContainerPointer; + using PointsContainer = typename OutputMeshType::PointsContainer; + using CellsContainerPointer = typename OutputMeshType::CellsContainerPointer; + using CellsContainer = typename OutputMeshType::CellsContainer; + using PointIdentifier = typename OutputMeshType::PointIdentifier; + using PointVectorType = typename std::vector; + using CellIdentifier = typename OutputMeshType::CellIdentifier; + using CellInterfaceType = CellInterface; + using TriangleCellType = TriangleCell; + using TriangleAutoPointer = typename TriangleCellType::SelfAutoPointer; + using TriangleCellAutoPointer = typename TriangleCellType::CellAutoPointer; + using QuadrilateralCellType = QuadrilateralCell; + using QuadrilateralAutoPointer = typename QuadrilateralCellType::SelfAutoPointer; + using QuadrilateralCellAutoPointer = typename QuadrilateralCellType::CellAutoPointer; + + using InputImageType = TInputImage; + using InputImagePointer = typename InputImageType::Pointer; + using InputImageConstPointer = typename InputImageType::ConstPointer; + using InputPixelType = typename InputImageType::PixelType; + using SizeType = typename InputImageType::SizeType; + using SpacingType = typename InputImageType::SpacingType; + using SpacingValueType = typename InputImageType::SpacingValueType; + using IndexType = typename InputImageType::IndexType; + using PointType = typename OutputMeshType::PointType; + + using InterpolatorType = TInterpolator; + using InterpolatorPointer = typename InterpolatorType::Pointer; + using InterpolatorOutputType = typename InterpolatorType::OutputType; + + /** Other convenient type alias. */ + using InputImageIteratorType = ConstShapedNeighborhoodIterator; + using GradientFilterType = std::conditional_t, + GradientImageFilter>; + using GradientFilterPointer = typename GradientFilterType::Pointer; + using GradientImageType = typename GradientFilterType::OutputImageType; + using GradientImagePointer = typename GradientImageType::Pointer; + using GradientPixelType = typename GradientFilterType::OutputPixelType; + using GradientInterpolatorType = itk::VectorLinearInterpolateImageFunction; + using GradientInterpolatorPointer = typename GradientInterpolatorType::Pointer; + + /** Get/Set the iso-surface value. + * This parameter specifies the value of the iso-surface for which to + * generate the mesh. Pixels equal to or less than this value are + * considered on the surface or inside the resultant mesh. + */ + itkGetMacro(IsoSurfaceValue, InputPixelType); + itkSetMacro(IsoSurfaceValue, InputPixelType); + + /** Accept the input image. */ + using Superclass::SetInput; + virtual void + SetInput(const InputImageType * inputImage); + + /** Get/Set interpolate function. */ + itkGetConstObjectMacro(Interpolator, InterpolatorType); + itkSetObjectMacro(Interpolator, InterpolatorType); + + /** Get/Set whether triangle or quadrilateral faces should be generated. + * True = triangle faces, False = quadrilateral faces. + Default = true (triangle faces). */ + itkGetMacro(GenerateTriangleFaces, bool); + itkSetMacro(GenerateTriangleFaces, bool); + itkBooleanMacro(GenerateTriangleFaces); + + /** Get/Set whether the vertices should be project to the iso-surface. + Default = true. */ + itkGetMacro(ProjectVerticesToIsoSurface, bool); + itkSetMacro(ProjectVerticesToIsoSurface, bool); + itkBooleanMacro(ProjectVerticesToIsoSurface); + + /** Get/Set whether the adjacent input pixel value should be saved as cell data in the output mesh. + Default = false. */ + itkGetMacro(SavePixelAsCellData, bool); + itkSetMacro(SavePixelAsCellData, bool); + itkBooleanMacro(SavePixelAsCellData); + + /** Get/Set the threshold for the "distance" from iso-surface during vertex projection. + Note that the distance is actually measured in pixel value units (not space). + The smaller this value, the closer the vertices will be to the iso-surface. + Small values result in longer convergence time (i.e. slower). + Values are clamped to the range [0.0, max pixel value]. + Default = 0.5. */ + itkGetMacro(ProjectVertexSurfaceDistanceThreshold, double); + itkSetClampMacro(ProjectVertexSurfaceDistanceThreshold, double, 0.0, NumericTraits::max()); + + /** Get/Set the the initial step length for vertex projection. + Values are clamped to the range [0.0, large]. + Default = max spacing * 0.25 (expressed in physical space). */ + itkGetMacro(ProjectVertexStepLength, double); + itkSetClampMacro(ProjectVertexStepLength, double, 0.0, 100000.0); + + /** Get/Set the step length relaxation factor during vertex projection. + The step length is multiplied by this factor each iteration to allow convergence. + Values are clamped to the range [0.0, 1.0]. + Default = 0.95. */ + itkGetMacro(ProjectVertexStepLengthRelaxationFactor, double); + itkSetClampMacro(ProjectVertexStepLengthRelaxationFactor, double, 0.0, 1.0); + + /** Get/Set the maximum number of steps used during vertex projection. + Default = 50. */ + itkGetMacro(ProjectVertexMaximumNumberOfSteps, unsigned int); + itkSetMacro(ProjectVertexMaximumNumberOfSteps, unsigned int); + + /** Calculate connected components labels for all possible 2x2x2 binary images. */ + void + CalculateLabelsArray(); + +protected: + CuberilleImageToMeshFilter(); + ~CuberilleImageToMeshFilter() override; + void + PrintSelf(std::ostream & os, Indent indent) const override; + + void + GenerateData() override; + void + GenerateOutputInformation() override {}; // do nothing + +private: + /** \class VertexLookupNode A private class containing lookup details for vertices. + * \ingroup Cuberille */ + class VertexLookupNode + { + public: + /** Convenient type alias */ + using Self = VertexLookupNode; + + /** Constructors */ + VertexLookupNode() = default; + VertexLookupNode(unsigned long x, unsigned long y) + : m_X(x) + , m_Y(y) + {} + + /** Parameters */ + unsigned long + GetX() + { + return m_X; + } + unsigned long + GetY() + { + return m_Y; + } + + /** Comparison operators for sorting */ + bool + operator>(const Self & node) const + { + return (m_Y > node.m_Y) || ((m_Y == node.m_Y) && (m_X > node.m_X)); + } + bool + operator>=(const Self & node) const + { + return (m_Y > node.m_Y) || ((m_Y == node.m_Y) && (m_X >= node.m_X)); + } + bool + operator<(const Self & node) const + { + return (m_Y < node.m_Y) || ((m_Y == node.m_Y) && (m_X < node.m_X)); + } + bool + operator<=(const Self & node) const + { + return (m_Y < node.m_Y) || ((m_Y == node.m_Y) && (m_X <= node.m_X)); + } + + private: + unsigned long m_X{ 0 }; + unsigned long m_Y{ 0 }; + }; + + /** \class VertexLookupMap A private class providing vertex lookup functionality. + * \ingroup Cuberille */ + template + class VertexLookupMap + { + public: + /** Convenient type alias */ + using Self = VertexLookupMap; + using MapType = std::map; + + /** Constructors */ + VertexLookupMap() = default; + + /** Clear the lookup map. */ + void + Clear() + { + m_Map.clear(); + } + + /** Add the given vertex identifer to the given [x,y] position. */ + void + AddVertex(unsigned int x, unsigned int y, PointVectorType ids) + { + VertexLookupNode node(x, y); + m_Map.insert(typename MapType::value_type(node, ids)); + } + + /** Get the vertex identifer for the given [x,y] position. + * Returns true if the vertex exists and id contains the identifer. + * Returns false if the vertex does not exist and id is undefined. */ + bool + GetVertex(unsigned int x, unsigned int y, const size_t component, PointIdentifier & id) + { + bool result = false; + VertexLookupNode node(x, y); + auto it = m_Map.find(node); + if (it != m_Map.end()) + { + result = true; + id = it->second.at(component); + } + return result; + } + + private: + MapType m_Map; + }; + + /** Some convenient type alias. */ + using VertexLookupMapType = VertexLookupMap; + + /** Private functions to implement the algorithm. */ + + /** Compute gradient image. */ + inline void + ComputeGradientImage(); + + /** Set a flag activating each vertex for the given face. */ + inline void + SetVerticesFromFace(unsigned int face, std::array & vertexHasQuad); + + /** Get the vertex lookup index from the given index and vertex number. */ + inline IndexType + GetVertexLookupIndex(unsigned int vertex, IndexType index); + + /** Project vertex to the iso-surface by stepping along normal. */ + inline void + ProjectVertexToIsoSurface(PointType & vertex); + + /** Add a vertex to the given mesh. Increments point identifier. */ + inline PointVectorType + AddVertex(PointIdentifier & id, + IndexType index, + const InputImageType * image, + OutputMeshType * mesh, + const size_t numComponents); + + /** Add quadrilateral face to the given mesh. Increments cell identifier. */ + inline void + AddQuadFace(CellIdentifier & id, + std::array f, + OutputMeshType * mesh, + const InputPixelType & pixel); + + /** Calculate the local 2x2x2 bitmask for a given vertex index. */ + size_t + CalculateBitmaskIDForVertexIndex(const IndexType & vindex); + + using TLabel = signed char; + using TLabels = std::array; + using TLabelsArray = std::array; + + TLabelsArray m_LabelsArray; + + InputPixelType m_IsoSurfaceValue; + InterpolatorPointer m_Interpolator; + GradientInterpolatorPointer m_GradientInterpolator; + SpacingValueType m_MaxSpacing; + bool m_GenerateTriangleFaces{ true }; + bool m_ProjectVerticesToIsoSurface{ true }; + bool m_SavePixelAsCellData{ false }; + double m_ProjectVertexSurfaceDistanceThreshold{ 0.5 }; + double m_ProjectVertexStepLength{ -1.0 }; + double m_ProjectVertexStepLengthRelaxationFactor{ 0.95 }; + unsigned int m_ProjectVertexMaximumNumberOfSteps{ 50 }; +}; + +} // end namespace itk + +#ifndef ITK_MANUAL_INSTANTIATION +# include "itkCuberilleImageToMeshFilter.hxx" +#endif + +#endif diff --git a/Modules/Filtering/Cuberille/include/itkCuberilleImageToMeshFilter.hxx b/Modules/Filtering/Cuberille/include/itkCuberilleImageToMeshFilter.hxx new file mode 100644 index 000000000000..e95cbde8bad2 --- /dev/null +++ b/Modules/Filtering/Cuberille/include/itkCuberilleImageToMeshFilter.hxx @@ -0,0 +1,1008 @@ +/*========================================================================= + * + * 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 itkCuberilleImageToMeshFilter_hxx +#define itkCuberilleImageToMeshFilter_hxx + +#include "itkMath.h" +#include "itkNumericTraits.h" +#include +#include +#include +#include + +namespace itk +{ + +template +CuberilleImageToMeshFilter::CuberilleImageToMeshFilter() + : m_IsoSurfaceValue(NumericTraits::One) + , m_MaxSpacing(NumericTraits::One) + +{ + this->SetNumberOfRequiredInputs(1); + this->CalculateLabelsArray(); +} + +template +CuberilleImageToMeshFilter::~CuberilleImageToMeshFilter() +{ + m_GradientInterpolator = nullptr; +} + +template +void +CuberilleImageToMeshFilter::SetInput(const InputImageType * image) +{ + this->ProcessObject::SetNthInput(0, const_cast(image)); +} + +template +void +CuberilleImageToMeshFilter::GenerateData() +{ + + // Get input/output + InputImageConstPointer image = Superclass::GetInput(0); + typename OutputMeshType::Pointer mesh = Superclass::GetOutput(); + + // Compute maximum spacing + m_MaxSpacing = image->GetSpacing().GetVnlVector().max_value(); + + // Set default step length + if (m_ProjectVertexStepLength < 0.0) + { + m_ProjectVertexStepLength = m_MaxSpacing * 0.25; + } + + // Create interpolator for pixel value + if (m_Interpolator.IsNull()) + { + m_Interpolator = LinearInterpolateImageFunction::New(); + } + m_Interpolator->SetInputImage(image); + + // Create interpolator for gradient image + ComputeGradientImage(); + + // Set up iterator + typename InputImageIteratorType::SizeType radius; + radius.Fill(1); + InputImageIteratorType it(radius, image, image->GetBufferedRegion()); + setConnectivity(&it, false); // Set face connectivity + + // TODO: Estimate number of vertices/faces for which to reserve + // TODO: There is an issue with quad edge mesh related to reserving, see + // https://www.itk.org/mailman/private/insight-developers/2010-June/014653.html + // SizeType size = image->GetLargestPossibleRegion().GetSize(); + // mesh->GetPoints()->Reserve( (size[0]*size[1]*size[2]) / 200 ); + // mesh->GetCells()->Reserve( (size[0]*size[1]*size[2]) / 100 ); + + // Set up helper structures + unsigned int look0 = 1; + unsigned int look1 = 0; + unsigned char numFaces = 0; + std::array faceHasQuad; + std::array vertexHasQuad; + std::array v; + std::array f; + PointIdentifier nextVertexId = 0; + CellIdentifier nextCellId = 0; + IndexType index; + typename IndexType::IndexValueType lastZ = -1; + InputPixelType center = 0; + std::array offset; + offset[0][0] = -1; + offset[0][1] = +0; + offset[0][2] = +0; // -X + offset[1][0] = +0; + offset[1][1] = -1; + offset[1][2] = +0; // -Y + offset[2][0] = +1; + offset[2][1] = +0; + offset[2][2] = +0; // +X + offset[3][0] = +0; + offset[3][1] = +1; + offset[3][2] = +0; // +Y + offset[4][0] = +0; + offset[4][1] = +0; + offset[4][2] = -1; // -Z + offset[5][0] = +0; + offset[5][1] = +0; + offset[5][2] = +1; // +Z + VertexLookupMapType lookup[2]; + lookup[look0].Clear(); + lookup[look1].Clear(); + + // TODO: Handle voxels on the edge of the image + + // Iterate input image + for (it.GoToBegin(); !it.IsAtEnd(); ++it) + { + // Determine if current pixel is suitable + center = it.GetCenterPixel(); + if (center < m_IsoSurfaceValue) + { + continue; + } + + // Re-initialize for new pixel + numFaces = 0; + faceHasQuad.fill(false); + vertexHasQuad.fill(false); + + // Re-initialize for new z plane + index = it.GetIndex(); + if (index[2] != lastZ) + { + std::swap(look0, look1); + lookup[look1].Clear(); + lastZ = index[2]; + } + + // Compute which faces (if any) have quads + for (unsigned int i = 0; i < 6; ++i) + { + // NOTE: suitability check above means center <= m_IsoSurfaceValue + faceHasQuad[i] = (it.GetPixel(offset[i]) < m_IsoSurfaceValue); + if (faceHasQuad[i]) + { + numFaces++; + SetVerticesFromFace(i, vertexHasQuad); + } + } + + // Process each face + if (numFaces > 0) + { + // Create vertices + for (unsigned int i = 0; i < 8; ++i) + { + if (!vertexHasQuad[i]) + { + continue; + } + + // Use the vertex lookup to get the vertex for the correct slice + IndexType vindex = GetVertexLookupIndex(i, index); + const auto bitmaskID = this->CalculateBitmaskIDForVertexIndex(vindex); + const auto diff = vindex - index; + size_t offsetID = 7 - (diff[0] + diff[1] * 2 + diff[2] * 4); + const auto components = this->m_LabelsArray.at(bitmaskID); + const auto component = components.at(offsetID); + unsigned int look = (i < 4) ? look0 : look1; // First four are first slice + if (!lookup[look].GetVertex(vindex[0], vindex[1], component, v[i])) + { + // Vertex was not in lookup, create and add to lookup + v[i] = nextVertexId; + const auto numComponents = (*std::max_element(components.begin(), components.end())) + 1; + const auto pv = AddVertex(nextVertexId, vindex, image, mesh, numComponents); + lookup[look].AddVertex(vindex[0], vindex[1], pv); + } + + } // end foreach vertex + + // Create faces + if (faceHasQuad[0]) + { + f[0] = v[0]; + f[1] = v[4]; + f[2] = v[7]; + f[3] = v[3]; + AddQuadFace(nextCellId, f, mesh, center); + } + if (faceHasQuad[1]) + { + f[0] = v[0]; + f[1] = v[1]; + f[2] = v[5]; + f[3] = v[4]; + AddQuadFace(nextCellId, f, mesh, center); + } + if (faceHasQuad[2]) + { + f[0] = v[1]; + f[1] = v[2]; + f[2] = v[6]; + f[3] = v[5]; + AddQuadFace(nextCellId, f, mesh, center); + } + if (faceHasQuad[3]) + { + f[0] = v[2]; + f[1] = v[3]; + f[2] = v[7]; + f[3] = v[6]; + AddQuadFace(nextCellId, f, mesh, center); + } + if (faceHasQuad[4]) + { + f[0] = v[0]; + f[1] = v[3]; + f[2] = v[2]; + f[3] = v[1]; + AddQuadFace(nextCellId, f, mesh, center); + } + if (faceHasQuad[5]) + { + f[0] = v[4]; + f[1] = v[5]; + f[2] = v[6]; + f[3] = v[7]; + AddQuadFace(nextCellId, f, mesh, center); + } + + } // end if num faces > 0 + } +} + +template +void +CuberilleImageToMeshFilter::SetVerticesFromFace(unsigned int face, + std::array & v) +{ + switch (face) + { + case 0: + v[0] = true; + v[4] = true; + v[7] = true; + v[3] = true; + break; + case 1: + v[0] = true; + v[1] = true; + v[5] = true; + v[4] = true; + break; + case 2: + v[1] = true; + v[2] = true; + v[6] = true; + v[5] = true; + break; + case 3: + v[2] = true; + v[3] = true; + v[7] = true; + v[6] = true; + break; + case 4: + v[0] = true; + v[3] = true; + v[2] = true; + v[1] = true; + break; + case 5: + v[4] = true; + v[5] = true; + v[6] = true; + v[7] = true; + break; + } +} + +template +typename TInputImage::IndexType +CuberilleImageToMeshFilter::GetVertexLookupIndex( + unsigned int vertex, + typename TInputImage::IndexType index) +{ + IndexType result(index); + switch (vertex) + { + case 0: + break; + case 1: + result[0] += 1; + break; + case 2: + result[0] += 1; + result[1] += 1; + break; + case 3: + result[1] += 1; + break; + case 4: + result[2] += 1; + break; + case 5: + result[0] += 1; + result[2] += 1; + break; + case 6: + result[0] += 1; + result[1] += 1; + result[2] += 1; + break; + case 7: + result[1] += 1; + result[2] += 1; + break; + } + return result; +} + +template +typename CuberilleImageToMeshFilter::PointVectorType +CuberilleImageToMeshFilter::AddVertex( + typename TOutputMesh::PointIdentifier & id, + typename TInputImage::IndexType index, + const TInputImage * image, + TOutputMesh * mesh, + const size_t numComponents) +{ + PointType vertex; + image->TransformIndexToPhysicalPoint(index, vertex); + const auto spacing = image->GetSpacing(); + const auto direction = image->GetDirection(); + const auto offset = direction * spacing * 0.5; + vertex -= offset; + if (m_ProjectVerticesToIsoSurface) + { + ProjectVertexToIsoSurface(vertex); + } + PointVectorType pointIDVector; + for (size_t i = 0; i < numComponents; ++i) + { + pointIDVector.push_back(id); + mesh->GetPoints()->InsertElement(id++, vertex); + } + return pointIDVector; +} + +template +void +CuberilleImageToMeshFilter::AddQuadFace( + typename TOutputMesh::CellIdentifier & id, + std::array face, + TOutputMesh * mesh, + const InputPixelType & pixel) +{ + if (m_GenerateTriangleFaces) + { + // Get vertices + PointType v[4]; + for (unsigned int i = 0; i < 4; ++i) + { + v[i] = mesh->GetPoints()->GetElement(face[i]); + } + + // Split the quad along the longest edge to avoid skinny triangles + PointIdentifier face1[3]; + PointIdentifier face2[3]; + if (v[0].SquaredEuclideanDistanceTo(v[2]) >= v[1].SquaredEuclideanDistanceTo(v[3])) + { + face1[0] = face[0]; + face1[1] = face[1]; + face1[2] = face[3]; + face2[0] = face[1]; + face2[1] = face[2]; + face2[2] = face[3]; + } + else + { + face1[0] = face[0]; + face1[1] = face[1]; + face1[2] = face[2]; + face2[0] = face[0]; + face2[1] = face[2]; + face2[2] = face[3]; + } + + // Add triangle 1 cell + TriangleCellAutoPointer tri1; + tri1.TakeOwnership(new TriangleCellType); + tri1->SetPointIds(face1); + mesh->SetCell(id++, tri1); + if (this->m_SavePixelAsCellData) + { + mesh->SetCellData((id - 1), pixel); + } + + // Add triangle 2 cell + TriangleCellAutoPointer tri2; + tri2.TakeOwnership(new TriangleCellType); + tri2->SetPointIds(face2); + mesh->SetCell(id++, tri2); + if (this->m_SavePixelAsCellData) + { + mesh->SetCellData((id - 1), pixel); + } + } + else + { + // Add quateraleral cell + QuadrilateralCellAutoPointer quad1; + quad1.TakeOwnership(new QuadrilateralCellType); + quad1->SetPointIds(face.data()); + mesh->SetCell(id++, quad1); + if (this->m_SavePixelAsCellData) + { + mesh->SetCellData((id - 1), pixel); + } + } +} + +template +void +CuberilleImageToMeshFilter::ProjectVertexToIsoSurface(PointType & vertex) +{ + if constexpr (UseAdvancedProjection) + { + // Set up + bool done = false; + double step = m_ProjectVertexStepLength; + unsigned int numberOfSteps = 0; + PointType temp[2]; + InterpolatorOutputType value[2]; + double diff[2]; + unsigned int i, swaps = 0; + int previousi = -1; + + while (!done) + { + // Compute normal vector + GradientPixelType normal = m_GradientInterpolator->Evaluate(vertex); + if (normal.Normalize() == 0.0) // old norm was zero + { + break; + } + + // Step along both directions of normal + for (i = 0; i < InputImageType::ImageDimension; ++i) + { + temp[0][i] = vertex[i] + (normal[i] * +1.0 * step); + temp[1][i] = vertex[i] + (normal[i] * -1.0 * step); + } + step *= m_ProjectVertexStepLengthRelaxationFactor; + + // Compute which direction moves vertex closer to iso-surface value + value[0] = m_Interpolator->Evaluate(temp[0]); + value[1] = m_Interpolator->Evaluate(temp[1]); + diff[0] = itk::Math::abs(value[0] - m_IsoSurfaceValue); + diff[1] = itk::Math::abs(value[1] - m_IsoSurfaceValue); + i = (diff[0] <= diff[1]) ? 0 : 1; + if (previousi < 0) + { + previousi = i; + } + swaps += (int)(previousi != i); + vertex = temp[i]; + + // Determine whether vertex is close enough to iso-surface value + done |= diff[i] < m_ProjectVertexSurfaceDistanceThreshold; + if (done) + { + break; + } + + // Determine whether we have done enough steps + done |= numberOfSteps++ > m_ProjectVertexMaximumNumberOfSteps; + if (done) + { + break; + } + + // Determine whether there has been too many sign swaps (oscillating) + done |= (swaps >= 5); + if (done) + { + break; + } + } + } + else if constexpr (UseLineSearchProjection) + { + + // Set up + GradientPixelType normal; + InterpolatorOutputType value; + PointType temp, bestVertex; + double sign, d, metric, bestMetric = 10000; + + // Compute normal vector + normal = m_GradientInterpolator->Evaluate(vertex); + if (normal.Normalize() == 0.0) // old norm was zero + { + return; + } + + // Search on both sides of the line + for (sign = -1.0; sign <= 1.0; sign += 2.0) + { + for (unsigned int j = 1; j < m_ProjectVertexMaximumNumberOfSteps / 2; ++j) + { + // Compute current location along line + d = (double)j / ((double)m_ProjectVertexMaximumNumberOfSteps / 2.0); + for (unsigned int i = 0; i < InputImageType::ImageDimension; ++i) + { + temp[i] = vertex[i] + (normal[i] * sign * m_ProjectVertexStepLength * d); + } + // Compute metric (combination of difference and distance) + value = m_Interpolator->Evaluate(temp); + metric = itk::Math::abs(value - m_IsoSurfaceValue); // Difference + + // Determine if current position is the "best" + if (metric < bestMetric) + { + bestMetric = metric; + bestVertex = temp; + } + } + } + vertex = bestVertex; + } + else + { + // Set up + bool done = false; + double sign = 1.0; + double step = m_ProjectVertexStepLength; + unsigned int numberOfSteps = 0; + GradientPixelType normal; + InterpolatorOutputType value; + + while (!done) + { + // Compute normal vector + normal = m_GradientInterpolator->Evaluate(vertex); + if (normal.Normalize() == 0.0) // old norm was zero + { + break; + } + + // Compute whether vertex is close enough to iso-surface value + value = m_Interpolator->Evaluate(vertex); + done |= itk::Math::abs(value - m_IsoSurfaceValue) < m_ProjectVertexSurfaceDistanceThreshold; + if (done) + { + break; + } + + // Step along the normal towards the iso-surface value + sign = (value < m_IsoSurfaceValue) ? +1.0 : -1.0; + for (unsigned int i = 0; i < InputImageType::ImageDimension; ++i) + { + vertex[i] += (normal[i] * sign * step); + } + step *= m_ProjectVertexStepLengthRelaxationFactor; + done |= numberOfSteps++ > m_ProjectVertexMaximumNumberOfSteps; + } + } +} + +template +void +CuberilleImageToMeshFilter::ComputeGradientImage() +{ + if (m_ProjectVerticesToIsoSurface && m_GradientInterpolator.IsNull()) + { + typename GradientFilterType::Pointer gradientFilter = GradientFilterType::New(); + gradientFilter->SetInput(Superclass::GetInput(0)); + if constexpr (UseGradientRecursiveGaussian) + { + gradientFilter->SetSigma(m_MaxSpacing * 1.0); + gradientFilter->SetNormalizeAcrossScale(true); + } + gradientFilter->Update(); + m_GradientInterpolator = GradientInterpolatorType::New(); + m_GradientInterpolator->SetInputImage(gradientFilter->GetOutput()); + gradientFilter->GetOutput()->DisconnectPipeline(); + gradientFilter = nullptr; + } +} + +template +size_t +CuberilleImageToMeshFilter::CalculateBitmaskIDForVertexIndex( + const IndexType & vindex) +{ + typename IndexType::OffsetType ones; + ones.Fill(1); + const IndexType localorigin = vindex - ones; + std::bitset<8> bitmask; + for (size_t i = 0; i < 8; ++i) + { + std::bitset<3> bits(i); + typename IndexType::OffsetType offset; + offset[0] = bits[0]; + offset[1] = bits[1]; + offset[2] = bits[2]; + const auto index = localorigin + offset; + const auto pixel = this->GetInput()->GetPixel(index); + bitmask[i] = (pixel >= this->m_IsoSurfaceValue); + } + return static_cast(bitmask.to_ulong()); +} + +template +void +CuberilleImageToMeshFilter::CalculateLabelsArray() +{ + + // The commented code below iterates through all possible binary 2x2x2 regions, + // and runs connected components on the result. The total number of foreground + // connected components is equal to the number of vertices which must be + // replicated--these vertices are stored in the VertexMap. One less than the + // number of the connected component is equal to the index of the corresponding + // vertex in the VertexMap. + // + // This method works for all but four of the 256 possible cases. The failing + // cases occur when there are two background pixels at opposite corners of the + // 2x2x2 region, and all remaining pixels are foreground. In these cases, + // there is only one connected component, but two vertices must be added. These + // cases are handled separately. + // + // While it would be possible to generate this at runtime, it could theoretically + // cause a significant performance hit in a case where many instances of this class + // must be instantiated. Therefore, the values are hardcoded below. The commented + // code was used to generate the hardcoded values. + // + // NOTE: There was previously a bug in itk::ConnectedComponentImageFilter which + // caused the hardcoded values below to be incorrect, though this has since been + // fixed. Therefore, the commented code below will fail if built prior to the + // following patch: + // + // Git Hash: c32a7846cb6502eaa79780de3ff7b0ce81597118 + // + // using TImage = itk::Image; + // using TConnected = itk::ConnectedComponentImageFilter< TImage, TImage >; + // + // for (size_t mask = 0; mask < std::pow(2, 8); ++mask) { + // + // std::bitset<8> bitmask(mask); + // + // const auto image = TImage::New(); + // TImage::SizeType size; + // size.Fill( 2 ); + // + // TImage::IndexType origin; + // origin.Fill( 0 ); + // + // TImage::RegionType region(origin, size); + // + // image->SetRegions( region ); + // image->Allocate(); + // image->FillBuffer( 0 ); + // + // for (size_t index = 0; index < std::pow(2, 3); ++index) { + // std::bitset<3> bitindex(index); + // image->SetPixel( {{bitindex[0], bitindex[1], bitindex[2]}}, bitmask[index] ); + // } + // + // const auto connected = TConnected::New(); + // connected->SetInput( image ); + // connected->FullyConnectedOff(); + // connected->Update(); + // + // TLabels labels; + // + // for (size_t index = 0; index < std::pow(2, 3); ++index) { + // std::bitset<3> bitindex(index); + // const auto component = connected->GetOutput()->GetPixel({{ + // bitindex[0], + // bitindex[1], + // bitindex[2] + // }}); + // labels[index] = component - 1; + // } + // + // this->m_LabelsArray[mask] = labels; + // + // } + // + // // Manually handle the corner cases, discussed above. + // + // this->m_LabelsArray[126] = { { -1, 0, 0, 1, 0, 1, 1, -1 } }; + // this->m_LabelsArray[189] = { { 0, -1, 1, 0, 1, 0, -1, 1 } }; + // this->m_LabelsArray[219] = { { 0, 1, -1, 0, 1, -1, 0, 1 } }; + // this->m_LabelsArray[231] = { { 0, 1, 1, -1, -1, 0, 0, 1 } }; + // + // for (size_t i = 0; i < std::pow(2,8); ++i) { + // std::cout << "this->m_LabelsArray[" << i << "] "; + // std::cout << "= { {"; + // for (size_t j = 0; j < std::pow(2,3); ++j) { + // const auto k = static_cast(this->m_LabelsArray[i][j]); + // std::cout << ' ' << k; + // if (j != 7) std::cout << ','; + // } + // std::cout << " } };\n"; + // } + + this->m_LabelsArray[0] = { { -1, -1, -1, -1, -1, -1, -1, -1 } }; + this->m_LabelsArray[1] = { { 0, -1, -1, -1, -1, -1, -1, -1 } }; + this->m_LabelsArray[2] = { { -1, 0, -1, -1, -1, -1, -1, -1 } }; + this->m_LabelsArray[3] = { { 0, 0, -1, -1, -1, -1, -1, -1 } }; + this->m_LabelsArray[4] = { { -1, -1, 0, -1, -1, -1, -1, -1 } }; + this->m_LabelsArray[5] = { { 0, -1, 0, -1, -1, -1, -1, -1 } }; + this->m_LabelsArray[6] = { { -1, 0, 1, -1, -1, -1, -1, -1 } }; + this->m_LabelsArray[7] = { { 0, 0, 0, -1, -1, -1, -1, -1 } }; + this->m_LabelsArray[8] = { { -1, -1, -1, 0, -1, -1, -1, -1 } }; + this->m_LabelsArray[9] = { { 0, -1, -1, 1, -1, -1, -1, -1 } }; + this->m_LabelsArray[10] = { { -1, 0, -1, 0, -1, -1, -1, -1 } }; + this->m_LabelsArray[11] = { { 0, 0, -1, 0, -1, -1, -1, -1 } }; + this->m_LabelsArray[12] = { { -1, -1, 0, 0, -1, -1, -1, -1 } }; + this->m_LabelsArray[13] = { { 0, -1, 0, 0, -1, -1, -1, -1 } }; + this->m_LabelsArray[14] = { { -1, 0, 0, 0, -1, -1, -1, -1 } }; + this->m_LabelsArray[15] = { { 0, 0, 0, 0, -1, -1, -1, -1 } }; + this->m_LabelsArray[16] = { { -1, -1, -1, -1, 0, -1, -1, -1 } }; + this->m_LabelsArray[17] = { { 0, -1, -1, -1, 0, -1, -1, -1 } }; + this->m_LabelsArray[18] = { { -1, 0, -1, -1, 1, -1, -1, -1 } }; + this->m_LabelsArray[19] = { { 0, 0, -1, -1, 0, -1, -1, -1 } }; + this->m_LabelsArray[20] = { { -1, -1, 0, -1, 1, -1, -1, -1 } }; + this->m_LabelsArray[21] = { { 0, -1, 0, -1, 0, -1, -1, -1 } }; + this->m_LabelsArray[22] = { { -1, 0, 1, -1, 2, -1, -1, -1 } }; + this->m_LabelsArray[23] = { { 0, 0, 0, -1, 0, -1, -1, -1 } }; + this->m_LabelsArray[24] = { { -1, -1, -1, 0, 1, -1, -1, -1 } }; + this->m_LabelsArray[25] = { { 0, -1, -1, 1, 0, -1, -1, -1 } }; + this->m_LabelsArray[26] = { { -1, 0, -1, 0, 1, -1, -1, -1 } }; + this->m_LabelsArray[27] = { { 0, 0, -1, 0, 0, -1, -1, -1 } }; + this->m_LabelsArray[28] = { { -1, -1, 0, 0, 1, -1, -1, -1 } }; + this->m_LabelsArray[29] = { { 0, -1, 0, 0, 0, -1, -1, -1 } }; + this->m_LabelsArray[30] = { { -1, 0, 0, 0, 1, -1, -1, -1 } }; + this->m_LabelsArray[31] = { { 0, 0, 0, 0, 0, -1, -1, -1 } }; + this->m_LabelsArray[32] = { { -1, -1, -1, -1, -1, 0, -1, -1 } }; + this->m_LabelsArray[33] = { { 0, -1, -1, -1, -1, 1, -1, -1 } }; + this->m_LabelsArray[34] = { { -1, 0, -1, -1, -1, 0, -1, -1 } }; + this->m_LabelsArray[35] = { { 0, 0, -1, -1, -1, 0, -1, -1 } }; + this->m_LabelsArray[36] = { { -1, -1, 0, -1, -1, 1, -1, -1 } }; + this->m_LabelsArray[37] = { { 0, -1, 0, -1, -1, 1, -1, -1 } }; + this->m_LabelsArray[38] = { { -1, 0, 1, -1, -1, 0, -1, -1 } }; + this->m_LabelsArray[39] = { { 0, 0, 0, -1, -1, 0, -1, -1 } }; + this->m_LabelsArray[40] = { { -1, -1, -1, 0, -1, 1, -1, -1 } }; + this->m_LabelsArray[41] = { { 0, -1, -1, 1, -1, 2, -1, -1 } }; + this->m_LabelsArray[42] = { { -1, 0, -1, 0, -1, 0, -1, -1 } }; + this->m_LabelsArray[43] = { { 0, 0, -1, 0, -1, 0, -1, -1 } }; + this->m_LabelsArray[44] = { { -1, -1, 0, 0, -1, 1, -1, -1 } }; + this->m_LabelsArray[45] = { { 0, -1, 0, 0, -1, 1, -1, -1 } }; + this->m_LabelsArray[46] = { { -1, 0, 0, 0, -1, 0, -1, -1 } }; + this->m_LabelsArray[47] = { { 0, 0, 0, 0, -1, 0, -1, -1 } }; + this->m_LabelsArray[48] = { { -1, -1, -1, -1, 0, 0, -1, -1 } }; + this->m_LabelsArray[49] = { { 0, -1, -1, -1, 0, 0, -1, -1 } }; + this->m_LabelsArray[50] = { { -1, 0, -1, -1, 0, 0, -1, -1 } }; + this->m_LabelsArray[51] = { { 0, 0, -1, -1, 0, 0, -1, -1 } }; + this->m_LabelsArray[52] = { { -1, -1, 0, -1, 1, 1, -1, -1 } }; + this->m_LabelsArray[53] = { { 0, -1, 0, -1, 0, 0, -1, -1 } }; + this->m_LabelsArray[54] = { { -1, 0, 1, -1, 0, 0, -1, -1 } }; + this->m_LabelsArray[55] = { { 0, 0, 0, -1, 0, 0, -1, -1 } }; + this->m_LabelsArray[56] = { { -1, -1, -1, 0, 1, 1, -1, -1 } }; + this->m_LabelsArray[57] = { { 0, -1, -1, 1, 0, 0, -1, -1 } }; + this->m_LabelsArray[58] = { { -1, 0, -1, 0, 0, 0, -1, -1 } }; + this->m_LabelsArray[59] = { { 0, 0, -1, 0, 0, 0, -1, -1 } }; + this->m_LabelsArray[60] = { { -1, -1, 0, 0, 1, 1, -1, -1 } }; + this->m_LabelsArray[61] = { { 0, -1, 0, 0, 0, 0, -1, -1 } }; + this->m_LabelsArray[62] = { { -1, 0, 0, 0, 0, 0, -1, -1 } }; + this->m_LabelsArray[63] = { { 0, 0, 0, 0, 0, 0, -1, -1 } }; + this->m_LabelsArray[64] = { { -1, -1, -1, -1, -1, -1, 0, -1 } }; + this->m_LabelsArray[65] = { { 0, -1, -1, -1, -1, -1, 1, -1 } }; + this->m_LabelsArray[66] = { { -1, 0, -1, -1, -1, -1, 1, -1 } }; + this->m_LabelsArray[67] = { { 0, 0, -1, -1, -1, -1, 1, -1 } }; + this->m_LabelsArray[68] = { { -1, -1, 0, -1, -1, -1, 0, -1 } }; + this->m_LabelsArray[69] = { { 0, -1, 0, -1, -1, -1, 0, -1 } }; + this->m_LabelsArray[70] = { { -1, 0, 1, -1, -1, -1, 1, -1 } }; + this->m_LabelsArray[71] = { { 0, 0, 0, -1, -1, -1, 0, -1 } }; + this->m_LabelsArray[72] = { { -1, -1, -1, 0, -1, -1, 1, -1 } }; + this->m_LabelsArray[73] = { { 0, -1, -1, 1, -1, -1, 2, -1 } }; + this->m_LabelsArray[74] = { { -1, 0, -1, 0, -1, -1, 1, -1 } }; + this->m_LabelsArray[75] = { { 0, 0, -1, 0, -1, -1, 1, -1 } }; + this->m_LabelsArray[76] = { { -1, -1, 0, 0, -1, -1, 0, -1 } }; + this->m_LabelsArray[77] = { { 0, -1, 0, 0, -1, -1, 0, -1 } }; + this->m_LabelsArray[78] = { { -1, 0, 0, 0, -1, -1, 0, -1 } }; + this->m_LabelsArray[79] = { { 0, 0, 0, 0, -1, -1, 0, -1 } }; + this->m_LabelsArray[80] = { { -1, -1, -1, -1, 0, -1, 0, -1 } }; + this->m_LabelsArray[81] = { { 0, -1, -1, -1, 0, -1, 0, -1 } }; + this->m_LabelsArray[82] = { { -1, 0, -1, -1, 1, -1, 1, -1 } }; + this->m_LabelsArray[83] = { { 0, 0, -1, -1, 0, -1, 0, -1 } }; + this->m_LabelsArray[84] = { { -1, -1, 0, -1, 0, -1, 0, -1 } }; + this->m_LabelsArray[85] = { { 0, -1, 0, -1, 0, -1, 0, -1 } }; + this->m_LabelsArray[86] = { { -1, 0, 1, -1, 1, -1, 1, -1 } }; + this->m_LabelsArray[87] = { { 0, 0, 0, -1, 0, -1, 0, -1 } }; + this->m_LabelsArray[88] = { { -1, -1, -1, 0, 1, -1, 1, -1 } }; + this->m_LabelsArray[89] = { { 0, -1, -1, 1, 0, -1, 0, -1 } }; + this->m_LabelsArray[90] = { { -1, 0, -1, 0, 1, -1, 1, -1 } }; + this->m_LabelsArray[91] = { { 0, 0, -1, 0, 0, -1, 0, -1 } }; + this->m_LabelsArray[92] = { { -1, -1, 0, 0, 0, -1, 0, -1 } }; + this->m_LabelsArray[93] = { { 0, -1, 0, 0, 0, -1, 0, -1 } }; + this->m_LabelsArray[94] = { { -1, 0, 0, 0, 0, -1, 0, -1 } }; + this->m_LabelsArray[95] = { { 0, 0, 0, 0, 0, -1, 0, -1 } }; + this->m_LabelsArray[96] = { { -1, -1, -1, -1, -1, 0, 1, -1 } }; + this->m_LabelsArray[97] = { { 0, -1, -1, -1, -1, 1, 2, -1 } }; + this->m_LabelsArray[98] = { { -1, 0, -1, -1, -1, 0, 1, -1 } }; + this->m_LabelsArray[99] = { { 0, 0, -1, -1, -1, 0, 1, -1 } }; + this->m_LabelsArray[100] = { { -1, -1, 0, -1, -1, 1, 0, -1 } }; + this->m_LabelsArray[101] = { { 0, -1, 0, -1, -1, 1, 0, -1 } }; + this->m_LabelsArray[102] = { { -1, 0, 1, -1, -1, 0, 1, -1 } }; + this->m_LabelsArray[103] = { { 0, 0, 0, -1, -1, 0, 0, -1 } }; + this->m_LabelsArray[104] = { { -1, -1, -1, 0, -1, 1, 2, -1 } }; + this->m_LabelsArray[105] = { { 0, -1, -1, 1, -1, 2, 3, -1 } }; + this->m_LabelsArray[106] = { { -1, 0, -1, 0, -1, 0, 1, -1 } }; + this->m_LabelsArray[107] = { { 0, 0, -1, 0, -1, 0, 1, -1 } }; + this->m_LabelsArray[108] = { { -1, -1, 0, 0, -1, 1, 0, -1 } }; + this->m_LabelsArray[109] = { { 0, -1, 0, 0, -1, 1, 0, -1 } }; + this->m_LabelsArray[110] = { { -1, 0, 0, 0, -1, 0, 0, -1 } }; + this->m_LabelsArray[111] = { { 0, 0, 0, 0, -1, 0, 0, -1 } }; + this->m_LabelsArray[112] = { { -1, -1, -1, -1, 0, 0, 0, -1 } }; + this->m_LabelsArray[113] = { { 0, -1, -1, -1, 0, 0, 0, -1 } }; + this->m_LabelsArray[114] = { { -1, 0, -1, -1, 0, 0, 0, -1 } }; + this->m_LabelsArray[115] = { { 0, 0, -1, -1, 0, 0, 0, -1 } }; + this->m_LabelsArray[116] = { { -1, -1, 0, -1, 0, 0, 0, -1 } }; + this->m_LabelsArray[117] = { { 0, -1, 0, -1, 0, 0, 0, -1 } }; + this->m_LabelsArray[118] = { { -1, 0, 0, -1, 0, 0, 0, -1 } }; + this->m_LabelsArray[119] = { { 0, 0, 0, -1, 0, 0, 0, -1 } }; + this->m_LabelsArray[120] = { { -1, -1, -1, 0, 1, 1, 1, -1 } }; + this->m_LabelsArray[121] = { { 0, -1, -1, 1, 0, 0, 0, -1 } }; + this->m_LabelsArray[122] = { { -1, 0, -1, 0, 0, 0, 0, -1 } }; + this->m_LabelsArray[123] = { { 0, 0, -1, 0, 0, 0, 0, -1 } }; + this->m_LabelsArray[124] = { { -1, -1, 0, 0, 0, 0, 0, -1 } }; + this->m_LabelsArray[125] = { { 0, -1, 0, 0, 0, 0, 0, -1 } }; + this->m_LabelsArray[126] = { { -1, 0, 0, 1, 0, 1, 1, -1 } }; + this->m_LabelsArray[127] = { { 0, 0, 0, 0, 0, 0, 0, -1 } }; + this->m_LabelsArray[128] = { { -1, -1, -1, -1, -1, -1, -1, 0 } }; + this->m_LabelsArray[129] = { { 0, -1, -1, -1, -1, -1, -1, 1 } }; + this->m_LabelsArray[130] = { { -1, 0, -1, -1, -1, -1, -1, 1 } }; + this->m_LabelsArray[131] = { { 0, 0, -1, -1, -1, -1, -1, 1 } }; + this->m_LabelsArray[132] = { { -1, -1, 0, -1, -1, -1, -1, 1 } }; + this->m_LabelsArray[133] = { { 0, -1, 0, -1, -1, -1, -1, 1 } }; + this->m_LabelsArray[134] = { { -1, 0, 1, -1, -1, -1, -1, 2 } }; + this->m_LabelsArray[135] = { { 0, 0, 0, -1, -1, -1, -1, 1 } }; + this->m_LabelsArray[136] = { { -1, -1, -1, 0, -1, -1, -1, 0 } }; + this->m_LabelsArray[137] = { { 0, -1, -1, 1, -1, -1, -1, 1 } }; + this->m_LabelsArray[138] = { { -1, 0, -1, 0, -1, -1, -1, 0 } }; + this->m_LabelsArray[139] = { { 0, 0, -1, 0, -1, -1, -1, 0 } }; + this->m_LabelsArray[140] = { { -1, -1, 0, 0, -1, -1, -1, 0 } }; + this->m_LabelsArray[141] = { { 0, -1, 0, 0, -1, -1, -1, 0 } }; + this->m_LabelsArray[142] = { { -1, 0, 0, 0, -1, -1, -1, 0 } }; + this->m_LabelsArray[143] = { { 0, 0, 0, 0, -1, -1, -1, 0 } }; + this->m_LabelsArray[144] = { { -1, -1, -1, -1, 0, -1, -1, 1 } }; + this->m_LabelsArray[145] = { { 0, -1, -1, -1, 0, -1, -1, 1 } }; + this->m_LabelsArray[146] = { { -1, 0, -1, -1, 1, -1, -1, 2 } }; + this->m_LabelsArray[147] = { { 0, 0, -1, -1, 0, -1, -1, 1 } }; + this->m_LabelsArray[148] = { { -1, -1, 0, -1, 1, -1, -1, 2 } }; + this->m_LabelsArray[149] = { { 0, -1, 0, -1, 0, -1, -1, 1 } }; + this->m_LabelsArray[150] = { { -1, 0, 1, -1, 2, -1, -1, 3 } }; + this->m_LabelsArray[151] = { { 0, 0, 0, -1, 0, -1, -1, 1 } }; + this->m_LabelsArray[152] = { { -1, -1, -1, 0, 1, -1, -1, 0 } }; + this->m_LabelsArray[153] = { { 0, -1, -1, 1, 0, -1, -1, 1 } }; + this->m_LabelsArray[154] = { { -1, 0, -1, 0, 1, -1, -1, 0 } }; + this->m_LabelsArray[155] = { { 0, 0, -1, 0, 0, -1, -1, 0 } }; + this->m_LabelsArray[156] = { { -1, -1, 0, 0, 1, -1, -1, 0 } }; + this->m_LabelsArray[157] = { { 0, -1, 0, 0, 0, -1, -1, 0 } }; + this->m_LabelsArray[158] = { { -1, 0, 0, 0, 1, -1, -1, 0 } }; + this->m_LabelsArray[159] = { { 0, 0, 0, 0, 0, -1, -1, 0 } }; + this->m_LabelsArray[160] = { { -1, -1, -1, -1, -1, 0, -1, 0 } }; + this->m_LabelsArray[161] = { { 0, -1, -1, -1, -1, 1, -1, 1 } }; + this->m_LabelsArray[162] = { { -1, 0, -1, -1, -1, 0, -1, 0 } }; + this->m_LabelsArray[163] = { { 0, 0, -1, -1, -1, 0, -1, 0 } }; + this->m_LabelsArray[164] = { { -1, -1, 0, -1, -1, 1, -1, 1 } }; + this->m_LabelsArray[165] = { { 0, -1, 0, -1, -1, 1, -1, 1 } }; + this->m_LabelsArray[166] = { { -1, 0, 1, -1, -1, 0, -1, 0 } }; + this->m_LabelsArray[167] = { { 0, 0, 0, -1, -1, 0, -1, 0 } }; + this->m_LabelsArray[168] = { { -1, -1, -1, 0, -1, 0, -1, 0 } }; + this->m_LabelsArray[169] = { { 0, -1, -1, 1, -1, 1, -1, 1 } }; + this->m_LabelsArray[170] = { { -1, 0, -1, 0, -1, 0, -1, 0 } }; + this->m_LabelsArray[171] = { { 0, 0, -1, 0, -1, 0, -1, 0 } }; + this->m_LabelsArray[172] = { { -1, -1, 0, 0, -1, 0, -1, 0 } }; + this->m_LabelsArray[173] = { { 0, -1, 0, 0, -1, 0, -1, 0 } }; + this->m_LabelsArray[174] = { { -1, 0, 0, 0, -1, 0, -1, 0 } }; + this->m_LabelsArray[175] = { { 0, 0, 0, 0, -1, 0, -1, 0 } }; + this->m_LabelsArray[176] = { { -1, -1, -1, -1, 0, 0, -1, 0 } }; + this->m_LabelsArray[177] = { { 0, -1, -1, -1, 0, 0, -1, 0 } }; + this->m_LabelsArray[178] = { { -1, 0, -1, -1, 0, 0, -1, 0 } }; + this->m_LabelsArray[179] = { { 0, 0, -1, -1, 0, 0, -1, 0 } }; + this->m_LabelsArray[180] = { { -1, -1, 0, -1, 1, 1, -1, 1 } }; + this->m_LabelsArray[181] = { { 0, -1, 0, -1, 0, 0, -1, 0 } }; + this->m_LabelsArray[182] = { { -1, 0, 1, -1, 0, 0, -1, 0 } }; + this->m_LabelsArray[183] = { { 0, 0, 0, -1, 0, 0, -1, 0 } }; + this->m_LabelsArray[184] = { { -1, -1, -1, 0, 0, 0, -1, 0 } }; + this->m_LabelsArray[185] = { { 0, -1, -1, 0, 0, 0, -1, 0 } }; + this->m_LabelsArray[186] = { { -1, 0, -1, 0, 0, 0, -1, 0 } }; + this->m_LabelsArray[187] = { { 0, 0, -1, 0, 0, 0, -1, 0 } }; + this->m_LabelsArray[188] = { { -1, -1, 0, 0, 0, 0, -1, 0 } }; + this->m_LabelsArray[189] = { { 0, -1, 1, 0, 1, 0, -1, 1 } }; + this->m_LabelsArray[190] = { { -1, 0, 0, 0, 0, 0, -1, 0 } }; + this->m_LabelsArray[191] = { { 0, 0, 0, 0, 0, 0, -1, 0 } }; + this->m_LabelsArray[192] = { { -1, -1, -1, -1, -1, -1, 0, 0 } }; + this->m_LabelsArray[193] = { { 0, -1, -1, -1, -1, -1, 1, 1 } }; + this->m_LabelsArray[194] = { { -1, 0, -1, -1, -1, -1, 1, 1 } }; + this->m_LabelsArray[195] = { { 0, 0, -1, -1, -1, -1, 1, 1 } }; + this->m_LabelsArray[196] = { { -1, -1, 0, -1, -1, -1, 0, 0 } }; + this->m_LabelsArray[197] = { { 0, -1, 0, -1, -1, -1, 0, 0 } }; + this->m_LabelsArray[198] = { { -1, 0, 1, -1, -1, -1, 1, 1 } }; + this->m_LabelsArray[199] = { { 0, 0, 0, -1, -1, -1, 0, 0 } }; + this->m_LabelsArray[200] = { { -1, -1, -1, 0, -1, -1, 0, 0 } }; + this->m_LabelsArray[201] = { { 0, -1, -1, 1, -1, -1, 1, 1 } }; + this->m_LabelsArray[202] = { { -1, 0, -1, 0, -1, -1, 0, 0 } }; + this->m_LabelsArray[203] = { { 0, 0, -1, 0, -1, -1, 0, 0 } }; + this->m_LabelsArray[204] = { { -1, -1, 0, 0, -1, -1, 0, 0 } }; + this->m_LabelsArray[205] = { { 0, -1, 0, 0, -1, -1, 0, 0 } }; + this->m_LabelsArray[206] = { { -1, 0, 0, 0, -1, -1, 0, 0 } }; + this->m_LabelsArray[207] = { { 0, 0, 0, 0, -1, -1, 0, 0 } }; + this->m_LabelsArray[208] = { { -1, -1, -1, -1, 0, -1, 0, 0 } }; + this->m_LabelsArray[209] = { { 0, -1, -1, -1, 0, -1, 0, 0 } }; + this->m_LabelsArray[210] = { { -1, 0, -1, -1, 1, -1, 1, 1 } }; + this->m_LabelsArray[211] = { { 0, 0, -1, -1, 0, -1, 0, 0 } }; + this->m_LabelsArray[212] = { { -1, -1, 0, -1, 0, -1, 0, 0 } }; + this->m_LabelsArray[213] = { { 0, -1, 0, -1, 0, -1, 0, 0 } }; + this->m_LabelsArray[214] = { { -1, 0, 1, -1, 1, -1, 1, 1 } }; + this->m_LabelsArray[215] = { { 0, 0, 0, -1, 0, -1, 0, 0 } }; + this->m_LabelsArray[216] = { { -1, -1, -1, 0, 0, -1, 0, 0 } }; + this->m_LabelsArray[217] = { { 0, -1, -1, 0, 0, -1, 0, 0 } }; + this->m_LabelsArray[218] = { { -1, 0, -1, 0, 0, -1, 0, 0 } }; + this->m_LabelsArray[219] = { { 0, 1, -1, 0, 1, -1, 0, 1 } }; + this->m_LabelsArray[220] = { { -1, -1, 0, 0, 0, -1, 0, 0 } }; + this->m_LabelsArray[221] = { { 0, -1, 0, 0, 0, -1, 0, 0 } }; + this->m_LabelsArray[222] = { { -1, 0, 0, 0, 0, -1, 0, 0 } }; + this->m_LabelsArray[223] = { { 0, 0, 0, 0, 0, -1, 0, 0 } }; + this->m_LabelsArray[224] = { { -1, -1, -1, -1, -1, 0, 0, 0 } }; + this->m_LabelsArray[225] = { { 0, -1, -1, -1, -1, 1, 1, 1 } }; + this->m_LabelsArray[226] = { { -1, 0, -1, -1, -1, 0, 0, 0 } }; + this->m_LabelsArray[227] = { { 0, 0, -1, -1, -1, 0, 0, 0 } }; + this->m_LabelsArray[228] = { { -1, -1, 0, -1, -1, 0, 0, 0 } }; + this->m_LabelsArray[229] = { { 0, -1, 0, -1, -1, 0, 0, 0 } }; + this->m_LabelsArray[230] = { { -1, 0, 0, -1, -1, 0, 0, 0 } }; + this->m_LabelsArray[231] = { { 0, 1, 1, -1, -1, 0, 0, 1 } }; + this->m_LabelsArray[232] = { { -1, -1, -1, 0, -1, 0, 0, 0 } }; + this->m_LabelsArray[233] = { { 0, -1, -1, 1, -1, 1, 1, 1 } }; + this->m_LabelsArray[234] = { { -1, 0, -1, 0, -1, 0, 0, 0 } }; + this->m_LabelsArray[235] = { { 0, 0, -1, 0, -1, 0, 0, 0 } }; + this->m_LabelsArray[236] = { { -1, -1, 0, 0, -1, 0, 0, 0 } }; + this->m_LabelsArray[237] = { { 0, -1, 0, 0, -1, 0, 0, 0 } }; + this->m_LabelsArray[238] = { { -1, 0, 0, 0, -1, 0, 0, 0 } }; + this->m_LabelsArray[239] = { { 0, 0, 0, 0, -1, 0, 0, 0 } }; + this->m_LabelsArray[240] = { { -1, -1, -1, -1, 0, 0, 0, 0 } }; + this->m_LabelsArray[241] = { { 0, -1, -1, -1, 0, 0, 0, 0 } }; + this->m_LabelsArray[242] = { { -1, 0, -1, -1, 0, 0, 0, 0 } }; + this->m_LabelsArray[243] = { { 0, 0, -1, -1, 0, 0, 0, 0 } }; + this->m_LabelsArray[244] = { { -1, -1, 0, -1, 0, 0, 0, 0 } }; + this->m_LabelsArray[245] = { { 0, -1, 0, -1, 0, 0, 0, 0 } }; + this->m_LabelsArray[246] = { { -1, 0, 0, -1, 0, 0, 0, 0 } }; + this->m_LabelsArray[247] = { { 0, 0, 0, -1, 0, 0, 0, 0 } }; + this->m_LabelsArray[248] = { { -1, -1, -1, 0, 0, 0, 0, 0 } }; + this->m_LabelsArray[249] = { { 0, -1, -1, 0, 0, 0, 0, 0 } }; + this->m_LabelsArray[250] = { { -1, 0, -1, 0, 0, 0, 0, 0 } }; + this->m_LabelsArray[251] = { { 0, 0, -1, 0, 0, 0, 0, 0 } }; + this->m_LabelsArray[252] = { { -1, -1, 0, 0, 0, 0, 0, 0 } }; + this->m_LabelsArray[253] = { { 0, -1, 0, 0, 0, 0, 0, 0 } }; + this->m_LabelsArray[254] = { { -1, 0, 0, 0, 0, 0, 0, 0 } }; + this->m_LabelsArray[255] = { { 0, 0, 0, 0, 0, 0, 0, 0 } }; +} + +template +void +CuberilleImageToMeshFilter::PrintSelf(std::ostream & os, Indent indent) const +{ + Superclass::PrintSelf(os, indent); + + os << indent + << "IsoSurfaceValue: " << static_cast::PrintType>(m_IsoSurfaceValue) + << std::endl; + os << indent << "MaxSpacing: " << static_cast::PrintType>(m_MaxSpacing) + << std::endl; + os << indent << "GenerateTriangleFaces: " << static_cast::PrintType>(m_GenerateTriangleFaces) + << std::endl; + os << indent + << "ProjectVerticesToIsoSurface: " << static_cast::PrintType>(m_ProjectVerticesToIsoSurface) + << std::endl; + os << indent << "ProjectVertexSurfaceDistanceThreshold: " << m_ProjectVertexSurfaceDistanceThreshold << std::endl; + os << indent << "ProjectVertexStepLength: " << m_ProjectVertexStepLength << std::endl; + os << indent << "ProjectVertexStepLengthRelaxationFactor: " << m_ProjectVertexStepLengthRelaxationFactor << std::endl; + os << indent << "ProjectVertexMaximumNumberOfSteps: " << m_ProjectVertexMaximumNumberOfSteps << std::endl; +} +} // end namespace itk +#endif diff --git a/Modules/Filtering/Cuberille/itk-module.cmake b/Modules/Filtering/Cuberille/itk-module.cmake new file mode 100644 index 000000000000..56a9b1eb2b95 --- /dev/null +++ b/Modules/Filtering/Cuberille/itk-module.cmake @@ -0,0 +1,41 @@ +set( + DOCUMENTATION + "This module implements cuberille implicit surface +polygonization for ITK. This method operates by diving the surface into a +number of small cubes called cuberilles. Each cuberille is centered at a +pixel lying on the iso-surface and then quadrilaterals are generated for each +face. The original approach is improved by projecting the vertices of each +cuberille onto the implicit surface, smoothing the typical block-like +resultant mesh. + +A more detailed description can be found in the Insight Journal article: + + Mueller, D. \"Cuberille Implicit Surface Polygonization for ITK\" + https://hdl.handle.net/10380/3186 + https://www.insight-journal.org/browse/publication/740 + July 20, 2010. +" +) + +itk_module( + Cuberille + DEPENDS + ITKCommon + ITKImageFunction + ITKImageGradient + ITKMesh + ITKConnectedComponents + TEST_DEPENDS + ITKTestKernel + ITKQuadEdgeMesh + ITKQuadEdgeMeshFiltering + ITKThresholding + ITKIOImageBase + ITKIOMeta + ITKIONRRD + ITKIOMeshBase + ITKIOVTK + DESCRIPTION "${DOCUMENTATION}" + EXCLUDE_FROM_DEFAULT + ENABLE_SHARED +) diff --git a/Modules/Filtering/Cuberille/test/CMakeLists.txt b/Modules/Filtering/Cuberille/test/CMakeLists.txt new file mode 100644 index 000000000000..a8373c9f7478 --- /dev/null +++ b/Modules/Filtering/Cuberille/test/CMakeLists.txt @@ -0,0 +1,420 @@ +itk_module_test() + +set(DATA_DIR ${CMAKE_CURRENT_SOURCE_DIR}/Input) + +set( + CuberilleTests + CuberilleTest01.cxx + CuberilleTest02.cxx + CuberilleTest03.cxx + CuberilleTest04.cxx + CuberilleTest_Issue66.cxx +) + +createtestdriver(Cuberille "${Cuberille-Test_LIBRARIES}" "${CuberilleTests}") + +# Configure the Tests +itk_add_test( + NAME Cuberille_Blob0_00 + COMMAND + CuberilleTestDriver + CuberilleTest01 + ${DATA_DIR}/blob0.mha + ${ITK_TEST_OUTPUT_DIR}/blob0-01.vtk + 200 # Iso-surface value + 8 # Expected number of points + 6 # Expected number of cells + 0 # Generate triangle faces + 0 # Project vertices to iso-surface + 0.2 # Surface distance threshold + 0.24 # Step length + 0.95 # Step length relaxation factor + 100 # Maximum number of steps +) + +itk_add_test( + NAME Cuberille_Blob1_01 + COMMAND + CuberilleTestDriver + CuberilleTest01 + ${DATA_DIR}/blob1.mha + ${ITK_TEST_OUTPUT_DIR}/blob1-01.vtk + 200 # Iso-surface value + 12 # Expected number of points + 10 # Expected number of cells + 0 # Generate triangle faces + 0 # Project vertices to iso-surface + 0.2 # Surface distance threshold + 0.24 # Step length + 0.95 # Step length relaxation factor + 100 # Maximum number of steps +) + +itk_add_test( + NAME Cuberille_Blob2_01 + COMMAND + CuberilleTestDriver + CuberilleTest01 + ${DATA_DIR}/blob2.mha + ${ITK_TEST_OUTPUT_DIR}/blob2-01.vtk + 200 # Iso-surface value + 16 # Expected number of points + 12 # Expected number of cells + 0 # Generate triangle faces + 0 # Project vertices to iso-surface + 0.2 # Surface distance threshold + 0.24 # Step length + 0.95 # Step length relaxation factor + 100 # Maximum number of steps +) + +itk_add_test( + NAME Cuberille_Blob3_01 + COMMAND + CuberilleTestDriver + CuberilleTest01 + ${DATA_DIR}/blob3.mha + ${ITK_TEST_OUTPUT_DIR}/blob3-01.vtk + 200 # Iso-surface value + 240 # Expected number of points + 180 # Expected number of cells + 0 # Generate triangle faces + 0 # Project vertices to iso-surface + 0.2 # Surface distance threshold + 0.24 # Step length + 0.95 # Step length relaxation factor + 100 # Maximum number of steps +) + +itk_add_test( + NAME Cuberille_Blob4_01 + COMMAND + CuberilleTestDriver + CuberilleTest01 + ${DATA_DIR}/blob4.mha + ${ITK_TEST_OUTPUT_DIR}/blob4-01.vtk + 200 # Iso-surface value + 2124 # Expected number of points + 2122 # Expected number of cells + 0 # Generate triangle faces + 1 # Project vertices to iso-surface + 0.2 # Surface distance threshold + 0.24 # Step length + 0.95 # Step length relaxation factor + 100 # Maximum number of steps +) + +itk_add_test( + NAME Cuberille_MarschnerLobb_01 + COMMAND + CuberilleTestDriver + CuberilleTest01 + ${DATA_DIR}/marschnerlobb.mha + ${ITK_TEST_OUTPUT_DIR}/marschnerlobb-01.vtk + 55 # Iso-surface value + 22106 # Expected number of points + 22104 # Expected number of cells + 0 # Generate triangle faces + 1 # Project vertices to iso-surface + 0.2 # Surface distance threshold + 0.24 # Step length + 0.95 # Step length relaxation factor + 200 # Maximum number of steps +) + +itk_add_test( + NAME Cuberille_Fuel_01 + COMMAND + CuberilleTestDriver + CuberilleTest01 + ${DATA_DIR}/fuel.mha + ${ITK_TEST_OUTPUT_DIR}/fuel-01.vtk + 15 # Iso-surface value + 5302 # Expected number of points + 5316 # Expected number of cells + 0 # Generate triangle faces + 0 # Project vertices to iso-surface + 0.2 # Surface distance threshold + 0.24 # Step length + 0.95 # Step length relaxation factor + 100 # Maximum number of steps +) + +itk_add_test( + NAME Cuberille_Fuel_02 + COMMAND + CuberilleTestDriver + CuberilleTest01 + ${DATA_DIR}/fuel.mha + ${ITK_TEST_OUTPUT_DIR}/fuel-02.vtk + 15 # Iso-surface value + 5302 # Expected number of points + 5316 # Expected number of cells + 0 # Generate triangle faces + 1 # Project vertices to iso-surface + 0.2 # Surface distance threshold + 0.24 # Step length + 0.95 # Step length relaxation factor + 100 # Maximum number of steps +) + +itk_add_test( + NAME Cuberille_Fuel_03 + COMMAND + CuberilleTestDriver + CuberilleTest01 + ${DATA_DIR}/fuel.mha + ${ITK_TEST_OUTPUT_DIR}/fuel-03.vtk + 15 # Iso-surface value + 5302 # Expected number of points + 10632 # Expected number of cells + 1 # Generate triangle faces + 1 # Project vertices to iso-surface + 0.2 # Surface distance threshold + 0.24 # Step length + 0.95 # Step length relaxation factor + 100 # Maximum number of steps +) + +itk_add_test( + NAME Cuberille_HydrogenAtom_01 + COMMAND + CuberilleTestDriver + CuberilleTest01 + ${DATA_DIR}/hydrogenAtom.mha + ${ITK_TEST_OUTPUT_DIR}/hydrogenAtom-01.vtk + 15 # Iso-surface value + 29880 # Expected number of points + 29874 # Expected number of cells + 0 # Generate triangle faces + 1 # Project vertices to iso-surface + 0.2 # Surface distance threshold + 0.24 # Step length + 0.95 # Step length relaxation factor + 100 # Maximum number of steps +) + +itk_add_test( + NAME Cuberille_Neghip_01 + COMMAND + CuberilleTestDriver + CuberilleTest01 + ${DATA_DIR}/neghip.mha + ${ITK_TEST_OUTPUT_DIR}/neghip-01.vtk + 55 # Iso-surface value + 15164 # Expected number of points + 15136 # Expected number of cells + 0 # Generate triangle faces + 0 # Project vertices to iso-surface + 0.2 # Surface distance threshold + 0.24 # Step length + 0.95 # Step length relaxation factor + 100 # Maximum number of steps +) + +itk_add_test( + NAME Cuberille_Neghip_02 + COMMAND + CuberilleTestDriver + CuberilleTest01 + ${DATA_DIR}/neghip.mha + ${ITK_TEST_OUTPUT_DIR}/neghip-02.vtk + 55 # Iso-surface value + 15164 # Expected number of points + 15136 # Expected number of cells + 0 # Generate triangle faces + 1 # Project vertices to iso-surface + 0.2 # Surface distance threshold + 0.24 # Step length + 0.95 # Step length relaxation factor + 100 # Maximum number of steps +) + +itk_add_test( + NAME Cuberille_Neghip_03 + COMMAND + CuberilleTestDriver + CuberilleTest01 + ${DATA_DIR}/neghip.mha + ${ITK_TEST_OUTPUT_DIR}/neghip-03.vtk + 55 # Iso-surface value + 15164 # Expected number of points + 30272 # Expected number of cells + 1 # Generate triangle faces + 1 # Project vertices to iso-surface + 0.2 # Surface distance threshold + 0.24 # Step length + 0.95 # Step length relaxation factor + 100 # Maximum number of steps +) + +itk_add_test( + NAME Cuberille_Nucleon_01 + COMMAND + CuberilleTestDriver + CuberilleTest01 + ${DATA_DIR}/nucleon.mha + ${ITK_TEST_OUTPUT_DIR}/nucleon-01.vtk + 140 # Iso-surface value + 3504 # Expected number of points + 3500 # Expected number of cells + 0 # Generate triangle faces + 0 # Project vertices to iso-surface + 0.2 # Surface distance threshold + 0.24 # Step length + 0.95 # Step length relaxation factor + 100 # Maximum number of steps +) + +itk_add_test( + NAME Cuberille_Nucleon_02 + COMMAND + CuberilleTestDriver + CuberilleTest01 + ${DATA_DIR}/nucleon.mha + ${ITK_TEST_OUTPUT_DIR}/nucleon-02.vtk + 140 # Iso-surface value + 3504 # Expected number of points + 3500 # Expected number of cells + 0 # Generate triangle faces + 1 # Project vertices to iso-surface + 0.2 # Surface distance threshold + 0.24 # Step length + 0.95 # Step length relaxation factor + 100 # Maximum number of steps +) + +itk_add_test( + NAME Cuberille_Nucleon_03 + COMMAND + CuberilleTestDriver + CuberilleTest01 + ${DATA_DIR}/nucleon.mha + ${ITK_TEST_OUTPUT_DIR}/nucleon-03.vtk + 140 # Iso-surface value + 3504 # Expected number of points + 7000 # Expected number of cells + 1 # Generate triangle faces + 1 # Project vertices to iso-surface + 0.2 # Surface distance threshold + 0.24 # Step length + 0.95 # Step length relaxation factor + 100 # Maximum number of steps +) + +itk_add_test( + NAME Cuberille_Silicium_01 + COMMAND + CuberilleTestDriver + CuberilleTest01 + ${DATA_DIR}/silicium.mha + ${ITK_TEST_OUTPUT_DIR}/silicium-01.vtk + 85 # Iso-surface value + 20036 # Expected number of points + 20024 # Expected number of cells + 0 # Generate triangle faces + 0 # Project vertices to iso-surface + 0.2 # Surface distance threshold + 0.24 # Step length + 0.95 # Step length relaxation factor + 100 # Maximum number of steps +) + +itk_add_test( + NAME Cuberille_Silicium_02 + COMMAND + CuberilleTestDriver + CuberilleTest01 + ${DATA_DIR}/silicium.mha + ${ITK_TEST_OUTPUT_DIR}/silicium-02.vtk + 85 # Iso-surface value + 20036 # Expected number of points + 20024 # Expected number of cells + 0 # Generate triangle faces + 1 # Project vertices to iso-surface + 0.2 # Surface distance threshold + 0.24 # Step length + 0.95 # Step length relaxation factor + 100 # Maximum number of steps +) + +itk_add_test( + NAME Cuberille_Silicium_03 + COMMAND + CuberilleTestDriver + CuberilleTest01 + ${DATA_DIR}/silicium.mha + ${ITK_TEST_OUTPUT_DIR}/silicium-03.vtk + 85 # Iso-surface value + 20036 # Expected number of points + 40048 # Expected number of cells + 1 # Generate triangle faces + 1 # Project vertices to iso-surface + 0.2 # Surface distance threshold + 0.24 # Step length + 0.95 # Step length relaxation factor + 100 # Maximum number of steps +) + +#add_test( +# Cuberille_Engine_01 +# CuberilleTest01 +# Test01 +# ${DATA_DIR}/engine2.mha +# ${ITK_TEST_OUTPUT_DIR}/engine2-01.vtk +# 115 # Iso-surface value +# 319918 # Expected number of points +# 640112 # Expected number of cells +# 0 # Generate triangle faces +# 1 # Project vertices to iso-surface +# 0.2 # Surface distance threshold +# 0.24 # Step length +# 0.95 # Step length relaxation factor +# 100 # Maximum number of steps +#) + +#add_test( +# Cuberille_Bunny_01 +# CuberilleTest01 +# Test01 +# ${DATA_DIR}/bunny3.mha +# ${ITK_TEST_OUTPUT_DIR}/bunny3-01.vtk +# 155 # unsigned char +# #1550 # signed short +# 1021438 # Expected number of points +# 2042892 # Expected number of cells +# 1 # Generate triangle faces +# 1 # Project vertices to iso-surface +# 0.10 # Surface distance threshold +# 0.15 # Step length +# 0.85 # Step length relaxation factor +# 25 # Maximum number of steps +#) + +itk_add_test( + NAME CuberilleIssue66 + COMMAND + ${itk-module}TestDriver + CuberilleTest_Issue66 +) + +itk_add_test( + NAME CuberilleTestDirectionMatrix + COMMAND + ${itk-module}TestDriver + CuberilleTest02 +) + +itk_add_test( + NAME CuberilleTestCellData + COMMAND + ${itk-module}TestDriver + CuberilleTest03 +) + +itk_add_test( + NAME CuberilleTestNonManifoldGeometry + COMMAND + ${itk-module}TestDriver + CuberilleTest04 +) diff --git a/Modules/Filtering/Cuberille/test/CuberilleTest01.cxx b/Modules/Filtering/Cuberille/test/CuberilleTest01.cxx new file mode 100644 index 000000000000..78253d3cc380 --- /dev/null +++ b/Modules/Filtering/Cuberille/test/CuberilleTest01.cxx @@ -0,0 +1,286 @@ +/*========================================================================= + * + * 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. + * + *=========================================================================*/ + +#define USE_BSPLINE_INTERPOLATOR 0 +#define USE_MARCHING_CUBES 0 +#define USE_DECIMATION 0 + +#include "itkTimeProbe.h" +#include "itkImage.h" +#include "itkMesh.h" +#include "itkQuadEdgeMesh.h" +#include "itkCuberilleImageToMeshFilter.h" +#include "itkBinaryThresholdImageFilter.h" +#include "itkBinaryMask3DMeshSource.h" +#include "itkImageFileReader.h" +#include "itkImageFileWriter.h" +#include "itkVTKPolyDataWriter.h" +#include "itkLinearInterpolateImageFunction.h" +#include "itkBSplineInterpolateImageFunction.h" +#include "itkQuadricDecimationQuadEdgeMeshFilter.h" +#include "itkQuadEdgeMeshDecimationCriteria.h" +#include "itkTestingMacros.h" + +template +int +CuberilleTest01Helper(int argc, char * argv[]) +{ + + if (argc < 6) + { + std::cout << "Usage: " << argv[0]; + std::cout << "InputImage OutputMesh IsoSurfaceValue ExpectedNumberOfPoints ExpectedNumberOfCells"; + std::cout << "[GenerateTriangleFaces] [ProjectToIsoSurface] "; + std::cout << "[SurfaceDistanceThreshold] [StepLength] [StepLengthRelax] [MaximumNumberOfSteps]"; + std::cout << std::endl; + return EXIT_FAILURE; + } + + using PixelType = typename ImageType::PixelType; + using ImageFileReaderType = itk::ImageFileReader; + using MeshFileWriterType = itk::VTKPolyDataWriter; +#if USE_BSPLINE_INTERPOLATOR + using InterpolatorType = itk::BSplineInterpolateImageFunction; +#else + using InterpolatorType = itk::LinearInterpolateImageFunction; +#endif + using CuberilleType = itk::CuberilleImageToMeshFilter; + + + // Read command-line parameters + int arg = 1; + char * filenameInputImage = argv[arg++]; + char * filenameOutputMesh = argv[arg++]; + PixelType isoSurfaceValue = atoi(argv[arg++]); + unsigned int expectedNumberOfPoints = atoi(argv[arg++]); + unsigned int expectedNumberOfCells = atoi(argv[arg++]); + + bool generateTriangleFaces = true; + if (argc > arg) + { + generateTriangleFaces = atoi(argv[arg++]); + } + + bool projectToIsoSurface = true; + if (argc > arg) + { + projectToIsoSurface = atoi(argv[arg++]); + } + + double surfaceDistanceThreshold = 0.5; + if (argc > arg) + { + surfaceDistanceThreshold = atof(argv[arg++]); + } + + double stepLength = 0.25; + if (argc > arg) + { + stepLength = atof(argv[arg++]); + } + + double stepLengthRelax = 0.95; + if (argc > arg) + { + stepLengthRelax = atof(argv[arg++]); + } + + unsigned int maximumNumberOfSteps = 50; + if (argc > arg) + { + maximumNumberOfSteps = atoi(argv[arg++]); + } + + // Read input image + const auto reader = ImageFileReaderType::New(); + reader->SetFileName(filenameInputImage); + + ITK_TRY_EXPECT_NO_EXCEPTION(reader->UpdateLargestPossibleRegion()); + + typename ImageType::Pointer input = reader->GetOutput(); + input->DisconnectPipeline(); + + // Create output mesh + typename MeshType::Pointer outputMesh = nullptr; + itk::TimeProbe time; +#if USE_MARCHING_CUBES + + // Create marching cubes mesh + const auto threshold = BinaryThresholdFilterType::New(); + threshold->SetInput(input); + threshold->SetLowerThreshold(IsoSurfaceValue); + threshold->SetUpperThreshold(itk::NumericTraits::max()); + threshold->SetInsideValue(itk::NumericTraits::One); + threshold->SetOutsideValue(itk::NumericTraits::Zero); + threshold->UpdateLargestPossibleRegion(); + + const auto marching = MarchingCubesType::New(); + marching->SetInput(threshold->GetOutput()); + + time.Start(); + + ITK_TRY_EXPECT_NO_EXCEPTION(marching->Update()); + + time.Stop(); + + outputMesh = marching->GetOutput(); + outputMesh->DisconnectPipeline(); +#else + + // Create cuberille mesh filter + const auto cuberille = CuberilleType::New(); + + // How long does it take to pre-calculate the array labels array? + + itk::TimeProbe labelsArrayTimer; + + labelsArrayTimer.Start(); + cuberille->CalculateLabelsArray(); + labelsArrayTimer.Stop(); + + std::cout << "Time to calculate labels array: " << labelsArrayTimer.GetMean() << std::endl; + + cuberille->SetInput(input); + + cuberille->SetIsoSurfaceValue(isoSurfaceValue); + ITK_TEST_SET_GET_VALUE(isoSurfaceValue, cuberille->GetIsoSurfaceValue()); + + const auto interpolator = InterpolatorType::New(); +# if USE_BSPLINE_INTERPOLATOR + unsigned int splineOrder = 3; + interpolator->SetSplineOrder(splineOrder); +# endif + cuberille->SetInterpolator(interpolator); + ITK_TEST_SET_GET_VALUE(interpolator, cuberille->GetInterpolator()); + + cuberille->SetGenerateTriangleFaces(generateTriangleFaces); + ITK_TEST_SET_GET_BOOLEAN(cuberille, GenerateTriangleFaces, generateTriangleFaces); + + cuberille->SetProjectVerticesToIsoSurface(projectToIsoSurface); + ITK_TEST_SET_GET_BOOLEAN(cuberille, ProjectVerticesToIsoSurface, projectToIsoSurface); + + cuberille->SetProjectVertexSurfaceDistanceThreshold(surfaceDistanceThreshold); + ITK_TEST_SET_GET_VALUE(surfaceDistanceThreshold, cuberille->GetProjectVertexSurfaceDistanceThreshold()); + + cuberille->SetProjectVertexStepLength(stepLength); + ITK_TEST_SET_GET_VALUE(stepLength, cuberille->GetProjectVertexStepLength()); + + cuberille->SetProjectVertexStepLengthRelaxationFactor(stepLengthRelax); + ITK_TEST_SET_GET_VALUE(stepLengthRelax, cuberille->GetProjectVertexStepLengthRelaxationFactor()); + + cuberille->SetProjectVertexMaximumNumberOfSteps(maximumNumberOfSteps); + ITK_TEST_SET_GET_VALUE(maximumNumberOfSteps, cuberille->GetProjectVertexMaximumNumberOfSteps()); + + ITK_TEST_SET_GET_BOOLEAN(cuberille, SavePixelAsCellData, false); + + time.Start(); + + ITK_TRY_EXPECT_NO_EXCEPTION(cuberille->Update()); + + time.Stop(); + + outputMesh = cuberille->GetOutput(); + + outputMesh->DisconnectPipeline(); + +#endif + +#if USE_DECIMATION + // Decimation + using DecimationCriterionType = itk::NumberOfFacesCriterion; + const auto decimateCriterion = DecimationCriterionType::New(); + decimateCriterion->SetTopologicalChange(false); + decimateCriterion->SetNumberOfElements(2000); + + using DecimationType = itk::QuadricDecimationQuadEdgeMeshFilter; + const auto decimate = DecimationType::New(); + decimate->SetCriterion(decimateCriterion); + + decimate->SetInput(outputMesh); + + ITK_TRY_EXPECT_NO_EXCEPTION(decimate->Update()); +#endif + + // Write mesh + const auto writer = MeshFileWriterType::New(); +#if USE_DECIMATION + writer->SetInput(decimate->GetOutput()); +#else + writer->SetInput(outputMesh); +#endif + writer->SetFileName(filenameOutputMesh); + + ITK_TRY_EXPECT_NO_EXCEPTION(writer->Update()); + + // Assert number of points/cells + std::cout << "Polygonization took " << time.GetMean() << " seconds" << std::endl; + std::cout << "Mesh has " << outputMesh->GetNumberOfPoints() << " vertices "; + std::cout << "and " << outputMesh->GetNumberOfCells() << " cells" << std::endl; + if (expectedNumberOfPoints > 0 && outputMesh->GetNumberOfPoints() != expectedNumberOfPoints) + { + std::cerr << "Test failed!" << std::endl; + std::cerr << "Error: Expected mesh with " << expectedNumberOfPoints << " points, but found " + << outputMesh->GetNumberOfPoints() << std::endl; + return EXIT_FAILURE; + } + if (expectedNumberOfCells > 0 && outputMesh->GetNumberOfCells() != expectedNumberOfCells) + { + std::cerr << "Test failed!" << std::endl; + std::cerr << "Error: Expected mesh with " << expectedNumberOfCells << " cells, but found " + << outputMesh->GetNumberOfCells() << std::endl; + return EXIT_FAILURE; + } + + std::cout << "Test finished" << std::endl; + return EXIT_SUCCESS; +} + +int +CuberilleTest01(int argc, char * argv[]) +{ + + constexpr unsigned int Dimension = 3; + using PixelType = unsigned char; + using CoordinateType = double; + using QEMeshType = itk::QuadEdgeMesh; + using MeshType = itk::Mesh; + using ImageType = itk::Image; + + { + using CuberilleType = itk::CuberilleImageToMeshFilter; + const auto cuberille = CuberilleType::New(); + ITK_EXERCISE_BASIC_OBJECT_METHODS(cuberille, CuberilleImageToMeshFilter, ImageToMeshFilter); + } + + { + using CuberilleType = itk::CuberilleImageToMeshFilter; + const auto cuberille = CuberilleType::New(); + ITK_EXERCISE_BASIC_OBJECT_METHODS(cuberille, CuberilleImageToMeshFilter, ImageToMeshFilter); + } + + if (EXIT_FAILURE == CuberilleTest01Helper(argc, argv)) + { + return EXIT_FAILURE; + } + if (EXIT_FAILURE == CuberilleTest01Helper(argc, argv)) + { + return EXIT_FAILURE; + } + + return EXIT_SUCCESS; +} diff --git a/Modules/Filtering/Cuberille/test/CuberilleTest02.cxx b/Modules/Filtering/Cuberille/test/CuberilleTest02.cxx new file mode 100644 index 000000000000..b2f76f64d166 --- /dev/null +++ b/Modules/Filtering/Cuberille/test/CuberilleTest02.cxx @@ -0,0 +1,126 @@ +/*========================================================================= + * + * 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. + * + *=========================================================================*/ + +// STD +#include + +// ITK +#include "itkTestingMacros.h" +#include "itkImage.h" +#include "itkMesh.h" +#include "itkCuberilleImageToMeshFilter.h" +#include "itkNearestNeighborInterpolateImageFunction.h" + +const unsigned int Dimension = 3; +using TPixel = unsigned char; +using TCoordinate = double; +using TImage = itk::Image; +using TMesh = itk::Mesh; +using TImageToMesh = itk::CuberilleImageToMeshFilter; +using TInterp = itk::NearestNeighborInterpolateImageFunction; + +TImage::Pointer +CuberilleTestCreateImage(const std::array & flips) +{ + const auto image = TImage::New(); + + TImage::IndexType start; + start.Fill(0); + + TImage::SizeType size; + size.Fill(3); + + TImage::RegionType region; + region.SetSize(size); + region.SetIndex(start); + + image->SetRegions(region); + image->Allocate(); + image->FillBuffer(0); + + TImage::IndexType index; + index.Fill(1); + image->SetPixel(index, 1); + + TImage::DirectionType direction; + direction.SetIdentity(); + + for (size_t i = 0; i < flips.size(); ++i) + { + if (flips[i]) + { + direction(i, i) *= -1; + } + } + + image->SetDirection(direction); + + return image; +} + +void +CuberilleTestHelper(TImage::Pointer input) +{ + + const auto image_to_mesh = TImageToMesh::New(); + image_to_mesh->SetInput(input); + image_to_mesh->ProjectVerticesToIsoSurfaceOff(); + image_to_mesh->Update(); + + const auto mesh = TMesh::New(); + mesh->Graft(image_to_mesh->GetOutput()); + mesh->DisconnectPipeline(); + + TMesh::PointType center; + center.Fill(0.0); + + for (auto it = mesh->GetPoints()->Begin(); it != mesh->GetPoints()->End(); ++it) + { + center[0] += it.Value()[0]; + center[1] += it.Value()[1]; + center[2] += it.Value()[2]; + } + + center[0] /= 8.0; + center[1] /= 8.0; + center[2] /= 8.0; + + const auto interp = TInterp::New(); + interp->SetInputImage(input); + + if (1 != interp->Evaluate(center)) + { + throw 0; + } +} + +int +CuberilleTest02(int itkNotUsed(argc), char * itkNotUsed(argv)[]) +{ + + CuberilleTestHelper(CuberilleTestCreateImage({ { false, false, false } })); + CuberilleTestHelper(CuberilleTestCreateImage({ { false, false, true } })); + CuberilleTestHelper(CuberilleTestCreateImage({ { false, true, false } })); + CuberilleTestHelper(CuberilleTestCreateImage({ { false, true, true } })); + CuberilleTestHelper(CuberilleTestCreateImage({ { true, false, false } })); + CuberilleTestHelper(CuberilleTestCreateImage({ { true, false, true } })); + CuberilleTestHelper(CuberilleTestCreateImage({ { true, true, false } })); + CuberilleTestHelper(CuberilleTestCreateImage({ { true, true, true } })); + + return EXIT_SUCCESS; +} diff --git a/Modules/Filtering/Cuberille/test/CuberilleTest03.cxx b/Modules/Filtering/Cuberille/test/CuberilleTest03.cxx new file mode 100644 index 000000000000..0283983f5a1c --- /dev/null +++ b/Modules/Filtering/Cuberille/test/CuberilleTest03.cxx @@ -0,0 +1,177 @@ +/*========================================================================= + * + * 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. + * + *=========================================================================*/ + +// STD +#include + +// ITK +#include "itkTestingMacros.h" +#include "itkImage.h" +#include "itkMesh.h" +#include "itkCuberilleImageToMeshFilter.h" +#include "itkNearestNeighborInterpolateImageFunction.h" +#include "itkTriangleHelper.h" + +const unsigned int Dimension = 3; +using TPixel = unsigned char; +using TCoordinate = double; +using TImage = itk::Image; +using TMesh = itk::Mesh; +using TImageToMesh = itk::CuberilleImageToMeshFilter; +using TInterp = itk::NearestNeighborInterpolateImageFunction; +using TTriangleHelper = itk::TriangleHelper; + +TImage::Pointer +CuberilleTest03CreateImage() +{ + const auto image = TImage::New(); + + TImage::IndexType start; + start.Fill(0); + + TImage::SizeType size; + size[0] = 3; + size[1] = 5; + size[2] = 3; + + TImage::RegionType region; + region.SetSize(size); + region.SetIndex(start); + + image->SetRegions(region); + image->Allocate(); + image->FillBuffer(0); + + TImage::IndexType index; + index.Fill(1); + image->SetPixel(index, 1); + index[1] = 2; + image->SetPixel(index, 2); + index[1] = 3; + image->SetPixel(index, 3); + + return image; +} + +void +CuberilleTest03Helper(TImage::Pointer image) +{ + + const auto image_to_mesh = TImageToMesh::New(); + image_to_mesh->SetInput(image); + image_to_mesh->SavePixelAsCellDataOn(); + image_to_mesh->GenerateTriangleFacesOff(); + image_to_mesh->ProjectVerticesToIsoSurfaceOff(); + image_to_mesh->Update(); + + const auto mesh = TMesh::New(); + mesh->Graft(image_to_mesh->GetOutput()); + mesh->DisconnectPipeline(); + + const auto interp = TInterp::New(); + interp->SetInputImage(image); + + const auto half_spacing = image->GetSpacing() * 0.5; + + for (auto it = mesh->GetCells()->Begin(); it != mesh->GetCells()->End(); ++it) + { + const auto cell = it.Value(); + + typename TImage::PointType centroid; + centroid.SetToMidPoint(mesh->GetPoint(cell->GetPointIds()[0]), mesh->GetPoint(cell->GetPointIds()[2])); + + auto normal = TTriangleHelper::ComputeNormal(mesh->GetPoint(cell->GetPointIds()[0]), + mesh->GetPoint(cell->GetPointIds()[1]), + mesh->GetPoint(cell->GetPointIds()[2])); + + normal.Normalize(); + + for (size_t i = 0; i < 3; ++i) + { + normal[i] *= half_spacing[i]; + } + + const auto resample = centroid + -1.0f * normal; + const auto label = interp->Evaluate(resample); + + const auto data = mesh->GetCellData()->ElementAt(it.Index()); + + if (label != data) + { + std::cerr << "Calculated Pixel (" << label << ") != Cell Data (" << data << ").\n"; + throw 0; + } + } +} + +int +CuberilleTest03Parameters(const bool triangles, const bool project) +{ + + const auto image = CuberilleTest03CreateImage(); + const auto image_to_mesh = TImageToMesh::New(); + image_to_mesh->SetInput(image); + image_to_mesh->SavePixelAsCellDataOn(); + image_to_mesh->SetGenerateTriangleFaces(triangles); + image_to_mesh->SetProjectVerticesToIsoSurface(project); + ITK_TRY_EXPECT_NO_EXCEPTION(image_to_mesh->Update()); + ITK_TEST_EXPECT_EQUAL(image_to_mesh->GetOutput()->GetCells()->Size(), + image_to_mesh->GetOutput()->GetCellData()->Size()); + + return EXIT_SUCCESS; +} + +// - Create a test image, 3x3x5, with zeros along the edges and [1 2 3] +// in the interior. +// - Extract a Cuberille mesh calling SavePixelAsCellDataOn(). +// - For each cell in the mesh, calculate the centroid and normal. +// - Walk one-half pixel width from the centroid along the normal into the mesh. +// - Sample that point in the input segmentation image. +// - Assert that that cell data is equal to the adjacent pixel value. +int +CuberilleTest03(int itkNotUsed(argc), char * itkNotUsed(argv)[]) +{ + + { // Test Set/Get Methods + const auto image_to_mesh = TImageToMesh::New(); + bool save = true; + ITK_TEST_SET_GET_BOOLEAN(image_to_mesh, SavePixelAsCellData, save); + } + + if (EXIT_FAILURE == CuberilleTest03Parameters(false, false)) + { + throw 0; + } + if (EXIT_FAILURE == CuberilleTest03Parameters(false, true)) + { + throw 0; + } + if (EXIT_FAILURE == CuberilleTest03Parameters(true, false)) + { + throw 0; + } + if (EXIT_FAILURE == CuberilleTest03Parameters(true, true)) + { + throw 0; + } + + const auto image = CuberilleTest03CreateImage(); + CuberilleTest03Helper(image); + + return EXIT_SUCCESS; +} diff --git a/Modules/Filtering/Cuberille/test/CuberilleTest04.cxx b/Modules/Filtering/Cuberille/test/CuberilleTest04.cxx new file mode 100644 index 000000000000..9ec59091e729 --- /dev/null +++ b/Modules/Filtering/Cuberille/test/CuberilleTest04.cxx @@ -0,0 +1,94 @@ +/*========================================================================= + * + * 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 +#include +#include +#include +#include + +int +CuberilleTest04(int itkNotUsed(argc), char * itkNotUsed(argv)[]) +{ + + using TPixel = unsigned char; + using TCoordinate = double; + const unsigned int Dimension = 3; + using TImage = itk::Image; + using TMesh = itk::QuadEdgeMesh; + using TPad = itk::ConstantPadImageFilter; + using TExtract = itk::CuberilleImageToMeshFilter; + + for (size_t mask = 1; mask < std::pow(2, 8); ++mask) + { + + std::bitset<8> bitmask(mask); + + std::cout << mask << ' ' << bitmask << std::endl; + + const auto image = TImage::New(); + TImage::SizeType size; + size.Fill(2); + + TImage::IndexType origin; + origin.Fill(0); + + TImage::RegionType region(origin, size); + + image->SetRegions(region); + image->Allocate(); + image->FillBuffer(0); + + for (size_t index = 0; index < std::pow(2, 3); ++index) + { + std::bitset<3> bitindex(index); + image->SetPixel({ { bitindex[0], bitindex[1], bitindex[2] } }, bitmask[index]); + } + + typename TImage::SizeType padding; + padding.Fill(1); + + const auto pad = TPad::New(); + pad->SetInput(image); + pad->SetPadUpperBound(padding); + pad->SetPadLowerBound(padding); + pad->SetConstant(static_cast(0)); + + const auto extract = TExtract::New(); + extract->SetInput(pad->GetOutput()); + extract->GenerateTriangleFacesOn(); + extract->ProjectVerticesToIsoSurfaceOff(); + extract->SavePixelAsCellDataOn(); + ITK_TRY_EXPECT_NO_EXCEPTION(extract->Update()); + + const auto out = TMesh::New(); + out->Graft(extract->GetOutput()); + + ITK_TEST_EXPECT_TRUE(out->GetNumberOfFaces() == out->GetCellData()->Size()); + + for (auto it = out->GetEdgeCells()->Begin(); it != out->GetEdgeCells()->End(); ++it) + { + using TEdge = TMesh::EdgeCellType; + const auto qe = dynamic_cast(it.Value())->GetQEGeom(); + ITK_TEST_EXPECT_TRUE(qe->IsLeftSet()); + ITK_TEST_EXPECT_TRUE(qe->IsRightSet()); + } + } + + return EXIT_SUCCESS; +} diff --git a/Modules/Filtering/Cuberille/test/CuberilleTest_Issue66.cxx b/Modules/Filtering/Cuberille/test/CuberilleTest_Issue66.cxx new file mode 100644 index 000000000000..51ba04fa5d29 --- /dev/null +++ b/Modules/Filtering/Cuberille/test/CuberilleTest_Issue66.cxx @@ -0,0 +1,150 @@ +/*========================================================================= + * + * 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 +#include +#include +#include +#include + +const unsigned int Dimension = 3; +using TPixel = unsigned char; +using TImage = itk::Image; +using TMesh = itk::Mesh; +// using TMesh = itk::QuadEdgeMesh; +using TExtract = itk::CuberilleImageToMeshFilter; +using TMeshWriter = itk::MeshFileWriter; + +/* This is description from https://github.com/InsightSoftwareConsortium/ITKCuberille/issues/66 + +Imagine an image with only two foreground pixels, offset by {1, 1, 0}. When extracting the surface of this image as a +mesh, it is necessary to duplicate the coincident edge, creating two disconnected cubes, in order to ensure that the +result is manifold. A slice through the z plane would show: + +*---*---* +| A | | +*---*---* +| | B | +*---*---* +Now imagine that these foreground pixels are still flanked by background pixels in the z plain, but are "connected" into +the same connected component by foreground pixels in the two adjacent z plains: + + Z-1 Z Z+1 +*---*---* *---*---* *---*---* +| X | X | | A | | | X | X | +*---*---* *---*---* *---*---* +| | X | | | B | | | X | +*---*---* *---*---* *---*---* +Here, the coincident edge of pixels A and B must again be duplicated; but unlike the first example, the two distinct +edges actually connect the same two vertices. Quite reasonably, when asked to create a face which would require an edge +between two vertices, itk::QuadEdgeMesh checks: +https://github.com/InsightSoftwareConsortium/ITK/blob/v5.4.5/Modules/Core/QuadEdgeMesh/include/itkQuadEdgeMesh.hxx#L1212-L1221 +to see whether there is an existing edge, and refuses to create a new +edge if one already exists. + +Most of the time, this is what we want--but in the case described above, it acts as a bug. It's not immediately obvious +to me how to detect when this check is and is not necessary; as I work through this problem, I thought I would post it +as an issue in case anyone else had any insight--pun intended. ;-D + +*/ + +TImage::Pointer +case1() +{ + const auto image = TImage::New(); + TImage::RegionType region({ { 0, 0, 0 } }, { { 5, 4, 4 } }); + image->SetBufferedRegion(region); + image->AllocateInitialized(); // zero filled image + + image->SetPixel({ { 1, 1, 1 } }, 1); // A + image->SetPixel({ { 2, 2, 1 } }, 1); // B + + image->SetPixel({ { 1, 1, 0 } }, 1); // C + image->SetPixel({ { 1, 2, 0 } }, 1); // D + image->SetPixel({ { 2, 2, 0 } }, 1); // E + + image->SetPixel({ { 1, 1, 2 } }, 1); // F + image->SetPixel({ { 1, 2, 2 } }, 1); // G + image->SetPixel({ { 2, 2, 2 } }, 1); // H + + return image; +} + +TImage::Pointer +case2() +{ + const auto image = TImage::New(); + TImage::RegionType region({ { 0, 0, 0 } }, { { 5, 4, 4 } }); + image->SetBufferedRegion(region); + image->AllocateInitialized(); // zero filled image + + image->SetPixel({ { 1, 1, 1 } }, 1); // A + image->SetPixel({ { 2, 1, 1 } }, 1); // B + + image->SetPixel({ { 3, 1, 1 } }, 1); // C + image->SetPixel({ { 1, 2, 1 } }, 1); // D + image->SetPixel({ { 3, 2, 1 } }, 1); // E + + image->SetPixel({ { 1, 2, 2 } }, 1); // F + image->SetPixel({ { 2, 2, 2 } }, 1); // G + image->SetPixel({ { 3, 2, 2 } }, 1); // H + + return image; +} + + +int +CuberilleTest_Issue66(int, char *[]) +{ + int countFailed = 0; + for (const auto image : { case1(), case2() }) + { + + const auto extract = TExtract::New(); + extract->SetInput(image); + extract->SavePixelAsCellDataOn(); + extract->GenerateTriangleFacesOff(); + extract->ProjectVerticesToIsoSurfaceOff(); + + const auto m_writer = TMeshWriter::New(); + m_writer->SetInput(extract->GetOutput()); + m_writer->SetFileName("mesh.obj"); + m_writer->Update(); + + const auto n_cell = extract->GetOutput()->GetNumberOfCells(); + const auto n_data = extract->GetOutput()->GetCellData()->Size(); + + if (n_cell != n_data) + { + std::cout << "This fails if using itk::QuadEdgeMesh." << std::endl; + std::cout << "n_cell: " << n_cell << std::endl; + std::cout << "n_data: " << n_data << std::endl; + ++countFailed; + } + else + { + std::cout << "This succeeds if using itk::Mesh." << std::endl; + } + } + + if (countFailed == 0) + { + return EXIT_SUCCESS; + } + return EXIT_FAILURE; +} diff --git a/Modules/Filtering/Cuberille/test/Input/blob0.mha b/Modules/Filtering/Cuberille/test/Input/blob0.mha new file mode 100644 index 000000000000..389a3d84475e Binary files /dev/null and b/Modules/Filtering/Cuberille/test/Input/blob0.mha differ diff --git a/Modules/Filtering/Cuberille/test/Input/blob1.mha b/Modules/Filtering/Cuberille/test/Input/blob1.mha new file mode 100644 index 000000000000..0d40cfc358c1 Binary files /dev/null and b/Modules/Filtering/Cuberille/test/Input/blob1.mha differ diff --git a/Modules/Filtering/Cuberille/test/Input/blob2.mha b/Modules/Filtering/Cuberille/test/Input/blob2.mha new file mode 100644 index 000000000000..1549ca98ef98 Binary files /dev/null and b/Modules/Filtering/Cuberille/test/Input/blob2.mha differ diff --git a/Modules/Filtering/Cuberille/test/Input/blob3.mha b/Modules/Filtering/Cuberille/test/Input/blob3.mha new file mode 100644 index 000000000000..df856c829cc9 Binary files /dev/null and b/Modules/Filtering/Cuberille/test/Input/blob3.mha differ diff --git a/Modules/Filtering/Cuberille/test/Input/blob4.mha b/Modules/Filtering/Cuberille/test/Input/blob4.mha new file mode 100644 index 000000000000..adf54d971747 Binary files /dev/null and b/Modules/Filtering/Cuberille/test/Input/blob4.mha differ diff --git a/Modules/Filtering/Cuberille/test/Input/fuel.mha b/Modules/Filtering/Cuberille/test/Input/fuel.mha new file mode 100644 index 000000000000..b129c23631ae Binary files /dev/null and b/Modules/Filtering/Cuberille/test/Input/fuel.mha differ diff --git a/Modules/Filtering/Cuberille/test/Input/hydrogenAtom.mha b/Modules/Filtering/Cuberille/test/Input/hydrogenAtom.mha new file mode 100644 index 000000000000..ebf06c19d1d6 Binary files /dev/null and b/Modules/Filtering/Cuberille/test/Input/hydrogenAtom.mha differ diff --git a/Modules/Filtering/Cuberille/test/Input/marschnerlobb.mha b/Modules/Filtering/Cuberille/test/Input/marschnerlobb.mha new file mode 100644 index 000000000000..ddae2378fc33 Binary files /dev/null and b/Modules/Filtering/Cuberille/test/Input/marschnerlobb.mha differ diff --git a/Modules/Filtering/Cuberille/test/Input/neghip.mha b/Modules/Filtering/Cuberille/test/Input/neghip.mha new file mode 100644 index 000000000000..f350b557bd02 Binary files /dev/null and b/Modules/Filtering/Cuberille/test/Input/neghip.mha differ diff --git a/Modules/Filtering/Cuberille/test/Input/nucleon.mha b/Modules/Filtering/Cuberille/test/Input/nucleon.mha new file mode 100644 index 000000000000..52cce95dfa74 Binary files /dev/null and b/Modules/Filtering/Cuberille/test/Input/nucleon.mha differ diff --git a/Modules/Filtering/Cuberille/test/Input/silicium.mha b/Modules/Filtering/Cuberille/test/Input/silicium.mha new file mode 100644 index 000000000000..d33093e0d522 Binary files /dev/null and b/Modules/Filtering/Cuberille/test/Input/silicium.mha differ diff --git a/Modules/Filtering/Cuberille/wrapping/CMakeLists.txt b/Modules/Filtering/Cuberille/wrapping/CMakeLists.txt new file mode 100644 index 000000000000..23dd508f7a58 --- /dev/null +++ b/Modules/Filtering/Cuberille/wrapping/CMakeLists.txt @@ -0,0 +1,3 @@ +itk_wrap_module(Cuberille) +itk_auto_load_submodules() +itk_end_wrap_module() diff --git a/Modules/Filtering/Cuberille/wrapping/itkCuberilleImageToMeshFilter.wrap b/Modules/Filtering/Cuberille/wrapping/itkCuberilleImageToMeshFilter.wrap new file mode 100644 index 000000000000..cedbcbc1ea19 --- /dev/null +++ b/Modules/Filtering/Cuberille/wrapping/itkCuberilleImageToMeshFilter.wrap @@ -0,0 +1,10 @@ +itk_wrap_class("itk::CuberilleImageToMeshFilter" POINTER_WITH_2_SUPERCLASSES) +unique(types "${WRAP_ITK_SCALAR}") +foreach(d ${ITK_WRAP_IMAGE_DIMS}) + foreach(t ${types}) + itk_wrap_template("${ITKM_I${t}${d}}M${ITKM_${t}}${d}" + "${ITKT_I${t}${d}}, itk::Mesh< ${ITKT_${t}},${d} >, itk::InterpolateImageFunction< ${ITKT_I${t}${d}}, ${ITKT_D} >" + ) + endforeach() +endforeach() +itk_end_wrap_class() diff --git a/Modules/Filtering/Cuberille/wrapping/test/CMakeLists.txt b/Modules/Filtering/Cuberille/wrapping/test/CMakeLists.txt new file mode 100644 index 000000000000..e61e8418fe79 --- /dev/null +++ b/Modules/Filtering/Cuberille/wrapping/test/CMakeLists.txt @@ -0,0 +1,49 @@ +itk_python_expression_add_test( + NAME itkCuberilleImageToMeshFilterPythonTest + EXPRESSION "filter = itk.CuberilleImageToMeshFilter.New()" +) + +set(test_input_dir ${CMAKE_CURRENT_SOURCE_DIR}/../../test/Input) +# cmake_path(CONVERT ${test_input_dir} TO_CMAKE_PATH_LIST test_input_dir NORMALIZE) # CMake 3.20+ + +list(FIND ITK_WRAP_IMAGE_DIMS 3 wrap_3_index) +if(wrap_3_index GREATER -1) + itk_python_add_test( + NAME Cuberille_blob0_python + TEST_DRIVER_ARGS + COMMAND + cuberille_test.py + DATA{${test_input_dir}/blob0.mha} + ${ITK_TEST_OUTPUT_DIR}/blob0.vtk + ) + + itk_python_add_test( + NAME Cuberille_blob0_100_python + TEST_DRIVER_ARGS + COMMAND + cuberille_test.py + DATA{${test_input_dir}/blob0.mha} + ${ITK_TEST_OUTPUT_DIR}/blob0_100.vtk + 100 # iso-surface value + ) + + itk_python_add_test( + NAME Cuberille_hydrogenAtom_35_python + TEST_DRIVER_ARGS + COMMAND + cuberille_test.py + DATA{${test_input_dir}/hydrogenAtom.mha} + ${ITK_TEST_OUTPUT_DIR}/hydrogenAtom35.vtk + 35 # iso-surface value + ) + + itk_python_add_test( + NAME Cuberille_hydrogenAtom_50_python + TEST_DRIVER_ARGS + COMMAND + cuberille_test.py + DATA{${test_input_dir}/hydrogenAtom.mha} + ${ITK_TEST_OUTPUT_DIR}/hydrogenAtom50.vtk + 50 # iso-surface value + ) +endif() diff --git a/Modules/Filtering/Cuberille/wrapping/test/cuberille_test.py b/Modules/Filtering/Cuberille/wrapping/test/cuberille_test.py new file mode 100644 index 000000000000..ee90a25daa73 --- /dev/null +++ b/Modules/Filtering/Cuberille/wrapping/test/cuberille_test.py @@ -0,0 +1,30 @@ +import sys + +import itk + +image = itk.imread(sys.argv[1]) + +if len(sys.argv) > 3: + isosurface = sys.argv[3] +else: + calculator = itk.MinimumMaximumImageCalculator[type(image)].New() + calculator.SetImage(image) + calculator.Compute() + isosurface = (calculator.GetMaximum() + calculator.GetMinimum()) / 2 + +# convert number type +(input_image_template, (input_pixel_type, input_image_dimension)) = itk.template(image) +if input_pixel_type == itk.F or input_pixel_type == itk.D: + isosurface = float(isosurface) +else: + isosurface = int(isosurface) + +mesh = itk.cuberille_image_to_mesh_filter( + image, + generate_triangle_faces=True, + iso_surface_value=isosurface, + project_vertices_to_iso_surface=True, + project_vertex_surface_distance_threshold=0.05, +) + +itk.meshwrite(mesh, sys.argv[2]) diff --git a/Modules/Remote/Cuberille.remote.cmake b/Modules/Remote/Cuberille.remote.cmake deleted file mode 100644 index 78046b2839e7..000000000000 --- a/Modules/Remote/Cuberille.remote.cmake +++ /dev/null @@ -1,72 +0,0 @@ -#-- # Grading Level Criteria Report -#-- EVALUATION DATE: 2020-03-24 -#-- EVALUATORS: [Dženan Zukić, Davis Vigneault] -#-- -#-- ## Compliance level 5 star (AKA ITK main modules, or remote modules that could become core modules) -#-- - [ ] Widespread community dependance -#-- - [X] Above 90% code coverage -#-- - [ ] CI dashboards and testing monitored rigorously -#-- - [X] Key API features are exposed in wrapping interface -#-- - [ ] All requirements of Levels 4,3,2,1 -#-- -#-- ## Compliance Level 4 star (Very high-quality code, perhaps small community dependance) -#-- - [X] Meets all ITK code style standards -#-- - [X] No external requirements beyond those needed by ITK proper -#-- - [ ] Builds and passes tests on all supported platforms within 1 month of each core tagged release -#-- - [ ] Windows Shared Library Build with Visual Studio -#-- - [ ] Mac with clang compiller -#-- - [ ] Linux with gcc compiler -#-- - [X] Active developer community dedicated to maintaining code-base -#-- - [X] 75% code coverage demonstrated for testing suite -#-- - [X] Continuous integration testing performed -#-- - [X] All requirements of Levels 3,2,1 -#-- -#-- ## Compliance Level 3 star (Quality beta code) -#-- - [X] API | executable interface is considered mostly stable and feature complete -#-- - [X] 10% C0-code coverage demonstrated for testing suite -#-- - [X] Some tests exist and pass on at least some platform -#-- - [X] All requirements of Levels 2,1 -#-- -#-- ## Compliance Level 2 star (Alpha code feature API development or niche community/execution environment dependance ) -#-- - [X] Compiles for at least 1 niche set of execution envirionments, and perhaps others -#-- (may depend on specific external tools like a java environment, or specific external libraries to work ) -#-- - [X] All requirements of Levels 1 -#-- -#-- ## Compliance Level 1 star (Pre-alpha features under development and code of unknown quality) -#-- - [X] Code complies on at least 1 platform -#-- -#-- ## Compliance Level 0 star ( Code/Feature of known poor-quality or deprecated status ) -#-- - [ ] Code reviewed and explicitly identified as not recommended for use -#-- -#-- ### Please document here any justification for the criteria above -# Code style enforced by clang-format on 2020-02-19, and clang-tidy modernizations completed -# 2020-05-25 Continuous integration and clang format linting now performed by Github Actions -# 2020-05-25 Line coverage is 99.8%, function coverage is 98.0% - -# Contact: Matt McCormick -itk_fetch_module( - Cuberille - "This module implements cuberille implicit surface -polygonization for ITK. This method operates by diving the surface into a -number of small cubes called cuberilles. Each cuberille is centered at a -pixel lying on the iso-surface and then quadrilaterals are generated for each -face. The original approach is improved by projecting the vertices of each -cuberille onto the implicit surface, smoothing the typical block-like -resultant mesh. - -A more detailed description can be found in the Insight Journal article: - - Mueller, D. \"Cuberille Implicit Surface Polygonization for ITK\" - https://doi.org/10.54294/df9mgw - July-December, 2010. - -And the related: - - Mueller, D. \"Fast Marching Minimal Path Extraction in ITK\" - https://doi.org/10.54294/z5zwhh - January-June, 2008. -" - MODULE_COMPLIANCE_LEVEL 3 - GIT_REPOSITORY https://github.com/InsightSoftwareConsortium/ITKCuberille.git - GIT_TAG a152f89e5d45f7262523446555e1f4da86071a71 - ) diff --git a/pyproject.toml b/pyproject.toml index d75121de752e..52c56027848e 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -53,6 +53,7 @@ cmd = '''cmake -DCMAKE_C_COMPILER_LAUNCHER:STRING=ccache -DCMAKE_CXX_COMPILER_LAUNCHER:STRING=ccache -DModule_AnisotropicDiffusionLBR:BOOL=ON + -DModule_Cuberille:BOOL=ON -DModule_Montage:BOOL=ON -DModule_FastBilateral:BOOL=ON -DModule_GenericLabelInterpolator:BOOL=ON