GridInterpolation.h

Go to the documentation of this file.

//     ____   ______ __
//    / __ \ / ____// /
//   / /_/ // /    / /
//  / ____// /___ / /___   PixInsight Class Library
// /_/     \____//_____/   PCL 2.6.5
// ----------------------------------------------------------------------------
// pcl/GridInterpolation.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_GridInterpolation_h
#define __PCL_GridInterpolation_h
 
 
#include <pcl/Defs.h>
#include <pcl/Diagnostics.h>
 
#include <pcl/AbstractImage.h>
#include <pcl/BicubicInterpolation.h>
#include <pcl/Matrix.h>
#include <pcl/ParallelProcess.h>
#include <pcl/Rectangle.h>
#include <pcl/ReferenceArray.h>
#include <pcl/Thread.h>
 
#ifndef __PCL_BUILDING_PIXINSIGHT_APPLICATION
#  include <pcl/Console.h>
#  include <pcl/StandardStatus.h>
#endif
 
namespace pcl
{
 
// ----------------------------------------------------------------------------
 
class PCL_CLASS GridInterpolation : public ParallelProcess
{
public:
 
   GridInterpolation() = default;
 
   GridInterpolation( const GridInterpolation& ) = default;
 
   GridInterpolation( GridInterpolation&& ) = default;
 
   GridInterpolation& operator =( const GridInterpolation& ) = default;
 
   GridInterpolation& operator =( GridInterpolation&& ) = default;
 
   template <class R, class SI>
   void Initialize( const R& rect, int delta, const SI& S, bool verbose = true )
   {
      PCL_PRECONDITION( rect.IsRect() )
      PCL_PRECONDITION( delta > 0 )
 
      m_rect = rect.Ordered();
      double width = m_rect.Width();
      double height = m_rect.Height();
      if ( !m_rect.IsRect() || 1 + width == 1 || 1 + height == 1 )
         throw Error( "GridInterpolation::Initialize(): Empty interpolation space." );
 
      m_delta = Abs( delta );
      if ( 1 + m_delta == 1 )
         throw Error( "GridInterpolation::Initialize(): Zero or insignificant grid distance." );
 
      int rows = 1 + int( Ceil( height/m_delta ) );
      int cols = 1 + int( Ceil( width/m_delta ) );
      m_G = DMatrix( rows, cols );
 
      StatusMonitor monitor;
#ifndef __PCL_BUILDING_PIXINSIGHT_APPLICATION
      StandardStatus status;
      if ( verbose )
      {
         monitor.SetCallback( &status );
         monitor.Initialize( "Building surface interpolation grid", rows );
      }
#endif
 
      Array<size_type> L = Thread::OptimalThreadLoads( rows,
                                                       1/*overheadLimit*/,
                                                       m_parallel ? m_maxProcessors : 1 );
      AbstractImage::ThreadData data( monitor, rows );
      ReferenceArray<GridInitializationThread<SI> > threads;
      for ( int i = 0, n = 0; i < int( L.Length() ); n += int( L[i++] ) )
         threads.Add( new GridInitializationThread<SI>( data, *this, S, n, n + int( L[i] ) ) );
      AbstractImage::RunThreads( threads, data );
      threads.Destroy();
 
      m_I.Initialize( m_G.Begin(), cols, rows );
   }
 
   template <class R>
   void Initialize( const R& rect, double delta, const DMatrix& G )
   {
      PCL_PRECONDITION( rect.IsRect() )
      PCL_PRECONDITION( delta > 0 )
 
      m_rect = DRect( rect.Ordered() );
      double width = m_rect.Width();
      double height = m_rect.Height();
      if ( !m_rect.IsRect() || 1 + width == 1 || 1 + height == 1 )
         throw Error( "GridInterpolation::Initialize(): Empty interpolation space." );
 
      m_delta = Abs( delta );
      if ( 1 + m_delta == 1 )
         throw Error( "GridInterpolation::Initialize(): Zero or insignificant grid distance." );
 
      int rows = 1 + int( Ceil( height/m_delta ) );
      int cols = 1 + int( Ceil( width/m_delta ) );
      if ( G.Rows() != rows || G.Cols() != cols )
         throw Error( "GridInterpolation::Initialize(): Invalid matrix dimensions." );
 
      m_G = G;
      m_I.Initialize( m_G.Begin(), cols, rows );
   }
 
   bool IsValid() const
   {
      return !m_G.IsEmpty();
   }
 
