itk/itkBinaryThinningImageFilter3D.txx

   1 #ifndef _itkBinaryThinningImageFilter3D_txx\r
   2 #define _itkBinaryThinningImageFilter3D_txx\r
   3 \r
   4 #include <iostream>\r
   5 \r
   6 #include "itkBinaryThinningImageFilter3D.h"\r
   7 #include "itkImageRegionConstIterator.h"\r
   8 #include "itkImageRegionIterator.h"\r
   9 #include "itkNeighborhoodIterator.h"\r
  10 #include <vector>\r
  11 \r
  12 namespace itk\r
  13 {\r
  14 \r
  15 /**\r
  16  *    Constructor\r
  17  */\r
  18 template <class TInputImage,class TOutputImage>\r
  19 BinaryThinningImageFilter3D<TInputImage,TOutputImage>\r
  20 ::BinaryThinningImageFilter3D()\r
  21 {\r
  22 \r
  23   this->SetNumberOfRequiredOutputs( 1 );\r
  24 \r
  25   OutputImagePointer thinImage = OutputImageType::New();\r
  26   this->SetNthOutput( 0, thinImage.GetPointer() );\r
  27 \r
  28 }\r
  29 \r
  30 /**\r
  31  *  Return the thinning Image pointer\r
  32  */\r
  33 template <class TInputImage,class TOutputImage>\r
  34 typename BinaryThinningImageFilter3D<\r
  35   TInputImage,TOutputImage>::OutputImageType * \r
  36 BinaryThinningImageFilter3D<TInputImage,TOutputImage>\r
  37 ::GetThinning(void)\r
  38 {\r
  39   return  dynamic_cast< OutputImageType * >(\r
  40     this->ProcessObject::GetOutput(0) );\r
  41 }\r
  42 \r
  43 \r
  44 /**\r
  45  *  Prepare data for computation\r
  46  *  Copy the input image to the output image, changing from the input\r
  47  *  type to the output type.\r
  48  */\r
  49 template <class TInputImage,class TOutputImage>\r
  50 void \r
  51 BinaryThinningImageFilter3D<TInputImage,TOutputImage>\r
  52 ::PrepareData(void) \r
  53 {\r
  54   \r
  55   itkDebugMacro(<< "PrepareData Start");\r
  56   OutputImagePointer thinImage = GetThinning();\r
  57 \r
  58   InputImagePointer  inputImage  = \r
  59     dynamic_cast<const TInputImage  *>( ProcessObject::GetInput(0) );\r
  60 \r
  61   thinImage->SetBufferedRegion( thinImage->GetRequestedRegion() );\r
  62   thinImage->Allocate();\r
  63 \r
  64   typename OutputImageType::RegionType region  = thinImage->GetRequestedRegion();\r
  65 \r
  66 \r
  67   ImageRegionConstIterator< TInputImage >  it( inputImage,  region );\r
  68   ImageRegionIterator< TOutputImage > ot( thinImage,  region );\r
  69 \r
  70   it.GoToBegin();\r
  71   ot.GoToBegin();\r
  72 \r
  73   itkDebugMacro(<< "PrepareData: Copy input to output");\r
  74  \r
  75   // Copy the input to the output, changing all foreground pixels to\r
  76   // have value 1 in the process.\r
  77   while( !ot.IsAtEnd() )\r
  78       {\r
  79       if ( it.Get() )\r
  80         {\r
  81         ot.Set( NumericTraits<OutputImagePixelType>::One );\r
  82         }\r
  83       else\r
  84         {\r
  85         ot.Set( NumericTraits<OutputImagePixelType>::Zero );\r
  86         }\r
  87       ++it;\r
  88       ++ot;\r
  89       }\r
  90   itkDebugMacro(<< "PrepareData End");    \r
  91 }\r
  92 \r
  93 /**\r
  94  *  Post processing for computing thinning\r
  95  */\r
  96 template <class TInputImage,class TOutputImage>\r
  97 void \r
  98 BinaryThinningImageFilter3D<TInputImage,TOutputImage>\r
  99 ::ComputeThinImage() \r
 100 {\r
 101   itkDebugMacro( << "ComputeThinImage Start");\r
 102   OutputImagePointer thinImage = GetThinning();\r
 103 \r
 104   typename OutputImageType::RegionType region = thinImage->GetRequestedRegion();\r
 105   \r
 106   ConstBoundaryConditionType boundaryCondition;\r
 107   boundaryCondition.SetConstant( 0 );\r
 108 \r
 109   typename NeighborhoodIteratorType::RadiusType radius;\r
 110   radius.Fill(1);\r
 111   NeighborhoodIteratorType ot( radius, thinImage, region );\r
 112   ot.SetBoundaryCondition( boundaryCondition );\r
 113 \r
 114   std::vector < IndexType > simpleBorderPoints;\r
 115   typename std::vector < IndexType >::iterator simpleBorderPointsIt;\r
 116 \r
 117   // Define offsets\r
 118   typedef typename NeighborhoodIteratorType::OffsetType OffsetType;\r
 119   OffsetType N   = {{ 0,-1, 0}};  // north\r
 120   OffsetType S   = {{ 0, 1, 0}};  // south\r
 121   OffsetType E   = {{ 1, 0, 0}};  // east\r
 122   OffsetType W   = {{-1, 0, 0}};  // west\r
 123   OffsetType U   = {{ 0, 0, 1}};  // up\r
 124   OffsetType B   = {{ 0, 0,-1}};  // bottom\r
 125 \r
 126   // prepare Euler LUT [Lee94]\r
 127   int eulerLUT[256]; \r
 128   fillEulerLUT( eulerLUT );\r
 129   // Loop through the image several times until there is no change.\r
 130   int unchangedBorders = 0;\r
 131   while( unchangedBorders < 6 )  // loop until no change for all the six border types\r
 132   {\r
 133     unchangedBorders = 0;\r
 134     for( int currentBorder = 1; currentBorder <= 6; currentBorder++)\r
 135     {\r
 136       // Loop through the image.\r
 137       for ( ot.GoToBegin(); !ot.IsAtEnd(); ++ot )\r
 138       { \r
 139         // check if point is foreground\r
 140         if ( ot.GetCenterPixel() != 1 )\r
 141         {\r
 142           continue;         // current point is already background \r
 143         }\r
 144         // check 6-neighbors if point is a border point of type currentBorder\r
 145         bool isBorderPoint = false;\r
 146         if( currentBorder == 1 && ot.GetPixel(N)<=0 )\r
 147           isBorderPoint = true;\r
 148         if( currentBorder == 2 && ot.GetPixel(S)<=0 )\r
 149           isBorderPoint = true;\r
 150         if( currentBorder == 3 && ot.GetPixel(E)<=0 )\r
 151           isBorderPoint = true;\r
 152         if( currentBorder == 4 && ot.GetPixel(W)<=0 )\r
 153           isBorderPoint = true;\r
 154         if( currentBorder == 5 && ot.GetPixel(U)<=0 )\r
 155           isBorderPoint = true;\r
 156         if( currentBorder == 6 && ot.GetPixel(B)<=0 )\r
 157           isBorderPoint = true;\r
 158         if( !isBorderPoint )\r
 159         {\r
 160           continue;         // current point is not deletable\r
 161         }        \r
 162         // check if point is the end of an arc\r
 163         int numberOfNeighbors = -1;   // -1 and not 0 because the center pixel will be counted as well  \r
 164         for( int i = 0; i < 27; i++ ) // i =  0..26\r
 165           if( ot.GetPixel(i)==1 )\r
 166             numberOfNeighbors++;\r
 167 \r
 168         if( numberOfNeighbors == 1 )\r
 169         {\r
 170           continue;         // current point is not deletable\r
 171         }\r
 172 \r
 173         // check if point is Euler invariant\r
 174         if( !isEulerInvariant( ot.GetNeighborhood(), eulerLUT ) )\r
 175         {\r
 176           continue;         // current point is not deletable\r
 177         }\r
 178 \r
 179         // check if point is simple (deletion does not change connectivity in the 3x3x3 neighborhood)\r
 180         if( !isSimplePoint( ot.GetNeighborhood() ) )\r
 181         {\r
 182           continue;         // current point is not deletable\r
 183         }\r
 184 \r
 185         // add all simple border points to a list for sequential re-checking\r
 186         simpleBorderPoints.push_back( ot.GetIndex() );\r
 187       } // end image iteration loop\r
 188 \r
 189       // sequential re-checking to preserve connectivity when\r
 190       // deleting in a parallel way\r
 191       bool noChange = true;\r
 192       for( simpleBorderPointsIt=simpleBorderPoints.begin(); simpleBorderPointsIt!=simpleBorderPoints.end(); simpleBorderPointsIt++)\r
 193       {\r
 194         // 1. Set simple border point to 0\r
 195         thinImage->SetPixel( *simpleBorderPointsIt, NumericTraits<OutputImagePixelType>::Zero);\r
 196         // 2. Check if neighborhood is still connected\r
 197         ot.SetLocation( *simpleBorderPointsIt );\r
 198         if( !isSimplePoint( ot.GetNeighborhood() ) )\r
 199         {\r
 200           // we cannot delete current point, so reset\r
 201           thinImage->SetPixel( *simpleBorderPointsIt, NumericTraits<OutputImagePixelType>::One );\r
 202         }\r
 203         else\r
 204         {\r
 205           noChange = false;\r
 206         }\r
 207       }\r
 208       if( noChange )\r
 209         unchangedBorders++;\r
 210 \r
 211       simpleBorderPoints.clear();\r
 212     } // end currentBorder for loop\r
 213   } // end unchangedBorders while loop\r
 214 \r
 215   itkDebugMacro( << "ComputeThinImage End");\r
 216 }\r
 217 \r
 218 /**\r
 219  *  Generate ThinImage\r
 220  */\r
 221 template <class TInputImage,class TOutputImage>\r
 222 void \r
 223 BinaryThinningImageFilter3D<TInputImage,TOutputImage>\r
 224 ::GenerateData() \r
 225 {\r
 226 \r
 227   this->PrepareData();\r
 228 \r
 229   itkDebugMacro(<< "GenerateData: Computing Thinning Image");\r
 230   this->ComputeThinImage();\r
 231 } // end GenerateData()\r
 232 \r
 233 /** \r
 234  * Fill the Euler look-up table (LUT) for later check of the Euler invariance. (see [Lee94])\r
 235  */\r
 236 template <class TInputImage,class TOutputImage>\r
 237 void \r
 238 BinaryThinningImageFilter3D<TInputImage,TOutputImage>\r
 239 ::fillEulerLUT(int *LUT)\r
 240 {\r
 241   LUT[1]  =  1;\r
 242   LUT[3]  = -1;\r
 243   LUT[5]  = -1;\r
 244   LUT[7]  =  1;\r
 245   LUT[9]  = -3;\r
 246   LUT[11] = -1;\r
 247   LUT[13] = -1;\r
 248   LUT[15] =  1;\r
 249   LUT[17] = -1;\r
 250   LUT[19] =  1;\r
 251   LUT[21] =  1;\r
 252   LUT[23] = -1;\r
 253   LUT[25] =  3;\r
 254   LUT[27] =  1;\r
 255   LUT[29] =  1;\r
 256   LUT[31] = -1;\r
 257   LUT[33] = -3;\r
 258   LUT[35] = -1;\r
 259   LUT[37] =  3;\r
 260   LUT[39] =  1;\r
 261   LUT[41] =  1;\r
 262   LUT[43] = -1;\r
 263   LUT[45] =  3;\r
 264   LUT[47] =  1;\r
 265   LUT[49] = -1;\r
 266   LUT[51] =  1;\r
 267 \r
 268   LUT[53] =  1;\r
 269   LUT[55] = -1;\r
 270   LUT[57] =  3;\r
 271   LUT[59] =  1;\r
 272   LUT[61] =  1;\r
 273   LUT[63] = -1;\r
 274   LUT[65] = -3;\r
 275   LUT[67] =  3;\r
 276   LUT[69] = -1;\r
 277   LUT[71] =  1;\r
 278   LUT[73] =  1;\r
 279   LUT[75] =  3;\r
 280   LUT[77] = -1;\r
 281   LUT[79] =  1;\r
 282   LUT[81] = -1;\r
 283   LUT[83] =  1;\r
 284   LUT[85] =  1;\r
 285   LUT[87] = -1;\r
 286   LUT[89] =  3;\r
 287   LUT[91] =  1;\r
 288   LUT[93] =  1;\r
 289   LUT[95] = -1;\r
 290   LUT[97] =  1;\r
 291   LUT[99] =  3;\r
 292   LUT[101] =  3;\r
 293   LUT[103] =  1;\r
 294 \r
 295   LUT[105] =  5;\r
 296   LUT[107] =  3;\r
 297   LUT[109] =  3;\r
 298   LUT[111] =  1;\r
 299   LUT[113] = -1;\r
 300   LUT[115] =  1;\r
 301   LUT[117] =  1;\r
 302   LUT[119] = -1;\r
 303   LUT[121] =  3;\r
 304   LUT[123] =  1;\r
 305   LUT[125] =  1;\r
 306   LUT[127] = -1;\r
 307   LUT[129] = -7;\r
 308   LUT[131] = -1;\r
 309   LUT[133] = -1;\r
 310   LUT[135] =  1;\r
 311   LUT[137] = -3;\r
 312   LUT[139] = -1;\r
 313   LUT[141] = -1;\r
 314   LUT[143] =  1;\r
 315   LUT[145] = -1;\r
 316   LUT[147] =  1;\r
 317   LUT[149] =  1;\r
 318   LUT[151] = -1;\r
 319   LUT[153] =  3;\r
 320   LUT[155] =  1;\r
 321 \r
 322   LUT[157] =  1;\r
 323   LUT[159] = -1;\r
 324   LUT[161] = -3;\r
 325   LUT[163] = -1;\r
 326   LUT[165] =  3;\r
 327   LUT[167] =  1;\r
 328   LUT[169] =  1;\r
 329   LUT[171] = -1;\r
 330   LUT[173] =  3;\r
 331   LUT[175] =  1;\r
 332   LUT[177] = -1;\r
 333   LUT[179] =  1;\r
 334   LUT[181] =  1;\r
 335   LUT[183] = -1;\r
 336   LUT[185] =  3;\r
 337   LUT[187] =  1;\r
 338   LUT[189] =  1;\r
 339   LUT[191] = -1;\r
 340   LUT[193] = -3;\r
 341   LUT[195] =  3;\r
 342   LUT[197] = -1;\r
 343   LUT[199] =  1;\r
 344   LUT[201] =  1;\r
 345   LUT[203] =  3;\r
 346   LUT[205] = -1;\r
 347   LUT[207] =  1;\r
 348 \r
 349   LUT[209] = -1;\r
 350   LUT[211] =  1;\r
 351   LUT[213] =  1;\r
 352   LUT[215] = -1;\r
 353   LUT[217] =  3;\r
 354   LUT[219] =  1;\r
 355   LUT[221] =  1;\r
 356   LUT[223] = -1;\r
 357   LUT[225] =  1;\r
 358   LUT[227] =  3;\r
 359   LUT[229] =  3;\r
 360   LUT[231] =  1;\r
 361   LUT[233] =  5;\r
 362   LUT[235] =  3;\r
 363   LUT[237] =  3;\r
 364   LUT[239] =  1;\r
 365   LUT[241] = -1;\r
 366   LUT[243] =  1;\r
 367   LUT[245] =  1;\r
 368   LUT[247] = -1;\r
 369   LUT[249] =  3;\r
 370   LUT[251] =  1;\r
 371   LUT[253] =  1;\r
 372   LUT[255] = -1;\r
 373 }\r
 374 \r
 375 /** \r
 376  * Check for Euler invariance. (see [Lee94])\r
 377  */\r
 378 template <class TInputImage,class TOutputImage>\r
 379 bool \r
 380 BinaryThinningImageFilter3D<TInputImage,TOutputImage>\r
 381 ::isEulerInvariant(NeighborhoodType neighbors, int *LUT)\r
 382 {\r
 383   // calculate Euler characteristic for each octant and sum up\r
 384   int EulerChar = 0;\r
 385   unsigned char n;\r
 386   // Octant SWU\r
 387   n = 1;\r
 388   if( neighbors[24]==1 )\r
 389     n |= 128;\r
 390   if( neighbors[25]==1 )\r
 391     n |=  64;\r
 392   if( neighbors[15]==1 )\r
 393     n |=  32;\r
 394   if( neighbors[16]==1 )\r
 395     n |=  16;\r
 396   if( neighbors[21]==1 )\r
 397     n |=   8;\r
 398   if( neighbors[22]==1 )\r
 399     n |=   4;\r
 400   if( neighbors[12]==1 )\r
 401     n |=   2;\r
 402   EulerChar += LUT[n];\r
 403   // Octant SEU\r
 404   n = 1;\r
 405   if( neighbors[26]==1 )\r
 406     n |= 128;\r
 407   if( neighbors[23]==1 )\r
 408     n |=  64;\r
 409   if( neighbors[17]==1 )\r
 410     n |=  32;\r
 411   if( neighbors[14]==1 )\r
 412     n |=  16;\r
 413   if( neighbors[25]==1 )\r
 414     n |=   8;\r
 415   if( neighbors[22]==1 )\r
 416     n |=   4;\r
 417   if( neighbors[16]==1 )\r
 418     n |=   2;\r
 419   EulerChar += LUT[n];\r
 420   // Octant NWU\r
 421   n = 1;\r
 422   if( neighbors[18]==1 )\r
 423     n |= 128;\r
 424   if( neighbors[21]==1 )\r
 425     n |=  64;\r
 426   if( neighbors[9]==1 )\r
 427     n |=  32;\r
 428   if( neighbors[12]==1 )\r
 429     n |=  16;\r
 430   if( neighbors[19]==1 )\r
 431     n |=   8;\r
 432   if( neighbors[22]==1 )\r
 433     n |=   4;\r
 434   if( neighbors[10]==1 )\r
 435     n |=   2;\r
 436   EulerChar += LUT[n];\r
 437   // Octant NEU\r
 438   n = 1;\r
 439   if( neighbors[20]==1 )\r
 440     n |= 128;\r
 441   if( neighbors[23]==1 )\r
 442     n |=  64;\r
 443   if( neighbors[19]==1 )\r
 444     n |=  32;\r
 445   if( neighbors[22]==1 )\r
 446     n |=  16;\r
 447   if( neighbors[11]==1 )\r
 448     n |=   8;\r
 449   if( neighbors[14]==1 )\r
 450     n |=   4;\r
 451   if( neighbors[10]==1 )\r
 452     n |=   2;\r
 453   EulerChar += LUT[n];\r
 454   // Octant SWB\r
 455   n = 1;\r
 456   if( neighbors[6]==1 )\r
 457     n |= 128;\r
 458   if( neighbors[15]==1 )\r
 459     n |=  64;\r
 460   if( neighbors[7]==1 )\r
 461     n |=  32;\r
 462   if( neighbors[16]==1 )\r
 463     n |=  16;\r
 464   if( neighbors[3]==1 )\r
 465     n |=   8;\r
 466   if( neighbors[12]==1 )\r
 467     n |=   4;\r
 468   if( neighbors[4]==1 )\r
 469     n |=   2;\r
 470   EulerChar += LUT[n];\r
 471   // Octant SEB\r
 472   n = 1;\r
 473   if( neighbors[8]==1 )\r
 474     n |= 128;\r
 475   if( neighbors[7]==1 )\r
 476     n |=  64;\r
 477   if( neighbors[17]==1 )\r
 478     n |=  32;\r
 479   if( neighbors[16]==1 )\r
 480     n |=  16;\r
 481   if( neighbors[5]==1 )\r
 482     n |=   8;\r
 483   if( neighbors[4]==1 )\r
 484     n |=   4;\r
 485   if( neighbors[14]==1 )\r
 486     n |=   2;\r
 487   EulerChar += LUT[n];\r
 488   // Octant NWB\r
 489   n = 1;\r
 490   if( neighbors[0]==1 )\r
 491     n |= 128;\r
 492   if( neighbors[9]==1 )\r
 493     n |=  64;\r
 494   if( neighbors[3]==1 )\r
 495     n |=  32;\r
 496   if( neighbors[12]==1 )\r
 497     n |=  16;\r
 498   if( neighbors[1]==1 )\r
 499     n |=   8;\r
 500   if( neighbors[10]==1 )\r
 501     n |=   4;\r
 502   if( neighbors[4]==1 )\r
 503     n |=   2;\r
 504   EulerChar += LUT[n];\r
 505   // Octant NEB\r
 506   n = 1;\r
 507   if( neighbors[2]==1 )\r
 508     n |= 128;\r
 509   if( neighbors[1]==1 )\r
 510     n |=  64;\r
 511   if( neighbors[11]==1 )\r
 512     n |=  32;\r
 513   if( neighbors[10]==1 )\r
 514     n |=  16;\r
 515   if( neighbors[5]==1 )\r
 516     n |=   8;\r
 517   if( neighbors[4]==1 )\r
 518     n |=   4;\r
 519   if( neighbors[14]==1 )\r
 520     n |=   2;\r
 521   EulerChar += LUT[n];\r
 522   if( EulerChar == 0 )\r
 523     return true;\r
 524   else\r
 525     return false;\r
 526 }\r
 527 \r
 528 /** \r
 529  * Check if current point is a Simple Point.\r
 530  * This method is named 'N(v)_labeling' in [Lee94].\r
 531  * Outputs the number of connected objects in a neighborhood of a point\r
 532  * after this point would have been removed.\r
 533  */\r
 534 template <class TInputImage,class TOutputImage>\r
 535 bool \r
 536 BinaryThinningImageFilter3D<TInputImage,TOutputImage>\r
 537 ::isSimplePoint(NeighborhoodType neighbors)\r
 538 {\r
 539   // copy neighbors for labeling\r
 540   int cube[26];\r
 541   int i;\r
 542   for( i = 0; i < 13; i++ )  // i =  0..12 -> cube[0..12]\r
 543     cube[i] = neighbors[i];\r
 544   // i != 13 : ignore center pixel when counting (see [Lee94])\r
 545   for( i = 14; i < 27; i++ ) // i = 14..26 -> cube[13..25]\r
 546     cube[i-1] = neighbors[i];\r
 547   // set initial label\r
 548   int label = 2;\r
 549   // for all points in the neighborhood\r
 550   for( int i = 0; i < 26; i++ )\r
 551   {\r
 552     if( cube[i]==1 )     // voxel has not been labelled yet\r
 553     {\r
 554       // start recursion with any octant that contains the point i\r
 555       switch( i )\r
 556       {\r
 557       case 0:\r
 558       case 1:\r
 559       case 3:\r
 560       case 4:\r
 561       case 9:\r
 562       case 10:\r
 563       case 12:\r
 564         Octree_labeling(1, label, cube );\r
 565         break;\r
 566       case 2:\r
 567       case 5:\r
 568       case 11:\r
 569       case 13:\r
 570         Octree_labeling(2, label, cube );\r
 571         break;\r
 572       case 6:\r
 573       case 7:\r
 574       case 14:\r
 575       case 15:\r
 576         Octree_labeling(3, label, cube );\r
 577         break;\r
 578       case 8:\r
 579       case 16:\r
 580         Octree_labeling(4, label, cube );\r
 581         break;\r
 582       case 17:\r
 583       case 18:\r
 584       case 20:\r
 585       case 21:\r
 586         Octree_labeling(5, label, cube );\r
 587         break;\r
 588       case 19:\r
 589       case 22:\r
 590         Octree_labeling(6, label, cube );\r
 591         break;\r
 592       case 23:\r
 593       case 24:\r
 594         Octree_labeling(7, label, cube );\r
 595         break;\r
 596       case 25:\r
 597         Octree_labeling(8, label, cube );\r
 598         break;\r
 599       }\r
 600       label++;\r
 601       if( label-2 >= 2 )\r
 602       {\r
 603         return false;\r
 604       }\r
 605     }\r
 606   }\r
 607   //return label-2; in [Lee94] if the number of connected compontents would be needed\r
 608   return true;\r
 609 }\r
 610 \r
 611 /** \r
 612  * Octree_labeling [Lee94]\r
 613  * This is a recursive method that calulates the number of connected\r
 614  * components in the 3D neighbourhood after the center pixel would\r
 615  * have been removed.\r
 616  */\r
 617 template <class TInputImage,class TOutputImage>\r
 618 void \r
 619 BinaryThinningImageFilter3D<TInputImage,TOutputImage>\r
 620 ::Octree_labeling(int octant, int label, int *cube)\r
 621 {\r
 622   // check if there are points in the octant with value 1\r
 623   if( octant==1 )\r
 624   {\r
 625         // set points in this octant to current label\r
 626         // and recurseive labeling of adjacent octants\r
 627     if( cube[0] == 1 )\r
 628       cube[0] = label;\r
 629     if( cube[1] == 1 )\r
 630     {\r
 631       cube[1] = label;        \r
 632       Octree_labeling( 2, label, cube);\r
 633     }\r
 634     if( cube[3] == 1 )\r
 635     {\r
 636       cube[3] = label;        \r
 637       Octree_labeling( 3, label, cube);\r
 638     }\r
 639     if( cube[4] == 1 )\r
 640     {\r
 641       cube[4] = label;        \r
 642       Octree_labeling( 2, label, cube);\r
 643       Octree_labeling( 3, label, cube);\r
 644       Octree_labeling( 4, label, cube);\r
 645     }\r
 646     if( cube[9] == 1 )\r
 647     {\r
 648       cube[9] = label;        \r
 649       Octree_labeling( 5, label, cube);\r
 650     }\r
 651     if( cube[10] == 1 )\r
 652     {\r
 653       cube[10] = label;        \r
 654       Octree_labeling( 2, label, cube);\r
 655       Octree_labeling( 5, label, cube);\r
 656       Octree_labeling( 6, label, cube);\r
 657     }\r
 658     if( cube[12] == 1 )\r
 659     {\r
 660       cube[12] = label;        \r
 661       Octree_labeling( 3, label, cube);\r
 662       Octree_labeling( 5, label, cube);\r
 663       Octree_labeling( 7, label, cube);\r
 664     }\r
 665   }\r
 666   if( octant==2 )\r
 667   {\r
 668     if( cube[1] == 1 )\r
 669     {\r
 670       cube[1] = label;\r
 671       Octree_labeling( 1, label, cube);\r
 672     }\r
 673     if( cube[4] == 1 )\r
 674     {\r
 675       cube[4] = label;        \r
 676       Octree_labeling( 1, label, cube);\r
 677       Octree_labeling( 3, label, cube);\r
 678       Octree_labeling( 4, label, cube);\r
 679     }\r
 680     if( cube[10] == 1 )\r
 681     {\r
 682       cube[10] = label;        \r
 683       Octree_labeling( 1, label, cube);\r
 684       Octree_labeling( 5, label, cube);\r
 685       Octree_labeling( 6, label, cube);\r
 686     }\r
 687     if( cube[2] == 1 )\r
 688       cube[2] = label;        \r
 689     if( cube[5] == 1 )\r
 690     {\r
 691       cube[5] = label;        \r
 692       Octree_labeling( 4, label, cube);\r
 693     }\r
 694     if( cube[11] == 1 )\r
 695     {\r
 696       cube[11] = label;        \r
 697       Octree_labeling( 6, label, cube);\r
 698     }\r
 699     if( cube[13] == 1 )\r
 700     {\r
 701       cube[13] = label;        \r
 702       Octree_labeling( 4, label, cube);\r
 703       Octree_labeling( 6, label, cube);\r
 704       Octree_labeling( 8, label, cube);\r
 705     }\r
 706   }\r
 707   if( octant==3 )\r
 708   {\r
 709     if( cube[3] == 1 )\r
 710     {\r
 711       cube[3] = label;        \r
 712       Octree_labeling( 1, label, cube);\r
 713     }\r
 714     if( cube[4] == 1 )\r
 715     {\r
 716       cube[4] = label;        \r
 717       Octree_labeling( 1, label, cube);\r
 718       Octree_labeling( 2, label, cube);\r
 719       Octree_labeling( 4, label, cube);\r
 720     }\r
 721     if( cube[12] == 1 )\r
 722     {\r
 723       cube[12] = label;        \r
 724       Octree_labeling( 1, label, cube);\r
 725       Octree_labeling( 5, label, cube);\r
 726       Octree_labeling( 7, label, cube);\r
 727     }\r
 728     if( cube[6] == 1 )\r
 729       cube[6] = label;        \r
 730     if( cube[7] == 1 )\r
 731     {\r
 732       cube[7] = label;        \r
 733       Octree_labeling( 4, label, cube);\r
 734     }\r
 735     if( cube[14] == 1 )\r
 736     {\r
 737       cube[14] = label;        \r
 738       Octree_labeling( 7, label, cube);\r
 739     }\r
 740     if( cube[15] == 1 )\r
 741     {\r
 742       cube[15] = label;        \r
 743       Octree_labeling( 4, label, cube);\r
 744       Octree_labeling( 7, label, cube);\r
 745       Octree_labeling( 8, label, cube);\r
 746     }\r
 747   }\r
 748   if( octant==4 )\r
 749   {\r
 750         if( cube[4] == 1 )\r
 751     {\r
 752       cube[4] = label;        \r
 753       Octree_labeling( 1, label, cube);\r
 754       Octree_labeling( 2, label, cube);\r
 755       Octree_labeling( 3, label, cube);\r
 756     }\r
 757         if( cube[5] == 1 )\r
 758     {\r
 759       cube[5] = label;        \r
 760       Octree_labeling( 2, label, cube);\r
 761     }\r
 762     if( cube[13] == 1 )\r
 763     {\r
 764       cube[13] = label;        \r
 765       Octree_labeling( 2, label, cube);\r
 766       Octree_labeling( 6, label, cube);\r
 767       Octree_labeling( 8, label, cube);\r
 768     }\r
 769     if( cube[7] == 1 )\r
 770     {\r
 771       cube[7] = label;        \r
 772       Octree_labeling( 3, label, cube);\r
 773     }\r
 774     if( cube[15] == 1 )\r
 775     {\r
 776       cube[15] = label;        \r
 777       Octree_labeling( 3, label, cube);\r
 778       Octree_labeling( 7, label, cube);\r
 779       Octree_labeling( 8, label, cube);\r
 780     }\r
 781     if( cube[8] == 1 )\r
 782       cube[8] = label;        \r
 783     if( cube[16] == 1 )\r
 784     {\r
 785       cube[16] = label;        \r
 786       Octree_labeling( 8, label, cube);\r
 787     }\r
 788   }\r
 789   if( octant==5 )\r
 790   {\r
 791         if( cube[9] == 1 )\r
 792     {\r
 793       cube[9] = label;        \r
 794       Octree_labeling( 1, label, cube);\r
 795     }\r
 796     if( cube[10] == 1 )\r
 797     {\r
 798       cube[10] = label;        \r
 799       Octree_labeling( 1, label, cube);\r
 800       Octree_labeling( 2, label, cube);\r
 801       Octree_labeling( 6, label, cube);\r
 802     }\r
 803     if( cube[12] == 1 )\r
 804     {\r
 805       cube[12] = label;        \r
 806       Octree_labeling( 1, label, cube);\r
 807       Octree_labeling( 3, label, cube);\r
 808       Octree_labeling( 7, label, cube);\r
 809     }\r
 810     if( cube[17] == 1 )\r
 811       cube[17] = label;        \r
 812     if( cube[18] == 1 )\r
 813     {\r
 814       cube[18] = label;        \r
 815       Octree_labeling( 6, label, cube);\r
 816     }\r
 817     if( cube[20] == 1 )\r
 818     {\r
 819       cube[20] = label;        \r
 820       Octree_labeling( 7, label, cube);\r
 821     }\r
 822     if( cube[21] == 1 )\r
 823     {\r
 824       cube[21] = label;        \r
 825       Octree_labeling( 6, label, cube);\r
 826       Octree_labeling( 7, label, cube);\r
 827       Octree_labeling( 8, label, cube);\r
 828     }\r
 829   }\r
 830   if( octant==6 )\r
 831   {\r
 832         if( cube[10] == 1 )\r
 833     {\r
 834       cube[10] = label;        \r
 835       Octree_labeling( 1, label, cube);\r
 836       Octree_labeling( 2, label, cube);\r
 837       Octree_labeling( 5, label, cube);\r
 838     }\r
 839     if( cube[11] == 1 )\r
 840     {\r
 841       cube[11] = label;        \r
 842       Octree_labeling( 2, label, cube);\r
 843     }\r
 844     if( cube[13] == 1 )\r
 845     {\r
 846       cube[13] = label;        \r
 847       Octree_labeling( 2, label, cube);\r
 848       Octree_labeling( 4, label, cube);\r
 849       Octree_labeling( 8, label, cube);\r
 850     }\r
 851     if( cube[18] == 1 )\r
 852     {\r
 853       cube[18] = label;        \r
 854       Octree_labeling( 5, label, cube);\r
 855     }\r
 856     if( cube[21] == 1 )\r
 857     {\r
 858       cube[21] = label;        \r
 859       Octree_labeling( 5, label, cube);\r
 860       Octree_labeling( 7, label, cube);\r
 861       Octree_labeling( 8, label, cube);\r
 862     }\r
 863     if( cube[19] == 1 )\r
 864       cube[19] = label;        \r
 865     if( cube[22] == 1 )\r
 866     {\r
 867       cube[22] = label;        \r
 868       Octree_labeling( 8, label, cube);\r
 869     }\r
 870   }\r
 871   if( octant==7 )\r
 872   {\r
 873         if( cube[12] == 1 )\r
 874     {\r
 875       cube[12] = label;        \r
 876       Octree_labeling( 1, label, cube);\r
 877       Octree_labeling( 3, label, cube);\r
 878       Octree_labeling( 5, label, cube);\r
 879     }\r
 880         if( cube[14] == 1 )\r
 881     {\r
 882       cube[14] = label;        \r
 883       Octree_labeling( 3, label, cube);\r
 884     }\r
 885     if( cube[15] == 1 )\r
 886     {\r
 887       cube[15] = label;        \r
 888       Octree_labeling( 3, label, cube);\r
 889       Octree_labeling( 4, label, cube);\r
 890       Octree_labeling( 8, label, cube);\r
 891     }\r
 892     if( cube[20] == 1 )\r
 893     {\r
 894       cube[20] = label;        \r
 895       Octree_labeling( 5, label, cube);\r
 896     }\r
 897     if( cube[21] == 1 )\r
 898     {\r
 899       cube[21] = label;        \r
 900       Octree_labeling( 5, label, cube);\r
 901       Octree_labeling( 6, label, cube);\r
 902       Octree_labeling( 8, label, cube);\r
 903     }\r
 904     if( cube[23] == 1 )\r
 905       cube[23] = label;        \r
 906     if( cube[24] == 1 )\r
 907     {\r
 908       cube[24] = label;        \r
 909       Octree_labeling( 8, label, cube);\r
 910     }\r
 911   }\r
 912   if( octant==8 )\r
 913   {\r
 914         if( cube[13] == 1 )\r
 915     {\r
 916       cube[13] = label;        \r
 917       Octree_labeling( 2, label, cube);\r
 918       Octree_labeling( 4, label, cube);\r
 919       Octree_labeling( 6, label, cube);\r
 920     }\r
 921         if( cube[15] == 1 )\r
 922     {\r
 923       cube[15] = label;        \r
 924       Octree_labeling( 3, label, cube);\r
 925       Octree_labeling( 4, label, cube);\r
 926       Octree_labeling( 7, label, cube);\r
 927     }\r
 928         if( cube[16] == 1 )\r
 929     {\r
 930       cube[16] = label;        \r
 931       Octree_labeling( 4, label, cube);\r
 932     }\r
 933         if( cube[21] == 1 )\r
 934     {\r
 935       cube[21] = label;        \r
 936       Octree_labeling( 5, label, cube);\r
 937       Octree_labeling( 6, label, cube);\r
 938       Octree_labeling( 7, label, cube);\r
 939     }\r
 940         if( cube[22] == 1 )\r
 941     {\r
 942       cube[22] = label;        \r
 943       Octree_labeling( 6, label, cube);\r
 944     }\r
 945         if( cube[24] == 1 )\r
 946     {\r
 947       cube[24] = label;        \r
 948       Octree_labeling( 7, label, cube);\r
 949     }\r
 950         if( cube[25] == 1 )\r
 951       cube[25] = label;        \r
 952   } \r
 953 }\r
 954 \r
 955 \r
 956 /**\r
 957  *  Print Self\r
 958  */\r
 959 template <class TInputImage,class TOutputImage>\r
 960 void \r
 961 BinaryThinningImageFilter3D<TInputImage,TOutputImage>\r
 962 ::PrintSelf(std::ostream& os, Indent indent) const\r
 963 {\r
 964   Superclass::PrintSelf(os,indent);\r
 965   \r
 966   os << indent << "Thinning image: " << std::endl;\r
 967 \r
 968 }\r
 969 \r
 970 } // end namespace itk\r
 971 \r
 972 #endif\r