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
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.