Medical Imaging Interaction Toolkit  2016.11.0
Medical Imaging Interaction Toolkit
itkOrientationDistributionFunction.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 
17 #ifndef _itkOrientationDistributionFunction_txx
18 #define _itkOrientationDistributionFunction_txx
19 
20 #include <stdio.h>
21 #include <stdlib.h>
22 #include <time.h>
23 #include <algorithm>
24 #include <vector>
25 #include <vnl/algo/vnl_matrix_inverse.h>
26 #include "itkPointShell.h"
27 
28 #define _USE_MATH_DEFINES
29 #include <cmath>
30 #ifdef _MSC_VER
31  #if _MSC_VER <= 1700
32  #define fmin(a,b) ((a<=b)?(a):(b))
33  #define fmax(a,b) ((a>=b)?(a):(b))
34  #define isnan(c) (c!=c)
35  #endif
36 #endif
37 
38 #include <itkMatrix.h>
39 #include <vnl/vnl_matrix.h>
40 #include <vnl/vnl_matrix_fixed.h>
41 #include <vnl/vnl_inverse.h>
42 #include <typeinfo>
43 #include <ciso646>
44 
45 namespace itk
46 {
47 
48  template<class T, unsigned int N>
49  vtkPolyData* itk::OrientationDistributionFunction<T,N>::m_BaseMesh = nullptr;
50 
51  template<class T, unsigned int N>
52  double itk::OrientationDistributionFunction<T,N>::m_MaxChordLength = -1.0;
53 
54  template<class T, unsigned int N>
55  vnl_matrix_fixed<double, 3, N>* itk::OrientationDistributionFunction<T,N>::m_Directions
56  = itk::PointShell<N, vnl_matrix_fixed<double, 3, N> >::DistributePointShell();
57 
58  template<class T, unsigned int N>
59  std::vector< std::vector<int>* >* itk::OrientationDistributionFunction<T,N>::m_NeighborIdxs = nullptr;
60 
61  template<class T, unsigned int N>
62  std::vector< std::vector<int>* >* itk::OrientationDistributionFunction<T,N>::m_AngularRangeIdxs = nullptr;
63 
64  template<class T, unsigned int N>
65  std::vector<int>* itk::OrientationDistributionFunction<T,N>::m_HalfSphereIdxs = nullptr;
66 
67  template<class T, unsigned int N>
68  itk::SimpleFastMutexLock itk::OrientationDistributionFunction<T,N>::m_MutexBaseMesh;
69  template<class T, unsigned int N>
70  itk::SimpleFastMutexLock itk::OrientationDistributionFunction<T,N>::m_MutexHalfSphereIdxs;
71  template<class T, unsigned int N>
72  itk::SimpleFastMutexLock itk::OrientationDistributionFunction<T,N>::m_MutexNeighbors;
73  template<class T, unsigned int N>
74  itk::SimpleFastMutexLock itk::OrientationDistributionFunction<T,N>::m_MutexAngularRange;
75 
76 
77 #define ODF_PI M_PI
78 
79  /**
80  * Assignment Operator
81  */
82  template<class T, unsigned int NOdfDirections>
83  OrientationDistributionFunction<T, NOdfDirections>&
84  OrientationDistributionFunction<T, NOdfDirections>
85  ::operator= (const Self& r)
86  {
87  BaseArray::operator=(r);
88  return *this;
89  }
90 
91  /**
92  * Assignment Operator from a scalar constant
93  */
94  template<class T, unsigned int NOdfDirections>
95  OrientationDistributionFunction<T, NOdfDirections>&
96  OrientationDistributionFunction<T, NOdfDirections>
97  ::operator= (const ComponentType & r)
98  {
99  BaseArray::operator=(&r);
100  return *this;
101  }
102 
103  /**
104  * Assigment from a plain array
105  */
106  template<class T, unsigned int NOdfDirections>
107  OrientationDistributionFunction<T, NOdfDirections>&
108  OrientationDistributionFunction<T, NOdfDirections>
109  ::operator= (const ComponentArrayType r )
110  {
111  BaseArray::operator=(r);
112  return *this;
113  }
114 
115  /**
116  * Returns a temporary copy of a vector
117  */
118  template<class T, unsigned int NOdfDirections>
119  OrientationDistributionFunction<T, NOdfDirections>
120  OrientationDistributionFunction<T, NOdfDirections>
121  ::operator+(const Self & r) const
122  {
123  Self result;
124  for( unsigned int i=0; i<InternalDimension; i++)
125  {
126  result[i] = (*this)[i] + r[i];
127  }
128  return result;
129  }
130 
131  /**
132  * Returns a temporary copy of a vector
133  */
134  template<class T, unsigned int NOdfDirections>
135  OrientationDistributionFunction<T, NOdfDirections>
136  OrientationDistributionFunction<T, NOdfDirections>
137  ::operator-(const Self & r) const
138  {
139  Self result;
140  for( unsigned int i=0; i<InternalDimension; i++)
141  {
142  result[i] = (*this)[i] - r[i];
143  }
144  return result;
145  }
146 
147  /**
148  * Performs addition in place
149  */
150  template<class T, unsigned int NOdfDirections>
151  const OrientationDistributionFunction<T, NOdfDirections> &
152  OrientationDistributionFunction<T, NOdfDirections>
153  ::operator+=(const Self & r)
154  {
155  for( unsigned int i=0; i<InternalDimension; i++)
156  {
157  (*this)[i] += r[i];
158  }
159  return *this;
160  }
161 
162  /**
163  * Performs subtraction in place
164  */
165  template<class T, unsigned int NOdfDirections>
166  const OrientationDistributionFunction<T, NOdfDirections> &
167  OrientationDistributionFunction<T, NOdfDirections>
168  ::operator-=(const Self & r)
169  {
170  for( unsigned int i=0; i<InternalDimension; i++)
171  {
172  (*this)[i] -= r[i];
173  }
174  return *this;
175  }
176 
177  /**
178  * Performs multiplication by a scalar, in place
179  */
180  template<class T, unsigned int NOdfDirections>
181  const OrientationDistributionFunction<T, NOdfDirections> &
182  OrientationDistributionFunction<T, NOdfDirections>
183  ::operator*=(const RealValueType & r)
184  {
185  for( unsigned int i=0; i<InternalDimension; i++)
186  {
187  (*this)[i] *= r;
188  }
189  return *this;
190  }
191 
192  /**
193  * Performs division by a scalar, in place
194  */
195  template<class T, unsigned int NOdfDirections>
196  const OrientationDistributionFunction<T, NOdfDirections> &
197  OrientationDistributionFunction<T, NOdfDirections>
198  ::operator/=(const RealValueType & r)
199  {
200  for( unsigned int i=0; i<InternalDimension; i++)
201  {
202  (*this)[i] /= r;
203  }
204  return *this;
205  }
206 
207  /**
208  * Performs multiplication with a scalar
209  */
210  template<class T, unsigned int NOdfDirections>
211  OrientationDistributionFunction<T, NOdfDirections>
212  OrientationDistributionFunction<T, NOdfDirections>
213  ::operator*(const RealValueType & r) const
214  {
215  Self result;
216  for( unsigned int i=0; i<InternalDimension; i++)
217  {
218  result[i] = (*this)[i] * r;
219  }
220  return result;
221  }
222 
223  /**
224  * Performs division by a scalar
225  */
226  template<class T, unsigned int NOdfDirections>
227  OrientationDistributionFunction<T, NOdfDirections>
228  OrientationDistributionFunction<T, NOdfDirections>
229  ::operator/(const RealValueType & r) const
230  {
231  Self result;
232  for( unsigned int i=0; i<InternalDimension; i++)
233  {
234  result[i] = (*this)[i] / r;
235  }
236  return result;
237  }
238 
239  /**
240  * Matrix notation access to elements
241  */
242  template<class T, unsigned int NOdfDirections>
243  const typename OrientationDistributionFunction<T, NOdfDirections>::ValueType &
244  OrientationDistributionFunction<T, NOdfDirections>
245  ::operator()(unsigned int row, unsigned int col) const
246  {
247  unsigned int k;
248 
249  if( row < col )
250  {
251  k = row * InternalDimension + col - row * ( row + 1 ) / 2;
252  }
253  else
254  {
255  k = col * InternalDimension + row - col * ( col + 1 ) / 2;
256  }
257 
258 
259  if( k >= InternalDimension )
260  {
261  k = 0;
262  }
263 
264  return (*this)[k];
265  }
266 
267  /**
268  * Matrix notation access to elements
269  */
270  template<class T, unsigned int NOdfDirections>
271  typename OrientationDistributionFunction<T, NOdfDirections>::ValueType &
272  OrientationDistributionFunction<T, NOdfDirections>
273  ::operator()(unsigned int row, unsigned int col)
274  {
275  unsigned int k;
276 
277  if( row < col )
278  {
279  k = row * InternalDimension + col - row * ( row + 1 ) / 2;
280  }
281  else
282  {
283  k = col * InternalDimension + row - col * ( col + 1 ) / 2;
284  }
285 
286 
287  if( k >= InternalDimension )
288  {
289  k = 0;
290  }
291 
292  return (*this)[k];
293  }
294 
295  /**
296  * Set the Tensor to an Identity.
297  * Set ones in the diagonal and zeroes every where else.
298  */
299  template<class T, unsigned int NOdfDirections>
300  void
301  OrientationDistributionFunction<T, NOdfDirections>
302  ::SetIsotropic()
303  {
304  this->Fill(NumericTraits< T >::One / NOdfDirections);
305  }
306 
307  /**
308  * InitFromTensor()
309  */
310  template<class T, unsigned int NOdfDirections>
311  void
312  OrientationDistributionFunction<T, NOdfDirections>
313  ::InitFromTensor(itk::DiffusionTensor3D<T> tensor)
314  {
315  for(unsigned int i=0; i<NOdfDirections; i++)
316  {
317  /*
318  * | t0 t1 t2 | g0
319  * g0 g1 g2 * | t1 t3 t4 | * g1
320  * | t2 t4 t5 | g2
321  *
322  * = g0 * (t0g0*t1g1*t2g2)
323  * + g1 * (t1g0+t3g1+t4g2)
324  * + g2 * (t2g0+t4g1+t5g2)
325  */
326  T g0 = (*m_Directions)(0,i);
327  T g1 = (*m_Directions)(1,i);
328  T g2 = (*m_Directions)(2,i);
329  T t0 = tensor[0];
330  T t1 = tensor[1];
331  T t2 = tensor[2];
332  T t3 = tensor[3];
333  T t4 = tensor[4];
334  T t5 = tensor[5];
335  (*this)[i] = g0 * (t0*g0+t1*g1+t2*g2)
336  + g1 * (t1*g0+t3*g1+t4*g2)
337  + g2 * (t2*g0+t4*g1+t5*g2);
338 
339  if ((*this)[i]<0 || (*this)[i]!=(*this)[i])
340  (*this)[i] = 0;
341  }
342  }
343 
344  /**
345  * InitFromEllipsoid()
346  */
347  template<class T, unsigned int NOdfDirections>
348  void OrientationDistributionFunction<T, NOdfDirections>::
349  InitFromEllipsoid( itk::DiffusionTensor3D<T> tensor )
350  {
351  FixedArray<T, 6> nulltensor;
352  nulltensor.Fill(0.0);
353  if( tensor == nulltensor )
354  {
355  for ( unsigned int it=0; it < NOdfDirections; ++it ){ (*this)[it] = (T)0; }
356  MITK_DEBUG << "OrientationDistributionFunction<" << typeid(T).name() << ", " << NOdfDirections
357  << ">::InitFromEllipsoid(" << typeid(tensor).name()
358  << ") encountered a nulltensor as dti input point and ignorend it.";
359  return;
360  }
361 
362  typename itk::DiffusionTensor3D<T>::EigenValuesArrayType eigenValues;
363  typename itk::DiffusionTensor3D<T>::EigenVectorsMatrixType eigenVectors;
364  tensor.ComputeEigenAnalysis( eigenValues, eigenVectors ); // gives normalized eigenvectors as lines i.e. rows.
365  double a = eigenValues[0]; // those eigenvalues are the 3 |axes of the ellipsoid|,
366  double b = eigenValues[1]; // ComputeEigenAnalysis gives eigenValues in ascending < order per default,
367  double c = eigenValues[2]; // therefor the third eigenVector is the main direction of diffusion.
368 
369  if( a <= 0.0 || b <= 0.0 || c <= 0.0 )
370  {
371  for ( unsigned int it=0; it < NOdfDirections; ++it ){ (*this)[it] = (T)0; }
372  MITK_DEBUG << "OrientationDistributionFunction<" << typeid(T).name() << ", " << NOdfDirections
373  << ">::InitFromEllipsoid(" << typeid(tensor).name()
374  << ") encountered an eigenvalue <= 0 and ignored this input point.";
375  return;
376  }
377 
378  // check magnitude and scale towards 1 to minimize numerical condition kappa:
379 #ifdef _MSC_VER
380  #if _MSC_VER <= 1700
381  int exponent_a = floor(std::log(a)/std::log(2));
382  int exponent_b = floor(std::log(b)/std::log(2));
383  int exponent_c = floor(std::log(c)/std::log(2));
384  #else
385  int exponent_a = std::ilogb(a);
386  int exponent_b = std::ilogb(b);
387  int exponent_c = std::ilogb(c);
388  #endif
389 #else
390  int exponent_a = std::ilogb(a);
391  int exponent_b = std::ilogb(b);
392  int exponent_c = std::ilogb(c);
393 #endif
394 
395  T min_exponent= fmin(exponent_a, fmin(exponent_b, exponent_c) );
396  T max_exponent= fmax(exponent_c, fmax(exponent_b, exponent_a) );
397  int scale_exponent = floor(0.5 * (min_exponent + max_exponent));
398  double scaling = pow(2, scale_exponent);
399  a= a/scaling;
400  b= b/scaling;
401  c= c/scaling;
402 
403  vnl_matrix_fixed<double, 3, 3> eigenBase; // for change of base system.
404  for (int row = 0 ; row < 3; ++row) // Transposing because ComputeEigenAnalysis(,) gave us _row_ vectors.
405  {
406  for (int col = 0; col < 3; ++col)
407  {
408  eigenBase(row, col) = eigenVectors(col, row);
409  }
410  }
411  eigenBase= vnl_inverse(eigenBase); // Assuming canonical orthonormal system x=[1;0;0];y=[0;1;0];z=[0;0;1] for original DT.
412  eigenBase.assert_finite();
413 
414 #ifndef NDEBUG
415  double kappa=1.0;
416  { // calculate numerical condition kappa= ||f(x)-f(x~)|| approximately:
417  double gxaa = pow( a, -2.0); double gybb = pow( b, -2.0); double gzcc = pow( c, -2.0);
418  kappa = sqrt( pow( a, 2.0)+ pow( b, 2.0)+ pow( c, 2.0) + 1 ) / (gxaa + gybb + gzcc)
419  * sqrt( pow( a, -6.0)+ pow( b, -6.0)+ pow( c, -6.0) + pow( a, -4.0)+ pow( b, -4.0)+ pow( c, -4.0) );
420  MITK_DEBUG <<"kappa= "<< kappa << ", eigenvalues= [" << a <<", "<< b <<", "<< c <<"], eigenbase= ["
421  <<eigenBase(0,0)<<", "<<eigenBase(1,0)<<", "<<eigenBase(2,0)<<"; "
422  <<eigenBase(0,1)<<", "<<eigenBase(1,1)<<", "<<eigenBase(2,1)<<"; "
423  <<eigenBase(0,2)<<", "<<eigenBase(1,2)<<", "<<eigenBase(2,2)<<"] ";
424  if( std::isnan(kappa) )
425  {
426  MITK_DEBUG << "oh noes: kappa was NaN, setting kappa to 1e5"; // typical value kappa=1e5 for 1e-6<=a,b,c<=1e-4.
427  kappa=1e5;
428  };
429  }
430 #endif
431 
432  for( unsigned int i=0; i < NOdfDirections; ++i )
433  {
434  /// calculate probability p(g)=r as ellipsoids magnitude in direction of g' = B^-1 * g:
435  /// (g0*r/evx)² + (g1*r/evy)² + (g2*r/evz)² = 1, |g'|=1.
436 
437  vnl_vector_fixed<double, 3> g( (*m_Directions)(0,i), (*m_Directions)(1,i), (*m_Directions)(2,i) );
438  g = eigenBase*g; // passive change of base system of g.
439  g = g.normalize(); // unit vectors necessary.
440  (*this)[i] = scaling / sqrt( (g[0]/a)*(g[0]/a) + (g[1]/b)*(g[1]/b) + (g[2]/c)*(g[2]/c) );
441 
442 #ifndef NDEBUG
443  { // boundary check for numerical stability, assuming sigma=6, ||f(x~)-f~(x~)|| <= eps*kappa*sigma.
444  T min_ev= fmin(a, fmin(b, c) ); T max_ev= fmax(c, fmax(b, a) );
445  double eps= std::numeric_limits<T>::epsilon();
446  assert( scaling*min_ev <= ((*this)[i] + eps*kappa*6.0) ); // we should be between smallest and
447  assert( (*this)[i] <= (scaling*max_ev + eps*kappa*6.0) ); // biggest eigenvalue.
448  }
449 #endif
450  if ( (*this)[i] < T(0) || (*this)[i] > T(1) || std::isnan((*this)[i]) ) // P∈[0;1] sanity check.
451  { // C: NaN != NaN, C++11: isnan((*this)[i]).
452  MITK_DEBUG << "OrientationDistributionFunction<" << typeid(T).name() << ", " << NOdfDirections
453  << ">::InitFromEllipsoid(" << typeid(tensor).name()
454  << ") encountered a probability value out of range [0;1] and set it to zero: (*this)["
455  << i <<"]= " << (*this)[i];
456  (*this)[i] = T(0);
457  }
458  }
459  }
460 
461  /**
462  * L2-Normalization
463  */
464  template<class T, unsigned int NOdfDirections>
465  void
466  OrientationDistributionFunction<T, NOdfDirections>
467  ::L2Normalize()
468  {
469  T sum = 0;
470  for( unsigned int i=0; i<InternalDimension; i++)
471  {
472  sum += (*this)[i]*(*this)[i];
473  }
474  sum = std::sqrt(sum);
475  for( unsigned int i=0; i<InternalDimension; i++)
476  {
477  (*this)[i] = (*this)[i] / sum;
478  }
479  }
480 
481  /**
482  * Normalization to PDF
483  */
484  template<class T, unsigned int NOdfDirections>
485  void
486  OrientationDistributionFunction<T, NOdfDirections>
487  ::Normalize()
488  {
489  T sum = 0;
490  for( unsigned int i=0; i<InternalDimension; i++)
491  {
492  sum += (*this)[i];
493  }
494  if (sum>0)
495  {
496  for( unsigned int i=0; i<InternalDimension; i++)
497  {
498  (*this)[i] = (*this)[i] / sum;
499  }
500  }
501  }
502 
503  /**
504  * Min/Max-Normalization
505  */
506  template<class T, unsigned int NOdfDirections>
507  OrientationDistributionFunction<T, NOdfDirections>
508  OrientationDistributionFunction<T, NOdfDirections>
509  ::MinMaxNormalize() const
510  {
511  T max = NumericTraits<T>::NonpositiveMin();
512  T min = NumericTraits<T>::max();
513  for( unsigned int i=0; i<InternalDimension; i++)
514  {
515  max = (*this)[i] > max ? (*this)[i] : max;
516  min = (*this)[i] < min ? (*this)[i] : min;
517  }
518  Self retval;
519  for( unsigned int i=0; i<InternalDimension; i++)
520  {
521  retval[i] = ((*this)[i] - min) / (max - min);
522  }
523  return retval;
524  }
525 
526  /**
527  * Max-Normalization
528  */
529  template<class T, unsigned int NOdfDirections>
530  OrientationDistributionFunction<T, NOdfDirections>
531  OrientationDistributionFunction<T, NOdfDirections>
532  ::MaxNormalize() const
533  {
534  T max = NumericTraits<T>::NonpositiveMin();
535  for( unsigned int i=0; i<InternalDimension; i++)
536  {
537  max = (*this)[i] > max ? (*this)[i] : max;
538  }
539  Self retval;
540  for( unsigned int i=0; i<InternalDimension; i++)
541  {
542  retval[i] = (*this)[i] / max;
543  }
544  return retval;
545  }
546 
547  template<class T, unsigned int NOdfDirections>
548  T
549  OrientationDistributionFunction<T, NOdfDirections>
550  ::GetMaxValue() const
551  {
552  T max = NumericTraits<T>::NonpositiveMin();
553  for( unsigned int i=0; i<InternalDimension; i++)
554  {
555  if((*this)[i] >= max )
556  {
557  max = (*this)[i];
558  }
559  }
560  return max;
561  }
562 
563  template<class T, unsigned int NOdfDirections>
564  T
565  OrientationDistributionFunction<T, NOdfDirections>
566  ::GetMinValue() const
567  {
568  T min = NumericTraits<T>::max();
569  for( unsigned int i=0; i<InternalDimension; i++)
570  {
571  if((*this)[i] >= min )
572  {
573  min = (*this)[i];
574  }
575  }
576  return min;
577  }
578 
579  template<class T, unsigned int NOdfDirections>
580  T
581  OrientationDistributionFunction<T, NOdfDirections>
582  ::GetMeanValue() const
583  {
584  T sum = 0;
585  for( unsigned int i=0; i<InternalDimension; i++)
586  sum += (*this)[i];
587  return sum/InternalDimension;
588  }
589 
590  template<class T, unsigned int NOdfDirections>
591  double
592  OrientationDistributionFunction<T, NOdfDirections>
593  ::GetMaxChordLength()
594  {
595  if(m_MaxChordLength<0.0)
596  {
597  ComputeBaseMesh();
598  double max_dist = -1;
599  vtkPoints* points = m_BaseMesh->GetPoints();
600  for(int i=0; i<NOdfDirections; i++)
601  {
602  double p[3];
603  points->GetPoint(i,p);
604  std::vector<int> neighbors = GetNeighbors(i);
605  for(std::size_t j=0; j<neighbors.size(); j++)
606  {
607  double n[3];
608  points->GetPoint(neighbors[j],n);
609  double d = sqrt(
610  (p[0]-n[0])*(p[0]-n[0]) +
611  (p[1]-n[1])*(p[1]-n[1]) +
612  (p[2]-n[2])*(p[2]-n[2]));
613  max_dist = d>max_dist ? d : max_dist;
614  }
615  }
616  m_MaxChordLength = max_dist;
617  }
618  return m_MaxChordLength;
619  }
620 
621  template<class T, unsigned int NOdfDirections>
622  void
623  OrientationDistributionFunction<T, NOdfDirections>
624  ::ComputeBaseMesh()
625  {
626 
627  m_MutexBaseMesh.Lock();
628  if(m_BaseMesh == nullptr)
629  {
630 
631  vtkPoints* points = vtkPoints::New();
632  for(unsigned int j=0; j<NOdfDirections; j++){
633  double x = (*m_Directions)(0,j);
634  double y = (*m_Directions)(1,j);
635  double z = (*m_Directions)(2,j);
636  double az = atan2(y,x);
637  double elev = atan2(z,sqrt(x*x+y*y));
638  double r = sqrt(x*x+y*y+z*z);
639  points->InsertNextPoint(az,elev,r);
640  }
641 
642  vtkPolyData* polydata = vtkPolyData::New();
643  polydata->SetPoints( points );
644  vtkDelaunay2D *delaunay = vtkDelaunay2D::New();
645  delaunay->SetInputData( polydata );
646  delaunay->Update();
647 
648  vtkCellArray* vtkpolys = delaunay->GetOutput()->GetPolys();
649  vtkCellArray* vtknewpolys = vtkCellArray::New();
650  vtkIdType npts; vtkIdType *pts;
651  while(vtkpolys->GetNextCell(npts,pts))
652  {
653  bool insert = true;
654  for(int i=0; i<npts; i++)
655  {
656  double *tmpPoint = points->GetPoint(pts[i]);
657  double az = tmpPoint[0];
658  double elev = tmpPoint[1];
659  if((std::abs(az)>ODF_PI-0.5) || (std::abs(elev)>ODF_PI/2-0.5))
660  insert = false;
661  }
662  if(insert)
663  vtknewpolys->InsertNextCell(npts, pts);
664  }
665 
666  vtkPoints* points2 = vtkPoints::New();
667  for(unsigned int j=0; j<NOdfDirections; j++){
668  double x = -(*m_Directions)(0,j);
669  double y = -(*m_Directions)(2,j);
670  double z = -(*m_Directions)(1,j);
671  double az = atan2(y,x);
672  double elev = atan2(z,sqrt(x*x+y*y));
673  double r = sqrt(x*x+y*y+z*z);
674  points2->InsertNextPoint(az,elev,r);
675  }
676 
677  vtkPolyData* polydata2 = vtkPolyData::New();
678  polydata2->SetPoints( points2 );
679  vtkDelaunay2D *delaunay2 = vtkDelaunay2D::New();
680  delaunay2->SetInputData( polydata2 );
681  delaunay2->Update();
682 
683  vtkpolys = delaunay2->GetOutput()->GetPolys();
684  while(vtkpolys->GetNextCell(npts,pts))
685  {
686  bool insert = true;
687  for(int i=0; i<npts; i++)
688  {
689  double *tmpPoint = points2->GetPoint(pts[i]);
690  double az = tmpPoint[0];
691  double elev = tmpPoint[1];
692  if((std::abs(az)>ODF_PI-0.5) || (std::abs(elev)>ODF_PI/2-0.5))
693  insert = false;
694  }
695  if(insert)
696  vtknewpolys->InsertNextCell(npts, pts);
697  }
698 
699  polydata->SetPolys(vtknewpolys);
700 
701  for (unsigned int p = 0; p < NOdfDirections; p++)
702  {
703  points->SetPoint(p,m_Directions->get_column(p).data_block());
704  }
705  polydata->SetPoints( points );
706 
707  m_BaseMesh = polydata;
708  }
709  m_MutexBaseMesh.Unlock();
710  }
711 
712  /**
713  * Extract the principle diffusion direction
714  */
715  template<class T, unsigned int NOdfDirections>
716  int
717  OrientationDistributionFunction<T, NOdfDirections>
718  ::GetPrincipleDiffusionDirection() const
719  {
720  T max = NumericTraits<T>::NonpositiveMin();
721  int maxidx = -1;
722  for( unsigned int i=0; i<InternalDimension; i++)
723  {
724  if((*this)[i] >= max )
725  {
726  max = (*this)[i];
727  maxidx = i;
728  }
729  }
730 
731  return maxidx;
732  }
733 
734  template<class T, unsigned int NOdfDirections>
735  std::vector<int>
736  OrientationDistributionFunction<T, NOdfDirections>
737  ::GetNeighbors(int idx)
738  {
739  ComputeBaseMesh();
740 
741  m_MutexNeighbors.Lock();
742  if(m_NeighborIdxs == nullptr)
743  {
744  m_NeighborIdxs = new std::vector< std::vector<int>* >();
745  vtkCellArray* polys = m_BaseMesh->GetPolys();
746 
747  for(unsigned int i=0; i<NOdfDirections; i++)
748  {
749  auto idxs = new std::vector<int>();
750  polys->InitTraversal();
751  vtkIdType npts; vtkIdType *pts;
752  while(polys->GetNextCell(npts,pts))
753  {
754  if( pts[0] == i )
755  {
756  idxs->push_back(pts[1]);
757  idxs->push_back(pts[2]);
758  }
759  else if( pts[1] == i )
760  {
761  idxs->push_back(pts[0]);
762  idxs->push_back(pts[2]);
763  }
764  else if( pts[2] == i )
765  {
766  idxs->push_back(pts[0]);
767  idxs->push_back(pts[1]);
768  }
769  }
770  std::sort(idxs->begin(), idxs->end());
771  std::vector< int >::iterator endLocation;
772  endLocation = std::unique( idxs->begin(), idxs->end() );
773  idxs->erase(endLocation, idxs->end());
774  m_NeighborIdxs->push_back(idxs);
775  }
776  }
777  m_MutexNeighbors.Unlock();
778 
779  return *m_NeighborIdxs->at(idx);
780  }
781 
782  /**
783  * Extract the n-th diffusion direction
784  */
785  template<class T, unsigned int NOdfDirections>
786  int
787  OrientationDistributionFunction<T, NOdfDirections>
788  ::GetNthDiffusionDirection(int n, vnl_vector_fixed<double,3> rndVec) const
789  {
790 
791  if( n == 0 )
792  return GetPrincipleDiffusionDirection();
793 
794  m_MutexHalfSphereIdxs.Lock();
795  if( !m_HalfSphereIdxs )
796  {
797  m_HalfSphereIdxs = new std::vector<int>();
798  for( unsigned int i=0; i<InternalDimension; i++)
799  {
800  if(dot_product(m_Directions->get_column(i),rndVec) > 0.0)
801  {
802  m_HalfSphereIdxs->push_back(i);
803  }
804  }
805  }
806  m_MutexHalfSphereIdxs.Unlock();
807 
808  // collect indices of directions
809  // that are local maxima
810  std::vector<int> localMaxima;
811  std::vector<int>::iterator it;
812  for( it=m_HalfSphereIdxs->begin();
813  it!=m_HalfSphereIdxs->end();
814  it++)
815  {
816  std::vector<int> nbs = GetNeighbors(*it);
817  std::vector<int>::iterator it2;
818  bool max = true;
819  for(it2 = nbs.begin();
820  it2 != nbs.end();
821  it2++)
822  {
823  if((*this)[*it2] > (*this)[*it])
824  {
825  max = false;
826  break;
827  }
828  }
829  if(max)
830  localMaxima.push_back(*it);
831  }
832 
833  // delete n highest local maxima from list
834  // and return remaining highest
835  int maxidx = -1;
836  std::vector<int>::iterator itMax;
837  for( int i=0; i<=n; i++ )
838  {
839  maxidx = -1;
840  T max = NumericTraits<T>::NonpositiveMin();
841  for(it = localMaxima.begin();
842  it != localMaxima.end();
843  it++)
844  {
845  if((*this)[*it]>max)
846  {
847  max = (*this)[*it];
848  maxidx = *it;
849  }
850  }
851  it = find(localMaxima.begin(), localMaxima.end(), maxidx);
852  if(it!=localMaxima.end())
853  localMaxima.erase(it);
854  }
855 
856  return maxidx;
857  }
858 
859  template < typename TComponent, unsigned int NOdfDirections >
860  vnl_vector_fixed<double,3> itk::OrientationDistributionFunction<TComponent, NOdfDirections>
861  ::GetDirection( int i )
862  {
863  return m_Directions->get_column(i);
864  }
865 
866  /**
867  * Interpolate a position between sampled directions
868  */
869  template<class T, unsigned int NOdfDirections>
870  T
871  OrientationDistributionFunction<T, NOdfDirections>
872  ::GetInterpolatedComponent(vnl_vector_fixed<double,3> dir, InterpolationMethods method) const
873  {
874 
875  ComputeBaseMesh();
876  double retval = -1.0;
877 
878  switch(method)
879  {
880  case ODF_NEAREST_NEIGHBOR_INTERP:
881  {
882 
883  vtkPoints* points = m_BaseMesh->GetPoints();
884  double current_min = NumericTraits<double>::max();
885  int current_min_idx = -1;
886  for(int i=0; i<NOdfDirections; i++)
887  {
888  vnl_vector_fixed<double,3> P(points->GetPoint(i));
889  double dist = (dir-P).two_norm();
890  current_min_idx = dist<current_min ? i : current_min_idx;
891  current_min = dist<current_min ? dist : current_min;
892  }
893  retval = this->GetNthComponent(current_min_idx);
894  break;
895 
896  }
897  case ODF_TRILINEAR_BARYCENTRIC_INTERP:
898  {
899 
900  double maxChordLength = GetMaxChordLength();
901  vtkCellArray* polys = m_BaseMesh->GetPolys();
902  vtkPoints* points = m_BaseMesh->GetPoints();
903  vtkIdType npts; vtkIdType *pts;
904  double current_min = NumericTraits<double>::max();
905  polys->InitTraversal();
906  while(polys->GetNextCell(npts,pts))
907  {
908  vnl_vector_fixed<double,3> A(points->GetPoint(pts[0]));
909  vnl_vector_fixed<double,3> B(points->GetPoint(pts[1]));
910  vnl_vector_fixed<double,3> C(points->GetPoint(pts[2]));
911 
912  vnl_vector_fixed<double,3> d1;
913  d1.put(0,(dir-A).two_norm());
914  d1.put(1,(dir-B).two_norm());
915  d1.put(2,(dir-C).two_norm());
916  double maxval = d1.max_value();
917  if(maxval>maxChordLength)
918  {
919  continue;
920  }
921 
922  // Compute vectors
923  vnl_vector_fixed<double,3> v0 = C - A;
924  vnl_vector_fixed<double,3> v1 = B - A;
925 
926  // Project direction to plane ABC
927  vnl_vector_fixed<double,3> v6 = dir;
928  vnl_vector_fixed<double,3> cross = vnl_cross_3d(v0, v1);
929  cross = cross.normalize();
930  vtkPlane::ProjectPoint(v6.data_block(),A.data_block(),cross.data_block(),v6.data_block());
931  v6 = v6-A;
932 
933  // Calculate barycentric coords
934  vnl_matrix_fixed<double,3,2> mat;
935  mat.set_column(0, v0);
936  mat.set_column(1, v1);
937  vnl_matrix_inverse<double> inv(mat);
938  vnl_matrix_fixed<double,2,3> inver = inv.pinverse();
939  vnl_vector<double> uv = inv.pinverse()*v6;
940 
941  // Check if point is in triangle
942  double eps = 0.01;
943  if( (uv(0) >= 0-eps) && (uv(1) >= 0-eps) && (uv(0) + uv(1) <= 1+eps) )
944  {
945  // check if minimum angle is the max so far
946  if(d1.two_norm() < current_min)
947  {
948  current_min = d1.two_norm();
949  vnl_vector<double> barycentricCoords(3);
950  barycentricCoords[2] = uv[0]<0 ? 0 : (uv[0]>1?1:uv[0]);
951  barycentricCoords[1] = uv[1]<0 ? 0 : (uv[1]>1?1:uv[1]);
952  barycentricCoords[0] = 1-(barycentricCoords[1]+barycentricCoords[2]);
953  retval = barycentricCoords[0]*this->GetNthComponent(pts[0]) +
954  barycentricCoords[1]*this->GetNthComponent(pts[1]) +
955  barycentricCoords[2]*this->GetNthComponent(pts[2]);
956  }
957  }
958  }
959 
960  break;
961  }
962  case ODF_SPHERICAL_GAUSSIAN_BASIS_FUNCTIONS:
963  {
964  double maxChordLength = GetMaxChordLength();
965  double sigma = asin(maxChordLength/2);
966 
967  // this is the contribution of each kernel to each sampling point on the
968  // equator
969  vnl_vector<double> contrib;
970  contrib.set_size(NOdfDirections);
971 
972  vtkPoints* points = m_BaseMesh->GetPoints();
973  double sum = 0;
974  for(int i=0; i<NOdfDirections; i++)
975  {
976  vnl_vector_fixed<double,3> P(points->GetPoint(i));
977  double stv = dir[0]*P[0]
978  + dir[1]*P[1]
979  + dir[2]*P[2];
980  stv = (stv<-1.0) ? -1.0 : ( (stv>1.0) ? 1.0 : stv);
981  double x = acos(stv);
982  contrib[i] = (1.0/(sigma*sqrt(2.0*ODF_PI)))
983  *exp((-x*x)/(2*sigma*sigma));
984  sum += contrib[i];
985  }
986 
987  retval = 0;
988  for(int i=0; i<NOdfDirections; i++)
989  {
990  retval += (contrib[i] / sum)*this->GetNthComponent(i);
991  }
992  break;
993  }
994 
995  }
996 
997  if(retval==-1)
998  {
999  std::cout << "Interpolation failed" << std::endl;
1000  return 0;
1001  }
1002 
1003  return retval;
1004 
1005  }
1006 
1007  /**
1008  * Calculate Generalized Fractional Anisotropy
1009  */
1010  template<class T, unsigned int NOdfDirections>
1011  T
1012  OrientationDistributionFunction<T, NOdfDirections>
1013  ::GetGeneralizedFractionalAnisotropy() const
1014  {
1015  double mean = 0;
1016  double std = 0;
1017  double rms = 0;
1018  for( unsigned int i=0; i<InternalDimension; i++)
1019  {
1020  T val = (*this)[i];
1021  mean += val;
1022  }
1023  mean /= NOdfDirections;
1024  for( unsigned int i=0; i<InternalDimension; i++)
1025  {
1026  T val = (*this)[i];
1027  std += (val - mean) * (val - mean);
1028  rms += val*val;
1029  }
1030  std *= NOdfDirections;
1031  rms *= NOdfDirections - 1;
1032 
1033  if(rms == 0)
1034  {
1035  return 0;
1036  }
1037  else
1038  {
1039  return sqrt(std/rms);
1040  }
1041  }
1042 
1043 
1044  template < typename T, unsigned int N>
1045  T itk::OrientationDistributionFunction<T, N>
1046  ::GetGeneralizedGFA( int k, int p ) const
1047  {
1048  double mean = 0;
1049  double std = 0;
1050  double rms = 0;
1051  double max = NumericTraits<double>::NonpositiveMin();
1052  for( unsigned int i=0; i<InternalDimension; i++)
1053  {
1054  double val = (double)(*this)[i];
1055  mean += pow(val,(double)p);
1056  max = val > max ? val : max;
1057  }
1058  max = pow(max,(double)p);
1059  mean /= N;
1060  for( unsigned int i=0; i<InternalDimension; i++)
1061  {
1062  double val = (double)(*this)[i];
1063  std += (pow(val,(double)p) - mean) * (pow(val,(double)p) - mean);
1064  if(k>0)
1065  {
1066  rms += pow(val,(double)(p*k));
1067  }
1068  }
1069  std /= N - 1;
1070  std = sqrt(std);
1071 
1072  if(k>0)
1073  {
1074  rms /= N;
1075  rms = pow(rms,(double)(1.0/k));
1076  }
1077  else if(k<0) // lim k->inf gives us the maximum
1078  {
1079  rms = max;
1080  }
1081  else // k==0 undefined, we define zeros root from 1 as 1
1082  {
1083  rms = 1;
1084  }
1085 
1086  if(rms == 0)
1087  {
1088  return 0;
1089  }
1090  else
1091  {
1092  return (T)(std/rms);
1093  }
1094  }
1095 
1096  /**
1097  * Calculate Nematic Order Parameter
1098  */
1099  template < typename T, unsigned int N >
1100  T itk::OrientationDistributionFunction<T, N>
1101  ::GetNematicOrderParameter() const
1102  {
1103  // not yet implemented
1104  return 0;
1105  }
1106 
1107  /**
1108  * Calculate StdDev by MaxValue
1109  */
1110  template < typename T, unsigned int N >
1111  T itk::OrientationDistributionFunction<T, N>
1112  ::GetStdDevByMaxValue() const
1113  {
1114  double mean = 0;
1115  double std = 0;
1116  T max = NumericTraits<T>::NonpositiveMin();
1117  for( unsigned int i=0; i<InternalDimension; i++)
1118  {
1119  T val = (*this)[i];
1120  mean += val;
1121  max = (*this)[i] > max ? (*this)[i] : max;
1122  }
1123  mean /= InternalDimension;
1124  for( unsigned int i=0; i<InternalDimension; i++)
1125  {
1126  T val = (*this)[i];
1127  std += (val - mean) * (val - mean);
1128  }
1129  std /= InternalDimension-1;
1130 
1131  if(max == 0)
1132  {
1133  return 0;
1134  }
1135  else
1136  {
1137  return (sqrt(std)/max);
1138  }
1139  }
1140 
1141  template < typename T, unsigned int N >
1142  T itk::OrientationDistributionFunction<T, N>
1143  ::GetPrincipleCurvature(double alphaMinDegree, double alphaMaxDegree, int invert) const
1144  {
1145  // following loop only performed once
1146  // (computing indices of each angular range)
1147  m_MutexAngularRange.Lock();
1148  if(m_AngularRangeIdxs == nullptr)
1149  {
1150  m_AngularRangeIdxs = new std::vector< std::vector<int>* >();
1151  for(unsigned int i=0; i<N; i++)
1152  {
1153  vnl_vector_fixed<double,3> pDir = GetDirection(i);
1154  auto idxs = new std::vector<int>();
1155  for(unsigned int j=0; j<N; j++)
1156  {
1157  vnl_vector_fixed<double,3> cDir = GetDirection(j);
1158  double angle = ( 180 / ODF_PI ) * acos( dot_product(pDir, cDir) );
1159  if( (angle < alphaMaxDegree) && (angle > alphaMinDegree) )
1160  {
1161  idxs->push_back(j);
1162  }
1163  }
1164  m_AngularRangeIdxs->push_back(idxs);
1165  }
1166  }
1167  m_MutexAngularRange.Unlock();
1168 
1169  // find the maximum (or minimum) direction (remember index and value)
1170  T mode;
1171  int pIdx = -1;
1172  if(invert == 0)
1173  {
1174  pIdx = GetPrincipleDiffusionDirection();
1175  mode = (*this)[pIdx];
1176  }
1177  else
1178  {
1179  mode = NumericTraits<T>::max();
1180  for( unsigned int i=0; i<N; i++)
1181  {
1182  if((*this)[i] < mode )
1183  {
1184  mode = (*this)[i];
1185  pIdx = i;
1186  }
1187  }
1188  }
1189  //////////////////////
1190  //////// compute median of mode and its neighbors to become more stable to noise
1191  //////// compared to simply using the mode
1192  //////////////////////
1193 
1194  //////// values of mode and its neighbors
1195  //////std::vector<int> nbs = GetNeighbors(pIdx);
1196  //////std::vector<T> modeAndNeighborVals;
1197  //////modeAndNeighborVals.push_back(mode);
1198  //////int numNeighbors = nbs.size();
1199  //////for(int i=0; i<numNeighbors; i++)
1200  //////{
1201  ////// modeAndNeighborVals.push_back((*this)[nbs[i]]);
1202  //////}
1203 
1204  //////// sort by value
1205  //////std::sort( modeAndNeighborVals.begin(), modeAndNeighborVals.end() );
1206 
1207  //////// median of mode and neighbors
1208  //////mode = modeAndNeighborVals[floor(0.5*(double)(numNeighbors+1)+0.5)];
1209 
1210  ////////////////
1211  // computing a quantile of the angular range
1212  ////////////////
1213 
1214  // define quantile
1215  double quantile = 0.00;
1216 
1217  // collect all values in angular range of mode
1218  std::vector<T> odfValuesInAngularRange;
1219  int numInRange = m_AngularRangeIdxs->at(pIdx)->size();
1220  for(int i=0; i<numInRange; i++)
1221  {
1222  odfValuesInAngularRange.push_back((*this)[(*m_AngularRangeIdxs->at(pIdx))[i] ]);
1223  }
1224 
1225  // sort them by value
1226  std::sort( odfValuesInAngularRange.begin(), odfValuesInAngularRange.end() );
1227 
1228  // median of angular range
1229  T median = odfValuesInAngularRange[floor(quantile*(double)numInRange+0.5)];
1230 
1231  // compute and return final value
1232  if(mode > median)
1233  {
1234  return mode/median - 1.0;
1235  }
1236  else
1237  {
1238  return median/mode - 1.0;
1239  }
1240  }
1241 
1242  /**
1243  * Calculate Normalized Entropy
1244  */
1245  template < typename T, unsigned int N >
1246  T itk::OrientationDistributionFunction<T, N>
1247  ::GetNormalizedEntropy() const
1248  {
1249  double mean = 0;
1250  for( unsigned int i=0; i<InternalDimension; i++)
1251  {
1252  T val = (*this)[i];
1253  if( val != 0 )
1254  {
1255  val = log(val);
1256  }
1257  else
1258  {
1259  val = log(0.0000001);
1260  }
1261  mean += val;
1262  }
1263  double _n = (double) InternalDimension;
1264  mean /= _n;
1265  return (T) (-_n / log(_n) * mean);
1266  }
1267 
1268  /**
1269  * Pre-multiply the Tensor by a Matrix
1270  */
1271  template<class T, unsigned int NOdfDirections>
1272  OrientationDistributionFunction<T, NOdfDirections>
1273  OrientationDistributionFunction<T, NOdfDirections>
1274  ::PreMultiply( const MatrixType & m ) const
1275  {
1276  Self result;
1277  typedef typename NumericTraits<T>::AccumulateType AccumulateType;
1278  for(unsigned int r=0; r<NOdfDirections; r++)
1279  {
1280  for(unsigned int c=0; c<NOdfDirections; c++)
1281  {
1282  AccumulateType sum = NumericTraits<AccumulateType>::ZeroValue();
1283  for(unsigned int t=0; t<NOdfDirections; t++)
1284  {
1285  sum += m(r,t) * (*this)(t,c);
1286  }
1287  result(r,c) = static_cast<T>( sum );
1288  }
1289  }
1290  return result;
1291  }
1292 
1293  /**
1294  * Post-multiply the Tensor by a Matrix
1295  */
1296  template<class T, unsigned int NOdfDirections>
1297  OrientationDistributionFunction<T, NOdfDirections>
1298  OrientationDistributionFunction<T, NOdfDirections>
1299  ::PostMultiply( const MatrixType & m ) const
1300  {
1301  Self result;
1302  typedef typename NumericTraits<T>::AccumulateType AccumulateType;
1303  for(unsigned int r=0; r<NOdfDirections; r++)
1304  {
1305  for(unsigned int c=0; c<NOdfDirections; c++)
1306  {
1307  AccumulateType sum = NumericTraits<AccumulateType>::ZeroValue();
1308  for(unsigned int t=0; t<NOdfDirections; t++)
1309  {
1310  sum += (*this)(r,t) * m(t,c);
1311  }
1312  result(r,c) = static_cast<T>( sum );
1313  }
1314  }
1315  return result;
1316  }
1317 
1318  /**
1319  * Print content to an ostream
1320  */
1321  template<class T, unsigned int NOdfDirections>
1322  std::ostream &
1323  operator<<(std::ostream& os,const OrientationDistributionFunction<T, NOdfDirections> & c )
1324  {
1325  for(unsigned int i=0; i<c.GetNumberOfComponents(); i++)
1326  {
1327  os << static_cast<typename NumericTraits<T>::PrintType>(c[i]) << " ";
1328  }
1329  return os;
1330  }
1331 
1332  /**
1333  * Read content from an istream
1334  */
1335  template<class T, unsigned int NOdfDirections>
1336  std::istream &
1337  operator>>(std::istream& is, OrientationDistributionFunction<T, NOdfDirections> & dt )
1338  {
1339  for(unsigned int i=0; i < dt.GetNumberOfComponents(); i++)
1340  {
1341  is >> dt[i];
1342  }
1343  return is;
1344  }
1345 
1346 } // end namespace itk
1347 
1348 #endif