PCL
BicubicFilterInterpolation.h
Go to the documentation of this file.
1 // ____ ______ __
2 // / __ \ / ____// /
3 // / /_/ // / / /
4 // / ____// /___ / /___ PixInsight Class Library
5 // /_/ \____//_____/ PCL 2.7.0
6 // ----------------------------------------------------------------------------
7 // pcl/BicubicFilterInterpolation.h - Released 2024-06-18T15:48:54Z
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-2024 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 (https://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_BicubicFilterInterpolation_h
53 #define __PCL_BicubicFilterInterpolation_h
54 
56 
57 #include <pcl/Defs.h>
58 #include <pcl/Diagnostics.h>
59 
61 #include <pcl/Math.h>
62 #include <pcl/Utility.h>
63 #include <pcl/Vector.h>
64 
65 namespace pcl
66 {
67 
68 // ----------------------------------------------------------------------------
69 
86 class PCL_CLASS CubicFilter
87 {
88 public:
89 
94  CubicFilter( double B, double C )
95  : m_B( B )
96  , m_C( C )
97  {
98  m_k31 = ( 12 - 9*m_B - 6*m_C)/6;
99  m_k21 = (-18 + 12*m_B + 6*m_C)/6;
100  m_k01 = ( 6 - 2*m_B )/6;
101  m_k32 = ( -m_B - 6*m_C)/6;
102  m_k22 = ( 6*m_B + 30*m_C)/6;
103  m_k12 = ( -12*m_B - 48*m_C)/6;
104  m_k02 = ( 8*m_B + 24*m_C)/6;
105  }
106 
110  CubicFilter( const CubicFilter& ) = default;
111 
115  virtual ~CubicFilter()
116  {
117  }
118 
122  CubicFilter& operator =( const CubicFilter& ) = default;
123 
134  PCL_HOT_FUNCTION
135  double operator()( double x ) const noexcept
136  {
137  if ( x < 0 )
138  x = -x;
139  return (x < 1) ? m_k01 + x*x*(m_k21 + x*m_k31) :
140  m_k02 + x*(m_k12 + x*(m_k22 + x*m_k32));
141  }
142 
151  double Width() const noexcept
152  {
153  return 2.0;
154  }
155 
159  virtual String Description() const
160  {
161  return String().Format( "Cubic filter (B=%.6f, C=%.6f)", m_B, m_C );
162  }
163 
167  virtual CubicFilter* Clone() const
168  {
169  return new CubicFilter( *this );
170  }
171 
172 protected:
173 
174  double m_k31, m_k21, m_k01, m_k32, m_k22, m_k12, m_k02;
175 
176 private:
177 
178  double m_B, m_C; // for reference only; not used in calculations
179 };
180 
191 {
192 public:
193 
198  : CubicFilter( 1.0/3, 1.0/3 )
199  {
200  }
201 
206 
211  {
212  }
213 
216  String Description() const override
217  {
218  return "Mitchell-Netravali cubic filter";
219  }
220 
223  CubicFilter* Clone() const override
224  {
225  return new MitchellNetravaliCubicFilter( *this );
226  }
227 };
228 
238 class PCL_CLASS CatmullRomSplineFilter : public CubicFilter
239 {
240 public:
241 
246  : CubicFilter( 0, 0.5 )
247  {
248  }
249 
254 
259  {
260  }
261 
264  String Description() const override
265  {
266  return "Catmull-Rom spline filter";
267  }
268 
271  CubicFilter* Clone() const override
272  {
273  return new CatmullRomSplineFilter( *this );
274  }
275 };
276 
286 class PCL_CLASS CubicBSplineFilter : public CubicFilter
287 {
288 public:
289 
294  : CubicFilter( 1, 0 )
295  {
296  }
297 
302 
307  {
308  }
309 
312  String Description() const override
313  {
314  return "Cubic B-spline filter";
315  }
316 
319  CubicFilter* Clone() const override
320  {
321  return new CubicBSplineFilter( *this );
322  }
323 };
324 
325 // ----------------------------------------------------------------------------
326 
327 #define m_width this->m_width
328 #define m_height this->m_height
329 #define m_fillBorder this->m_fillBorder
330 #define m_fillValue this->m_fillValue
331 #define m_data this->m_data
332 
333 // ----------------------------------------------------------------------------
334 
350 template <typename T>
352 {
353 public:
354 
366  BicubicFilterInterpolation( int rx, int ry, const CubicFilter& filter )
367  : m_rx( Max( 1, rx ) )
368  , m_ry( Max( 1, ry ) )
369  , m_filter( filter )
370  {
371  PCL_PRECONDITION( rx >= 1 )
372  PCL_PRECONDITION( ry >= 1 )
373  Initialize();
374  }
375 
380 
385  {
386  }
387 
390  void Initialize( const T* data, int dataWidth, int dataHeight ) override
391  {
392  BidimensionalInterpolation<T>::Initialize( data, dataWidth, dataHeight );
393 
394  if ( m_rx > m_width || m_ry > m_height )
395  {
396  m_rx = Min( m_rx, m_width );
397  m_ry = Min( m_ry, m_height );
398  Initialize();
399  }
400  }
401 
412  PCL_HOT_FUNCTION
413  double operator()( double x, double y ) const override
414  {
415  PCL_PRECONDITION( m_data != nullptr )
416  PCL_PRECONDITION( m_width > 0 && m_height > 0 )
417  PCL_PRECONDITION( x >= 0 && x < m_width )
418  PCL_PRECONDITION( y >= 0 && y < m_height )
419  PCL_CHECK( m_rx >= 1 && m_rx <= m_width )
420  PCL_CHECK( m_ry >= 1 && m_ry <= m_height )
421  PCL_CHECK( !m_k.IsEmpty() )
422 
423  int x0 = Range( RoundInt( x ), 0, m_width-1 );
424  int y0 = Range( RoundInt( y ), 0, m_height-1 );
425 
426  double dx = x - x0;
427  double dy = y - y0;
428 
429  int64 d = int64( y0 - m_ry )*m_width + x0 - m_rx;
430 
431  if ( m_rx >= m_ry )
432  {
433  double wx = 0;
434  for ( int xi = -m_rx, c = 0; xi <= m_rx; ++xi, ++c )
435  wx += (m_k[c] = m_filter( m_sx*(xi - dx) ));
436 
437  double sy = 0, wy = 0;
438  for ( int yi = -m_ry; yi <= m_ry; ++yi )
439  {
440  double ky = m_filter( m_sy*(yi - dy) );
441  wy += ky;
442 
443  const T* p;
444  int y = y0 + yi;
445 
446  if ( y < 0 )
447  {
448  d += m_width;
449 
450  if ( m_fillBorder )
451  {
452  sy += m_fillValue * ky;
453  continue;
454  }
455 
456  p = m_data + (d - 2*int64( m_width )*(y + 1));
457  }
458  else if ( y >= m_height )
459  {
460  if ( m_fillBorder )
461  {
462  sy += m_fillValue * ky;
463  continue;
464  }
465 
466  p = m_data + (d -= m_width);
467  }
468  else
469  {
470  p = m_data + d;
471  d += m_width;
472  }
473 
474  int x = x0 - m_rx;
475  int x2 = x0 + m_rx;
476  int x1 = Min( x2, m_width-1 );
477  double sx = 0;
478  const double* kx = m_k.Begin();
479 
480  while ( x < 0 )
481  {
482  ++x;
483  ++p;
484  sx += (m_fillBorder ? m_fillValue : double( *(p - (x + x)) )) * *kx++;
485  }
486 
487  for ( ;; )
488  {
489  sx += *p++ * *kx++;
490  if ( ++x > x1 )
491  break;
492  }
493 
494  while ( x <= x2 )
495  {
496  sx += (m_fillBorder ? m_fillValue : double( *--p )) * *kx++;
497  ++x;
498  }
499 
500  sy += sx/wx * ky;
501  }
502 
503  return sy/wy;
504  }
505  else
506  {
507  double wy = 0;
508  for ( int yi = -m_ry, r = 0; yi <= m_ry; ++yi, ++r )
509  wy += (m_k[r] = m_filter( m_sy*(yi - dy) ));
510 
511  double sx = 0, wx = 0;
512  for ( int xi = -m_rx; xi <= m_rx; ++xi )
513  {
514  double kx = m_filter( m_sx*(xi - dx) );
515  wx += kx;
516 
517  const T* p;
518  int x = x0 + xi;
519 
520  if ( x < 0 )
521  {
522  ++d;
523 
524  if ( m_fillBorder )
525  {
526  sx += m_fillValue * kx;
527  continue;
528  }
529 
530  p = m_data + (d - 2*(x + 1));
531  }
532  else if ( x >= m_width )
533  {
534  if ( m_fillBorder )
535  {
536  sx += m_fillValue * kx;
537  continue;
538  }
539 
540  p = m_data + --d;
541  }
542  else
543  p = m_data + d++;
544 
545  int y = y0 - m_ry;
546  int y2 = y0 + m_ry;
547  int y1 = Min( y2, m_height-1 );
548  double sy = 0;
549  const double* ky = m_k.Begin();
550 
551  while ( y < 0 )
552  {
553  ++y;
554  p += m_width;
555  sy += (m_fillBorder ? m_fillValue : double( *(p - (y + y)*m_width) )) * *ky++;
556  }
557 
558  for ( ;; )
559  {
560  sy += *p * *ky++;
561  p += m_width;
562  if ( ++y > y1 )
563  break;
564  }
565 
566  while ( y <= y2 )
567  {
568  sy += (m_fillBorder ? m_fillValue : double( *(p -= m_width) )) * *ky++;
569  ++y;
570  }
571 
572  sx += sy/wy * kx;
573  }
574 
575  return sx/wx;
576  }
577  }
578 
582  int HorizontalRadius() const noexcept
583  {
584  return m_rx;
585  }
586 
590  int VerticalRadius() const noexcept
591  {
592  return m_ry;
593  }
594 
603  void SetRadii( int rx, int ry )
604  {
605  if ( rx != m_rx || ry != m_ry )
606  {
607  m_rx = Max( 1, rx );
608  m_ry = Max( 1, ry );
609 
610  if ( this->data != nullptr )
611  {
612  if ( m_rx > m_width )
613  m_rx = m_width;
614  if ( m_ry > m_height )
615  m_ry = m_height;
616  }
617 
618  Initialize();
619  }
620  }
621 
626  const CubicFilter& Filter() const noexcept
627  {
628  return m_filter;
629  }
630 
634  void SetFilter( const CubicFilter& filter )
635  {
636  m_filter = filter;
637  Initialize();
638  }
639 
640 protected:
641 
642  int m_rx, m_ry; // horizontal and vertical radii
643  double m_sx, m_sy; // filter scaling ratios
644  mutable DVector m_k; // workspace for interpolated values
645  CubicFilter m_filter;
646 
647  void Initialize()
648  {
649  m_sx = m_filter.Width()/(m_rx + 0.5);
650  m_sy = m_filter.Width()/(m_ry + 0.5);
651  m_k = DVector( 1 + (Max( m_rx, m_ry ) << 1) );
652  }
653 };
654 
655 // ----------------------------------------------------------------------------
656 
657 #undef m_width
658 #undef m_height
659 #undef m_fillBorder
660 #undef m_fillValue
661 #undef m_data
662 
663 // ----------------------------------------------------------------------------
664 
665 } // pcl
666 
667 #endif // __PCL_BicubicFilterInterpolation_h
668 
669 // ----------------------------------------------------------------------------
670 // EOF pcl/BicubicFilterInterpolation.h - Released 2024-06-18T15:48:54Z
Bicubic filter interpolation algorithms.
const CubicFilter & Filter() const noexcept
void Initialize(const T *data, int dataWidth, int dataHeight) override
BicubicFilterInterpolation(const BicubicFilterInterpolation &)=default
PCL_HOT_FUNCTION double operator()(double x, double y) const override
BicubicFilterInterpolation(int rx, int ry, const CubicFilter &filter)
void SetFilter(const CubicFilter &filter)
A generic interface to two-dimensional interpolation algorithms.
virtual void Initialize(const T *data, int width, int height)
CubicFilter * Clone() const override
CatmullRomSplineFilter(const CatmullRomSplineFilter &)=default
CubicBSplineFilter(const CubicBSplineFilter &)=default
CubicFilter * Clone() const override
String Description() const override
Mitchell-Netravali parameterized cubic filters.
virtual CubicFilter * Clone() const
CubicFilter(double B, double C)
PCL_HOT_FUNCTION double operator()(double x) const noexcept
CubicFilter(const CubicFilter &)=default
virtual String Description() const
double Width() const noexcept
Generic vector of arbitrary length.
Definition: Vector.h:107
Mitchell-Netravali cubic filter with B=C=1/3.
CubicFilter * Clone() const override
MitchellNetravaliCubicFilter(const MitchellNetravaliCubicFilter &)=default
Unicode (UTF-16) string.
Definition: String.h:8113
String & Format(const_c_string8 fmt,...)
Definition: String.h:11081
int RoundInt(T x) noexcept
Definition: Math.h:1503
signed long long int64
Definition: Defs.h:676
constexpr const T & Min(const T &a, const T &b) noexcept
Definition: Utility.h:90
constexpr const T & Range(const T &x, const T &a, const T &b) noexcept
Definition: Utility.h:190
constexpr const T & Max(const T &a, const T &b) noexcept
Definition: Utility.h:119
PCL root namespace.
Definition: AbstractImage.h:77