-#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
+#ifndef _itkBinaryThinningImageFilter3D_txx
+#define _itkBinaryThinningImageFilter3D_txx
+
+#include <iostream>
+
+#include "itkBinaryThinningImageFilter3D.h"
+#include "itkImageRegionConstIterator.h"
+#include "itkImageRegionIterator.h"
+#include "itkNeighborhoodIterator.h"
+#include <vector>
+
+namespace itk
+{
+
+/**
+ * Constructor
+ */
+template <class TInputImage,class TOutputImage>
+BinaryThinningImageFilter3D<TInputImage,TOutputImage>
+::BinaryThinningImageFilter3D()
+{
+
+ this->SetNumberOfRequiredOutputs( 1 );
+
+ OutputImagePointer thinImage = OutputImageType::New();
+ this->SetNthOutput( 0, thinImage.GetPointer() );
+
+}
+
+/**
+ * Return the thinning Image pointer
+ */
+template <class TInputImage,class TOutputImage>
+typename BinaryThinningImageFilter3D<
+TInputImage,TOutputImage>::OutputImageType *
+BinaryThinningImageFilter3D<TInputImage,TOutputImage>
+::GetThinning(void)
+{
+ return dynamic_cast< OutputImageType * >(
+ this->ProcessObject::GetOutput(0) );
+}
+
+
+/**
+ * Prepare data for computation
+ * Copy the input image to the output image, changing from the input
+ * type to the output type.
+ */
+template <class TInputImage,class TOutputImage>
+void
+BinaryThinningImageFilter3D<TInputImage,TOutputImage>
+::PrepareData(void)
+{
+
+ itkDebugMacro(<< "PrepareData Start");
+ OutputImagePointer thinImage = GetThinning();
+
+ InputImagePointer inputImage =
+ dynamic_cast<const TInputImage *>( ProcessObject::GetInput(0) );
+
+ thinImage->SetBufferedRegion( thinImage->GetRequestedRegion() );
+ thinImage->Allocate();
+
+ typename OutputImageType::RegionType region = thinImage->GetRequestedRegion();
+
+
+ ImageRegionConstIterator< TInputImage > it( inputImage, region );
+ ImageRegionIterator< TOutputImage > ot( thinImage, region );
+
+ it.GoToBegin();
+ ot.GoToBegin();
+
+ itkDebugMacro(<< "PrepareData: Copy input to output");
+
+ // Copy the input to the output, changing all foreground pixels to
+ // have value 1 in the process.
+ while( !ot.IsAtEnd() ) {
+ if ( it.Get() ) {
+ ot.Set( NumericTraits<OutputImagePixelType>::One );
+ } else {
+ ot.Set( NumericTraits<OutputImagePixelType>::Zero );
+ }
+ ++it;
+ ++ot;
+ }
+ itkDebugMacro(<< "PrepareData End");
+}
+
+/**
+ * Post processing for computing thinning
+ */
+template <class TInputImage,class TOutputImage>
+void
+BinaryThinningImageFilter3D<TInputImage,TOutputImage>
+::ComputeThinImage()
+{
+ itkDebugMacro( << "ComputeThinImage Start");
+ OutputImagePointer thinImage = GetThinning();
+
+ typename OutputImageType::RegionType region = thinImage->GetRequestedRegion();
+
+ ConstBoundaryConditionType boundaryCondition;
+ boundaryCondition.SetConstant( 0 );
+
+ typename NeighborhoodIteratorType::RadiusType radius;
+ radius.Fill(1);
+ NeighborhoodIteratorType ot( radius, thinImage, region );
+ ot.SetBoundaryCondition( boundaryCondition );
+
+ std::vector < IndexType > simpleBorderPoints;
+ typename std::vector < IndexType >::iterator simpleBorderPointsIt;
+
+ // Define offsets
+ typedef typename NeighborhoodIteratorType::OffsetType OffsetType;
+ OffsetType N = {{ 0,-1, 0}}; // north
+ OffsetType S = {{ 0, 1, 0}}; // south
+ OffsetType E = {{ 1, 0, 0}}; // east
+ OffsetType W = {{-1, 0, 0}}; // west
+ OffsetType U = {{ 0, 0, 1}}; // up
+ OffsetType B = {{ 0, 0,-1}}; // bottom
+
+ // prepare Euler LUT [Lee94]
+ int eulerLUT[256];
+ fillEulerLUT( eulerLUT );
+ // Loop through the image several times until there is no change.
+ int unchangedBorders = 0;
+ while( unchangedBorders < 6 ) { // loop until no change for all the six border types
+ unchangedBorders = 0;
+ for( int currentBorder = 1; currentBorder <= 6; currentBorder++) {
+ // Loop through the image.
+ for ( ot.GoToBegin(); !ot.IsAtEnd(); ++ot ) {
+ // check if point is foreground
+ if ( ot.GetCenterPixel() != 1 ) {
+ continue; // current point is already background
+ }
+ // check 6-neighbors if point is a border point of type currentBorder
+ bool isBorderPoint = false;
+ if( currentBorder == 1 && ot.GetPixel(N)<=0 )
+ isBorderPoint = true;
+ if( currentBorder == 2 && ot.GetPixel(S)<=0 )
+ isBorderPoint = true;
+ if( currentBorder == 3 && ot.GetPixel(E)<=0 )
+ isBorderPoint = true;
+ if( currentBorder == 4 && ot.GetPixel(W)<=0 )
+ isBorderPoint = true;
+ if( currentBorder == 5 && ot.GetPixel(U)<=0 )
+ isBorderPoint = true;
+ if( currentBorder == 6 && ot.GetPixel(B)<=0 )
+ isBorderPoint = true;
+ if( !isBorderPoint ) {
+ continue; // current point is not deletable
+ }
+ // check if point is the end of an arc
+ int numberOfNeighbors = -1; // -1 and not 0 because the center pixel will be counted as well
+ for( int i = 0; i < 27; i++ ) // i = 0..26
+ if( ot.GetPixel(i)==1 )
+ numberOfNeighbors++;
+
+ if( numberOfNeighbors == 1 ) {
+ continue; // current point is not deletable
+ }
+
+ // check if point is Euler invariant
+ if( !isEulerInvariant( ot.GetNeighborhood(), eulerLUT ) ) {
+ continue; // current point is not deletable
+ }
+
+ // check if point is simple (deletion does not change connectivity in the 3x3x3 neighborhood)
+ if( !isSimplePoint( ot.GetNeighborhood() ) ) {
+ continue; // current point is not deletable
+ }
+
+ // add all simple border points to a list for sequential re-checking
+ simpleBorderPoints.push_back( ot.GetIndex() );
+ } // end image iteration loop
+
+ // sequential re-checking to preserve connectivity when
+ // deleting in a parallel way
+ bool noChange = true;
+ for( simpleBorderPointsIt=simpleBorderPoints.begin(); simpleBorderPointsIt!=simpleBorderPoints.end(); simpleBorderPointsIt++) {
+ // 1. Set simple border point to 0
+ thinImage->SetPixel( *simpleBorderPointsIt, NumericTraits<OutputImagePixelType>::Zero);
+ // 2. Check if neighborhood is still connected
+ ot.SetLocation( *simpleBorderPointsIt );
+ if( !isSimplePoint( ot.GetNeighborhood() ) ) {
+ // we cannot delete current point, so reset
+ thinImage->SetPixel( *simpleBorderPointsIt, NumericTraits<OutputImagePixelType>::One );
+ } else {
+ noChange = false;
+ }
+ }
+ if( noChange )
+ unchangedBorders++;
+
+ simpleBorderPoints.clear();
+ } // end currentBorder for loop
+ } // end unchangedBorders while loop
+
+ itkDebugMacro( << "ComputeThinImage End");
+}
+
+/**
+ * Generate ThinImage
+ */
+template <class TInputImage,class TOutputImage>
+void
+BinaryThinningImageFilter3D<TInputImage,TOutputImage>
+::GenerateData()
+{
+
+ this->PrepareData();
+
+ itkDebugMacro(<< "GenerateData: Computing Thinning Image");
+ this->ComputeThinImage();
+} // end GenerateData()
+
+/**
+ * Fill the Euler look-up table (LUT) for later check of the Euler invariance. (see [Lee94])
+ */
+template <class TInputImage,class TOutputImage>
+void
+BinaryThinningImageFilter3D<TInputImage,TOutputImage>
+::fillEulerLUT(int *LUT)
+{
+ LUT[1] = 1;
+ LUT[3] = -1;
+ LUT[5] = -1;
+ LUT[7] = 1;
+ LUT[9] = -3;
+ LUT[11] = -1;
+ LUT[13] = -1;
+ LUT[15] = 1;
+ LUT[17] = -1;
+ LUT[19] = 1;
+ LUT[21] = 1;
+ LUT[23] = -1;
+ LUT[25] = 3;
+ LUT[27] = 1;
+ LUT[29] = 1;
+ LUT[31] = -1;
+ LUT[33] = -3;
+ LUT[35] = -1;
+ LUT[37] = 3;
+ LUT[39] = 1;
+ LUT[41] = 1;
+ LUT[43] = -1;
+ LUT[45] = 3;
+ LUT[47] = 1;
+ LUT[49] = -1;
+ LUT[51] = 1;
+
+ LUT[53] = 1;
+ LUT[55] = -1;
+ LUT[57] = 3;
+ LUT[59] = 1;
+ LUT[61] = 1;
+ LUT[63] = -1;
+ LUT[65] = -3;
+ LUT[67] = 3;
+ LUT[69] = -1;
+ LUT[71] = 1;
+ LUT[73] = 1;
+ LUT[75] = 3;
+ LUT[77] = -1;
+ LUT[79] = 1;
+ LUT[81] = -1;
+ LUT[83] = 1;
+ LUT[85] = 1;
+ LUT[87] = -1;
+ LUT[89] = 3;
+ LUT[91] = 1;
+ LUT[93] = 1;
+ LUT[95] = -1;
+ LUT[97] = 1;
+ LUT[99] = 3;
+ LUT[101] = 3;
+ LUT[103] = 1;
+
+ LUT[105] = 5;
+ LUT[107] = 3;
+ LUT[109] = 3;
+ LUT[111] = 1;
+ LUT[113] = -1;
+ LUT[115] = 1;
+ LUT[117] = 1;
+ LUT[119] = -1;
+ LUT[121] = 3;
+ LUT[123] = 1;
+ LUT[125] = 1;
+ LUT[127] = -1;
+ LUT[129] = -7;
+ LUT[131] = -1;
+ LUT[133] = -1;
+ LUT[135] = 1;
+ LUT[137] = -3;
+ LUT[139] = -1;
+ LUT[141] = -1;
+ LUT[143] = 1;
+ LUT[145] = -1;
+ LUT[147] = 1;
+ LUT[149] = 1;
+ LUT[151] = -1;
+ LUT[153] = 3;
+ LUT[155] = 1;
+
+ LUT[157] = 1;
+ LUT[159] = -1;
+ LUT[161] = -3;
+ LUT[163] = -1;
+ LUT[165] = 3;
+ LUT[167] = 1;
+ LUT[169] = 1;
+ LUT[171] = -1;
+ LUT[173] = 3;
+ LUT[175] = 1;
+ LUT[177] = -1;
+ LUT[179] = 1;
+ LUT[181] = 1;
+ LUT[183] = -1;
+ LUT[185] = 3;
+ LUT[187] = 1;
+ LUT[189] = 1;
+ LUT[191] = -1;
+ LUT[193] = -3;
+ LUT[195] = 3;
+ LUT[197] = -1;
+ LUT[199] = 1;
+ LUT[201] = 1;
+ LUT[203] = 3;
+ LUT[205] = -1;
+ LUT[207] = 1;
+
+ LUT[209] = -1;
+ LUT[211] = 1;
+ LUT[213] = 1;
+ LUT[215] = -1;
+ LUT[217] = 3;
+ LUT[219] = 1;
+ LUT[221] = 1;
+ LUT[223] = -1;
+ LUT[225] = 1;
+ LUT[227] = 3;
+ LUT[229] = 3;
+ LUT[231] = 1;
+ LUT[233] = 5;
+ LUT[235] = 3;
+ LUT[237] = 3;
+ LUT[239] = 1;
+ LUT[241] = -1;
+ LUT[243] = 1;
+ LUT[245] = 1;
+ LUT[247] = -1;
+ LUT[249] = 3;
+ LUT[251] = 1;
+ LUT[253] = 1;
+ LUT[255] = -1;
+}
+
+/**
+ * Check for Euler invariance. (see [Lee94])
+ */
+template <class TInputImage,class TOutputImage>
+bool
+BinaryThinningImageFilter3D<TInputImage,TOutputImage>
+::isEulerInvariant(NeighborhoodType neighbors, int *LUT)
+{
+ // calculate Euler characteristic for each octant and sum up
+ int EulerChar = 0;
+ unsigned char n;
+ // Octant SWU
+ n = 1;
+ if( neighbors[24]==1 )
+ n |= 128;
+ if( neighbors[25]==1 )
+ n |= 64;
+ if( neighbors[15]==1 )
+ n |= 32;
+ if( neighbors[16]==1 )
+ n |= 16;
+ if( neighbors[21]==1 )
+ n |= 8;
+ if( neighbors[22]==1 )
+ n |= 4;
+ if( neighbors[12]==1 )
+ n |= 2;
+ EulerChar += LUT[n];
+ // Octant SEU
+ n = 1;
+ if( neighbors[26]==1 )
+ n |= 128;
+ if( neighbors[23]==1 )
+ n |= 64;
+ if( neighbors[17]==1 )
+ n |= 32;
+ if( neighbors[14]==1 )
+ n |= 16;
+ if( neighbors[25]==1 )
+ n |= 8;
+ if( neighbors[22]==1 )
+ n |= 4;
+ if( neighbors[16]==1 )
+ n |= 2;
+ EulerChar += LUT[n];
+ // Octant NWU
+ n = 1;
+ if( neighbors[18]==1 )
+ n |= 128;
+ if( neighbors[21]==1 )
+ n |= 64;
+ if( neighbors[9]==1 )
+ n |= 32;
+ if( neighbors[12]==1 )
+ n |= 16;
+ if( neighbors[19]==1 )
+ n |= 8;
+ if( neighbors[22]==1 )
+ n |= 4;
+ if( neighbors[10]==1 )
+ n |= 2;
+ EulerChar += LUT[n];
+ // Octant NEU
+ n = 1;
+ if( neighbors[20]==1 )
+ n |= 128;
+ if( neighbors[23]==1 )
+ n |= 64;
+ if( neighbors[19]==1 )
+ n |= 32;
+ if( neighbors[22]==1 )
+ n |= 16;
+ if( neighbors[11]==1 )
+ n |= 8;
+ if( neighbors[14]==1 )
+ n |= 4;
+ if( neighbors[10]==1 )
+ n |= 2;
+ EulerChar += LUT[n];
+ // Octant SWB
+ n = 1;
+ if( neighbors[6]==1 )
+ n |= 128;
+ if( neighbors[15]==1 )
+ n |= 64;
+ if( neighbors[7]==1 )
+ n |= 32;
+ if( neighbors[16]==1 )
+ n |= 16;
+ if( neighbors[3]==1 )
+ n |= 8;
+ if( neighbors[12]==1 )
+ n |= 4;
+ if( neighbors[4]==1 )
+ n |= 2;
+ EulerChar += LUT[n];
+ // Octant SEB
+ n = 1;
+ if( neighbors[8]==1 )
+ n |= 128;
+ if( neighbors[7]==1 )
+ n |= 64;
+ if( neighbors[17]==1 )
+ n |= 32;
+ if( neighbors[16]==1 )
+ n |= 16;
+ if( neighbors[5]==1 )
+ n |= 8;
+ if( neighbors[4]==1 )
+ n |= 4;
+ if( neighbors[14]==1 )
+ n |= 2;
+ EulerChar += LUT[n];
+ // Octant NWB
+ n = 1;
+ if( neighbors[0]==1 )
+ n |= 128;
+ if( neighbors[9]==1 )
+ n |= 64;
+ if( neighbors[3]==1 )
+ n |= 32;
+ if( neighbors[12]==1 )
+ n |= 16;
+ if( neighbors[1]==1 )
+ n |= 8;
+ if( neighbors[10]==1 )
+ n |= 4;
+ if( neighbors[4]==1 )
+ n |= 2;
+ EulerChar += LUT[n];
+ // Octant NEB
+ n = 1;
+ if( neighbors[2]==1 )
+ n |= 128;
+ if( neighbors[1]==1 )
+ n |= 64;
+ if( neighbors[11]==1 )
+ n |= 32;
+ if( neighbors[10]==1 )
+ n |= 16;
+ if( neighbors[5]==1 )
+ n |= 8;
+ if( neighbors[4]==1 )
+ n |= 4;
+ if( neighbors[14]==1 )
+ n |= 2;
+ EulerChar += LUT[n];
+ if( EulerChar == 0 )
+ return true;
+ else
+ return false;
+}
+
+/**
+ * Check if current point is a Simple Point.
+ * This method is named 'N(v)_labeling' in [Lee94].
+ * Outputs the number of connected objects in a neighborhood of a point
+ * after this point would have been removed.
+ */
+template <class TInputImage,class TOutputImage>
+bool
+BinaryThinningImageFilter3D<TInputImage,TOutputImage>
+::isSimplePoint(NeighborhoodType neighbors)
+{
+ // copy neighbors for labeling
+ int cube[26];
+ int i;
+ for( i = 0; i < 13; i++ ) // i = 0..12 -> cube[0..12]
+ cube[i] = neighbors[i];
+ // i != 13 : ignore center pixel when counting (see [Lee94])
+ for( i = 14; i < 27; i++ ) // i = 14..26 -> cube[13..25]
+ cube[i-1] = neighbors[i];
+ // set initial label
+ int label = 2;
+ // for all points in the neighborhood
+ for( int i = 0; i < 26; i++ ) {
+ if( cube[i]==1 ) { // voxel has not been labelled yet
+ // start recursion with any octant that contains the point i
+ switch( i ) {
+ case 0:
+ case 1:
+ case 3:
+ case 4:
+ case 9:
+ case 10:
+ case 12:
+ Octree_labeling(1, label, cube );
+ break;
+ case 2:
+ case 5:
+ case 11:
+ case 13:
+ Octree_labeling(2, label, cube );
+ break;
+ case 6:
+ case 7:
+ case 14:
+ case 15:
+ Octree_labeling(3, label, cube );
+ break;
+ case 8:
+ case 16:
+ Octree_labeling(4, label, cube );
+ break;
+ case 17:
+ case 18:
+ case 20:
+ case 21:
+ Octree_labeling(5, label, cube );
+ break;
+ case 19:
+ case 22:
+ Octree_labeling(6, label, cube );
+ break;
+ case 23:
+ case 24:
+ Octree_labeling(7, label, cube );
+ break;
+ case 25:
+ Octree_labeling(8, label, cube );
+ break;
+ }
+ label++;
+ if( label-2 >= 2 ) {
+ return false;
+ }
+ }
+ }
+ //return label-2; in [Lee94] if the number of connected compontents would be needed
+ return true;
+}
+
+/**
+ * Octree_labeling [Lee94]
+ * This is a recursive method that calulates the number of connected
+ * components in the 3D neighbourhood after the center pixel would
+ * have been removed.
+ */
+template <class TInputImage,class TOutputImage>
+void
+BinaryThinningImageFilter3D<TInputImage,TOutputImage>
+::Octree_labeling(int octant, int label, int *cube)
+{
+ // check if there are points in the octant with value 1
+ if( octant==1 ) {
+ // set points in this octant to current label
+ // and recurseive labeling of adjacent octants
+ if( cube[0] == 1 )
+ cube[0] = label;
+ if( cube[1] == 1 ) {
+ cube[1] = label;
+ Octree_labeling( 2, label, cube);
+ }
+ if( cube[3] == 1 ) {
+ cube[3] = label;
+ Octree_labeling( 3, label, cube);
+ }
+ if( cube[4] == 1 ) {
+ cube[4] = label;
+ Octree_labeling( 2, label, cube);
+ Octree_labeling( 3, label, cube);
+ Octree_labeling( 4, label, cube);
+ }
+ if( cube[9] == 1 ) {
+ cube[9] = label;
+ Octree_labeling( 5, label, cube);
+ }
+ if( cube[10] == 1 ) {
+ cube[10] = label;
+ Octree_labeling( 2, label, cube);
+ Octree_labeling( 5, label, cube);
+ Octree_labeling( 6, label, cube);
+ }
+ if( cube[12] == 1 ) {
+ cube[12] = label;
+ Octree_labeling( 3, label, cube);
+ Octree_labeling( 5, label, cube);
+ Octree_labeling( 7, label, cube);
+ }
+ }
+ if( octant==2 ) {
+ if( cube[1] == 1 ) {
+ cube[1] = label;
+ Octree_labeling( 1, label, cube);
+ }
+ if( cube[4] == 1 ) {
+ cube[4] = label;
+ Octree_labeling( 1, label, cube);
+ Octree_labeling( 3, label, cube);
+ Octree_labeling( 4, label, cube);
+ }
+ if( cube[10] == 1 ) {
+ cube[10] = label;
+ Octree_labeling( 1, label, cube);
+ Octree_labeling( 5, label, cube);
+ Octree_labeling( 6, label, cube);
+ }
+ if( cube[2] == 1 )
+ cube[2] = label;
+ if( cube[5] == 1 ) {
+ cube[5] = label;
+ Octree_labeling( 4, label, cube);
+ }
+ if( cube[11] == 1 ) {
+ cube[11] = label;
+ Octree_labeling( 6, label, cube);
+ }
+ if( cube[13] == 1 ) {
+ cube[13] = label;
+ Octree_labeling( 4, label, cube);
+ Octree_labeling( 6, label, cube);
+ Octree_labeling( 8, label, cube);
+ }
+ }
+ if( octant==3 ) {
+ if( cube[3] == 1 ) {
+ cube[3] = label;
+ Octree_labeling( 1, label, cube);
+ }
+ if( cube[4] == 1 ) {
+ cube[4] = label;
+ Octree_labeling( 1, label, cube);
+ Octree_labeling( 2, label, cube);
+ Octree_labeling( 4, label, cube);
+ }
+ if( cube[12] == 1 ) {
+ cube[12] = label;
+ Octree_labeling( 1, label, cube);
+ Octree_labeling( 5, label, cube);
+ Octree_labeling( 7, label, cube);
+ }
+ if( cube[6] == 1 )
+ cube[6] = label;
+ if( cube[7] == 1 ) {
+ cube[7] = label;
+ Octree_labeling( 4, label, cube);
+ }
+ if( cube[14] == 1 ) {
+ cube[14] = label;
+ Octree_labeling( 7, label, cube);
+ }
+ if( cube[15] == 1 ) {
+ cube[15] = label;
+ Octree_labeling( 4, label, cube);
+ Octree_labeling( 7, label, cube);
+ Octree_labeling( 8, label, cube);
+ }
+ }
+ if( octant==4 ) {
+ if( cube[4] == 1 ) {
+ cube[4] = label;
+ Octree_labeling( 1, label, cube);
+ Octree_labeling( 2, label, cube);
+ Octree_labeling( 3, label, cube);
+ }
+ if( cube[5] == 1 ) {
+ cube[5] = label;
+ Octree_labeling( 2, label, cube);
+ }
+ if( cube[13] == 1 ) {
+ cube[13] = label;
+ Octree_labeling( 2, label, cube);
+ Octree_labeling( 6, label, cube);
+ Octree_labeling( 8, label, cube);
+ }
+ if( cube[7] == 1 ) {
+ cube[7] = label;
+ Octree_labeling( 3, label, cube);
+ }
+ if( cube[15] == 1 ) {
+ cube[15] = label;
+ Octree_labeling( 3, label, cube);
+ Octree_labeling( 7, label, cube);
+ Octree_labeling( 8, label, cube);
+ }
+ if( cube[8] == 1 )
+ cube[8] = label;
+ if( cube[16] == 1 ) {
+ cube[16] = label;
+ Octree_labeling( 8, label, cube);
+ }
+ }
+ if( octant==5 ) {
+ if( cube[9] == 1 ) {
+ cube[9] = label;
+ Octree_labeling( 1, label, cube);
+ }
+ if( cube[10] == 1 ) {
+ cube[10] = label;
+ Octree_labeling( 1, label, cube);
+ Octree_labeling( 2, label, cube);
+ Octree_labeling( 6, label, cube);
+ }
+ if( cube[12] == 1 ) {
+ cube[12] = label;
+ Octree_labeling( 1, label, cube);
+ Octree_labeling( 3, label, cube);
+ Octree_labeling( 7, label, cube);
+ }
+ if( cube[17] == 1 )
+ cube[17] = label;
+ if( cube[18] == 1 ) {
+ cube[18] = label;
+ Octree_labeling( 6, label, cube);
+ }
+ if( cube[20] == 1 ) {
+ cube[20] = label;
+ Octree_labeling( 7, label, cube);
+ }
+ if( cube[21] == 1 ) {
+ cube[21] = label;
+ Octree_labeling( 6, label, cube);
+ Octree_labeling( 7, label, cube);
+ Octree_labeling( 8, label, cube);
+ }
+ }
+ if( octant==6 ) {
+ if( cube[10] == 1 ) {
+ cube[10] = label;
+ Octree_labeling( 1, label, cube);
+ Octree_labeling( 2, label, cube);
+ Octree_labeling( 5, label, cube);
+ }
+ if( cube[11] == 1 ) {
+ cube[11] = label;
+ Octree_labeling( 2, label, cube);
+ }
+ if( cube[13] == 1 ) {
+ cube[13] = label;
+ Octree_labeling( 2, label, cube);
+ Octree_labeling( 4, label, cube);
+ Octree_labeling( 8, label, cube);
+ }
+ if( cube[18] == 1 ) {
+ cube[18] = label;
+ Octree_labeling( 5, label, cube);
+ }
+ if( cube[21] == 1 ) {
+ cube[21] = label;
+ Octree_labeling( 5, label, cube);
+ Octree_labeling( 7, label, cube);
+ Octree_labeling( 8, label, cube);
+ }
+ if( cube[19] == 1 )
+ cube[19] = label;
+ if( cube[22] == 1 ) {
+ cube[22] = label;
+ Octree_labeling( 8, label, cube);
+ }
+ }
+ if( octant==7 ) {
+ if( cube[12] == 1 ) {
+ cube[12] = label;
+ Octree_labeling( 1, label, cube);
+ Octree_labeling( 3, label, cube);
+ Octree_labeling( 5, label, cube);
+ }
+ if( cube[14] == 1 ) {
+ cube[14] = label;
+ Octree_labeling( 3, label, cube);
+ }
+ if( cube[15] == 1 ) {
+ cube[15] = label;
+ Octree_labeling( 3, label, cube);
+ Octree_labeling( 4, label, cube);
+ Octree_labeling( 8, label, cube);
+ }
+ if( cube[20] == 1 ) {
+ cube[20] = label;
+ Octree_labeling( 5, label, cube);
+ }
+ if( cube[21] == 1 ) {
+ cube[21] = label;
+ Octree_labeling( 5, label, cube);
+ Octree_labeling( 6, label, cube);
+ Octree_labeling( 8, label, cube);
+ }
+ if( cube[23] == 1 )
+ cube[23] = label;
+ if( cube[24] == 1 ) {
+ cube[24] = label;
+ Octree_labeling( 8, label, cube);
+ }
+ }
+ if( octant==8 ) {
+ if( cube[13] == 1 ) {
+ cube[13] = label;
+ Octree_labeling( 2, label, cube);
+ Octree_labeling( 4, label, cube);
+ Octree_labeling( 6, label, cube);
+ }
+ if( cube[15] == 1 ) {
+ cube[15] = label;
+ Octree_labeling( 3, label, cube);
+ Octree_labeling( 4, label, cube);
+ Octree_labeling( 7, label, cube);
+ }
+ if( cube[16] == 1 ) {
+ cube[16] = label;
+ Octree_labeling( 4, label, cube);
+ }
+ if( cube[21] == 1 ) {
+ cube[21] = label;
+ Octree_labeling( 5, label, cube);
+ Octree_labeling( 6, label, cube);
+ Octree_labeling( 7, label, cube);
+ }
+ if( cube[22] == 1 ) {
+ cube[22] = label;
+ Octree_labeling( 6, label, cube);
+ }
+ if( cube[24] == 1 ) {
+ cube[24] = label;
+ Octree_labeling( 7, label, cube);
+ }
+ if( cube[25] == 1 )
+ cube[25] = label;
+ }
+}
+
+
+/**
+ * Print Self
+ */
+template <class TInputImage,class TOutputImage>
+void
+BinaryThinningImageFilter3D<TInputImage,TOutputImage>
+::PrintSelf(std::ostream& os, Indent indent) const
+{
+ Superclass::PrintSelf(os,indent);
+
+ os << indent << "Thinning image: " << std::endl;
+
+}
+
+} // end namespace itk
+
+#endif
git://git.creatis.insa-lyon.fr / clitk.git / commitdiff