   void Clear()
   {
      m_I.Clear();
      m_G.Clear();
   }
 
   const DRect& ReferenceRect() const
   {
      return m_rect;
   }
 
   double Delta() const
   {
      return m_delta;
   }
 
   const DMatrix& InterpolationMatrix() const
   {
      return m_G;
   }
 
   template <typename T>
   double operator ()( T x, T y ) const
   {
      double fx = (double( x ) - m_rect.x0)/m_delta;
      double fy = (double( y ) - m_rect.y0)/m_delta;
      return m_I( fx, fy );
   }
 
   template <typename T>
   double operator ()( const GenericPoint<T>& p ) const
   {
      return operator ()( p.x, p.y );
   }
 
private:
 
   using grid_interpolation = BicubicBSplineInterpolation<double>;
 
   DRect              m_rect;
   double             m_delta;
   DMatrix            m_G;
   grid_interpolation m_I;
 
   template <class SI>
   class GridInitializationThread : public Thread
   {
   public:
 
      GridInitializationThread( const AbstractImage::ThreadData& data,
                                GridInterpolation& grid, const SI& surface, int startRow, int endRow )
         : m_data( data )
         , m_grid( grid )
         , m_surface( surface )
         , m_startRow( startRow )
         , m_endRow( endRow )
      {
      }
 
      PCL_HOT_FUNCTION void Run() override
      {
         INIT_THREAD_MONITOR()
 
         double y = m_grid.m_rect.y0 + m_startRow*m_grid.m_delta;
         for ( int i = m_startRow; i < m_endRow; ++i, y += m_grid.m_delta )
         {
            double x = m_grid.m_rect.x0;
            for ( int j = 0; j < m_grid.m_G.Cols(); ++j, x += m_grid.m_delta )
               m_grid.m_G[i][j] = m_surface( x, y );
 
            UPDATE_THREAD_MONITOR( 32 )
         }
      }
 
   private:
 
      const AbstractImage::ThreadData& m_data;
            GridInterpolation& m_grid;
      const SI&                m_surface;
            int                m_startRow, m_endRow;
   };
};
 
// ----------------------------------------------------------------------------
 
class PCL_CLASS PointGridInterpolation : public ParallelProcess
{
public:
 
   PointGridInterpolation() = default;
 
   PointGridInterpolation( const PointGridInterpolation& ) = default;
 
   PointGridInterpolation( PointGridInterpolation&& ) = default;
 
   PointGridInterpolation& operator =( const PointGridInterpolation& ) = default;
 
   PointGridInterpolation& operator =( PointGridInterpolation&& ) = default;
 
   template <class R, class PSI>
   void Initialize( const R& rect, double delta, const PSI& PS, bool verbose = true )
   {
      PCL_PRECONDITION( rect.IsRect() )
      PCL_PRECONDITION( delta > 0 )
 
      m_rect = DRect( rect.Ordered() );
      double width = m_rect.Width();
      double height = m_rect.Height();
      if ( !m_rect.IsRect() || 1 + width == 1 || 1 + height == 1 )
         throw Error( "PointGridInterpolation::Initialize(): Empty interpolation space." );
 
      m_delta = Abs( delta );
      if ( 1 + m_delta == 1 )
         throw Error( "PointGridInterpolation::Initialize(): Zero or insignificant grid distance." );
 
      int rows = 1 + int( Ceil( height/m_delta ) );
      int cols = 1 + int( Ceil( width/m_delta ) );
      m_Gx = DMatrix( rows, cols );
      m_Gy = DMatrix( rows, cols );
 
      StatusMonitor monitor;
#ifndef __PCL_BUILDING_PIXINSIGHT_APPLICATION
      StandardStatus status;
      if ( verbose )
      {
         monitor.SetCallback( &status );
         monitor.Initialize( "Building surface interpolation grid", rows );
      }
#endif
 
      Array<size_type> L = Thread::OptimalThreadLoads( rows,
                                                       1/*overheadLimit*/,
                                                       m_parallel ? m_maxProcessors : 1 );
      AbstractImage::ThreadData data( monitor, rows );
      ReferenceArray<GridInitializationThread<PSI> > threads;
      for ( int i = 0, n = 0; i < int( L.Length() ); n += int( L[i++] ) )
         threads.Add( new GridInitializationThread<PSI>( data, *this, PS, n, n + int( L[i] ) ) );
      AbstractImage::RunThreads( threads, data );
      threads.Destroy();
 
      m_Ix.Initialize( m_Gx.Begin(), cols, rows );
      m_Iy.Initialize( m_Gy.Begin(), cols, rows );
   }
 
