SurfacePolynomial.h

Go to the documentation of this file.

//     ____   ______ __
//    / __ \ / ____// /
//   / /_/ // /    / /
//  / ____// /___ / /___   PixInsight Class Library
// /_/     \____//_____/   PCL 2.8.5
// ----------------------------------------------------------------------------
// pcl/SurfacePolynomial.h - Released 2024-12-28T16:53:48Z
// ----------------------------------------------------------------------------
// This file is part of the PixInsight Class Library (PCL).
// PCL is a multiplatform C++ framework for development of PixInsight modules.
//
// Copyright (c) 2003-2024 Pleiades Astrophoto S.L. All Rights Reserved.
//
// Redistribution and use in both source and binary forms, with or without
// modification, is permitted provided that the following conditions are met:
//
// 1. All redistributions of source code must retain the above copyright
//    notice, this list of conditions and the following disclaimer.
//
// 2. All redistributions in binary form must reproduce the above copyright
//    notice, this list of conditions and the following disclaimer in the
//    documentation and/or other materials provided with the distribution.
//
// 3. Neither the names "PixInsight" and "Pleiades Astrophoto", nor the names
//    of their contributors, may be used to endorse or promote products derived
//    from this software without specific prior written permission. For written
//    permission, please contact info@pixinsight.com.
//
// 4. All products derived from this software, in any form whatsoever, must
//    reproduce the following acknowledgment in the end-user documentation
//    and/or other materials provided with the product:
//
//    "This product is based on software from the PixInsight project, developed
//    by Pleiades Astrophoto and its contributors (https://pixinsight.com/)."
//
//    Alternatively, if that is where third-party acknowledgments normally
//    appear, this acknowledgment must be reproduced in the product itself.
//
// THIS SOFTWARE IS PROVIDED BY PLEIADES ASTROPHOTO AND ITS CONTRIBUTORS
// "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
// TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL PLEIADES ASTROPHOTO OR ITS
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, BUSINESS
// INTERRUPTION; PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; AND LOSS OF USE,
// DATA OR PROFITS) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
// POSSIBILITY OF SUCH DAMAGE.
// ----------------------------------------------------------------------------
 
#ifndef __PCL_SurfacePolynomial_h
#define __PCL_SurfacePolynomial_h
 
 
#include <pcl/Defs.h>
#include <pcl/Diagnostics.h>
 
#include <pcl/Matrix.h>
#include <pcl/Point.h>
#include <pcl/Vector.h>
 
namespace pcl
{
 
// ----------------------------------------------------------------------------
 
template <typename T>
class PCL_CLASS SurfacePolynomial
{
public:
 
   using vector_type = GenericVector<T>;
 
   using scalar = typename vector_type::scalar;
 
   SurfacePolynomial() = default;
 
   SurfacePolynomial( const SurfacePolynomial& ) = default;
 
   SurfacePolynomial( SurfacePolynomial&& ) = default;
 
   virtual ~SurfacePolynomial()
   {
   }
 
   SurfacePolynomial& operator =( const SurfacePolynomial& ) = default;
 
   SurfacePolynomial& operator =( SurfacePolynomial&& ) = default;
 
   bool IsValid() const
   {
      return !m_polynomial.IsEmpty();
   }
 
   int Degree() const
   {
      return m_degree;
   }
 
   void SetDegree( int degree )
   {
      PCL_PRECONDITION( degree >= 1 )
      Clear();
      m_degree = pcl::Max( 1, degree );
   }
 
