1 /*===================================================================
3 The Medical Imaging Interaction Toolkit (MITK)
5 Copyright (c) German Cancer Research Center,
6 Division of Medical and Biological Informatics.
9 This software is distributed WITHOUT ANY WARRANTY; without
10 even the implied warranty of MERCHANTABILITY or FITNESS FOR
13 See LICENSE.txt or http://www.mitk.org for details.
15 ===================================================================*/
17 #ifndef _itkOrientationDistributionFunction_txx
18 #define _itkOrientationDistributionFunction_txx
25 #include <vnl/algo/vnl_matrix_inverse.h>
26 #include "itkPointShell.h"
28 #define _USE_MATH_DEFINES
32 #define fmin(a,b) ((a<=b)?(a):(b))
33 #define fmax(a,b) ((a>=b)?(a):(b))
34 #define isnan(c) (c!=c)
38 #include <itkMatrix.h>
39 #include <vnl/vnl_matrix.h>
40 #include <vnl/vnl_matrix_fixed.h>
41 #include <vnl/vnl_inverse.h>
48 template<class T, unsigned int N>
49 vtkPolyData* itk::OrientationDistributionFunction<T,N>::m_BaseMesh = nullptr;
51 template<class T, unsigned int N>
52 double itk::OrientationDistributionFunction<T,N>::m_MaxChordLength = -1.0;
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();
58 template<class T, unsigned int N>
59 std::vector< std::vector<int>* >* itk::OrientationDistributionFunction<T,N>::m_NeighborIdxs = nullptr;
61 template<class T, unsigned int N>
62 std::vector< std::vector<int>* >* itk::OrientationDistributionFunction<T,N>::m_AngularRangeIdxs = nullptr;
64 template<class T, unsigned int N>
65 std::vector<int>* itk::OrientationDistributionFunction<T,N>::m_HalfSphereIdxs = nullptr;
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;
82 template<class T, unsigned int NOdfDirections>
83 OrientationDistributionFunction<T, NOdfDirections>&
84 OrientationDistributionFunction<T, NOdfDirections>
85 ::operator= (const Self& r)
87 BaseArray::operator=(r);
92 * Assignment Operator from a scalar constant
94 template<class T, unsigned int NOdfDirections>
95 OrientationDistributionFunction<T, NOdfDirections>&
96 OrientationDistributionFunction<T, NOdfDirections>
97 ::operator= (const ComponentType & r)
99 BaseArray::operator=(&r);
104 * Assigment from a plain array
106 template<class T, unsigned int NOdfDirections>
107 OrientationDistributionFunction<T, NOdfDirections>&
108 OrientationDistributionFunction<T, NOdfDirections>
109 ::operator= (const ComponentArrayType r )
111 BaseArray::operator=(r);
116 * Returns a temporary copy of a vector
118 template<class T, unsigned int NOdfDirections>
119 OrientationDistributionFunction<T, NOdfDirections>
120 OrientationDistributionFunction<T, NOdfDirections>
121 ::operator+(const Self & r) const
124 for( unsigned int i=0; i<InternalDimension; i++)
126 result[i] = (*this)[i] + r[i];
132 * Returns a temporary copy of a vector
134 template<class T, unsigned int NOdfDirections>
135 OrientationDistributionFunction<T, NOdfDirections>
136 OrientationDistributionFunction<T, NOdfDirections>
137 ::operator-(const Self & r) const
140 for( unsigned int i=0; i<InternalDimension; i++)
142 result[i] = (*this)[i] - r[i];
148 * Performs addition in place
150 template<class T, unsigned int NOdfDirections>
151 const OrientationDistributionFunction<T, NOdfDirections> &
152 OrientationDistributionFunction<T, NOdfDirections>
153 ::operator+=(const Self & r)
155 for( unsigned int i=0; i<InternalDimension; i++)
163 * Performs subtraction in place
165 template<class T, unsigned int NOdfDirections>
166 const OrientationDistributionFunction<T, NOdfDirections> &
167 OrientationDistributionFunction<T, NOdfDirections>
168 ::operator-=(const Self & r)
170 for( unsigned int i=0; i<InternalDimension; i++)
178 * Performs multiplication by a scalar, in place
180 template<class T, unsigned int NOdfDirections>
181 const OrientationDistributionFunction<T, NOdfDirections> &
182 OrientationDistributionFunction<T, NOdfDirections>
183 ::operator*=(const RealValueType & r)
185 for( unsigned int i=0; i<InternalDimension; i++)
193 * Performs division by a scalar, in place
195 template<class T, unsigned int NOdfDirections>
196 const OrientationDistributionFunction<T, NOdfDirections> &
197 OrientationDistributionFunction<T, NOdfDirections>
198 ::operator/=(const RealValueType & r)
200 for( unsigned int i=0; i<InternalDimension; i++)
208 * Performs multiplication with a scalar
210 template<class T, unsigned int NOdfDirections>
211 OrientationDistributionFunction<T, NOdfDirections>
212 OrientationDistributionFunction<T, NOdfDirections>
213 ::operator*(const RealValueType & r) const
216 for( unsigned int i=0; i<InternalDimension; i++)
218 result[i] = (*this)[i] * r;
224 * Performs division by a scalar
226 template<class T, unsigned int NOdfDirections>
227 OrientationDistributionFunction<T, NOdfDirections>
228 OrientationDistributionFunction<T, NOdfDirections>
229 ::operator/(const RealValueType & r) const
232 for( unsigned int i=0; i<InternalDimension; i++)
234 result[i] = (*this)[i] / r;
240 * Matrix notation access to elements
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
251 k = row * InternalDimension + col - row * ( row + 1 ) / 2;
255 k = col * InternalDimension + row - col * ( col + 1 ) / 2;
259 if( k >= InternalDimension )
268 * Matrix notation access to elements
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)
279 k = row * InternalDimension + col - row * ( row + 1 ) / 2;
283 k = col * InternalDimension + row - col * ( col + 1 ) / 2;
287 if( k >= InternalDimension )
296 * Set the Tensor to an Identity.
297 * Set ones in the diagonal and zeroes every where else.
299 template<class T, unsigned int NOdfDirections>
301 OrientationDistributionFunction<T, NOdfDirections>
304 this->Fill(NumericTraits< T >::One / NOdfDirections);
310 template<class T, unsigned int NOdfDirections>
312 OrientationDistributionFunction<T, NOdfDirections>
313 ::InitFromTensor(itk::DiffusionTensor3D<T> tensor)
315 for(unsigned int i=0; i<NOdfDirections; i++)
319 * g0 g1 g2 * | t1 t3 t4 | * g1
322 * = g0 * (t0g0*t1g1*t2g2)
323 * + g1 * (t1g0+t3g1+t4g2)
324 * + g2 * (t2g0+t4g1+t5g2)
326 T g0 = (*m_Directions)(0,i);
327 T g1 = (*m_Directions)(1,i);
328 T g2 = (*m_Directions)(2,i);
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);
339 if ((*this)[i]<0 || (*this)[i]!=(*this)[i])
345 * InitFromEllipsoid()
347 template<class T, unsigned int NOdfDirections>
348 void OrientationDistributionFunction<T, NOdfDirections>::
349 InitFromEllipsoid( itk::DiffusionTensor3D<T> tensor )
351 FixedArray<T, 6> nulltensor;
352 nulltensor.Fill(0.0);
353 if( tensor == nulltensor )
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.";
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.
369 if( a <= 0.0 || b <= 0.0 || c <= 0.0 )
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.";
378 // check magnitude and scale towards 1 to minimize numerical condition kappa:
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));
385 int exponent_a = std::ilogb(a);
386 int exponent_b = std::ilogb(b);
387 int exponent_c = std::ilogb(c);
390 int exponent_a = std::ilogb(a);
391 int exponent_b = std::ilogb(b);
392 int exponent_c = std::ilogb(c);
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);
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.
406 for (int col = 0; col < 3; ++col)
408 eigenBase(row, col) = eigenVectors(col, row);
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();
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) )
426 MITK_DEBUG << "oh noes: kappa was NaN, setting kappa to 1e5"; // typical value kappa=1e5 for 1e-6<=a,b,c<=1e-4.
432 for( unsigned int i=0; i < NOdfDirections; ++i )
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.
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) );
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.
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];
464 template<class T, unsigned int NOdfDirections>
466 OrientationDistributionFunction<T, NOdfDirections>
470 for( unsigned int i=0; i<InternalDimension; i++)
472 sum += (*this)[i]*(*this)[i];
474 sum = std::sqrt(sum);
475 for( unsigned int i=0; i<InternalDimension; i++)
477 (*this)[i] = (*this)[i] / sum;
482 * Normalization to PDF
484 template<class T, unsigned int NOdfDirections>
486 OrientationDistributionFunction<T, NOdfDirections>
490 for( unsigned int i=0; i<InternalDimension; i++)
496 for( unsigned int i=0; i<InternalDimension; i++)
498 (*this)[i] = (*this)[i] / sum;
504 * Min/Max-Normalization
506 template<class T, unsigned int NOdfDirections>
507 OrientationDistributionFunction<T, NOdfDirections>
508 OrientationDistributionFunction<T, NOdfDirections>
509 ::MinMaxNormalize() const
511 T max = NumericTraits<T>::NonpositiveMin();
512 T min = NumericTraits<T>::max();
513 for( unsigned int i=0; i<InternalDimension; i++)
515 max = (*this)[i] > max ? (*this)[i] : max;
516 min = (*this)[i] < min ? (*this)[i] : min;
519 for( unsigned int i=0; i<InternalDimension; i++)
521 retval[i] = ((*this)[i] - min) / (max - min);
529 template<class T, unsigned int NOdfDirections>
530 OrientationDistributionFunction<T, NOdfDirections>
531 OrientationDistributionFunction<T, NOdfDirections>
532 ::MaxNormalize() const
534 T max = NumericTraits<T>::NonpositiveMin();
535 for( unsigned int i=0; i<InternalDimension; i++)
537 max = (*this)[i] > max ? (*this)[i] : max;
540 for( unsigned int i=0; i<InternalDimension; i++)
542 retval[i] = (*this)[i] / max;
547 template<class T, unsigned int NOdfDirections>
549 OrientationDistributionFunction<T, NOdfDirections>
550 ::GetMaxValue() const
552 T max = NumericTraits<T>::NonpositiveMin();
553 for( unsigned int i=0; i<InternalDimension; i++)
555 if((*this)[i] >= max )
563 template<class T, unsigned int NOdfDirections>
565 OrientationDistributionFunction<T, NOdfDirections>
566 ::GetMinValue() const
568 T min = NumericTraits<T>::max();
569 for( unsigned int i=0; i<InternalDimension; i++)
571 if((*this)[i] >= min )
579 template<class T, unsigned int NOdfDirections>
581 OrientationDistributionFunction<T, NOdfDirections>
582 ::GetMeanValue() const
585 for( unsigned int i=0; i<InternalDimension; i++)
587 return sum/InternalDimension;
590 template<class T, unsigned int NOdfDirections>
592 OrientationDistributionFunction<T, NOdfDirections>
593 ::GetMaxChordLength()
595 if(m_MaxChordLength<0.0)
598 double max_dist = -1;
599 vtkPoints* points = m_BaseMesh->GetPoints();
600 for(int i=0; i<NOdfDirections; i++)
603 points->GetPoint(i,p);
604 std::vector<int> neighbors = GetNeighbors(i);
605 for(std::size_t j=0; j<neighbors.size(); j++)
608 points->GetPoint(neighbors[j],n);
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;
616 m_MaxChordLength = max_dist;
618 return m_MaxChordLength;
621 template<class T, unsigned int NOdfDirections>
623 OrientationDistributionFunction<T, NOdfDirections>
627 m_MutexBaseMesh.Lock();
628 if(m_BaseMesh == nullptr)
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);
642 vtkPolyData* polydata = vtkPolyData::New();
643 polydata->SetPoints( points );
644 vtkDelaunay2D *delaunay = vtkDelaunay2D::New();
645 delaunay->SetInputData( polydata );
648 vtkCellArray* vtkpolys = delaunay->GetOutput()->GetPolys();
649 vtkCellArray* vtknewpolys = vtkCellArray::New();
650 vtkIdType npts; vtkIdType *pts;
651 while(vtkpolys->GetNextCell(npts,pts))
654 for(int i=0; i<npts; i++)
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))
663 vtknewpolys->InsertNextCell(npts, pts);
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);
677 vtkPolyData* polydata2 = vtkPolyData::New();
678 polydata2->SetPoints( points2 );
679 vtkDelaunay2D *delaunay2 = vtkDelaunay2D::New();
680 delaunay2->SetInputData( polydata2 );
683 vtkpolys = delaunay2->GetOutput()->GetPolys();
684 while(vtkpolys->GetNextCell(npts,pts))
687 for(int i=0; i<npts; i++)
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))
696 vtknewpolys->InsertNextCell(npts, pts);
699 polydata->SetPolys(vtknewpolys);
701 for (unsigned int p = 0; p < NOdfDirections; p++)
703 points->SetPoint(p,m_Directions->get_column(p).data_block());
705 polydata->SetPoints( points );
707 m_BaseMesh = polydata;
709 m_MutexBaseMesh.Unlock();
713 * Extract the principle diffusion direction
715 template<class T, unsigned int NOdfDirections>
717 OrientationDistributionFunction<T, NOdfDirections>
718 ::GetPrincipleDiffusionDirection() const
720 T max = NumericTraits<T>::NonpositiveMin();
722 for( unsigned int i=0; i<InternalDimension; i++)
724 if((*this)[i] >= max )
734 template<class T, unsigned int NOdfDirections>
736 OrientationDistributionFunction<T, NOdfDirections>
737 ::GetNeighbors(int idx)
741 m_MutexNeighbors.Lock();
742 if(m_NeighborIdxs == nullptr)
744 m_NeighborIdxs = new std::vector< std::vector<int>* >();
745 vtkCellArray* polys = m_BaseMesh->GetPolys();
747 for(unsigned int i=0; i<NOdfDirections; i++)
749 auto idxs = new std::vector<int>();
750 polys->InitTraversal();
751 vtkIdType npts; vtkIdType *pts;
752 while(polys->GetNextCell(npts,pts))
756 idxs->push_back(pts[1]);
757 idxs->push_back(pts[2]);
759 else if( pts[1] == i )
761 idxs->push_back(pts[0]);
762 idxs->push_back(pts[2]);
764 else if( pts[2] == i )
766 idxs->push_back(pts[0]);
767 idxs->push_back(pts[1]);
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);
777 m_MutexNeighbors.Unlock();
779 return *m_NeighborIdxs->at(idx);
783 * Extract the n-th diffusion direction
785 template<class T, unsigned int NOdfDirections>
787 OrientationDistributionFunction<T, NOdfDirections>
788 ::GetNthDiffusionDirection(int n, vnl_vector_fixed<double,3> rndVec) const
792 return GetPrincipleDiffusionDirection();
794 m_MutexHalfSphereIdxs.Lock();
795 if( !m_HalfSphereIdxs )
797 m_HalfSphereIdxs = new std::vector<int>();
798 for( unsigned int i=0; i<InternalDimension; i++)
800 if(dot_product(m_Directions->get_column(i),rndVec) > 0.0)
802 m_HalfSphereIdxs->push_back(i);
806 m_MutexHalfSphereIdxs.Unlock();
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();
816 std::vector<int> nbs = GetNeighbors(*it);
817 std::vector<int>::iterator it2;
819 for(it2 = nbs.begin();
823 if((*this)[*it2] > (*this)[*it])
830 localMaxima.push_back(*it);
833 // delete n highest local maxima from list
834 // and return remaining highest
836 std::vector<int>::iterator itMax;
837 for( int i=0; i<=n; i++ )
840 T max = NumericTraits<T>::NonpositiveMin();
841 for(it = localMaxima.begin();
842 it != localMaxima.end();
851 it = find(localMaxima.begin(), localMaxima.end(), maxidx);
852 if(it!=localMaxima.end())
853 localMaxima.erase(it);
859 template < typename TComponent, unsigned int NOdfDirections >
860 vnl_vector_fixed<double,3> itk::OrientationDistributionFunction<TComponent, NOdfDirections>
861 ::GetDirection( int i )
863 return m_Directions->get_column(i);
867 * Interpolate a position between sampled directions
869 template<class T, unsigned int NOdfDirections>
871 OrientationDistributionFunction<T, NOdfDirections>
872 ::GetInterpolatedComponent(vnl_vector_fixed<double,3> dir, InterpolationMethods method) const
876 double retval = -1.0;
880 case ODF_NEAREST_NEIGHBOR_INTERP:
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++)
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;
893 retval = this->GetNthComponent(current_min_idx);
897 case ODF_TRILINEAR_BARYCENTRIC_INTERP:
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))
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]));
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)
923 vnl_vector_fixed<double,3> v0 = C - A;
924 vnl_vector_fixed<double,3> v1 = B - A;
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());
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;
941 // Check if point is in triangle
943 if( (uv(0) >= 0-eps) && (uv(1) >= 0-eps) && (uv(0) + uv(1) <= 1+eps) )
945 // check if minimum angle is the max so far
946 if(d1.two_norm() < current_min)
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]);
962 case ODF_SPHERICAL_GAUSSIAN_BASIS_FUNCTIONS:
964 double maxChordLength = GetMaxChordLength();
965 double sigma = asin(maxChordLength/2);
967 // this is the contribution of each kernel to each sampling point on the
969 vnl_vector<double> contrib;
970 contrib.set_size(NOdfDirections);
972 vtkPoints* points = m_BaseMesh->GetPoints();
974 for(int i=0; i<NOdfDirections; i++)
976 vnl_vector_fixed<double,3> P(points->GetPoint(i));
977 double stv = dir[0]*P[0]
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));
988 for(int i=0; i<NOdfDirections; i++)
990 retval += (contrib[i] / sum)*this->GetNthComponent(i);
999 std::cout << "Interpolation failed" << std::endl;
1008 * Calculate Generalized Fractional Anisotropy
1010 template<class T, unsigned int NOdfDirections>
1012 OrientationDistributionFunction<T, NOdfDirections>
1013 ::GetGeneralizedFractionalAnisotropy() const
1018 for( unsigned int i=0; i<InternalDimension; i++)
1023 mean /= NOdfDirections;
1024 for( unsigned int i=0; i<InternalDimension; i++)
1027 std += (val - mean) * (val - mean);
1030 std *= NOdfDirections;
1031 rms *= NOdfDirections - 1;
1039 return sqrt(std/rms);
1044 template < typename T, unsigned int N>
1045 T itk::OrientationDistributionFunction<T, N>
1046 ::GetGeneralizedGFA( int k, int p ) const
1051 double max = NumericTraits<double>::NonpositiveMin();
1052 for( unsigned int i=0; i<InternalDimension; i++)
1054 double val = (double)(*this)[i];
1055 mean += pow(val,(double)p);
1056 max = val > max ? val : max;
1058 max = pow(max,(double)p);
1060 for( unsigned int i=0; i<InternalDimension; i++)
1062 double val = (double)(*this)[i];
1063 std += (pow(val,(double)p) - mean) * (pow(val,(double)p) - mean);
1066 rms += pow(val,(double)(p*k));
1075 rms = pow(rms,(double)(1.0/k));
1077 else if(k<0) // lim k->inf gives us the maximum
1081 else // k==0 undefined, we define zeros root from 1 as 1
1092 return (T)(std/rms);
1097 * Calculate Nematic Order Parameter
1099 template < typename T, unsigned int N >
1100 T itk::OrientationDistributionFunction<T, N>
1101 ::GetNematicOrderParameter() const
1103 // not yet implemented
1108 * Calculate StdDev by MaxValue
1110 template < typename T, unsigned int N >
1111 T itk::OrientationDistributionFunction<T, N>
1112 ::GetStdDevByMaxValue() const
1116 T max = NumericTraits<T>::NonpositiveMin();
1117 for( unsigned int i=0; i<InternalDimension; i++)
1121 max = (*this)[i] > max ? (*this)[i] : max;
1123 mean /= InternalDimension;
1124 for( unsigned int i=0; i<InternalDimension; i++)
1127 std += (val - mean) * (val - mean);
1129 std /= InternalDimension-1;
1137 return (sqrt(std)/max);
1141 template < typename T, unsigned int N >
1142 T itk::OrientationDistributionFunction<T, N>
1143 ::GetPrincipleCurvature(double alphaMinDegree, double alphaMaxDegree, int invert) const
1145 // following loop only performed once
1146 // (computing indices of each angular range)
1147 m_MutexAngularRange.Lock();
1148 if(m_AngularRangeIdxs == nullptr)
1150 m_AngularRangeIdxs = new std::vector< std::vector<int>* >();
1151 for(unsigned int i=0; i<N; i++)
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++)
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) )
1164 m_AngularRangeIdxs->push_back(idxs);
1167 m_MutexAngularRange.Unlock();
1169 // find the maximum (or minimum) direction (remember index and value)
1174 pIdx = GetPrincipleDiffusionDirection();
1175 mode = (*this)[pIdx];
1179 mode = NumericTraits<T>::max();
1180 for( unsigned int i=0; i<N; i++)
1182 if((*this)[i] < mode )
1189 //////////////////////
1190 //////// compute median of mode and its neighbors to become more stable to noise
1191 //////// compared to simply using the mode
1192 //////////////////////
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++)
1201 ////// modeAndNeighborVals.push_back((*this)[nbs[i]]);
1204 //////// sort by value
1205 //////std::sort( modeAndNeighborVals.begin(), modeAndNeighborVals.end() );
1207 //////// median of mode and neighbors
1208 //////mode = modeAndNeighborVals[floor(0.5*(double)(numNeighbors+1)+0.5)];
1211 // computing a quantile of the angular range
1215 double quantile = 0.00;
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++)
1222 odfValuesInAngularRange.push_back((*this)[(*m_AngularRangeIdxs->at(pIdx))[i] ]);
1225 // sort them by value
1226 std::sort( odfValuesInAngularRange.begin(), odfValuesInAngularRange.end() );
1228 // median of angular range
1229 T median = odfValuesInAngularRange[floor(quantile*(double)numInRange+0.5)];
1231 // compute and return final value
1234 return mode/median - 1.0;
1238 return median/mode - 1.0;
1243 * Calculate Normalized Entropy
1245 template < typename T, unsigned int N >
1246 T itk::OrientationDistributionFunction<T, N>
1247 ::GetNormalizedEntropy() const
1250 for( unsigned int i=0; i<InternalDimension; i++)
1259 val = log(0.0000001);
1263 double _n = (double) InternalDimension;
1265 return (T) (-_n / log(_n) * mean);
1269 * Pre-multiply the Tensor by a Matrix
1271 template<class T, unsigned int NOdfDirections>
1272 OrientationDistributionFunction<T, NOdfDirections>
1273 OrientationDistributionFunction<T, NOdfDirections>
1274 ::PreMultiply( const MatrixType & m ) const
1277 typedef typename NumericTraits<T>::AccumulateType AccumulateType;
1278 for(unsigned int r=0; r<NOdfDirections; r++)
1280 for(unsigned int c=0; c<NOdfDirections; c++)
1282 AccumulateType sum = NumericTraits<AccumulateType>::ZeroValue();
1283 for(unsigned int t=0; t<NOdfDirections; t++)
1285 sum += m(r,t) * (*this)(t,c);
1287 result(r,c) = static_cast<T>( sum );
1294 * Post-multiply the Tensor by a Matrix
1296 template<class T, unsigned int NOdfDirections>
1297 OrientationDistributionFunction<T, NOdfDirections>
1298 OrientationDistributionFunction<T, NOdfDirections>
1299 ::PostMultiply( const MatrixType & m ) const
1302 typedef typename NumericTraits<T>::AccumulateType AccumulateType;
1303 for(unsigned int r=0; r<NOdfDirections; r++)
1305 for(unsigned int c=0; c<NOdfDirections; c++)
1307 AccumulateType sum = NumericTraits<AccumulateType>::ZeroValue();
1308 for(unsigned int t=0; t<NOdfDirections; t++)
1310 sum += (*this)(r,t) * m(t,c);
1312 result(r,c) = static_cast<T>( sum );
1319 * Print content to an ostream
1321 template<class T, unsigned int NOdfDirections>
1323 operator<<(std::ostream& os,const OrientationDistributionFunction<T, NOdfDirections> & c )
1325 for(unsigned int i=0; i<c.GetNumberOfComponents(); i++)
1327 os << static_cast<typename NumericTraits<T>::PrintType>(c[i]) << " ";
1333 * Read content from an istream
1335 template<class T, unsigned int NOdfDirections>
1337 operator>>(std::istream& is, OrientationDistributionFunction<T, NOdfDirections> & dt )
1339 for(unsigned int i=0; i < dt.GetNumberOfComponents(); i++)
1346 } // end namespace itk