itk/itkBinaryThinningImageFilter3D.txx

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