DMRITool  v0.1.1-139-g860d86b4
Diffusion MRI Tool
itkDWIGenerator.hxx
Go to the documentation of this file.
1 /*=========================================================================
2 
3  Program: DWI Generator
4 
5  Copyright (c) Pew-Thian Yap, Jian Cheng. All rights reserved.
6  See http://www.unc.edu/~ptyap/ for details.
7 
8  This software is distributed WITHOUT ANY WARRANTY; without even
9  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
10  PURPOSE. See the above copyright notices for more information.
11 
12  =========================================================================*/
13 
14 #ifndef __itkDWIGenerator_hxx
15 #define __itkDWIGenerator_hxx
16 
17 #include "itkDWIGenerator.h"
18 #include "itkImageRegionIteratorWithIndex.h"
19 #include "itkImage.h"
20 #include "itkDiffusionTensor.h"
22 #include "itksys/SystemTools.hxx"
23 
24 #include "utl.h"
25 
26 namespace itk
27 {
28 
29 template <class TOutputImage, class TScalarImage>
32 {
33 }
34 
35 template <class TOutputImage, class TScalarImage>
38 {
39 }
40 
41 template <class TOutputImage, class TScalarImage>
43 ::PrintSelf(std::ostream& os, Indent indent) const
44 {
45  Superclass::PrintSelf(os, indent);
46  os << indent << "FileName: " << m_FileName << std::endl;
47 }
48 
49 template <class TOutputImage, class TScalarImage>
50 typename LightObject::Pointer
53 {
54  typename LightObject::Pointer loPtr = Superclass::InternalClone();
55 
56  typename Self::Pointer rval = dynamic_cast<Self *>(loPtr.GetPointer());
57  if(rval.IsNull())
58  {
59  itkExceptionMacro(<< "downcast to type " << this->GetNameOfClass()<< " failed.");
60  }
61 
62  rval->m_BackgroundDiffusionParameterValues = m_BackgroundDiffusionParameterValues;
63  rval->m_FileName = m_FileName;
64  return loPtr;
65 }
66 
67 template <class TOutputImage, class TScalarImage>
70 {
71  utlShowPosition(this->GetDebug());
72 
73  this->Initialization();
74 
75  bool isOutputPeak = this->m_MaxNumberOfPeaks>0;
76 
77  MatrixPointer qSpaceOrientationMatrix = this->m_SamplingSchemeQSpace->GetOrientationsCartesian();
78  STDVectorPointer bVector = this->m_SamplingSchemeQSpace->GetBVector();
79  MatrixPointer rSpaceOrientationMatrix = this->m_SamplingSchemeRSpace->GetOrientationsCartesian();
80  STDVectorPointer rVector = this->m_SamplingSchemeRSpace->GetRadiusVector();
81  double tau = this->m_SamplingSchemeQSpace->GetTau();
82  if (this->GetDebug())
83  std::cout << "tau = " << tau << std::endl << std::flush;
84 
85  int numberOfComponentsPerPixel4Peak = isOutputPeak?PeakContainerHelper::GetDimension(this->m_PeakType,this->m_MaxNumberOfPeaks):-1;
86  // utlSAGlobalException(isOutputPeak && numberOfParameterSets>PeakContainerHelper::GetNumberOfPeaks(this->m_PeakType,numberOfComponentsPerPixel4Peak))(numberOfParameterSets)(numberOfComponentsPerPixel4Peak).msg("wrong peak size");
87  OutputImagePixelType pixelPeak;
88  if (isOutputPeak)
89  {
90  pixelPeak.SetSize(numberOfComponentsPerPixel4Peak);
91  }
92 
93  // Parse parameter file
94  // This section should be rewritten to support generic image dimension
95 // char * pch;
96 // unsigned int idx;
97  unsigned int ndims = 0;
98 // char tempLine[256];
99  std::string line, extractedLine;
100  double scale = 1;
101  double sigma = -1;
102  double snr = -1;
103  double backgroundScale = -1;
104 
105  DiffusionParameterValuesType diffusionParameterValues;
106  DiffusionParameterContainerType diffusionParameterContainer;
107  diffusionParameterContainer.clear();
108  DiffusionParameterValuesType backgroundDiffusionParameterValues;
109  backgroundDiffusionParameterValues.clear();
110 
111  // Read Header
112  utlGlobalException(!utl::IsFileExist(m_FileName), "The file does not exist! m_FileName="<<m_FileName);
113 
114  std::vector<std::vector<std::string> > stringMatrix;
115  utl::ReadLines(m_FileName, stringMatrix, " ");
116 
117  for ( int i = 0; i < stringMatrix.size(); i += 1 )
118  {
119  if (stringMatrix[i][0]=="NDims" && stringMatrix[i][1]=="=")
120  std::istringstream ( stringMatrix[i][2] ) >> ndims;
121 
122  if (stringMatrix[i][0]=="DimSize" && stringMatrix[i][1]=="=")
123  {
124  for ( int kk = 0; kk < stringMatrix[i].size()-2; kk += 1 )
125  std::istringstream ( stringMatrix[i][kk+2] ) >> this->m_OutputSize[kk];
126  }
127 
128  if (stringMatrix[i][0]=="ElementSpacing" && stringMatrix[i][1]=="=")
129  {
130  for ( int kk = 0; kk < stringMatrix[i].size()-2; kk += 1 )
131  std::istringstream ( stringMatrix[i][kk+2] ) >> this->m_OutputSpacing[kk];
132  }
133 
134  // if (stringMatrix[i][0]=="BValues" && stringMatrix[i][1]=="=")
135  // {
136  // if (utl::IsNumber(stringMatrix[i][2])) // float numbers
137  // {
138  // bVector->SetSize(stringMatrix[i].size()-2);
139  // for ( int kk = 0; kk < stringMatrix[i].size()-2; kk += 1 )
140  // std::istringstream ( stringMatrix[i][kk+2] ) >> (*bVector)[kk];
141  // }
142  // else // bVector in txt file
143  // {
144  // utlException(stringMatrix[i].size()!=3, "need to set only one txt file");
145  // std::string ext, file;
146  // utl::getFileExtension(stringMatrix[i][2], ext, file);
147  // utlException(ext!="txt", "wrong bVector file");
148  // }
149  // }
150 
151  if (stringMatrix[i][0]=="Scale" && stringMatrix[i][1]=="=")
152  std::istringstream ( stringMatrix[i][2] ) >> this->m_B0Scale;
153 
154  if (stringMatrix[i][0]=="ModelType" && stringMatrix[i][1]=="=")
155  {
156  if (stringMatrix[i][2]=="SYMMETRICAL_TENSOR_IN_CARTESIAN_COORDS")
157  this->m_ModelType = Superclass::SYMMETRICAL_TENSOR_IN_CARTESIAN_COORDS;
158  else if (stringMatrix[i][2]=="SYMMETRICAL_TENSOR_IN_SPHERICAL_COORDS")
159  this->m_ModelType = Superclass::SYMMETRICAL_TENSOR_IN_SPHERICAL_COORDS;
160  else if (stringMatrix[i][2]=="TENSOR_IN_EULER_ANGLES")
161  {
162  this->m_ModelType = Superclass::TENSOR_IN_EULER_ANGLES;
163  utlGlobalException(true, "TODO");
164  }
165  else if (stringMatrix[i][2]=="CYLINDER_SPHERICAL_MODEL")
166  this->m_ModelType = Superclass::CYLINDER_SPHERICAL_MODEL;
167  else
168  utlGlobalException(true, "wrong model type");
169  }
170 
171  if (stringMatrix[i][0]=="BackgroundScale" && stringMatrix[i][1]=="=")
172  std::istringstream ( stringMatrix[i][2] ) >> backgroundScale;
173 
174  if (stringMatrix[i][0]=="DiffusionParameters" && stringMatrix[i][1]=="=")
175  {
176  diffusionParameterValues.clear();
177  for ( int kk = 0; kk < stringMatrix[i].size()-2; kk += 1 )
178  diffusionParameterValues.push_back( atof(stringMatrix[i][kk+2].c_str()) );
179 
180  diffusionParameterContainer.push_back( diffusionParameterValues );
181  }
182 
183  if (stringMatrix[i][0]=="BackgroundDiffusionParameters" && stringMatrix[i][1]=="=")
184  {
185  backgroundDiffusionParameterValues.clear();
186  for ( int kk = 0; kk < stringMatrix[i].size()-2; kk += 1 )
187  backgroundDiffusionParameterValues.push_back( atof(stringMatrix[i][kk+2].c_str()) );
188  }
189 
190  if (stringMatrix[i][0]=="RicianNoiseSigma" && stringMatrix[i][1]=="=")
191  std::istringstream ( stringMatrix[i][2] ) >> sigma;
192 
193  if (stringMatrix[i][0]=="SNR" && stringMatrix[i][1]=="=")
194  std::istringstream ( stringMatrix[i][2] ) >> snr;
195  }
196 
197  if ( this->m_NoiseSigma >= 0 )
198  {
199  sigma = this->m_NoiseSigma;
200  }
201 
202  if ( this->m_SNR >= 0 )
203  {
204  snr = this->m_SNR;
205  }
206 
207  utlException(snr>0 && sigma>0, "snr and sigma can not be set simultaneously. snr="<<snr<<", sigma="<<sigma);
208 
209  scale = this->m_B0Scale;
210  if (backgroundScale<=0)
211  backgroundScale = this->m_B0Scale;
212 
213  // NOTE: for VectorImage, the spacing is determined by the size
214  for ( int i = this->m_OutputSize.GetSizeDimension()-1; i >= 0; i -= 1 )
215  {
216  if (this->m_OutputSize[i]==1)
217  this->m_OutputSpacing[i] = 1;
218  else
219  break;
220  }
221 
222  // allocate all outputs
223  this->AllocateOutputs();
224 
225  // typename OutputImageType::Pointer outputPtr = this->GetOutput();
226  OutputImagePointer outputDWI=this->GetDWIImage(), outputODF=this->GetODFImage(), outputEAP=this->GetEAPImage(), outputPeak=this->GetPeakImage();
227  ScalarImagePointer b0Image = this->GetB0Image(), outputRTO = this->GetRTOImage(), outputMSD = this->GetMSDImage();
228  OutputImagePixelType pixelDWI, tempPixelDWI, pixelEAP, tempPixelEAP, pixelODF, tempPixelODF;
229  double pixelRTO=0, tempPixelRTO=0, pixelMSD=0, tempPixelMSD=0;
230  unsigned int numberOfComponentsPerPixel4DWI = this->GetNumberOfQSpaceSamples();
231  unsigned int numberOfComponentsPerPixel4EAP = this->GetNumberOfRSpaceSamples();
232  unsigned int numberOfComponentsPerPixel4ODF = this->GetNumberOfRSpaceSamples();
233  OutputImagePixelType zeroPixel;
234  if (this->m_IsOutputDWI)
235  {
236  pixelDWI.SetSize( numberOfComponentsPerPixel4DWI );
237  pixelDWI.Fill( 0 );
238  tempPixelDWI.SetSize( numberOfComponentsPerPixel4DWI );
239  tempPixelDWI.Fill( 0 );
240  }
241  if (this->m_IsOutputEAP)
242  {
243  pixelEAP.SetSize( numberOfComponentsPerPixel4EAP );
244  pixelEAP.Fill( 0 );
245  tempPixelEAP.SetSize( numberOfComponentsPerPixel4EAP );
246  tempPixelEAP.Fill( 0 );
247  }
248  if (this->m_IsOutputODF)
249  {
250  pixelODF.SetSize( numberOfComponentsPerPixel4ODF );
251  pixelODF.Fill( 0 );
252  tempPixelODF.SetSize( numberOfComponentsPerPixel4ODF );
253  tempPixelODF.Fill( 0 );
254  }
255 
256  OutputImageIndexType index;
257 
258 
259  // Convenient variables
260  int numberOfParameterSets = -1;
261  unsigned int numberOfOrientationComponentsPerParameterSet = 3;
262  int numberOfComponentsPerParameterSet = -1; // Added partial volume fractions
263  int numberOfDiffusivityComponentsPerParameterSet = -1;
264 
265  bool isTensorModel=false;
267 
268  bool isCylinderModel=false;
269  typename CylinderModelType::PointType cylinderAxis;
270  typename CylinderModelType::VectorPointer tempPixelDWIVec, tempPixelODFVec, tempPixelEAPVec;
271 
272  if (this->m_ModelType==Superclass::SYMMETRICAL_TENSOR_IN_CARTESIAN_COORDS)
273  {
274  numberOfDiffusivityComponentsPerParameterSet=2;
275  numberOfOrientationComponentsPerParameterSet = 3;
276  numberOfComponentsPerParameterSet = 6;
277  isTensorModel = true;
278  }
279  else if (this->m_ModelType==Superclass::SYMMETRICAL_TENSOR_IN_SPHERICAL_COORDS)
280  {
281  numberOfDiffusivityComponentsPerParameterSet=2;
282  numberOfOrientationComponentsPerParameterSet = 2;
283  numberOfComponentsPerParameterSet = 5;
284  isTensorModel = true;
285  }
286  else if (this->m_ModelType==Superclass::TENSOR_IN_EULER_ANGLES)
287  {
288  numberOfDiffusivityComponentsPerParameterSet=3;
289  numberOfOrientationComponentsPerParameterSet = 3;
290  numberOfComponentsPerParameterSet = 7;
291  isTensorModel = true;
292  utlGlobalException(true, "TODO");
293  }
294  else if (this->m_ModelType==Superclass::CYLINDER_SPHERICAL_MODEL)
295  {
296  numberOfDiffusivityComponentsPerParameterSet=0;
297  numberOfOrientationComponentsPerParameterSet = 2;
298  numberOfComponentsPerParameterSet = 3;
299  isCylinderModel = true;
300  }
301 
302  utlGlobalException(!isTensorModel && !isCylinderModel, "wrong model type");
303 
304  if (isCylinderModel)
305  {
306  this->m_CylinderModel->SetDebug(this->GetDebug());
307  this->m_CylinderModel->BuildTable();
308  this->m_CylinderModel->SetSamplingSchemeQSpace(this->m_SamplingSchemeQSpace);
309  this->m_CylinderModel->SetSamplingSchemeRSpace(this->m_SamplingSchemeRSpace);
310  }
311 
312  Vector<PrecisionType, 3> principalDirection(0.0);
313  STDVectorType diffusivities;
314  if (numberOfDiffusivityComponentsPerParameterSet>0)
315  diffusivities.resize(numberOfDiffusivityComponentsPerParameterSet);
316  Vector<PrecisionType, 3> e1;
317  e1[0]=1.0, e1[1]=0.0, e1[2]=0.0;
318  double partialVolumeWeight = 0;
319 
320  double norm;
321 // double temp;
322 
323  // std::cout << "scale = " << scale << std::endl;
324  // std::cout << "backgroundScale = " << backgroundScale << std::endl;
325  // utlPrintVar1(true, numberOfDiffusivityComponentsPerParameterSet);
326 
327  // Simulated DWI data
328  for ( unsigned int k=0; k<diffusionParameterContainer.size(); k++ )
329  {
330  diffusionParameterValues = diffusionParameterContainer[k];
331 
332  for ( unsigned int dim=0; dim<ndims; dim++ )
333  {
334  index[dim] = diffusionParameterValues[dim];
335  }
336  if (this->GetDebug())
337  {
338  std::cout << "index = " << index << std::endl << std::flush;
339  utl::PrintVector(diffusionParameterValues, "diffusionParameterValues");
340  }
341 
342  // 3 elements for directions, 2 elements for diffusivities
343  utlSAGlobalException(!utl::IsInt( (1.0*diffusionParameterValues.size()-ndims)/numberOfComponentsPerParameterSet) )(diffusionParameterValues.size())(ndims)(numberOfComponentsPerParameterSet).msg("wrong diffusion parameter size");
344  numberOfParameterSets = ( diffusionParameterValues.size() - ndims ) /numberOfComponentsPerParameterSet;
345  if (this->m_IsOutputDWI)
346  pixelDWI.Fill( 0 );
347  if (this->m_IsOutputEAP)
348  pixelEAP.Fill( 0 );
349  if (this->m_IsOutputODF)
350  pixelODF.Fill( 0 );
351  if (isOutputPeak)
352  pixelPeak.Fill( 0 );
353  pixelMSD=0, pixelRTO=0;
354 
355  // std::cout << "k = "<< k << ", pixel = " << pixel << std::endl << std::flush;
356  // utlPrintVar4(true, numberOfParameterSets, diffusionParameterValues.size(), ndims, numberOfComponentsPerParameterSet);
357  PrecisionType sumPartialVolumeWeight=0.0;
358  int numberOfPeaks=0;
359  for ( unsigned int s=0; s<numberOfParameterSets; s++ )
360  {
361  for ( unsigned int d=0; d<numberOfOrientationComponentsPerParameterSet; d++ )
362  principalDirection[d] = diffusionParameterValues[ndims + s* numberOfComponentsPerParameterSet + d];
363 
364  norm = principalDirection.GetNorm();
365  if (this->m_ModelType==Superclass::SYMMETRICAL_TENSOR_IN_SPHERICAL_COORDS || this->m_ModelType==Superclass::CYLINDER_SPHERICAL_MODEL)
366  {
367  principalDirection *= vnl_math::pi/180.0;
368  double theta=principalDirection[0], phi=principalDirection[1];
369  utl::spherical2Cartesian(1.0,theta,phi,principalDirection[0],principalDirection[1],principalDirection[2]);
370  }
371  else
372  {
373  for ( unsigned int d=0; d<3; d++ )
374  principalDirection[d] /= norm;
375  }
376 
377  for ( unsigned int d=0; d<numberOfDiffusivityComponentsPerParameterSet; d++ )
378  {
379  diffusivities[d] = diffusionParameterValues[ndims + s
380  * numberOfComponentsPerParameterSet
381  + numberOfOrientationComponentsPerParameterSet + d ];
382  }
383  if (diffusivities.size()>0)
384  std::sort(diffusivities.begin(), diffusivities.end(), std::greater<double>());
385 
386  partialVolumeWeight = diffusionParameterValues[ndims + s
387  * numberOfComponentsPerParameterSet
388  + numberOfOrientationComponentsPerParameterSet
389  + numberOfDiffusivityComponentsPerParameterSet];
390 
391  sumPartialVolumeWeight += partialVolumeWeight;
392 
393  // utlPrintVar4(true, k, index[0], index[1], index[2]);
394  // utl::PrintVector(diffusivities, "diffusivities");
395 
396  if (isOutputPeak)
397  {
398  PeakContainerHelper::SetPeak<Vector<double,3>, OutputImagePixelType >(principalDirection, pixelPeak, s, this->m_PeakType);
399  if (this->m_PeakType==NXYZV || this->m_PeakType==XYZV)
400  PeakContainerHelper::SetPeakValue<OutputImagePixelType >(partialVolumeWeight, pixelPeak, s, this->m_PeakType);
401  numberOfPeaks++;
402  }
403 
404  if (isTensorModel)
405  {
406  std::vector<double> vec(3,1);
407  if (this->m_ModelType==Superclass::SYMMETRICAL_TENSOR_IN_CARTESIAN_COORDS || this->m_ModelType==Superclass::SYMMETRICAL_TENSOR_IN_SPHERICAL_COORDS)
408  vec[0]=diffusivities[0], vec[1]=diffusivities[1], vec[2]=diffusivities[1];
409  else if (this->m_ModelType==Superclass::TENSOR_IN_EULER_ANGLES)
410  vec[0]=diffusivities[0], vec[1]=diffusivities[1], vec[2]=diffusivities[2];
411  tensor.Fill(0.0);
412  tensor.SetEigenValues(vec);
413  // utlPrintVar2(true, "before", tensor);
414  Matrix<double,3,3> rotation;
415  utl::RotationMatrixFromVectors<Vector<double,3>, Matrix<double,3> >(e1, principalDirection, rotation);
416  tensor.Rotate(rotation);
417  // utlPrintVar2(true, "after", tensor);
418 
419  if ( this->m_IsOutputDWI )
420  {
421  tensor.GetDWISamples(tempPixelDWI, *qSpaceOrientationMatrix, *bVector);
422  // utlPrintVar1(true, tempPixel);
423  // utlPrintVar1(true, tempPixel*scale*partialVolumeWeight);
424  pixelDWI += tempPixelDWI*scale*partialVolumeWeight;
425  }
426  if ( this->m_IsOutputODF )
427  {
428  tensor.GetODFSamples(tempPixelODF, *rSpaceOrientationMatrix,this->m_ODFOrder, false);
429  pixelODF += tempPixelODF*partialVolumeWeight;
430  }
431  if ( this->m_IsOutputEAP)
432  {
433  tensor.GetEAPSamples(tempPixelEAP, *rSpaceOrientationMatrix, *rVector,tau);
434  pixelEAP += tempPixelEAP*partialVolumeWeight;
435  // std::cout << "\n\nindex = " << index << std::endl << std::flush;
436  }
437  if ( this->m_IsOutputRTO)
438  {
439  tempPixelRTO = tensor.GetReturnToOrigin(tau);
440  pixelRTO += tempPixelRTO*partialVolumeWeight;
441  }
442  if ( this->m_IsOutputMSD)
443  {
444  tempPixelMSD = tensor.GetMeanSquaredDisplacement(tau);
445  pixelMSD += tempPixelMSD*partialVolumeWeight;
446  }
447  }
448  else if (isCylinderModel)
449  {
450  cylinderAxis[0] = principalDirection[0];
451  cylinderAxis[1] = principalDirection[1];
452  cylinderAxis[2] = principalDirection[2];
453  this->m_CylinderModel->SetCylinderAxis(cylinderAxis);
454  if (this->m_IsOutputDWI)
455  {
456  if (this->GetDebug())
457  this->m_CylinderModel->Print(std::cout<<"this->m_CylinderModel=");
458  this->m_CylinderModel->ComputeDWISamples();
459  tempPixelDWIVec = this->m_CylinderModel->GetDWISamples();
460  for ( int i = 0; i < tempPixelDWIVec->size(); i += 1 )
461  pixelDWI[i] += (*tempPixelDWIVec)[i]*scale*partialVolumeWeight;
462  }
463  if ( this->m_IsOutputODF )
464  {
465  utlGlobalException(true, "TODO: to be implemented");
466  }
467  if ( this->m_IsOutputEAP )
468  {
469  utlGlobalException(true, "TODO: to be implemented");
470  }
471  if ( this->m_IsOutputRTO )
472  {
473  utlGlobalException(true, "TODO: to be implemented");
474  }
475  if ( this->m_IsOutputMSD )
476  {
477  utlGlobalException(true, "TODO: to be implemented");
478  }
479  }
480 
481  }
482  if (isOutputPeak)
483  {
484  if (this->m_PeakType==NXYZ || this->m_PeakType==NXYZV)
485  pixelPeak[0] = numberOfPeaks;
486  }
487 
488  utlException(sumPartialVolumeWeight<1e-9, "wrong partialVolumeWeight");
489  if ( this->m_IsOutputDWI )
490  pixelDWI /= sumPartialVolumeWeight;
491  if ( this->m_IsOutputODF )
492  pixelODF /= sumPartialVolumeWeight;
493  if ( this->m_IsOutputEAP )
494  pixelEAP /= sumPartialVolumeWeight;
495  if ( this->m_IsOutputRTO )
496  pixelRTO /= sumPartialVolumeWeight;
497  if ( this->m_IsOutputMSD )
498  pixelMSD /= sumPartialVolumeWeight;
499  if (this->m_IsOutputODF && this->m_ODFOrder!=2)
500  {
501  PrecisionType normFactor = 4*vnl_math::pi*utl::GetSumOfVector<OutputImagePixelType>(pixelODF,pixelODF.Size()) / pixelODF.Size();
502  if (normFactor!=0)
503  pixelODF /= normFactor;
504  }
505 
506  // Add Rician noise
507  if ( this->m_IsOutputDWI && (sigma>0 || snr>0) && pixelDWI.GetSquaredNorm()>0 )
508  {
509  PrecisionType sigma_real = -1;
510  if ( snr > 0 && sigma <= 0)
511  {
512  double mean = 0;
513  for ( unsigned int k=0; k<numberOfComponentsPerPixel4DWI; k++ )
514  mean += pixelDWI[k];
515  mean /= numberOfComponentsPerPixel4DWI;
516  sigma_real = mean/snr;
517  }
518  else if ( snr <= 0 && sigma > 0)
519 // sigma_real = scale * sigma;
520  sigma_real = sigma;
521  pixelDWI = utl::AddNoise<OutputImagePixelType>(pixelDWI, pixelDWI.GetSize(), sigma_real, true);
522  }
523 
524  if ( isOutputPeak )
525  {
526  outputPeak->SetPixel(index, pixelPeak);
527  }
528  if ( this->m_IsOutputDWI )
529  {
530  outputDWI->SetPixel(index, pixelDWI);
531  b0Image->SetPixel(index, scale);
532  }
533  if ( this->m_IsOutputODF )
534  outputODF->SetPixel(index, pixelODF);
535  if ( this->m_IsOutputEAP )
536  outputEAP->SetPixel(index, pixelEAP);
537  if ( this->m_IsOutputRTO )
538  outputRTO->SetPixel(index, pixelRTO);
539  if ( this->m_IsOutputMSD )
540  outputMSD->SetPixel(index, pixelMSD);
541  }
542 
543  ImageRegionIteratorWithIndex<OutputImageType> outputDWIIt, outputODFIt, outputEAPIt, outputPeakIt;
544  ImageRegionIteratorWithIndex<ScalarImageType> b0It, outputRTOIt, outputMSDIt;
545  if ( this->m_IsOutputDWI )
546  {
547  outputDWIIt = ImageRegionIteratorWithIndex<OutputImageType>( outputDWI, outputDWI->GetRequestedRegion() );
548  b0It = ImageRegionIteratorWithIndex<ScalarImageType>( b0Image, outputDWI->GetRequestedRegion() );
549  }
550  if ( this->m_IsOutputODF )
551  outputODFIt = ImageRegionIteratorWithIndex<OutputImageType>( outputODF, outputODF->GetRequestedRegion() );
552  if ( this->m_IsOutputEAP )
553  outputEAPIt = ImageRegionIteratorWithIndex<OutputImageType>( outputEAP, outputEAP->GetRequestedRegion() );
554  if ( isOutputPeak )
555  outputPeakIt = ImageRegionIteratorWithIndex<OutputImageType>( outputPeak, outputPeak->GetRequestedRegion() );
556  if ( this->m_IsOutputRTO )
557  outputRTOIt = ImageRegionIteratorWithIndex<ScalarImageType>( outputRTO, outputRTO->GetRequestedRegion() );
558  if ( this->m_IsOutputMSD )
559  outputMSDIt = ImageRegionIteratorWithIndex<ScalarImageType>( outputMSD, outputMSD->GetRequestedRegion() );
560 
561 
562  if (m_BackgroundDiffusionParameterValues.size()>0)
563  backgroundDiffusionParameterValues = m_BackgroundDiffusionParameterValues;
564 
565  // Background
566  if ( backgroundDiffusionParameterValues.size() > 0 )
567  {
568  numberOfParameterSets = backgroundDiffusionParameterValues.size() /numberOfComponentsPerParameterSet;
569  // utl::PrintVector(backgroundDiffusionParameterValues, "backgroundDiffusionParameterValues");
570  // utlPrintVar3(true, numberOfComponentsPerParameterSet, numberOfParameterSets, numberOfOrientationComponentsPerParameterSet);
571  if (this->m_IsOutputDWI)
572  pixelDWI.Fill( 0 );
573  if (this->m_IsOutputEAP)
574  pixelEAP.Fill( 0 );
575  if (this->m_IsOutputODF)
576  pixelODF.Fill( 0 );
577  pixelRTO=0, pixelMSD=0;
578  PrecisionType sumPartialVolumeWeight=0.0;
579 
580  for ( unsigned int s=0; s<numberOfParameterSets; s++ )
581  {
582  for ( unsigned int d=0; d<numberOfOrientationComponentsPerParameterSet; d++ )
583  principalDirection[d] = backgroundDiffusionParameterValues[s* numberOfComponentsPerParameterSet + d];
584 
585  // std::cout << "principalDirection = " << principalDirection << std::endl << std::flush;
586  // utlException (std::abs(norm)<1e-8, "directions have zero norm");
587  norm = principalDirection.GetNorm();
588  if (this->m_ModelType==Superclass::SYMMETRICAL_TENSOR_IN_SPHERICAL_COORDS || this->m_ModelType==Superclass::CYLINDER_SPHERICAL_MODEL)
589  {
590  principalDirection *= vnl_math::pi/180.0;
591  double theta=principalDirection[0], phi=principalDirection[1];
592  utl::spherical2Cartesian(1.0,theta,phi,principalDirection[0],principalDirection[1],principalDirection[2]);
593  }
594  else
595  {
596  for ( unsigned int d=0; d<3; d++ )
597  principalDirection[d] /= norm;
598  }
599 
600  for ( unsigned int d=0; d<numberOfDiffusivityComponentsPerParameterSet; d++ )
601  {
602  diffusivities[d] = backgroundDiffusionParameterValues[s* numberOfComponentsPerParameterSet
603  + numberOfOrientationComponentsPerParameterSet + d ];
604  }
605  if (diffusivities.size()>0)
606  std::sort(diffusivities.begin(), diffusivities.end(), std::greater<double>());
607 
608  partialVolumeWeight = backgroundDiffusionParameterValues[s* numberOfComponentsPerParameterSet
609  + numberOfOrientationComponentsPerParameterSet
610  + numberOfDiffusivityComponentsPerParameterSet];
611 
612  sumPartialVolumeWeight += partialVolumeWeight;
613 
614  if (isTensorModel)
615  {
616  std::vector<double> vec(3,1);
617  if (this->m_ModelType==Superclass::SYMMETRICAL_TENSOR_IN_CARTESIAN_COORDS || this->m_ModelType==Superclass::SYMMETRICAL_TENSOR_IN_SPHERICAL_COORDS)
618  vec[0]=diffusivities[0], vec[1]=diffusivities[1], vec[2]=diffusivities[1];
619  else if (this->m_ModelType==Superclass::TENSOR_IN_EULER_ANGLES)
620  vec[0]=diffusivities[0], vec[1]=diffusivities[1], vec[2]=diffusivities[2];
621  tensor.Fill(0.0);
622  tensor.SetEigenValues(vec);
623  // tensor.Print(std::cout<<"tensor = ");
624  Matrix<double,3,3> rotation;
625  utl::RotationMatrixFromVectors<Vector<double,3>, Matrix<double,3> >(e1, principalDirection, rotation);
626  tensor.Rotate(rotation);
627  // tensor.Print(std::cout<<"tensor = ");
628 
629  if ( this->m_IsOutputDWI )
630  {
631  tensor.GetDWISamples(tempPixelDWI, *qSpaceOrientationMatrix, *bVector);
632  pixelDWI += tempPixelDWI*backgroundScale*partialVolumeWeight;
633  }
634  if ( this->m_IsOutputODF )
635  {
636  tensor.GetODFSamples(tempPixelODF, *rSpaceOrientationMatrix,this->m_ODFOrder, false);
637  pixelODF += tempPixelODF*partialVolumeWeight;
638  }
639  if ( this->m_IsOutputEAP)
640  {
641  tensor.GetEAPSamples(tempPixelEAP, *rSpaceOrientationMatrix, *rVector,tau);
642  pixelEAP += tempPixelEAP*partialVolumeWeight;
643  }
644  if ( this->m_IsOutputRTO)
645  {
646  tempPixelRTO = tensor.GetReturnToOrigin(tau);
647  pixelRTO += tempPixelRTO*partialVolumeWeight;
648  }
649  if ( this->m_IsOutputMSD)
650  {
651  tempPixelMSD = tensor.GetMeanSquaredDisplacement(tau);
652  pixelMSD += tempPixelMSD*partialVolumeWeight;
653  }
654  }
655  else if (isCylinderModel)
656  {
657  cylinderAxis[0] = principalDirection[0];
658  cylinderAxis[1] = principalDirection[1];
659  cylinderAxis[2] = principalDirection[2];
660  this->m_CylinderModel->SetCylinderAxis(cylinderAxis);
661  if (this->m_IsOutputDWI)
662  {
663  // cyl->Print(std::cout<<"cy1=");
664  this->m_CylinderModel->ComputeDWISamples();
665  tempPixelDWIVec = this->m_CylinderModel->GetDWISamples();
666  for ( int i = 0; i < tempPixelDWIVec->size(); i += 1 )
667  pixelDWI[i] += (*tempPixelDWIVec)[i]*backgroundScale*partialVolumeWeight;
668  }
669 
670  }
671 
672  }
673  utlException(sumPartialVolumeWeight<1e-9, "wrong partialVolumeWeight");
674  if ( this->m_IsOutputDWI )
675  pixelDWI /= sumPartialVolumeWeight;
676  if ( this->m_IsOutputODF )
677  pixelODF /= sumPartialVolumeWeight;
678  if ( this->m_IsOutputEAP )
679  pixelEAP /= sumPartialVolumeWeight;
680  if ( this->m_IsOutputRTO )
681  pixelRTO /= sumPartialVolumeWeight;
682  if ( this->m_IsOutputMSD )
683  pixelMSD /= sumPartialVolumeWeight;
684  if (this->m_IsOutputODF && this->m_ODFOrder!=2)
685  {
686  PrecisionType normFactor = 4*vnl_math::pi*utl::GetSumOfVector<OutputImagePixelType>(pixelODF,pixelODF.Size()) / pixelODF.Size();
687  if (normFactor!=0)
688  pixelODF /= normFactor;
689  }
690 
691  // Add Rician noise
692  PrecisionType sigma_real = -1;
693  if ( this->m_IsOutputDWI && (sigma>0 || snr>0) && pixelDWI.GetSquaredNorm()>0 )
694  {
695  if ( snr > 0 && sigma <= 0)
696  {
697  double mean = 0;
698  for ( unsigned int k=0; k<numberOfComponentsPerPixel4DWI; k++ )
699  mean += pixelDWI[k];
700  mean /= numberOfComponentsPerPixel4DWI;
701  sigma_real = mean/snr;
702  }
703  else if ( snr <= 0 && sigma > 0)
704 // sigma_real = backgroundScale * sigma;
705  sigma_real = sigma;
706  }
707 
708  // Set pixels
709  if (this->m_IsOutputDWI)
710  {
711  b0It.GoToBegin();
712  outputDWIIt.GoToBegin();
713  OutputImagePixelType pixel_noise=pixelDWI;
714  while ( !b0It.IsAtEnd() )
715  {
716  if ( outputDWIIt.Get().GetNorm() == 0 )
717  {
718  if (sigma_real>0)
719  pixel_noise = utl::AddNoise<OutputImagePixelType>(pixelDWI, pixelDWI.GetSize(), sigma_real, true);
720  outputDWIIt.Set( pixel_noise );
721  b0It.Set( backgroundScale );
722  }
723  ++outputDWIIt;
724  ++b0It;
725  }
726  }
727 
728  if (this->m_IsOutputODF)
729  {
730  outputODFIt.GoToBegin();
731  while ( !outputODFIt.IsAtEnd() )
732  {
733  if ( outputODFIt.Get().GetNorm() == 0 )
734  outputODFIt.Set( pixelODF );
735  ++outputODFIt;
736  }
737  }
738 
739  if (this->m_IsOutputEAP)
740  {
741  outputEAPIt.GoToBegin();
742  while ( !outputEAPIt.IsAtEnd() )
743  {
744  if ( outputEAPIt.Get().GetNorm() == 0 )
745  outputEAPIt.Set( pixelEAP );
746  ++outputEAPIt;
747  }
748  }
749 
750  if (this->m_IsOutputRTO)
751  {
752  outputRTOIt.GoToBegin();
753  while ( !outputRTOIt.IsAtEnd() )
754  {
755  if ( outputRTOIt.Get() == 0 )
756  outputRTOIt.Set( pixelRTO );
757  ++outputRTOIt;
758  }
759  }
760  if (this->m_IsOutputMSD)
761  {
762  outputMSDIt.GoToBegin();
763  while ( !outputMSDIt.IsAtEnd() )
764  {
765  if ( outputMSDIt.Get() == 0 )
766  outputMSDIt.Set( pixelMSD );
767  ++outputMSDIt;
768  }
769  }
770  }
771 
772 }
773 
774 
775 } //namespace ITK
776 
777 #endif
Superclass::DiffusionParameterValuesType DiffusionParameterValuesType
void GetEAPSamples(TVectorType &vec, const utl::NDArray< TPrecision, 2 > &gradients, const std::vector< TPrecision > &rValues, const TPrecision &tau) const
Superclass::STDVectorPointer STDVectorPointer
Superclass::OutputImagePointer OutputImagePointer
SmartPointer< Self > Pointer
helper functions specifically used in dmritool
void GetDWISamples(TVectorType &vec, const utl::NDArray< TPrecision, 2 > &gradients, const std::vector< TPrecision > &bValues) const
static int GetDimension(const PeakType peakType, const int numberOfPeaks)
#define utlException(cond, expout)
Definition: utlCoreMacro.h:548
bool IsFileExist(const std::string &file)
Definition: utlCore.h:529
void spherical2Cartesian(const T r, const T theta, const T phi, T &x, T &y, T &z)
Definition: utlCore.h:1602
RealValueType GetReturnToOrigin(const TPrecision tau) const
#define utlSAGlobalException(expr)
Definition: utlCoreMacro.h:362
Superclass::OutputImageIndexType OutputImageIndexType
static void sort(int *irOut, T *prOut, int beg, int end)
Definition: misc.h:141
void PrintSelf(std::ostream &os, Indent indent) const ITK_OVERRIDE
Self & Rotate(const typename Self::MatrixType &mat)
Superclass::ScalarImagePointer ScalarImagePointer
Superclass::DiffusionParameterContainerType DiffusionParameterContainerType
#define utlGlobalException(cond, expout)
Definition: utlCoreMacro.h:372
Superclass::VectorPointer VectorPointer
bool IsInt(const std::string &input, const double epss=1e-10)
Definition: utlCore.h:792
Superclass::MatrixPointer MatrixPointer
LightObject::Pointer InternalClone() const ITK_OVERRIDE
void SetEigenValues(const TArrayType &array)
Generate DWI data based on provided parameter file.
Superclass::STDVectorType STDVectorType
void PrintVector(const std::vector< T > &vec, const std::string &str="", const char *separate=" ", std::ostream &os=std::cout, bool showStats=true)
Definition: utlCore.h:1002
Generate DWI data based on provided parameter file.
void ReadLines(const std::string &filename, std::vector< std::vector< std::string > > &strVec, const char *cc=" ")
Definition: utlCore.h:862
Superclass::OutputImagePixelType OutputImagePixelType
#define utlShowPosition(cond)
Definition: utlCoreMacro.h:554
Superclass::PrecisionType PrecisionType
RealValueType GetMeanSquaredDisplacement(const TPrecision tau) const
void GetODFSamples(TVectorType &vec, const utl::NDArray< TPrecision, 2 > &gradients, const int &odfOrder=2, const bool &isNormalize=false) const
tensor with some useful functions