ShepardInterpolation.h

Go to the documentation of this file.

//     ____   ______ __
//    / __ \ / ____// /
//   / /_/ // /    / /
//  / ____// /___ / /___   PixInsight Class Library
// /_/     \____//_____/   PCL 2.6.5
// ----------------------------------------------------------------------------
// pcl/ShepardInterpolation.h - Released 2024-01-13T15:47:58Z
// ----------------------------------------------------------------------------
// 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_ShepardInterpolation_h
#define __PCL_ShepardInterpolation_h
 
 
#include <pcl/Defs.h>
#include <pcl/Diagnostics.h>
 
#include <pcl/QuadTree.h>
#include <pcl/Vector.h>
 
#ifndef __PCL_BUILDING_PIXINSIGHT_APPLICATION
#  include <pcl/Console.h>
#endif
 
namespace pcl
{
 
// ----------------------------------------------------------------------------
 
/*
 * Default normalized search radius for local Shepard interpolation. This is an
 * initial search radius relative to the unit circle for the adaptive quadtree
 * search algorithm.
 */
#define __PCL_SHEPARD_DEFAULT_SEARCH_RADIUS  0.10
 
/*
 * Default power parameter for local Shepard interpolation. Larger values tend
 * to yield more accurate interpolation devices. Small powers lead to more
 * approximating (smoothing) devices. The chosen default value is intermediate.
 */
#define __PCL_SHEPARD_DEFAULT_POWER          4
 
/*
 * Default regularization (smoothing) factor for local Shepard interpolation,
 * in the range [0,1). This is a clipping fraction for Winsorization of nearby
 * function values in the point interpolation routine.
 */
#define __PCL_SHEPARD_DEFAULT_REGULARIZATION 0
 
// ----------------------------------------------------------------------------
 
template <typename T>
class PCL_CLASS ShepardInterpolation
{
public:
 
   using vector_type = GenericVector<T>;
 
   using scalar = typename vector_type::scalar;
 
   using search_tree = QuadTree<vector_type>;
 
   using search_rect = typename search_tree::rectangle;
 
   constexpr static int BucketCapacity = 16;
 
   ShepardInterpolation() = default;
 
   ShepardInterpolation( const ShepardInterpolation& ) = delete;
 
   ShepardInterpolation( ShepardInterpolation&& ) = default;
 
   virtual ~ShepardInterpolation()
   {
   }
 
   ShepardInterpolation& operator =( const ShepardInterpolation& ) = delete;
 
   ShepardInterpolation& operator =( ShepardInterpolation&& ) = default;
 
   bool IsValid() const
   {
      return !m_T.IsEmpty();
   }
 
   void SetPower( int m )
   {
      PCL_PRECONDITION( m > 0 )
      m_mu = (m > 0) ? m : __PCL_SHEPARD_DEFAULT_POWER;
   }
 
   int Power() const
   {
      return m_mu;
   }
 
   void SetRadius( double R )
   {
      PCL_PRECONDITION( R > 0 )
      m_R = (R < 0 || 1 + R == 1) ? __PCL_SHEPARD_DEFAULT_SEARCH_RADIUS : R;
   }
 
   double Radius() const
   {
      return m_R;
   }
 
   void SetSmoothing( float r )
   {
      PCL_PRECONDITION( r >= 0 && r < 1 )
      m_r = (r < 0 || r >= 1) ? __PCL_SHEPARD_DEFAULT_REGULARIZATION : r;
   }
 
   float Smoothing() const
   {
      return m_r;
   }
 
   void Initialize( const T* x, const T* y, const T* z, int n )
   {
      DoInitialize( nullptr/*rect*/, x, y, z, n );
   }
 
   void Initialize( const search_rect& rect, const T* x, const T* y, const T* z, int n )
   {
      DoInitialize( &rect, x, y, z, n );
   }
 