   template <class R, class SI>
   void Initialize( const R& rect, double delta, const SI& Sx, const SI& Sy, bool verbose = true )
   {
      PCL_PRECONDITION( rect.IsRect() )
      PCL_PRECONDITION( delta > 0 )
 
      m_rect = DRect( rect.Ordered() );
      double width = m_rect.Width();
      double height = m_rect.Height();
      if ( !m_rect.IsRect() || 1 + width == 1 || 1 + height == 1 )
         throw Error( "PointGridInterpolation::Initialize(): Empty interpolation space." );
 
      m_delta = Abs( delta );
      if ( 1 + m_delta == 1 )
         throw Error( "PointGridInterpolation::Initialize(): Zero or insignificant grid distance." );
 
      int rows = 1 + int( Ceil( height/m_delta ) );
      int cols = 1 + int( Ceil( width/m_delta ) );
      m_Gx = DMatrix( rows, cols );
      m_Gy = DMatrix( rows, cols );
 
      StatusMonitor monitor;
#ifndef __PCL_BUILDING_PIXINSIGHT_APPLICATION
      StandardStatus status;
      if ( verbose )
      {
         monitor.SetCallback( &status );
         monitor.Initialize( "Building surface interpolation grid", rows );
      }
#endif
 
      Array<size_type> L = Thread::OptimalThreadLoads( rows,
                                                       1/*overheadLimit*/,
                                                       m_parallel ? m_maxProcessors : 1 );
      AbstractImage::ThreadData data( monitor, rows );
      ReferenceArray<GridInitializationXYThread<SI> > threads;
      for ( int i = 0, n = 0; i < int( L.Length() ); n += int( L[i++] ) )
         threads.Add( new GridInitializationXYThread<SI>( data, *this, Sx, Sy, n, n + int( L[i] ) ) );
      AbstractImage::RunThreads( threads, data );
      threads.Destroy();
 
      m_Ix.Initialize( m_Gx.Begin(), cols, rows );
      m_Iy.Initialize( m_Gy.Begin(), cols, rows );
   }
 
   template <class R>
   void Initialize( const R& rect, int delta, const DMatrix& Gx, const DMatrix& Gy )
   {
      PCL_PRECONDITION( rect.IsRect() )
      PCL_PRECONDITION( delta > 0 )
 
      m_rect = DRect( rect.Ordered() );
      double width = m_rect.Width();
      double height = m_rect.Height();
      if ( !m_rect.IsRect() || 1 + width == 1 || 1 + height == 1 )
         throw Error( "PointGridInterpolation::Initialize(): Empty interpolation region." );
 
      m_delta = Abs( delta );
      if ( 1 + m_delta == 1 )
         throw Error( "PointGridInterpolation::Initialize(): Zero or insignificant grid distance." );
 
      int rows = 1 + int( Ceil( height/m_delta ) );
      int cols = 1 + int( Ceil( width/m_delta ) );
      if ( Gx.Rows() != rows || Gx.Cols() != cols || Gy.Rows() != rows || Gy.Cols() != cols )
         throw Error( "PointGridInterpolation::Initialize(): Invalid matrix dimensions." );
 
      m_Gx = Gx;
      m_Gy = Gy;
      m_Ix.Initialize( m_Gx.Begin(), cols, rows );
      m_Iy.Initialize( m_Gy.Begin(), cols, rows );
   }
 
   template <class PSI>
   void ApplyLocalModel( const PSI& PS,
                         const String& message = "Applying local interpolation model",
                         bool verbose = true )
   {
      PCL_PRECONDITION( IsValid() )
      if ( !IsValid() )
         throw Error( "PointGridInterpolation::ApplyLocalModel(): Uninitialized interpolation." );
 
      StatusMonitor monitor;
#ifndef __PCL_BUILDING_PIXINSIGHT_APPLICATION
      StandardStatus status;
      if ( verbose )
      {
         monitor.SetCallback( &status );
         monitor.Initialize( message, m_Gx.Rows() );
      }
#endif
 
      Array<size_type> L = Thread::OptimalThreadLoads( m_Gx.Rows(),
                                                       1/*overheadLimit*/,
                                                       m_parallel ? m_maxProcessors : 1 );
      AbstractImage::ThreadData data( monitor, m_Gx.Rows() );
      ReferenceArray<LocalModelThread<PSI> > threads;
      for ( int i = 0, n = 0; i < int( L.Length() ); n += int( L[i++] ) )
         threads.Add( new LocalModelThread<PSI>( data, *this, PS, n, n + int( L[i] ) ) );
      AbstractImage::RunThreads( threads, data );
      threads.Destroy();
   }
 
