itk/itkBinaryThinningImageFilter3D.txx

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