BicubicFilterInterpolation.h

Go to the documentation of this file.

//     ____   ______ __
//    / __ \ / ____// /
//   / /_/ // /    / /
//  / ____// /___ / /___   PixInsight Class Library
// /_/     \____//_____/   PCL 2.6.5
// ----------------------------------------------------------------------------
// pcl/BicubicFilterInterpolation.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_BicubicFilterInterpolation_h
#define __PCL_BicubicFilterInterpolation_h
 
 
#include <pcl/Defs.h>
#include <pcl/Diagnostics.h>
 
#include <pcl/BidimensionalInterpolation.h>
#include <pcl/Math.h>
#include <pcl/Utility.h>
#include <pcl/Vector.h>
 
namespace pcl
{
 
// ----------------------------------------------------------------------------
 
class PCL_CLASS CubicFilter
{
public:
 
   CubicFilter( double B, double C )
      : m_B( B )
      , m_C( C )
   {
      m_k31 = ( 12 -  9*m_B -  6*m_C)/6;
      m_k21 = (-18 + 12*m_B +  6*m_C)/6;
      m_k01 = (  6 -  2*m_B         )/6;
      m_k32 = (        -m_B -  6*m_C)/6;
      m_k22 = (       6*m_B + 30*m_C)/6;
      m_k12 = (     -12*m_B - 48*m_C)/6;
      m_k02 = (       8*m_B + 24*m_C)/6;
   }
 
   CubicFilter( const CubicFilter& ) = default;
 
   virtual ~CubicFilter()
   {
   }
 
   CubicFilter& operator =( const CubicFilter& ) = default;
 
   PCL_HOT_FUNCTION
   double operator()( double x ) const noexcept
   {
      if ( x < 0 )
         x = -x;
      return (x < 1) ? m_k01 + x*x*(m_k21 + x*m_k31) :
                       m_k02 + x*(m_k12 + x*(m_k22 + x*m_k32));
   }
 
   double Width() const noexcept
   {
      return 2.0;
   }
 
   virtual String Description() const
   {
      return String().Format( "Cubic filter (B=%.6f, C=%.6f)", m_B, m_C );
   }
 
   virtual CubicFilter* Clone() const
   {
      return new CubicFilter( *this );
   }
 
protected:
 
   double m_k31, m_k21, m_k01, m_k32, m_k22, m_k12, m_k02;
 
private:
 
   double m_B, m_C; // for reference only; not used in calculations
};
 
class PCL_CLASS MitchellNetravaliCubicFilter : public CubicFilter
{
public:
 
   MitchellNetravaliCubicFilter()
      : CubicFilter( 1.0/3, 1.0/3 )
   {
   }
 
   MitchellNetravaliCubicFilter( const MitchellNetravaliCubicFilter& ) = default;
 
   ~MitchellNetravaliCubicFilter() override
   {
   }
 
   String Description() const override
   {
      return "Mitchell-Netravali cubic filter";
   }
 
   CubicFilter* Clone() const override
   {
      return new MitchellNetravaliCubicFilter( *this );
   }
};
 
class PCL_CLASS CatmullRomSplineFilter : public CubicFilter
{
public:
 
   CatmullRomSplineFilter()
      : CubicFilter( 0, 0.5 )
   {
   }
 
   CatmullRomSplineFilter( const CatmullRomSplineFilter& ) = default;
 
   ~CatmullRomSplineFilter() override
   {
   }
 
   String Description() const override
   {
      return "Catmull-Rom spline filter";
   }
 
   CubicFilter* Clone() const override
   {
      return new CatmullRomSplineFilter( *this );
   }
};
 
class PCL_CLASS CubicBSplineFilter : public CubicFilter
{
public:
 
   CubicBSplineFilter()
      : CubicFilter( 1, 0 )
   {
   }
 
   CubicBSplineFilter( const CubicBSplineFilter& ) = default;
 
   ~CubicBSplineFilter() override
   {
   }
 
   String Description() const override
   {
      return "Cubic B-spline filter";
   }
 
   CubicFilter* Clone() const override
   {
      return new CubicBSplineFilter( *this );
   }
};
 
// ----------------------------------------------------------------------------
 
#define m_width      this->m_width
#define m_height     this->m_height
#define m_fillBorder this->m_fillBorder
#define m_fillValue  this->m_fillValue
#define m_data       this->m_data
 
// ----------------------------------------------------------------------------
 
template <typename T>
class PCL_CLASS BicubicFilterInterpolation : public BidimensionalInterpolation<T>
{
public:
 
   BicubicFilterInterpolation( int rx, int ry, const CubicFilter& filter )
      : m_rx( Max( 1, rx ) )
      , m_ry( Max( 1, ry ) )
      , m_filter( filter )
   {
      PCL_PRECONDITION( rx >= 1 )
      PCL_PRECONDITION( ry >= 1 )
      Initialize();
   }
 
   BicubicFilterInterpolation( const BicubicFilterInterpolation& ) = default;
 
   ~BicubicFilterInterpolation() override
   {
   }
 
   void Initialize( const T* data, int dataWidth, int dataHeight ) override
   {
      BidimensionalInterpolation<T>::Initialize( data, dataWidth, dataHeight );
 
      if ( m_rx > m_width || m_ry > m_height )
      {
         m_rx = Min( m_rx, m_width );
         m_ry = Min( m_ry, m_height );
         Initialize();
      }
   }
 