   template <class SI>
   void ApplyLocalModel( const SI& Sx, const SI& Sy,
                         const String& message = "Applying local interpolation model",
                         bool verbose = true )
   {
      PCL_PRECONDITION( IsValid() )
      if ( !IsValid() )
         throw Error( "PointGridInterpolation::ApplyLocalModel(): Uninitialized interpolation." );
 
      StatusMonitor monitor;
#ifndef __PCL_BUILDING_PIXINSIGHT_APPLICATION
      StandardStatus status;
      if ( verbose )
      {
         monitor.SetCallback( &status );
         monitor.Initialize( message, m_Gx.Rows() );
      }
#endif
 
      Array<size_type> L = Thread::OptimalThreadLoads( m_Gx.Rows(),
                                                       1/*overheadLimit*/,
                                                       m_parallel ? m_maxProcessors : 1 );
      AbstractImage::ThreadData data( monitor, m_Gx.Rows() );
      ReferenceArray<LocalModelThread2<SI> > threads;
      for ( int i = 0, n = 0; i < int( L.Length() ); n += int( L[i++] ) )
         threads.Add( new LocalModelThread2<SI>( data, *this, Sx, Sy, n, n + int( L[i] ) ) );
      AbstractImage::RunThreads( threads, data );
      threads.Destroy();
   }
 
   bool IsValid() const
   {
      return !m_Gx.IsEmpty() && !m_Gy.IsEmpty();
   }
 
   void Clear()
   {
      m_Ix.Clear();
      m_Iy.Clear();
      m_Gx.Clear();
      m_Gy.Clear();
   }
 
   const DRect& ReferenceRect() const
   {
      return m_rect;
   }
 
   double Delta() const
   {
      return m_delta;
   }
 
   const DMatrix& XInterpolationMatrix() const
   {
      return m_Gx;
   }
 
   const DMatrix& YInterpolationMatrix() const
   {
      return m_Gy;
   }
 
   template <typename T>
   DPoint operator ()( T x, T y ) const
   {
      PCL_PRECONDITION( IsValid() )
      double fx = (double( x ) - m_rect.x0)/m_delta;
      double fy = (double( y ) - m_rect.y0)/m_delta;
      return DPoint( m_Ix( fx, fy ), m_Iy( fx, fy ) );
   }
 
   template <typename T>
   DPoint operator ()( const GenericPoint<T>& p ) const
   {
      return operator ()( p.x, p.y );
   }
 
private:
 
   using grid_interpolation = BicubicBSplineInterpolation<double>;
 
   DRect              m_rect;
   double             m_delta;
   DMatrix            m_Gx, m_Gy;
   grid_interpolation m_Ix, m_Iy;
 
   template <class PSI>
   class GridInitializationThread : public Thread
   {
   public:
 
      GridInitializationThread( const AbstractImage::ThreadData& data,
                                PointGridInterpolation& grid, const PSI& surface, int startRow, int endRow )
         : m_data( data )
         , m_grid( grid )
         , m_surface( surface )
         , m_startRow( startRow )
         , m_endRow( endRow )
      {
      }
 
      PCL_HOT_FUNCTION void Run() override
      {
         INIT_THREAD_MONITOR()
 
         double y = m_grid.m_rect.y0 + m_startRow*m_grid.m_delta;
         for ( int i = m_startRow; i < m_endRow; ++i, y += m_grid.m_delta )
         {
            double x = m_grid.m_rect.x0;
            for ( int j = 0; j < m_grid.m_Gx.Cols(); ++j, x += m_grid.m_delta )
            {
               DPoint p = m_surface( x, y );
               m_grid.m_Gx[i][j] = p.x;
               m_grid.m_Gy[i][j] = p.y;
            }
 
            UPDATE_THREAD_MONITOR( 32 )
         }
      }
 
   private:
 
      const AbstractImage::ThreadData& m_data;
            PointGridInterpolation&    m_grid;
      const PSI&                       m_surface;
            int                        m_startRow, m_endRow;
   };
 