   void Initialize( const T* x, const T* y, const T* z, int n )
   {
      PCL_PRECONDITION( x != nullptr && y != nullptr && z != nullptr )
      PCL_PRECONDITION( n > 2 )
 
      if ( n < 3 )
         throw Error( "At least three input nodes are required in SurfacePolynomial::Initialize()" );
 
      Clear();
 
      // Find mean coordinate values
      m_x0 = m_y0 = 0;
      for ( int i = 0; i < n; ++i )
      {
         m_x0 += x[i];
         m_y0 += y[i];
      }
      m_x0 /= n;
      m_y0 /= n;
 
      // Find radius of unit circle
      m_r0 = 0;
      for ( int i = 0; i < n; ++i )
      {
         double dx = m_x0 - x[i];
         double dy = m_y0 - y[i];
         double r = Sqrt( dx*dx + dy*dy );
         if ( r > m_r0 )
            m_r0 = r;
      }
      if ( 1 + m_r0 == 1 )
         throw Error( "SurfacePolynomial::Initialize(): Empty or insignificant interpolation space" );
 
      m_r0 = 1/m_r0;
 
      const int size = (m_degree + 1)*(m_degree + 2) >> 1;
 
      DMatrix M( 0.0, size, size );
      DVector R( 0.0, size );
      {
         // Transform coordinates to unit circle
         DVector X( n ), Y( n );
         for ( int i = 0; i < n; ++i )
         {
            X[i] = m_r0*(x[i] - m_x0);
            Y[i] = m_r0*(y[i] - m_y0);
         }
 
         GenericVector<DVector> Z( n );
         for ( int k = 0; k < n; ++k )
         {
            Z[k] = DVector( size );
            for ( int i = 0, l = 0; i <= m_degree; ++i )
               for ( int j = 0; j <= m_degree-i; ++j, ++l )
                  Z[k][l] = PowI( X[k], i ) * PowI( Y[k], j );
         }
 
         int n2 = n*n;
         for ( int i = 0; i < size; ++i )
         {
            for ( int j = 0; j < size; ++j )
            {
               for ( int k = 0; k < n; ++k )
                  M[i][j] += Z[k][i] * Z[k][j];
               M[i][j] /= n2;
            }
 
            for ( int k = 0; k < n; ++k )
               R[i] += z[k] * Z[k][i];
            R[i] /= n2;
         }
      }
 
      for ( int i = 0; i < size; ++i )
      {
         double pMi = M[i][i];
         if ( M[i][i] != 0 )
         {
            for ( int j = i; j < size; ++j )
               M[i][j] /= pMi;
            R[i] /= pMi;
         }
 
         for ( int k = 1; i+k < size; ++k )
         {
            double pMk = M[i+k][i];
            if ( M[i+k][i] != 0 )
            {
               for ( int j = i; j < size; ++j )
               {
                  M[i+k][j] /= pMk;
                  M[i+k][j] -= M[i][j];
               }
               R[i+k] /= pMk;
               R[i+k] -= R[i];
            }
         }
      }
 
      m_polynomial = vector_type( size );
      for ( int i = size; --i >= 0; )
      {
         m_polynomial[i] = scalar( R[i] );
         for ( int j = i; --j >= 0; )
            R[j] -= M[j][i]*R[i];
      }
   }
 
   T operator ()( double x, double y ) const
   {
      PCL_PRECONDITION( !m_polynomial.IsEmpty() )
 
      double dx = m_r0*(x - m_x0);
      double dy = m_r0*(y - m_y0);
      double z = 0;
      double px = 1;
      for ( int i = 0, l = 0; ; px *= dx )
      {
         double py = 1;
         for ( int j = 0; j <= m_degree-i; py *= dy, ++j, ++l )
            z += m_polynomial[l]*px*py;
         if ( ++i > m_degree )
            break;
      }
      return T( z );
   }
 
   void Clear()
   {
      m_polynomial.Clear();
   }
 
   template <typename T1>
   void Evaluate( T1* Z, const T1* X, const T1* Y, size_type n ) const
   {
      PCL_PRECONDITION( !m_polynomial.IsEmpty() )
      PCL_PRECONDITION( X != nullptr && Y != nullptr && Z != nullptr )
      PCL_PRECONDITION( n > 0 )
      for ( size_type i = 0; i < n; ++i )
         Z[i] = T1( operator()( double( X[i] ), double( Y[i] ) ) );
   }
 
   bool HasFastVectorEvaluation() const
   {
      return false;
   }
 
protected:
 
