4 #include <itkImageRegionIterator.h> 5 #include <itkImageRegionIteratorWithIndex.h> 6 #include <itkImageRegionConstIterator.h> 8 #include <vnl/algo/vnl_determinant.h> 9 #include <vnl/algo/vnl_matrix_inverse.h> 10 #include <vnl/algo/vnl_symmetric_eigensystem.h> 20 template <
class InputPixelScalarType,
class OutputPixelScalarType>
25 if (i == m_GradientDirections.size())
26 m_GradientDirections.push_back(grad);
27 else if (i > m_GradientDirections.size())
28 std::cerr <<
"Trying to add a direction not contiguous... Add directions contiguously (0,1,2,3,...)..." << std::endl;
30 m_GradientDirections[i] = grad;
33 template <
class InputPixelScalarType,
class OutputPixelScalarType>
38 if (this->GetComputationMask())
41 typedef itk::ImageRegionConstIterator <InputImageType> B0IteratorType;
42 typedef itk::ImageRegionIterator <MaskImageType> MaskIteratorType;
44 unsigned int firstB0Index = 0;
45 while (m_BValuesList[firstB0Index] > 10)
48 B0IteratorType b0Itr(this->GetInput(firstB0Index),this->GetOutput()->GetLargestPossibleRegion());
50 typename MaskImageType::Pointer maskImage = MaskImageType::New();
51 maskImage->Initialize();
52 maskImage->SetRegions(this->GetInput(0)->GetLargestPossibleRegion());
53 maskImage->SetSpacing (this->GetInput(0)->GetSpacing());
54 maskImage->SetOrigin (this->GetInput(0)->GetOrigin());
55 maskImage->SetDirection (this->GetInput(0)->GetDirection());
56 maskImage->Allocate();
58 MaskIteratorType maskItr (maskImage,this->GetOutput()->GetLargestPossibleRegion());
59 while (!b0Itr.IsAtEnd())
61 if (b0Itr.Get() > m_B0Threshold)
70 this->SetComputationMask(maskImage);
73 template <
class InputPixelScalarType,
class OutputPixelScalarType>
81 this->Superclass::GenerateOutputInformation();
84 output->SetVectorLength(m_NumberOfComponents);
87 template <
class InputPixelScalarType,
class OutputPixelScalarType>
92 if (m_BValuesList.size() != this->GetNumberOfIndexedInputs())
94 std::string error(
"There should be the same number of input images and input b-values... ");
95 error += std::to_string (m_BValuesList.size());
97 error += std::to_string (this->GetNumberOfIndexedInputs());
99 throw itk::ExceptionObject(__FILE__, __LINE__,error,ITK_LOCATION);
102 itk::ImageRegionIterator <OutputImageType> fillOut(this->GetOutput(),this->GetOutput()->GetLargestPossibleRegion());
103 typename OutputImageType::PixelType emptyModelVec(m_NumberOfComponents);
104 emptyModelVec.Fill(0);
106 while (!fillOut.IsAtEnd())
108 fillOut.Set(emptyModelVec);
112 m_EstimatedB0Image = OutputB0ImageType::New();
113 m_EstimatedB0Image->Initialize();
114 m_EstimatedB0Image->SetOrigin(this->GetOutput()->GetOrigin());
115 m_EstimatedB0Image->SetSpacing(this->GetOutput()->GetSpacing());
116 m_EstimatedB0Image->SetDirection(this->GetOutput()->GetDirection());
117 m_EstimatedB0Image->SetRegions(this->GetOutput()->GetLargestPossibleRegion());
119 m_EstimatedB0Image->Allocate();
120 m_EstimatedB0Image->FillBuffer(0.0);
122 m_EstimatedVarianceImage = OutputB0ImageType::New();
123 m_EstimatedVarianceImage->Initialize();
124 m_EstimatedVarianceImage->SetOrigin(this->GetOutput()->GetOrigin());
125 m_EstimatedVarianceImage->SetSpacing(this->GetOutput()->GetSpacing());
126 m_EstimatedVarianceImage->SetDirection(this->GetOutput()->GetDirection());
127 m_EstimatedVarianceImage->SetRegions(this->GetOutput()->GetLargestPossibleRegion());
129 m_EstimatedVarianceImage->Allocate();
130 m_EstimatedVarianceImage->FillBuffer(0.0);
132 vnl_matrix <double> initSolverSystem(this->GetNumberOfIndexedInputs(),m_NumberOfComponents + 1);
133 initSolverSystem.fill(0.0);
134 for (
unsigned int i = 0;i < this->GetNumberOfIndexedInputs();++i)
136 initSolverSystem(i,0) = 1.0;
137 unsigned int pos = 1;
138 for (
unsigned int j = 0;j < 3;++j)
140 for (
unsigned int k = 0;k < j;++k)
142 initSolverSystem(i,pos) = - 2.0 * m_BValuesList[i] * m_GradientDirections[i][j] * m_GradientDirections[i][k];
146 initSolverSystem(i,pos) = - m_BValuesList[i] * m_GradientDirections[i][j] * m_GradientDirections[i][j];
151 vnl_matrix_inverse <double> inverter (initSolverSystem);
152 m_InitialMatrixSolver = inverter.pinverse();
154 Superclass::BeforeThreadedGenerateData();
157 template <
class InputPixelScalarType,
class OutputPixelScalarType>
162 typedef itk::ImageRegionConstIterator <InputImageType> ImageIteratorType;
164 unsigned int numInputs = this->GetNumberOfIndexedInputs();
165 std::vector <ImageIteratorType> inIterators;
166 for (
unsigned int i = 0;i < numInputs;++i)
167 inIterators.push_back(ImageIteratorType(this->GetInput(i),outputRegionForThread));
169 typedef itk::ImageRegionConstIterator <MaskImageType> MaskIteratorType;
170 MaskIteratorType maskIterator(this->GetComputationMask(),outputRegionForThread);
172 typedef itk::ImageRegionIterator <OutputImageType> OutImageIteratorType;
173 OutImageIteratorType outIterator(this->GetOutput(),outputRegionForThread);
175 typedef itk::ImageRegionIterator <OutputB0ImageType> OutB0ImageIteratorType;
176 OutB0ImageIteratorType outB0Iterator(m_EstimatedB0Image,outputRegionForThread);
177 OutB0ImageIteratorType outVarianceIterator(m_EstimatedVarianceImage,outputRegionForThread);
179 typedef typename OutputImageType::PixelType OutputPixelType;
180 std::vector <double> dwi (numInputs,0);
181 std::vector <double> lnDwi (numInputs,0);
182 OutputPixelType resVec(m_NumberOfComponents);
184 std::vector <double> predictedValues(numInputs,0);
190 typedef itk::SymmetricEigenAnalysis < vnl_matrix <double>, vnl_diag_matrix<double>, vnl_matrix <double> > EigenAnalysisType;
191 EigenAnalysisType eigen(3);
193 while (!outIterator.IsAtEnd())
197 if (maskIterator.Get() == 0)
199 outIterator.Set(resVec);
200 outB0Iterator.Set(0);
201 outVarianceIterator.Set(0);
203 for (
unsigned int i = 0;i < numInputs;++i)
209 ++outVarianceIterator;
214 for (
unsigned int i = 0;i < numInputs;++i)
216 dwi[i] = inIterators[i].Get();
217 lnDwi[i] = std::log(std::max(1.0e-6,dwi[i]));
220 for (
unsigned int i = 0;i < m_NumberOfComponents;++i)
222 for (
unsigned int j = 0;j < numInputs;++j)
223 resVec[i] += m_InitialMatrixSolver(i+1,j) * lnDwi[j];
232 std::vector <double> optimizedValue(m_NumberOfComponents, 0.0);
234 double cThetaControl = 0;
235 for (
unsigned int i = 0;i < 3;++i)
238 if (std::abs(cThetaControl + 1.0) > 1.0e-5)
241 double minValue = 1.0e-7;
242 double maxValue = 1.0e-2;
243 for (
unsigned int i = 0;i < 3;++i)
245 optimizedValue[i] += M_PI;
246 int num2Pi = std::floor(optimizedValue[i] / (2.0 * M_PI));
247 optimizedValue[i] -= 2.0 * M_PI * num2Pi + M_PI;
249 optimizedValue[i + 3] = std::min(maxValue, std::max(data.
workEigenValues[i], minValue));
253 nlopt::opt opt(nlopt::LN_BOBYQA, m_NumberOfComponents);
255 std::vector <double> lowerBounds(m_NumberOfComponents, - M_PI);
256 for (
unsigned int i = 0;i < 3;++i)
257 lowerBounds[i + 3] = minValue;
259 opt.set_lower_bounds(lowerBounds);
261 std::vector <double> upperBounds(m_NumberOfComponents, M_PI);
262 for (
unsigned int i = 0;i < 3;++i)
263 upperBounds[i + 3] = maxValue;
265 opt.set_upper_bounds(upperBounds);
266 opt.set_xtol_rel(1e-4);
267 opt.set_ftol_rel(1e-4);
268 opt.set_maxeval(2500);
276 opt.set_min_objective(OptimizationFunction, &data);
280 opt.optimize(optimizedValue, minf);
282 catch(nlopt::roundoff_limited& e)
284 bool failedOpt =
false;
285 for (
unsigned int i = 0;i < 6;++i)
287 if (!std::isfinite(optimizedValue[i]))
297 outIterator.Set(resVec);
298 outB0Iterator.Set(0);
299 outVarianceIterator.Set(0);
301 for (
unsigned int i = 0;i < numInputs;++i)
307 ++outVarianceIterator;
308 this->IncrementNumberOfProcessedPoints();
315 for (
unsigned int i = 0;i < 3;++i)
321 double outVarianceValue;
322 double outB0Value = this->ComputeB0AndVarianceFromTensorVector(data.
workTensor,dwi,outVarianceValue);
324 outIterator.Set(resVec);
325 outB0Iterator.Set(outB0Value);
326 outVarianceIterator.Set(outVarianceValue);
328 for (
unsigned int i = 0;i < numInputs;++i)
331 this->IncrementNumberOfProcessedPoints();
335 ++outVarianceIterator;
339 template <
class InputPixelScalarType,
class OutputPixelScalarType>
344 double sumSquaredObservedSignals = 0;
345 double sumPredictedPerObservedSignals = 0;
346 double sumSquaredPredictedSignals = 0;
347 unsigned int numInputs = dwiSignal.size();
349 for (
unsigned int i = 0;i < numInputs;++i)
351 double observedValue = dwiSignal[i];
352 double predictedValue = 0;
353 double bValue = m_BValuesList[i];
359 for (
unsigned int j = 0;j < 3;++j)
361 predictedValue += tensorValue(j,j) * m_GradientDirections[i][j] * m_GradientDirections[i][j];
362 for (
unsigned int k = j + 1;k < 3;++k)
363 predictedValue += 2.0 * tensorValue(j,k) * m_GradientDirections[i][j] * m_GradientDirections[i][k];
366 predictedValue = std::exp(- bValue * predictedValue);
369 sumSquaredObservedSignals += observedValue * observedValue;
370 sumSquaredPredictedSignals += predictedValue * predictedValue;
371 sumPredictedPerObservedSignals += predictedValue * observedValue;
374 double outB0Value = sumPredictedPerObservedSignals / sumSquaredPredictedSignals;
375 outVarianceValue = (sumSquaredObservedSignals - 2.0 * outB0Value * sumPredictedPerObservedSignals + outB0Value * outB0Value * sumSquaredPredictedSignals) / (numInputs - 1.0);
379 template <
class InputPixelScalarType,
class OutputPixelScalarType>
389 template <
class InputPixelScalarType,
class OutputPixelScalarType>
393 std::vector <double> &predictedValues, vnl_matrix <double> &rotationMatrix,
394 vnl_matrix<double> &workTensor, vnl_diag_matrix <double> &workEigenValues)
397 workEigenValues.set_size(3);
398 for (
unsigned int i = 0;i < 3;++i)
399 workEigenValues[i] = x[3 + i];
401 workTensor.set_size(3,3);
404 for (
unsigned int i = 0;i < this->GetNumberOfIndexedInputs();++i)
406 predictedValues[i] = 0;
407 double bValue = m_BValuesList[i];
410 predictedValues[i] = 1;
414 for (
unsigned int j = 0;j < 3;++j)
416 predictedValues[i] += workTensor(j,j) * m_GradientDirections[i][j] * m_GradientDirections[i][j];
417 for (
unsigned int k = j + 1;k < 3;++k)
418 predictedValues[i] += 2.0 * workTensor(j,k) * m_GradientDirections[i][j] * m_GradientDirections[i][k];
421 predictedValues[i] = std::exp(- bValue * predictedValues[i]);
425 double normConstant = 0;
427 for (
unsigned int i = 0;i < this->GetNumberOfIndexedInputs();++i)
429 b0Val += predictedValues[i] * observedData[i];
430 normConstant += predictedValues[i] * predictedValues[i];
433 b0Val /= normConstant;
435 double costValue = 0;
436 for (
unsigned int i = 0;i < this->GetNumberOfIndexedInputs();++i)
438 double predVal = b0Val * predictedValues[i];
439 costValue += (predVal - observedData[i]) * (predVal - observedData[i]);
void AddGradientDirection(unsigned int i, vnl_vector_fixed< double, 3 > &grad)
double ComputeB0AndVarianceFromTensorVector(const vnl_matrix< double > &tensorValue, const std::vector< double > &dwiSignal, double &outVarianceValue)
void Get3DRotationLogarithm(const vnl_matrix< T > &rotationMatrix, std::vector< T > &outputAngles)
Computation of a 3D rotation matrix logarithm. Rodrigues' explicit formula.
Superclass::OutputImageRegionType OutputImageRegionType
itk::VectorImage< OutputPixelScalarType, 3 > OutputImageType
void BeforeThreadedGenerateData() ITK_OVERRIDE
vnl_matrix< double > workTensor
static double OptimizationFunction(const std::vector< double > &x, std::vector< double > &grad, void *func_data)
void DynamicThreadedGenerateData(const OutputImageRegionType &outputRegionForThread) ITK_OVERRIDE
void GenerateOutputInformation() ITK_OVERRIDE
vnl_diag_matrix< double > workEigenValues
void GetVectorRepresentation(const vnl_matrix< T1 > &tensor, itk::VariableLengthVector< T2 > &vector, unsigned int vecDim=0, bool scale=false)
void CheckComputationMask() ITK_OVERRIDE
vnl_matrix< double > rotationMatrix
std::vector< double > predictedValues
void GetTensorFromVectorRepresentation(const itk::VariableLengthVector< T1 > &vector, vnl_matrix< T2 > &tensor, unsigned int tensDim=0, bool scale=false)
void RecomposeTensor(vnl_diag_matrix< T1 > &eigs, vnl_matrix< T1 > &eigVecs, vnl_matrix< T2 > &resMatrix)
Recompose tensor from values extracted using SymmetricEigenAnalysis (vnl_symmetric_eigensystem transp...
void Get3DRotationExponential(const std::vector< T > &angles, vnl_matrix< T > &outputRotation)
Computation of a 3D rotation matrix exponential. Rodrigues' explicit formula.
std::vector< double > dwi
double ComputeCostAtPosition(const std::vector< double > &x, const std::vector< double > &observedData, std::vector< double > &predictedValues, vnl_matrix< double > &rotationMatrix, vnl_matrix< double > &workTensor, vnl_diag_matrix< double > &workEigenValues)