Medical Imaging Interaction Toolkit  2016.11.0
Medical Imaging Interaction Toolkit
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
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