1 /*===================================================================
3 The Medical Imaging Interaction Toolkit (MITK)
5 Copyright (c) German Cancer Research Center,
6 Division of Medical and Biological Informatics.
9 This software is distributed WITHOUT ANY WARRANTY; without
10 even the implied warranty of MERCHANTABILITY or FITNESS FOR
13 See LICENSE.txt or http://www.mitk.org for details.
15 ===================================================================*/
16 #ifndef _itkSkeletonizationFilter_txx
17 #define _itkSkeletonizationFilter_txx
19 #include "itkSkeletonizationFilter.h"
20 #include "mitkProgressBar.h"
27 template< class TInputImage, class TOutputImage >
28 SkeletonizationFilter<TInputImage, TOutputImage>::SkeletonizationFilter()
30 m_DirectionImage = VectorImageType::New();
33 template< class TInputImage, class TOutputImage >
34 SkeletonizationFilter<TInputImage, TOutputImage>::~SkeletonizationFilter()
42 template< class TInputImage, class TOutputImage >
43 void SkeletonizationFilter<TInputImage, TOutputImage>::GenerateData()
45 //----------------------------------------------------------------------//
47 //----------------------------------------------------------------------//
48 mitk::ProgressBar::GetInstance()->AddStepsToDo( 3 );
50 std::cout << "Skeletonize" << std::endl;
53 const InputImageType* faImage = this->GetInput();
54 typename InputImageType::SizeType size = faImage->GetRequestedRegion().GetSize();
56 //typename RealImageType::SizeType size = m_FaImage->GetRequestedRegion().GetSize();
58 m_DirectionImage->SetRegions(faImage->GetRequestedRegion());
59 m_DirectionImage->SetDirection(faImage->GetDirection());
60 m_DirectionImage->SetSpacing(faImage->GetSpacing());
61 m_DirectionImage->SetOrigin(faImage->GetOrigin());
62 m_DirectionImage->Allocate();
63 m_DirectionImage->FillBuffer(0.0);
66 for(unsigned int z=1; z<size[2]-1; z++) for(unsigned int y=1; y<size[1]-1; y++) for(unsigned int x=1; x<size[0]-1; x++)
68 typename InputImageType::IndexType ix;
69 ix[0]=x; ix[1]=y; ix[2]=z;
70 float theval = faImage->GetPixel(ix);
75 /* Calculate point of gravity. We will consider each 3x3x3 neighbourhood as a unit cube. The center
76 * point of each voxel will be a multiplicative of 1/6. The center of the unit cube is 3/6 = 1/2/
79 float cogX = 0.0; float cogY = 0.0; float cogZ = 0.0; float sum = 0.0; float l;
80 int vecX = 0; int vecY = 0; int vecZ = 0;
82 for(int dz=-1; dz<=1; dz++) for(int dy=-1; dy<=1; dy++) for(int dx=-1; dx<=1;dx++)
84 typename InputImageType::IndexType p;
85 p[0] = x+dx; p[1] = y+dy; p[2] = z+dz;
86 float mass = faImage->GetPixel(p);
89 cogX += (float)dx*mass; cogY += (float)dy*mass; cogZ += (float)dz*mass;
92 cogX /= sum; cogY /= sum; cogZ /= sum;
93 l = sqrt(cogX*cogX + cogY*cogY + cogZ*cogZ);
95 if (l > 0.1) /* is CofG far enough away from centre voxel? */
97 vecX = std::max(std::min(round(cogX/l),1),-1);
98 vecY = std::max(std::min(round(cogY/l),1),-1);
99 vecZ = std::max(std::min(round(cogZ/l),1),-1);
102 // Find direction of max curvature
105 float maxcost=0, centreval=2*theval;
106 for(int zz=0; zz<=1; zz++) // note - starts at zero as we're only searching half the voxels
108 for(int yy=-1; yy<=1; yy++)
110 for(int xx=-1; xx<=1; xx++)
112 if ( (zz==1) || (yy==1) || ((yy==0)&&(xx==1)) )
114 float weighting = pow( (float)(xx*xx+yy*yy+zz*zz) , (float)-0.7 ); // power is arbitrary: maybe test other functions here
117 typename InputImageType::IndexType i,j;
118 i[0] = x+xx; i[1] = y+yy; i[2] = z+zz;
119 j[0] = x-xx; j[1] = y-yy; j[2] = z-zz;
120 float cost = weighting * ( centreval
121 - (float)faImage->GetPixel(i)
122 - (float)faImage->GetPixel(j));
138 vec[0] = vecX; vec[1] = vecY; vec[2]=vecZ;
139 m_DirectionImage->SetPixel(ix, vec);
144 mitk::ProgressBar::GetInstance()->Progress();
147 // Smooth m_DirectionImage and store in directionSmoothed by finding the
148 // mode in a 3*3 neighbourhoud
149 VectorImageType::Pointer directionSmoothed = VectorImageType::New();
150 directionSmoothed->SetRegions(faImage->GetRequestedRegion());
151 directionSmoothed->SetDirection(faImage->GetDirection());
152 directionSmoothed->SetSpacing(faImage->GetSpacing());
153 directionSmoothed->SetOrigin(faImage->GetOrigin());
154 directionSmoothed->Allocate();
156 VectorImageType::PixelType p;
157 p[0]=0; p[1]=0; p[2]=0;
158 directionSmoothed->FillBuffer(p);
162 for(unsigned int z=1; z<size[2]-1; z++) for(unsigned int y=1; y<size[1]-1; y++) for(unsigned int x=1; x<size[0]-1; x++)
164 VectorImageType::IndexType ix;
165 ix[0]=x; ix[1]=y; ix[2]=z;
167 // Find the vector that occured most
168 auto localsum = new int[27];
169 int localmax=0, xxx, yyy, zzz;
171 for(int zz=0; zz<27; zz++) localsum[zz]=0;
173 for(int zz=-1; zz<=1; zz++) for(int yy=-1; yy<=1; yy++) for(int xx=-1; xx<=1; xx++)
175 VectorImageType::IndexType i;
176 i[0] = x+xx; i[1] = y+yy; i[2] = z+zz;
177 VectorType v = m_DirectionImage->GetPixel(i);
182 localsum[(1+zzz)*9+(1+yyy)*3+1+xxx]++;
183 localsum[(1-zzz)*9+(1-yyy)*3+1-xxx]++;
186 for(int zz=-1; zz<=1; zz++) for(int yy=-1; yy<=1; yy++) for(int xx=-1; xx<=1; xx++)
188 if (localsum[(1+zz)*9+(1+yy)*3+1+xx]>localmax)
190 localmax=localsum[(1+zz)*9+(1+yy)*3+1+xx];
192 v[0] = xx; v[1] = yy; v[2] = zz;
193 directionSmoothed->SetPixel(ix, v);
201 m_DirectionImage = directionSmoothed;
203 mitk::ProgressBar::GetInstance()->Progress();
205 // Do non-max-suppression in the direction of perp and set as output of the filter
206 typename OutputImageType::Pointer outputImg = OutputImageType::New();
207 outputImg->SetRegions(faImage->GetRequestedRegion());
208 outputImg->SetDirection(faImage->GetDirection());
209 outputImg->SetSpacing(faImage->GetSpacing());
210 outputImg->SetOrigin(faImage->GetOrigin());
211 outputImg->Allocate();
212 outputImg->FillBuffer(0.0);
214 for(unsigned int z=1; z<size[2]-1; z++) for(unsigned int y=1; y<size[1]-1; y++) for(unsigned int x=1; x<size[0]-1; x++)
217 typename InputImageType::IndexType ix;
218 ix[0]=x; ix[1]=y; ix[2]=z;
220 float theval = faImage->GetPixel(ix);
221 VectorType v = directionSmoothed->GetPixel(ix);
223 typename VectorImageType::IndexType i;
225 i[0] = x-v[0]; i[1] = y-v[1]; i[2] = z-v[2];
226 float min = faImage->GetPixel(i);
228 i[0] = x+v[0]; i[1] = y+v[1]; i[2] = z+v[2];
229 float plus = faImage->GetPixel(i);
231 i[0] = x-2*v[0]; i[1] = y-2*v[1]; i[2] = z-2*v[2];
232 float minmin = faImage->GetPixel(i);
234 i[0] = x+2*v[0]; i[1] = y+2*v[1]; i[2] = z+2*v[2];
235 float plusplus = faImage->GetPixel(i);
237 if( ((v[0]!=0) || (v[1]!=0) || (v[2]!=0)) &&
238 theval >= plus && theval > min && theval >= plusplus && theval > minmin )
240 outputImg->SetPixel(ix, theval);
247 Superclass::SetNthOutput( 0, outputImg );
248 mitk::ProgressBar::GetInstance()->Progress();
254 // Can provide a vector image to visualize the gradient image used in the search for local maxima.
255 template< class TInputImage, class TOutputImage >
256 itk::VectorImage<int, 3>::Pointer SkeletonizationFilter<TInputImage, TOutputImage>::GetGradientImage()
258 GradientImageType::Pointer gradImg = GradientImageType::New();
260 if(m_DirectionImage.IsNotNull())
262 gradImg->SetSpacing(m_DirectionImage->GetSpacing());
263 gradImg->SetOrigin(m_DirectionImage->GetOrigin());
264 gradImg->SetDirection(m_DirectionImage->GetDirection());
265 gradImg->SetRegions(m_DirectionImage->GetLargestPossibleRegion().GetSize());
266 gradImg->SetVectorLength(3);
270 VectorImageType::SizeType size = m_DirectionImage->GetLargestPossibleRegion().GetSize();
272 for(std::size_t i=0; i<size[0]; i++)
274 for(std::size_t j=0; j<size[1]; j++)
276 for(std::size_t k=0; k<size[2]; k++)
283 VectorType vec = m_DirectionImage->GetPixel(ix);
285 itk::VariableLengthVector<int> pixel;
287 pixel.SetElement(0, vec.GetElement(0));
288 pixel.SetElement(1, vec.GetElement(1));
289 pixel.SetElement(2, vec.GetElement(2));
291 gradImg->SetPixel(ix, pixel);
306 #endif // _itkSkeletonizationFilter_txx