Medical Imaging Interaction Toolkit  2016.11.0
Medical Imaging Interaction Toolkit
itkFiniteDiffOdfMaximaExtractionFilter.cpp
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 __itkFiniteDiffOdfMaximaExtractionFilter_cpp
17 #define __itkFiniteDiffOdfMaximaExtractionFilter_cpp
18 
19 #include "itkFiniteDiffOdfMaximaExtractionFilter.h"
20 #include <itkImageRegionConstIterator.h>
21 #include <itkImageRegionConstIteratorWithIndex.h>
22 #include <itkImageRegionIterator.h>
23 #include <vnl/vnl_vector.h>
24 #include <itkOrientationDistributionFunction.h>
25 #include <itkContinuousIndex.h>
26 
27 #include <vtkSmartPointer.h>
28 #include <vtkPolyData.h>
29 #include <vtkCellArray.h>
30 #include <vtkPoints.h>
31 #include <vtkPolyLine.h>
32 
33 #include <boost/math/special_functions.hpp>
34 #include <boost/progress.hpp>
35 
36 
37 using namespace boost::math;
38 
39 namespace itk {
40 
41 
42 static bool CompareVectors(const vnl_vector_fixed< double, 3 >& v1, const vnl_vector_fixed< double, 3 >& v2)
43 {
44  return (v1.magnitude()>v2.magnitude());
45 }
46 
47 template< class PixelType, int ShOrder, int NrOdfDirections >
48 FiniteDiffOdfMaximaExtractionFilter< PixelType, ShOrder, NrOdfDirections>
49 ::FiniteDiffOdfMaximaExtractionFilter()
50  : m_NormalizationMethod(MAX_VEC_NORM)
51  , m_MaxNumPeaks(2)
52  , m_PeakThreshold(0.4)
53  , m_AbsolutePeakThreshold(0)
54  , m_ClusteringThreshold(0.9)
55  , m_AngularThreshold(0.7)
56  , m_NumCoeffs((ShOrder*ShOrder + ShOrder + 2)/2 + ShOrder)
57  , m_Toolkit(FSL)
58 {
59  this->SetNumberOfRequiredInputs(1);
60 }
61 
62 template< class PixelType, int ShOrder, int NrOdfDirections >
63 void FiniteDiffOdfMaximaExtractionFilter< PixelType, ShOrder, NrOdfDirections>
64 ::FindCandidatePeaks(OdfType& odf, double thr, std::vector< DirectionType >& container)
65 {
66  double gfa = odf.GetGeneralizedFractionalAnisotropy();
67  //Find the peaks using a finite difference method
68  bool flag = true;
69  vnl_vector_fixed< bool, NrOdfDirections > used; used.fill(false);
70  //Find the peaks
71  for (int i=0; i<NrOdfDirections; i++)
72  {
73  if (used[i])
74  continue;
75 
76  double val = odf.GetElement(i);
77  if (val>thr && val*gfa>m_AbsolutePeakThreshold) // limit to one hemisphere ???
78  {
79  flag = true;
80  std::vector< int > neighbours = odf.GetNeighbors(i);
81  for (unsigned int j=0; j<neighbours.size(); j++)
82  if (val<=odf.GetElement(neighbours.at(j)))
83  {
84  flag = false;
85  break;
86  }
87  if (flag) // point is a peak
88  {
89  container.push_back(odf.GetDirection(i).normalize());
90  used[i] = true;
91  for (unsigned int j=0; j<neighbours.size(); j++)
92  used[neighbours.at(j)] = true;
93  }
94  }
95  }
96 }
97 
98 template< class PixelType, int ShOrder, int NrOdfDirections >
99 std::vector< vnl_vector_fixed< double, 3 > > FiniteDiffOdfMaximaExtractionFilter< PixelType, ShOrder, NrOdfDirections>::MeanShiftClustering(std::vector< vnl_vector_fixed< double, 3 > >& inDirs)
100 {
101  std::vector< DirectionType > outDirs;
102  if (inDirs.empty())
103  return inDirs;
104 
105  DirectionType oldMean, currentMean, workingMean;
106  std::vector< int > touched;
107 
108  // initialize
109  touched.resize(inDirs.size(), 0);
110  bool free = true;
111  currentMean = inDirs[0]; // initialize first seed
112  while (free)
113  {
114  oldMean.fill(0.0);
115 
116  // start mean-shift clustering
117  float angle = 0.0;
118  int counter = 0;
119  while ((currentMean-oldMean).magnitude()>0.0001)
120  {
121  counter = 0;
122  oldMean = currentMean;
123  workingMean = oldMean;
124  workingMean.normalize();
125  currentMean.fill(0.0);
126  for (unsigned int i=0; i<inDirs.size(); i++)
127  {
128  angle = dot_product(workingMean, inDirs[i]);
129  if (angle>=m_ClusteringThreshold)
130  {
131  currentMean += inDirs[i];
132  touched[i] = 1;
133  counter++;
134  }
135  else if (-angle>=m_ClusteringThreshold)
136  {
137  currentMean -= inDirs[i];
138  touched[i] = 1;
139  counter++;
140  }
141  }
142  }
143 
144  // found stable mean
145  if (counter>0)
146  {
147  float mag = currentMean.magnitude();
148  if (mag>0)
149  {
150  currentMean /= mag;
151  outDirs.push_back(currentMean);
152  }
153  }
154 
155  // find next unused seed
156  free = false;
157  for (unsigned int i=0; i<touched.size(); i++)
158  if (touched[i]==0)
159  {
160  currentMean = inDirs[i];
161  free = true;
162  }
163  }
164 
165  if (inDirs.size()==outDirs.size())
166  {
167  return outDirs;
168  }
169  else
170  return MeanShiftClustering(outDirs);
171 }
172 
173 template< class PixelType, int ShOrder, int NrOdfDirections >
174 void FiniteDiffOdfMaximaExtractionFilter< PixelType, ShOrder, NrOdfDirections>
175 ::BeforeThreadedGenerateData()
176 {
177  typename CoefficientImageType::Pointer ShCoeffImage = static_cast< CoefficientImageType* >( this->ProcessObject::GetInput(0) );
178  itk::Vector<double,3> spacing = ShCoeffImage->GetSpacing();
179  double minSpacing = spacing[0];
180  if (spacing[1]<minSpacing)
181  minSpacing = spacing[1];
182  if (spacing[2]<minSpacing)
183  minSpacing = spacing[2];
184 
185  mitk::Point3D origin = ShCoeffImage->GetOrigin();
186  itk::Matrix<double, 3, 3> direction = ShCoeffImage->GetDirection();
187  ImageRegion<3> imageRegion = ShCoeffImage->GetLargestPossibleRegion();
188 
189  if (m_MaskImage.IsNotNull())
190  {
191  origin = m_MaskImage->GetOrigin();
192  direction = m_MaskImage->GetDirection();
193  imageRegion = m_MaskImage->GetLargestPossibleRegion();
194  }
195 
196  m_DirectionImageContainer = ItkDirectionImageContainer::New();
197  for (unsigned int i=0; i<m_MaxNumPeaks; i++)
198  {
199  itk::Vector< float, 3 > nullVec; nullVec.Fill(0.0);
200  ItkDirectionImage::Pointer img = ItkDirectionImage::New();
201  img->SetSpacing( spacing );
202  img->SetOrigin( origin );
203  img->SetDirection( direction );
204  img->SetRegions( imageRegion );
205  img->Allocate();
206  img->FillBuffer(nullVec);
207  m_DirectionImageContainer->InsertElement(m_DirectionImageContainer->Size(), img);
208  }
209  if (m_MaskImage.IsNull())
210  {
211  m_MaskImage = ItkUcharImgType::New();
212  m_MaskImage->SetSpacing( spacing );
213  m_MaskImage->SetOrigin( origin );
214  m_MaskImage->SetDirection( direction );
215  m_MaskImage->SetRegions( imageRegion );
216  m_MaskImage->Allocate();
217  m_MaskImage->FillBuffer(1);
218  }
219  m_NumDirectionsImage = ItkUcharImgType::New();
220  m_NumDirectionsImage->SetSpacing( spacing );
221  m_NumDirectionsImage->SetOrigin( origin );
222  m_NumDirectionsImage->SetDirection( direction );
223  m_NumDirectionsImage->SetRegions( imageRegion );
224  m_NumDirectionsImage->Allocate();
225  m_NumDirectionsImage->FillBuffer(0);
226 
227  this->SetNumberOfIndexedOutputs(m_MaxNumPeaks);
228 
229  // calculate SH basis
230  OdfType odf;
231  vnl_matrix_fixed<double, 3, NrOdfDirections>* directions = odf.GetDirections();
232  vnl_matrix< double > sphCoords;
233  std::vector< DirectionType > dirs;
234  for (int i=0; i<NrOdfDirections; i++)
235  dirs.push_back(directions->get_column(i));
236  Cart2Sph(dirs, sphCoords); // convert candidate peaks to spherical angles
237  m_ShBasis = CalcShBasis(sphCoords); // evaluate spherical harmonics at each peak
238 
239  MITK_INFO << "Starting finite differences maximum extraction";
240  MITK_INFO << "ODF sampling points: " << NrOdfDirections;
241  MITK_INFO << "SH order: " << ShOrder;
242  MITK_INFO << "Maximum peaks: " << m_MaxNumPeaks;
243  MITK_INFO << "Relative threshold: " << m_PeakThreshold;
244  MITK_INFO << "Absolute threshold: " << m_AbsolutePeakThreshold;
245  MITK_INFO << "Clustering threshold: " << m_ClusteringThreshold;
246  MITK_INFO << "Angular threshold: " << m_AngularThreshold;
247 }
248 
249 template< class PixelType, int ShOrder, int NrOdfDirections >
250 void FiniteDiffOdfMaximaExtractionFilter< PixelType, ShOrder, NrOdfDirections>
251 ::AfterThreadedGenerateData()
252 {
253  MITK_INFO << "Generating vector field";
254  vtkSmartPointer<vtkCellArray> m_VtkCellArray = vtkSmartPointer<vtkCellArray>::New();
255  vtkSmartPointer<vtkPoints> m_VtkPoints = vtkSmartPointer<vtkPoints>::New();
256 
257  typename CoefficientImageType::Pointer ShCoeffImage = static_cast< CoefficientImageType* >( this->ProcessObject::GetInput(0) );
258  ImageRegionConstIterator< CoefficientImageType > cit(ShCoeffImage, ShCoeffImage->GetLargestPossibleRegion() );
259 
260  mitk::Vector3D spacing = ShCoeffImage->GetSpacing();
261  double minSpacing = spacing[0];
262  if (spacing[1]<minSpacing)
263  minSpacing = spacing[1];
264  if (spacing[2]<minSpacing)
265  minSpacing = spacing[2];
266 
267  int maxProgress = ShCoeffImage->GetLargestPossibleRegion().GetSize()[0]*ShCoeffImage->GetLargestPossibleRegion().GetSize()[1]*ShCoeffImage->GetLargestPossibleRegion().GetSize()[2];
268  boost::progress_display disp(maxProgress);
269 
270  while( !cit.IsAtEnd() )
271  {
272  ++disp;
273 
274  typename CoefficientImageType::IndexType index = cit.GetIndex();
275  if (m_MaskImage->GetPixel(index)==0)
276  {
277  ++cit;
278  continue;
279  }
280 
281  for (unsigned int i=0; i<m_DirectionImageContainer->Size(); i++)
282  {
283  ItkDirectionImage::Pointer img = m_DirectionImageContainer->GetElement(i);
284  itk::Vector< float, 3 > pixel = img->GetPixel(index);
285  DirectionType dir; dir[0] = pixel[0]; dir[1] = pixel[1]; dir[2] = pixel[2];
286 
287  vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
288  itk::ContinuousIndex<double, 3> center;
289  center[0] = index[0];
290  center[1] = index[1];
291  center[2] = index[2];
292  itk::Point<double> worldCenter;
293  m_MaskImage->TransformContinuousIndexToPhysicalPoint( center, worldCenter );
294 
295  itk::Point<double> worldStart;
296  worldStart[0] = worldCenter[0]-dir[0]/2 * minSpacing;
297  worldStart[1] = worldCenter[1]-dir[1]/2 * minSpacing;
298  worldStart[2] = worldCenter[2]-dir[2]/2 * minSpacing;
299  vtkIdType id = m_VtkPoints->InsertNextPoint(worldStart.GetDataPointer());
300  container->GetPointIds()->InsertNextId(id);
301  itk::Point<double> worldEnd;
302  worldEnd[0] = worldCenter[0]+dir[0]/2 * minSpacing;
303  worldEnd[1] = worldCenter[1]+dir[1]/2 * minSpacing;
304  worldEnd[2] = worldCenter[2]+dir[2]/2 * minSpacing;
305  id = m_VtkPoints->InsertNextPoint(worldEnd.GetDataPointer());
306  container->GetPointIds()->InsertNextId(id);
307  m_VtkCellArray->InsertNextCell(container);
308  }
309  ++cit;
310  }
311 
312  vtkSmartPointer<vtkPolyData> directionsPolyData = vtkSmartPointer<vtkPolyData>::New();
313  directionsPolyData->SetPoints(m_VtkPoints);
314  directionsPolyData->SetLines(m_VtkCellArray);
315  m_OutputFiberBundle = mitk::FiberBundle::New(directionsPolyData);
316 
317  for (unsigned int i=0; i<m_DirectionImageContainer->Size(); i++)
318  {
319  ItkDirectionImage::Pointer img = m_DirectionImageContainer->GetElement(i);
320  this->SetNthOutput(i, img);
321  }
322 }
323 
324 template< class PixelType, int ShOrder, int NrOdfDirections >
325 void FiniteDiffOdfMaximaExtractionFilter< PixelType, ShOrder, NrOdfDirections>
326 ::ThreadedGenerateData( const OutputImageRegionType& outputRegionForThread, ThreadIdType threadID )
327 {
328  typename CoefficientImageType::Pointer ShCoeffImage = static_cast< CoefficientImageType* >( this->ProcessObject::GetInput(0) );
329 
330  ImageRegionConstIterator< CoefficientImageType > cit(ShCoeffImage, outputRegionForThread );
331 
332  OdfType odf;
333  while( !cit.IsAtEnd() )
334  {
335  typename CoefficientImageType::IndexType index = cit.GetIndex();
336  if (m_MaskImage->GetPixel(index)==0)
337  {
338  ++cit;
339  continue;
340  }
341 
342  CoefficientPixelType c = cit.Get();
343 
344  // calculate ODF
345  double max = 0;
346  odf.Fill(0.0);
347  for (int i=0; i<NrOdfDirections; i++)
348  {
349  for (int j=0; j<m_NumCoeffs; j++)
350  odf[i] += c[j]*m_ShBasis(i,j);
351  if (odf[i]>max)
352  max = odf[i];
353  }
354  if (max<0.0001)
355  {
356  ++cit;
357  continue;
358  }
359 
360  std::vector< DirectionType > candidates, peaks, temp;
361  peaks.clear();
362  max *= m_PeakThreshold; // relative threshold
363  FindCandidatePeaks(odf, max, candidates); // find all local maxima
364  candidates = MeanShiftClustering(candidates); // cluster maxima
365 
366  vnl_matrix< double > shBasis, sphCoords;
367  Cart2Sph(candidates, sphCoords); // convert candidate peaks to spherical angles
368  shBasis = CalcShBasis(sphCoords); // evaluate spherical harmonics at each peak
369  max = 0.0;
370  for (unsigned int i=0; i<candidates.size(); i++) // scale peaks according to ODF value
371  {
372  double val = 0;
373  for (int j=0; j<m_NumCoeffs; j++)
374  val += c[j]*shBasis(i,j);
375  if (val>max)
376  max = val;
377  peaks.push_back(candidates[i]*val);
378  }
379  std::sort( peaks.begin(), peaks.end(), CompareVectors ); // sort peaks
380 
381  // kick out directions to close to a larger direction (too far away to cluster but too close to keep)
382  unsigned int m = peaks.size();
383  if ( m>m_DirectionImageContainer->Size() )
384  m = m_DirectionImageContainer->Size();
385  for (unsigned int i=0; i<m; i++)
386  {
387  DirectionType v1 = peaks.at(i);
388  double val = v1.magnitude();
389  if (val<max*m_PeakThreshold || val<m_AbsolutePeakThreshold)
390  break;
391 
392  bool flag = true;
393  for (unsigned int j=0; j<peaks.size(); j++)
394  if (i!=j)
395  {
396  DirectionType v2 = peaks.at(j);
397  double val2 = v2.magnitude();
398  double angle = fabs(dot_product(v1,v2)/(val*val2));
399  if (angle>m_AngularThreshold && val<val2)
400  {
401  flag = false;
402  break;
403  }
404  }
405 
406  if (flag)
407  temp.push_back(v1);
408  }
409  peaks = temp;
410 
411  // fill output image
412  unsigned int num = peaks.size();
413  if ( num>m_DirectionImageContainer->Size() )
414  num = m_DirectionImageContainer->Size();
415  for (unsigned int i=0; i<num; i++)
416  {
417  vnl_vector<double> dir = peaks.at(i);
418 
419  ItkDirectionImage::Pointer img = m_DirectionImageContainer->GetElement(i);
420 
421  switch (m_NormalizationMethod)
422  {
423  case NO_NORM:
424  break;
425  case SINGLE_VEC_NORM:
426  dir.normalize();
427  break;
428  case MAX_VEC_NORM:
429  dir /= max;
430  break;
431  }
432 
433  dir = m_MaskImage->GetDirection()*dir;
434 
435  if (m_FlipX)
436  dir[0] = -dir[0];
437  if (m_FlipY)
438  dir[1] = -dir[1];
439  if (m_FlipZ)
440  dir[2] = -dir[2];
441 
442  itk::Vector< float, 3 > pixel;
443  pixel.SetElement(0, dir[0]);
444  pixel.SetElement(1, dir[1]);
445  pixel.SetElement(2, dir[2]);
446  img->SetPixel(index, pixel);
447  }
448  m_NumDirectionsImage->SetPixel(index, num);
449  ++cit;
450  }
451  MITK_INFO << "Thread " << threadID << " finished extraction";
452 }
453 
454 // convert cartesian to spherical coordinates
455 template< class PixelType, int ShOrder, int NrOdfDirections >
456 void FiniteDiffOdfMaximaExtractionFilter< PixelType, ShOrder, NrOdfDirections>
457 ::Cart2Sph(const std::vector< DirectionType >& dir, vnl_matrix<double>& sphCoords)
458 {
459  sphCoords.set_size(dir.size(), 2);
460 
461  for (unsigned int i=0; i<dir.size(); i++)
462  {
463  double mag = dir[i].magnitude();
464 
465  if( mag<0.0001 )
466  {
467  sphCoords(i,0) = M_PI/2; // theta
468  sphCoords(i,1) = M_PI/2; // phi
469  }
470  else
471  {
472  sphCoords(i,0) = acos(dir[i](2)/mag); // theta
473  sphCoords(i,1) = atan2(dir[i](1), dir[i](0)); // phi
474  }
475  }
476 }
477 
478 // generate spherical harmonic values of the desired order for each input direction
479 template< class PixelType, int ShOrder, int NrOdfDirections >
480 vnl_matrix<double> FiniteDiffOdfMaximaExtractionFilter< PixelType, ShOrder, NrOdfDirections>
481 ::CalcShBasis(vnl_matrix<double>& sphCoords)
482 {
483  int M = sphCoords.rows();
484  int j, m; double mag, plm;
485  vnl_matrix<double> shBasis;
486  shBasis.set_size(M, m_NumCoeffs);
487 
488  for (int p=0; p<M; p++)
489  {
490  j=0;
491  for (int l=0; l<=ShOrder; l=l+2)
492  for (m=-l; m<=l; m++)
493  {
494  switch (m_Toolkit)
495  {
496  case FSL:
497  plm = legendre_p<double>(l,abs(m),cos(sphCoords(p,0)));
498  mag = sqrt((double)(2*l+1)/(4.0*M_PI)*factorial<double>(l-abs(m))/factorial<double>(l+abs(m)))*plm;
499 
500  if (m<0)
501  shBasis(p,j) = sqrt(2.0)*mag*cos(fabs((double)m)*sphCoords(p,1));
502  else if (m==0)
503  shBasis(p,j) = mag;
504  else
505  shBasis(p,j) = pow(-1.0, m)*sqrt(2.0)*mag*sin(m*sphCoords(p,1));
506  break;
507  case MRTRIX:
508 
509  plm = legendre_p<double>(l,abs(m),-cos(sphCoords(p,0)));
510  mag = sqrt((double)(2*l+1)/(4.0*M_PI)*factorial<double>(l-abs(m))/factorial<double>(l+abs(m)))*plm;
511  if (m>0)
512  shBasis(p,j) = mag*cos(m*sphCoords(p,1));
513  else if (m==0)
514  shBasis(p,j) = mag;
515  else
516  shBasis(p,j) = mag*sin(-m*sphCoords(p,1));
517  break;
518  }
519 
520  j++;
521  }
522  }
523  return shBasis;
524 }
525 
526 }
527 
528 #endif // __itkFiniteDiffOdfMaximaExtractionFilter_cpp
itk::FiniteDiffOdfMaximaExtractionFilter::CalcShBasis
vnl_matrix< double > CalcShBasis(vnl_matrix< double > &sphCoords)
Definition: itkFiniteDiffOdfMaximaExtractionFilter.cpp:481
mitk::Pointer
itk::SmartPointer< Self > Pointer
Definition: mitkRenderingManager.h:389
MITK_INFO
#define MITK_INFO
Definition: mitkLogMacros.h:22
itk::FiniteDiffOdfMaximaExtractionFilter::MeanShiftClustering
std::vector< DirectionType > MeanShiftClustering(std::vector< DirectionType > &inDirs)
Definition: itkFiniteDiffOdfMaximaExtractionFilter.cpp:99
itk::FiniteDiffOdfMaximaExtractionFilter::FindCandidatePeaks
void FindCandidatePeaks(OdfType &odf, double odfMax, std::vector< DirectionType > &inDirs)
Definition: itkFiniteDiffOdfMaximaExtractionFilter.cpp:64
itk::OrientationDistributionFunction
Represents an ODF for Q-Ball imaging.
Definition: itkOrientationDistributionFunction.h:50
itk::OrientationDistributionFunction::GetDirection
static vnl_vector_fixed< double, 3 > GetDirection(int i)
itk::FiniteDiffOdfMaximaExtractionFilter::ThreadedGenerateData
void ThreadedGenerateData(const OutputImageRegionType &outputRegionForThread, ThreadIdType threadID)
Definition: itkFiniteDiffOdfMaximaExtractionFilter.cpp:326
itk::OrientationDistributionFunction::GetNeighbors
static std::vector< int > GetNeighbors(int idx)
itk::FiniteDiffOdfMaximaExtractionFilter::Cart2Sph
void Cart2Sph(const std::vector< DirectionType > &dir, vnl_matrix< double > &sphCoords)
Definition: itkFiniteDiffOdfMaximaExtractionFilter.cpp:457
max
static T max(T x, T y)
Definition: svm.cpp:70
itk::CompareVectors
static bool CompareVectors(const vnl_vector_fixed< double, 3 > &v1, const vnl_vector_fixed< double, 3 > &v2)
Definition: itkFiniteDiffOdfMaximaExtractionFilter.cpp:42
itk::OrientationDistributionFunction::GetGeneralizedFractionalAnisotropy
TComponent GetGeneralizedFractionalAnisotropy() const
mitk::New
static itkEventMacro(BoundingShapeInteractionEvent, itk::AnyEvent) class MITKBOUNDINGSHAPE_EXPORT BoundingShapeInteractor Pointer New()
Basic interaction methods for mitk::GeometryData.