ANIMA  4.0
animaDistortionCorrectionImageFilter.hxx
Go to the documentation of this file.
1 #pragma once
3 
4 #include <itkVariableLengthVector.h>
5 #include <itkImageRegionIterator.h>
6 #include <itkImageRegionConstIterator.h>
7 #include <itkImageLinearIteratorWithIndex.h>
8 
11 
12 #include <vector>
13 #include <cstddef>
14 #include <iterator>
15 #include <iostream>
16 
17 namespace anima
18 {
19 
20 template< typename TInputImage >
21 DistortionCorrectionImageFilter < TInputImage >
23 {
24  m_Direction = 0;
25  this->SetNumberOfRequiredInputs( 2 );
26 }
27 
28 template< typename TInputImage >
30 ::GenerateOutputInformation(void)
31 {
32  // Override the method in itkImageSource, so we can set the vector length of
33  // the output itk::VectorImage
34 
35  this->Superclass::GenerateOutputInformation();
36 
37  OutputImageType *output = this->GetOutput();
38  output->SetVectorLength( 3 );
39 }
40 
41 template< typename TInputImage >
43 ::GenerateData()
44 {
45  // Call a method that can be overridden by a subclass to allocate
46  // memory for the filter's outputs
47  this->AllocateOutputs();
48 
49  // Call a method that can be overridden by a subclass to perform
50  // some calculations prior to splitting the main computations into
51  // separate threads
52  this->BeforeThreadedGenerateData();
53 
54  using RegionType = itk::ImageRegion <TInputImage::ImageDimension>;
55  typename OutputImageType::Pointer outputImage(this->GetOutput());
56  const RegionType region = outputImage->GetRequestedRegion();
57 
58  this->GetMultiThreader()->SetNumberOfWorkUnits(this->GetNumberOfWorkUnits());
59  this->GetMultiThreader()->template ParallelizeImageRegionRestrictDirection<OutputImageType::ImageDimension>(
60  this->m_Direction, region, [this](const RegionType & lambdaRegion) { this->DynamicThreadedGenerateData(lambdaRegion); }, this);
61 }
62 
63 template< typename TInputImage >
65 ::BeforeThreadedGenerateData()
66 {
67  Superclass::BeforeThreadedGenerateData();
68 
69  m_ReferenceGeometry = this->GetInput(0)->GetDirection();
70 
71  for (unsigned int i = 0;i < 3;++i)
72  for (unsigned int j = 0;j < 3;++j)
73  m_ReferenceGeometry(i,j) *= this->GetInput(0)->GetSpacing()[j];
74 }
75 
76 template< typename TInputImage >
78 ::DynamicThreadedGenerateData(const OutputImageRegionType& outputRegionForThread)
79 {
80  InputImagePointer forwardImage = this->GetInput(0);
81  InputImagePointer backwardImage = this->GetInput(1);
82 
83  typedef itk::ImageLinearConstIteratorWithIndex< TInputImage > ConstIteratorType;
84  typedef itk::ImageLinearIteratorWithIndex< OutputImageType > IteratorType;
85 
86  ConstIteratorType forwardIt (forwardImage, outputRegionForThread);
87  ConstIteratorType backwardIt (backwardImage, outputRegionForThread);
88  IteratorType outputIt (this->GetOutput(), outputRegionForThread);
89 
90  forwardIt.SetDirection(m_Direction);
91  backwardIt.SetDirection(m_Direction);
92  outputIt.SetDirection(m_Direction);
93 
94  unsigned int lengthLine = outputRegionForThread.GetSize()[ this->m_Direction ];
95 
96  double step = 0;
97  double epsilon = 0.0001;
98  itk::VariableLengthVector< PixelType > fillVectorField(3);
99  itk::VariableLengthVector< PixelType > fillScaledVectorField(3);
100  fillVectorField.Fill(0);
101  fillScaledVectorField.Fill(0);
102 
103  std::vector<double> forwardScale;
104  std::vector<double> backwardScale;
105  std::vector<double> middleScale;
106  std::vector<double> differenceForwardBackward;
107  std::vector<double> differenceForwardBackwardResampled;
108 
109  while (!forwardIt.IsAtEnd())
110  {
111  // Compute Cumulative sum
112  std::vector<double> forwardCumulatedLine;
113  std::vector<double> backwardCumulatedLine;
114 
115  forwardCumulatedLine.push_back(forwardIt.Get());
116  backwardCumulatedLine.push_back(backwardIt.Get());
117 
118  ++forwardIt;
119  ++backwardIt;
120 
121  while (!forwardIt.IsAtEndOfLine())
122  {
123 
124  forwardCumulatedLine.push_back(forwardCumulatedLine.back() + forwardIt.Get() + epsilon);
125  backwardCumulatedLine.push_back(backwardCumulatedLine.back() + backwardIt.Get() + epsilon);
126 
127  ++forwardIt;
128  ++backwardIt;
129 
130  }
131 
132  double forwMax = forwardCumulatedLine.back();
133  double forwMin = forwardCumulatedLine.front();
134  double backMax = backwardCumulatedLine.back();
135  double backMin = backwardCumulatedLine.front();
136 
137  for (int i=0; i<lengthLine; i++)
138  backwardCumulatedLine[i]-= backMin;
139  for (int i=0; i<lengthLine; i++)
140  backwardCumulatedLine[i]*= (forwMax - forwMin)/(backMax - backMin);
141  for (int i=0; i<lengthLine; i++)
142  backwardCumulatedLine[i]+= forwMin;
143 
144  step = (forwMax - forwMin)/(10*forwardCumulatedLine.size());
145 
146  InverseCubicInterpolator(forwardCumulatedLine, forwardScale, step);
147  InverseCubicInterpolator(backwardCumulatedLine, backwardScale, step);
148 
149  differenceForwardBackward.resize(forwardScale.size());
150  middleScale.resize(forwardScale.size());
151 
152  for (int i=0; i< differenceForwardBackward.size(); i++)
153  {
154  differenceForwardBackward[i] = (forwardScale[i] - backwardScale[i])/2;
155  middleScale[i] = (forwardScale[i] + backwardScale[i])/2;
156  }
157 
158  CubicInterpolator<double>(differenceForwardBackward, middleScale, differenceForwardBackwardResampled, lengthLine);
159 
160  unsigned int pos = 0;
161 
162  while (!outputIt.IsAtEndOfLine())
163  {
164  fillVectorField[m_Direction] = differenceForwardBackwardResampled[pos];
165 
166  for (unsigned int i = 0;i < 3;++i)
167  {
168  fillScaledVectorField[i] = m_ReferenceGeometry(i,m_Direction) * fillVectorField[m_Direction];
169  }
170 
171  outputIt.Set(fillScaledVectorField);
172  ++outputIt;
173  pos++;
174  }
175 
176  forwardIt.NextLine();
177  backwardIt.NextLine();
178  outputIt.NextLine();
179  }
180 }
181 
182 template< typename TInputImage >
184 ::AfterThreadedGenerateData()
185 {
186  // Final step to smooth the obtained vector field
187  if (m_FieldSmoothingSigma > 0)
188  {
190 
191  typename SmoothingFilterType::Pointer fieldSmoother = SmoothingFilterType::New();
192  fieldSmoother->SetInput(this->GetOutput());
193  fieldSmoother->SetNumberOfWorkUnits(this->GetNumberOfWorkUnits());
194  fieldSmoother->SetSigma(m_FieldSmoothingSigma);
195 
196  fieldSmoother->Update();
197 
198  this->SetNthOutput(0,fieldSmoother->GetOutput());
199  }
200 
201  Superclass::AfterThreadedGenerateData();
202 }
203 
204 } // end namespace anima
itk::VectorImage< typename TInputImage::InternalPixelType, TInputImage::ImageDimension > OutputImageType
void InverseCubicInterpolator(std::vector< T > &inputVect, std::vector< T > &outputVect, T step)