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 ===================================================================*/
16 #include "itkStochasticTractographyFilter.h"
17 #include "vnl/vnl_math.h"
18 #include "vnl/vnl_matrix_fixed.h"
19 #include "vnl/vnl_vector_fixed.h"
20 #include "vnl/vnl_matrix.h"
21 #include "vnl/vnl_sym_matrix.h"
22 #include "vnl/vnl_vector.h"
23 #include "vnl/vnl_diag_matrix.h"
24 #include "vnl/algo/vnl_qr.h"
25 //#include "vnl/algo/vnl_svd.h"
26 #include "vnl/algo/vnl_matrix_inverse.h"
27 //#include "vnl/algo/vnl_symmetric_eigensystem.h"
28 #include "itkSymmetricEigenAnalysis.h"
29 #include "vnl/vnl_transpose.h"
30 #include "itkVariableSizeMatrix.h"
31 #include "itkPathIterator.h"
32 #include "itkImageRegionIterator.h"
33 #include "itkImageRegionConstIterator.h"
34 #include "itkImageRegionConstIteratorWithIndex.h"
38 template< class TInputDWIImage, class TInputWhiteMatterProbabilityImage, class TOutputConnectivityImage >
39 StochasticTractographyFilter< TInputDWIImage, TInputWhiteMatterProbabilityImage, TOutputConnectivityImage >
40 ::StochasticTractographyFilter():
41 m_TotalTracts(0),m_MaxTractLength(0),m_Gradients(NULL), m_TransformedGradients(NULL),m_bValues(NULL),
42 m_SampleDirections(NULL), m_A(NULL), m_Aqr(NULL), m_LikelihoodCachePtr(NULL),
43 m_MaxLikelihoodCacheSize(0), m_CurrentLikelihoodCacheElements(0),
44 m_ClockPtr(NULL), m_TotalDelegatedTracts(0), m_OutputTractContainer(NULL){
45 this->m_SeedIndex[0]=0;
46 this->m_SeedIndex[1]=0;
47 this->m_SeedIndex[2]=0;
48 this->m_MeasurementFrame.set_identity();
49 this->SetNumberOfRequiredInputs(1); //Filter needs a DWI image and a Mask Image, but check will be outside itkProcessObject and inside GenerateData()
52 m_ClockPtr = RealTimeClock::New();
53 this->m_RandomGenerator.reseed( ((unsigned long) this->m_ClockPtr->GetTimeInSeconds()) );
54 //load in default sample directions
55 this->LoadDefaultSampleDirections();
58 template< class TInputDWIImage, class TInputWhiteMatterProbabilityImage, class TOutputConnectivityImage >
59 StochasticTractographyFilter< TInputDWIImage, TInputWhiteMatterProbabilityImage, TOutputConnectivityImage >
60 ::~StochasticTractographyFilter(){
65 template< class TInputDWIImage, class TInputWhiteMatterProbabilityImage, class TOutputConnectivityImage >
67 StochasticTractographyFilter< TInputDWIImage, TInputWhiteMatterProbabilityImage, TOutputConnectivityImage >
68 ::ProbabilisticallyInterpolate( vnl_random& randomgenerator,
69 const TractType::ContinuousIndexType& cindex,
70 typename InputDWIImageType::IndexType& index){
72 for(int i=0; i<3; i++){
73 if ((vcl_ceil(cindex[i]+vnl_math::eps)-cindex[i]) < randomgenerator.drand64())
74 index[i]=(int)vcl_ceil(cindex[i]);
75 else index[i]=(int)vcl_floor(cindex[i]);
79 template< class TInputDWIImage, class TInputWhiteMatterProbabilityImage, class TOutputConnectivityImage >
81 StochasticTractographyFilter< TInputDWIImage, TInputWhiteMatterProbabilityImage, TOutputConnectivityImage >
82 ::UpdateGradientDirections(void){
83 //the gradient direction is transformed into IJK space
84 //by moving into the image space and then to IJK space
86 this->m_TransformedGradients = GradientDirectionContainerType::New();
87 unsigned int N = this->m_Gradients->Size();
88 for(unsigned int i=0; i<N; i++){
89 GradientDirectionContainerType::Element g_i =
90 this->m_MeasurementFrame *
91 this->m_Gradients->GetElement(i);
93 /** The correction to LPS space is not neccessary as of itk 3.2 **/
96 g_i = this->GetInput()->GetDirection().GetInverse() * g_i;
97 this->m_TransformedGradients->InsertElement(i, g_i);
101 template< class TInputDWIImage, class TInputWhiteMatterProbabilityImage, class TOutputConnectivityImage >
103 StochasticTractographyFilter< TInputDWIImage, TInputWhiteMatterProbabilityImage, TOutputConnectivityImage >
104 ::UpdateTensorModelFittingMatrices( void ){
105 //std::cout<<"UpdateTensorFittingMatrix\n";
106 //estimate the parameters using linear LS estimation
107 //using convention specified by Salvador
108 //solve for Beta in: logPhi=X*Beta
109 //number of rows of the matrix depends on the number of inputs,
110 //i.e. the number of measurements of the voxel (n)
111 unsigned int N = this->m_TransformedGradients->Size();
115 this->m_A = new vnl_matrix< double >(N, 7); //potential memory leak here
116 vnl_matrix< double >& A = *(this->m_A);
118 for(unsigned int j=0; j< N ; j++){
119 GradientDirectionContainerType::Element g = m_TransformedGradients->GetElement(j);
120 const bValueType& b_i = m_bValues->GetElement(j);
123 A(j,1)=-1*b_i*(g[0]*g[0]);
124 A(j,2)=-1*b_i*(g[1]*g[1]);
125 A(j,3)=-1*b_i*(g[2]*g[2]);
126 A(j,4)=-1*b_i*(2*g[0]*g[1]);
127 A(j,5)=-1*b_i*(2*g[0]*g[2]);
128 A(j,6)=-1*b_i*(2*g[1]*g[2]);
131 //Store a QR decomposition to quickly estimate
132 //the weighing matrix for each voxel
133 if(this->m_Aqr!=NULL)
135 this->m_Aqr = new vnl_qr< double >(A); //potential memory leak here
138 template< class TInputDWIImage, class TInputWhiteMatterProbabilityImage, class TOutputConnectivityImage >
140 StochasticTractographyFilter< TInputDWIImage, TInputWhiteMatterProbabilityImage, TOutputConnectivityImage >
141 ::CalculateTensorModelParameters( const DWIVectorImageType::PixelType& dwivalues,
142 vnl_diag_matrix<double>& W,
143 TensorModelParamType& tensormodelparams){
145 unsigned int N = this->m_TransformedGradients->Size();
147 //setup const references for code clarity
148 const vnl_matrix< double >& A = *(this->m_A);
149 const vnl_qr< double >& Aqr = *(this->m_Aqr);
151 //vnl_vector is used because the itk vector is limited in its methods and does not
152 //contain an internal vnl class like VariableSizematrix
153 //also itk_matrix has methods which are compatible with vnl_vectors
154 vnl_vector< double > logPhi( N );
156 for(unsigned int j=0; j< N ; j++){
157 //fill up the logPhi vector using log(dwi) values
158 logPhi.put(j, vcl_log(static_cast<double>(dwivalues[j]) + vnl_math::eps));
161 /** Find WLS estimate of the parameters of the Tensor model **/
163 // First estimate W by LS estimation of the intensities
164 //vnl_matrix< double > Q = Aqr.Q();
165 //vnl_vector< double > QtB = Aqr.Q().transpose()*logPhi;
166 //vnl_vector< double > QTB = Aqr.QtB(logPhi);
167 //vnl_matrix< double > R = Aqr.R();
168 W = A* vnl_qr< double >(Aqr.R()).solve(Aqr.QtB(logPhi));
169 //W = A * Aqr.solve(logPhi);
170 for(vnl_diag_matrix< double >::iterator i = W.begin();i!=W.end(); i++){
174 // Now solve for parameters using the estimated weighing matrix
175 tensormodelparams = vnl_qr< double >((W*A).transpose()*W*A).solve(
176 (W*A).transpose()*W*logPhi);
178 //tensormodelparams = vnl_qr< double >((W*A)).solve(W*logPhi);
181 template< class TInputDWIImage, class TInputWhiteMatterProbabilityImage, class TOutputConnectivityImage >
183 StochasticTractographyFilter< TInputDWIImage, TInputWhiteMatterProbabilityImage, TOutputConnectivityImage >
184 ::CalculateConstrainedModelParameters( const TensorModelParamType& tensormodelparams,
185 ConstrainedModelParamType& constrainedmodelparams){
187 vnl_sym_matrix< double > D( 3, 0 );
190 //set the tensor model parameters into a Diffusion tensor
191 D(0,0) = tensormodelparams[1];
192 D(0,1) = tensormodelparams[4];
193 D(0,2) = tensormodelparams[5];
194 D(1,0) = tensormodelparams[4];
195 D(1,1) = tensormodelparams[2];
196 D(1,2) = tensormodelparams[6];
197 D(2,0) = tensormodelparams[5];
198 D(2,1) = tensormodelparams[6];
199 D(2,2) = tensormodelparams[3];
201 //pass through the no gradient intensity Z_0 and
202 //calculate alpha, beta and v hat (the eigenvector
203 //associated with the largest eigenvalue)
204 vnl_matrix_fixed< double, 3, 3 > S(0.0);
205 vnl_vector_fixed< double, 3 > Lambda(0.0);
206 SymmetricEigenAnalysis< vnl_sym_matrix< double >,
207 vnl_vector_fixed< double, 3 >, vnl_matrix_fixed< double, 3, 3 > >
209 eigensystem.ComputeEigenValuesAndVectors( D, Lambda, S );
211 //need to take abs to get rid of negative eigenvalues
212 alpha = (vcl_abs(Lambda[0]) + vcl_abs(Lambda[1])) / 2;
213 beta = vcl_abs(Lambda[2]) - alpha;
215 constrainedmodelparams[0] = tensormodelparams[0];
216 constrainedmodelparams[1] = alpha;
217 constrainedmodelparams[2] = beta;
218 constrainedmodelparams[3] = S[2][0];
219 constrainedmodelparams[4] = S[2][1];
220 constrainedmodelparams[5] = S[2][2];
223 template< class TInputDWIImage, class TInputWhiteMatterProbabilityImage, class TOutputConnectivityImage >
225 StochasticTractographyFilter< TInputDWIImage, TInputWhiteMatterProbabilityImage, TOutputConnectivityImage >
226 ::CalculateNoiseFreeDWIFromConstrainedModel( const ConstrainedModelParamType& constrainedmodelparams,
227 DWIVectorImageType::PixelType& noisefreedwi){
229 unsigned int N = this->m_TransformedGradients->Size();
230 const double& z_0 = constrainedmodelparams[0];
231 const double& alpha = constrainedmodelparams[1];
232 const double& beta = constrainedmodelparams[2];
233 TractOrientationContainerType::Element v_hat( constrainedmodelparams[3],
234 constrainedmodelparams[4],
235 constrainedmodelparams[5]);
237 for(unsigned int i=0; i < N ; i++ ){
238 const double& b_i = this->m_bValues->GetElement(i);
239 const GradientDirectionContainerType::Element& g_i =
240 this->m_TransformedGradients->GetElement(i);
242 noisefreedwi.SetElement(i,
243 vcl_exp(z_0-(alpha*b_i+beta*b_i*vnl_math_sqr(dot_product(g_i, v_hat)))));
247 template< class TInputDWIImage, class TInputWhiteMatterProbabilityImage, class TOutputConnectivityImage >
249 StochasticTractographyFilter< TInputDWIImage, TInputWhiteMatterProbabilityImage, TOutputConnectivityImage >
250 ::CalculateResidualVariance( const DWIVectorImageType::PixelType& noisydwi,
251 const DWIVectorImageType::PixelType& noisefreedwi,
252 const vnl_diag_matrix< double >& W,
253 const unsigned int numberofparameters,
254 double& residualvariance){
256 unsigned int N = this->m_TransformedGradients->Size();
260 /** Not sure if we should be taking difference of log or nonlog intensities **/
261 /** residual variance is too low if we take the difference of log intensities **/
262 /** perhaps using WLS will correct this problem **/
263 for(unsigned int i=0; i<N; i++)
264 residualvariance+=vnl_math_sqr(W(i,i) * (vcl_log(noisydwi[i]/noisefreedwi[i])));
265 residualvariance/=(N-numberofparameters);
268 template< class TInputDWIImage, class TInputWhiteMatterProbabilityImage, class TOutputConnectivityImage >
270 StochasticTractographyFilter< TInputDWIImage, TInputWhiteMatterProbabilityImage, TOutputConnectivityImage >
271 ::CalculateLikelihood( const DWIVectorImageType::PixelType &dwipixel,
272 TractOrientationContainerType::ConstPointer orientations,
273 ProbabilityDistributionImageType::PixelType& likelihood){
275 unsigned int N = this->m_TransformedGradients->Size();
276 TensorModelParamType tensorparams( 0.0 );
277 vnl_diag_matrix< double > W(N,0);
278 ConstrainedModelParamType constrainedparams( 0.0 );
279 DWIVectorImageType::PixelType noisefreedwi(N);
280 double residualvariance=0;
281 double jointlikelihood=1;
283 CalculateTensorModelParameters( dwipixel, W, tensorparams );
284 CalculateConstrainedModelParameters( tensorparams, constrainedparams );
285 CalculateNoiseFreeDWIFromConstrainedModel( constrainedparams, noisefreedwi );
286 CalculateResidualVariance( dwipixel, noisefreedwi, W, 6, residualvariance );
288 for(unsigned int i=0; i < orientations->Size(); i++){
289 /** Vary the entry corresponding to the estimated
290 Tract orientation over the selected sample directions,
291 while preserving the best estimate for the other parameters **/
292 TractOrientationContainerType::Element currentdir = orientations->GetElement(i);
294 /** Incorporate the current sample direction into the secondary parameters **/
295 constrainedparams[3]=currentdir[0];
296 constrainedparams[4]=currentdir[1];
297 constrainedparams[5]=currentdir[2];
299 /** Obtain the estimated
300 intensity for this choice of Tract direction **/
301 CalculateNoiseFreeDWIFromConstrainedModel(constrainedparams, noisefreedwi);
303 jointlikelihood = 1.0;
304 for(unsigned int j=0; j<N; j++){
305 /** Calculate the likelihood given the residualvariance,
306 estimated intensity and the actual intensity (refer to Friman) **/
308 (noisefreedwi[j]/vcl_sqrt(2*vnl_math::pi*residualvariance))*
309 vcl_exp(-vnl_math_sqr(noisefreedwi[j]*vcl_log(dwipixel[j]/noisefreedwi[j]))/
310 (2*residualvariance));
312 likelihood[i]=jointlikelihood;
316 template< class TInputDWIImage, class TInputWhiteMatterProbabilityImage, class TOutputConnectivityImage >
318 StochasticTractographyFilter< TInputDWIImage, TInputWhiteMatterProbabilityImage, TOutputConnectivityImage >
319 ::CalculatePrior( TractOrientationContainerType::Element v_prev,
320 TractOrientationContainerType::ConstPointer orientations,
321 ProbabilityDistributionImageType::PixelType& prior ){
323 const double gamma = 1;
325 for(unsigned int i=0; i < orientations->Size(); i++){
326 if(v_prev.squared_magnitude()==0){
330 prior[i] = dot_product(orientations->GetElement(i),v_prev);;
335 prior[i]=vcl_pow(prior[i],gamma);
341 template< class TInputDWIImage, class TInputWhiteMatterProbabilityImage, class TOutputConnectivityImage >
343 StochasticTractographyFilter< TInputDWIImage, TInputWhiteMatterProbabilityImage, TOutputConnectivityImage >
344 ::CalculatePosterior( const ProbabilityDistributionImageType::PixelType& likelihood,
345 const ProbabilityDistributionImageType::PixelType& prior,
346 ProbabilityDistributionImageType::PixelType& posterior){
349 for(unsigned int i=0; i<likelihood.Size(); i++){
350 sum+=likelihood[i]*prior[i];
352 for(unsigned int i=0; i<likelihood.Size(); i++){
353 posterior[i] = (likelihood[i]*prior[i])/sum;
357 template< class TInputDWIImage, class TInputWhiteMatterProbabilityImage, class TOutputConnectivityImage >
359 StochasticTractographyFilter< TInputDWIImage, TInputWhiteMatterProbabilityImage, TOutputConnectivityImage >
360 ::SampleTractOrientation( vnl_random& randomgenerator,
361 const ProbabilityDistributionImageType::PixelType& posterior,
362 TractOrientationContainerType::ConstPointer orientations,
363 TractOrientationContainerType::Element& choosendirection ){
365 double randomnum = randomgenerator.drand64();
369 //will crash in the unlikely case that 0 was choosen as the randomnum
370 while(cumsum < randomnum){
371 cumsum+=posterior[i];
374 choosendirection = orientations->GetElement(i-1);
376 //std::cout<< "cumsum: " << cumsum<<std::endl;
377 //std::cout<<"selected orientation:( " << (i-1)
378 // <<") "<<choosendirection<< std::endl;
380 template< class TInputDWIImage, class TInputWhiteMatterProbabilityImage, class TOutputConnectivityImage >
382 StochasticTractographyFilter< TInputDWIImage, TInputWhiteMatterProbabilityImage, TOutputConnectivityImage >
383 ::FiberExistenceTest( vnl_random& randomgenerator,
384 typename InputWhiteMatterProbabilityImageType::ConstPointer wmpimage,
385 typename InputWhiteMatterProbabilityImageType::IndexType index ){
386 double randomnum = randomgenerator.drand64();
387 if( randomnum < wmpimage->GetPixel( index ) )
392 //the seedindex is in continuous IJK coordinates
393 template< class TInputDWIImage, class TInputWhiteMatterProbabilityImage, class TOutputConnectivityImage >
395 StochasticTractographyFilter< TInputDWIImage, TInputWhiteMatterProbabilityImage, TOutputConnectivityImage >
396 ::StochasticTractGeneration( typename InputDWIImageType::ConstPointer dwiimagePtr,
397 typename InputWhiteMatterProbabilityImageType::ConstPointer wmpimagePtr,
398 typename InputDWIImageType::IndexType seedindex,
399 unsigned long randomseed,
400 TractType::Pointer tract){
402 TractType::ContinuousIndexType cindex_curr = seedindex;
403 typename InputDWIImageType::IndexType index_curr = {{0,0,0}};
404 ProbabilityDistributionImageType::PixelType
405 prior_curr(this->m_SampleDirections->Size());
406 ProbabilityDistributionImageType::PixelType
407 posterior_curr(this->m_SampleDirections->Size());
408 TractOrientationContainerType::Element v_curr(0,0,0);
409 TractOrientationContainerType::Element v_prev(0,0,0);
412 vnl_random randomgenerator(randomseed);
413 //std::cout<<randomseed<<std::endl;
415 for(unsigned int j=0; j<this->m_MaxTractLength; j++){
416 this->ProbabilisticallyInterpolate( randomgenerator, cindex_curr, index_curr );
418 if(!dwiimagePtr->GetLargestPossibleRegion().IsInside(index_curr)){
422 if( FiberExistenceTest( randomgenerator, wmpimagePtr, index_curr ) ){
423 tract->AddVertex(cindex_curr);
425 this->CalculatePrior( v_prev, this->m_SampleDirections, prior_curr);
427 const ProbabilityDistributionImageType::PixelType&
428 cachelikelihood_curr = this->AccessLikelihoodCache(index_curr);
430 if( cachelikelihood_curr.GetSize() != 0){
431 //use the cached direction
432 this->CalculatePosterior( cachelikelihood_curr, prior_curr, posterior_curr);
435 //do the likelihood calculation and discard
436 //std::cout<<"Cache Miss!\n";
437 ProbabilityDistributionImageType::PixelType
438 likelihood_curr_temp(this->m_SampleDirections->Size());
440 this->CalculateLikelihood(static_cast< DWIVectorImageType::PixelType >(
441 dwiimagePtr->GetPixel(index_curr)),
442 this->m_SampleDirections,
443 likelihood_curr_temp);
444 this->CalculatePosterior( likelihood_curr_temp, prior_curr, posterior_curr);
446 this->SampleTractOrientation(randomgenerator, posterior_curr,
447 this->m_SampleDirections, v_curr);
449 //takes into account voxels of different sizes
450 //converts from a step length of 1 mm to the corresponding length in IJK space
451 const typename InputDWIImageType::SpacingType& spacing = dwiimagePtr->GetSpacing();
452 cindex_curr[0]+=v_curr[0]/spacing[0];
453 cindex_curr[1]+=v_curr[1]/spacing[1];
454 cindex_curr[2]+=v_curr[2]/spacing[2];
458 //fiber doesn't exist in this voxel
459 //std::cout<<"Stopped Tracking: No Fiber in this Voxel\n";
465 template< class TInputDWIImage, class TInputWhiteMatterProbabilityImage, class TOutputConnectivityImage >
467 StochasticTractographyFilter< TInputDWIImage, TInputWhiteMatterProbabilityImage, TOutputConnectivityImage >
468 ::GenerateTractContainerOutput(){
469 //allocate tractcontainer
470 this->m_OutputTractContainer = TractContainerType::New();
472 this->UpdateGradientDirections();
473 this->UpdateTensorModelFittingMatrices();
474 this->m_TotalDelegatedTracts = 0;
476 //calculate the number of voxels to cache from Megabyte memory size limit
477 ProbabilityDistributionImageType::PixelType
478 element(this->GetSampleDirections()->Size());
479 unsigned long elementsize = sizeof(ProbabilityDistributionImageType::PixelType) +
480 sizeof(double)*element.Size();
481 this->m_MaxLikelihoodCacheElements =
482 (this->m_MaxLikelihoodCacheSize*1048576)/elementsize;
483 std::cout << "MaxLikelhoodCacheElements: "
484 << this->m_MaxLikelihoodCacheElements
487 //setup the multithreader
488 StochasticTractGenerationCallbackStruct data;
490 this->GetMultiThreader()->SetSingleMethod( StochasticTractGenerationCallback,
492 this->GetMultiThreader()->SetNumberOfThreads(this->GetNumberOfThreads());
493 std::cout<<"Number of Threads: " << this->GetMultiThreader()->GetNumberOfThreads() << std::endl;
494 //start the multithreaded execution
495 this->GetMultiThreader()->SingleMethodExecute();
496 std::cout<< "CurrentLikelihoodCacheElements: " <<
497 this->m_CurrentLikelihoodCacheElements << std::endl;
500 template< class TInputDWIImage, class TInputWhiteMatterProbabilityImage, class TOutputConnectivityImage >
502 StochasticTractographyFilter< TInputDWIImage, TInputWhiteMatterProbabilityImage, TOutputConnectivityImage >
504 //Generate the tracts
505 this->GenerateTractContainerOutput();
508 this->AllocateOutputs();
510 //write tracts to output image
511 this->TractContainerToConnectivityMap(this->m_OutputTractContainer);
515 template< class TInputDWIImage, class TInputWhiteMatterProbabilityImage, class TOutputConnectivityImage >
516 ITK_THREAD_RETURN_TYPE
517 StochasticTractographyFilter< TInputDWIImage, TInputWhiteMatterProbabilityImage, TOutputConnectivityImage >
518 ::StochasticTractGenerationCallback( void *arg )
520 StochasticTractGenerationCallbackStruct* str=
521 (StochasticTractGenerationCallbackStruct *)
522 (((MultiThreader::ThreadInfoStruct *)(arg))->UserData);
524 typename InputDWIImageType::ConstPointer inputDWIImagePtr = str->Filter->GetInput();
525 typename InputWhiteMatterProbabilityImageType::ConstPointer inputWMPImage =
526 str->Filter->GetWhiteMatterProbabilityImage();
528 unsigned long randomseed=0;
530 while(str->Filter->DelegateTract(randomseed)){
531 //std::cout<<randomseed<<std::endl;
533 TractType::Pointer tract = TractType::New();
535 str->Filter->StochasticTractGeneration( inputDWIImagePtr,
537 str->Filter->GetSeedIndex(),
541 //only store tract if it is of nonzero length
542 if( tract->GetVertexList()->Size() > 4 ){
543 //std::cout<<"Storing tract\n";
544 str->Filter->StoreTract(tract);
547 //std::cout<<"Not Storing Tract\n";
550 return ITK_THREAD_RETURN_VALUE;
553 template< class TInputDWIImage, class TInputWhiteMatterProbabilityImage, class TOutputConnectivityImage >
554 ProbabilityDistributionImageType::PixelType&
555 StochasticTractographyFilter< TInputDWIImage, TInputWhiteMatterProbabilityImage, TOutputConnectivityImage >
556 ::AccessLikelihoodCache( typename InputDWIImageType::IndexType index )
558 this->m_LikelihoodCacheMutexImagePtr->GetPixel(index).Lock();
560 ProbabilityDistributionImageType::PixelType& likelihood =
561 m_LikelihoodCachePtr->GetPixel( index );
562 typename InputDWIImageType::ConstPointer inputDWIImagePtr = this->GetInput();
564 if( likelihood.GetSize() !=0){
565 //entry found in cache
566 this->m_LikelihoodCacheMutexImagePtr->GetPixel(index).Unlock();
569 //we need to lock m_CurrentLikelihoodCacheElements as well but not crucial right now
570 else if( this->m_CurrentLikelihoodCacheElements < this->m_MaxLikelihoodCacheElements ){
571 //entry not found in cache but we have space to store it
572 likelihood.SetSize(this->m_SampleDirections->Size());
574 this->CalculateLikelihood(static_cast< DWIVectorImageType::PixelType >(
575 inputDWIImagePtr->GetPixel(index)),
576 this->m_SampleDirections,
578 this->m_CurrentLikelihoodCacheElements++;
580 this->m_LikelihoodCacheMutexImagePtr->GetPixel(index).Unlock();
584 //entry not found in cache and no space to store it
585 this->m_LikelihoodCacheMutexImagePtr->GetPixel(index).Unlock();
590 this->m_LikelihoodCacheMutexImagePtr->GetPixel(index).Unlock();
597 template< class TInputDWIImage, class TInputWhiteMatterProbabilityImage, class TOutputConnectivityImage >
599 StochasticTractographyFilter< TInputDWIImage, TInputWhiteMatterProbabilityImage, TOutputConnectivityImage >
600 ::DelegateTract(unsigned long& randomseed){
601 bool success = false;
602 this->m_TotalDelegatedTractsMutex.Lock();
603 if(this->m_TotalDelegatedTracts < this->m_TotalTracts){
604 randomseed = this->m_RandomGenerator.lrand32();
605 this->m_TotalDelegatedTracts++;
607 //a tract was successfully delegated
609 else success = false; //all tracts have been delegated
610 this->m_TotalDelegatedTractsMutex.Unlock();
615 template< class TInputDWIImage, class TInputWhiteMatterProbabilityImage, class TOutputConnectivityImage >
617 StochasticTractographyFilter< TInputDWIImage, TInputWhiteMatterProbabilityImage, TOutputConnectivityImage >
618 ::TractContainerToConnectivityMap(TractContainerType::Pointer tractcontainer){
619 //zero the output image
620 typename OutputConnectivityImageType::Pointer outputPtr = this->GetOutput();
621 outputPtr->FillBuffer(0);
623 typedef PathIterator< OutputConnectivityImageType, TractType > OutputTractIteratorType;
625 for(unsigned int i=0; i<tractcontainer->Size(); i++ ){
626 TractType::Pointer tract = tractcontainer->GetElement(i);
627 //std::cout<< tract->EndOfInput() <<std::endl;
628 OutputTractIteratorType outputtractIt( outputPtr,
631 for(outputtractIt.GoToBegin(); !outputtractIt.IsAtEnd(); ++outputtractIt){
632 /* there is an issue using outputtractIt.Value() */
633 // outputtractIt.Set(outputtractIt.Get()+1);
640 template< class TInputDWIImage, class TInputWhiteMatterProbabilityImage, class TOutputConnectivityImage >
642 StochasticTractographyFilter< TInputDWIImage, TInputWhiteMatterProbabilityImage, TOutputConnectivityImage >
643 ::StoreTract(TractType::Pointer tract){
644 this->m_OutputTractContainerMutex.Lock();
645 this->m_OutputTractContainer->InsertElement(
646 this->m_OutputTractContainer->Size(),
648 this->m_OutputTractContainerMutex.Unlock();
651 template< class TInputDWIImage, class TInputWhiteMatterProbabilityImage, class TOutputConnectivityImage >
653 StochasticTractographyFilter< TInputDWIImage, TInputWhiteMatterProbabilityImage, TOutputConnectivityImage >
654 ::GenerateTensorImageOutput(void){
655 this->UpdateGradientDirections();
656 this->UpdateTensorModelFittingMatrices();
658 //allocate the tensor image
659 this->m_OutputTensorImage = OutputTensorImageType::New();
660 m_OutputTensorImage->CopyInformation( this->GetInput() );
661 m_OutputTensorImage->SetBufferedRegion( this->GetInput()->GetBufferedRegion() );
662 m_OutputTensorImage->SetRequestedRegion( this->GetInput()->GetRequestedRegion() );
663 m_OutputTensorImage->Allocate();
665 //define an iterator for the input and output images
666 typedef itk::ImageRegionConstIterator< InputDWIImageType > DWIImageIteratorType;
667 typedef itk::ImageRegionIterator< OutputTensorImageType > TensorImageIteratorType;
670 inputDWIit( this->GetInput(), m_OutputTensorImage->GetRequestedRegion() );
672 TensorImageIteratorType outputtensorit
673 ( m_OutputTensorImage, m_OutputTensorImage->GetRequestedRegion() );
675 unsigned int N = this->m_TransformedGradients->Size();
676 TensorModelParamType tensormodelparams( 0.0 );
677 vnl_diag_matrix< double > W(N,0);
679 for(inputDWIit.GoToBegin(), outputtensorit.GoToBegin();
680 !outputtensorit.IsAtEnd(); ++inputDWIit, ++outputtensorit){
681 CalculateTensorModelParameters( inputDWIit.Get(),
682 W, tensormodelparams);
684 OutputTensorImageType::PixelType& D = outputtensorit.Value();
685 //set the tensor model parameters into a Diffusion tensor
686 D(0,0) = tensormodelparams[1];
687 D(0,1) = tensormodelparams[4];
688 D(0,2) = tensormodelparams[5];
689 D(1,0) = tensormodelparams[4];
690 D(1,1) = tensormodelparams[2];
691 D(1,2) = tensormodelparams[6];
692 D(2,0) = tensormodelparams[5];
693 D(2,1) = tensormodelparams[6];
694 D(2,2) = tensormodelparams[3];