   T operator ()( double x, double y ) const
   {
      PCL_PRECONDITION( !m_T.IsEmpty() )
 
      const double dx = m_r0*(x - m_x0);
      const double dy = m_r0*(y - m_y0);
 
      for ( double R = m_R; ; R += m_R )
      {
         int m = 0;
         double R2 = R*R;
         Array<DPoint> V;
         m_T.Search( DRect( dx-R, dy-R, dx+R, dy+R ),
                     [&]( const vector_type& v, void* )
                     {
                        double x = dx - v[0];
                        double y = dy - v[1];
                        double r2 = x*x + y*y;
                        if ( r2 < R2 )
                        {
                           ++m;
                           double w = PowI( 1 - Sqrt( r2 )/R, m_mu );
                           /*
                            * N.B. The equivalent code below is about 400 times
                            * slower than the above call to PowI() for m_mu=16.
                            * Measured on a TR 2990WX.
                            *
                           double w = 1 - Sqrt( r2 )/R;
                           for ( int i = 1; i < m_mu; ++i )
                              w *= w;
                            */
                           V << DPoint( w, w*v[2] );
                        }
                     },
                     nullptr );
         if ( m >= 3 )
         {
            if ( m_r > 0 )
            {
               /*
                * Regularization by Winsorization of the weighted sample.
                */
               int r = Min( TruncInt( m_r * m ), (m >> 1) - (m & 1)^1 );
               if ( r > 0 )
               {
                  V.Sort( []( const DPoint& v1, const DPoint& v2 )
                     {
                        return v1.y < v2.y;
                     } );
                  for ( int i = 0; i < r; ++i )
                  {
                     V[i].y = V[r].y;
                     V[m-i-1].y = V[m-r-1].y;
                  }
               }
            }
            DPoint Wz( 0 );
            for ( const DPoint& v : V )
               Wz += v;
            if ( 1 + Wz.x != 1 )
               return T( Wz.y/Wz.x );
            if ( R >= 1 )
               break; // degenerate!
         }
      }
 
      return 0; // empty!?
   }
 
   template <typename Tp>
   T operator ()( const GenericPoint<Tp>& p ) const
   {
      return operator ()( double( p.x ), double( p.y ) );
   }
 
   void Clear()
   {
      m_T.Clear();
   }
 
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_mu = __PCL_SHEPARD_DEFAULT_POWER;          // power parameter (leveling factor)
   double      m_R  = __PCL_SHEPARD_DEFAULT_SEARCH_RADIUS;  // initial search radius
   float       m_r  = __PCL_SHEPARD_DEFAULT_REGULARIZATION; // regularization (clipping fraction)
   search_tree m_T;        // tree points store input coordinates and function values
 
   void DoInitialize( const search_rect* rect, 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( "ShepardInterpolation::Initialize(): At least three input nodes must be specified." );
 
      Clear();
 
      try
      {
         if ( rect == nullptr )
         {
            /*
             * 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 = x[i] - m_x0;
               double dy = y[i] - m_y0;
               double r = Sqrt( dx*dx + dy*dy );
               if ( r > m_r0 )
                  m_r0 = r;
            }
         }
         else
         {
            m_x0 = rect->CenterX();
            m_y0 = rect->CenterY();
            m_r0 = rect->Diagonal()/2;
         }
 
         if ( 1 + m_r0 == 1 )
            throw Error( "ShepardInterpolation::Initialize(): Empty or insignificant interpolation space." );
         m_r0 = 1/m_r0;
 
         /*
          * Build working vector. Transform coordinates to the unit circle.
          */
         Array<vector_type> V;
         for ( int i = 0; i < n; ++i )
            V << vector_type( m_r0*(x[i] - m_x0), m_r0*(y[i] - m_y0), z[i] );
 
         /*
          * Find and remove duplicate input nodes. Two nodes are equal if their
          * coordinates don't differ more than the machine epsilon for the
          * floating point type T.
          */
         V.Sort(  []( const vector_type& p, const vector_type& q )
                  {
                     return (p[0] != q[0]) ? p[0] < q[0] : p[1] < q[1];
                  } );
         Array<int> remove;
         for ( int i = 0, j = 1; j < n; ++i, ++j )
            if ( Abs( V[i][0] - V[j][0] ) <= std::numeric_limits<T>::epsilon() )
               if ( Abs( V[i][1] - V[j][1] ) <= std::numeric_limits<T>::epsilon() )
                  remove << i;
         if ( !remove.IsEmpty() )
         {
            Array<vector_type> U;
            int i = 0;
            for ( int j : remove )
            {
               for ( ; i < j; ++i )
                  U << V[i];
               ++i;
            }
            for ( ; i < n; ++i )
               U << V[i];
            if ( U.Length() < 3 )
               throw Error( "ShepardInterpolation::Initialize(): Less than three input nodes left after sanitization." );
            V = U;
         }
 
         /*
          * Build the point search tree.
          */
         if ( rect == nullptr )
            m_T.Build( V, BucketCapacity );
         else
         {
            m_T.Build( *rect, V, BucketCapacity );
            if ( m_T.Length() < 3 )
               throw Error( "ShepardInterpolation::Initialize(): Less than three input nodes in the specified search region." );
         }
      }
      catch ( ... )
      {
         Clear();
         throw;
      }
   }
};
 
// ----------------------------------------------------------------------------
 
template <class P = DPoint>
class PCL_CLASS PointShepardInterpolation
{
public:
 
