ANIMA  4.0
animaRecursiveLineYvvGaussianImageFilter.h
Go to the documentation of this file.
1 #pragma once
2 
3 #include <itkInPlaceImageFilter.h>
4 #include <itkNumericTraits.h>
5 #include <itkVector.h>
6 #include <itkVariableLengthVector.h>
7 
8 namespace anima
9 {
10 template <typename TInputImage, typename TOutputImage=TInputImage>
12 public itk::InPlaceImageFilter<TInputImage,TOutputImage>
13 {
14 public:
17  typedef itk::InPlaceImageFilter<TInputImage,TOutputImage> Superclass;
18  typedef itk::SmartPointer<Self> Pointer;
19  typedef itk::SmartPointer<const Self> ConstPointer;
20 
22  itkNewMacro(Self)
23 
24 
25  itkTypeMacro(RecursiveLineYvvGaussianImageFilter, InPlaceImageFilter)
26 
27 
28  typedef typename TInputImage::Pointer InputImagePointer;
29  typedef typename TInputImage::ConstPointer InputImageConstPointer;
30 
36  typedef typename TInputImage::PixelType InputPixelType;
37  typedef typename itk::NumericTraits<InputPixelType>::RealType RealType;
38  typedef typename itk::NumericTraits<InputPixelType>::ScalarRealType ScalarRealType;
39 
40  typedef typename TOutputImage::RegionType OutputImageRegionType;
41 
43  typedef TInputImage InputImageType;
44 
46  typedef TOutputImage OutputImageType;
47 
49  itkGetConstMacro(Direction, unsigned int)
50 
52  itkSetMacro(Direction, unsigned int)
53 
55  void SetInputImage( const TInputImage * );
56 
58  const TInputImage * GetInputImage( void );
59 
75  itkSetMacro(NormalizeAcrossScale, bool)
76  itkGetConstMacro(NormalizeAcrossScale, bool)
77 
80  itkGetConstMacro(Sigma, ScalarRealType)
81  itkSetMacro(Sigma, ScalarRealType)
82 
83 protected:
85  virtual ~RecursiveLineYvvGaussianImageFilter() {}
86  void PrintSelf(std::ostream& os, itk::Indent indent) const ITK_OVERRIDE;
87 
89  void BeforeThreadedGenerateData() ITK_OVERRIDE;
90  void GenerateData() ITK_OVERRIDE;
91  void DynamicThreadedGenerateData(const OutputImageRegionType& outputRegionForThread) ITK_OVERRIDE;
92 
101  void EnlargeOutputRequestedRegion(itk::DataObject *output) ITK_OVERRIDE;
102 
107  virtual void SetUp(ScalarRealType spacing);
108 
115  void FilterDataArray(RealType *outs, const RealType *data, unsigned int ln,
116  RealType &sV0, RealType &sV1, RealType &sV2);
117 
119  template <class T>
120  inline void ComputeCausalPart(T &out, const T &data, T &V0, T &V1, T &V2)
121  {
122  out = data + V0 * m_B1 + V1 * m_B2 + V2 * m_B3;
123  V2 = V1;V1 = V0;V0 = out;
124  }
125 
127  template <class T>
128  inline void ComputeCausalPart(itk::VariableLengthVector<T> &out, const itk::VariableLengthVector<T> &data,
129  itk::VariableLengthVector<T> &V0, itk::VariableLengthVector<T> &V1,
130  itk::VariableLengthVector<T> &V2)
131  {
132  unsigned int vSize = data.GetSize();
133  if (out.GetSize() != vSize)
134  out.SetSize(vSize);
135 
136  for (unsigned int i = 0;i < vSize;++i)
137  {
138  out[i] = data[i] + V0[i] * m_B1 + V1[i] * m_B2 + V2[i] * m_B3;
139  V2[i] = V1[i];V1[i] = V0[i];V0[i] = out[i];
140  }
141  }
142 
144  template <class T>
145  inline void ComputeAntiCausalPart(T &out, T &data, T &V0, T &V1, T &V2)
146  {
147  out = data * m_B + V0 * m_B1 + V1 * m_B2 + V2 * m_B3;
148  V2 = V1;V1 = V0;V0 = out;
149  }
150 
152  template <class T>
153  inline void ComputeAntiCausalPart(itk::VariableLengthVector<T> &out, itk::VariableLengthVector<T> &data,
154  itk::VariableLengthVector<T> &V0, itk::VariableLengthVector<T> &V1,
155  itk::VariableLengthVector<T> &V2)
156  {
157  unsigned int vSize = data.GetSize();
158  for (unsigned int i = 0;i < vSize;++i)
159  {
160  out[i] = data[i] * m_B + V0[i] * m_B1 + V1[i] * m_B2 + V2[i] * m_B3;
161  V2[i] = V1[i];V1[i] = V0[i];V0[i] = out[i];
162  }
163  }
164 
165  template <class T>
166  inline void ComputeCausalBase(const T &data, T &V0, T &V1, T &V2)
167  {
168  V0 = data / (1.0 - m_B1 - m_B2 - m_B3);V1 = V0;V2 = V1;
169  }
170 
171  template <class T>
172  inline void ComputeCausalBase(const itk::VariableLengthVector<T> &data, itk::VariableLengthVector<T> &V0,
173  itk::VariableLengthVector<T> &V1, itk::VariableLengthVector<T> &V2)
174  {
175  unsigned int vSize = data.GetSize();
176  if (V0.GetSize() != vSize)
177  V0.SetSize(vSize);
178  if (V1.GetSize() != vSize)
179  V1.SetSize(vSize);
180  if (V2.GetSize() != vSize)
181  V2.SetSize(vSize);
182 
183  for (unsigned int i = 0;i < vSize;++i)
184  {
185  V0[i] = data[i] / (1.0 - m_B1 - m_B2 - m_B3);V1[i] = V0[i];V2[i] = V1[i];
186  }
187  }
188 
189  template <class T>
190  inline void ComputeAntiCausalBase(const T &data, T *outs, T &V0, T &V1, T &V2, unsigned int ln)
191  {
192  // Handle outside values according to Triggs and Sdika
193  const T u_p = data / (1.0 - m_B1 - m_B2 - m_B3);
194  const T v_p = u_p / (1.0 - m_B1 - m_B2 - m_B3);
195 
196  V0 = v_p;
197  V1 = v_p;
198  V2 = v_p;
199 
200  for (unsigned int i = 0;i < 3;++i)
201  {
202  V0 += (outs[ln - 1 - i] - u_p) * m_MMatrix(0,i);
203  V1 += (outs[ln - 1 - i] - u_p) * m_MMatrix(1,i);
204  V2 += (outs[ln - 1 - i] - u_p) * m_MMatrix(2,i);
205  }
206 
207  // This was not in the 2006 Triggs paper but sounds quite logical since m_B is not one
208  V0 *= m_B;
209  V1 *= m_B;
210  V2 *= m_B;
211  }
212 
213  template <class T>
214  inline void ComputeAntiCausalBase(const itk::VariableLengthVector<T> &data, itk::VariableLengthVector<T> *outs,
215  itk::VariableLengthVector<T> &V0, itk::VariableLengthVector<T> &V1,
216  itk::VariableLengthVector<T> &V2, unsigned int ln)
217  {
218  unsigned int vSize = data.GetSize();
219  // Handle outside values according to Triggs and Sdika
220  double factor = 1.0 / (1.0 - m_B1 - m_B2 - m_B3);
221 
222  for (unsigned int i = 0;i < vSize;++i)
223  {
224  V0[i] = data[i] * factor * factor;
225  V1[i] = data[i] * factor * factor;
226  V2[i] = data[i] * factor * factor;
227  }
228 
229  for (unsigned int i = 0;i < 3;++i)
230  for (unsigned int j = 0;j < vSize;++j)
231  {
232  V0[j] += (outs[ln - 1 - i][j] - data[j] * factor) * m_MMatrix(0,i);
233  V1[j] += (outs[ln - 1 - i][j] - data[j] * factor) * m_MMatrix(1,i);
234  V2[j] += (outs[ln - 1 - i][j] - data[j] * factor) * m_MMatrix(2,i);
235  }
236 
237  // This was not in the 2006 Triggs paper but sounds quite logical since m_B is not one
238  for (unsigned int i = 0;i < vSize;++i)
239  {
240  V0[i] *= m_B;
241  V1[i] *= m_B;
242  V2[i] *= m_B;
243  }
244  }
245 
251 
252  // Initialization matrix for anti-causal pass
253  vnl_matrix <ScalarRealType> m_MMatrix;
254 
255 private:
256  ITK_DISALLOW_COPY_AND_ASSIGN(RecursiveLineYvvGaussianImageFilter);
257 
260  unsigned int m_Direction;
261 
263  ScalarRealType m_Sigma;
264 
266  bool m_NormalizeAcrossScale;
267 };
268 
269 
270 } // end of namespace anima
271 
void EnlargeOutputRequestedRegion(itk::DataObject *output) ITK_OVERRIDE
itk::InPlaceImageFilter< TInputImage, TOutputImage > Superclass
void ComputeAntiCausalPart(itk::VariableLengthVector< T > &out, itk::VariableLengthVector< T > &data, itk::VariableLengthVector< T > &V0, itk::VariableLengthVector< T > &V1, itk::VariableLengthVector< T > &V2)
itk::NumericTraits< InputPixelType >::ScalarRealType ScalarRealType
void ComputeCausalPart(itk::VariableLengthVector< T > &out, const itk::VariableLengthVector< T > &data, itk::VariableLengthVector< T > &V0, itk::VariableLengthVector< T > &V1, itk::VariableLengthVector< T > &V2)
void ComputeAntiCausalPart(T &out, T &data, T &V0, T &V1, T &V2)
itk::NumericTraits< InputPixelType >::RealType RealType
void ComputeAntiCausalBase(const T &data, T *outs, T &V0, T &V1, T &V2, unsigned int ln)
void FilterDataArray(RealType *outs, const RealType *data, unsigned int ln, RealType &sV0, RealType &sV1, RealType &sV2)
void PrintSelf(std::ostream &os, itk::Indent indent) const ITK_OVERRIDE
void ComputeCausalPart(T &out, const T &data, T &V0, T &V1, T &V2)
void ComputeAntiCausalBase(const itk::VariableLengthVector< T > &data, itk::VariableLengthVector< T > *outs, itk::VariableLengthVector< T > &V0, itk::VariableLengthVector< T > &V1, itk::VariableLengthVector< T > &V2, unsigned int ln)
void DynamicThreadedGenerateData(const OutputImageRegionType &outputRegionForThread) ITK_OVERRIDE
void ComputeCausalBase(const itk::VariableLengthVector< T > &data, itk::VariableLengthVector< T > &V0, itk::VariableLengthVector< T > &V1, itk::VariableLengthVector< T > &V2)