   double      m_r0 = 1;     // scaling factor for normalization of node coordinates
   double      m_x0 = 0;     // zero offset for normalization of X node coordinates
   double      m_y0 = 0;     // zero offset for normalization of Y node coordinates
   int         m_degree = 3; // polynomial degree > 0
   vector_type m_polynomial; // coefficients of the 2-D surface polynomial
};
 
// ----------------------------------------------------------------------------
 
template <class P = DPoint>
class PCL_CLASS PointSurfacePolynomial
{
public:
 
   using point = P;
 
   using point_list = Array<point>;
 
   using surface = SurfacePolynomial<double>;
 
   PointSurfacePolynomial() = default;
 
   PointSurfacePolynomial( const PointSurfacePolynomial& ) = default;
 
   PointSurfacePolynomial( PointSurfacePolynomial&& ) = default;
 
   PointSurfacePolynomial( const point_list& P1, const point_list& P2, int degree = 3 )
   {
      Initialize( P1, P2, degree );
   }
 
   PointSurfacePolynomial( const surface& Sx, const surface& Sy )
   {
      Initialize( Sx, Sy );
   }
 
   PointSurfacePolynomial& operator =( const PointSurfacePolynomial& ) = default;
 
   PointSurfacePolynomial& operator =( PointSurfacePolynomial&& ) = default;
 
   void Initialize( const point_list& P1, const point_list& P2, int degree = 3 )
   {
      PCL_PRECONDITION( P1.Length() >= 3 )
      PCL_PRECONDITION( P1.Length() <= P2.Length() )
      PCL_PRECONDITION( degree > 0 )
 
      m_Sx.Clear();
      m_Sy.Clear();
 
      m_Sx.SetDegree( degree );
      m_Sy.SetDegree( degree );
 
      if ( P1.Length() < 3 || P2.Length() < 3 )
         throw Error( "PointSurfacePolynomial::Initialize(): At least three input nodes must be specified." );
 
      if ( P2.Length() < P1.Length() )
         throw Error( "PointSurfacePolynomial::Initialize(): Incompatible point array lengths." );
 
      DVector X( int( P1.Length() ) ),
              Y( int( P1.Length() ) ),
              Zx( int( P1.Length() ) ),
              Zy( int( P1.Length() ) );
      for ( int i = 0; i < X.Length(); ++i )
      {
          X[i] = P1[i].x;
          Y[i] = P1[i].y;
         Zx[i] = P2[i].x;
         Zy[i] = P2[i].y;
      }
      m_Sx.Initialize( X.Begin(), Y.Begin(), Zx.Begin(), X.Length() );
      m_Sy.Initialize( X.Begin(), Y.Begin(), Zy.Begin(), X.Length() );
   }
 
   void Clear()
   {
      m_Sx.Clear();
      m_Sy.Clear();
   }
 
   bool IsValid() const
   {
      return m_Sx.IsValid() && m_Sy.IsValid();
   }
 
   const surface& SurfaceX() const
   {
      return m_Sx;
   }
 
   const surface& SurfaceY() const
   {
      return m_Sy;
   }
 
   template <typename T>
   DPoint operator ()( T x, T y ) const
   {
      return DPoint( m_Sx( x, y ), m_Sy( x, y ) );
   }
 
   template <typename T>
   DPoint operator ()( const GenericPoint<T>& p ) const
   {
      return operator ()( p.x, p.y );
   }
 
   template <typename T>
   void Evaluate( T* ZX, T* ZY, const T* X, const T* Y, size_type n ) const
   {
      PCL_PRECONDITION( ISValid() )
      m_Sx.Evaluate( ZX, X, Y, n );
      m_Sy.Evaluate( ZY, X, Y, n );
   }
 
   bool HasFastVectorEvaluation() const
   {
      return false;
   }
 
private:
 
   /*
    * The surface interpolations in the X and Y plane directions.
    */
   surface m_Sx, m_Sy;
};
 
// ----------------------------------------------------------------------------
 
}  // pcl
 
#endif   // __PCL_SurfacePolynomial_h
 
// ----------------------------------------------------------------------------
// EOF pcl/SurfacePolynomial.h - Released 2024-12-28T16:53:48Z