PCL
ChebyshevFit.h
Go to the documentation of this file.
1 // ____ ______ __
2 // / __ \ / ____// /
3 // / /_/ // / / /
4 // / ____// /___ / /___ PixInsight Class Library
5 // /_/ \____//_____/ PCL 2.1.19
6 // ----------------------------------------------------------------------------
7 // pcl/ChebyshevFit.h - Released 2019-11-07T10:59:34Z
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-2019 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 (http://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_ChebyshevFit_h
53 #define __PCL_ChebyshevFit_h
54 
56 
57 #include <pcl/Defs.h>
58 
59 #include <pcl/Constants.h>
60 #include <pcl/ErrorHandler.h>
61 #include <pcl/Exception.h>
62 #include <pcl/MultiVector.h>
63 
64 namespace pcl
65 {
66 
67 // ----------------------------------------------------------------------------
68 
69 template <typename Tx, typename Ty> class GenericScalarChebyshevFit;
70 
100 template <typename Tx, typename Ty>
102 {
103 public:
104 
109 
114 
119 
124 
159  template <class F>
160  GenericChebyshevFit( F f, Tx x1, Tx x2, int N, int n ) :
161  dx( Abs( x2 - x1 ) ),
162  x0( (x1 + x2)/2 ),
163  m( Max( 2, n ), Max( 1, N ) )
164  {
165  PCL_PRECONDITION( N > 0 )
166  PCL_PRECONDITION( n > 1 )
167  if ( 1 + dx == 1 )
168  throw Error( "GenericChebyshevFit: Empty or insignificant function evaluation interval." );
169  N = m.Length();
170  n = m[0];
171  Tx dx2 = dx/2;
172  function_values y( n, N );
173  for ( int j = 0; j < n; ++j )
174  y[j] = f( x0 + Cos( Const<Ty>::pi()*(j + 0.5)/n )*dx2 );
175  c = coefficient_series( N, n );
176  Ty k = 2.0/n;
177  for ( int i = 0; i < N; ++i )
178  for ( int j = 0; j < n; ++j )
179  {
180  Ty s = Ty( 0 );
181  for ( int k = 0; k < n; ++k )
182  s += y[k][i] * Cos( Const<Ty>::pi()*j*(k + 0.5)/n );
183  c[i][j] = k*s;
184  }
185  }
186 
200  GenericChebyshevFit( const coefficient_series& ck, Tx x1, Tx x2 ) :
201  dx( Abs( x2 - x1 ) ),
202  x0( (x1 + x2)/2 ),
203  c( ck ),
204  m( int( ck.Length() ) )
205  {
206  if ( 1 + dx == 1 )
207  throw Error( "GenericChebyshevFit: Empty or insignificant function evaluation interval." );
208  if ( c.IsEmpty() )
209  throw Error( "GenericChebyshevFit: Empty polynomial expansion." );
210  for ( int i = 0; i < m.Length(); ++i )
211  if ( (m[i] = c[i].Length()) < 1 )
212  throw Error( "GenericChebyshevFit: Invalid coefficient series for dimension " + String( i ) + '.' );
213  }
214 
225  GenericChebyshevFit() = default;
226 
230  GenericChebyshevFit( const GenericChebyshevFit& ) = default;
231 
236 
241 
246 
252  bool IsValid() const
253  {
254  return !c.IsEmpty();
255  }
256 
265  Tx LowerBound() const
266  {
267  return x0 - dx/2;
268  }
269 
278  Tx UpperBound() const
279  {
280  return x0 + dx/2;
281  }
282 
288  int NumberOfComponents() const
289  {
290  return m.Length();
291  }
292 
301  int Length( int i = 0 ) const
302  {
303  PCL_PRECONDITION( i >= 0 && i < NumberOfComponents() )
304  return c[i].Length();
305  }
306 
316  int TruncatedLength( int i = -1 ) const
317  {
318  PCL_PRECONDITION( i < 0 || i < NumberOfComponents() )
319  if ( i < 0 )
320  return m.MaxComponent();
321  return m[i];
322  }
323 
331  {
332  int N = 0;
333  for ( int i = 0; i < NumberOfComponents(); ++i )
334  N += c[i].Length();
335  return N;
336  }
337 
346  {
347  return m.Sum();
348  }
349 
359  bool IsTruncated( int i = -1 ) const
360  {
361  PCL_PRECONDITION( i < 0 || i < NumberOfComponents() )
362  if ( i < 0 )
363  return m.MaxComponent() < Length();
364  return m[i] < c[i].Length();
365  }
366 
377  Ty TruncationError( int i = -1 ) const
378  {
379  PCL_PRECONDITION( i < 0 || i < NumberOfComponents() )
380  if ( i < 0 )
381  {
382  int N = NumberOfComponents();
383  function_value e( N );
384  for ( int j = 0; j < N; ++j )
385  e[j] = TruncationError( j );
386  return e.MaxComponent();
387  }
388  if ( m[i] < c[i].Length() )
389  {
390  Ty e = Ty( 0 );
391  for ( int j = c[i].Length(); --j >= m[i]; )
392  e += Abs( c[i][j] );
393  return e;
394  }
395  return Abs( c[i][c[i].Length()-1] );
396  }
397 
428  bool Truncate( Ty e )
429  {
430  e = Abs( e );
431  int N = NumberOfComponents();
432  int tc = 0;
433  for ( int i = 0; i < N; ++i )
434  {
435  Ty s = Ty( 0 );
436  for ( m[i] = c[i].Length(); m[i] > 2; --m[i] )
437  if ( (s += Abs( c[i][m[i]-1] )) >= e )
438  break;
439  if ( m[i] < c[i].Length() )
440  ++tc;
441  }
442  return tc == N;
443  }
444 
449  const coefficient_series& Coefficients() const
450  {
451  return c;
452  }
453 
472  function_value Evaluate( Tx x ) const
473  {
474  PCL_PRECONDITION( x >= LowerBound() )
475  PCL_PRECONDITION( x <= UpperBound() )
476  const Ty y0 = Ty( 2*(x - x0)/dx );
477  const Ty y2 = 2*y0;
478  function_value y( NumberOfComponents() );
479  for ( int i = 0; i < y.Length(); ++i )
480  {
481  Ty d0 = Ty( 0 );
482  Ty d1 = Ty( 0 );
483  const Ty* k = c[i].At( m[i] );
484  for ( int j = m[i]; --j > 0; )
485  {
486  Ty d = d1;
487  d1 = y2*d1 - d0 + *--k;
488  d0 = d;
489  }
490  y[i] = y0*d1 - d0 + *--k/2;
491  }
492  return y;
493  }
494 
498  function_value operator ()( Tx x ) const
499  {
500  return Evaluate( x );
501  }
502 
518  {
519  int N = NumberOfComponents();
521  ch1.dx = dx;
522  ch1.x0 = x0;
523  ch1.c = coefficient_series( N );
524  ch1.m = IVector( N );
525  for ( int i = 0; i < N; ++i )
526  {
527  int n = Max( 1, c[i].Length()-1 );
528  ch1.c[i] = coefficients( n );
529  ch1.m[i] = n;
530  if ( n > 1 )
531  {
532  ch1.c[i][n-1] = 2*n*c[i][n];
533  ch1.c[i][n-2] = 2*(n-1)*c[i][n-1];
534  for ( int j = n-3; j >= 0; --j )
535  ch1.c[i][j] = ch1.c[i][j+2] + 2*(j+1)*c[i][j+1];
536  ch1.c[i] *= Ty( 2 )/dx;
537  }
538  else
539  ch1.c[i] = Ty( 0 );
540  }
541  return ch1;
542  }
543 
560  {
561  int N = NumberOfComponents();
563  ch.dx = dx;
564  ch.x0 = x0;
565  ch.c = coefficient_series( N );
566  ch.m = IVector( N );
567  for ( int i = 0; i < N; ++i )
568  {
569  int n = c[i].Length();
570  ch.c[i] = coefficients( n );
571  ch.m[i] = n;
572  if ( n > 1 )
573  {
574  const Ty k = Ty( dx )/4;
575  Ty s = Ty( 0 );
576  int f = 1;
577  for ( int j = 1; j < n-1; ++j, f = -f )
578  s += f*(ch.c[i][j] = k*(c[i][j-1] - c[i][j+1])/j);
579  ch.c[i][0] = 2*(s + f*(ch.c[i][n-1] = k*c[i][n-2]/(n-1)));
580  }
581  else
582  ch.c[i] = Ty( 0 );
583  }
584  return ch;
585  }
586 
587 private:
588 
589  Tx dx; // x2 - x1
590  Tx x0; // (x1 + x2)/2
591  coefficient_series c; // Chebyshev polynomial coefficients
592  IVector m; // length of the truncated coefficient series
593 
594  friend class GenericScalarChebyshevFit<Tx,Ty>;
595 };
596 
597 // ----------------------------------------------------------------------------
598 
610 template <typename Tx, typename Ty>
612 {
613 public:
614 
622  template <class F>
623  GenericScalarChebyshevFit( F f, Tx x1, Tx x2, int n ) :
624  GenericChebyshevFit<Tx,Ty>( [f]( Tx x ){ return GenericVector<Ty>( f( x ), 1 ); }, x1, x2, 1, n )
625  {
626  }
627 
638  GenericScalarChebyshevFit() = default;
639 
644 
649 
654 
659 
666  int TruncatedLength() const
667  {
669  }
670 
676  bool IsTruncated() const
677  {
679  }
680 
687  Ty TruncationError() const
688  {
690  }
691 
699  Ty Evaluate( Tx x ) const
700  {
702  }
703 
707  Ty operator ()( Tx x ) const
708  {
709  return Evaluate( x );
710  }
711 
721  {
723  }
724 
734  {
736  }
737 
738 private:
739 
746  {
747  this->dx = T.dx;
748  this->x0 = T.x0;
749  this->c = typename GenericChebyshevFit<Tx,Ty>::coefficient_series( 1 );
750  this->c[0] = T.c[0];
751  this->m = IVector( 1 );
752  this->m[0] = T.m[0];
753  }
754 };
755 
756 // ----------------------------------------------------------------------------
757 
758 #ifndef __PCL_NO_CHEBYSHEV_FIT_INSTANTIATE
759 
773 
782 typedef F32ChebyshevFit FChebyshevFit;
783 
793 
802 typedef F64ChebyshevFit DChebyshevFit;
803 
812 typedef DChebyshevFit ChebyshevFit;
813 
814 #ifndef _MSC_VER
815 
829 
842 typedef F80ChebyshevFit LDChebyshevFit;
843 
844 #endif // !_MSC_VER
845 
855 
864 typedef F32ScalarChebyshevFit FScalarChebyshevFit;
865 
875 
884 typedef F64ScalarChebyshevFit DScalarChebyshevFit;
885 
894 typedef DScalarChebyshevFit ScalarChebyshevFit;
895 
896 #ifndef _MSC_VER
897 
911 
925 typedef F80ScalarChebyshevFit LDScalarChebyshevFit;
926 
927 #endif // !_MSC_VER
928 
929 #endif // !__PCL_NO_CHEBYSHEV_FIT_INSTANTIATE
930 
931 // ----------------------------------------------------------------------------
932 
933 } // pcl
934 
935 #endif // __PCL_ChebyshevFit_h
936 
937 // ----------------------------------------------------------------------------
938 // EOF pcl/ChebyshevFit.h - Released 2019-11-07T10:59:34Z
bool IsEmpty() const
Definition: Array.h:311
Fundamental numeric constants.
Definition: Constants.h:72
Approximation of vector-valued functions by Chebyshev polynomial expansions.
Definition: ChebyshevFit.h:101
function_value operator()(Tx x) const
Definition: ChebyshevFit.h:498
80-bit extended precision floating point scalar Chebyshev function approximation. ...
64-bit floating point Chebyshev function approximation.
80-bit extended precision floating point scalar Chebyshev function approximation. ...
64-bit floating point Chebyshev function approximation.
PCL root namespace.
Definition: AbstractImage.h:76
64-bit floating point scalar Chebyshev function approximation.
GenericVector< Ty > function_value
Definition: ChebyshevFit.h:118
const coefficient_series & Coefficients() const
Definition: ChebyshevFit.h:449
int Length() const
Definition: Vector.h:1502
64-bit floating point Chebyshev function approximation.
GenericChebyshevFit & operator=(const GenericChebyshevFit &)=default
32-bit floating point scalar Chebyshev function approximation.
GenericVector< Ty > coefficients
Definition: ChebyshevFit.h:108
Generic vector of arbitrary length.
Definition: Vector.h:106
int TruncatedLength(int i=-1) const
Definition: ChebyshevFit.h:316
function_value Evaluate(Tx x) const
Definition: ChebyshevFit.h:472
T Abs(const Complex< T > &c)
Definition: Complex.h:420
Unicode (UTF-16) string.
Definition: String.h:7911
Approximation of scalar-valued functions by Chebyshev polynomial expansion.
Definition: ChebyshevFit.h:69
constexpr const T & Max(const T &a, const T &b)
Definition: Utility.h:119
GenericChebyshevFit(const coefficient_series &ck, Tx x1, Tx x2)
Definition: ChebyshevFit.h:200
int NumberOfTruncatedCoefficients() const
Definition: ChebyshevFit.h:345
int NumberOfCoefficients() const
Definition: ChebyshevFit.h:330
64-bit floating point scalar Chebyshev function approximation.
GenericMultiVector< Ty > function_values
Definition: ChebyshevFit.h:123
component MaxComponent() const
Definition: Vector.h:974
int Length(int i=0) const
Definition: ChebyshevFit.h:301
A simple exception with an associated error message.
Definition: Exception.h:213
80-bit extended precision floating point Chebyshev function approximation.
GenericChebyshevFit Integral() const
Definition: ChebyshevFit.h:559
size_type Length() const
Definition: Array.h:265
bool IsTruncated(int i=-1) const
Definition: ChebyshevFit.h:359
GenericScalarChebyshevFit Integral() const
Definition: ChebyshevFit.h:733
GenericScalarChebyshevFit(F f, Tx x1, Tx x2, int n)
Definition: ChebyshevFit.h:623
Ty TruncationError(int i=-1) const
Definition: ChebyshevFit.h:377
Complex< T > Cos(const Complex< T > &c)
Definition: Complex.h:797
32-bit signed integer vector.
32-bit floating point Chebyshev function approximation.
int NumberOfComponents() const
Definition: ChebyshevFit.h:288
64-bit floating point scalar Chebyshev function approximation.
80-bit extended precision floating point Chebyshev function approximation.
GenericChebyshevFit Derivative() const
Definition: ChebyshevFit.h:517
GenericChebyshevFit(F f, Tx x1, Tx x2, int N, int n)
Definition: ChebyshevFit.h:160
iterator At(size_type i)
Definition: Array.h:354
32-bit floating point Chebyshev function approximation.
32-bit floating point scalar Chebyshev function approximation.
GenericMultiVector< Ty > coefficient_series
Definition: ChebyshevFit.h:113
GenericScalarChebyshevFit Derivative() const
Definition: ChebyshevFit.h:720