itkOrientationDistributionFunction.txx
Go to the documentation of this file.
00001 /*========================================================================= 00002 00003 Program: Insight Segmentation & Registration Toolkit 00004 Language: C++ 00005 Date: $Date: 2008-03-10 22:48:13 $ 00006 Version: $Revision: 1.14 $ 00007 00008 Copyright (c) Insight Software Consortium. All rights reserved. 00009 See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details. 00010 00011 This software is distributed WITHOUT ANY WARRANTY; without even 00012 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 00013 PURPOSE. See the above copyright notices for more information. 00014 00015 =========================================================================*/ 00016 #ifndef _itkOrientationDistributionFunction_txx 00017 #define _itkOrientationDistributionFunction_txx 00018 00019 #include <stdio.h> 00020 #include <stdlib.h> 00021 #include <time.h> 00022 #include <algorithm> 00023 #include <vector> 00024 #include <vnl/algo/vnl_matrix_inverse.h> 00025 #include "itkPointShell.h" 00026 00027 //#include "itkNumericTraitsTensorPixel.h" 00028 00029 namespace itk 00030 { 00031 00032 /* 00033 00034 #define INIT_STATIC_ODF_VARS(N_DIRS) \ 00035 _INIT_STATIC_ODF_VARS(float,N_DIRS) \ 00036 _INIT_STATIC_ODF_VARS(double,N_DIRS) \ 00037 00038 #define _INIT_STATIC_ODF_VARS(PIXTYPE,N_DIRS) \ 00039 vtkPolyData* itk::OrientationDistributionFunction<PIXTYPE,N_DIRS>::m_BaseMesh = NULL; \ 00040 vnl_matrix_fixed<double, 3, N_DIRS>* itk::OrientationDistributionFunction<PIXTYPE,N_DIRS>::m_Directions \ 00041 = itk::PointShell<N_DIRS, vnl_matrix_fixed<double, 3, N_DIRS> >::DistributePointShell(); \ 00042 std::vector< std::vector<int>* >* itk::OrientationDistributionFunction<PIXTYPE,N_DIRS>::m_NeighborIdxs = NULL; \ 00043 std::vector<int>* itk::OrientationDistributionFunction<PIXTYPE,N_DIRS>::m_HalfSphereIdxs = NULL; \ 00044 bool itk::OrientationDistributionFunction<PIXTYPE,N_DIRS>::m_Mutex = false; \ 00045 00046 00047 INIT_STATIC_ODF_VARS(40) 00048 INIT_STATIC_ODF_VARS(60) 00049 INIT_STATIC_ODF_VARS(80) 00050 INIT_STATIC_ODF_VARS(100) 00051 INIT_STATIC_ODF_VARS(150) 00052 INIT_STATIC_ODF_VARS(200) 00053 INIT_STATIC_ODF_VARS(250) 00054 */ 00055 00056 template<class T, unsigned int N> 00057 vtkPolyData* itk::OrientationDistributionFunction<T,N>::m_BaseMesh = NULL; 00058 00059 template<class T, unsigned int N> 00060 double itk::OrientationDistributionFunction<T,N>::m_MaxChordLength = -1.0; 00061 00062 template<class T, unsigned int N> 00063 vnl_matrix_fixed<double, 3, N>* itk::OrientationDistributionFunction<T,N>::m_Directions 00064 = itk::PointShell<N, vnl_matrix_fixed<double, 3, N> >::DistributePointShell(); 00065 00066 template<class T, unsigned int N> 00067 std::vector< std::vector<int>* >* itk::OrientationDistributionFunction<T,N>::m_NeighborIdxs = NULL; 00068 00069 template<class T, unsigned int N> 00070 std::vector< std::vector<int>* >* itk::OrientationDistributionFunction<T,N>::m_AngularRangeIdxs = NULL; 00071 00072 template<class T, unsigned int N> 00073 std::vector<int>* itk::OrientationDistributionFunction<T,N>::m_HalfSphereIdxs = NULL; 00074 00075 template<class T, unsigned int N> 00076 itk::SimpleFastMutexLock itk::OrientationDistributionFunction<T,N>::m_MutexBaseMesh; 00077 template<class T, unsigned int N> 00078 itk::SimpleFastMutexLock itk::OrientationDistributionFunction<T,N>::m_MutexHalfSphereIdxs; 00079 template<class T, unsigned int N> 00080 itk::SimpleFastMutexLock itk::OrientationDistributionFunction<T,N>::m_MutexNeighbors; 00081 template<class T, unsigned int N> 00082 itk::SimpleFastMutexLock itk::OrientationDistributionFunction<T,N>::m_MutexAngularRange; 00083 00084 00085 #define ODF_PI 3.14159265358979323846 00086 00087 /* 00088 * Assignment Operator 00089 */ 00090 template<class T, unsigned int NOdfDirections> 00091 OrientationDistributionFunction<T, NOdfDirections>& 00092 OrientationDistributionFunction<T, NOdfDirections> 00093 ::operator= (const Self& r) 00094 { 00095 BaseArray::operator=(r); 00096 return *this; 00097 } 00098 00099 00100 /* 00101 * Assignment Operator from a scalar constant 00102 */ 00103 template<class T, unsigned int NOdfDirections> 00104 OrientationDistributionFunction<T, NOdfDirections>& 00105 OrientationDistributionFunction<T, NOdfDirections> 00106 ::operator= (const ComponentType & r) 00107 { 00108 BaseArray::operator=(&r); 00109 return *this; 00110 } 00111 00112 00113 /* 00114 * Assigment from a plain array 00115 */ 00116 template<class T, unsigned int NOdfDirections> 00117 OrientationDistributionFunction<T, NOdfDirections>& 00118 OrientationDistributionFunction<T, NOdfDirections> 00119 ::operator= (const ComponentArrayType r ) 00120 { 00121 BaseArray::operator=(r); 00122 return *this; 00123 } 00124 00125 00126 00130 template<class T, unsigned int NOdfDirections> 00131 OrientationDistributionFunction<T, NOdfDirections> 00132 OrientationDistributionFunction<T, NOdfDirections> 00133 ::operator+(const Self & r) const 00134 { 00135 Self result; 00136 for( unsigned int i=0; i<InternalDimension; i++) 00137 { 00138 result[i] = (*this)[i] + r[i]; 00139 } 00140 return result; 00141 } 00142 00143 00144 00145 00149 template<class T, unsigned int NOdfDirections> 00150 OrientationDistributionFunction<T, NOdfDirections> 00151 OrientationDistributionFunction<T, NOdfDirections> 00152 ::operator-(const Self & r) const 00153 { 00154 Self result; 00155 for( unsigned int i=0; i<InternalDimension; i++) 00156 { 00157 result[i] = (*this)[i] - r[i]; 00158 } 00159 return result; 00160 } 00161 00162 00163 00167 template<class T, unsigned int NOdfDirections> 00168 const OrientationDistributionFunction<T, NOdfDirections> & 00169 OrientationDistributionFunction<T, NOdfDirections> 00170 ::operator+=(const Self & r) 00171 { 00172 for( unsigned int i=0; i<InternalDimension; i++) 00173 { 00174 (*this)[i] += r[i]; 00175 } 00176 return *this; 00177 } 00178 00179 00180 00181 00185 template<class T, unsigned int NOdfDirections> 00186 const OrientationDistributionFunction<T, NOdfDirections> & 00187 OrientationDistributionFunction<T, NOdfDirections> 00188 ::operator-=(const Self & r) 00189 { 00190 for( unsigned int i=0; i<InternalDimension; i++) 00191 { 00192 (*this)[i] -= r[i]; 00193 } 00194 return *this; 00195 } 00196 00197 00198 00202 template<class T, unsigned int NOdfDirections> 00203 const OrientationDistributionFunction<T, NOdfDirections> & 00204 OrientationDistributionFunction<T, NOdfDirections> 00205 ::operator*=(const RealValueType & r) 00206 { 00207 for( unsigned int i=0; i<InternalDimension; i++) 00208 { 00209 (*this)[i] *= r; 00210 } 00211 return *this; 00212 } 00213 00214 00215 00219 template<class T, unsigned int NOdfDirections> 00220 const OrientationDistributionFunction<T, NOdfDirections> & 00221 OrientationDistributionFunction<T, NOdfDirections> 00222 ::operator/=(const RealValueType & r) 00223 { 00224 for( unsigned int i=0; i<InternalDimension; i++) 00225 { 00226 (*this)[i] /= r; 00227 } 00228 return *this; 00229 } 00230 00231 00232 00233 00237 template<class T, unsigned int NOdfDirections> 00238 OrientationDistributionFunction<T, NOdfDirections> 00239 OrientationDistributionFunction<T, NOdfDirections> 00240 ::operator*(const RealValueType & r) const 00241 { 00242 Self result; 00243 for( unsigned int i=0; i<InternalDimension; i++) 00244 { 00245 result[i] = (*this)[i] * r; 00246 } 00247 return result; 00248 } 00249 00250 00254 template<class T, unsigned int NOdfDirections> 00255 OrientationDistributionFunction<T, NOdfDirections> 00256 OrientationDistributionFunction<T, NOdfDirections> 00257 ::operator/(const RealValueType & r) const 00258 { 00259 Self result; 00260 for( unsigned int i=0; i<InternalDimension; i++) 00261 { 00262 result[i] = (*this)[i] / r; 00263 } 00264 return result; 00265 } 00266 00267 00268 /* 00269 * Matrix notation access to elements 00270 */ 00271 template<class T, unsigned int NOdfDirections> 00272 const typename OrientationDistributionFunction<T, NOdfDirections>::ValueType & 00273 OrientationDistributionFunction<T, NOdfDirections> 00274 ::operator()(unsigned int row, unsigned int col) const 00275 { 00276 unsigned int k; 00277 00278 if( row < col ) 00279 { 00280 k = row * InternalDimension + col - row * ( row + 1 ) / 2; 00281 } 00282 else 00283 { 00284 k = col * InternalDimension + row - col * ( col + 1 ) / 2; 00285 } 00286 00287 00288 if( k >= InternalDimension ) 00289 { 00290 k = 0; 00291 } 00292 00293 return (*this)[k]; 00294 } 00295 00296 00297 00298 /* 00299 * Matrix notation access to elements 00300 */ 00301 template<class T, unsigned int NOdfDirections> 00302 typename OrientationDistributionFunction<T, NOdfDirections>::ValueType & 00303 OrientationDistributionFunction<T, NOdfDirections> 00304 ::operator()(unsigned int row, unsigned int col) 00305 { 00306 unsigned int k; 00307 00308 if( row < col ) 00309 { 00310 k = row * InternalDimension + col - row * ( row + 1 ) / 2; 00311 } 00312 else 00313 { 00314 k = col * InternalDimension + row - col * ( col + 1 ) / 2; 00315 } 00316 00317 00318 if( k >= InternalDimension ) 00319 { 00320 k = 0; 00321 } 00322 00323 return (*this)[k]; 00324 } 00325 00326 00327 /* 00328 * Set the Tensor to an Identity. 00329 * Set ones in the diagonal and zeroes every where else. 00330 */ 00331 template<class T, unsigned int NOdfDirections> 00332 void 00333 OrientationDistributionFunction<T, NOdfDirections> 00334 ::SetIsotropic() 00335 { 00336 this->Fill(NumericTraits< T >::One / NOdfDirections); 00337 } 00338 00339 /* 00340 * Set the Tensor to an Identity. 00341 * Set ones in the diagonal and zeroes every where else. 00342 */ 00343 template<class T, unsigned int NOdfDirections> 00344 void 00345 OrientationDistributionFunction<T, NOdfDirections> 00346 ::InitFromTensor(itk::DiffusionTensor3D<T> tensor) 00347 { 00348 for(unsigned int i=0; i<NOdfDirections; i++) 00349 { 00350 /* 00351 * | t0 t1 t2 | g0 00352 * g0 g1 g2 * | t1 t3 t4 | * g1 00353 * | t2 t4 t5 | g2 00354 * 00355 * = g0 * (t0g0*t1g1*t2g2) 00356 * + g1 * (t1g0+t3g1+t4g2) 00357 * + g2 * (t2g0+t4g1+t5g2) 00358 */ 00359 T g0 = (*m_Directions)(0,i); 00360 T g1 = (*m_Directions)(1,i); 00361 T g2 = (*m_Directions)(2,i); 00362 T t0 = tensor[0]; 00363 T t1 = tensor[1]; 00364 T t2 = tensor[2]; 00365 T t3 = tensor[3]; 00366 T t4 = tensor[4]; 00367 T t5 = tensor[5]; 00368 (*this)[i] = g0 * (t0*g0+t1*g1+t2*g2) 00369 + g1 * (t1*g0+t3*g1+t4*g2) 00370 + g2 * (t2*g0+t4*g1+t5*g2); 00371 00372 (*this)[i] = (*this)[i] < 0 ? 0 : (*this)[i]; 00373 } 00374 } 00375 00376 /* 00377 * Normalization to PDF 00378 */ 00379 template<class T, unsigned int NOdfDirections> 00380 void 00381 OrientationDistributionFunction<T, NOdfDirections> 00382 ::Normalize() 00383 { 00384 T sum = 0; 00385 for( unsigned int i=0; i<InternalDimension; i++) 00386 { 00387 sum += (*this)[i]; 00388 } 00389 for( unsigned int i=0; i<InternalDimension; i++) 00390 { 00391 (*this)[i] = (*this)[i] / sum; 00392 } 00393 } 00394 00395 /* 00396 * Min/Max-Normalization 00397 */ 00398 template<class T, unsigned int NOdfDirections> 00399 OrientationDistributionFunction<T, NOdfDirections> 00400 OrientationDistributionFunction<T, NOdfDirections> 00401 ::MinMaxNormalize() const 00402 { 00403 T max = NumericTraits<T>::NonpositiveMin(); 00404 T min = NumericTraits<T>::max(); 00405 for( unsigned int i=0; i<InternalDimension; i++) 00406 { 00407 max = (*this)[i] > max ? (*this)[i] : max; 00408 min = (*this)[i] < min ? (*this)[i] : min; 00409 } 00410 Self retval; 00411 for( unsigned int i=0; i<InternalDimension; i++) 00412 { 00413 retval[i] = ((*this)[i] - min) / (max - min); 00414 } 00415 return retval; 00416 } 00417 00418 /* 00419 * Max-Normalization 00420 */ 00421 template<class T, unsigned int NOdfDirections> 00422 OrientationDistributionFunction<T, NOdfDirections> 00423 OrientationDistributionFunction<T, NOdfDirections> 00424 ::MaxNormalize() const 00425 { 00426 T max = NumericTraits<T>::NonpositiveMin(); 00427 for( unsigned int i=0; i<InternalDimension; i++) 00428 { 00429 max = (*this)[i] > max ? (*this)[i] : max; 00430 } 00431 Self retval; 00432 for( unsigned int i=0; i<InternalDimension; i++) 00433 { 00434 retval[i] = (*this)[i] / max; 00435 } 00436 return retval; 00437 } 00438 00439 template<class T, unsigned int NOdfDirections> 00440 double 00441 OrientationDistributionFunction<T, NOdfDirections> 00442 ::GetMaxChordLength() 00443 { 00444 if(m_MaxChordLength<0.0) 00445 { 00446 ComputeBaseMesh(); 00447 double max_dist = -1; 00448 vtkPoints* points = m_BaseMesh->GetPoints(); 00449 for(int i=0; i<NOdfDirections; i++) 00450 { 00451 double p[3]; 00452 points->GetPoint(i,p); 00453 std::vector<int> neighbors = GetNeighbors(i); 00454 for(int j=0; j<neighbors.size(); j++) 00455 { 00456 double n[3]; 00457 points->GetPoint(neighbors[j],n); 00458 double d = sqrt( 00459 (p[0]-n[0])*(p[0]-n[0]) + 00460 (p[1]-n[1])*(p[1]-n[1]) + 00461 (p[2]-n[2])*(p[2]-n[2])); 00462 max_dist = d>max_dist ? d : max_dist; 00463 } 00464 } 00465 m_MaxChordLength = max_dist; 00466 } 00467 return m_MaxChordLength; 00468 } 00469 00470 template<class T, unsigned int NOdfDirections> 00471 void 00472 OrientationDistributionFunction<T, NOdfDirections> 00473 ::ComputeBaseMesh() 00474 { 00475 00476 m_MutexBaseMesh.Lock(); 00477 if(m_BaseMesh == NULL) 00478 { 00479 00480 vtkPoints* points = vtkPoints::New(); 00481 for(unsigned int j=0; j<NOdfDirections; j++){ 00482 double x = (*m_Directions)(0,j); 00483 double y = (*m_Directions)(1,j); 00484 double z = (*m_Directions)(2,j); 00485 double az = atan2(y,x); 00486 double elev = atan2(z,sqrt(x*x+y*y)); 00487 double r = sqrt(x*x+y*y+z*z); 00488 points->InsertNextPoint(az,elev,r); 00489 } 00490 00491 vtkPolyData* polydata = vtkPolyData::New(); 00492 polydata->SetPoints( points ); 00493 vtkDelaunay2D *delaunay = vtkDelaunay2D::New(); 00494 delaunay->SetInput( polydata ); 00495 delaunay->Update(); 00496 00497 vtkCellArray* vtkpolys = delaunay->GetOutput()->GetPolys(); 00498 vtkCellArray* vtknewpolys = vtkCellArray::New(); 00499 vtkIdType npts; vtkIdType *pts; 00500 while(vtkpolys->GetNextCell(npts,pts)) 00501 { 00502 bool insert = true; 00503 for(int i=0; i<npts; i++) 00504 { 00505 double *tmpPoint = points->GetPoint(pts[i]); 00506 double az = tmpPoint[0]; 00507 double elev = tmpPoint[1]; 00508 if((abs(az)>ODF_PI-0.5) || (abs(elev)>ODF_PI/2-0.5)) 00509 insert = false; 00510 } 00511 if(insert) 00512 vtknewpolys->InsertNextCell(npts, pts); 00513 } 00514 00515 vtkPoints* points2 = vtkPoints::New(); 00516 for(unsigned int j=0; j<NOdfDirections; j++){ 00517 double x = -(*m_Directions)(0,j); 00518 double y = -(*m_Directions)(2,j); 00519 double z = -(*m_Directions)(1,j); 00520 double az = atan2(y,x); 00521 double elev = atan2(z,sqrt(x*x+y*y)); 00522 double r = sqrt(x*x+y*y+z*z); 00523 points2->InsertNextPoint(az,elev,r); 00524 } 00525 00526 vtkPolyData* polydata2 = vtkPolyData::New(); 00527 polydata2->SetPoints( points2 ); 00528 vtkDelaunay2D *delaunay2 = vtkDelaunay2D::New(); 00529 delaunay2->SetInput( polydata2 ); 00530 delaunay2->Update(); 00531 00532 vtkpolys = delaunay2->GetOutput()->GetPolys(); 00533 while(vtkpolys->GetNextCell(npts,pts)) 00534 { 00535 bool insert = true; 00536 for(int i=0; i<npts; i++) 00537 { 00538 double *tmpPoint = points2->GetPoint(pts[i]); 00539 double az = tmpPoint[0]; 00540 double elev = tmpPoint[1]; 00541 if((abs(az)>ODF_PI-0.5) || (abs(elev)>ODF_PI/2-0.5)) 00542 insert = false; 00543 } 00544 if(insert) 00545 vtknewpolys->InsertNextCell(npts, pts); 00546 } 00547 00548 polydata->SetPolys(vtknewpolys); 00549 00550 for (vtkIdType p = 0; p < NOdfDirections; p++) 00551 { 00552 points->SetPoint(p,m_Directions->get_column(p).data_block()); 00553 } 00554 polydata->SetPoints( points ); 00555 00557 //vtkCleanPolyData* cleaner = vtkCleanPolyData::New(); 00558 //cleaner->SetInput( polydata ); 00559 //cleaner->SetAbsoluteTolerance( 0.0 ); 00560 //cleaner->Update(); 00561 00562 //vtkAppendPolyData* append = vtkAppendPolyData::New(); 00563 //append->AddInput(cleaner->GetOutput()); 00564 //append->Update(); 00565 //m_BaseMesh = append->GetOutput(); 00566 00567 m_BaseMesh = polydata; 00568 } 00569 m_MutexBaseMesh.Unlock(); 00570 } 00571 00572 /* 00573 * Extract the principle diffusion direction 00574 */ 00575 template<class T, unsigned int NOdfDirections> 00576 int 00577 OrientationDistributionFunction<T, NOdfDirections> 00578 ::GetPrincipleDiffusionDirection() const 00579 { 00580 T max = NumericTraits<T>::NonpositiveMin(); 00581 int maxidx = -1; 00582 for( unsigned int i=0; i<InternalDimension; i++) 00583 { 00584 if((*this)[i] > max ) 00585 { 00586 max = (*this)[i]; 00587 maxidx = i; 00588 } 00589 } 00590 00591 return maxidx; 00592 } 00593 00594 template<class T, unsigned int NOdfDirections> 00595 std::vector<int> 00596 OrientationDistributionFunction<T, NOdfDirections> 00597 ::GetNeighbors(int idx) 00598 { 00599 ComputeBaseMesh(); 00600 00601 m_MutexNeighbors.Lock(); 00602 if(m_NeighborIdxs == NULL) 00603 { 00604 m_NeighborIdxs = new std::vector< std::vector<int>* >(); 00605 vtkCellArray* polys = m_BaseMesh->GetPolys(); 00606 00607 for(unsigned int i=0; i<NOdfDirections; i++) 00608 { 00609 std::vector<int> *idxs = new std::vector<int>(); 00610 polys->InitTraversal(); 00611 vtkIdType npts; vtkIdType *pts; 00612 while(polys->GetNextCell(npts,pts)) 00613 { 00614 if( pts[0] == i ) 00615 { 00616 idxs->push_back(pts[1]); 00617 idxs->push_back(pts[2]); 00618 } 00619 else if( pts[1] == i ) 00620 { 00621 idxs->push_back(pts[0]); 00622 idxs->push_back(pts[2]); 00623 } 00624 else if( pts[2] == i ) 00625 { 00626 idxs->push_back(pts[0]); 00627 idxs->push_back(pts[1]); 00628 } 00629 } 00630 std::sort(idxs->begin(), idxs->end()); 00631 std::vector< int >::iterator endLocation; 00632 endLocation = std::unique( idxs->begin(), idxs->end() ); 00633 idxs->erase(endLocation, idxs->end()); 00634 m_NeighborIdxs->push_back(idxs); 00635 } 00636 } 00637 m_MutexNeighbors.Unlock(); 00638 00639 return *m_NeighborIdxs->at(idx); 00640 } 00641 00642 /* 00643 * Extract the n-th diffusion direction 00644 */ 00645 template<class T, unsigned int NOdfDirections> 00646 int 00647 OrientationDistributionFunction<T, NOdfDirections> 00648 ::GetNthDiffusionDirection(int n, vnl_vector_fixed<double,3> rndVec) const 00649 { 00650 00651 if( n == 0 ) 00652 return GetPrincipleDiffusionDirection(); 00653 00654 m_MutexHalfSphereIdxs.Lock(); 00655 if( !m_HalfSphereIdxs ) 00656 { 00657 m_HalfSphereIdxs = new std::vector<int>(); 00658 for( unsigned int i=0; i<InternalDimension; i++) 00659 { 00660 if(dot_product(m_Directions->get_column(i),rndVec) > 0.0) 00661 { 00662 m_HalfSphereIdxs->push_back(i); 00663 } 00664 } 00665 } 00666 m_MutexHalfSphereIdxs.Unlock(); 00667 00668 // collect indices of directions 00669 // that are local maxima 00670 std::vector<int> localMaxima; 00671 std::vector<int>::iterator it; 00672 for( it=m_HalfSphereIdxs->begin(); 00673 it!=m_HalfSphereIdxs->end(); 00674 it++) 00675 { 00676 std::vector<int> nbs = GetNeighbors(*it); 00677 std::vector<int>::iterator it2; 00678 bool max = true; 00679 for(it2 = nbs.begin(); 00680 it2 != nbs.end(); 00681 it2++) 00682 { 00683 if((*this)[*it2] > (*this)[*it]) 00684 { 00685 max = false; 00686 break; 00687 } 00688 } 00689 if(max) 00690 localMaxima.push_back(*it); 00691 } 00692 00693 // delete n highest local maxima from list 00694 // and return remaining highest 00695 int maxidx = -1; 00696 std::vector<int>::iterator itMax; 00697 for( int i=0; i<=n; i++ ) 00698 { 00699 maxidx = -1; 00700 T max = NumericTraits<T>::NonpositiveMin(); 00701 for(it = localMaxima.begin(); 00702 it != localMaxima.end(); 00703 it++) 00704 { 00705 if((*this)[*it]>max) 00706 { 00707 max = (*this)[*it]; 00708 maxidx = *it; 00709 } 00710 } 00711 it = find(localMaxima.begin(), localMaxima.end(), maxidx); 00712 if(it!=localMaxima.end()) 00713 localMaxima.erase(it); 00714 } 00715 00716 return maxidx; 00717 } 00718 00719 template < typename TComponent, unsigned int NOdfDirections > 00720 vnl_vector_fixed<double,3> itk::OrientationDistributionFunction<TComponent, NOdfDirections> 00721 ::GetDirection( int i ) 00722 { 00723 return m_Directions->get_column(i); 00724 } 00725 00726 /* 00727 * Interpolate a position between sampled directions 00728 */ 00729 template<class T, unsigned int NOdfDirections> 00730 T 00731 OrientationDistributionFunction<T, NOdfDirections> 00732 ::GetInterpolatedComponent(vnl_vector_fixed<double,3> dir, InterpolationMethods method) const 00733 { 00734 00735 ComputeBaseMesh(); 00736 double retval = -1.0; 00737 00738 switch(method) 00739 { 00740 case ODF_NEAREST_NEIGHBOR_INTERP: 00741 { 00742 00743 vtkPoints* points = m_BaseMesh->GetPoints(); 00744 double current_min = NumericTraits<double>::max(); 00745 int current_min_idx = -1; 00746 for(int i=0; i<NOdfDirections; i++) 00747 { 00748 vnl_vector_fixed<double,3> P(points->GetPoint(i)); 00749 double dist = (dir-P).two_norm(); 00750 current_min_idx = dist<current_min ? i : current_min_idx; 00751 current_min = dist<current_min ? dist : current_min; 00752 } 00753 retval = this->GetNthComponent(current_min_idx); 00754 break; 00755 00756 } 00757 case ODF_TRILINEAR_BARYCENTRIC_INTERP: 00758 { 00759 00760 double maxChordLength = GetMaxChordLength(); 00761 vtkCellArray* polys = m_BaseMesh->GetPolys(); 00762 vtkPoints* points = m_BaseMesh->GetPoints(); 00763 vtkIdType npts; vtkIdType *pts; 00764 double current_min = NumericTraits<double>::max(); 00765 polys->InitTraversal(); 00766 while(polys->GetNextCell(npts,pts)) 00767 { 00768 vnl_vector_fixed<double,3> A(points->GetPoint(pts[0])); 00769 vnl_vector_fixed<double,3> B(points->GetPoint(pts[1])); 00770 vnl_vector_fixed<double,3> C(points->GetPoint(pts[2])); 00771 00772 vnl_vector_fixed<double,3> d1; 00773 d1.put(0,(dir-A).two_norm()); 00774 d1.put(1,(dir-B).two_norm()); 00775 d1.put(2,(dir-C).two_norm()); 00776 double maxval = d1.max_value(); 00777 if(maxval>maxChordLength) 00778 { 00779 continue; 00780 } 00781 00782 // Compute vectors 00783 vnl_vector_fixed<double,3> v0 = C - A; 00784 vnl_vector_fixed<double,3> v1 = B - A; 00785 00786 // Project direction to plane ABC 00787 vnl_vector_fixed<double,3> v6 = dir; 00788 vnl_vector_fixed<double,3> cross = vnl_cross_3d(v0, v1); 00789 cross = cross.normalize(); 00790 vtkPlane::ProjectPoint(v6.data_block(),A.data_block(),cross.data_block(),v6.data_block()); 00791 v6 = v6-A; 00792 00793 // Calculate barycentric coords 00794 vnl_matrix_fixed<double,3,2> mat; 00795 mat.set_column(0, v0); 00796 mat.set_column(1, v1); 00797 vnl_matrix_inverse<double> inv(mat); 00798 vnl_matrix_fixed<double,2,3> inver = inv.pinverse(); 00799 vnl_vector<double> uv = inv.pinverse()*v6; 00800 00801 // Check if point is in triangle 00802 double eps = 0.01; 00803 if( (uv(0) >= 0-eps) && (uv(1) >= 0-eps) && (uv(0) + uv(1) <= 1+eps) ) 00804 { 00805 // check if minimum angle is the max so far 00806 if(d1.two_norm() < current_min) 00807 { 00808 current_min = d1.two_norm(); 00809 vnl_vector<double> barycentricCoords(3); 00810 barycentricCoords[2] = uv[0]<0 ? 0 : (uv[0]>1?1:uv[0]); 00811 barycentricCoords[1] = uv[1]<0 ? 0 : (uv[1]>1?1:uv[1]); 00812 barycentricCoords[0] = 1-(barycentricCoords[1]+barycentricCoords[2]); 00813 retval = barycentricCoords[0]*this->GetNthComponent(pts[0]) + 00814 barycentricCoords[1]*this->GetNthComponent(pts[1]) + 00815 barycentricCoords[2]*this->GetNthComponent(pts[2]); 00816 } 00817 } 00818 } 00819 00820 break; 00821 } 00822 case ODF_SPHERICAL_GAUSSIAN_BASIS_FUNCTIONS: 00823 { 00824 double maxChordLength = GetMaxChordLength(); 00825 double sigma = asin(maxChordLength/2); 00826 00827 // this is the contribution of each kernel to each sampling point on the 00828 // equator 00829 vnl_vector<double> contrib; 00830 contrib.set_size(NOdfDirections); 00831 00832 vtkPoints* points = m_BaseMesh->GetPoints(); 00833 double sum = 0; 00834 for(int i=0; i<NOdfDirections; i++) 00835 { 00836 vnl_vector_fixed<double,3> P(points->GetPoint(i)); 00837 double stv = dir[0]*P[0] 00838 + dir[1]*P[1] 00839 + dir[2]*P[2]; 00840 stv = (stv<-1.0) ? -1.0 : ( (stv>1.0) ? 1.0 : stv); 00841 double x = acos(stv); 00842 contrib[i] = (1.0/(sigma*sqrt(2.0*ODF_PI))) 00843 *exp((-x*x)/(2*sigma*sigma)); 00844 sum += contrib[i]; 00845 } 00846 00847 retval = 0; 00848 for(int i=0; i<NOdfDirections; i++) 00849 { 00850 retval += (contrib[i] / sum)*this->GetNthComponent(i); 00851 } 00852 break; 00853 } 00854 00855 } 00856 00857 if(retval==-1) 00858 { 00859 std::cout << "Interpolation failed" << std::endl; 00860 return 0; 00861 } 00862 00863 return retval; 00864 00865 } 00866 00867 /* 00868 * Calculate Generalized Fractional Anisotropy 00869 */ 00870 template<class T, unsigned int NOdfDirections> 00871 T 00872 OrientationDistributionFunction<T, NOdfDirections> 00873 ::GetGeneralizedFractionalAnisotropy() const 00874 { 00875 double mean = 0; 00876 double std = 0; 00877 double rms = 0; 00878 for( unsigned int i=0; i<InternalDimension; i++) 00879 { 00880 T val = (*this)[i]; 00881 mean += val; 00882 } 00883 mean /= NOdfDirections; 00884 for( unsigned int i=0; i<InternalDimension; i++) 00885 { 00886 T val = (*this)[i]; 00887 std += (val - mean) * (val - mean); 00888 rms += val*val; 00889 } 00890 std *= NOdfDirections; 00891 rms *= NOdfDirections - 1; 00892 00893 if(rms == 0) 00894 { 00895 return 0; 00896 } 00897 else 00898 { 00899 return sqrt(std/rms); 00900 } 00901 } 00902 00903 00904 template < typename T, unsigned int N> 00905 T itk::OrientationDistributionFunction<T, N> 00906 ::GetGeneralizedGFA( int k, int p ) const 00907 { 00908 double mean = 0; 00909 double std = 0; 00910 double rms = 0; 00911 double max = NumericTraits<double>::NonpositiveMin(); 00912 for( unsigned int i=0; i<InternalDimension; i++) 00913 { 00914 double val = (double)(*this)[i]; 00915 mean += pow(val,(double)p); 00916 max = val > max ? val : max; 00917 } 00918 max = pow(max,(double)p); 00919 mean /= N; 00920 for( unsigned int i=0; i<InternalDimension; i++) 00921 { 00922 double val = (double)(*this)[i]; 00923 std += (pow(val,(double)p) - mean) * (pow(val,(double)p) - mean); 00924 if(k>0) 00925 { 00926 rms += pow(val,(double)(p*k)); 00927 } 00928 } 00929 std /= N - 1; 00930 std = sqrt(std); 00931 00932 if(k>0) 00933 { 00934 rms /= N; 00935 rms = pow(rms,(double)(1.0/k)); 00936 } 00937 else if(k<0) // lim k->inf gives us the maximum 00938 { 00939 rms = max; 00940 } 00941 else // k==0 undefined, we define zeros root from 1 as 1 00942 { 00943 rms = 1; 00944 } 00945 00946 if(rms == 0) 00947 { 00948 return 0; 00949 } 00950 else 00951 { 00952 return (T)(std/rms); 00953 } 00954 } 00955 00956 /* 00957 * Calculate Nematic Order Parameter 00958 */ 00959 template < typename T, unsigned int N > 00960 T itk::OrientationDistributionFunction<T, N> 00961 ::GetNematicOrderParameter() const 00962 { 00963 // not yet implemented 00964 return 0; 00965 } 00966 00967 /* 00968 * Calculate StdDev by MaxValue 00969 */ 00970 template < typename T, unsigned int N > 00971 T itk::OrientationDistributionFunction<T, N> 00972 ::GetStdDevByMaxValue() const 00973 { 00974 double mean = 0; 00975 double std = 0; 00976 T max = NumericTraits<T>::NonpositiveMin(); 00977 for( unsigned int i=0; i<InternalDimension; i++) 00978 { 00979 T val = (*this)[i]; 00980 mean += val; 00981 max = (*this)[i] > max ? (*this)[i] : max; 00982 } 00983 mean /= InternalDimension; 00984 for( unsigned int i=0; i<InternalDimension; i++) 00985 { 00986 T val = (*this)[i]; 00987 std += (val - mean) * (val - mean); 00988 } 00989 std /= InternalDimension-1; 00990 00991 if(max == 0) 00992 { 00993 return 0; 00994 } 00995 else 00996 { 00997 return (sqrt(std)/max); 00998 } 00999 } 01000 01001 template < typename T, unsigned int N > 01002 T itk::OrientationDistributionFunction<T, N> 01003 ::GetPrincipleCurvature(double alphaMinDegree, double alphaMaxDegree, int invert) const 01004 { 01005 // following loop only performed once 01006 // (computing indices of each angular range) 01007 m_MutexAngularRange.Lock(); 01008 if(m_AngularRangeIdxs == NULL) 01009 { 01010 m_AngularRangeIdxs = new std::vector< std::vector<int>* >(); 01011 for(unsigned int i=0; i<N; i++) 01012 { 01013 vnl_vector_fixed<double,3> pDir = GetDirection(i); 01014 std::vector<int> *idxs = new std::vector<int>(); 01015 for(unsigned int j=0; j<N; j++) 01016 { 01017 vnl_vector_fixed<double,3> cDir = GetDirection(j); 01018 double angle = ( 180 / ODF_PI ) * acos( dot_product(pDir, cDir) ); 01019 if( (angle < alphaMaxDegree) && (angle > alphaMinDegree) ) 01020 { 01021 idxs->push_back(j); 01022 } 01023 } 01024 m_AngularRangeIdxs->push_back(idxs); 01025 } 01026 } 01027 m_MutexAngularRange.Unlock(); 01028 01029 // find the maximum (or minimum) direction (remember index and value) 01030 T mode; 01031 int pIdx = -1; 01032 if(invert == 0) 01033 { 01034 pIdx = GetPrincipleDiffusionDirection(); 01035 mode = (*this)[pIdx]; 01036 } 01037 else 01038 { 01039 mode = NumericTraits<T>::max(); 01040 for( unsigned int i=0; i<N; i++) 01041 { 01042 if((*this)[i] < mode ) 01043 { 01044 mode = (*this)[i]; 01045 pIdx = i; 01046 } 01047 } 01048 } 01053 01063 01066 01069 01071 // computing a quantile of the angular range 01073 01074 // define quantile 01075 double quantile = 0.00; 01076 01077 // collect all values in angular range of mode 01078 std::vector<T> odfValuesInAngularRange; 01079 int numInRange = m_AngularRangeIdxs->at(pIdx)->size(); 01080 for(int i=0; i<numInRange; i++) 01081 { 01082 odfValuesInAngularRange.push_back((*this)[(*m_AngularRangeIdxs->at(pIdx))[i] ]); 01083 } 01084 01085 // sort them by value 01086 std::sort( odfValuesInAngularRange.begin(), odfValuesInAngularRange.end() ); 01087 01088 // median of angular range 01089 T median = odfValuesInAngularRange[floor(quantile*(double)numInRange+0.5)]; 01090 01091 // compute and return final value 01092 if(mode > median) 01093 { 01094 return mode/median - 1.0; 01095 } 01096 else 01097 { 01098 return median/mode - 1.0; 01099 } 01100 } 01101 01102 /* 01103 * Calculate Normalized Entropy 01104 */ 01105 template < typename T, unsigned int N > 01106 T itk::OrientationDistributionFunction<T, N> 01107 ::GetNormalizedEntropy() const 01108 { 01109 double mean = 0; 01110 for( unsigned int i=0; i<InternalDimension; i++) 01111 { 01112 T val = (*this)[i]; 01113 if( val != 0 ) 01114 { 01115 val = log(val); 01116 } 01117 else 01118 { 01119 val = log(0.0000001); 01120 } 01121 mean += val; 01122 } 01123 double _n = (double) InternalDimension; 01124 mean /= _n; 01125 return (T) (-_n / log(_n) * mean); 01126 } 01127 01128 /* 01129 * Pre-multiply the Tensor by a Matrix 01130 */ 01131 template<class T, unsigned int NOdfDirections> 01132 OrientationDistributionFunction<T, NOdfDirections> 01133 OrientationDistributionFunction<T, NOdfDirections> 01134 ::PreMultiply( const MatrixType & m ) const 01135 { 01136 Self result; 01137 typedef typename NumericTraits<T>::AccumulateType AccumulateType; 01138 for(unsigned int r=0; r<NOdfDirections; r++) 01139 { 01140 for(unsigned int c=0; c<NOdfDirections; c++) 01141 { 01142 AccumulateType sum = NumericTraits<AccumulateType>::ZeroValue(); 01143 for(unsigned int t=0; t<NOdfDirections; t++) 01144 { 01145 sum += m(r,t) * (*this)(t,c); 01146 } 01147 result(r,c) = static_cast<T>( sum ); 01148 } 01149 } 01150 return result; 01151 } 01152 01153 /* 01154 * Post-multiply the Tensor by a Matrix 01155 */ 01156 template<class T, unsigned int NOdfDirections> 01157 OrientationDistributionFunction<T, NOdfDirections> 01158 OrientationDistributionFunction<T, NOdfDirections> 01159 ::PostMultiply( const MatrixType & m ) const 01160 { 01161 Self result; 01162 typedef typename NumericTraits<T>::AccumulateType AccumulateType; 01163 for(unsigned int r=0; r<NOdfDirections; r++) 01164 { 01165 for(unsigned int c=0; c<NOdfDirections; c++) 01166 { 01167 AccumulateType sum = NumericTraits<AccumulateType>::ZeroValue(); 01168 for(unsigned int t=0; t<NOdfDirections; t++) 01169 { 01170 sum += (*this)(r,t) * m(t,c); 01171 } 01172 result(r,c) = static_cast<T>( sum ); 01173 } 01174 } 01175 return result; 01176 } 01177 01178 01179 01180 01184 template<class T, unsigned int NOdfDirections> 01185 std::ostream & 01186 operator<<(std::ostream& os,const OrientationDistributionFunction<T, NOdfDirections> & c ) 01187 { 01188 for(unsigned int i=0; i<c.GetNumberOfComponents(); i++) 01189 { 01190 os << static_cast<typename NumericTraits<T>::PrintType>(c[i]) << " "; 01191 } 01192 return os; 01193 } 01194 01195 01199 template<class T, unsigned int NOdfDirections> 01200 std::istream & 01201 operator>>(std::istream& is, OrientationDistributionFunction<T, NOdfDirections> & dt ) 01202 { 01203 for(unsigned int i=0; i < dt.GetNumberOfComponents(); i++) 01204 { 01205 is >> dt[i]; 01206 } 01207 return is; 01208 } 01209 01210 01211 01212 } // end namespace itk 01213 01214 #endif
