18 #ifndef __itkDWISingleVoxelGenerator_hxx 19 #define __itkDWISingleVoxelGenerator_hxx 28 template <
class TOutputImage,
class TScalarImage>
37 template <
class TOutputImage,
class TScalarImage>
41 Superclass::PrintSelf(os, indent);
42 PrintVar1(
true, m_RandomType, os<<indent);
43 if (m_StoredOrientationMatrix->Rows()>0)
47 template <
class TOutputImage,
class TScalarImage>
48 typename LightObject::Pointer
52 typename LightObject::Pointer loPtr = Superclass::InternalClone();
57 itkExceptionMacro(<<
"downcast to type " << this->GetNameOfClass()<<
" failed.");
60 rval->m_DiffusionParameterValues = m_DiffusionParameterValues;
61 rval->m_RandomType = m_RandomType;
62 rval->m_StoredOrientationMatrix = m_StoredOrientationMatrix;
66 template <
class TOutputImage,
class TScalarImage>
70 Superclass::Initialization();
71 if (m_RandomType==UNIFORM && m_StoredOrientationMatrix->Rows())
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();
79 template <
class TOutputImage,
class TScalarImage>
84 this->Initialization();
86 bool isOutputPeak = this->m_MaxNumberOfPeaks>0;
88 MatrixPointer qSpaceOrientationMatrix = this->m_SamplingSchemeQSpace->GetOrientationsCartesian();
90 MatrixPointer rSpaceOrientationMatrix = this->m_SamplingSchemeRSpace->GetOrientationsCartesian();
92 double tau = this->m_SamplingSchemeQSpace->GetTau();
95 this->AllocateOutputs();
97 int numberOfDiffusivityComponentsPerParameterSet = -1;
98 unsigned int numberOfOrientationComponentsPerParameterSet = 3;
99 if (this->m_ModelType==Superclass::SYMMETRICAL_TENSOR_IN_CARTESIAN_COORDS)
101 numberOfDiffusivityComponentsPerParameterSet=2;
102 numberOfOrientationComponentsPerParameterSet = 3;
104 else if (this->m_ModelType==Superclass::SYMMETRICAL_TENSOR_IN_SPHERICAL_COORDS)
106 numberOfDiffusivityComponentsPerParameterSet=2;
107 numberOfOrientationComponentsPerParameterSet = 2;
109 else if (this->m_ModelType==Superclass::TENSOR_IN_EULER_ANGLES)
111 numberOfDiffusivityComponentsPerParameterSet=3;
112 numberOfOrientationComponentsPerParameterSet = 3;
115 else if (this->m_ModelType==Superclass::CYLINDER_SPHERICAL_MODEL)
117 numberOfOrientationComponentsPerParameterSet = 2;
121 int numberOfComponentsPerParameterSet = numberOfDiffusivityComponentsPerParameterSet + numberOfOrientationComponentsPerParameterSet+1;
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;
128 std::vector<DiffusionTensor<double> > tensorVec;
129 std::vector<double> weightVec;
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;
137 STDVectorType diffusivities(numberOfDiffusivityComponentsPerParameterSet);
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);
151 for (
unsigned int s=0; s<numberOfParameterSets; s++ )
153 for (
unsigned int d=0; d<numberOfOrientationComponentsPerParameterSet; d++ )
154 principalDirection[d] = m_DiffusionParameterValues[s* numberOfComponentsPerParameterSet + d];
156 norm = principalDirection.GetNorm();
157 if (this->m_ModelType==Superclass::SYMMETRICAL_TENSOR_IN_SPHERICAL_COORDS || this->m_ModelType==Superclass::CYLINDER_SPHERICAL_MODEL)
159 principalDirection *= vnl_math::pi/180.0;
160 double theta=principalDirection[0], phi=principalDirection[1];
165 for (
unsigned int d=0; d<numberOfOrientationComponentsPerParameterSet; d++ )
166 principalDirection[d] /= norm;
169 for (
unsigned int d=0; d<numberOfDiffusivityComponentsPerParameterSet; d++ )
171 diffusivities[d] = m_DiffusionParameterValues[s*numberOfComponentsPerParameterSet + numberOfOrientationComponentsPerParameterSet + d ];
173 std::sort(diffusivities.begin(), diffusivities.end(), std::greater<double>());
175 partialVolumeWeight = m_DiffusionParameterValues[s*numberOfComponentsPerParameterSet + numberOfOrientationComponentsPerParameterSet + numberOfDiffusivityComponentsPerParameterSet];
176 sumPartialVolumeWeight += partialVolumeWeight;
177 weightVec.push_back(partialVolumeWeight);
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];
191 Matrix<double,3,3> rotation;
192 utl::RotationMatrixFromVectors<Vector<double,3>, Matrix<double,3> >(e1, principalDirection, rotation);
196 tensorVec.push_back(tensor);
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);
206 if (isOutputPeak && (this->m_PeakType==
NXYZ || this->m_PeakType==
NXYZV))
207 pixelPeakInput[0] = numberOfPeaks;
209 utlSAException(sumPartialVolumeWeight<1e-9)(sumPartialVolumeWeight).msg(
"wrong partialVolumeWeight");
210 for (
int s = 0; s < numberOfParameterSets; s += 1 )
211 weightVec[s] /= sumPartialVolumeWeight;
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();
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();
223 if (this->m_IsOutputDWI)
225 pixelDWI.SetSize( numberOfComponentsPerPixel4DWI );
227 tempPixelDWI.SetSize( numberOfComponentsPerPixel4DWI );
228 tempPixelDWI.Fill( 0 );
230 if (this->m_IsOutputEAP)
232 pixelEAP.SetSize( numberOfComponentsPerPixel4EAP );
234 tempPixelEAP.SetSize( numberOfComponentsPerPixel4EAP );
235 tempPixelEAP.Fill( 0 );
237 if (this->m_IsOutputODF)
239 pixelODF.SetSize( numberOfComponentsPerPixel4ODF );
241 tempPixelODF.SetSize( numberOfComponentsPerPixel4ODF );
242 tempPixelODF.Fill( 0 );
245 ImageRegionIteratorWithIndex<OutputImageType> outputDWIIt, outputODFIt, outputEAPIt, outputPeakIt;
246 ImageRegionIteratorWithIndex<ScalarImageType> b0It, outputRTOIt, outputMSDIt;
247 if ( this->m_IsOutputDWI )
250 outputDWIIt = ImageRegionIteratorWithIndex<OutputImageType>( outputDWI, outputDWI->GetRequestedRegion() );
251 b0It = ImageRegionIteratorWithIndex<ScalarImageType>( b0Image, b0Image->GetRequestedRegion() );
253 if ( this->m_IsOutputODF )
254 outputODFIt = ImageRegionIteratorWithIndex<OutputImageType>( outputODF, outputODF->GetRequestedRegion() );
255 if ( this->m_IsOutputEAP )
256 outputEAPIt = ImageRegionIteratorWithIndex<OutputImageType>( outputEAP, outputEAP->GetRequestedRegion() );
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() );
266 outputDWIIt.GoToBegin();
267 outputODFIt.GoToBegin();
268 outputEAPIt.GoToBegin();
269 outputPeakIt.GoToBegin();
271 Matrix<double,3,3> rotation;
273 std::vector<double> peak;
274 while ( !outputDWIIt.IsAtEnd() )
276 if (this->m_IsOutputDWI)
278 if (this->m_IsOutputODF)
280 if (this->m_IsOutputEAP)
282 pixelMSD=0, pixelRTO=0;
284 pixelPeak = pixelPeakInput;
286 if (m_RandomType==UNIFORM)
288 std::vector<double> vv(3,0);
289 if (m_StoredOrientationMatrix->Rows()>0)
291 vv[0] = (*m_StoredOrientationMatrix)(indexGrad,0);
292 vv[1] = (*m_StoredOrientationMatrix)(indexGrad,1);
293 vv[2] = (*m_StoredOrientationMatrix)(indexGrad,2);
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);
306 for (
int s = 0; s < numberOfParameterSets; s += 1 )
308 tensor = tensorVec[s];
309 partialVolumeWeight = weightVec[s];
314 if (m_RandomType==UNIFORM)
321 vnl_vector<double> peakRotated = rotation.GetVnlMatrix()*peakVnl;
322 PeakContainerHelper::SetPeak<vnl_vector<double>,
OutputImagePixelType >(peakRotated, pixelPeak, s, this->m_PeakType);
327 if ( this->m_IsOutputDWI )
329 tensor.
GetDWISamples(tempPixelDWI, *qSpaceOrientationMatrix, *bVector);
330 pixelDWI += tempPixelDWI*this->m_B0Scale*partialVolumeWeight;
332 if ( this->m_IsOutputODF )
334 tensor.
GetODFSamples(tempPixelODF, *rSpaceOrientationMatrix,this->m_ODFOrder,
false);
335 pixelODF += tempPixelODF*partialVolumeWeight;
337 if ( this->m_IsOutputEAP)
339 tensor.
GetEAPSamples(tempPixelEAP, *rSpaceOrientationMatrix, *rVector, tau);
340 pixelEAP += tempPixelEAP*this->m_B0Scale*partialVolumeWeight;
342 if ( this->m_IsOutputRTO)
345 pixelRTO += tempPixelRTO*partialVolumeWeight;
347 if ( this->m_IsOutputMSD)
350 pixelMSD += tempPixelMSD*partialVolumeWeight;
354 if (this->m_IsOutputODF && this->m_ODFOrder!=2)
356 PrecisionType normFactor = 4*vnl_math::pi*utl::GetSumOfVector<OutputImagePixelType>(pixelODF,pixelODF.Size()) / pixelODF.Size();
358 pixelODF /= normFactor;
361 if ( this->m_IsOutputDWI )
363 if ( (this->m_NoiseSigma>0 || this->m_SNR>0) && pixelDWI.GetSquaredNorm()>0 )
366 if ( this->m_SNR > 0 && this->m_NoiseSigma <= 0)
369 for (
unsigned int k=0; k<numberOfComponentsPerPixel4DWI; k++ )
371 mean /= numberOfComponentsPerPixel4DWI;
372 sigma_real = mean/this->m_SNR;
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);
378 outputDWIIt.Set(pixelDWI);
379 b0It.Set(this->m_B0Scale);
384 if ( this->m_IsOutputODF )
386 outputODFIt.Set(pixelODF);
390 if ( this->m_IsOutputEAP )
392 outputEAPIt.Set(pixelEAP);
397 outputPeakIt.Set(pixelPeak);
400 if ( this->m_IsOutputRTO )
402 outputRTOIt.Set(pixelRTO);
405 if ( this->m_IsOutputMSD )
407 outputMSDIt.Set(pixelMSD);
NDArray is a N-Dimensional array class (row-major, c version)
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)
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)
RealValueType GetReturnToOrigin(const TPrecision tau) const
#define utlSAGlobalException(expr)
static void sort(int *irOut, T *prOut, int beg, int end)
Self & Rotate(const typename Self::MatrixType &mat)
LightObject::Pointer InternalClone() const
utl_shared_ptr< STDVectorType > STDVectorPointer
#define utlGlobalException(cond, expout)
utl_shared_ptr< MatrixType > MatrixPointer
bool IsInt(const std::string &input, const double epss=1e-10)
static VectorType GetPeak(const ContainerType &vec, const int index, const PeakType peakType)
std::vector< double > RandomPointInSphere(const bool hemis)
OutputImageType::PixelType OutputImagePixelType
ScalarImageType::Pointer ScalarImagePointer
static int GetNumberOfPeaks(const PeakType peakType, const int dimension)
void SetEigenValues(const TArrayType &array)
#define utlSAException(expr)
Generate DWI data based on provided parameter file.
DWISingleVoxelGenerator()
SmartPointer< Self > Pointer
void PrintSelf(std::ostream &os, Indent indent) const
OutputImageType::Pointer OutputImagePointer
vnl_vector< T > StdVectorToVnlVector(const std::vector< T > &vec)
#define utlShowPosition(cond)
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