   PCL_HOT_FUNCTION
   double operator()( double x, double y ) const override
   {
      PCL_PRECONDITION( m_data != nullptr )
      PCL_PRECONDITION( m_width > 0 && m_height > 0 )
      PCL_PRECONDITION( x >= 0 && x < m_width )
      PCL_PRECONDITION( y >= 0 && y < m_height )
      PCL_CHECK( m_rx >= 1 && m_rx <= m_width )
      PCL_CHECK( m_ry >= 1 && m_ry <= m_height )
      PCL_CHECK( !m_k.IsEmpty() )
 
      int x0 = Range( RoundInt( x ), 0, m_width-1 );
      int y0 = Range( RoundInt( y ), 0, m_height-1 );
 
      double dx = x - x0;
      double dy = y - y0;
 
      int64 d = int64( y0 - m_ry )*m_width + x0 - m_rx;
 
      if ( m_rx >= m_ry )
      {
         double wx = 0;
         for ( int xi = -m_rx, c = 0; xi <= m_rx; ++xi, ++c )
            wx += (m_k[c] = m_filter( m_sx*(xi - dx) ));
 
         double sy = 0, wy = 0;
         for ( int yi = -m_ry; yi <= m_ry; ++yi )
         {
            double ky = m_filter( m_sy*(yi - dy) );
            wy += ky;
 
            const T* p;
            int y = y0 + yi;
 
            if ( y < 0 )
            {
               d += m_width;
 
               if ( m_fillBorder )
               {
                  sy += m_fillValue * ky;
                  continue;
               }
 
               p = m_data + (d - 2*int64( m_width )*(y + 1));
            }
            else if ( y >= m_height )
            {
               if ( m_fillBorder )
               {
                  sy += m_fillValue * ky;
                  continue;
               }
 
               p = m_data + (d -= m_width);
            }
            else
            {
               p = m_data + d;
               d += m_width;
            }
 
            int x  = x0 - m_rx;
            int x2 = x0 + m_rx;
            int x1 = Min( x2, m_width-1 );
            double sx = 0;
            const double* kx = m_k.Begin();
 
            while ( x < 0 )
            {
               ++x;
               ++p;
               sx += (m_fillBorder ? m_fillValue : double( *(p - (x + x)) )) * *kx++;
            }
 
            for ( ;; )
            {
               sx += *p++ * *kx++;
               if ( ++x > x1 )
                  break;
            }
 
            while ( x <= x2 )
            {
               sx += (m_fillBorder ? m_fillValue : double( *--p )) * *kx++;
               ++x;
            }
 
            sy += sx/wx * ky;
         }
 
         return sy/wy;
      }
      else
      {
         double wy = 0;
         for ( int yi = -m_ry, r = 0; yi <= m_ry; ++yi, ++r )
            wy += (m_k[r] = m_filter( m_sy*(yi - dy) ));
 
         double sx = 0, wx = 0;
         for ( int xi = -m_rx; xi <= m_rx; ++xi )
         {
            double kx = m_filter( m_sx*(xi - dx) );
            wx += kx;
 
            const T* p;
            int x = x0 + xi;
 
            if ( x < 0 )
            {
               ++d;
 
               if ( m_fillBorder )
               {
                  sx += m_fillValue * kx;
                  continue;
               }
 
               p = m_data + (d - 2*(x + 1));
            }
            else if ( x >= m_width )
            {
               if ( m_fillBorder )
               {
                  sx += m_fillValue * kx;
                  continue;
               }
 
               p = m_data + --d;
            }
            else
               p = m_data + d++;
 
            int y  = y0 - m_ry;
            int y2 = y0 + m_ry;
            int y1 = Min( y2, m_height-1 );
            double sy = 0;
            const double* ky = m_k.Begin();
 
            while ( y < 0 )
            {
               ++y;
               p += m_width;
               sy += (m_fillBorder ? m_fillValue : double( *(p - (y + y)*m_width) )) * *ky++;
            }
 
            for ( ;; )
            {
               sy += *p * *ky++;
               p += m_width;
               if ( ++y > y1 )
                  break;
            }
 
            while ( y <= y2 )
            {
               sy += (m_fillBorder ? m_fillValue : double( *(p -= m_width) )) * *ky++;
               ++y;
            }
 
            sx += sy/wy * kx;
         }
 
         return sx/wx;
      }
   }
 
   int HorizontalRadius() const noexcept
   {
      return m_rx;
   }
 
   int VerticalRadius() const noexcept
   {
      return m_ry;
   }
 
   void SetRadii( int rx, int ry )
   {
      if ( rx != m_rx || ry != m_ry )
      {
         m_rx = Max( 1, rx );
         m_ry = Max( 1, ry );
 
         if ( this->data != nullptr )
         {
            if ( m_rx > m_width )
               m_rx = m_width;
            if ( m_ry > m_height )
               m_ry = m_height;
         }
 
         Initialize();
      }
   }
 
   const CubicFilter& Filter() const noexcept
   {
      return m_filter;
   }
 
   void SetFilter( const CubicFilter& filter )
   {
      m_filter = filter;
      Initialize();
   }
 
protected:
 
           int         m_rx, m_ry; // horizontal and vertical radii
           double      m_sx, m_sy; // filter scaling ratios
   mutable DVector     m_k;        // workspace for interpolated values
           CubicFilter m_filter;
 
   void Initialize()
   {
      m_sx = m_filter.Width()/(m_rx + 0.5);
      m_sy = m_filter.Width()/(m_ry + 0.5);
      m_k = DVector( 1 + (Max( m_rx, m_ry ) << 1) );
   }
};
 
// ----------------------------------------------------------------------------
 
#undef m_width
#undef m_height
#undef m_fillBorder
#undef m_fillValue
#undef m_data
 
// ----------------------------------------------------------------------------
 
} // pcl
 
#endif   // __PCL_BicubicFilterInterpolation_h
 
// ----------------------------------------------------------------------------
// EOF pcl/BicubicFilterInterpolation.h - Released 2024-01-13T15:47:58Z