21 #ifndef otbBSplineDecompositionImageFilter_hxx
22 #define otbBSplineDecompositionImageFilter_hxx
24 #include "itkImageRegionConstIteratorWithIndex.h"
25 #include "itkImageRegionIterator.h"
26 #include "itkProgressReporter.h"
27 #include "itkVector.h"
35 template <
class TInputImage,
class TOutputImage>
41 m_IteratorDirection = 0;
42 this->SetSplineOrder(SplineOrder);
49 template <
class TInputImage,
class TOutputImage>
52 Superclass::PrintSelf(os, indent);
53 os << indent <<
"Spline Order: " << m_SplineOrder << std::endl;
57 template <
class TInputImage,
class TOutputImage>
66 if (m_DataLength[m_IteratorDirection] == 1)
72 for (
int k = 0; k < m_NumberOfPoles; ++k)
75 c0 = c0 * (1.0 - m_SplinePoles[k]) * (1.0 - 1.0 / m_SplinePoles[k]);
79 for (
unsigned int n = 0; n < m_DataLength[m_IteratorDirection]; ++n)
85 for (
int k = 0; k < m_NumberOfPoles; ++k)
88 this->SetInitialCausalCoefficient(m_SplinePoles[k]);
90 for (
unsigned int n = 1; n < m_DataLength[m_IteratorDirection]; ++n)
92 m_Scratch[n] += m_SplinePoles[k] * m_Scratch[n - 1];
96 this->SetInitialAntiCausalCoefficient(m_SplinePoles[k]);
98 for (
int n = m_DataLength[m_IteratorDirection] - 2; 0 <= n; n--)
100 m_Scratch[n] = m_SplinePoles[k] * (m_Scratch[n + 1] - m_Scratch[n]);
106 template <
class TInputImage,
class TOutputImage>
109 if (SplineOrder == m_SplineOrder)
113 m_SplineOrder = SplineOrder;
118 template <
class TInputImage,
class TOutputImage>
124 switch (m_SplineOrder)
128 m_SplinePoles[0] = std::sqrt(3.0) - 2.0;
138 m_SplinePoles[0] = std::sqrt(8.0) - 3.0;
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;
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;
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.");
160 template <
class TInputImage,
class TOutputImage>
166 double sum, zn, z2n, iz;
167 unsigned long horizon;
170 horizon = m_DataLength[m_IteratorDirection];
172 if (m_Tolerance > 0.0)
174 horizon = (long)std::ceil(log(m_Tolerance) / std::log(fabs(z)));
176 if (horizon < m_DataLength[m_IteratorDirection])
180 for (
unsigned int n = 1; n < horizon; ++n)
182 sum += zn * m_Scratch[n];
191 z2n = std::pow(z, (
double)(m_DataLength[m_IteratorDirection] - 1L));
192 sum = m_Scratch[0] + z2n * m_Scratch[m_DataLength[m_IteratorDirection] - 1L];
194 for (
unsigned int n = 1; n <= (m_DataLength[m_IteratorDirection] - 2); ++n)
196 sum += (zn + z2n) * m_Scratch[n];
200 m_Scratch[0] = sum / (1.0 - zn * zn);
204 template <
class TInputImage,
class TOutputImage>
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]);
214 template <
class TInputImage,
class TOutputImage>
219 itk::Size<ImageDimension> size = output->GetBufferedRegion().GetSize();
221 unsigned int count = output->GetBufferedRegion().GetNumberOfPixels() / size[0] * ImageDimension;
223 itk::ProgressReporter progress(
this, 0, count, 10);
226 this->CopyImageToImage();
228 for (
unsigned int n = 0; n < ImageDimension; ++n)
230 m_IteratorDirection = n;
235 CIterator.SetDirection(m_IteratorDirection);
237 while (!CIterator.IsAtEnd())
240 this->CopyCoefficientsToScratch(CIterator);
243 this->DataToCoefficients1D();
247 CIterator.GoToBeginOfLine();
248 this->CopyScratchToCoefficients(CIterator);
249 CIterator.NextLine();
250 progress.CompletedPixel();
258 template <
class TInputImage,
class TOutputImage>
262 typedef itk::ImageRegionConstIteratorWithIndex<TInputImage> InputIterator;
263 typedef itk::ImageRegionIterator<TOutputImage> OutputIterator;
264 typedef typename TOutputImage::PixelType OutputPixelType;
266 InputIterator inIt(this->GetInput(), this->GetInput()->GetBufferedRegion());
267 OutputIterator outIt(this->GetOutput(), this->GetOutput()->GetBufferedRegion());
272 while (!outIt.IsAtEnd())
274 outIt.Set(
static_cast<OutputPixelType
>(inIt.Get()));
283 template <
class TInputImage,
class TOutputImage>
286 typedef typename TOutputImage::PixelType OutputPixelType;
288 while (!Iter.IsAtEndOfLine())
290 Iter.Set(
static_cast<OutputPixelType
>(m_Scratch[j]));
300 template <
class TInputImage,
class TOutputImage>
304 while (!Iter.IsAtEndOfLine())
306 m_Scratch[j] =
static_cast<double>(Iter.Get());
316 template <
class TInputImage,
class TOutputImage>
322 m_DataLength = inputPtr->GetBufferedRegion().GetSize();
324 unsigned long maxLength = 0;
325 for (
unsigned int n = 0; n < ImageDimension; ++n)
327 if (m_DataLength[n] > maxLength)
329 maxLength = m_DataLength[n];
332 m_Scratch.resize(maxLength);
336 outputPtr->SetBufferedRegion(outputPtr->GetRequestedRegion());
337 outputPtr->Allocate();
340 this->DataToCoefficientsND();