OTB  9.0.0
Orfeo Toolbox
otbBSplineDecompositionImageFilter.hxx
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2005-2022 Centre National d'Etudes Spatiales (CNES)
3  *
4  * This file is part of Orfeo Toolbox
5  *
6  * https://www.orfeo-toolbox.org/
7  *
8  * Licensed under the Apache License, Version 2.0 (the "License");
9  * you may not use this file except in compliance with the License.
10  * You may obtain a copy of the License at
11  *
12  * http://www.apache.org/licenses/LICENSE-2.0
13  *
14  * Unless required by applicable law or agreed to in writing, software
15  * distributed under the License is distributed on an "AS IS" BASIS,
16  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
17  * See the License for the specific language governing permissions and
18  * limitations under the License.
19  */
20 
21 #ifndef otbBSplineDecompositionImageFilter_hxx
22 #define otbBSplineDecompositionImageFilter_hxx
24 #include "itkImageRegionConstIteratorWithIndex.h"
25 #include "itkImageRegionIterator.h"
26 #include "itkProgressReporter.h"
27 #include "itkVector.h"
28 
29 namespace otb
30 {
31 
35 template <class TInputImage, class TOutputImage>
37 {
38  m_SplineOrder = 0;
39  int SplineOrder = 3;
40  m_Tolerance = 1e-10; // Need some guidance on this one...what is reasonable?
41  m_IteratorDirection = 0;
42  this->SetSplineOrder(SplineOrder);
43 }
45 
49 template <class TInputImage, class TOutputImage>
50 void BSplineDecompositionImageFilter<TInputImage, TOutputImage>::PrintSelf(std::ostream& os, itk::Indent indent) const
51 {
52  Superclass::PrintSelf(os, indent);
53  os << indent << "Spline Order: " << m_SplineOrder << std::endl;
54 }
56 
57 template <class TInputImage, class TOutputImage>
59 {
60 
61  // See Unser, 1993, Part II, Equation 2.5,
62  // or Unser, 1999, Box 2. for an explanation.
63 
64  double c0 = 1.0;
65 
66  if (m_DataLength[m_IteratorDirection] == 1) // Required by mirror boundaries
67  {
68  return false;
69  }
70 
71  // Compute overall gain
72  for (int k = 0; k < m_NumberOfPoles; ++k)
73  {
74  // Note for cubic splines lambda = 6
75  c0 = c0 * (1.0 - m_SplinePoles[k]) * (1.0 - 1.0 / m_SplinePoles[k]);
76  }
77 
78  // apply the gain
79  for (unsigned int n = 0; n < m_DataLength[m_IteratorDirection]; ++n)
80  {
81  m_Scratch[n] *= c0;
82  }
83 
84  // loop over all poles
85  for (int k = 0; k < m_NumberOfPoles; ++k)
86  {
87  // causal initialization
88  this->SetInitialCausalCoefficient(m_SplinePoles[k]);
89  // causal recursion
90  for (unsigned int n = 1; n < m_DataLength[m_IteratorDirection]; ++n)
91  {
92  m_Scratch[n] += m_SplinePoles[k] * m_Scratch[n - 1];
93  }
94 
95  // anticausal initialization
96  this->SetInitialAntiCausalCoefficient(m_SplinePoles[k]);
97  // anticausal recursion
98  for (int n = m_DataLength[m_IteratorDirection] - 2; 0 <= n; n--)
99  {
100  m_Scratch[n] = m_SplinePoles[k] * (m_Scratch[n + 1] - m_Scratch[n]);
101  }
102  }
103  return true;
104 }
105 
106 template <class TInputImage, class TOutputImage>
108 {
109  if (SplineOrder == m_SplineOrder)
110  {
111  return;
112  }
113  m_SplineOrder = SplineOrder;
114  this->SetPoles();
115  this->Modified();
116 }
117 
118 template <class TInputImage, class TOutputImage>
120 {
121  /* See Unser, 1997. Part II, Table I for Pole values */
122  // See also, Handbook of Medical Imaging, Processing and Analysis, Ed. Isaac N. Bankman,
123  // 2000, pg. 416.
124  switch (m_SplineOrder)
125  {
126  case 3:
127  m_NumberOfPoles = 1;
128  m_SplinePoles[0] = std::sqrt(3.0) - 2.0;
129  break;
130  case 0:
131  m_NumberOfPoles = 0;
132  break;
133  case 1:
134  m_NumberOfPoles = 0;
135  break;
136  case 2:
137  m_NumberOfPoles = 1;
138  m_SplinePoles[0] = std::sqrt(8.0) - 3.0;
139  break;
140  case 4:
141  m_NumberOfPoles = 2;
142  m_SplinePoles[0] = std::sqrt(664.0 - std::sqrt(438976.0)) + std::sqrt(304.0) - 19.0;
143  m_SplinePoles[1] = std::sqrt(664.0 + std::sqrt(438976.0)) - std::sqrt(304.0) - 19.0;
144  break;
145  case 5:
146  m_NumberOfPoles = 2;
147  m_SplinePoles[0] = std::sqrt(135.0 / 2.0 - std::sqrt(17745.0 / 4.0)) + std::sqrt(105.0 / 4.0) - 13.0 / 2.0;
148  m_SplinePoles[1] = std::sqrt(135.0 / 2.0 + std::sqrt(17745.0 / 4.0)) - std::sqrt(105.0 / 4.0) - 13.0 / 2.0;
149  break;
150  default:
151  // SplineOrder not implemented yet.
152  itk::ExceptionObject err(__FILE__, __LINE__);
153  err.SetLocation(ITK_LOCATION);
154  err.SetDescription("SplineOrder must be between 0 and 5. Requested spline order has not been implemented yet.");
155  throw err;
156  break;
157  }
158 }
159 
160 template <class TInputImage, class TOutputImage>
162 {
163  /* beginning InitialCausalCoefficient */
164  /* See Unser, 1999, Box 2 for explanation */
165 
166  double sum, zn, z2n, iz;
167  unsigned long horizon;
168 
169  /* this initialization corresponds to mirror boundaries */
170  horizon = m_DataLength[m_IteratorDirection];
171  zn = z;
172  if (m_Tolerance > 0.0)
173  {
174  horizon = (long)std::ceil(log(m_Tolerance) / std::log(fabs(z)));
175  }
176  if (horizon < m_DataLength[m_IteratorDirection])
177  {
178  /* accelerated loop */
179  sum = m_Scratch[0]; // verify this
180  for (unsigned int n = 1; n < horizon; ++n)
181  {
182  sum += zn * m_Scratch[n];
183  zn *= z;
184  }
185  m_Scratch[0] = sum;
186  }
187  else
188  {
189  /* full loop */
190  iz = 1.0 / z;
191  z2n = std::pow(z, (double)(m_DataLength[m_IteratorDirection] - 1L));
192  sum = m_Scratch[0] + z2n * m_Scratch[m_DataLength[m_IteratorDirection] - 1L];
193  z2n *= z2n * iz;
194  for (unsigned int n = 1; n <= (m_DataLength[m_IteratorDirection] - 2); ++n)
195  {
196  sum += (zn + z2n) * m_Scratch[n];
197  zn *= z;
198  z2n *= iz;
199  }
200  m_Scratch[0] = sum / (1.0 - zn * zn);
201  }
202 }
203 
204 template <class TInputImage, class TOutputImage>
206 {
207  // this initialization corresponds to mirror boundaries
208  /* See Unser, 1999, Box 2 for explanation */
209  // Also see erratum at http://bigwww.epfl.ch/publications/unser9902.html
210  m_Scratch[m_DataLength[m_IteratorDirection] - 1] =
211  (z / (z * z - 1.0)) * (z * m_Scratch[m_DataLength[m_IteratorDirection] - 2] + m_Scratch[m_DataLength[m_IteratorDirection] - 1]);
212 }
213 
214 template <class TInputImage, class TOutputImage>
216 {
217  OutputImagePointer output = this->GetOutput();
218 
219  itk::Size<ImageDimension> size = output->GetBufferedRegion().GetSize();
220 
221  unsigned int count = output->GetBufferedRegion().GetNumberOfPixels() / size[0] * ImageDimension;
222 
223  itk::ProgressReporter progress(this, 0, count, 10);
224 
225  // Initialize coefficient array
226  this->CopyImageToImage(); // Coefficients are initialized to the input data
227 
228  for (unsigned int n = 0; n < ImageDimension; ++n)
229  {
230  m_IteratorDirection = n;
231  // Loop through each dimension
232 
233  // Initialize iterators
234  OutputLinearIterator CIterator(output, output->GetBufferedRegion());
235  CIterator.SetDirection(m_IteratorDirection);
236  // For each data vector
237  while (!CIterator.IsAtEnd())
238  {
239  // Copy coefficients to scratch
240  this->CopyCoefficientsToScratch(CIterator);
241 
242  // Perform 1D BSpline calculations
243  this->DataToCoefficients1D();
244 
245  // Copy scratch back to coefficients.
246  // Brings us back to the end of the line we were working on.
247  CIterator.GoToBeginOfLine();
248  this->CopyScratchToCoefficients(CIterator); // m_Scratch = m_Image;
249  CIterator.NextLine();
250  progress.CompletedPixel();
251  }
252  }
253 }
254 
258 template <class TInputImage, class TOutputImage>
260 {
261 
262  typedef itk::ImageRegionConstIteratorWithIndex<TInputImage> InputIterator;
263  typedef itk::ImageRegionIterator<TOutputImage> OutputIterator;
264  typedef typename TOutputImage::PixelType OutputPixelType;
265 
266  InputIterator inIt(this->GetInput(), this->GetInput()->GetBufferedRegion());
267  OutputIterator outIt(this->GetOutput(), this->GetOutput()->GetBufferedRegion());
268 
269  inIt.GoToBegin();
270  outIt.GoToBegin();
271 
272  while (!outIt.IsAtEnd())
273  {
274  outIt.Set(static_cast<OutputPixelType>(inIt.Get()));
275  ++inIt;
276  ++outIt;
277  }
278 }
279 
283 template <class TInputImage, class TOutputImage>
285 {
286  typedef typename TOutputImage::PixelType OutputPixelType;
287  unsigned long j = 0;
288  while (!Iter.IsAtEndOfLine())
289  {
290  Iter.Set(static_cast<OutputPixelType>(m_Scratch[j]));
291  ++Iter;
292  ++j;
293  }
294 }
296 
300 template <class TInputImage, class TOutputImage>
302 {
303  unsigned long j = 0;
304  while (!Iter.IsAtEndOfLine())
305  {
306  m_Scratch[j] = static_cast<double>(Iter.Get());
307  ++Iter;
308  ++j;
309  }
310 }
312 
316 template <class TInputImage, class TOutputImage>
318 {
319 
320  // Allocate scratch memory
321  InputImageConstPointer inputPtr = this->GetInput();
322  m_DataLength = inputPtr->GetBufferedRegion().GetSize();
323 
324  unsigned long maxLength = 0;
325  for (unsigned int n = 0; n < ImageDimension; ++n)
326  {
327  if (m_DataLength[n] > maxLength)
328  {
329  maxLength = m_DataLength[n];
330  }
331  }
332  m_Scratch.resize(maxLength);
333 
334  // Allocate memory for output image
335  OutputImagePointer outputPtr = this->GetOutput();
336  outputPtr->SetBufferedRegion(outputPtr->GetRequestedRegion());
337  outputPtr->Allocate();
338 
339  // Calculate actual output
340  this->DataToCoefficientsND();
341 
342  // Clean up
343  m_Scratch.clear();
344 }
345 
346 } // namespace otb
347 
348 #endif
otb::BSplineDecompositionImageFilter::PrintSelf
void PrintSelf(std::ostream &os, itk::Indent indent) const override
Definition: otbBSplineDecompositionImageFilter.hxx:50
otb::BSplineDecompositionImageFilter::SetSplineOrder
void SetSplineOrder(unsigned int SplineOrder)
Definition: otbBSplineDecompositionImageFilter.hxx:107
otb::BSplineDecompositionImageFilter::OutputLinearIterator
itk::ImageLinearIteratorWithIndex< TOutputImage > OutputLinearIterator
Definition: otbBSplineDecompositionImageFilter.h:70
otb::BSplineDecompositionImageFilter::CopyScratchToCoefficients
void CopyScratchToCoefficients(OutputLinearIterator &)
Definition: otbBSplineDecompositionImageFilter.hxx:284
otb
The "otb" namespace contains all Orfeo Toolbox (OTB) classes.
Definition: otbJoinContainer.h:32
otb::BSplineDecompositionImageFilter::DataToCoefficientsND
void DataToCoefficientsND()
Definition: otbBSplineDecompositionImageFilter.hxx:215
otbBSplineDecompositionImageFilter.h
otb::BSplineDecompositionImageFilter::CopyImageToImage
void CopyImageToImage()
Definition: otbBSplineDecompositionImageFilter.hxx:259
otb::BSplineDecompositionImageFilter::InputImageConstPointer
Superclass::InputImageConstPointer InputImageConstPointer
Definition: otbBSplineDecompositionImageFilter.h:61
otb::BSplineDecompositionImageFilter::SetInitialCausalCoefficient
virtual void SetInitialCausalCoefficient(double z)
Definition: otbBSplineDecompositionImageFilter.hxx:161
otb::BSplineDecompositionImageFilter::SetPoles
virtual void SetPoles()
Definition: otbBSplineDecompositionImageFilter.hxx:119
otb::BSplineDecompositionImageFilter::DataToCoefficients1D
virtual bool DataToCoefficients1D()
Definition: otbBSplineDecompositionImageFilter.hxx:58
otb::BSplineDecompositionImageFilter::CopyCoefficientsToScratch
void CopyCoefficientsToScratch(OutputLinearIterator &)
Definition: otbBSplineDecompositionImageFilter.hxx:301
otb::BSplineDecompositionImageFilter::BSplineDecompositionImageFilter
BSplineDecompositionImageFilter()
Definition: otbBSplineDecompositionImageFilter.hxx:36
otb::BSplineDecompositionImageFilter::GenerateData
void GenerateData() override
Definition: otbBSplineDecompositionImageFilter.hxx:317
otb::BSplineDecompositionImageFilter::SetInitialAntiCausalCoefficient
virtual void SetInitialAntiCausalCoefficient(double z)
Definition: otbBSplineDecompositionImageFilter.hxx:205
otb::BSplineDecompositionImageFilter::OutputImagePointer
Superclass::OutputImagePointer OutputImagePointer
Definition: otbBSplineDecompositionImageFilter.h:62