1 #ifndef _itkBinaryThinningImageFilter3D_txx
2 #define _itkBinaryThinningImageFilter3D_txx
6 #include "itkBinaryThinningImageFilter3D.h"
7 #include "itkImageRegionConstIterator.h"
8 #include "itkImageRegionIterator.h"
9 #include "itkNeighborhoodIterator.h"
18 template <class TInputImage,class TOutputImage>
19 BinaryThinningImageFilter3D<TInputImage,TOutputImage>
20 ::BinaryThinningImageFilter3D()
23 this->SetNumberOfRequiredOutputs( 1 );
25 OutputImagePointer thinImage = OutputImageType::New();
26 this->SetNthOutput( 0, thinImage.GetPointer() );
31 * Return the thinning Image pointer
33 template <class TInputImage,class TOutputImage>
34 typename BinaryThinningImageFilter3D<
35 TInputImage,TOutputImage>::OutputImageType *
36 BinaryThinningImageFilter3D<TInputImage,TOutputImage>
39 return dynamic_cast< OutputImageType * >(
40 this->ProcessObject::GetOutput(0) );
45 * Prepare data for computation
46 * Copy the input image to the output image, changing from the input
47 * type to the output type.
49 template <class TInputImage,class TOutputImage>
51 BinaryThinningImageFilter3D<TInputImage,TOutputImage>
55 itkDebugMacro(<< "PrepareData Start");
56 OutputImagePointer thinImage = GetThinning();
58 InputImagePointer inputImage =
59 dynamic_cast<const TInputImage *>( ProcessObject::GetInput(0) );
61 thinImage->SetBufferedRegion( thinImage->GetRequestedRegion() );
62 thinImage->Allocate();
64 typename OutputImageType::RegionType region = thinImage->GetRequestedRegion();
67 ImageRegionConstIterator< TInputImage > it( inputImage, region );
68 ImageRegionIterator< TOutputImage > ot( thinImage, region );
73 itkDebugMacro(<< "PrepareData: Copy input to output");
75 // Copy the input to the output, changing all foreground pixels to
76 // have value 1 in the process.
77 while( !ot.IsAtEnd() ) {
79 ot.Set( NumericTraits<OutputImagePixelType>::One );
81 ot.Set( NumericTraits<OutputImagePixelType>::Zero );
86 itkDebugMacro(<< "PrepareData End");
90 * Post processing for computing thinning
92 template <class TInputImage,class TOutputImage>
94 BinaryThinningImageFilter3D<TInputImage,TOutputImage>
97 itkDebugMacro( << "ComputeThinImage Start");
98 OutputImagePointer thinImage = GetThinning();
100 typename OutputImageType::RegionType region = thinImage->GetRequestedRegion();
102 ConstBoundaryConditionType boundaryCondition;
103 boundaryCondition.SetConstant( 0 );
105 typename NeighborhoodIteratorType::RadiusType radius;
107 NeighborhoodIteratorType ot( radius, thinImage, region );
108 ot.SetBoundaryCondition( boundaryCondition );
110 std::vector < IndexType > simpleBorderPoints;
111 typename std::vector < IndexType >::iterator simpleBorderPointsIt;
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
122 // prepare Euler LUT [Lee94]
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
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
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 )
159 if( numberOfNeighbors == 1 ) {
160 continue; // current point is not deletable
163 // check if point is Euler invariant
164 if( !isEulerInvariant( ot.GetNeighborhood(), eulerLUT ) ) {
165 continue; // current point is not deletable
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
173 // add all simple border points to a list for sequential re-checking
174 simpleBorderPoints.push_back( ot.GetIndex() );
175 } // end image iteration loop
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 );
195 simpleBorderPoints.clear();
196 } // end currentBorder for loop
197 } // end unchangedBorders while loop
199 itkDebugMacro( << "ComputeThinImage End");
205 template <class TInputImage,class TOutputImage>
207 BinaryThinningImageFilter3D<TInputImage,TOutputImage>
213 itkDebugMacro(<< "GenerateData: Computing Thinning Image");
214 this->ComputeThinImage();
215 } // end GenerateData()
218 * Fill the Euler look-up table (LUT) for later check of the Euler invariance. (see [Lee94])
220 template <class TInputImage,class TOutputImage>
222 BinaryThinningImageFilter3D<TInputImage,TOutputImage>
223 ::fillEulerLUT(int *LUT)
360 * Check for Euler invariance. (see [Lee94])
362 template <class TInputImage,class TOutputImage>
364 BinaryThinningImageFilter3D<TInputImage,TOutputImage>
365 ::isEulerInvariant(NeighborhoodType neighbors, int *LUT)
367 // calculate Euler characteristic for each octant and sum up
372 if( neighbors[24]==1 )
374 if( neighbors[25]==1 )
376 if( neighbors[15]==1 )
378 if( neighbors[16]==1 )
380 if( neighbors[21]==1 )
382 if( neighbors[22]==1 )
384 if( neighbors[12]==1 )
389 if( neighbors[26]==1 )
391 if( neighbors[23]==1 )
393 if( neighbors[17]==1 )
395 if( neighbors[14]==1 )
397 if( neighbors[25]==1 )
399 if( neighbors[22]==1 )
401 if( neighbors[16]==1 )
406 if( neighbors[18]==1 )
408 if( neighbors[21]==1 )
410 if( neighbors[9]==1 )
412 if( neighbors[12]==1 )
414 if( neighbors[19]==1 )
416 if( neighbors[22]==1 )
418 if( neighbors[10]==1 )
423 if( neighbors[20]==1 )
425 if( neighbors[23]==1 )
427 if( neighbors[19]==1 )
429 if( neighbors[22]==1 )
431 if( neighbors[11]==1 )
433 if( neighbors[14]==1 )
435 if( neighbors[10]==1 )
440 if( neighbors[6]==1 )
442 if( neighbors[15]==1 )
444 if( neighbors[7]==1 )
446 if( neighbors[16]==1 )
448 if( neighbors[3]==1 )
450 if( neighbors[12]==1 )
452 if( neighbors[4]==1 )
457 if( neighbors[8]==1 )
459 if( neighbors[7]==1 )
461 if( neighbors[17]==1 )
463 if( neighbors[16]==1 )
465 if( neighbors[5]==1 )
467 if( neighbors[4]==1 )
469 if( neighbors[14]==1 )
474 if( neighbors[0]==1 )
476 if( neighbors[9]==1 )
478 if( neighbors[3]==1 )
480 if( neighbors[12]==1 )
482 if( neighbors[1]==1 )
484 if( neighbors[10]==1 )
486 if( neighbors[4]==1 )
491 if( neighbors[2]==1 )
493 if( neighbors[1]==1 )
495 if( neighbors[11]==1 )
497 if( neighbors[10]==1 )
499 if( neighbors[5]==1 )
501 if( neighbors[4]==1 )
503 if( neighbors[14]==1 )
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.
518 template <class TInputImage,class TOutputImage>
520 BinaryThinningImageFilter3D<TInputImage,TOutputImage>
521 ::isSimplePoint(NeighborhoodType neighbors)
523 // copy neighbors for labeling
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];
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
545 Octree_labeling(1, label, cube );
551 Octree_labeling(2, label, cube );
557 Octree_labeling(3, label, cube );
561 Octree_labeling(4, label, cube );
567 Octree_labeling(5, label, cube );
571 Octree_labeling(6, label, cube );
575 Octree_labeling(7, label, cube );
578 Octree_labeling(8, label, cube );
587 //return label-2; in [Lee94] if the number of connected compontents would be needed
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
597 template <class TInputImage,class TOutputImage>
599 BinaryThinningImageFilter3D<TInputImage,TOutputImage>
600 ::Octree_labeling(int octant, int label, int *cube)
602 // check if there are points in the octant with value 1
604 // set points in this octant to current label
605 // and recurseive labeling of adjacent octants
610 Octree_labeling( 2, label, cube);
614 Octree_labeling( 3, label, cube);
618 Octree_labeling( 2, label, cube);
619 Octree_labeling( 3, label, cube);
620 Octree_labeling( 4, label, cube);
624 Octree_labeling( 5, label, cube);
626 if( cube[10] == 1 ) {
628 Octree_labeling( 2, label, cube);
629 Octree_labeling( 5, label, cube);
630 Octree_labeling( 6, label, cube);
632 if( cube[12] == 1 ) {
634 Octree_labeling( 3, label, cube);
635 Octree_labeling( 5, label, cube);
636 Octree_labeling( 7, label, cube);
642 Octree_labeling( 1, label, cube);
646 Octree_labeling( 1, label, cube);
647 Octree_labeling( 3, label, cube);
648 Octree_labeling( 4, label, cube);
650 if( cube[10] == 1 ) {
652 Octree_labeling( 1, label, cube);
653 Octree_labeling( 5, label, cube);
654 Octree_labeling( 6, label, cube);
660 Octree_labeling( 4, label, cube);
662 if( cube[11] == 1 ) {
664 Octree_labeling( 6, label, cube);
666 if( cube[13] == 1 ) {
668 Octree_labeling( 4, label, cube);
669 Octree_labeling( 6, label, cube);
670 Octree_labeling( 8, label, cube);
676 Octree_labeling( 1, label, cube);
680 Octree_labeling( 1, label, cube);
681 Octree_labeling( 2, label, cube);
682 Octree_labeling( 4, label, cube);
684 if( cube[12] == 1 ) {
686 Octree_labeling( 1, label, cube);
687 Octree_labeling( 5, label, cube);
688 Octree_labeling( 7, label, cube);
694 Octree_labeling( 4, label, cube);
696 if( cube[14] == 1 ) {
698 Octree_labeling( 7, label, cube);
700 if( cube[15] == 1 ) {
702 Octree_labeling( 4, label, cube);
703 Octree_labeling( 7, label, cube);
704 Octree_labeling( 8, label, cube);
710 Octree_labeling( 1, label, cube);
711 Octree_labeling( 2, label, cube);
712 Octree_labeling( 3, label, cube);
716 Octree_labeling( 2, label, cube);
718 if( cube[13] == 1 ) {
720 Octree_labeling( 2, label, cube);
721 Octree_labeling( 6, label, cube);
722 Octree_labeling( 8, label, cube);
726 Octree_labeling( 3, label, cube);
728 if( cube[15] == 1 ) {
730 Octree_labeling( 3, label, cube);
731 Octree_labeling( 7, label, cube);
732 Octree_labeling( 8, label, cube);
736 if( cube[16] == 1 ) {
738 Octree_labeling( 8, label, cube);
744 Octree_labeling( 1, label, cube);
746 if( cube[10] == 1 ) {
748 Octree_labeling( 1, label, cube);
749 Octree_labeling( 2, label, cube);
750 Octree_labeling( 6, label, cube);
752 if( cube[12] == 1 ) {
754 Octree_labeling( 1, label, cube);
755 Octree_labeling( 3, label, cube);
756 Octree_labeling( 7, label, cube);
760 if( cube[18] == 1 ) {
762 Octree_labeling( 6, label, cube);
764 if( cube[20] == 1 ) {
766 Octree_labeling( 7, label, cube);
768 if( cube[21] == 1 ) {
770 Octree_labeling( 6, label, cube);
771 Octree_labeling( 7, label, cube);
772 Octree_labeling( 8, label, cube);
776 if( cube[10] == 1 ) {
778 Octree_labeling( 1, label, cube);
779 Octree_labeling( 2, label, cube);
780 Octree_labeling( 5, label, cube);
782 if( cube[11] == 1 ) {
784 Octree_labeling( 2, label, cube);
786 if( cube[13] == 1 ) {
788 Octree_labeling( 2, label, cube);
789 Octree_labeling( 4, label, cube);
790 Octree_labeling( 8, label, cube);
792 if( cube[18] == 1 ) {
794 Octree_labeling( 5, label, cube);
796 if( cube[21] == 1 ) {
798 Octree_labeling( 5, label, cube);
799 Octree_labeling( 7, label, cube);
800 Octree_labeling( 8, label, cube);
804 if( cube[22] == 1 ) {
806 Octree_labeling( 8, label, cube);
810 if( cube[12] == 1 ) {
812 Octree_labeling( 1, label, cube);
813 Octree_labeling( 3, label, cube);
814 Octree_labeling( 5, label, cube);
816 if( cube[14] == 1 ) {
818 Octree_labeling( 3, label, cube);
820 if( cube[15] == 1 ) {
822 Octree_labeling( 3, label, cube);
823 Octree_labeling( 4, label, cube);
824 Octree_labeling( 8, label, cube);
826 if( cube[20] == 1 ) {
828 Octree_labeling( 5, label, cube);
830 if( cube[21] == 1 ) {
832 Octree_labeling( 5, label, cube);
833 Octree_labeling( 6, label, cube);
834 Octree_labeling( 8, label, cube);
838 if( cube[24] == 1 ) {
840 Octree_labeling( 8, label, cube);
844 if( cube[13] == 1 ) {
846 Octree_labeling( 2, label, cube);
847 Octree_labeling( 4, label, cube);
848 Octree_labeling( 6, label, cube);
850 if( cube[15] == 1 ) {
852 Octree_labeling( 3, label, cube);
853 Octree_labeling( 4, label, cube);
854 Octree_labeling( 7, label, cube);
856 if( cube[16] == 1 ) {
858 Octree_labeling( 4, label, cube);
860 if( cube[21] == 1 ) {
862 Octree_labeling( 5, label, cube);
863 Octree_labeling( 6, label, cube);
864 Octree_labeling( 7, label, cube);
866 if( cube[22] == 1 ) {
868 Octree_labeling( 6, label, cube);
870 if( cube[24] == 1 ) {
872 Octree_labeling( 7, label, cube);
883 template <class TInputImage,class TOutputImage>
885 BinaryThinningImageFilter3D<TInputImage,TOutputImage>
886 ::PrintSelf(std::ostream& os, Indent indent) const
888 Superclass::PrintSelf(os,indent);
890 os << indent << "Thinning image: " << std::endl;
894 } // end namespace itk