--- /dev/null
+#ifndef __itkBinaryThinningImageFilter3D_h\r
+#define __itkBinaryThinningImageFilter3D_h\r
+\r
+#include <itkNeighborhoodIterator.h>\r
+#include <itkImageToImageFilter.h>\r
+#include <itkImageRegionIteratorWithIndex.h>\r
+#include <itkConstantBoundaryCondition.h>\r
+\r
+namespace itk\r
+{\r
+/** \class BinaryThinningImageFilter3D\r
+*\r
+* \brief This filter computes one-pixel-wide skeleton of a 3D input image.\r
+*\r
+* This class is parametrized over the type of the input image\r
+* and the type of the output image.\r
+* \r
+* The input is assumed to be a binary image. All non-zero valued voxels\r
+* are set to 1 internally to simplify the computation. The filter will\r
+* produce a skeleton of the object. The output background values are 0,\r
+* and the foreground values are 1.\r
+* \r
+* A 26-neighbourhood configuration is used for the foreground and a\r
+* 6-neighbourhood configuration for the background. Thinning is performed\r
+* symmetrically in order to guarantee that the skeleton lies medial within\r
+* the object.\r
+*\r
+* This filter is a parallel thinning algorithm and is an implementation\r
+* of the algorithm described in:\r
+* \r
+* T.C. Lee, R.L. Kashyap, and C.N. Chu.\r
+* Building skeleton models via 3-D medial surface/axis thinning algorithms.\r
+* Computer Vision, Graphics, and Image Processing, 56(6):462--478, 1994.\r
+* \r
+* To do: Make use of multi-threading.\r
+*\r
+* \author Hanno Homann, Oxford University, Wolfson Medical Vision Lab, UK.\r
+* \r
+* \sa MorphologyImageFilter\r
+* \ingroup ImageEnhancement MathematicalMorphologyImageFilters\r
+*/\r
+\r
+template <class TInputImage,class TOutputImage>\r
+class BinaryThinningImageFilter3D :\r
+ public ImageToImageFilter<TInputImage,TOutputImage>\r
+{\r
+public:\r
+ /** Standard class typedefs. */\r
+ typedef BinaryThinningImageFilter3D Self;\r
+ typedef ImageToImageFilter<TInputImage,TOutputImage> Superclass;\r
+ typedef SmartPointer<Self> Pointer;\r
+ typedef SmartPointer<const Self> ConstPointer;\r
+\r
+ /** Method for creation through the object factory */\r
+ itkNewMacro(Self);\r
+\r
+ /** Run-time type information (and related methods). */\r
+ itkTypeMacro( BinaryThinningImageFilter3D, ImageToImageFilter );\r
+\r
+ /** Type for input image. */\r
+ typedef TInputImage InputImageType;\r
+\r
+ /** Type for output image: Skelenton of the object. */\r
+ typedef TOutputImage OutputImageType;\r
+\r
+ /** Type for the region of the input image. */\r
+ typedef typename InputImageType::RegionType RegionType;\r
+\r
+ /** Type for the index of the input image. */\r
+ typedef typename RegionType::IndexType IndexType;\r
+\r
+ /** Type for the pixel type of the input image. */\r
+ typedef typename InputImageType::PixelType InputImagePixelType ;\r
+\r
+ /** Type for the pixel type of the input image. */\r
+ typedef typename OutputImageType::PixelType OutputImagePixelType ;\r
+\r
+ /** Type for the size of the input image. */\r
+ typedef typename RegionType::SizeType SizeType;\r
+\r
+ /** Pointer Type for input image. */\r
+ typedef typename InputImageType::ConstPointer InputImagePointer;\r
+\r
+ /** Pointer Type for the output image. */\r
+ typedef typename OutputImageType::Pointer OutputImagePointer;\r
+ \r
+ /** Boundary condition type for the neighborhood iterator */\r
+ typedef ConstantBoundaryCondition< TInputImage > ConstBoundaryConditionType;\r
+ \r
+ /** Neighborhood iterator type */\r
+ typedef NeighborhoodIterator<TInputImage, ConstBoundaryConditionType> NeighborhoodIteratorType;\r
+ \r
+ /** Neighborhood type */\r
+ typedef typename NeighborhoodIteratorType::NeighborhoodType NeighborhoodType;\r
+\r
+ /** Get Skelenton by thinning image. */\r
+ OutputImageType * GetThinning(void);\r
+\r
+ /** ImageDimension enumeration */\r
+ itkStaticConstMacro(InputImageDimension, unsigned int,\r
+ TInputImage::ImageDimension );\r
+ itkStaticConstMacro(OutputImageDimension, unsigned int,\r
+ TOutputImage::ImageDimension );\r
+\r
+#ifdef ITK_USE_CONCEPT_CHECKING\r
+ /** Begin concept checking */\r
+ itkConceptMacro(SameDimensionCheck,\r
+ (Concept::SameDimension<InputImageDimension, 3>));\r
+ itkConceptMacro(SameTypeCheck,\r
+ (Concept::SameType<InputImagePixelType, OutputImagePixelType>));\r
+ itkConceptMacro(InputAdditiveOperatorsCheck,\r
+ (Concept::AdditiveOperators<InputImagePixelType>));\r
+ itkConceptMacro(InputConvertibleToIntCheck,\r
+ (Concept::Convertible<InputImagePixelType, int>));\r
+ itkConceptMacro(IntConvertibleToInputCheck,\r
+ (Concept::Convertible<int, InputImagePixelType>));\r
+ itkConceptMacro(InputIntComparableCheck,\r
+ (Concept::Comparable<InputImagePixelType, int>));\r
+ /** End concept checking */\r
+#endif\r
+\r
+protected:\r
+ BinaryThinningImageFilter3D();\r
+ virtual ~BinaryThinningImageFilter3D() {};\r
+ void PrintSelf(std::ostream& os, Indent indent) const;\r
+\r
+ /** Compute thinning Image. */\r
+ void GenerateData();\r
+\r
+ /** Prepare data. */\r
+ void PrepareData();\r
+\r
+ /** Compute thinning Image. */\r
+ void ComputeThinImage();\r
+ \r
+ /** isEulerInvariant [Lee94] */\r
+ bool isEulerInvariant(NeighborhoodType neighbors, int *LUT);\r
+ void fillEulerLUT(int *LUT); \r
+ /** isSimplePoint [Lee94] */\r
+ bool isSimplePoint(NeighborhoodType neighbors);\r
+ /** Octree_labeling [Lee94] */\r
+ void Octree_labeling(int octant, int label, int *cube);\r
+\r
+\r
+private: \r
+ BinaryThinningImageFilter3D(const Self&); //purposely not implemented\r
+ void operator=(const Self&); //purposely not implemented\r
+\r
+}; // end of BinaryThinningImageFilter3D class\r
+\r
+} //end namespace itk\r
+\r
+#ifndef ITK_MANUAL_INSTANTIATION\r
+#include "itkBinaryThinningImageFilter3D.txx"\r
+#endif\r
+\r
+#endif\r
--- /dev/null
+#ifndef _itkBinaryThinningImageFilter3D_txx\r
+#define _itkBinaryThinningImageFilter3D_txx\r
+\r
+#include <iostream>\r
+\r
+#include "itkBinaryThinningImageFilter3D.h"\r
+#include "itkImageRegionConstIterator.h"\r
+#include "itkImageRegionIterator.h"\r
+#include "itkNeighborhoodIterator.h"\r
+#include <vector>\r
+\r
+namespace itk\r
+{\r
+\r
+/**\r
+ * Constructor\r
+ */\r
+template <class TInputImage,class TOutputImage>\r
+BinaryThinningImageFilter3D<TInputImage,TOutputImage>\r
+::BinaryThinningImageFilter3D()\r
+{\r
+\r
+ this->SetNumberOfRequiredOutputs( 1 );\r
+\r
+ OutputImagePointer thinImage = OutputImageType::New();\r
+ this->SetNthOutput( 0, thinImage.GetPointer() );\r
+\r
+}\r
+\r
+/**\r
+ * Return the thinning Image pointer\r
+ */\r
+template <class TInputImage,class TOutputImage>\r
+typename BinaryThinningImageFilter3D<\r
+ TInputImage,TOutputImage>::OutputImageType * \r
+BinaryThinningImageFilter3D<TInputImage,TOutputImage>\r
+::GetThinning(void)\r
+{\r
+ return dynamic_cast< OutputImageType * >(\r
+ this->ProcessObject::GetOutput(0) );\r
+}\r
+\r
+\r
+/**\r
+ * Prepare data for computation\r
+ * Copy the input image to the output image, changing from the input\r
+ * type to the output type.\r
+ */\r
+template <class TInputImage,class TOutputImage>\r
+void \r
+BinaryThinningImageFilter3D<TInputImage,TOutputImage>\r
+::PrepareData(void) \r
+{\r
+ \r
+ itkDebugMacro(<< "PrepareData Start");\r
+ OutputImagePointer thinImage = GetThinning();\r
+\r
+ InputImagePointer inputImage = \r
+ dynamic_cast<const TInputImage *>( ProcessObject::GetInput(0) );\r
+\r
+ thinImage->SetBufferedRegion( thinImage->GetRequestedRegion() );\r
+ thinImage->Allocate();\r
+\r
+ typename OutputImageType::RegionType region = thinImage->GetRequestedRegion();\r
+\r
+\r
+ ImageRegionConstIterator< TInputImage > it( inputImage, region );\r
+ ImageRegionIterator< TOutputImage > ot( thinImage, region );\r
+\r
+ it.GoToBegin();\r
+ ot.GoToBegin();\r
+\r
+ itkDebugMacro(<< "PrepareData: Copy input to output");\r
+ \r
+ // Copy the input to the output, changing all foreground pixels to\r
+ // have value 1 in the process.\r
+ while( !ot.IsAtEnd() )\r
+ {\r
+ if ( it.Get() )\r
+ {\r
+ ot.Set( NumericTraits<OutputImagePixelType>::One );\r
+ }\r
+ else\r
+ {\r
+ ot.Set( NumericTraits<OutputImagePixelType>::Zero );\r
+ }\r
+ ++it;\r
+ ++ot;\r
+ }\r
+ itkDebugMacro(<< "PrepareData End"); \r
+}\r
+\r
+/**\r
+ * Post processing for computing thinning\r
+ */\r
+template <class TInputImage,class TOutputImage>\r
+void \r
+BinaryThinningImageFilter3D<TInputImage,TOutputImage>\r
+::ComputeThinImage() \r
+{\r
+ itkDebugMacro( << "ComputeThinImage Start");\r
+ OutputImagePointer thinImage = GetThinning();\r
+\r
+ typename OutputImageType::RegionType region = thinImage->GetRequestedRegion();\r
+ \r
+ ConstBoundaryConditionType boundaryCondition;\r
+ boundaryCondition.SetConstant( 0 );\r
+\r
+ typename NeighborhoodIteratorType::RadiusType radius;\r
+ radius.Fill(1);\r
+ NeighborhoodIteratorType ot( radius, thinImage, region );\r
+ ot.SetBoundaryCondition( boundaryCondition );\r
+\r
+ std::vector < IndexType > simpleBorderPoints;\r
+ typename std::vector < IndexType >::iterator simpleBorderPointsIt;\r
+\r
+ // Define offsets\r
+ typedef typename NeighborhoodIteratorType::OffsetType OffsetType;\r
+ OffsetType N = {{ 0,-1, 0}}; // north\r
+ OffsetType S = {{ 0, 1, 0}}; // south\r
+ OffsetType E = {{ 1, 0, 0}}; // east\r
+ OffsetType W = {{-1, 0, 0}}; // west\r
+ OffsetType U = {{ 0, 0, 1}}; // up\r
+ OffsetType B = {{ 0, 0,-1}}; // bottom\r
+\r
+ // prepare Euler LUT [Lee94]\r
+ int eulerLUT[256]; \r
+ fillEulerLUT( eulerLUT );\r
+ // Loop through the image several times until there is no change.\r
+ int unchangedBorders = 0;\r
+ while( unchangedBorders < 6 ) // loop until no change for all the six border types\r
+ {\r
+ unchangedBorders = 0;\r
+ for( int currentBorder = 1; currentBorder <= 6; currentBorder++)\r
+ {\r
+ // Loop through the image.\r
+ for ( ot.GoToBegin(); !ot.IsAtEnd(); ++ot )\r
+ { \r
+ // check if point is foreground\r
+ if ( ot.GetCenterPixel() != 1 )\r
+ {\r
+ continue; // current point is already background \r
+ }\r
+ // check 6-neighbors if point is a border point of type currentBorder\r
+ bool isBorderPoint = false;\r
+ if( currentBorder == 1 && ot.GetPixel(N)<=0 )\r
+ isBorderPoint = true;\r
+ if( currentBorder == 2 && ot.GetPixel(S)<=0 )\r
+ isBorderPoint = true;\r
+ if( currentBorder == 3 && ot.GetPixel(E)<=0 )\r
+ isBorderPoint = true;\r
+ if( currentBorder == 4 && ot.GetPixel(W)<=0 )\r
+ isBorderPoint = true;\r
+ if( currentBorder == 5 && ot.GetPixel(U)<=0 )\r
+ isBorderPoint = true;\r
+ if( currentBorder == 6 && ot.GetPixel(B)<=0 )\r
+ isBorderPoint = true;\r
+ if( !isBorderPoint )\r
+ {\r
+ continue; // current point is not deletable\r
+ } \r
+ // check if point is the end of an arc\r
+ int numberOfNeighbors = -1; // -1 and not 0 because the center pixel will be counted as well \r
+ for( int i = 0; i < 27; i++ ) // i = 0..26\r
+ if( ot.GetPixel(i)==1 )\r
+ numberOfNeighbors++;\r
+\r
+ if( numberOfNeighbors == 1 )\r
+ {\r
+ continue; // current point is not deletable\r
+ }\r
+\r
+ // check if point is Euler invariant\r
+ if( !isEulerInvariant( ot.GetNeighborhood(), eulerLUT ) )\r
+ {\r
+ continue; // current point is not deletable\r
+ }\r
+\r
+ // check if point is simple (deletion does not change connectivity in the 3x3x3 neighborhood)\r
+ if( !isSimplePoint( ot.GetNeighborhood() ) )\r
+ {\r
+ continue; // current point is not deletable\r
+ }\r
+\r
+ // add all simple border points to a list for sequential re-checking\r
+ simpleBorderPoints.push_back( ot.GetIndex() );\r
+ } // end image iteration loop\r
+\r
+ // sequential re-checking to preserve connectivity when\r
+ // deleting in a parallel way\r
+ bool noChange = true;\r
+ for( simpleBorderPointsIt=simpleBorderPoints.begin(); simpleBorderPointsIt!=simpleBorderPoints.end(); simpleBorderPointsIt++)\r
+ {\r
+ // 1. Set simple border point to 0\r
+ thinImage->SetPixel( *simpleBorderPointsIt, NumericTraits<OutputImagePixelType>::Zero);\r
+ // 2. Check if neighborhood is still connected\r
+ ot.SetLocation( *simpleBorderPointsIt );\r
+ if( !isSimplePoint( ot.GetNeighborhood() ) )\r
+ {\r
+ // we cannot delete current point, so reset\r
+ thinImage->SetPixel( *simpleBorderPointsIt, NumericTraits<OutputImagePixelType>::One );\r
+ }\r
+ else\r
+ {\r
+ noChange = false;\r
+ }\r
+ }\r
+ if( noChange )\r
+ unchangedBorders++;\r
+\r
+ simpleBorderPoints.clear();\r
+ } // end currentBorder for loop\r
+ } // end unchangedBorders while loop\r
+\r
+ itkDebugMacro( << "ComputeThinImage End");\r
+}\r
+\r
+/**\r
+ * Generate ThinImage\r
+ */\r
+template <class TInputImage,class TOutputImage>\r
+void \r
+BinaryThinningImageFilter3D<TInputImage,TOutputImage>\r
+::GenerateData() \r
+{\r
+\r
+ this->PrepareData();\r
+\r
+ itkDebugMacro(<< "GenerateData: Computing Thinning Image");\r
+ this->ComputeThinImage();\r
+} // end GenerateData()\r
+\r
+/** \r
+ * Fill the Euler look-up table (LUT) for later check of the Euler invariance. (see [Lee94])\r
+ */\r
+template <class TInputImage,class TOutputImage>\r
+void \r
+BinaryThinningImageFilter3D<TInputImage,TOutputImage>\r
+::fillEulerLUT(int *LUT)\r
+{\r
+ LUT[1] = 1;\r
+ LUT[3] = -1;\r
+ LUT[5] = -1;\r
+ LUT[7] = 1;\r
+ LUT[9] = -3;\r
+ LUT[11] = -1;\r
+ LUT[13] = -1;\r
+ LUT[15] = 1;\r
+ LUT[17] = -1;\r
+ LUT[19] = 1;\r
+ LUT[21] = 1;\r
+ LUT[23] = -1;\r
+ LUT[25] = 3;\r
+ LUT[27] = 1;\r
+ LUT[29] = 1;\r
+ LUT[31] = -1;\r
+ LUT[33] = -3;\r
+ LUT[35] = -1;\r
+ LUT[37] = 3;\r
+ LUT[39] = 1;\r
+ LUT[41] = 1;\r
+ LUT[43] = -1;\r
+ LUT[45] = 3;\r
+ LUT[47] = 1;\r
+ LUT[49] = -1;\r
+ LUT[51] = 1;\r
+\r
+ LUT[53] = 1;\r
+ LUT[55] = -1;\r
+ LUT[57] = 3;\r
+ LUT[59] = 1;\r
+ LUT[61] = 1;\r
+ LUT[63] = -1;\r
+ LUT[65] = -3;\r
+ LUT[67] = 3;\r
+ LUT[69] = -1;\r
+ LUT[71] = 1;\r
+ LUT[73] = 1;\r
+ LUT[75] = 3;\r
+ LUT[77] = -1;\r
+ LUT[79] = 1;\r
+ LUT[81] = -1;\r
+ LUT[83] = 1;\r
+ LUT[85] = 1;\r
+ LUT[87] = -1;\r
+ LUT[89] = 3;\r
+ LUT[91] = 1;\r
+ LUT[93] = 1;\r
+ LUT[95] = -1;\r
+ LUT[97] = 1;\r
+ LUT[99] = 3;\r
+ LUT[101] = 3;\r
+ LUT[103] = 1;\r
+\r
+ LUT[105] = 5;\r
+ LUT[107] = 3;\r
+ LUT[109] = 3;\r
+ LUT[111] = 1;\r
+ LUT[113] = -1;\r
+ LUT[115] = 1;\r
+ LUT[117] = 1;\r
+ LUT[119] = -1;\r
+ LUT[121] = 3;\r
+ LUT[123] = 1;\r
+ LUT[125] = 1;\r
+ LUT[127] = -1;\r
+ LUT[129] = -7;\r
+ LUT[131] = -1;\r
+ LUT[133] = -1;\r
+ LUT[135] = 1;\r
+ LUT[137] = -3;\r
+ LUT[139] = -1;\r
+ LUT[141] = -1;\r
+ LUT[143] = 1;\r
+ LUT[145] = -1;\r
+ LUT[147] = 1;\r
+ LUT[149] = 1;\r
+ LUT[151] = -1;\r
+ LUT[153] = 3;\r
+ LUT[155] = 1;\r
+\r
+ LUT[157] = 1;\r
+ LUT[159] = -1;\r
+ LUT[161] = -3;\r
+ LUT[163] = -1;\r
+ LUT[165] = 3;\r
+ LUT[167] = 1;\r
+ LUT[169] = 1;\r
+ LUT[171] = -1;\r
+ LUT[173] = 3;\r
+ LUT[175] = 1;\r
+ LUT[177] = -1;\r
+ LUT[179] = 1;\r
+ LUT[181] = 1;\r
+ LUT[183] = -1;\r
+ LUT[185] = 3;\r
+ LUT[187] = 1;\r
+ LUT[189] = 1;\r
+ LUT[191] = -1;\r
+ LUT[193] = -3;\r
+ LUT[195] = 3;\r
+ LUT[197] = -1;\r
+ LUT[199] = 1;\r
+ LUT[201] = 1;\r
+ LUT[203] = 3;\r
+ LUT[205] = -1;\r
+ LUT[207] = 1;\r
+\r
+ LUT[209] = -1;\r
+ LUT[211] = 1;\r
+ LUT[213] = 1;\r
+ LUT[215] = -1;\r
+ LUT[217] = 3;\r
+ LUT[219] = 1;\r
+ LUT[221] = 1;\r
+ LUT[223] = -1;\r
+ LUT[225] = 1;\r
+ LUT[227] = 3;\r
+ LUT[229] = 3;\r
+ LUT[231] = 1;\r
+ LUT[233] = 5;\r
+ LUT[235] = 3;\r
+ LUT[237] = 3;\r
+ LUT[239] = 1;\r
+ LUT[241] = -1;\r
+ LUT[243] = 1;\r
+ LUT[245] = 1;\r
+ LUT[247] = -1;\r
+ LUT[249] = 3;\r
+ LUT[251] = 1;\r
+ LUT[253] = 1;\r
+ LUT[255] = -1;\r
+}\r
+\r
+/** \r
+ * Check for Euler invariance. (see [Lee94])\r
+ */\r
+template <class TInputImage,class TOutputImage>\r
+bool \r
+BinaryThinningImageFilter3D<TInputImage,TOutputImage>\r
+::isEulerInvariant(NeighborhoodType neighbors, int *LUT)\r
+{\r
+ // calculate Euler characteristic for each octant and sum up\r
+ int EulerChar = 0;\r
+ unsigned char n;\r
+ // Octant SWU\r
+ n = 1;\r
+ if( neighbors[24]==1 )\r
+ n |= 128;\r
+ if( neighbors[25]==1 )\r
+ n |= 64;\r
+ if( neighbors[15]==1 )\r
+ n |= 32;\r
+ if( neighbors[16]==1 )\r
+ n |= 16;\r
+ if( neighbors[21]==1 )\r
+ n |= 8;\r
+ if( neighbors[22]==1 )\r
+ n |= 4;\r
+ if( neighbors[12]==1 )\r
+ n |= 2;\r
+ EulerChar += LUT[n];\r
+ // Octant SEU\r
+ n = 1;\r
+ if( neighbors[26]==1 )\r
+ n |= 128;\r
+ if( neighbors[23]==1 )\r
+ n |= 64;\r
+ if( neighbors[17]==1 )\r
+ n |= 32;\r
+ if( neighbors[14]==1 )\r
+ n |= 16;\r
+ if( neighbors[25]==1 )\r
+ n |= 8;\r
+ if( neighbors[22]==1 )\r
+ n |= 4;\r
+ if( neighbors[16]==1 )\r
+ n |= 2;\r
+ EulerChar += LUT[n];\r
+ // Octant NWU\r
+ n = 1;\r
+ if( neighbors[18]==1 )\r
+ n |= 128;\r
+ if( neighbors[21]==1 )\r
+ n |= 64;\r
+ if( neighbors[9]==1 )\r
+ n |= 32;\r
+ if( neighbors[12]==1 )\r
+ n |= 16;\r
+ if( neighbors[19]==1 )\r
+ n |= 8;\r
+ if( neighbors[22]==1 )\r
+ n |= 4;\r
+ if( neighbors[10]==1 )\r
+ n |= 2;\r
+ EulerChar += LUT[n];\r
+ // Octant NEU\r
+ n = 1;\r
+ if( neighbors[20]==1 )\r
+ n |= 128;\r
+ if( neighbors[23]==1 )\r
+ n |= 64;\r
+ if( neighbors[19]==1 )\r
+ n |= 32;\r
+ if( neighbors[22]==1 )\r
+ n |= 16;\r
+ if( neighbors[11]==1 )\r
+ n |= 8;\r
+ if( neighbors[14]==1 )\r
+ n |= 4;\r
+ if( neighbors[10]==1 )\r
+ n |= 2;\r
+ EulerChar += LUT[n];\r
+ // Octant SWB\r
+ n = 1;\r
+ if( neighbors[6]==1 )\r
+ n |= 128;\r
+ if( neighbors[15]==1 )\r
+ n |= 64;\r
+ if( neighbors[7]==1 )\r
+ n |= 32;\r
+ if( neighbors[16]==1 )\r
+ n |= 16;\r
+ if( neighbors[3]==1 )\r
+ n |= 8;\r
+ if( neighbors[12]==1 )\r
+ n |= 4;\r
+ if( neighbors[4]==1 )\r
+ n |= 2;\r
+ EulerChar += LUT[n];\r
+ // Octant SEB\r
+ n = 1;\r
+ if( neighbors[8]==1 )\r
+ n |= 128;\r
+ if( neighbors[7]==1 )\r
+ n |= 64;\r
+ if( neighbors[17]==1 )\r
+ n |= 32;\r
+ if( neighbors[16]==1 )\r
+ n |= 16;\r
+ if( neighbors[5]==1 )\r
+ n |= 8;\r
+ if( neighbors[4]==1 )\r
+ n |= 4;\r
+ if( neighbors[14]==1 )\r
+ n |= 2;\r
+ EulerChar += LUT[n];\r
+ // Octant NWB\r
+ n = 1;\r
+ if( neighbors[0]==1 )\r
+ n |= 128;\r
+ if( neighbors[9]==1 )\r
+ n |= 64;\r
+ if( neighbors[3]==1 )\r
+ n |= 32;\r
+ if( neighbors[12]==1 )\r
+ n |= 16;\r
+ if( neighbors[1]==1 )\r
+ n |= 8;\r
+ if( neighbors[10]==1 )\r
+ n |= 4;\r
+ if( neighbors[4]==1 )\r
+ n |= 2;\r
+ EulerChar += LUT[n];\r
+ // Octant NEB\r
+ n = 1;\r
+ if( neighbors[2]==1 )\r
+ n |= 128;\r
+ if( neighbors[1]==1 )\r
+ n |= 64;\r
+ if( neighbors[11]==1 )\r
+ n |= 32;\r
+ if( neighbors[10]==1 )\r
+ n |= 16;\r
+ if( neighbors[5]==1 )\r
+ n |= 8;\r
+ if( neighbors[4]==1 )\r
+ n |= 4;\r
+ if( neighbors[14]==1 )\r
+ n |= 2;\r
+ EulerChar += LUT[n];\r
+ if( EulerChar == 0 )\r
+ return true;\r
+ else\r
+ return false;\r
+}\r
+\r
+/** \r
+ * Check if current point is a Simple Point.\r
+ * This method is named 'N(v)_labeling' in [Lee94].\r
+ * Outputs the number of connected objects in a neighborhood of a point\r
+ * after this point would have been removed.\r
+ */\r
+template <class TInputImage,class TOutputImage>\r
+bool \r
+BinaryThinningImageFilter3D<TInputImage,TOutputImage>\r
+::isSimplePoint(NeighborhoodType neighbors)\r
+{\r
+ // copy neighbors for labeling\r
+ int cube[26];\r
+ int i;\r
+ for( i = 0; i < 13; i++ ) // i = 0..12 -> cube[0..12]\r
+ cube[i] = neighbors[i];\r
+ // i != 13 : ignore center pixel when counting (see [Lee94])\r
+ for( i = 14; i < 27; i++ ) // i = 14..26 -> cube[13..25]\r
+ cube[i-1] = neighbors[i];\r
+ // set initial label\r
+ int label = 2;\r
+ // for all points in the neighborhood\r
+ for( int i = 0; i < 26; i++ )\r
+ {\r
+ if( cube[i]==1 ) // voxel has not been labelled yet\r
+ {\r
+ // start recursion with any octant that contains the point i\r
+ switch( i )\r
+ {\r
+ case 0:\r
+ case 1:\r
+ case 3:\r
+ case 4:\r
+ case 9:\r
+ case 10:\r
+ case 12:\r
+ Octree_labeling(1, label, cube );\r
+ break;\r
+ case 2:\r
+ case 5:\r
+ case 11:\r
+ case 13:\r
+ Octree_labeling(2, label, cube );\r
+ break;\r
+ case 6:\r
+ case 7:\r
+ case 14:\r
+ case 15:\r
+ Octree_labeling(3, label, cube );\r
+ break;\r
+ case 8:\r
+ case 16:\r
+ Octree_labeling(4, label, cube );\r
+ break;\r
+ case 17:\r
+ case 18:\r
+ case 20:\r
+ case 21:\r
+ Octree_labeling(5, label, cube );\r
+ break;\r
+ case 19:\r
+ case 22:\r
+ Octree_labeling(6, label, cube );\r
+ break;\r
+ case 23:\r
+ case 24:\r
+ Octree_labeling(7, label, cube );\r
+ break;\r
+ case 25:\r
+ Octree_labeling(8, label, cube );\r
+ break;\r
+ }\r
+ label++;\r
+ if( label-2 >= 2 )\r
+ {\r
+ return false;\r
+ }\r
+ }\r
+ }\r
+ //return label-2; in [Lee94] if the number of connected compontents would be needed\r
+ return true;\r
+}\r
+\r
+/** \r
+ * Octree_labeling [Lee94]\r
+ * This is a recursive method that calulates the number of connected\r
+ * components in the 3D neighbourhood after the center pixel would\r
+ * have been removed.\r
+ */\r
+template <class TInputImage,class TOutputImage>\r
+void \r
+BinaryThinningImageFilter3D<TInputImage,TOutputImage>\r
+::Octree_labeling(int octant, int label, int *cube)\r
+{\r
+ // check if there are points in the octant with value 1\r
+ if( octant==1 )\r
+ {\r
+ // set points in this octant to current label\r
+ // and recurseive labeling of adjacent octants\r
+ if( cube[0] == 1 )\r
+ cube[0] = label;\r
+ if( cube[1] == 1 )\r
+ {\r
+ cube[1] = label; \r
+ Octree_labeling( 2, label, cube);\r
+ }\r
+ if( cube[3] == 1 )\r
+ {\r
+ cube[3] = label; \r
+ Octree_labeling( 3, label, cube);\r
+ }\r
+ if( cube[4] == 1 )\r
+ {\r
+ cube[4] = label; \r
+ Octree_labeling( 2, label, cube);\r
+ Octree_labeling( 3, label, cube);\r
+ Octree_labeling( 4, label, cube);\r
+ }\r
+ if( cube[9] == 1 )\r
+ {\r
+ cube[9] = label; \r
+ Octree_labeling( 5, label, cube);\r
+ }\r
+ if( cube[10] == 1 )\r
+ {\r
+ cube[10] = label; \r
+ Octree_labeling( 2, label, cube);\r
+ Octree_labeling( 5, label, cube);\r
+ Octree_labeling( 6, label, cube);\r
+ }\r
+ if( cube[12] == 1 )\r
+ {\r
+ cube[12] = label; \r
+ Octree_labeling( 3, label, cube);\r
+ Octree_labeling( 5, label, cube);\r
+ Octree_labeling( 7, label, cube);\r
+ }\r
+ }\r
+ if( octant==2 )\r
+ {\r
+ if( cube[1] == 1 )\r
+ {\r
+ cube[1] = label;\r
+ Octree_labeling( 1, label, cube);\r
+ }\r
+ if( cube[4] == 1 )\r
+ {\r
+ cube[4] = label; \r
+ Octree_labeling( 1, label, cube);\r
+ Octree_labeling( 3, label, cube);\r
+ Octree_labeling( 4, label, cube);\r
+ }\r
+ if( cube[10] == 1 )\r
+ {\r
+ cube[10] = label; \r
+ Octree_labeling( 1, label, cube);\r
+ Octree_labeling( 5, label, cube);\r
+ Octree_labeling( 6, label, cube);\r
+ }\r
+ if( cube[2] == 1 )\r
+ cube[2] = label; \r
+ if( cube[5] == 1 )\r
+ {\r
+ cube[5] = label; \r
+ Octree_labeling( 4, label, cube);\r
+ }\r
+ if( cube[11] == 1 )\r
+ {\r
+ cube[11] = label; \r
+ Octree_labeling( 6, label, cube);\r
+ }\r
+ if( cube[13] == 1 )\r
+ {\r
+ cube[13] = label; \r
+ Octree_labeling( 4, label, cube);\r
+ Octree_labeling( 6, label, cube);\r
+ Octree_labeling( 8, label, cube);\r
+ }\r
+ }\r
+ if( octant==3 )\r
+ {\r
+ if( cube[3] == 1 )\r
+ {\r
+ cube[3] = label; \r
+ Octree_labeling( 1, label, cube);\r
+ }\r
+ if( cube[4] == 1 )\r
+ {\r
+ cube[4] = label; \r
+ Octree_labeling( 1, label, cube);\r
+ Octree_labeling( 2, label, cube);\r
+ Octree_labeling( 4, label, cube);\r
+ }\r
+ if( cube[12] == 1 )\r
+ {\r
+ cube[12] = label; \r
+ Octree_labeling( 1, label, cube);\r
+ Octree_labeling( 5, label, cube);\r
+ Octree_labeling( 7, label, cube);\r
+ }\r
+ if( cube[6] == 1 )\r
+ cube[6] = label; \r
+ if( cube[7] == 1 )\r
+ {\r
+ cube[7] = label; \r
+ Octree_labeling( 4, label, cube);\r
+ }\r
+ if( cube[14] == 1 )\r
+ {\r
+ cube[14] = label; \r
+ Octree_labeling( 7, label, cube);\r
+ }\r
+ if( cube[15] == 1 )\r
+ {\r
+ cube[15] = label; \r
+ Octree_labeling( 4, label, cube);\r
+ Octree_labeling( 7, label, cube);\r
+ Octree_labeling( 8, label, cube);\r
+ }\r
+ }\r
+ if( octant==4 )\r
+ {\r
+ if( cube[4] == 1 )\r
+ {\r
+ cube[4] = label; \r
+ Octree_labeling( 1, label, cube);\r
+ Octree_labeling( 2, label, cube);\r
+ Octree_labeling( 3, label, cube);\r
+ }\r
+ if( cube[5] == 1 )\r
+ {\r
+ cube[5] = label; \r
+ Octree_labeling( 2, label, cube);\r
+ }\r
+ if( cube[13] == 1 )\r
+ {\r
+ cube[13] = label; \r
+ Octree_labeling( 2, label, cube);\r
+ Octree_labeling( 6, label, cube);\r
+ Octree_labeling( 8, label, cube);\r
+ }\r
+ if( cube[7] == 1 )\r
+ {\r
+ cube[7] = label; \r
+ Octree_labeling( 3, label, cube);\r
+ }\r
+ if( cube[15] == 1 )\r
+ {\r
+ cube[15] = label; \r
+ Octree_labeling( 3, label, cube);\r
+ Octree_labeling( 7, label, cube);\r
+ Octree_labeling( 8, label, cube);\r
+ }\r
+ if( cube[8] == 1 )\r
+ cube[8] = label; \r
+ if( cube[16] == 1 )\r
+ {\r
+ cube[16] = label; \r
+ Octree_labeling( 8, label, cube);\r
+ }\r
+ }\r
+ if( octant==5 )\r
+ {\r
+ if( cube[9] == 1 )\r
+ {\r
+ cube[9] = label; \r
+ Octree_labeling( 1, label, cube);\r
+ }\r
+ if( cube[10] == 1 )\r
+ {\r
+ cube[10] = label; \r
+ Octree_labeling( 1, label, cube);\r
+ Octree_labeling( 2, label, cube);\r
+ Octree_labeling( 6, label, cube);\r
+ }\r
+ if( cube[12] == 1 )\r
+ {\r
+ cube[12] = label; \r
+ Octree_labeling( 1, label, cube);\r
+ Octree_labeling( 3, label, cube);\r
+ Octree_labeling( 7, label, cube);\r
+ }\r
+ if( cube[17] == 1 )\r
+ cube[17] = label; \r
+ if( cube[18] == 1 )\r
+ {\r
+ cube[18] = label; \r
+ Octree_labeling( 6, label, cube);\r
+ }\r
+ if( cube[20] == 1 )\r
+ {\r
+ cube[20] = label; \r
+ Octree_labeling( 7, label, cube);\r
+ }\r
+ if( cube[21] == 1 )\r
+ {\r
+ cube[21] = label; \r
+ Octree_labeling( 6, label, cube);\r
+ Octree_labeling( 7, label, cube);\r
+ Octree_labeling( 8, label, cube);\r
+ }\r
+ }\r
+ if( octant==6 )\r
+ {\r
+ if( cube[10] == 1 )\r
+ {\r
+ cube[10] = label; \r
+ Octree_labeling( 1, label, cube);\r
+ Octree_labeling( 2, label, cube);\r
+ Octree_labeling( 5, label, cube);\r
+ }\r
+ if( cube[11] == 1 )\r
+ {\r
+ cube[11] = label; \r
+ Octree_labeling( 2, label, cube);\r
+ }\r
+ if( cube[13] == 1 )\r
+ {\r
+ cube[13] = label; \r
+ Octree_labeling( 2, label, cube);\r
+ Octree_labeling( 4, label, cube);\r
+ Octree_labeling( 8, label, cube);\r
+ }\r
+ if( cube[18] == 1 )\r
+ {\r
+ cube[18] = label; \r
+ Octree_labeling( 5, label, cube);\r
+ }\r
+ if( cube[21] == 1 )\r
+ {\r
+ cube[21] = label; \r
+ Octree_labeling( 5, label, cube);\r
+ Octree_labeling( 7, label, cube);\r
+ Octree_labeling( 8, label, cube);\r
+ }\r
+ if( cube[19] == 1 )\r
+ cube[19] = label; \r
+ if( cube[22] == 1 )\r
+ {\r
+ cube[22] = label; \r
+ Octree_labeling( 8, label, cube);\r
+ }\r
+ }\r
+ if( octant==7 )\r
+ {\r
+ if( cube[12] == 1 )\r
+ {\r
+ cube[12] = label; \r
+ Octree_labeling( 1, label, cube);\r
+ Octree_labeling( 3, label, cube);\r
+ Octree_labeling( 5, label, cube);\r
+ }\r
+ if( cube[14] == 1 )\r
+ {\r
+ cube[14] = label; \r
+ Octree_labeling( 3, label, cube);\r
+ }\r
+ if( cube[15] == 1 )\r
+ {\r
+ cube[15] = label; \r
+ Octree_labeling( 3, label, cube);\r
+ Octree_labeling( 4, label, cube);\r
+ Octree_labeling( 8, label, cube);\r
+ }\r
+ if( cube[20] == 1 )\r
+ {\r
+ cube[20] = label; \r
+ Octree_labeling( 5, label, cube);\r
+ }\r
+ if( cube[21] == 1 )\r
+ {\r
+ cube[21] = label; \r
+ Octree_labeling( 5, label, cube);\r
+ Octree_labeling( 6, label, cube);\r
+ Octree_labeling( 8, label, cube);\r
+ }\r
+ if( cube[23] == 1 )\r
+ cube[23] = label; \r
+ if( cube[24] == 1 )\r
+ {\r
+ cube[24] = label; \r
+ Octree_labeling( 8, label, cube);\r
+ }\r
+ }\r
+ if( octant==8 )\r
+ {\r
+ if( cube[13] == 1 )\r
+ {\r
+ cube[13] = label; \r
+ Octree_labeling( 2, label, cube);\r
+ Octree_labeling( 4, label, cube);\r
+ Octree_labeling( 6, label, cube);\r
+ }\r
+ if( cube[15] == 1 )\r
+ {\r
+ cube[15] = label; \r
+ Octree_labeling( 3, label, cube);\r
+ Octree_labeling( 4, label, cube);\r
+ Octree_labeling( 7, label, cube);\r
+ }\r
+ if( cube[16] == 1 )\r
+ {\r
+ cube[16] = label; \r
+ Octree_labeling( 4, label, cube);\r
+ }\r
+ if( cube[21] == 1 )\r
+ {\r
+ cube[21] = label; \r
+ Octree_labeling( 5, label, cube);\r
+ Octree_labeling( 6, label, cube);\r
+ Octree_labeling( 7, label, cube);\r
+ }\r
+ if( cube[22] == 1 )\r
+ {\r
+ cube[22] = label; \r
+ Octree_labeling( 6, label, cube);\r
+ }\r
+ if( cube[24] == 1 )\r
+ {\r
+ cube[24] = label; \r
+ Octree_labeling( 7, label, cube);\r
+ }\r
+ if( cube[25] == 1 )\r
+ cube[25] = label; \r
+ } \r
+}\r
+\r
+\r
+/**\r
+ * Print Self\r
+ */\r
+template <class TInputImage,class TOutputImage>\r
+void \r
+BinaryThinningImageFilter3D<TInputImage,TOutputImage>\r
+::PrintSelf(std::ostream& os, Indent indent) const\r
+{\r
+ Superclass::PrintSelf(os,indent);\r
+ \r
+ os << indent << "Thinning image: " << std::endl;\r
+\r
+}\r
+\r
+} // end namespace itk\r
+\r
+#endif\r
git://git.creatis.insa-lyon.fr / clitk.git / commitdiff