   template <class SI>
   class GridInitializationXYThread : public Thread
   {
   public:
 
      GridInitializationXYThread( const AbstractImage::ThreadData& data,
                                  PointGridInterpolation& grid,
                                  const SI& surfaceX, const SI& surfaceY, int startRow, int endRow )
         : m_data( data )
         , m_grid( grid )
         , m_surfaceX( surfaceX )
         , m_surfaceY( surfaceY )
         , m_startRow( startRow )
         , m_endRow( endRow )
      {
      }
 
      PCL_HOT_FUNCTION void Run() override
      {
         INIT_THREAD_MONITOR()
 
         double y = m_grid.m_rect.y0 + m_startRow*m_grid.m_delta;
         for ( int i = m_startRow; i < m_endRow; ++i, y += m_grid.m_delta )
         {
            double x = m_grid.m_rect.x0;
            for ( int j = 0; j < m_grid.m_Gx.Cols(); ++j, x += m_grid.m_delta )
            {
               m_grid.m_Gx[i][j] = m_surfaceX( x, y );
               m_grid.m_Gy[i][j] = m_surfaceY( x, y );
            }
 
            UPDATE_THREAD_MONITOR( 32 )
         }
      }
 
   private:
 
      const AbstractImage::ThreadData& m_data;
            PointGridInterpolation&    m_grid;
      const SI&                        m_surfaceX;
      const SI&                        m_surfaceY;
            int                        m_startRow, m_endRow;
   };
 
   template <class PSI>
   class LocalModelThread : public Thread
   {
   public:
 
      LocalModelThread( const AbstractImage::ThreadData& data,
                        PointGridInterpolation& grid, const PSI& model, int startRow, int endRow )
         : m_data( data )
         , m_grid( grid )
         , m_model( model )
         , m_startRow( startRow )
         , m_endRow( endRow )
      {
      }
 
      PCL_HOT_FUNCTION void Run() override
      {
         INIT_THREAD_MONITOR()
 
         double y = m_grid.m_rect.y0 + m_startRow*m_grid.m_delta;
         for ( int i = m_startRow; i < m_endRow; ++i, y += m_grid.m_delta )
         {
            double x = m_grid.m_rect.x0;
            for ( int j = 0; j < m_grid.m_Gx.Cols(); ++j, x += m_grid.m_delta )
            {
               DPoint d = m_model( x, y );
               m_grid.m_Gx[i][j] += d.x;
               m_grid.m_Gy[i][j] += d.y;
            }
 
            UPDATE_THREAD_MONITOR( 32 )
         }
      }
 
   private:
 
      const AbstractImage::ThreadData& m_data;
            PointGridInterpolation&    m_grid;
      const PSI&                       m_model;
            int                        m_startRow, m_endRow;
   };
 
   template <class SI>
   class LocalModelThread2 : public Thread
   {
   public:
 
      LocalModelThread2( const AbstractImage::ThreadData& data,
                         PointGridInterpolation& grid, const SI& modelX, const SI& modelY, int startRow, int endRow )
         : m_data( data )
         , m_grid( grid )
         , m_modelX( modelX )
         , m_modelY( modelY )
         , m_startRow( startRow )
         , m_endRow( endRow )
      {
      }
 
      PCL_HOT_FUNCTION void Run() override
      {
         INIT_THREAD_MONITOR()
 
         double y = m_grid.m_rect.y0 + m_startRow*m_grid.m_delta;
         for ( int i = m_startRow; i < m_endRow; ++i, y += m_grid.m_delta )
         {
            double x = m_grid.m_rect.x0;
            for ( int j = 0; j < m_grid.m_Gx.Cols(); ++j, x += m_grid.m_delta )
            {
               m_grid.m_Gx[i][j] += m_modelX( x, y );
               m_grid.m_Gy[i][j] += m_modelY( x, y );
            }
 
            UPDATE_THREAD_MONITOR( 32 )
         }
      }
 
   private:
 
      const AbstractImage::ThreadData& m_data;
            PointGridInterpolation&    m_grid;
      const SI&                        m_modelX, m_modelY;
            int                        m_startRow, m_endRow;
   };
 
   friend class SplineWorldTransformation;
};
 
// ----------------------------------------------------------------------------
 
} // pcl
 
#endif   // __PCL_GridInterpolation_h
 
// ----------------------------------------------------------------------------
// EOF pcl/GridInterpolation.h - Released 2024-01-13T15:47:58Z