Medical Imaging Interaction Toolkit  2016.11.0
Medical Imaging Interaction Toolkit
itkStochasticTractographyFilter.txx
Go to the documentation of this file.
1 /*===================================================================
2 
3 The Medical Imaging Interaction Toolkit (MITK)
4 
5 Copyright (c) German Cancer Research Center,
6 Division of Medical and Biological Informatics.
7 All rights reserved.
8 
9 This software is distributed WITHOUT ANY WARRANTY; without
10 even the implied warranty of MERCHANTABILITY or FITNESS FOR
11 A PARTICULAR PURPOSE.
12 
13 See LICENSE.txt or http://www.mitk.org for details.
14 
15 ===================================================================*/
16 #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"
35 
36 namespace itk{
37 
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()
50 
51 
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();
56 }
57 
58 template< class TInputDWIImage, class TInputWhiteMatterProbabilityImage, class TOutputConnectivityImage >
59 StochasticTractographyFilter< TInputDWIImage, TInputWhiteMatterProbabilityImage, TOutputConnectivityImage >
60 ::~StochasticTractographyFilter(){
61  delete this->m_A;
62  delete this->m_Aqr;
63 }
64 
65 template< class TInputDWIImage, class TInputWhiteMatterProbabilityImage, class TOutputConnectivityImage >
66 void
67 StochasticTractographyFilter< TInputDWIImage, TInputWhiteMatterProbabilityImage, TOutputConnectivityImage >
68 ::ProbabilisticallyInterpolate( vnl_random& randomgenerator,
69  const TractType::ContinuousIndexType& cindex,
70  typename InputDWIImageType::IndexType& index){
71 
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]);
76  }
77 }
78 
79 template< class TInputDWIImage, class TInputWhiteMatterProbabilityImage, class TOutputConnectivityImage >
80 void
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
85 
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);
92 
93  /** The correction to LPS space is not neccessary as of itk 3.2 **/
94  //g_i[0] = -g_i[0];
95  //g_i[1] = -g_i[1];
96  g_i = this->GetInput()->GetDirection().GetInverse() * g_i;
97  this->m_TransformedGradients->InsertElement(i, g_i);
98  }
99 }
100 
101 template< class TInputDWIImage, class TInputWhiteMatterProbabilityImage, class TOutputConnectivityImage >
102 void
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();
112 
113  if(this->m_A!=NULL)
114  delete this->m_A;
115  this->m_A = new vnl_matrix< double >(N, 7); //potential memory leak here
116  vnl_matrix< double >& A = *(this->m_A);
117 
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);
121 
122  A(j,0)=1.0;
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]);
129  }
130 
131  //Store a QR decomposition to quickly estimate
132  //the weighing matrix for each voxel
133  if(this->m_Aqr!=NULL)
134  delete this->m_Aqr;
135  this->m_Aqr = new vnl_qr< double >(A); //potential memory leak here
136 }
137 
138 template< class TInputDWIImage, class TInputWhiteMatterProbabilityImage, class TOutputConnectivityImage >
139 void
140 StochasticTractographyFilter< TInputDWIImage, TInputWhiteMatterProbabilityImage, TOutputConnectivityImage >
141 ::CalculateTensorModelParameters( const DWIVectorImageType::PixelType& dwivalues,
142  vnl_diag_matrix<double>& W,
143  TensorModelParamType& tensormodelparams){
144 
145  unsigned int N = this->m_TransformedGradients->Size();
146 
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);
150 
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 );
155 
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));
159  }
160 
161  /** Find WLS estimate of the parameters of the Tensor model **/
162 
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++){
171  *i = vcl_exp( *i );
172  }
173 
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);
177  //int a;
178  //tensormodelparams = vnl_qr< double >((W*A)).solve(W*logPhi);
179 }
180 
181 template< class TInputDWIImage, class TInputWhiteMatterProbabilityImage, class TOutputConnectivityImage >
182 void
183 StochasticTractographyFilter< TInputDWIImage, TInputWhiteMatterProbabilityImage, TOutputConnectivityImage >
184 ::CalculateConstrainedModelParameters( const TensorModelParamType& tensormodelparams,
185  ConstrainedModelParamType& constrainedmodelparams){
186 
187  vnl_sym_matrix< double > D( 3, 0 );
188  double alpha =0;
189  double beta=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];
200 
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 > >
208  eigensystem( 3 );
209  eigensystem.ComputeEigenValuesAndVectors( D, Lambda, S );
210 
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;
214 
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];
221 }
222 
223 template< class TInputDWIImage, class TInputWhiteMatterProbabilityImage, class TOutputConnectivityImage >
224 void
225 StochasticTractographyFilter< TInputDWIImage, TInputWhiteMatterProbabilityImage, TOutputConnectivityImage >
226 ::CalculateNoiseFreeDWIFromConstrainedModel( const ConstrainedModelParamType& constrainedmodelparams,
227  DWIVectorImageType::PixelType& noisefreedwi){
228 
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]);
236 
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);
241 
242  noisefreedwi.SetElement(i,
243  vcl_exp(z_0-(alpha*b_i+beta*b_i*vnl_math_sqr(dot_product(g_i, v_hat)))));
244  }
245 }
246 
247 template< class TInputDWIImage, class TInputWhiteMatterProbabilityImage, class TOutputConnectivityImage >
248 void
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){
255 
256  unsigned int N = this->m_TransformedGradients->Size();
257 
258  residualvariance=0;
259 
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);
266 }
267 
268 template< class TInputDWIImage, class TInputWhiteMatterProbabilityImage, class TOutputConnectivityImage >
269 void
270 StochasticTractographyFilter< TInputDWIImage, TInputWhiteMatterProbabilityImage, TOutputConnectivityImage >
271 ::CalculateLikelihood( const DWIVectorImageType::PixelType &dwipixel,
272  TractOrientationContainerType::ConstPointer orientations,
273  ProbabilityDistributionImageType::PixelType& likelihood){
274 
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;
282 
283  CalculateTensorModelParameters( dwipixel, W, tensorparams );
284  CalculateConstrainedModelParameters( tensorparams, constrainedparams );
285  CalculateNoiseFreeDWIFromConstrainedModel( constrainedparams, noisefreedwi );
286  CalculateResidualVariance( dwipixel, noisefreedwi, W, 6, residualvariance );
287 
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);
293 
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];
298 
299  /** Obtain the estimated
300  intensity for this choice of Tract direction **/
301  CalculateNoiseFreeDWIFromConstrainedModel(constrainedparams, noisefreedwi);
302 
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) **/
307  jointlikelihood *=
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));
311  }
312  likelihood[i]=jointlikelihood;
313  }
314 }
315 
316 template< class TInputDWIImage, class TInputWhiteMatterProbabilityImage, class TOutputConnectivityImage >
317 void
318 StochasticTractographyFilter< TInputDWIImage, TInputWhiteMatterProbabilityImage, TOutputConnectivityImage >
319 ::CalculatePrior( TractOrientationContainerType::Element v_prev,
320  TractOrientationContainerType::ConstPointer orientations,
321  ProbabilityDistributionImageType::PixelType& prior ){
322 
323  const double gamma = 1;
324 
325  for(unsigned int i=0; i < orientations->Size(); i++){
326  if(v_prev.squared_magnitude()==0){
327  prior[i]=1.0;
328  }
329  else{
330  prior[i] = dot_product(orientations->GetElement(i),v_prev);;
331  if(prior[i]<0){
332  prior[i]=0;
333  }
334  else{
335  prior[i]=vcl_pow(prior[i],gamma);
336  }
337  }
338  }
339 }
340 
341 template< class TInputDWIImage, class TInputWhiteMatterProbabilityImage, class TOutputConnectivityImage >
342 void
343 StochasticTractographyFilter< TInputDWIImage, TInputWhiteMatterProbabilityImage, TOutputConnectivityImage >
344 ::CalculatePosterior( const ProbabilityDistributionImageType::PixelType& likelihood,
345  const ProbabilityDistributionImageType::PixelType& prior,
346  ProbabilityDistributionImageType::PixelType& posterior){
347 
348  double sum=0;
349  for(unsigned int i=0; i<likelihood.Size(); i++){
350  sum+=likelihood[i]*prior[i];
351  }
352  for(unsigned int i=0; i<likelihood.Size(); i++){
353  posterior[i] = (likelihood[i]*prior[i])/sum;
354  }
355 }
356 
357 template< class TInputDWIImage, class TInputWhiteMatterProbabilityImage, class TOutputConnectivityImage >
358 void
359 StochasticTractographyFilter< TInputDWIImage, TInputWhiteMatterProbabilityImage, TOutputConnectivityImage >
360 ::SampleTractOrientation( vnl_random& randomgenerator,
361  const ProbabilityDistributionImageType::PixelType& posterior,
362  TractOrientationContainerType::ConstPointer orientations,
363  TractOrientationContainerType::Element& choosendirection ){
364 
365  double randomnum = randomgenerator.drand64();
366  int i=0;
367  double cumsum=0;
368 
369  //will crash in the unlikely case that 0 was choosen as the randomnum
370  while(cumsum < randomnum){
371  cumsum+=posterior[i];
372  i++;
373  }
374  choosendirection = orientations->GetElement(i-1);
375 
376  //std::cout<< "cumsum: " << cumsum<<std::endl;
377  //std::cout<<"selected orientation:( " << (i-1)
378  // <<") "<<choosendirection<< std::endl;
379 }
380 template< class TInputDWIImage, class TInputWhiteMatterProbabilityImage, class TOutputConnectivityImage >
381 bool
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 ) )
388  return true;
389  else
390  return false;
391 }
392 //the seedindex is in continuous IJK coordinates
393 template< class TInputDWIImage, class TInputWhiteMatterProbabilityImage, class TOutputConnectivityImage >
394 void
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){
401 
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);
410 
411  tract->Initialize();
412  vnl_random randomgenerator(randomseed);
413  //std::cout<<randomseed<<std::endl;
414 
415  for(unsigned int j=0; j<this->m_MaxTractLength; j++){
416  this->ProbabilisticallyInterpolate( randomgenerator, cindex_curr, index_curr );
417 
418  if(!dwiimagePtr->GetLargestPossibleRegion().IsInside(index_curr)){
419  break;
420  }
421 
422  if( FiberExistenceTest( randomgenerator, wmpimagePtr, index_curr ) ){
423  tract->AddVertex(cindex_curr);
424 
425  this->CalculatePrior( v_prev, this->m_SampleDirections, prior_curr);
426 
427  const ProbabilityDistributionImageType::PixelType&
428  cachelikelihood_curr = this->AccessLikelihoodCache(index_curr);
429 
430  if( cachelikelihood_curr.GetSize() != 0){
431  //use the cached direction
432  this->CalculatePosterior( cachelikelihood_curr, prior_curr, posterior_curr);
433  }
434  else{
435  //do the likelihood calculation and discard
436  //std::cout<<"Cache Miss!\n";
437  ProbabilityDistributionImageType::PixelType
438  likelihood_curr_temp(this->m_SampleDirections->Size());
439 
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);
445  }
446  this->SampleTractOrientation(randomgenerator, posterior_curr,
447  this->m_SampleDirections, v_curr);
448 
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];
455  v_prev=v_curr;
456  }
457  else{
458  //fiber doesn't exist in this voxel
459  //std::cout<<"Stopped Tracking: No Fiber in this Voxel\n";
460  break;
461  }
462  }
463 }
464 
465 template< class TInputDWIImage, class TInputWhiteMatterProbabilityImage, class TOutputConnectivityImage >
466 void
467 StochasticTractographyFilter< TInputDWIImage, TInputWhiteMatterProbabilityImage, TOutputConnectivityImage >
468 ::GenerateTractContainerOutput(){
469  //allocate tractcontainer
470  this->m_OutputTractContainer = TractContainerType::New();
471 
472  this->UpdateGradientDirections();
473  this->UpdateTensorModelFittingMatrices();
474  this->m_TotalDelegatedTracts = 0;
475 
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
485  << std::endl;
486 
487  //setup the multithreader
488  StochasticTractGenerationCallbackStruct data;
489  data.Filter = this;
490  this->GetMultiThreader()->SetSingleMethod( StochasticTractGenerationCallback,
491  &data );
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;
498 }
499 
500 template< class TInputDWIImage, class TInputWhiteMatterProbabilityImage, class TOutputConnectivityImage >
501 void
502 StochasticTractographyFilter< TInputDWIImage, TInputWhiteMatterProbabilityImage, TOutputConnectivityImage >
503 ::GenerateData(){
504  //Generate the tracts
505  this->GenerateTractContainerOutput();
506 
507  //allocate outputs
508  this->AllocateOutputs();
509 
510  //write tracts to output image
511  this->TractContainerToConnectivityMap(this->m_OutputTractContainer);
512 
513 }
514 
515 template< class TInputDWIImage, class TInputWhiteMatterProbabilityImage, class TOutputConnectivityImage >
516 ITK_THREAD_RETURN_TYPE
517 StochasticTractographyFilter< TInputDWIImage, TInputWhiteMatterProbabilityImage, TOutputConnectivityImage >
518 ::StochasticTractGenerationCallback( void *arg )
519 {
520  StochasticTractGenerationCallbackStruct* str=
521  (StochasticTractGenerationCallbackStruct *)
522  (((MultiThreader::ThreadInfoStruct *)(arg))->UserData);
523 
524  typename InputDWIImageType::ConstPointer inputDWIImagePtr = str->Filter->GetInput();
525  typename InputWhiteMatterProbabilityImageType::ConstPointer inputWMPImage =
526  str->Filter->GetWhiteMatterProbabilityImage();
527 
528  unsigned long randomseed=0;
529 
530  while(str->Filter->DelegateTract(randomseed)){
531  //std::cout<<randomseed<<std::endl;
532  //generate the tract
533  TractType::Pointer tract = TractType::New();
534 
535  str->Filter->StochasticTractGeneration( inputDWIImagePtr,
536  inputWMPImage,
537  str->Filter->GetSeedIndex(),
538  randomseed,
539  tract);
540 
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);
545  }
546  else{
547  //std::cout<<"Not Storing Tract\n";
548  }
549  }
550  return ITK_THREAD_RETURN_VALUE;
551 }
552 
553 template< class TInputDWIImage, class TInputWhiteMatterProbabilityImage, class TOutputConnectivityImage >
554 ProbabilityDistributionImageType::PixelType&
555 StochasticTractographyFilter< TInputDWIImage, TInputWhiteMatterProbabilityImage, TOutputConnectivityImage >
556 ::AccessLikelihoodCache( typename InputDWIImageType::IndexType index )
557 {
558  this->m_LikelihoodCacheMutexImagePtr->GetPixel(index).Lock();
559 
560  ProbabilityDistributionImageType::PixelType& likelihood =
561  m_LikelihoodCachePtr->GetPixel( index );
562  typename InputDWIImageType::ConstPointer inputDWIImagePtr = this->GetInput();
563 
564  if( likelihood.GetSize() !=0){
565  //entry found in cache
566  this->m_LikelihoodCacheMutexImagePtr->GetPixel(index).Unlock();
567  return likelihood;
568  }
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());
573 
574  this->CalculateLikelihood(static_cast< DWIVectorImageType::PixelType >(
575  inputDWIImagePtr->GetPixel(index)),
576  this->m_SampleDirections,
577  likelihood);
578  this->m_CurrentLikelihoodCacheElements++;
579 
580  this->m_LikelihoodCacheMutexImagePtr->GetPixel(index).Unlock();
581  return likelihood;
582  }
583  else{
584  //entry not found in cache and no space to store it
585  this->m_LikelihoodCacheMutexImagePtr->GetPixel(index).Unlock();
586  return likelihood;
587  }
588 
589 
590  this->m_LikelihoodCacheMutexImagePtr->GetPixel(index).Unlock();
591 
592  // dummy
593  return likelihood;
594 
595 }
596 
597 template< class TInputDWIImage, class TInputWhiteMatterProbabilityImage, class TOutputConnectivityImage >
598 bool
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++;
606  success = true;
607  //a tract was successfully delegated
608  }
609  else success = false; //all tracts have been delegated
610  this->m_TotalDelegatedTractsMutex.Unlock();
611 
612  return success;
613 }
614 
615 template< class TInputDWIImage, class TInputWhiteMatterProbabilityImage, class TOutputConnectivityImage >
616 void
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);
622 
623  typedef PathIterator< OutputConnectivityImageType, TractType > OutputTractIteratorType;
624 
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,
629  tract );
630 
631  for(outputtractIt.GoToBegin(); !outputtractIt.IsAtEnd(); ++outputtractIt){
632  /* there is an issue using outputtractIt.Value() */
633 // outputtractIt.Set(outputtractIt.Get()+1);
634 
635 
636  }
637  }
638 }
639 
640 template< class TInputDWIImage, class TInputWhiteMatterProbabilityImage, class TOutputConnectivityImage >
641 void
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(),
647  tract);
648  this->m_OutputTractContainerMutex.Unlock();
649 }
650 
651 template< class TInputDWIImage, class TInputWhiteMatterProbabilityImage, class TOutputConnectivityImage >
652 void
653 StochasticTractographyFilter< TInputDWIImage, TInputWhiteMatterProbabilityImage, TOutputConnectivityImage >
654 ::GenerateTensorImageOutput(void){
655  this->UpdateGradientDirections();
656  this->UpdateTensorModelFittingMatrices();
657 
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();
664 
665  //define an iterator for the input and output images
666  typedef itk::ImageRegionConstIterator< InputDWIImageType > DWIImageIteratorType;
667  typedef itk::ImageRegionIterator< OutputTensorImageType > TensorImageIteratorType;
668 
669  DWIImageIteratorType
670  inputDWIit( this->GetInput(), m_OutputTensorImage->GetRequestedRegion() );
671 
672  TensorImageIteratorType outputtensorit
673  ( m_OutputTensorImage, m_OutputTensorImage->GetRequestedRegion() );
674 
675  unsigned int N = this->m_TransformedGradients->Size();
676  TensorModelParamType tensormodelparams( 0.0 );
677  vnl_diag_matrix< double > W(N,0);
678 
679  for(inputDWIit.GoToBegin(), outputtensorit.GoToBegin();
680  !outputtensorit.IsAtEnd(); ++inputDWIit, ++outputtensorit){
681  CalculateTensorModelParameters( inputDWIit.Get(),
682  W, tensormodelparams);
683 
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];
695 
696  //std::cout<<D;
697  }
698 }
699 
700 }