1 /*=========================================================================
2 Program: vv http://www.creatis.insa-lyon.fr/rio/vv
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
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.
13 It is distributed under dual licence
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
23 #include "itkBinaryThinningImageFilter3D.h"
24 #include "itkImageRegionConstIterator.h"
25 #include "itkImageRegionIterator.h"
26 #include "itkNeighborhoodIterator.h"
35 template <class TInputImage,class TOutputImage>
36 BinaryThinningImageFilter3D<TInputImage,TOutputImage>
37 ::BinaryThinningImageFilter3D()
40 this->SetNumberOfRequiredOutputs( 1 );
42 OutputImagePointer thinImage = OutputImageType::New();
43 this->SetNthOutput( 0, thinImage.GetPointer() );
48 * Return the thinning Image pointer
50 template <class TInputImage,class TOutputImage>
51 typename BinaryThinningImageFilter3D<
52 TInputImage,TOutputImage>::OutputImageType *
53 BinaryThinningImageFilter3D<TInputImage,TOutputImage>
56 return dynamic_cast< OutputImageType * >(
57 this->ProcessObject::GetOutput(0) );
62 * Prepare data for computation
63 * Copy the input image to the output image, changing from the input
64 * type to the output type.
66 template <class TInputImage,class TOutputImage>
68 BinaryThinningImageFilter3D<TInputImage,TOutputImage>
72 itkDebugMacro(<< "PrepareData Start");
73 OutputImagePointer thinImage = GetThinning();
75 InputImagePointer inputImage =
76 dynamic_cast<const TInputImage *>( ProcessObject::GetInput(0) );
78 thinImage->SetBufferedRegion( thinImage->GetRequestedRegion() );
79 thinImage->Allocate();
81 typename OutputImageType::RegionType region = thinImage->GetRequestedRegion();
84 ImageRegionConstIterator< TInputImage > it( inputImage, region );
85 ImageRegionIterator< TOutputImage > ot( thinImage, region );
90 itkDebugMacro(<< "PrepareData: Copy input to output");
92 // Copy the input to the output, changing all foreground pixels to
93 // have value 1 in the process.
94 while( !ot.IsAtEnd() ) {
96 ot.Set( NumericTraits<OutputImagePixelType>::One );
98 ot.Set( NumericTraits<OutputImagePixelType>::Zero );
103 itkDebugMacro(<< "PrepareData End");
107 * Post processing for computing thinning
109 template <class TInputImage,class TOutputImage>
111 BinaryThinningImageFilter3D<TInputImage,TOutputImage>
114 itkDebugMacro( << "ComputeThinImage Start");
115 OutputImagePointer thinImage = GetThinning();
117 typename OutputImageType::RegionType region = thinImage->GetRequestedRegion();
119 ConstBoundaryConditionType boundaryCondition;
120 boundaryCondition.SetConstant( 0 );
122 typename NeighborhoodIteratorType::RadiusType radius;
124 NeighborhoodIteratorType ot( radius, thinImage, region );
125 ot.SetBoundaryCondition( boundaryCondition );
127 std::vector < IndexType > simpleBorderPoints;
128 typename std::vector < IndexType >::iterator simpleBorderPointsIt;
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
139 // prepare Euler LUT [Lee94]
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
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
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 )
176 if( numberOfNeighbors == 1 ) {
177 continue; // current point is not deletable
180 // check if point is Euler invariant
181 if( !isEulerInvariant( ot.GetNeighborhood(), eulerLUT ) ) {
182 continue; // current point is not deletable
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
190 // add all simple border points to a list for sequential re-checking
191 simpleBorderPoints.push_back( ot.GetIndex() );
192 } // end image iteration loop
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 );
212 simpleBorderPoints.clear();
213 } // end currentBorder for loop
214 } // end unchangedBorders while loop
216 itkDebugMacro( << "ComputeThinImage End");
222 template <class TInputImage,class TOutputImage>
224 BinaryThinningImageFilter3D<TInputImage,TOutputImage>
230 itkDebugMacro(<< "GenerateData: Computing Thinning Image");
231 this->ComputeThinImage();
232 } // end GenerateData()
235 * Fill the Euler look-up table (LUT) for later check of the Euler invariance. (see [Lee94])
237 template <class TInputImage,class TOutputImage>
239 BinaryThinningImageFilter3D<TInputImage,TOutputImage>
240 ::fillEulerLUT(int *LUT)
377 * Check for Euler invariance. (see [Lee94])
379 template <class TInputImage,class TOutputImage>
381 BinaryThinningImageFilter3D<TInputImage,TOutputImage>
382 ::isEulerInvariant(NeighborhoodType neighbors, int *LUT)
384 // calculate Euler characteristic for each octant and sum up
389 if( neighbors[24]==1 )
391 if( neighbors[25]==1 )
393 if( neighbors[15]==1 )
395 if( neighbors[16]==1 )
397 if( neighbors[21]==1 )
399 if( neighbors[22]==1 )
401 if( neighbors[12]==1 )
406 if( neighbors[26]==1 )
408 if( neighbors[23]==1 )
410 if( neighbors[17]==1 )
412 if( neighbors[14]==1 )
414 if( neighbors[25]==1 )
416 if( neighbors[22]==1 )
418 if( neighbors[16]==1 )
423 if( neighbors[18]==1 )
425 if( neighbors[21]==1 )
427 if( neighbors[9]==1 )
429 if( neighbors[12]==1 )
431 if( neighbors[19]==1 )
433 if( neighbors[22]==1 )
435 if( neighbors[10]==1 )
440 if( neighbors[20]==1 )
442 if( neighbors[23]==1 )
444 if( neighbors[19]==1 )
446 if( neighbors[22]==1 )
448 if( neighbors[11]==1 )
450 if( neighbors[14]==1 )
452 if( neighbors[10]==1 )
457 if( neighbors[6]==1 )
459 if( neighbors[15]==1 )
461 if( neighbors[7]==1 )
463 if( neighbors[16]==1 )
465 if( neighbors[3]==1 )
467 if( neighbors[12]==1 )
469 if( neighbors[4]==1 )
474 if( neighbors[8]==1 )
476 if( neighbors[7]==1 )
478 if( neighbors[17]==1 )
480 if( neighbors[16]==1 )
482 if( neighbors[5]==1 )
484 if( neighbors[4]==1 )
486 if( neighbors[14]==1 )
491 if( neighbors[0]==1 )
493 if( neighbors[9]==1 )
495 if( neighbors[3]==1 )
497 if( neighbors[12]==1 )
499 if( neighbors[1]==1 )
501 if( neighbors[10]==1 )
503 if( neighbors[4]==1 )
508 if( neighbors[2]==1 )
510 if( neighbors[1]==1 )
512 if( neighbors[11]==1 )
514 if( neighbors[10]==1 )
516 if( neighbors[5]==1 )
518 if( neighbors[4]==1 )
520 if( neighbors[14]==1 )
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.
535 template <class TInputImage,class TOutputImage>
537 BinaryThinningImageFilter3D<TInputImage,TOutputImage>
538 ::isSimplePoint(NeighborhoodType neighbors)
540 // copy neighbors for labeling
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];
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
562 Octree_labeling(1, label, cube );
568 Octree_labeling(2, label, cube );
574 Octree_labeling(3, label, cube );
578 Octree_labeling(4, label, cube );
584 Octree_labeling(5, label, cube );
588 Octree_labeling(6, label, cube );
592 Octree_labeling(7, label, cube );
595 Octree_labeling(8, label, cube );
604 //return label-2; in [Lee94] if the number of connected compontents would be needed
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
614 template <class TInputImage,class TOutputImage>
616 BinaryThinningImageFilter3D<TInputImage,TOutputImage>
617 ::Octree_labeling(int octant, int label, int *cube)
619 // check if there are points in the octant with value 1
621 // set points in this octant to current label
622 // and recurseive labeling of adjacent octants
627 Octree_labeling( 2, label, cube);
631 Octree_labeling( 3, label, cube);
635 Octree_labeling( 2, label, cube);
636 Octree_labeling( 3, label, cube);
637 Octree_labeling( 4, label, cube);
641 Octree_labeling( 5, label, cube);
643 if( cube[10] == 1 ) {
645 Octree_labeling( 2, label, cube);
646 Octree_labeling( 5, label, cube);
647 Octree_labeling( 6, label, cube);
649 if( cube[12] == 1 ) {
651 Octree_labeling( 3, label, cube);
652 Octree_labeling( 5, label, cube);
653 Octree_labeling( 7, label, cube);
659 Octree_labeling( 1, label, cube);
663 Octree_labeling( 1, label, cube);
664 Octree_labeling( 3, label, cube);
665 Octree_labeling( 4, label, cube);
667 if( cube[10] == 1 ) {
669 Octree_labeling( 1, label, cube);
670 Octree_labeling( 5, label, cube);
671 Octree_labeling( 6, label, cube);
677 Octree_labeling( 4, label, cube);
679 if( cube[11] == 1 ) {
681 Octree_labeling( 6, label, cube);
683 if( cube[13] == 1 ) {
685 Octree_labeling( 4, label, cube);
686 Octree_labeling( 6, label, cube);
687 Octree_labeling( 8, label, cube);
693 Octree_labeling( 1, label, cube);
697 Octree_labeling( 1, label, cube);
698 Octree_labeling( 2, label, cube);
699 Octree_labeling( 4, label, cube);
701 if( cube[12] == 1 ) {
703 Octree_labeling( 1, label, cube);
704 Octree_labeling( 5, label, cube);
705 Octree_labeling( 7, label, cube);
711 Octree_labeling( 4, label, cube);
713 if( cube[14] == 1 ) {
715 Octree_labeling( 7, label, cube);
717 if( cube[15] == 1 ) {
719 Octree_labeling( 4, label, cube);
720 Octree_labeling( 7, label, cube);
721 Octree_labeling( 8, label, cube);
727 Octree_labeling( 1, label, cube);
728 Octree_labeling( 2, label, cube);
729 Octree_labeling( 3, label, cube);
733 Octree_labeling( 2, label, cube);
735 if( cube[13] == 1 ) {
737 Octree_labeling( 2, label, cube);
738 Octree_labeling( 6, label, cube);
739 Octree_labeling( 8, label, cube);
743 Octree_labeling( 3, label, cube);
745 if( cube[15] == 1 ) {
747 Octree_labeling( 3, label, cube);
748 Octree_labeling( 7, label, cube);
749 Octree_labeling( 8, label, cube);
753 if( cube[16] == 1 ) {
755 Octree_labeling( 8, label, cube);
761 Octree_labeling( 1, label, cube);
763 if( cube[10] == 1 ) {
765 Octree_labeling( 1, label, cube);
766 Octree_labeling( 2, label, cube);
767 Octree_labeling( 6, label, cube);
769 if( cube[12] == 1 ) {
771 Octree_labeling( 1, label, cube);
772 Octree_labeling( 3, label, cube);
773 Octree_labeling( 7, label, cube);
777 if( cube[18] == 1 ) {
779 Octree_labeling( 6, label, cube);
781 if( cube[20] == 1 ) {
783 Octree_labeling( 7, label, cube);
785 if( cube[21] == 1 ) {
787 Octree_labeling( 6, label, cube);
788 Octree_labeling( 7, label, cube);
789 Octree_labeling( 8, label, cube);
793 if( cube[10] == 1 ) {
795 Octree_labeling( 1, label, cube);
796 Octree_labeling( 2, label, cube);
797 Octree_labeling( 5, label, cube);
799 if( cube[11] == 1 ) {
801 Octree_labeling( 2, label, cube);
803 if( cube[13] == 1 ) {
805 Octree_labeling( 2, label, cube);
806 Octree_labeling( 4, label, cube);
807 Octree_labeling( 8, label, cube);
809 if( cube[18] == 1 ) {
811 Octree_labeling( 5, label, cube);
813 if( cube[21] == 1 ) {
815 Octree_labeling( 5, label, cube);
816 Octree_labeling( 7, label, cube);
817 Octree_labeling( 8, label, cube);
821 if( cube[22] == 1 ) {
823 Octree_labeling( 8, label, cube);
827 if( cube[12] == 1 ) {
829 Octree_labeling( 1, label, cube);
830 Octree_labeling( 3, label, cube);
831 Octree_labeling( 5, label, cube);
833 if( cube[14] == 1 ) {
835 Octree_labeling( 3, label, cube);
837 if( cube[15] == 1 ) {
839 Octree_labeling( 3, label, cube);
840 Octree_labeling( 4, label, cube);
841 Octree_labeling( 8, label, cube);
843 if( cube[20] == 1 ) {
845 Octree_labeling( 5, label, cube);
847 if( cube[21] == 1 ) {
849 Octree_labeling( 5, label, cube);
850 Octree_labeling( 6, label, cube);
851 Octree_labeling( 8, label, cube);
855 if( cube[24] == 1 ) {
857 Octree_labeling( 8, label, cube);
861 if( cube[13] == 1 ) {
863 Octree_labeling( 2, label, cube);
864 Octree_labeling( 4, label, cube);
865 Octree_labeling( 6, label, cube);
867 if( cube[15] == 1 ) {
869 Octree_labeling( 3, label, cube);
870 Octree_labeling( 4, label, cube);
871 Octree_labeling( 7, label, cube);
873 if( cube[16] == 1 ) {
875 Octree_labeling( 4, label, cube);
877 if( cube[21] == 1 ) {
879 Octree_labeling( 5, label, cube);
880 Octree_labeling( 6, label, cube);
881 Octree_labeling( 7, label, cube);
883 if( cube[22] == 1 ) {
885 Octree_labeling( 6, label, cube);
887 if( cube[24] == 1 ) {
889 Octree_labeling( 7, label, cube);
900 template <class TInputImage,class TOutputImage>
902 BinaryThinningImageFilter3D<TInputImage,TOutputImage>
903 ::PrintSelf(std::ostream& os, Indent indent) const
905 Superclass::PrintSelf(os,indent);
907 os << indent << "Thinning image: " << std::endl;
911 } // end namespace itk