00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017 #ifndef __itkAnalyticalDiffusionQballReconstructionImageFilter_cpp
00018 #define __itkAnalyticalDiffusionQballReconstructionImageFilter_cpp
00019
00020 #include "itkAnalyticalDiffusionQballReconstructionImageFilter.h"
00021 #include "itkImageRegionConstIterator.h"
00022 #include "itkImageRegionConstIteratorWithIndex.h"
00023 #include "itkImageRegionIterator.h"
00024 #include "itkArray.h"
00025 #include "vnl/vnl_vector.h"
00026
00027 #include <boost/version.hpp>
00028
00029 #if BOOST_VERSION / 100000 > 0
00030 #if BOOST_VERSION / 100 % 1000 > 34
00031 #include <boost/math/special_functions/legendre.hpp>
00032 #endif
00033 #endif
00034
00035 #include "itkPointShell.h"
00036
00037 namespace itk {
00038
00039 #define QBALL_ANAL_RECON_PI 3.14159265358979323846
00040
00041 template< class T, class TG, class TO, int L, int NODF>
00042 AnalyticalDiffusionQballReconstructionImageFilter<T,TG,TO,L,NODF>
00043 ::AnalyticalDiffusionQballReconstructionImageFilter() :
00044 m_GradientDirectionContainer(NULL),
00045 m_NumberOfGradientDirections(0),
00046 m_NumberOfBaselineImages(1),
00047 m_Threshold(NumericTraits< ReferencePixelType >::NonpositiveMin()),
00048 m_BValue(1.0),
00049 m_Lambda(0.0),
00050 m_DirectionsDuplicated(false)
00051 {
00052
00053
00054 this->SetNumberOfRequiredInputs( 1 );
00055 }
00056
00057
00058 template<
00059 class TReferenceImagePixelType,
00060 class TGradientImagePixelType,
00061 class TOdfPixelType,
00062 int NOrderL,
00063 int NrOdfDirections>
00064 typename itk::AnalyticalDiffusionQballReconstructionImageFilter<
00065 TReferenceImagePixelType,TGradientImagePixelType,TOdfPixelType,
00066 NOrderL,NrOdfDirections>::OdfPixelType
00067 itk::AnalyticalDiffusionQballReconstructionImageFilter
00068 <TReferenceImagePixelType, TGradientImagePixelType, TOdfPixelType,
00069 NOrderL, NrOdfDirections>
00070 ::Normalize( OdfPixelType odf,
00071 typename NumericTraits<ReferencePixelType>::AccumulateType b0 )
00072 {
00073 switch( m_NormalizationMethod )
00074 {
00075 case QBAR_STANDARD:
00076 {
00077 odf.Normalize();
00078 return odf;
00079 break;
00080 }
00081 case QBAR_B_ZERO_B_VALUE:
00082 {
00083 for(int i=0; i<NrOdfDirections; i++)
00084 {
00085 odf[i] = ((TOdfPixelType)log((TOdfPixelType)b0)-odf[i])/m_BValue;
00086 }
00087 return odf;
00088 break;
00089 }
00090 case QBAR_B_ZERO:
00091 {
00092 odf *= 1.0/b0;
00093 return odf;
00094 break;
00095 }
00096 case QBAR_NONE:
00097 {
00098 return odf;
00099 break;
00100 }
00101 case QBAR_ADC_ONLY:
00102 {
00103 for(int i=0; i<NrOdfDirections; i++)
00104 {
00105 odf[i] = ((TOdfPixelType)log((TOdfPixelType)b0)-odf[i])/m_BValue;
00106 }
00107 return odf;
00108 break;
00109 }
00110 case QBAR_RAW_SIGNAL:
00111 {
00112 return odf;
00113 break;
00114 }
00115 case QBAR_SOLID_ANGLE:
00116 {
00117 for(int i=0; i<NrOdfDirections; i++)
00118 {
00119 odf[i] = odf[i] < 0 ? 0 : odf[i];
00120 }
00121 break;
00122 }
00123 case QBAR_NONNEG_SOLID_ANGLE:
00124 {
00125 break;
00126 }
00127 }
00128
00129 return odf;
00130 }
00131
00132
00133 template<
00134 class TReferenceImagePixelType,
00135 class TGradientImagePixelType,
00136 class TOdfPixelType,
00137 int NOrderL,
00138 int NrOdfDirections>
00139 vnl_vector<TOdfPixelType>
00140 itk::AnalyticalDiffusionQballReconstructionImageFilter
00141 <TReferenceImagePixelType, TGradientImagePixelType, TOdfPixelType,
00142 NOrderL, NrOdfDirections>
00143 ::PreNormalize( vnl_vector<TOdfPixelType> vec,
00144 typename NumericTraits<ReferencePixelType>::AccumulateType b0 )
00145 {
00146 switch( m_NormalizationMethod )
00147 {
00148 case QBAR_STANDARD:
00149 {
00150 return vec;
00151 break;
00152 }
00153 case QBAR_B_ZERO_B_VALUE:
00154 {
00155 int n = vec.size();
00156 for(int i=0; i<n; i++)
00157 {
00158 vec[i] = log(vec[i]);
00159 }
00160 return vec;
00161 break;
00162 }
00163 case QBAR_B_ZERO:
00164 {
00165 return vec;
00166 break;
00167 }
00168 case QBAR_NONE:
00169 {
00170 return vec;
00171 break;
00172 }
00173 case QBAR_ADC_ONLY:
00174 {
00175 int n = vec.size();
00176 for(int i=0; i<n; i++)
00177 {
00178 vec[i] = log(vec[i]);
00179 }
00180 return vec;
00181 break;
00182 }
00183 case QBAR_RAW_SIGNAL:
00184 {
00185 return vec;
00186 break;
00187 }
00188 case QBAR_SOLID_ANGLE:
00189 case QBAR_NONNEG_SOLID_ANGLE:
00190 {
00191 int n = vec.size();
00192 for(int i=0; i<n; i++)
00193 {
00194 float meas = vec[i];
00195 if(meas >= b0)
00196 {
00197 meas = b0 - 0.0001;
00198 }
00199 vec[i] = log(-log(meas/b0));
00200 }
00201 return vec;
00202 break;
00203 }
00204 }
00205
00206 return vec;
00207 }
00208
00209 template< class T, class TG, class TO, int L, int NODF>
00210 void AnalyticalDiffusionQballReconstructionImageFilter<T,TG,TO,L,NODF>
00211 ::BeforeThreadedGenerateData()
00212 {
00213
00214
00215
00216
00217
00218
00219
00220 if( m_NumberOfGradientDirections < 6 )
00221 {
00222 itkExceptionMacro( << "At least 6 gradient directions are required" );
00223 }
00224
00225
00226 std::string gradientImageClassName(
00227 this->ProcessObject::GetInput(0)->GetNameOfClass());
00228 if ( strcmp(gradientImageClassName.c_str(),"VectorImage") != 0 )
00229 {
00230 itkExceptionMacro( <<
00231 "There is only one Gradient image. I expect that to be a VectorImage. "
00232 << "But its of type: " << gradientImageClassName );
00233 }
00234
00235 this->ComputeReconstructionMatrix();
00236
00237 m_BZeroImage = BZeroImageType::New();
00238 typename GradientImagesType::Pointer img = static_cast< GradientImagesType * >(
00239 this->ProcessObject::GetInput(0) );
00240 m_BZeroImage->SetSpacing( img->GetSpacing() );
00241 m_BZeroImage->SetOrigin( img->GetOrigin() );
00242 m_BZeroImage->SetDirection( img->GetDirection() );
00243 m_BZeroImage->SetLargestPossibleRegion( img->GetLargestPossibleRegion());
00244 m_BZeroImage->SetBufferedRegion( img->GetLargestPossibleRegion() );
00245 m_BZeroImage->Allocate();
00246
00247 if(m_NormalizationMethod == QBAR_SOLID_ANGLE ||
00248 m_NormalizationMethod == QBAR_NONNEG_SOLID_ANGLE)
00249 {
00250 m_Lambda = 0.0;
00251 }
00252
00253 }
00254
00255 template< class T, class TG, class TO, int L, int NODF>
00256 void AnalyticalDiffusionQballReconstructionImageFilter<T,TG,TO,L,NODF>
00257 ::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread,
00258 int )
00259 {
00260 typename OutputImageType::Pointer outputImage =
00261 static_cast< OutputImageType * >(this->ProcessObject::GetOutput(0));
00262
00263 ImageRegionIterator< OutputImageType > oit(outputImage, outputRegionForThread);
00264 oit.GoToBegin();
00265
00266 ImageRegionIterator< BZeroImageType > oit2(m_BZeroImage, outputRegionForThread);
00267 oit2.GoToBegin();
00268
00269 typedef ImageRegionConstIterator< GradientImagesType > GradientIteratorType;
00270 typedef typename GradientImagesType::PixelType GradientVectorType;
00271 typename GradientImagesType::Pointer gradientImagePointer = NULL;
00272
00273
00274
00275 gradientImagePointer = static_cast< GradientImagesType * >(
00276 this->ProcessObject::GetInput(0) );
00277
00278 GradientIteratorType git(gradientImagePointer, outputRegionForThread );
00279 git.GoToBegin();
00280
00281
00282 std::vector<unsigned int> baselineind;
00283
00284 std::vector<unsigned int> gradientind;
00285
00286
00287 for(GradientDirectionContainerType::ConstIterator gdcit = this->m_GradientDirectionContainer->Begin();
00288 gdcit != this->m_GradientDirectionContainer->End(); ++gdcit)
00289 {
00290 if(gdcit.Value().one_norm() <= 0.0)
00291 {
00292 baselineind.push_back(gdcit.Index());
00293 }
00294 else
00295 {
00296 gradientind.push_back(gdcit.Index());
00297 }
00298 }
00299
00300 if( m_DirectionsDuplicated )
00301 {
00302 int gradIndSize = gradientind.size();
00303 for(int i=0; i<gradIndSize; i++)
00304 gradientind.push_back(gradientind[i]);
00305 }
00306
00307 while( !git.IsAtEnd() )
00308 {
00309 GradientVectorType b = git.Get();
00310
00311 typename NumericTraits<ReferencePixelType>::AccumulateType b0 = NumericTraits<ReferencePixelType>::Zero;
00312
00313
00314 for(unsigned int i = 0; i < baselineind.size(); ++i)
00315 {
00316 b0 += b[baselineind[i]];
00317 }
00318 b0 /= this->m_NumberOfBaselineImages;
00319
00320 OdfPixelType odf(0.0);
00321 vnl_vector<TO> B(m_NumberOfGradientDirections);
00322
00323 if( (b0 != 0) && (b0 >= m_Threshold) )
00324 {
00325 for( unsigned int i = 0; i< m_NumberOfGradientDirections; i++ )
00326 {
00327 B[i] = static_cast<TO>(b[gradientind[i]]);
00328 }
00329
00330 B = PreNormalize(B, b0);
00331 if(m_NormalizationMethod == QBAR_SOLID_ANGLE)
00332 {
00333 vnl_vector<TO> coeffs(m_NumberCoefficients);
00334 coeffs = ( (*m_CoeffReconstructionMatrix) * B );
00335 coeffs[0] += 1.0/(2.0*sqrt(QBALL_ANAL_RECON_PI));
00336 odf = ( (*m_SphericalHarmonicBasisMatrix) * coeffs ).data_block();
00337 }
00338 else if(m_NormalizationMethod == QBAR_NONNEG_SOLID_ANGLE)
00339 {
00347 QuadProgPP::Matrix<double>& G(m_G);
00348 int lenb = B.size();
00349 vnl_vector<double>* s = new vnl_vector<double>(lenb);
00350 for (int ii=0; ii<lenb; ii++) (*s)(ii) = B(ii);
00351 vnl_vector<double> g0_ = -1.0 * (*m_B_t) * (*s);
00352 QuadProgPP::Vector<double> g0(g0_.data_block(),m_NumberCoefficients);
00353 try
00354 {
00355 QuadProgPP::QuadProg::solve_quadprog(G,g0,m_CE,m_ce0,m_CI,m_ci0,m_x);
00356 }
00357 catch(...)
00358 {
00359 m_x = 0;
00360 }
00361 vnl_vector<TO> coeffs(m_NumberCoefficients);
00362 for (int ii=0; ii<m_NumberCoefficients; ii++) coeffs(ii) = (*m_LP)(ii) * m_x[ii];
00363 coeffs[0] += 1.0/(2.0*sqrt(QBALL_ANAL_RECON_PI));
00364 odf = ( (*m_SphericalHarmonicBasisMatrix) * coeffs ).data_block();
00365
00366 }
00367 else
00368 {
00369 odf = ( (*m_ReconstructionMatrix) * B ).data_block();
00370 }
00371 odf = Normalize(odf, b0);
00372
00373 }
00374
00375 oit.Set( odf );
00376 ++oit;
00377 oit2.Set( b0 );
00378 ++oit2;
00379 ++git;
00380 }
00381
00382 std::cout << "One Thread finished reconstruction" << std::endl;
00383 }
00384
00385 template< class T, class TG, class TO, int L, int NODF>
00386 void AnalyticalDiffusionQballReconstructionImageFilter<T,TG,TO,L,NODF>
00387 ::tofile2(vnl_matrix<double> *pA, std::string fname)
00388 {
00389 vnl_matrix<double> A = (*pA);
00390 ofstream myfile;
00391 myfile.open (fname.c_str());
00392 myfile << "A1=[";
00393 for(int i=0; i<A.rows(); i++)
00394 {
00395 for(int j=0; j<A.columns(); j++)
00396 {
00397 myfile << A(i,j) << " ";
00398 if(j==A.columns()-1 && i!=A.rows()-1)
00399 myfile << ";";
00400 }
00401 }
00402 myfile << "];";
00403 myfile.close();
00404
00405 }
00406
00407 template< class T, class TG, class TO, int L, int NODF>
00408 double AnalyticalDiffusionQballReconstructionImageFilter<T,TG,TO,L,NODF>
00409 ::factorial(int number) {
00410 if(number <= 1) return 1;
00411 double result = 1.0;
00412 for(int i=1; i<=number; i++)
00413 result *= i;
00414 return result;
00415 }
00416
00417 template< class T, class TG, class TO, int L, int NODF>
00418 void AnalyticalDiffusionQballReconstructionImageFilter<T,TG,TO,L,NODF>
00419 ::Cart2Sph(double x, double y, double z, double *cart)
00420 {
00421 double phi, th, rad;
00422 rad = sqrt(x*x+y*y+z*z);
00423 th = atan2(z,sqrt(x*x+y*y));
00424 phi = atan2(y,x);
00425 th = -th+QBALL_ANAL_RECON_PI/2;
00426 phi = -phi+QBALL_ANAL_RECON_PI;
00427 cart[0] = phi;
00428 cart[1] = th;
00429 cart[2] = rad;
00430 }
00431
00432 template< class T, class TG, class TO, int L, int NODF>
00433 double AnalyticalDiffusionQballReconstructionImageFilter<T,TG,TO,L,NODF>
00434 ::legendre0(int l)
00435 {
00436 if( l%2 != 0 )
00437 {
00438 return 0;
00439 }
00440 else
00441 {
00442 double prod1 = 1.0;
00443 for(int i=1;i<l;i+=2) prod1 *= i;
00444 double prod2 = 1.0;
00445 for(int i=2;i<=l;i+=2) prod2 *= i;
00446 return pow(-1.0,l/2.0)*(prod1/prod2);
00447 }
00448 }
00449
00450 template< class T, class TG, class TO, int L, int NODF>
00451 double AnalyticalDiffusionQballReconstructionImageFilter<T,TG,TO,L,NODF>
00452 ::spherical_harmonic(int m,int l,double theta,double phi, bool complexPart)
00453 {
00454 if( theta < 0 || theta > QBALL_ANAL_RECON_PI )
00455 {
00456 std::cout << "bad range" << std::endl;
00457 return 0;
00458 }
00459
00460 if( phi < 0 || phi > 2*QBALL_ANAL_RECON_PI )
00461 {
00462 std::cout << "bad range" << std::endl;
00463 return 0;
00464 }
00465
00466 double pml = 0;
00467 double fac1 = factorial(l+m);
00468 double fac2 = factorial(l-m);
00469
00470 if( m<0 )
00471 {
00472 #if BOOST_VERSION / 100000 > 0
00473 #if BOOST_VERSION / 100 % 1000 > 34
00474 pml = ::boost::math::legendre_p<double>(l, -m, cos(theta));
00475 #else
00476 std::cout << "ERROR: Boost 1.35 minimum required" << std::endl;
00477 #endif
00478 #else
00479 std::cout << "ERROR: Boost 1.35 minimum required" << std::endl;
00480 #endif
00481 double mypow = pow(-1.0,-m);
00482 double myfac = (fac1/fac2);
00483 pml *= mypow*myfac;
00484 }
00485 else
00486 {
00487 #if BOOST_VERSION / 100000 > 0
00488 #if BOOST_VERSION / 100 % 1000 > 34
00489 pml = ::boost::math::legendre_p<double>(l, m, cos(theta));
00490 #endif
00491 #endif
00492 }
00493
00494
00495
00496 double retval = sqrt(((2.0*(double)l+1.0)/(4.0*QBALL_ANAL_RECON_PI))*(fac2/fac1)) * pml;
00497 if( !complexPart )
00498 {
00499 retval *= cos(m*phi);
00500 }
00501 else
00502 {
00503 retval *= sin(m*phi);
00504 }
00505
00506 return retval;
00507 }
00508
00509 template< class T, class TG, class TO, int L, int NODF>
00510 double AnalyticalDiffusionQballReconstructionImageFilter<T,TG,TO,L,NODF>
00511 ::Yj(int m, int k, double theta, double phi)
00512 {
00513 if( -k <= m && m < 0 )
00514 {
00515 return sqrt(2.0) * spherical_harmonic(m,k,theta,phi,false);
00516 }
00517
00518 if( m == 0 )
00519 return spherical_harmonic(0,k,theta,phi,false);
00520
00521 if( 0 < m && m <= k )
00522 {
00523 return sqrt(2.0) * spherical_harmonic(m,k,theta,phi,true);
00524 }
00525
00526 return 0;
00527 }
00528
00529
00530
00531 template< class T, class TG, class TO, int L, int NODF>
00532 void AnalyticalDiffusionQballReconstructionImageFilter<T,TG,TO,L,NODF>
00533 ::ComputeReconstructionMatrix()
00534 {
00535
00536
00537
00538
00539 if( m_NumberOfGradientDirections < 6 )
00540 {
00541 itkExceptionMacro( << "Not enough gradient directions supplied. Need to supply at least 6" );
00542 }
00543
00544 {
00545
00546 bool warning = false;
00547 for(GradientDirectionContainerType::ConstIterator gdcit1 = this->m_GradientDirectionContainer->Begin();
00548 gdcit1 != this->m_GradientDirectionContainer->End(); ++gdcit1)
00549 {
00550 for(GradientDirectionContainerType::ConstIterator gdcit2 = this->m_GradientDirectionContainer->Begin();
00551 gdcit2 != this->m_GradientDirectionContainer->End(); ++gdcit2)
00552 {
00553 if(gdcit1.Value() == gdcit2.Value() && gdcit1.Index() != gdcit2.Index())
00554 {
00555 itkWarningMacro( << "Some of the Diffusion Gradients equal each other. Corresponding image data should be averaged before calling this filter." );
00556 warning = true;
00557 break;
00558 }
00559 }
00560 if (warning) break;
00561 }
00562
00563
00564
00565
00566 vnl_vector<double> centerMass(3);
00567 centerMass.fill(0.0);
00568 int count = 0;
00569 for(GradientDirectionContainerType::ConstIterator gdcit1 = this->m_GradientDirectionContainer->Begin();
00570 gdcit1 != this->m_GradientDirectionContainer->End(); ++gdcit1)
00571 {
00572 if(gdcit1.Value().one_norm() > 0.0)
00573 {
00574 centerMass += gdcit1.Value();
00575 count ++;
00576 }
00577 }
00578 centerMass /= count;
00579 if(centerMass.two_norm() > 0.1)
00580 {
00581 m_DirectionsDuplicated = true;
00582 m_NumberOfGradientDirections *= 2;
00583 }
00584 }
00585
00586 vnl_matrix<double> *Q
00587 = new vnl_matrix<double>(3, m_NumberOfGradientDirections);
00588
00589 {
00590 int i = 0;
00591 for(GradientDirectionContainerType::ConstIterator gdcit = this->m_GradientDirectionContainer->Begin();
00592 gdcit != this->m_GradientDirectionContainer->End(); ++gdcit)
00593 {
00594 if(gdcit.Value().one_norm() > 0.0)
00595 {
00596 double x = gdcit.Value().get(0);
00597 double y = gdcit.Value().get(1);
00598 double z = gdcit.Value().get(2);
00599 double cart[3];
00600 Cart2Sph(x,y,z,cart);
00601 (*Q)(0,i) = cart[0];
00602 (*Q)(1,i) = cart[1];
00603 (*Q)(2,i++) = cart[2];
00604 }
00605 }
00606 if(m_DirectionsDuplicated)
00607 {
00608 for(GradientDirectionContainerType::ConstIterator gdcit = this->m_GradientDirectionContainer->Begin();
00609 gdcit != this->m_GradientDirectionContainer->End(); ++gdcit)
00610 {
00611 if(gdcit.Value().one_norm() > 0.0)
00612 {
00613 double x = gdcit.Value().get(0);
00614 double y = gdcit.Value().get(1);
00615 double z = gdcit.Value().get(2);
00616 double cart[3];
00617 Cart2Sph(x,y,z,cart);
00618 (*Q)(0,i) = cart[0];
00619 (*Q)(1,i) = cart[1];
00620 (*Q)(2,i++) = cart[2];
00621 }
00622 }
00623 }
00624 }
00625
00626 int l = L;
00627 m_NumberCoefficients = (int)(l*l + l + 2.0)/2.0 + l;
00628 vnl_matrix<double>* B = new vnl_matrix<double>(m_NumberOfGradientDirections,m_NumberCoefficients);
00629 vnl_matrix<double>* _L = new vnl_matrix<double>(m_NumberCoefficients,m_NumberCoefficients);
00630 _L->fill(0.0);
00631 vnl_matrix<double>* LL = new vnl_matrix<double>(m_NumberCoefficients,m_NumberCoefficients);
00632 LL->fill(0.0);
00633 vnl_matrix<double>* P = new vnl_matrix<double>(m_NumberCoefficients,m_NumberCoefficients);
00634 P->fill(0.0);
00635 vnl_matrix<double>* Inv = new vnl_matrix<double>(m_NumberCoefficients,m_NumberCoefficients);
00636 P->fill(0.0);
00637 vnl_vector<int>* lj = new vnl_vector<int>(m_NumberCoefficients);
00638 m_LP = new vnl_vector<double>(m_NumberCoefficients);
00639
00640 for(unsigned int i=0; i<m_NumberOfGradientDirections; i++)
00641 {
00642 for(int k=0; k<=l; k+=2)
00643 {
00644 for(int m=-k; m<=k; m++)
00645 {
00646 int j = (k*k + k + 2)/2 + m - 1;
00647 (*_L)(j,j) = -k*(k+1);
00648 (*m_LP)(j) = -k*(k+1);
00649 (*lj)[j] = k;
00650 double phi = (*Q)(0,i);
00651 double th = (*Q)(1,i);
00652 (*B)(i,j) = Yj(m,k,th,phi);
00653
00654 }
00655 }
00656 }
00657
00658
00659 LL->update(*_L);
00660 *LL *= *_L;
00661
00662
00663 for(int i=0; i<m_NumberCoefficients; i++)
00664 {
00665
00666
00667
00668
00669 if(m_NormalizationMethod == QBAR_SOLID_ANGLE ||
00670 m_NormalizationMethod == QBAR_NONNEG_SOLID_ANGLE)
00671 {
00672 (*P)(i,i) = 2.0*QBALL_ANAL_RECON_PI*legendre0((*lj)[i]);
00673 (*m_LP)(i) *= (*P)(i,i);
00674 }
00675 else
00676 {
00677 (*P)(i,i) = legendre0((*lj)[i]);
00678 }
00679 }
00680 m_B_t = new vnl_matrix<double>(B->transpose());
00681
00682 vnl_matrix<double> B_t_B = (*m_B_t) * (*B);
00683
00684 vnl_matrix<double> lambdaLL(m_NumberCoefficients,m_NumberCoefficients);
00685 lambdaLL.update((*LL));
00686 lambdaLL *= m_Lambda;
00687
00688
00689 vnl_matrix<double> tmp( B_t_B + lambdaLL);
00690 vnl_matrix_inverse<double> *pseudoInverse
00691 = new vnl_matrix_inverse<double>( tmp );
00692
00693 (*Inv) = pseudoInverse->pinverse();
00694
00695 vnl_matrix<double> temp((*Inv) * (*m_B_t));
00696 double fac1 = (1.0/(16.0*QBALL_ANAL_RECON_PI*QBALL_ANAL_RECON_PI));
00697 switch(m_NormalizationMethod)
00698 {
00699 case QBAR_ADC_ONLY:
00700 case QBAR_RAW_SIGNAL:
00701 break;
00702 case QBAR_STANDARD:
00703 case QBAR_B_ZERO_B_VALUE:
00704 case QBAR_B_ZERO:
00705 case QBAR_NONE:
00706 temp = (*P) * temp;
00707 break;
00708 case QBAR_SOLID_ANGLE:
00709 temp = fac1 * (*P) * (*_L) * temp;
00710 break;
00711 case QBAR_NONNEG_SOLID_ANGLE:
00712 m_G = QuadProgPP::Matrix<double>(B_t_B.data_block(), B_t_B.rows(), B_t_B.cols());
00713 m_CE = QuadProgPP::Matrix<double>((double)0,m_NumberCoefficients,0);
00714 m_ce0 = QuadProgPP::Vector<double>((double)0,0);
00715 m_ci0 = QuadProgPP::Vector<double>(4*QBALL_ANAL_RECON_PI, NODF);
00716 m_x = QuadProgPP::Vector<double>(m_NumberCoefficients);
00717 (*m_LP) *= fac1;
00718 break;
00719 }
00720
00721
00722
00723 m_CoeffReconstructionMatrix = new vnl_matrix<TO>(m_NumberCoefficients,m_NumberOfGradientDirections);
00724 for(int i=0; i<m_NumberCoefficients; i++)
00725 {
00726 for(unsigned int j=0; j<m_NumberOfGradientDirections; j++)
00727 {
00728 (*m_CoeffReconstructionMatrix)(i,j) = (float) temp(i,j);
00729 }
00730 }
00731
00732
00733
00734 int NOdfDirections = NODF;
00735 vnl_matrix_fixed<double, 3, NODF>* U =
00736 itk::PointShell<NODF, vnl_matrix_fixed<double, 3, NODF> >::DistributePointShell();
00737
00738 m_SphericalHarmonicBasisMatrix = new vnl_matrix<TO>(NOdfDirections,m_NumberCoefficients);
00739 vnl_matrix<double>* sphericalHarmonicBasisMatrix2
00740 = new vnl_matrix<double>(NOdfDirections,m_NumberCoefficients);
00741 for(int i=0; i<NOdfDirections; i++)
00742 {
00743 double x = (*U)(0,i);
00744 double y = (*U)(1,i);
00745 double z = (*U)(2,i);
00746 double cart[3];
00747 Cart2Sph(x,y,z,cart);
00748 (*U)(0,i) = cart[0];
00749 (*U)(1,i) = cart[1];
00750 (*U)(2,i) = cart[2];
00751 }
00752
00753 for(int i=0; i<NOdfDirections; i++)
00754 {
00755 for(int k=0; k<=l; k+=2)
00756 {
00757 for(int m=-k; m<=k; m++)
00758 {
00759 int j = (k*k + k + 2)/2 + m - 1;
00760 double phi = (*U)(0,i);
00761 double th = (*U)(1,i);
00762 (*m_SphericalHarmonicBasisMatrix)(i,j) = Yj(m,k,th,phi);
00763 (*sphericalHarmonicBasisMatrix2)(i,j) = (*m_SphericalHarmonicBasisMatrix)(i,j);
00764 }
00765 }
00766 }
00767
00768 if(m_NormalizationMethod == QBAR_NONNEG_SOLID_ANGLE)
00769 {
00770 vnl_matrix<double> CI_t =
00771 (*sphericalHarmonicBasisMatrix2) * (*P) * (*_L);
00772 m_CI = QuadProgPP::Matrix<double>(CI_t.transpose().data_block(), m_NumberCoefficients, NOdfDirections);
00773 }
00774
00775 m_ReconstructionMatrix = new vnl_matrix<TO>(NOdfDirections,m_NumberOfGradientDirections);
00776 *m_ReconstructionMatrix = (*m_SphericalHarmonicBasisMatrix) * (*m_CoeffReconstructionMatrix);
00777
00778 }
00779
00780 template< class T, class TG, class TO, int L, int NODF>
00781 void AnalyticalDiffusionQballReconstructionImageFilter<T,TG,TO,L,NODF>
00782 ::SetGradientImage( GradientDirectionContainerType *gradientDirection,
00783 const GradientImagesType *gradientImage )
00784 {
00785 this->m_GradientDirectionContainer = gradientDirection;
00786
00787 unsigned int numImages = gradientDirection->Size();
00788 this->m_NumberOfBaselineImages = 0;
00789 for(GradientDirectionContainerType::Iterator it = this->m_GradientDirectionContainer->Begin();
00790 it != this->m_GradientDirectionContainer->End(); it++)
00791 {
00792 if(it.Value().one_norm() <= 0.0)
00793 {
00794 this->m_NumberOfBaselineImages++;
00795 }
00796 else
00797 {
00798 it.Value() = it.Value() / it.Value().two_norm();
00799 }
00800 }
00801
00802 this->m_NumberOfGradientDirections = numImages - this->m_NumberOfBaselineImages;
00803
00804
00805
00806 if( gradientImage->GetVectorLength() != this->m_NumberOfBaselineImages + m_NumberOfGradientDirections )
00807 {
00808 itkExceptionMacro( << m_NumberOfGradientDirections << " gradients + " << this->m_NumberOfBaselineImages
00809 << "baselines = " << m_NumberOfGradientDirections + this->m_NumberOfBaselineImages
00810 << " directions specified but image has " << gradientImage->GetVectorLength()
00811 << " components.");
00812 }
00813
00814 this->ProcessObject::SetNthInput( 0,
00815 const_cast< GradientImagesType* >(gradientImage) );
00816
00817 }
00818
00819 template< class T, class TG, class TO, int L, int NODF>
00820 void AnalyticalDiffusionQballReconstructionImageFilter<T,TG,TO,L,NODF>
00821 ::PrintSelf(std::ostream& os, Indent indent) const
00822 {
00823 Superclass::PrintSelf(os,indent);
00824
00825 os << indent << "OdfReconstructionMatrix: " << m_ReconstructionMatrix << std::endl;
00826 if ( m_GradientDirectionContainer )
00827 {
00828 os << indent << "GradientDirectionContainer: "
00829 << m_GradientDirectionContainer << std::endl;
00830 }
00831 else
00832 {
00833 os << indent <<
00834 "GradientDirectionContainer: (Gradient directions not set)" << std::endl;
00835 }
00836 os << indent << "NumberOfGradientDirections: " <<
00837 m_NumberOfGradientDirections << std::endl;
00838 os << indent << "NumberOfBaselineImages: " <<
00839 m_NumberOfBaselineImages << std::endl;
00840 os << indent << "Threshold for reference B0 image: " << m_Threshold << std::endl;
00841 os << indent << "BValue: " << m_BValue << std::endl;
00842 }
00843
00844 }
00845
00846 #endif // __itkAnalyticalDiffusionQballReconstructionImageFilter_cpp