Medical Imaging Interaction Toolkit  2016.11.0
Medical Imaging Interaction Toolkit
itkSkeletonizationFilter.txx
Go to the documentation of this file.
1 /*===================================================================
2 
3 The Medical Imaging Interaction Toolkit (MITK)
4 
5 Copyright (c) German Cancer Research Center,
6 Division of Medical and Biological Informatics.
7 All rights reserved.
8 
9 This software is distributed WITHOUT ANY WARRANTY; without
10 even the implied warranty of MERCHANTABILITY or FITNESS FOR
11 A PARTICULAR PURPOSE.
12 
13 See LICENSE.txt or http://www.mitk.org for details.
14 
15 ===================================================================*/
16 #ifndef _itkSkeletonizationFilter_txx
17 #define _itkSkeletonizationFilter_txx
18 
19 #include "itkSkeletonizationFilter.h"
20 #include "mitkProgressBar.h"
21 #include <limits>
22 #include <time.h>
23 
24 namespace itk
25 {
26 
27  template< class TInputImage, class TOutputImage >
28  SkeletonizationFilter<TInputImage, TOutputImage>::SkeletonizationFilter()
29  {
30  m_DirectionImage = VectorImageType::New();
31  }
32 
33  template< class TInputImage, class TOutputImage >
34  SkeletonizationFilter<TInputImage, TOutputImage>::~SkeletonizationFilter()
35  {
36 
37  }
38 
39 
40 
41 
42  template< class TInputImage, class TOutputImage >
43  void SkeletonizationFilter<TInputImage, TOutputImage>::GenerateData()
44  {
45  //----------------------------------------------------------------------//
46  // Progress bar //
47  //----------------------------------------------------------------------//
48  mitk::ProgressBar::GetInstance()->AddStepsToDo( 3 );
49 
50  std::cout << "Skeletonize" << std::endl;
51 
52 
53  const InputImageType* faImage = this->GetInput();
54  typename InputImageType::SizeType size = faImage->GetRequestedRegion().GetSize();
55 
56  //typename RealImageType::SizeType size = m_FaImage->GetRequestedRegion().GetSize();
57 
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);
64 
65 
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++)
67  {
68  typename InputImageType::IndexType ix;
69  ix[0]=x; ix[1]=y; ix[2]=z;
70  float theval = faImage->GetPixel(ix);
71 
72  if(theval != 0)
73  {
74 
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/
77  */
78 
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;
81 
82  for(int dz=-1; dz<=1; dz++) for(int dy=-1; dy<=1; dy++) for(int dx=-1; dx<=1;dx++)
83  {
84  typename InputImageType::IndexType p;
85  p[0] = x+dx; p[1] = y+dy; p[2] = z+dz;
86  float mass = faImage->GetPixel(p);
87 
88  sum += mass;
89  cogX += (float)dx*mass; cogY += (float)dy*mass; cogZ += (float)dz*mass;
90  }
91 
92  cogX /= sum; cogY /= sum; cogZ /= sum;
93  l = sqrt(cogX*cogX + cogY*cogY + cogZ*cogZ);
94 
95  if (l > 0.1) /* is CofG far enough away from centre voxel? */
96  {
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);
100  }
101  else
102  // Find direction of max curvature
103  {
104 
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
107  {
108  for(int yy=-1; yy<=1; yy++)
109  {
110  for(int xx=-1; xx<=1; xx++)
111  {
112  if ( (zz==1) || (yy==1) || ((yy==0)&&(xx==1)) )
113  {
114  float weighting = pow( (float)(xx*xx+yy*yy+zz*zz) , (float)-0.7 ); // power is arbitrary: maybe test other functions here
115 
116 
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));
123 
124  if (cost>maxcost)
125  {
126  maxcost=cost;
127  vecX=xx;
128  vecY=yy;
129  vecZ=zz;
130  }
131  }
132  }
133  }
134  }
135  }
136 
137  VectorType vec;
138  vec[0] = vecX; vec[1] = vecY; vec[2]=vecZ;
139  m_DirectionImage->SetPixel(ix, vec);
140 
141  }
142  }
143 
144  mitk::ProgressBar::GetInstance()->Progress();
145 
146 
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();
155 
156  VectorImageType::PixelType p;
157  p[0]=0; p[1]=0; p[2]=0;
158  directionSmoothed->FillBuffer(p);
159 
160 
161 
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++)
163  {
164  VectorImageType::IndexType ix;
165  ix[0]=x; ix[1]=y; ix[2]=z;
166 
167  // Find the vector that occured most
168  auto localsum = new int[27];
169  int localmax=0, xxx, yyy, zzz;
170 
171  for(int zz=0; zz<27; zz++) localsum[zz]=0;
172 
173  for(int zz=-1; zz<=1; zz++) for(int yy=-1; yy<=1; yy++) for(int xx=-1; xx<=1; xx++)
174  {
175  VectorImageType::IndexType i;
176  i[0] = x+xx; i[1] = y+yy; i[2] = z+zz;
177  VectorType v = m_DirectionImage->GetPixel(i);
178  xxx = v[0];
179  yyy = v[1];
180  zzz = v[2];
181 
182  localsum[(1+zzz)*9+(1+yyy)*3+1+xxx]++;
183  localsum[(1-zzz)*9+(1-yyy)*3+1-xxx]++;
184  }
185 
186  for(int zz=-1; zz<=1; zz++) for(int yy=-1; yy<=1; yy++) for(int xx=-1; xx<=1; xx++)
187  {
188  if (localsum[(1+zz)*9+(1+yy)*3+1+xx]>localmax)
189  {
190  localmax=localsum[(1+zz)*9+(1+yy)*3+1+xx];
191  VectorType v;
192  v[0] = xx; v[1] = yy; v[2] = zz;
193  directionSmoothed->SetPixel(ix, v);
194  }
195  }
196 
197  delete localsum;
198 
199  }
200 
201  m_DirectionImage = directionSmoothed;
202 
203  mitk::ProgressBar::GetInstance()->Progress();
204 
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);
213 
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++)
215  {
216 
217  typename InputImageType::IndexType ix;
218  ix[0]=x; ix[1]=y; ix[2]=z;
219 
220  float theval = faImage->GetPixel(ix);
221  VectorType v = directionSmoothed->GetPixel(ix);
222 
223  typename VectorImageType::IndexType i;
224 
225  i[0] = x-v[0]; i[1] = y-v[1]; i[2] = z-v[2];
226  float min = faImage->GetPixel(i);
227 
228  i[0] = x+v[0]; i[1] = y+v[1]; i[2] = z+v[2];
229  float plus = faImage->GetPixel(i);
230 
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);
233 
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);
236 
237  if( ((v[0]!=0) || (v[1]!=0) || (v[2]!=0)) &&
238  theval >= plus && theval > min && theval >= plusplus && theval > minmin )
239  {
240  outputImg->SetPixel(ix, theval);
241  }
242 
243  }
244 
245 
246 
247  Superclass::SetNthOutput( 0, outputImg );
248  mitk::ProgressBar::GetInstance()->Progress();
249 
250  }
251 
252 
253 
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()
257  {
258  GradientImageType::Pointer gradImg = GradientImageType::New();
259 
260  if(m_DirectionImage.IsNotNull())
261  {
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);
267  gradImg->Allocate();
268 
269 
270  VectorImageType::SizeType size = m_DirectionImage->GetLargestPossibleRegion().GetSize();
271 
272  for(std::size_t i=0; i<size[0]; i++)
273  {
274  for(std::size_t j=0; j<size[1]; j++)
275  {
276  for(std::size_t k=0; k<size[2]; k++)
277  {
278  itk::Index<3> ix;
279  ix[0] = i;
280  ix[1] = j;
281  ix[2] = k;
282 
283  VectorType vec = m_DirectionImage->GetPixel(ix);
284 
285  itk::VariableLengthVector<int> pixel;
286  pixel.SetSize(3);
287  pixel.SetElement(0, vec.GetElement(0));
288  pixel.SetElement(1, vec.GetElement(1));
289  pixel.SetElement(2, vec.GetElement(2));
290 
291  gradImg->SetPixel(ix, pixel);
292 
293  }
294  }
295  }
296 
297  }
298 
299 
300 
301 
302  return gradImg;
303  }
304 
305 }
306 #endif // _itkSkeletonizationFilter_txx