PCL
AkimaInterpolation.h
Go to the documentation of this file.
1 // ____ ______ __
2 // / __ \ / ____// /
3 // / /_/ // / / /
4 // / ____// /___ / /___ PixInsight Class Library
5 // /_/ \____//_____/ PCL 2.7.0
6 // ----------------------------------------------------------------------------
7 // pcl/AkimaInterpolation.h - Released 2024-06-18T15:48:54Z
8 // ----------------------------------------------------------------------------
9 // This file is part of the PixInsight Class Library (PCL).
10 // PCL is a multiplatform C++ framework for development of PixInsight modules.
11 //
12 // Copyright (c) 2003-2024 Pleiades Astrophoto S.L. All Rights Reserved.
13 //
14 // Redistribution and use in both source and binary forms, with or without
15 // modification, is permitted provided that the following conditions are met:
16 //
17 // 1. All redistributions of source code must retain the above copyright
18 // notice, this list of conditions and the following disclaimer.
19 //
20 // 2. All redistributions in binary form must reproduce the above copyright
21 // notice, this list of conditions and the following disclaimer in the
22 // documentation and/or other materials provided with the distribution.
23 //
24 // 3. Neither the names "PixInsight" and "Pleiades Astrophoto", nor the names
25 // of their contributors, may be used to endorse or promote products derived
26 // from this software without specific prior written permission. For written
27 // permission, please contact info@pixinsight.com.
28 //
29 // 4. All products derived from this software, in any form whatsoever, must
30 // reproduce the following acknowledgment in the end-user documentation
31 // and/or other materials provided with the product:
32 //
33 // "This product is based on software from the PixInsight project, developed
34 // by Pleiades Astrophoto and its contributors (https://pixinsight.com/)."
35 //
36 // Alternatively, if that is where third-party acknowledgments normally
37 // appear, this acknowledgment must be reproduced in the product itself.
38 //
39 // THIS SOFTWARE IS PROVIDED BY PLEIADES ASTROPHOTO AND ITS CONTRIBUTORS
40 // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
41 // TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
42 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL PLEIADES ASTROPHOTO OR ITS
43 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
44 // EXEMPLARY OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, BUSINESS
45 // INTERRUPTION; PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; AND LOSS OF USE,
46 // DATA OR PROFITS) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
47 // CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
48 // ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
49 // POSSIBILITY OF SUCH DAMAGE.
50 // ----------------------------------------------------------------------------
51 
52 #ifndef __PCL_AkimaInterpolation_h
53 #define __PCL_AkimaInterpolation_h
54 
56 
57 #include <pcl/Defs.h>
58 #include <pcl/Diagnostics.h>
59 
61 
62 namespace pcl
63 {
64 
65 // ----------------------------------------------------------------------------
66 
67 #define m_x this->m_x
68 #define m_y this->m_y
69 
95 template <typename T = double>
97 {
98 public:
99 
104 
109 
113  AkimaInterpolation() = default;
114 
119 
124 
128  AkimaInterpolation& operator =( const AkimaInterpolation& ) = default;
129 
133  AkimaInterpolation& operator =( AkimaInterpolation&& ) = default;
134 
139  {
140  }
141 
158  void Initialize( const vector_type& x, const vector_type& y ) override
159  {
160  if ( y.Length() < 5 )
161  throw Error( "AkimaInterpolation::Initialize(): Less than five data points specified." );
162 
163  try
164  {
165  Clear();
167 
168  int n = m_y.Length();
169  int N = n-1; // Number of subintervals
170 
171  m_b = coefficient_vector( N );
172  m_c = coefficient_vector( N );
173  m_d = coefficient_vector( N );
174 
175  // Chordal slopes
176  coefficient_vector m0( N+4 ); // room for 4 additional prescribed slopes
177  T* m = m0.At( 2 ); // allow negative subscripts m[-1] and m[-2]
178 
179  // Akima left-hand slopes to support corners
180  coefficient_vector tL( n );
181 
182  // Calculate chordal slopes for each subinterval
183  if ( m_x )
184  for ( int i = 0; i < N; ++i )
185  {
186  T h = m_x[i+1] - m_x[i];
187  if ( 1 + h*h == 1 )
188  throw Error( "AkimaInterpolation::Initialize(): Empty interpolation subinterval(s)." );
189  m[i] = (m_y[i+1] - m_y[i])/h;
190  }
191  else
192  for ( int i = 0; i < N; ++i )
193  m[i] = m_y[i+1] - m_y[i];
194 
195  // Prescribed slopes at ending locations
196  m[-2 ] = 3*m[ 0] - 2*m[ 1];
197  m[-1 ] = 2*m[ 0] - m[ 1];
198  m[ N ] = 2*m[N-1] - m[N-2];
199  m[N+1] = 3*m[N-1] - 2*m[N-2];
200 
201  /*
202  * Akima left-hand and right-hand slopes.
203  * Right-hand slopes are just the interpolation coefficients bi.
204  */
205  for ( int i = 0; i < n; ++i )
206  {
207  T f = Abs( m[i-1] - m[i-2] );
208  T e = Abs( m[i+1] - m[i] ) + f;
209  if ( 1 + e != 1 )
210  {
211  tL[i] = m[i-1] + f*(m[i] - m[i-1])/e;
212  if ( i < N )
213  m_b[i] = tL[i];
214  }
215  else
216  {
217  tL[i] = m[i-1];
218  if ( i < N )
219  m_b[i] = m[i];
220  }
221  }
222 
223  /*
224  * Interpolation coefficients m_b[i], m_c[i], m_d[i]. 'ai'
225  * coefficients are the m_y[i] ordinate values.
226  */
227  for ( int i = 0; i < N; ++i )
228  {
229  m_c[i] = 3*m[i] - 2*m_b[i] - tL[i+1];
230  m_d[i] = m_b[i] + tL[i+1] - 2*m[i];
231 
232  if ( m_x )
233  {
234  T h = m_x[i+1] - m_x[i];
235  m_c[i] /= h;
236  m_d[i] /= h*h;
237  }
238  }
239  }
240  catch ( ... )
241  {
242  Clear();
243  throw;
244  }
245  }
246 
254  PCL_HOT_FUNCTION
255  double operator()( double x ) const override
256  {
257  PCL_PRECONDITION( IsValid() )
258 
259  /*
260  * Find the subinterval i0 such that m_x[i0] <= x < m_x[i0+1].
261  * Find the distance dx = x - m_x[i], or dx = x - i for implicit x = {0,1,...n-1}.
262  */
263  int i0;
264  double dx;
265  if ( m_x )
266  {
267  i0 = 0;
268  int i1 = m_x.Length() - 1;
269  // if trying to extrapolate, return zero.
270  if ( unlikely( x < m_x[0] || x > m_x[i1] ) )
271  return 0;
272  while ( i1-i0 > 1 )
273  {
274  int im = (i0 + i1) >> 1;
275  if ( x < m_x[im] )
276  i1 = im;
277  else
278  i0 = im;
279  }
280  dx = x - double( m_x[i0] );
281  }
282  else
283  {
284  if ( x <= 0 )
285  return m_y[0];
286  if ( x >= m_y.Length()-1 )
287  return m_y[m_y.Length()-1];
288  i0 = TruncInt( x );
289  dx = x - i0;
290  }
291 
292  /*
293  * Use a Horner scheme to calculate b[i]*dx + c[i]*dx^2 + d[i]*dx^3.
294  */
295  return m_y[i0] + dx*(m_b[i0] + dx*(m_c[i0] + dx*m_d[i0]));
296  }
297 
301  void Clear() override
302  {
303  m_b.Clear();
304  m_c.Clear();
305  m_d.Clear();
307  }
308 
313  bool IsValid() const override
314  {
315  return m_b && m_c && m_d;
316  }
317 
318 protected:
319 
320  /*
321  * Interpolating coefficients for each subinterval.
322  * The coefficients for dx^0 are the input ordinate values in the m_y vector.
323  */
324  coefficient_vector m_b; // coefficients for dx^1
325  coefficient_vector m_c; // coefficients for dx^2
326  coefficient_vector m_d; // coefficients for dx^3
327 };
328 
329 #undef m_x
330 #undef m_y
331 
332 // ----------------------------------------------------------------------------
333 
334 } // pcl
335 
336 #endif // __PCL_AkimaInterpolation_h
337 
338 // ----------------------------------------------------------------------------
339 // EOF pcl/AkimaInterpolation.h - Released 2024-06-18T15:48:54Z
Akima subspline interpolation algorithm.
typename UnidimensionalInterpolation< T >::vector_type vector_type
PCL_HOT_FUNCTION double operator()(double x) const override
bool IsValid() const override
void Initialize(const vector_type &x, const vector_type &y) override
AkimaInterpolation(const AkimaInterpolation &)=default
AkimaInterpolation(AkimaInterpolation &&)=default
A simple exception with an associated error message.
Definition: Exception.h:239
Generic vector of arbitrary length.
Definition: Vector.h:107
A generic interface to one-dimensional interpolation algorithms.
virtual void Initialize(const vector_type &x, const vector_type &y)
T Abs(const Complex< T > &c) noexcept
Definition: Complex.h:429
int TruncInt(T x) noexcept
Definition: Math.h:1132
PCL root namespace.
Definition: AbstractImage.h:77