4 #include <itkImageRegionIterator.h> 5 #include <itkImageRegionConstIterator.h> 9 template <
typename TInputImage,
typename TOutputImage>
11 T1RelaxometryEstimationImageFilter <TInputImage,TOutputImage>
12 ::BeforeThreadedGenerateData()
14 if (this->GetNumberOfIndexedInputs() != m_FlipAngles.size())
15 itkExceptionMacro(
"There should be the same number of inputs and flip angles");
17 if (this->GetNumberOfIndexedInputs() != 2)
18 itkExceptionMacro(
"There should be only two inputs for DESPOT1");
20 Superclass::BeforeThreadedGenerateData();
23 template <
typename TInputImage,
typename TOutputImage>
26 ::CheckComputationMask()
28 if (this->GetComputationMask())
31 typedef itk::ImageRegionConstIterator <InputImageType> IteratorType;
32 typedef itk::ImageRegionIterator <MaskImageType> MaskIteratorType;
34 std::vector <IteratorType> inItrs(this->GetNumberOfIndexedInputs());
35 for (
unsigned int i = 0;i < this->GetNumberOfIndexedInputs();++i)
36 inItrs[i] = IteratorType(this->GetInput(i),this->GetOutput()->GetLargestPossibleRegion());
38 typename MaskImageType::Pointer maskImage = MaskImageType::New();
39 maskImage->Initialize();
40 maskImage->SetRegions(this->GetInput(0)->GetLargestPossibleRegion());
41 maskImage->SetSpacing (this->GetInput(0)->GetSpacing());
42 maskImage->SetOrigin (this->GetInput(0)->GetOrigin());
43 maskImage->SetDirection (this->GetInput(0)->GetDirection());
44 maskImage->Allocate();
46 MaskIteratorType maskItr (maskImage,this->GetOutput()->GetLargestPossibleRegion());
47 while (!maskItr.IsAtEnd())
49 double averageVal = 0;
50 for (
unsigned int i = 0;i < this->GetNumberOfIndexedInputs();++i)
51 averageVal += inItrs[i].Get();
53 averageVal /= this->GetNumberOfIndexedInputs();
55 bool maskPoint = (averageVal <= m_AverageSignalThreshold);
62 for (
unsigned int i = 0;i < this->GetNumberOfIndexedInputs();++i)
68 this->SetComputationMask(maskImage);
71 template <
typename TInputImage,
typename TOutputImage>
76 typedef itk::ImageRegionConstIterator <InputImageType> ImageIteratorType;
77 typedef itk::ImageRegionIterator <OutputImageType> OutImageIteratorType;
79 OutImageIteratorType outT1Iterator(this->GetOutput(0),outputRegionForThread);
80 OutImageIteratorType outM0Iterator(this->GetOutput(1),outputRegionForThread);
82 ImageIteratorType firstInputIterator(this->GetInput(0),outputRegionForThread);
83 ImageIteratorType secondInputIterator(this->GetInput(1),outputRegionForThread);
85 OutImageIteratorType b1MapItr;
87 b1MapItr = OutImageIteratorType (m_B1Map,outputRegionForThread);
89 typedef itk::ImageRegionIterator <MaskImageType> MaskIteratorType;
90 MaskIteratorType maskItr(this->GetComputationMask(),outputRegionForThread);
92 while (!maskItr.IsAtEnd())
96 b1Value = b1MapItr.Get();
98 if ((maskItr.Get() == 0)||(b1Value == 0))
100 outT1Iterator.Set(0);
101 outM0Iterator.Set(0);
107 ++firstInputIterator;
108 ++secondInputIterator;
116 double basicT1Value = m_TRValue;
118 double lnValue = secondInputIterator.Get() * std::sin(m_FlipAngles[0] * b1Value) * std::cos(m_FlipAngles[1] * b1Value) - firstInputIterator.Get() * std::cos(m_FlipAngles[0] * b1Value) * std::sin(m_FlipAngles[1] * b1Value);
119 lnValue /= secondInputIterator.Get() * std::sin(m_FlipAngles[0] * b1Value) - firstInputIterator.Get() * std::sin(m_FlipAngles[1] * b1Value);
121 if (lnValue != lnValue)
123 outT1Iterator.Set(1e-4);
124 outM0Iterator.Set(0);
130 ++firstInputIterator;
131 ++secondInputIterator;
139 lnValue = std::log(std::abs(lnValue));
141 basicT1Value /= lnValue;
143 if ((basicT1Value < 0)||(basicT1Value > m_T1UpperBoundValue))
144 basicT1Value = m_T1UpperBoundValue;
146 outT1Iterator.Set(basicT1Value);
149 double sumSignalData = 0;
150 double sumSquaredData = 0;
152 double dataValue = (1.0 - std::exp(- m_TRValue / basicT1Value)) * std::sin(m_FlipAngles[0] * b1Value) / (1.0 - std::exp(-m_TRValue / basicT1Value) * std::cos(m_FlipAngles[0] * b1Value));
154 sumSignalData += firstInputIterator.Get() * dataValue;
155 sumSquaredData += dataValue * dataValue;
157 dataValue = (1.0 - std::exp(- m_TRValue / basicT1Value)) * std::sin(m_FlipAngles[1] * b1Value) / (1.0 - std::exp(- m_TRValue / basicT1Value) * std::cos(m_FlipAngles[1] * b1Value));
159 sumSignalData += secondInputIterator.Get() * dataValue;
160 sumSquaredData += dataValue * dataValue;
162 double m0Value = sumSignalData / sumSquaredData;
164 if (m0Value <= 1.0e-4)
166 else if (m0Value > m_M0UpperBoundValue)
167 m0Value = m_M0UpperBoundValue;
169 outM0Iterator.Set(m0Value);
175 ++firstInputIterator;
176 ++secondInputIterator;
Superclass::OutputImageRegionType OutputImageRegionType