   using point = P;
 
   using point_list = Array<point>;
 
   using surface = ShepardInterpolation<double>;
 
   PointShepardInterpolation() = default;
 
   PointShepardInterpolation( const PointShepardInterpolation& ) = default;
 
   PointShepardInterpolation( PointShepardInterpolation&& ) = default;
 
   PointShepardInterpolation( const point_list& P1, const point_list& P2,
                              int power = __PCL_SHEPARD_DEFAULT_POWER,
                              double radius = __PCL_SHEPARD_DEFAULT_SEARCH_RADIUS )
   {
      Initialize( P1, P2, power, radius );
   }
 
   PointShepardInterpolation( const surface& Sx, const surface& Sy )
   {
      Initialize( Sx, Sy );
   }
 
   PointShepardInterpolation& operator =( const PointShepardInterpolation& ) = delete;
 
   PointShepardInterpolation& operator =( PointShepardInterpolation&& ) = default;
 
   void Initialize( const point_list& P1, const point_list& P2,
                    int power = __PCL_SHEPARD_DEFAULT_POWER,
                    double radius = __PCL_SHEPARD_DEFAULT_SEARCH_RADIUS,
                    float smoothing = __PCL_SHEPARD_DEFAULT_REGULARIZATION )
   {
      PCL_PRECONDITION( P1.Length() >= 3 )
      PCL_PRECONDITION( P1.Length() <= P2.Length() )
      PCL_PRECONDITION( power > 0 )
      PCL_PRECONDITION( radius > 0 )
      PCL_PRECONDITION( smoothing >= 0 && smoothing < 1 )
 
      m_Sx.Clear();
      m_Sy.Clear();
 
      m_Sx.SetPower( power );
      m_Sy.SetPower( power );
 
      m_Sx.SetRadius( radius );
      m_Sy.SetRadius( radius );
 
      m_Sx.SetSmoothing( smoothing );
      m_Sy.SetSmoothing( smoothing );
 
      if ( P1.Length() < 3 || P2.Length() < 3 )
         throw Error( "PointShepardInterpolation::Initialize(): At least three input nodes must be specified." );
 
      if ( P2.Length() < P1.Length() )
         throw Error( "PointShepardInterpolation::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 );
   }
 
private:
 
   /*
    * The surface interpolations in the X and Y plane directions.
    */
   surface m_Sx, m_Sy;
 
   friend class DrizzleData;
   friend class DrizzleDataDecoder;
};
 
// ----------------------------------------------------------------------------
 
}  // pcl
 
#endif   // __PCL_ShepardInterpolation_h
 
// ----------------------------------------------------------------------------
// EOF pcl/ShepardInterpolation.h - Released 2024-01-13T15:47:58Z