DMRITool  v0.1.1-139-g860d86b4
Diffusion MRI Tool
itkDWISingleVoxelGenerator.hxx
Go to the documentation of this file.
1 
18 #ifndef __itkDWISingleVoxelGenerator_hxx
19 #define __itkDWISingleVoxelGenerator_hxx
20 
22 #include "itkDiffusionTensor.h"
24 
25 namespace itk
26 {
27 
28 template <class TOutputImage, class TScalarImage>
31 : Superclass(),
32  m_StoredOrientationMatrix(new MatrixType())
33 {
34  m_RandomType = FIXED;
35 }
36 
37 template <class TOutputImage, class TScalarImage>
39 ::PrintSelf(std::ostream& os, Indent indent) const
40 {
41  Superclass::PrintSelf(os, indent);
42  PrintVar1(true, m_RandomType, os<<indent);
43  if (m_StoredOrientationMatrix->Rows()>0)
44  utl::PrintUtlMatrix(*m_StoredOrientationMatrix, "m_StoredOrientationMatrix");
45 }
46 
47 template <class TOutputImage, class TScalarImage>
48 typename LightObject::Pointer
51 {
52  typename LightObject::Pointer loPtr = Superclass::InternalClone();
53 
54  typename Self::Pointer rval = dynamic_cast<Self *>(loPtr.GetPointer());
55  if(rval.IsNull())
56  {
57  itkExceptionMacro(<< "downcast to type " << this->GetNameOfClass()<< " failed.");
58  }
59 
60  rval->m_DiffusionParameterValues = m_DiffusionParameterValues;
61  rval->m_RandomType = m_RandomType;
62  rval->m_StoredOrientationMatrix = m_StoredOrientationMatrix;
63  return loPtr;
64 }
65 
66 template <class TOutputImage, class TScalarImage>
69 {
70  Superclass::Initialization();
71  if (m_RandomType==UNIFORM && m_StoredOrientationMatrix->Rows())
72  {
73  for ( int i = 0; i < TOutputImage::ImageDimension; i += 1 )
74  this->m_OutputSize[i] = 1.0;
75  this->m_OutputSize[0] = m_StoredOrientationMatrix->Rows();
76  }
77 }
78 
79 template <class TOutputImage, class TScalarImage>
82 {
83  utlShowPosition(this->GetDebug());
84  this->Initialization();
85 
86  bool isOutputPeak = this->m_MaxNumberOfPeaks>0;
87 
88  MatrixPointer qSpaceOrientationMatrix = this->m_SamplingSchemeQSpace->GetOrientationsCartesian();
89  STDVectorPointer bVector = this->m_SamplingSchemeQSpace->GetBVector();
90  MatrixPointer rSpaceOrientationMatrix = this->m_SamplingSchemeRSpace->GetOrientationsCartesian();
91  STDVectorPointer rVector = this->m_SamplingSchemeRSpace->GetRadiusVector();
92  double tau = this->m_SamplingSchemeQSpace->GetTau();
93 
94  // allocate all outputs
95  this->AllocateOutputs();
96 
97  int numberOfDiffusivityComponentsPerParameterSet = -1;
98  unsigned int numberOfOrientationComponentsPerParameterSet = 3;
99  if (this->m_ModelType==Superclass::SYMMETRICAL_TENSOR_IN_CARTESIAN_COORDS)
100  {
101  numberOfDiffusivityComponentsPerParameterSet=2;
102  numberOfOrientationComponentsPerParameterSet = 3;
103  }
104  else if (this->m_ModelType==Superclass::SYMMETRICAL_TENSOR_IN_SPHERICAL_COORDS)
105  {
106  numberOfDiffusivityComponentsPerParameterSet=2;
107  numberOfOrientationComponentsPerParameterSet = 2;
108  }
109  else if (this->m_ModelType==Superclass::TENSOR_IN_EULER_ANGLES)
110  {
111  numberOfDiffusivityComponentsPerParameterSet=3;
112  numberOfOrientationComponentsPerParameterSet = 3;
113  utlGlobalException(true, "TODO");
114  }
115  else if (this->m_ModelType==Superclass::CYLINDER_SPHERICAL_MODEL)
116  {
117  numberOfOrientationComponentsPerParameterSet = 2;
118  utlGlobalException(true, "TODO");
119  }
120 
121  int numberOfComponentsPerParameterSet = numberOfDiffusivityComponentsPerParameterSet + numberOfOrientationComponentsPerParameterSet+1;
122  // utlSAPrint(m_DiffusionParameterValues.size())(numberOfDiffusivityComponentsPerParameterSet)(numberOfOrientationComponentsPerParameterSet)(numberOfComponentsPerParameterSet);
123  utlSAGlobalException(!utl::IsInt(m_DiffusionParameterValues.size()*1.0 / numberOfComponentsPerParameterSet))(m_DiffusionParameterValues.size())(numberOfComponentsPerParameterSet).msg("wrong diffusion parameter size");
124  int numberOfParameterSets = m_DiffusionParameterValues.size() /numberOfComponentsPerParameterSet;
125  int numberOfComponentsPerPixel4Peak = isOutputPeak?PeakContainerHelper::GetDimension(this->m_PeakType,this->m_MaxNumberOfPeaks):-1;
126  utlSAGlobalException(isOutputPeak && numberOfParameterSets>PeakContainerHelper::GetNumberOfPeaks(this->m_PeakType,numberOfComponentsPerPixel4Peak))(numberOfParameterSets)(numberOfComponentsPerPixel4Peak).msg("wrong peak size");
127 
128  std::vector<DiffusionTensor<double> > tensorVec;
129  std::vector<double> weightVec;
130 
132  Vector<PrecisionType, 3> principalDirection(0.0);
133  double norm, partialVolumeWeight;
134  Vector<PrecisionType, 3> e1;
135  e1[0]=1.0, e1[1]=0.0, e1[2]=0.0;
136  PrecisionType sumPartialVolumeWeight=0.0;
137  STDVectorType diffusivities(numberOfDiffusivityComponentsPerParameterSet);
138  // utl::PrintVector(m_DiffusionParameterValues, "m_DiffusionParameterValues");
139  // utlSAPrint(numberOfParameterSets)(numberOfOrientationComponentsPerParameterSet)(numberOfComponentsPerParameterSet);
140 
141  OutputImagePixelType pixelPeakInput, pixelPeak;
142  int numberOfPeaks=0;
143  if (isOutputPeak)
144  {
145  pixelPeakInput.SetSize(numberOfComponentsPerPixel4Peak);
146  pixelPeakInput.Fill(0.0);
147  if (this->m_PeakType==NXYZ || this->m_PeakType==NXYZV)
148  pixelPeakInput[0] = numberOfParameterSets;
149  pixelPeak.SetSize(numberOfComponentsPerPixel4Peak);
150  }
151  for ( unsigned int s=0; s<numberOfParameterSets; s++ )
152  {
153  for ( unsigned int d=0; d<numberOfOrientationComponentsPerParameterSet; d++ )
154  principalDirection[d] = m_DiffusionParameterValues[s* numberOfComponentsPerParameterSet + d];
155 
156  norm = principalDirection.GetNorm();
157  if (this->m_ModelType==Superclass::SYMMETRICAL_TENSOR_IN_SPHERICAL_COORDS || this->m_ModelType==Superclass::CYLINDER_SPHERICAL_MODEL)
158  {
159  principalDirection *= vnl_math::pi/180.0;
160  double theta=principalDirection[0], phi=principalDirection[1];
161  utl::spherical2Cartesian(1.0,theta,phi,principalDirection[0],principalDirection[1],principalDirection[2]);
162  }
163  else
164  {
165  for ( unsigned int d=0; d<numberOfOrientationComponentsPerParameterSet; d++ )
166  principalDirection[d] /= norm;
167  }
168 
169  for ( unsigned int d=0; d<numberOfDiffusivityComponentsPerParameterSet; d++ )
170  {
171  diffusivities[d] = m_DiffusionParameterValues[s*numberOfComponentsPerParameterSet + numberOfOrientationComponentsPerParameterSet + d ];
172  }
173  std::sort(diffusivities.begin(), diffusivities.end(), std::greater<double>());
174 
175  partialVolumeWeight = m_DiffusionParameterValues[s*numberOfComponentsPerParameterSet + numberOfOrientationComponentsPerParameterSet + numberOfDiffusivityComponentsPerParameterSet];
176  sumPartialVolumeWeight += partialVolumeWeight;
177  weightVec.push_back(partialVolumeWeight);
178 
179  // utlPrintVar2(true, s, PeakContainerHelper::GetString(this->m_PeakType));
180  // std::cout << "principalDirection = " << principalDirection << std::endl << std::flush;
181  // utl::PrintVector(diffusivities, "diffusivities");
182 
183  std::vector<double> vec(3,1);
184  if (this->m_ModelType==Superclass::SYMMETRICAL_TENSOR_IN_CARTESIAN_COORDS || this->m_ModelType==Superclass::SYMMETRICAL_TENSOR_IN_SPHERICAL_COORDS)
185  vec[0]=diffusivities[0], vec[1]=diffusivities[1], vec[2]=diffusivities[1];
186  else if (this->m_ModelType==Superclass::TENSOR_IN_EULER_ANGLES)
187  vec[0]=diffusivities[0], vec[1]=diffusivities[1], vec[2]=diffusivities[2];
188  tensor.Fill(0.0);
189  tensor.SetEigenValues(vec);
190  // utlPrintVar2(true, "before", tensor);
191  Matrix<double,3,3> rotation;
192  utl::RotationMatrixFromVectors<Vector<double,3>, Matrix<double,3> >(e1, principalDirection, rotation);
193  tensor.Rotate(rotation);
194 
195  // utlPrintVar2(true, "after", tensor);
196  tensorVec.push_back(tensor);
197  if (isOutputPeak)
198  {
199  PeakContainerHelper::SetPeak<Vector<double,3>, OutputImagePixelType >(principalDirection, pixelPeakInput, s, this->m_PeakType);
200  if (this->m_PeakType==NXYZV || this->m_PeakType==XYZV)
201  PeakContainerHelper::SetPeakValue<OutputImagePixelType >(partialVolumeWeight, pixelPeakInput, s, this->m_PeakType);
202  numberOfPeaks++;
203  }
204  }
205 
206  if (isOutputPeak && (this->m_PeakType==NXYZ || this->m_PeakType==NXYZV))
207  pixelPeakInput[0] = numberOfPeaks;
208 
209  utlSAException(sumPartialVolumeWeight<1e-9)(sumPartialVolumeWeight).msg("wrong partialVolumeWeight");
210  for ( int s = 0; s < numberOfParameterSets; s += 1 )
211  weightVec[s] /= sumPartialVolumeWeight;
212 
213 
214  // typename OutputImageType::Pointer outputPtr = this->GetOutput();
215  OutputImagePointer outputDWI=this->GetDWIImage(), outputODF=this->GetODFImage(), outputEAP=this->GetEAPImage(), outputPeak=this->GetPeakImage();
216  ScalarImagePointer b0Image = this->GetB0Image(), outputRTO = this->GetRTOImage(), outputMSD = this->GetMSDImage();
217  OutputImagePixelType pixelDWI, tempPixelDWI, pixelEAP, tempPixelEAP, pixelODF, tempPixelODF;
218  double pixelRTO=0, tempPixelRTO=0, pixelMSD=0, tempPixelMSD=0;
219  int numberOfComponentsPerPixel4DWI = this->GetNumberOfQSpaceSamples();
220  int numberOfComponentsPerPixel4EAP = this->GetNumberOfRSpaceSamples();
221  int numberOfComponentsPerPixel4ODF = this->GetNumberOfRSpaceSamples();
222  OutputImagePixelType zeroPixel;
223  if (this->m_IsOutputDWI)
224  {
225  pixelDWI.SetSize( numberOfComponentsPerPixel4DWI );
226  pixelDWI.Fill( 0 );
227  tempPixelDWI.SetSize( numberOfComponentsPerPixel4DWI );
228  tempPixelDWI.Fill( 0 );
229  }
230  if (this->m_IsOutputEAP)
231  {
232  pixelEAP.SetSize( numberOfComponentsPerPixel4EAP );
233  pixelEAP.Fill( 0 );
234  tempPixelEAP.SetSize( numberOfComponentsPerPixel4EAP );
235  tempPixelEAP.Fill( 0 );
236  }
237  if (this->m_IsOutputODF)
238  {
239  pixelODF.SetSize( numberOfComponentsPerPixel4ODF );
240  pixelODF.Fill( 0 );
241  tempPixelODF.SetSize( numberOfComponentsPerPixel4ODF );
242  tempPixelODF.Fill( 0 );
243  }
244 
245  ImageRegionIteratorWithIndex<OutputImageType> outputDWIIt, outputODFIt, outputEAPIt, outputPeakIt;
246  ImageRegionIteratorWithIndex<ScalarImageType> b0It, outputRTOIt, outputMSDIt;
247  if ( this->m_IsOutputDWI )
248  {
249  // outputDWI->Print(std::cout<<"outputDWI");
250  outputDWIIt = ImageRegionIteratorWithIndex<OutputImageType>( outputDWI, outputDWI->GetRequestedRegion() );
251  b0It = ImageRegionIteratorWithIndex<ScalarImageType>( b0Image, b0Image->GetRequestedRegion() );
252  }
253  if ( this->m_IsOutputODF )
254  outputODFIt = ImageRegionIteratorWithIndex<OutputImageType>( outputODF, outputODF->GetRequestedRegion() );
255  if ( this->m_IsOutputEAP )
256  outputEAPIt = ImageRegionIteratorWithIndex<OutputImageType>( outputEAP, outputEAP->GetRequestedRegion() );
257  if ( isOutputPeak )
258  outputPeakIt = ImageRegionIteratorWithIndex<OutputImageType>( outputPeak, outputPeak->GetRequestedRegion() );
259  if ( this->m_IsOutputRTO )
260  outputRTOIt = ImageRegionIteratorWithIndex<ScalarImageType>( outputRTO, outputRTO->GetRequestedRegion() );
261  if ( this->m_IsOutputMSD )
262  outputMSDIt = ImageRegionIteratorWithIndex<ScalarImageType>( outputMSD, outputMSD->GetRequestedRegion() );
263 
264 
265  b0It.GoToBegin();
266  outputDWIIt.GoToBegin();
267  outputODFIt.GoToBegin();
268  outputEAPIt.GoToBegin();
269  outputPeakIt.GoToBegin();
270  OutputImagePixelType pixel_noise=pixelDWI;
271  Matrix<double,3,3> rotation;
272  int indexGrad = 0;
273  std::vector<double> peak;
274  while ( !outputDWIIt.IsAtEnd() )
275  {
276  if (this->m_IsOutputDWI)
277  pixelDWI.Fill( 0 );
278  if (this->m_IsOutputODF)
279  pixelODF.Fill( 0 );
280  if (this->m_IsOutputEAP)
281  pixelEAP.Fill( 0 );
282  pixelMSD=0, pixelRTO=0;
283  if (isOutputPeak)
284  pixelPeak = pixelPeakInput;
285 
286  if (m_RandomType==UNIFORM)
287  {
288  std::vector<double> vv(3,0);
289  if (m_StoredOrientationMatrix->Rows()>0)
290  {
291  vv[0] = (*m_StoredOrientationMatrix)(indexGrad,0);
292  vv[1] = (*m_StoredOrientationMatrix)(indexGrad,1);
293  vv[2] = (*m_StoredOrientationMatrix)(indexGrad,2);
294  indexGrad++;
295  }
296  else
297  vv = utl::RandomPointInSphere(false);
298  for ( int i = 0; i < 3; i += 1 )
299  principalDirection[i] = vv[i];
300  utl::RotationMatrixFromVectors<Vector<double,3>, Matrix<double,3> >(e1, principalDirection, rotation);
301  // utl::PrintContainer(&e1[0], &e1[3], "e1");
302  // utl::PrintContainer(&principalDirection[0], &principalDirection[3], "principalDirection");
303  // utl::PrintMatrix(rotation, 3,3, "rotation");
304  }
305 
306  for ( int s = 0; s < numberOfParameterSets; s += 1 )
307  {
308  tensor = tensorVec[s];
309  partialVolumeWeight = weightVec[s];
310  if (isOutputPeak)
311  peak = PeakContainerHelper::GetPeak(pixelPeakInput, s, this->m_PeakType);
312  // tensor.Print(std::cout<<"tensor 0=");
313 
314  if (m_RandomType==UNIFORM)
315  {
316  tensor.Rotate(rotation);
317  if (isOutputPeak)
318  {
319  peak = PeakContainerHelper::GetPeak(pixelPeakInput, s, this->m_PeakType);
320  vnl_vector<double> peakVnl = utl::StdVectorToVnlVector(peak);
321  vnl_vector<double> peakRotated = rotation.GetVnlMatrix()*peakVnl;
322  PeakContainerHelper::SetPeak<vnl_vector<double>, OutputImagePixelType >(peakRotated, pixelPeak, s, this->m_PeakType);
323  }
324  }
325  // tensor.Print(std::cout<<"tensor 1=");
326 
327  if ( this->m_IsOutputDWI )
328  {
329  tensor.GetDWISamples(tempPixelDWI, *qSpaceOrientationMatrix, *bVector);
330  pixelDWI += tempPixelDWI*this->m_B0Scale*partialVolumeWeight;
331  }
332  if ( this->m_IsOutputODF )
333  {
334  tensor.GetODFSamples(tempPixelODF, *rSpaceOrientationMatrix,this->m_ODFOrder, false);
335  pixelODF += tempPixelODF*partialVolumeWeight;
336  }
337  if ( this->m_IsOutputEAP)
338  {
339  tensor.GetEAPSamples(tempPixelEAP, *rSpaceOrientationMatrix, *rVector, tau);
340  pixelEAP += tempPixelEAP*this->m_B0Scale*partialVolumeWeight;
341  }
342  if ( this->m_IsOutputRTO)
343  {
344  tempPixelRTO = tensor.GetReturnToOrigin(tau);
345  pixelRTO += tempPixelRTO*partialVolumeWeight;
346  }
347  if ( this->m_IsOutputMSD)
348  {
349  tempPixelMSD = tensor.GetMeanSquaredDisplacement(tau);
350  pixelMSD += tempPixelMSD*partialVolumeWeight;
351  }
352  }
353 
354  if (this->m_IsOutputODF && this->m_ODFOrder!=2)
355  {
356  PrecisionType normFactor = 4*vnl_math::pi*utl::GetSumOfVector<OutputImagePixelType>(pixelODF,pixelODF.Size()) / pixelODF.Size();
357  if (normFactor!=0)
358  pixelODF /= normFactor;
359  }
360 
361  if ( this->m_IsOutputDWI )
362  {
363  if ( (this->m_NoiseSigma>0 || this->m_SNR>0) && pixelDWI.GetSquaredNorm()>0 )
364  {
365  PrecisionType sigma_real = -1;
366  if ( this->m_SNR > 0 && this->m_NoiseSigma <= 0)
367  {
368  double mean = 0;
369  for ( unsigned int k=0; k<numberOfComponentsPerPixel4DWI; k++ )
370  mean += pixelDWI[k];
371  mean /= numberOfComponentsPerPixel4DWI;
372  sigma_real = mean/this->m_SNR;
373  }
374  else if ( this->m_SNR <= 0 && this->m_NoiseSigma > 0)
375  sigma_real = this->m_NoiseSigma;
376  pixelDWI = utl::AddNoise<OutputImagePixelType>(pixelDWI, pixelDWI.GetSize(), sigma_real, true);
377  }
378  outputDWIIt.Set(pixelDWI);
379  b0It.Set(this->m_B0Scale);
380  ++b0It;
381  ++outputDWIIt;
382  }
383 
384  if ( this->m_IsOutputODF )
385  {
386  outputODFIt.Set(pixelODF);
387  ++outputODFIt;
388  }
389 
390  if ( this->m_IsOutputEAP )
391  {
392  outputEAPIt.Set(pixelEAP);
393  ++outputEAPIt;
394  }
395  if ( isOutputPeak )
396  {
397  outputPeakIt.Set(pixelPeak);
398  ++outputPeakIt;
399  }
400  if ( this->m_IsOutputRTO )
401  {
402  outputRTOIt.Set(pixelRTO);
403  ++outputRTOIt;
404  }
405  if ( this->m_IsOutputMSD )
406  {
407  outputMSDIt.Set(pixelMSD);
408  ++outputMSDIt;
409  }
410  }
411 }
412 
413 }
414 
415 
416 #endif
417 
418 
419 
NDArray is a N-Dimensional array class (row-major, c version)
Definition: utlFunctors.h:131
void GetEAPSamples(TVectorType &vec, const utl::NDArray< TPrecision, 2 > &gradients, const std::vector< TPrecision > &rValues, const TPrecision &tau) const
ImageSource< TOutputImage > Superclass
#define PrintVar1(cond, var, os)
Definition: utlCoreMacro.h:447
std::vector< PrecisionType > STDVectorType
void PrintUtlMatrix(const NDArray< T, 2 > &mat, const std::string &str="", const char *separate=" ", std::ostream &os=std::cout)
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)
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
static void sort(int *irOut, T *prOut, int beg, int end)
Definition: misc.h:141
Self & Rotate(const typename Self::MatrixType &mat)
LightObject::Pointer InternalClone() const
utl_shared_ptr< STDVectorType > STDVectorPointer
#define utlGlobalException(cond, expout)
Definition: utlCoreMacro.h:372
utl_shared_ptr< MatrixType > MatrixPointer
bool IsInt(const std::string &input, const double epss=1e-10)
Definition: utlCore.h:792
static VectorType GetPeak(const ContainerType &vec, const int index, const PeakType peakType)
std::vector< double > RandomPointInSphere(const bool hemis)
Definition: utlCore.h:1489
OutputImageType::PixelType OutputImagePixelType
ScalarImageType::Pointer ScalarImagePointer
static int GetNumberOfPeaks(const PeakType peakType, const int dimension)
void SetEigenValues(const TArrayType &array)
#define utlSAException(expr)
Definition: utlCoreMacro.h:543
Generate DWI data based on provided parameter file.
SmartPointer< Self > Pointer
void PrintSelf(std::ostream &os, Indent indent) const
OutputImageType::Pointer OutputImagePointer
vnl_vector< T > StdVectorToVnlVector(const std::vector< T > &vec)
Definition: utlVNL.h:355
#define utlShowPosition(cond)
Definition: utlCoreMacro.h:554
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