4 #include <itkImageRegionIterator.h> 5 #include <itkImageRegionConstIterator.h> 15 template <
class TPixelScalarType>
17 MultiT2RelaxometryEstimationImageFilter <TPixelScalarType>
18 ::BeforeThreadedGenerateData()
20 if (m_NLEstimation && (!m_InitialT2Map || !m_InitialM0Map))
21 itkExceptionMacro(
"Missing inputs for non-local estimation");
23 Superclass::BeforeThreadedGenerateData();
25 m_T2CompartmentValues.resize(m_NumberOfT2Compartments);
26 double logStart = std::log(m_LowerT2Bound);
27 double logEnd = std::log(m_UpperT2Bound);
28 double step = (logEnd - logStart) / (m_NumberOfT2Compartments - 1.0);
30 double logValue = logStart;
31 m_T2CompartmentValues[0] = std::exp(logValue);
32 for (
unsigned int i = 1;i < m_NumberOfT2Compartments;++i)
35 m_T2CompartmentValues[i] = std::exp(logValue);
38 this->GetB1OutputImage()->FillBuffer(1.0);
40 m_T2OutputImage = VectorOutputImageType::New();
41 m_T2OutputImage->Initialize();
42 m_T2OutputImage->SetRegions(this->GetInput(0)->GetLargestPossibleRegion());
43 m_T2OutputImage->SetSpacing (this->GetInput(0)->GetSpacing());
44 m_T2OutputImage->SetOrigin (this->GetInput(0)->GetOrigin());
45 m_T2OutputImage->SetDirection (this->GetInput(0)->GetDirection());
46 m_T2OutputImage->SetVectorLength(m_NumberOfT2Compartments);
47 m_T2OutputImage->Allocate();
51 m_T2OutputImage->FillBuffer(zero);
55 this->PrepareNLPatchSearchers();
57 typedef itk::ImageRegionConstIterator <InputImageType> InIteratorType;
58 typedef itk::ImageRegionIterator <VectorOutputImageType> VectorIteratorType;
60 InIteratorType m0Iterator(m_InitialM0Map, m_InitialM0Map->GetLargestPossibleRegion());
61 VectorIteratorType t2Iterator(m_InitialT2Map, m_InitialT2Map->GetLargestPossibleRegion());
64 while (!m0Iterator.IsAtEnd())
66 tmpData = t2Iterator.Get();
67 double m0Value = m0Iterator.Get();
69 for (
unsigned int i = 0;i < m_NumberOfT2Compartments;++i)
70 tmpData[i] *= m0Value;
72 t2Iterator.Set(tmpData);
79 template <
class TPixelScalarType>
82 ::CheckComputationMask()
84 if (this->GetComputationMask())
87 typedef itk::ImageRegionConstIterator <InputImageType> IteratorType;
88 typedef itk::ImageRegionIterator <MaskImageType> MaskIteratorType;
90 std::vector <IteratorType> inItrs(this->GetNumberOfIndexedInputs());
91 for (
unsigned int i = 0;i < this->GetNumberOfIndexedInputs();++i)
92 inItrs[i] = IteratorType(this->GetInput(i),this->GetOutput(0)->GetLargestPossibleRegion());
94 typename MaskImageType::Pointer maskImage = MaskImageType::New();
95 maskImage->Initialize();
96 maskImage->SetRegions(this->GetInput(0)->GetLargestPossibleRegion());
97 maskImage->SetSpacing (this->GetInput(0)->GetSpacing());
98 maskImage->SetOrigin (this->GetInput(0)->GetOrigin());
99 maskImage->SetDirection (this->GetInput(0)->GetDirection());
100 maskImage->Allocate();
102 MaskIteratorType maskItr (maskImage,this->GetOutput(0)->GetLargestPossibleRegion());
103 while (!maskItr.IsAtEnd())
105 double averageVal = 0;
106 for (
unsigned int i = 0;i < this->GetNumberOfIndexedInputs();++i)
107 averageVal += inItrs[i].Get();
109 averageVal /= this->GetNumberOfIndexedInputs();
111 bool maskPoint = (averageVal <= m_AverageSignalThreshold);
118 for (
unsigned int i = 0;i < this->GetNumberOfIndexedInputs();++i)
124 this->SetComputationMask(maskImage);
127 template <
class TPixelScalarType>
132 typedef itk::ImageRegionConstIterator <InputImageType> ImageConstIteratorType;
133 typedef itk::ImageRegionIterator <InputImageType> ImageIteratorType;
134 typedef itk::ImageRegionIterator <VectorOutputImageType> VectorImageIteratorType;
136 VectorImageIteratorType outT2Iterator(this->GetT2OutputImage(),outputRegionForThread);
137 ImageIteratorType outM0Iterator(this->GetM0OutputImage(),outputRegionForThread);
138 ImageIteratorType outMWFIterator(this->GetMWFOutputImage(),outputRegionForThread);
139 ImageIteratorType outB1Iterator(this->GetB1OutputImage(),outputRegionForThread);
140 ImageIteratorType outCostIterator(this->GetCostOutputImage(),outputRegionForThread);
142 VectorImageIteratorType initT2Iterator;
143 ImageIteratorType initM0Iterator;
144 ImageIteratorType initB1Iterator;
148 initT2Iterator = VectorImageIteratorType(m_InitialT2Map, outputRegionForThread);
149 initM0Iterator = ImageIteratorType(m_InitialM0Map, outputRegionForThread);
150 initB1Iterator = ImageIteratorType(m_InitialB1Map, outputRegionForThread);
153 unsigned int numInputs = this->GetNumberOfIndexedInputs();
154 std::vector <ImageConstIteratorType> inIterators;
155 for (
unsigned int i = 0;i < numInputs;++i)
156 inIterators.push_back(ImageConstIteratorType(this->GetInput(i),outputRegionForThread));
158 typedef itk::ImageRegionIterator <MaskImageType> MaskIteratorType;
159 MaskIteratorType maskItr(this->GetComputationMask(),outputRegionForThread);
161 ImageIteratorType t1MapItr;
163 t1MapItr = ImageIteratorType(m_T1Map,outputRegionForThread);
166 T2VectorType signalValuesExtended(numInputs + m_NumberOfT2Compartments);
167 T2VectorType t2OptimizedWeights(m_NumberOfT2Compartments);
168 T2VectorType priorDistribution(m_NumberOfT2Compartments);
169 priorDistribution.fill(0);
174 DataMatrixType AMatrixExtended(numInputs + m_NumberOfT2Compartments,m_NumberOfT2Compartments,0);
188 typename B1CostFunctionType::Pointer cost = B1CostFunctionType::New();
189 cost->SetEchoSpacing(m_EchoSpacing);
190 cost->SetExcitationFlipAngle(m_T2ExcitationFlipAngle);
191 cost->SetT2FlipAngles(m_T2FlipAngles);
192 cost->SetM0Value(1.0);
194 unsigned int dimension = cost->GetNumberOfParameters();
195 itk::Array<double> lowerBounds(dimension);
196 itk::Array<double> upperBounds(dimension);
197 B1OptimizerType::ParametersType p(dimension);
202 std::vector <double> workDataWeights;
203 std::vector <OutputVectorType> workDataSamples;
205 unsigned int threadId = this->GetSafeThreadId();
207 while (!maskItr.IsAtEnd())
209 outputT2Weights.Fill(0);
211 if (maskItr.Get() == 0)
213 outT2Iterator.Set(outputT2Weights);
214 outM0Iterator.Set(0);
215 outMWFIterator.Set(0.0);
216 outB1Iterator.Set(0.0);
217 outCostIterator.Set(0.0);
226 for (
unsigned int i = 0;i < numInputs;++i)
242 double previousB1Value = -1.0;
243 double t1Value = 1000;
244 double m0Value = 0.0;
246 for (
unsigned int i = 0;i < numInputs;++i)
248 signalValues[i] = inIterators[i].Get();
249 signalValuesExtended[i] = signalValues[i];
254 t1Value = t1MapItr.Get();
259 cost->SetT1Value(t1Value);
260 cost->SetT2Values(m_T2CompartmentValues);
261 cost->SetT2RelaxometrySignals(signalValues);
266 b1Value = initB1Iterator.Get();
268 b1Value = outB1Iterator.Get();
272 outputT2Weights = initT2Iterator.Get();
273 this->ComputeTikhonovPrior(maskItr.GetIndex(),outputT2Weights,m_NLPatchSearchers[threadId],priorDistribution,
274 workDataWeights, workDataSamples);
277 unsigned int numGlobalIterations = 0;
280 while ((std::abs(b1Value - previousB1Value) > m_B1Tolerance)&&(numGlobalIterations < 100))
282 ++numGlobalIterations;
283 previousB1Value = b1Value;
286 AMatrixExtended.fill(0);
287 for (
unsigned int i = 0;i < m_NumberOfT2Compartments;++i)
289 epgSignalValues = epgSimulator.
GetValue(t1Value,m_T2CompartmentValues[i],b1Value,1.0);
290 for (
unsigned int j = 0;j < numInputs;++j)
292 AMatrix(j,i) = epgSignalValues[j];
293 AMatrixExtended(j,i) = epgSignalValues[j];
298 nnlsOpt->SetDataMatrix(AMatrix);
299 nnlsOpt->SetPoints(signalValues);
301 nnlsOpt->StartOptimization();
302 t2OptimizedWeights = nnlsOpt->GetCurrentPosition();
303 residual = nnlsOpt->GetCurrentResidual();
306 double normT2Weights = 0;
307 for (
unsigned int i = 0;i < m_NumberOfT2Compartments;++i)
308 normT2Weights += (t2OptimizedWeights[i] - priorDistribution[i]) * (t2OptimizedWeights[i] - priorDistribution[i]);
310 if (m_RegularizationIntensity > 1.0)
312 double lambdaSq = (m_RegularizationIntensity - 1.0) * residual / normT2Weights;
313 residual = this->ComputeTikhonovRegularizedSolution(nnlsOpt,AMatrixExtended,signalValuesExtended,
314 lambdaSq,priorDistribution,t2OptimizedWeights);
317 outCostIterator.Set(residual);
320 for (
unsigned int i = 0;i < m_NumberOfT2Compartments;++i)
321 m0Value += t2OptimizedWeights[i];
323 for (
unsigned int i = 0;i < m_NumberOfT2Compartments;++i)
324 outputT2Weights[i] = t2OptimizedWeights[i] / m0Value;
327 cost->SetT2Weights(t2OptimizedWeights);
329 B1OptimizerType::Pointer b1Optimizer = B1OptimizerType::New();
330 b1Optimizer->SetAlgorithm(NLOPT_LN_BOBYQA);
331 b1Optimizer->SetXTolRel(1.0e-4);
332 b1Optimizer->SetFTolRel(1.0e-6);
333 b1Optimizer->SetMaxEval(500);
334 b1Optimizer->SetVectorStorageSize(2000);
336 b1Optimizer->SetLowerBoundParameters(lowerBounds);
337 b1Optimizer->SetUpperBoundParameters(upperBounds);
339 b1Optimizer->SetInitialPosition(p);
340 b1Optimizer->SetMaximize(
false);
341 b1Optimizer->SetCostFunction(cost);
343 b1Optimizer->StartOptimization();
344 p = b1Optimizer->GetCurrentPosition();
348 outM0Iterator.Set(m0Value);
349 outT2Iterator.Set(outputT2Weights);
350 outCostIterator.Set(residual);
352 for (
unsigned int i = 0;i < m_NumberOfT2Compartments;++i)
354 if (m_T2CompartmentValues[i] > m_MyelinThreshold)
357 mwfValue += outputT2Weights[i];
360 outMWFIterator.Set(mwfValue);
361 outB1Iterator.Set(b1Value);
363 this->IncrementNumberOfProcessedPoints();
371 for (
unsigned int i = 0;i < numInputs;++i)
385 this->SafeReleaseThreadId(threadId);
388 template <
class TPixelScalarType>
391 ::PrepareNLPatchSearchers()
393 unsigned int numThreads = this->GetNumberOfWorkUnits();
394 m_NLPatchSearchers.resize(numThreads);
395 unsigned int maxAbsDisp = floor((
double)(m_SearchNeighborhood / m_SearchStepSize)) * m_SearchStepSize;
398 unsigned int numInputs = this->GetNumberOfIndexedInputs();
401 typename InputImageType::SizeType radius;
402 for (
unsigned int j = 0;j < InputImageType::ImageDimension;++j)
403 radius[j] = m_PatchHalfSize;
405 for (
unsigned int i = 0;i < numInputs;++i)
407 typename MeanAndVarianceImagesFilterType::Pointer filter = MeanAndVarianceImagesFilterType::New();
408 filter->SetInput(this->GetInput(i));
410 filter->SetRadius(radius);
411 filter->SetNumberOfWorkUnits(this->GetNumberOfWorkUnits());
415 meanImage->DisconnectPipeline();
417 varImage->DisconnectPipeline();
419 for (
unsigned int j = 0;j < numThreads;++j)
421 m_NLPatchSearchers[j].AddMeanImage(meanImage);
422 m_NLPatchSearchers[j].AddVarImage(varImage);
423 m_NLPatchSearchers[j].AddPatchTestImage(this->GetInput(i));
428 typename InputImageRegionType::IndexType baseIndex;
429 typename InputImageRegionType::IndexType valueIndex;
430 std::vector <double> noiseCovariances(numInputs,0);
432 for (
unsigned int i = 0;i < numInputs;++i)
434 typedef itk::ImageRegionConstIteratorWithIndex< InputImageType > InIteratorType;
435 typename InputImageType::RegionType largestRegion = this->GetInput(i)->GetLargestPossibleRegion();
436 InIteratorType dataIterator (this->GetInput(i), largestRegion);
438 double averageLocalSignal, diffSignal;
441 double averageCovariance = 0;
442 unsigned int numEstimations = 0;
443 unsigned int numLocalPixels = 2 * InputImageType::ImageDimension;
445 while (!dataIterator.IsAtEnd())
447 baseSignal =
static_cast<double>(dataIterator.Get());
448 baseIndex = dataIterator.GetIndex();
449 averageLocalSignal = 0;
451 for (
unsigned int d = 0; d < InputImageType::ImageDimension; ++d)
453 valueIndex = baseIndex;
454 int tmpIndex = baseIndex[d] - m_LocalNeighborhood;
455 valueIndex[d] = std::max(tmpIndex,0);
456 averageLocalSignal +=
static_cast<double> (this->GetInput(i)->GetPixel(valueIndex));
458 valueIndex = baseIndex;
459 tmpIndex = baseIndex[d] + m_LocalNeighborhood;
460 int maxIndex = largestRegion.GetSize()[d] - 1;
461 valueIndex[d] = std::min(tmpIndex, maxIndex);
462 averageLocalSignal +=
static_cast<double> (this->GetInput(i)->GetPixel(valueIndex));
465 averageLocalSignal /= numLocalPixels;
466 diffSignal = sqrt(numLocalPixels / (numLocalPixels + 1.0)) * (baseSignal - averageLocalSignal);
468 averageCovariance += diffSignal * diffSignal;
475 noiseCovariances[i] = averageCovariance / numEstimations;
479 for (
unsigned int i = 0;i < numThreads;++i)
481 m_NLPatchSearchers[i].SetPatchHalfSize(m_PatchHalfSize);
482 m_NLPatchSearchers[i].SetSearchStepSize(m_SearchStepSize);
483 m_NLPatchSearchers[i].SetMaxAbsDisp(maxAbsDisp);
484 m_NLPatchSearchers[i].SetInputImage(m_InitialT2Map);
485 m_NLPatchSearchers[i].SetBetaParameter(m_BetaParameter);
486 m_NLPatchSearchers[i].SetWeightThreshold(m_WeightThreshold);
487 m_NLPatchSearchers[i].SetMeanMinThreshold(m_MeanMinThreshold);
488 m_NLPatchSearchers[i].SetVarMinThreshold(m_VarMinThreshold);
489 m_NLPatchSearchers[i].SetNoiseCovariances(noiseCovariances);
493 template <
class TPixelScalarType>
498 std::vector <double> &workDataWeights, std::vector <OutputVectorType> &workDataSamples)
505 priorDistribution.fill(0);
507 double sumWeights = 0;
508 double maxWeight = -1;
509 for (
unsigned int i = 0;i < workDataSamples.size();++i)
512 for (
unsigned int j = 0;j < m_NumberOfT2Compartments;++j)
514 if (workDataSamples[i][j] != 0)
524 for (
unsigned int j = 0;j < m_NumberOfT2Compartments;++j)
525 priorDistribution[j] += workDataWeights[i] * workDataSamples[i][j];
527 if (maxWeight < workDataWeights[i])
528 maxWeight = workDataWeights[i];
530 sumWeights += workDataWeights[i];
536 for (
unsigned int j = 0;j < m_NumberOfT2Compartments;++j)
537 priorDistribution[j] += maxWeight * refDistribution[j];
539 sumWeights += maxWeight;
540 priorDistribution /= sumWeights;
543 template <
class TPixelScalarType>
550 unsigned int rowSize = AMatrix.rows() - AMatrix.cols();
551 unsigned int colSize = AMatrix.cols();
553 double lambda = std::sqrt(lambdaSq);
555 for (
unsigned int i = 0;i < colSize;++i)
557 signalValues[rowSize + i] = lambda * priorDistribution[i];
558 for (
unsigned int j = 0;j < colSize;++j)
561 AMatrix(rowSize + i,j) = 0;
563 AMatrix(rowSize + i,j) = lambda;
572 t2OptimizedWeights = nnlsOpt->GetCurrentPosition();
Applies an variance filter to an image.
OutputImageType::Pointer OutputImagePointer
virtual const std::vector< double > & GetDatabaseWeights() const
Implements an ITK wrapper for the NLOPT library.
void UpdateAtPosition(const IndexType &dataIndex)
void SetPoints(const ParametersType &points)
Cost function for estimating B1 from T2 relaxometry acquisition, following a multi-T2 EPG decay model...
VectorOutputImageType::PixelType OutputVectorType
Superclass::OutputImageRegionType OutputImageRegionType
void SetDataMatrix(const MatrixType &data)
std::vector< double > RealVectorType
Implements multi-peak T2 relaxometry estimation (with or without regularization)
NNLSOptimizerType::ParametersType T2VectorType
void StartOptimization() ITK_OVERRIDE
NNLSOptimizerType::MatrixType DataMatrixType
Non negative least squares optimizer. Implements Lawson et al method, of squared problem is activated...
void SetEchoSpacing(double val)
void SetExcitationFlipAngle(double val)
InputImageType::IndexType IndexType
NNLSOptimizerType::Pointer NNLSOptimizerPointer
virtual const std::vector< PixelType > & GetDatabaseSamples() const
RealVectorType & GetValue(double t1Value, double t2Value, double flipAngle, double m0Value)
Get EPG values at given point.
void SetNumberOfEchoes(unsigned int val)
double GetCurrentResidual()