PCL
BicubicFilterInterpolation.h
Go to the documentation of this file.
1 // ____ ______ __
2 // / __ \ / ____// /
3 // / /_/ // / / /
4 // / ____// /___ / /___ PixInsight Class Library
5 // /_/ \____//_____/ PCL 2.1.19
6 // ----------------------------------------------------------------------------
7 // pcl/BicubicFilterInterpolation.h - Released 2019-11-07T10:59:34Z
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-2019 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 (http://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 ) : m_B( B ), m_C( C )
95  {
96  m_k31 = ( 12 - 9*m_B - 6*m_C)/6;
97  m_k21 = (-18 + 12*m_B + 6*m_C)/6;
98  m_k01 = ( 6 - 2*m_B )/6;
99  m_k32 = ( -m_B - 6*m_C)/6;
100  m_k22 = ( 6*m_B + 30*m_C)/6;
101  m_k12 = ( -12*m_B - 48*m_C)/6;
102  m_k02 = ( 8*m_B + 24*m_C)/6;
103  }
104 
108  CubicFilter( const CubicFilter& ) = default;
109 
113  virtual ~CubicFilter()
114  {
115  }
116 
120  CubicFilter& operator =( const CubicFilter& ) = default;
121 
132  PCL_HOT_FUNCTION
133  double operator()( double x ) const
134  {
135  if ( x < 0 )
136  x = -x;
137  return (x < 1) ? m_k01 + x*x*(m_k21 + x*m_k31) :
138  m_k02 + x*(m_k12 + x*(m_k22 + x*m_k32));
139  }
140 
149  double Width() const
150  {
151  return 2.0;
152  }
153 
157  virtual String Description() const
158  {
159  return String().Format( "Cubic filter (B=%.6f, C=%.6f)", m_B, m_C );
160  }
161 
165  virtual CubicFilter* Clone() const
166  {
167  return new CubicFilter( *this );
168  }
169 
170 protected:
171 
172  double m_k31, m_k21, m_k01, m_k32, m_k22, m_k12, m_k02;
173 
174 private:
175 
176  double m_B, m_C; // for reference only; not used in calculations
177 };
178 
189 {
190 public:
191 
196  {
197  }
198 
203 
208  {
209  }
210 
213  String Description() const override
214  {
215  return "Mitchell-Netravali cubic filter";
216  }
217 
220  CubicFilter* Clone() const override
221  {
222  return new MitchellNetravaliCubicFilter( *this );
223  }
224 };
225 
235 class PCL_CLASS CatmullRomSplineFilter : public CubicFilter
236 {
237 public:
238 
243  {
244  }
245 
250 
255  {
256  }
257 
260  String Description() const override
261  {
262  return "Catmull-Rom spline filter";
263  }
264 
267  CubicFilter* Clone() const override
268  {
269  return new CatmullRomSplineFilter( *this );
270  }
271 };
272 
282 class PCL_CLASS CubicBSplineFilter : public CubicFilter
283 {
284 public:
285 
290  {
291  }
292 
296  CubicBSplineFilter( const CubicBSplineFilter& ) = default;
297 
302  {
303  }
304 
307  String Description() const override
308  {
309  return "Cubic B-spline filter";
310  }
311 
314  CubicFilter* Clone() const override
315  {
316  return new CubicBSplineFilter( *this );
317  }
318 };
319 
320 // ----------------------------------------------------------------------------
321 
322 #define m_width this->m_width
323 #define m_height this->m_height
324 #define m_fillBorder this->m_fillBorder
325 #define m_fillValue this->m_fillValue
326 #define m_data this->m_data
327 
328 // ----------------------------------------------------------------------------
329 
345 template <typename T>
347 {
348 public:
349 
361  BicubicFilterInterpolation( int rx, int ry, const CubicFilter& filter ) :
362  m_rx( Max( 1, rx ) ),
363  m_ry( Max( 1, ry ) ),
364  m_filter( filter )
365  {
366  PCL_PRECONDITION( rx >= 1 )
367  PCL_PRECONDITION( ry >= 1 )
368  Initialize();
369  }
370 
375 
380  {
381  }
382 
385  void Initialize( const T* data, int dataWidth, int dataHeight ) override
386  {
387  BidimensionalInterpolation<T>::Initialize( data, dataWidth, dataHeight );
388 
389  if ( m_rx > m_width || m_ry > m_height )
390  {
391  m_rx = Min( m_rx, m_width );
392  m_ry = Min( m_ry, m_height );
393  Initialize();
394  }
395  }
396 
407  PCL_HOT_FUNCTION
408  double operator()( double x, double y ) const override
409  {
410  PCL_PRECONDITION( m_data != nullptr )
411  PCL_PRECONDITION( m_width > 0 && m_height > 0 )
412  PCL_PRECONDITION( x >= 0 && x < m_width )
413  PCL_PRECONDITION( y >= 0 && y < m_height )
414  PCL_CHECK( m_rx >= 1 && m_rx <= m_width )
415  PCL_CHECK( m_ry >= 1 && m_ry <= m_height )
416  PCL_CHECK( !m_k.IsEmpty() )
417 
418  int x0 = Range( RoundInt( x ), 0, m_width-1 );
419  int y0 = Range( RoundInt( y ), 0, m_height-1 );
420 
421  double dx = x - x0;
422  double dy = y - y0;
423 
424  int64 d = int64( y0 - m_ry )*m_width + x0 - m_rx;
425 
426  if ( m_rx >= m_ry )
427  {
428  double wx = 0;
429  for ( int xi = -m_rx, c = 0; xi <= m_rx; ++xi, ++c )
430  wx += (m_k[c] = m_filter( m_sx*(xi - dx) ));
431 
432  double sy = 0, wy = 0;
433  for ( int yi = -m_ry; yi <= m_ry; ++yi )
434  {
435  double ky = m_filter( m_sy*(yi - dy) );
436  wy += ky;
437 
438  const T* p;
439  int y = y0 + yi;
440 
441  if ( y < 0 )
442  {
443  d += m_width;
444 
445  if ( m_fillBorder )
446  {
447  sy += m_fillValue * ky;
448  continue;
449  }
450 
451  p = m_data + (d - 2*int64( m_width )*(y + 1));
452  }
453  else if ( y >= m_height )
454  {
455  if ( m_fillBorder )
456  {
457  sy += m_fillValue * ky;
458  continue;
459  }
460 
461  p = m_data + (d -= m_width);
462  }
463  else
464  {
465  p = m_data + d;
466  d += m_width;
467  }
468 
469  int x = x0 - m_rx;
470  int x2 = x0 + m_rx;
471  int x1 = Min( x2, m_width-1 );
472  double sx = 0;
473  const double* kx = m_k.Begin();
474 
475  while ( x < 0 )
476  {
477  ++x;
478  ++p;
479  sx += (m_fillBorder ? m_fillValue : double( *(p - (x + x)) )) * *kx++;
480  }
481 
482  for ( ;; )
483  {
484  sx += *p++ * *kx++;
485  if ( ++x > x1 )
486  break;
487  }
488 
489  while ( x <= x2 )
490  {
491  sx += (m_fillBorder ? m_fillValue : double( *--p )) * *kx++;
492  ++x;
493  }
494 
495  sy += sx/wx * ky;
496  }
497 
498  return sy/wy;
499  }
500  else
501  {
502  double wy = 0;
503  for ( int yi = -m_ry, r = 0; yi <= m_ry; ++yi, ++r )
504  wy += (m_k[r] = m_filter( m_sy*(yi - dy) ));
505 
506  double sx = 0, wx = 0;
507  for ( int xi = -m_rx; xi <= m_rx; ++xi )
508  {
509  double kx = m_filter( m_sx*(xi - dx) );
510  wx += kx;
511 
512  const T* p;
513  int x = x0 + xi;
514 
515  if ( x < 0 )
516  {
517  ++d;
518 
519  if ( m_fillBorder )
520  {
521  sx += m_fillValue * kx;
522  continue;
523  }
524 
525  p = m_data + (d - 2*(x + 1));
526  }
527  else if ( x >= m_width )
528  {
529  if ( m_fillBorder )
530  {
531  sx += m_fillValue * kx;
532  continue;
533  }
534 
535  p = m_data + --d;
536  }
537  else
538  p = m_data + d++;
539 
540  int y = y0 - m_ry;
541  int y2 = y0 + m_ry;
542  int y1 = Min( y2, m_height-1 );
543  double sy = 0;
544  const double* ky = m_k.Begin();
545 
546  while ( y < 0 )
547  {
548  ++y;
549  p += m_width;
550  sy += (m_fillBorder ? m_fillValue : double( *(p - (y + y)*m_width) )) * *ky++;
551  }
552 
553  for ( ;; )
554  {
555  sy += *p * *ky++;
556  p += m_width;
557  if ( ++y > y1 )
558  break;
559  }
560 
561  while ( y <= y2 )
562  {
563  sy += (m_fillBorder ? m_fillValue : double( *(p -= m_width) )) * *ky++;
564  ++y;
565  }
566 
567  sx += sy/wy * kx;
568  }
569 
570  return sx/wx;
571  }
572  }
573 
577  int HorizontalRadius() const
578  {
579  return m_rx;
580  }
581 
585  int VerticalRadius() const
586  {
587  return m_ry;
588  }
589 
598  void SetRadii( int rx, int ry )
599  {
600  if ( rx != m_rx || ry != m_ry )
601  {
602  m_rx = Max( 1, rx );
603  m_ry = Max( 1, ry );
604 
605  if ( this->data != nullptr )
606  {
607  if ( m_rx > m_width )
608  m_rx = m_width;
609  if ( m_ry > m_height )
610  m_ry = m_height;
611  }
612 
613  Initialize();
614  }
615  }
616 
621  const CubicFilter& Filter() const
622  {
623  return m_filter;
624  }
625 
629  void SetFilter( const CubicFilter& filter )
630  {
631  m_filter = filter;
632  Initialize();
633  }
634 
635 protected:
636 
637  int m_rx, m_ry; // horizontal and vertical radii
638  double m_sx, m_sy; // filter scaling ratios
639  mutable DVector m_k; // workspace for interpolated values
640  CubicFilter m_filter;
641 
642  void Initialize()
643  {
644  m_sx = m_filter.Width()/(m_rx + 0.5);
645  m_sy = m_filter.Width()/(m_ry + 0.5);
646  m_k = DVector( 1 + (Max( m_rx, m_ry ) << 1) );
647  }
648 };
649 
650 // ----------------------------------------------------------------------------
651 
652 #undef m_width
653 #undef m_height
654 #undef m_fillBorder
655 #undef m_fillValue
656 #undef m_data
657 
658 // ----------------------------------------------------------------------------
659 
660 } // pcl
661 
662 #endif // __PCL_BicubicFilterInterpolation_h
663 
664 // ----------------------------------------------------------------------------
665 // EOF pcl/BicubicFilterInterpolation.h - Released 2019-11-07T10:59:34Z
Bicubic filter interpolation algorithms.
int RoundInt(T x)
Definition: Math.h:1397
PCL root namespace.
Definition: AbstractImage.h:76
A generic interface to two-dimensional interpolation algorithms.
void SetFilter(const CubicFilter &filter)
CubicFilter * Clone() const override
PCL_HOT_FUNCTION double operator()(double x, double y) const override
void Initialize(const T *data, int dataWidth, int dataHeight) override
constexpr const T & Range(const T &x, const T &a, const T &b)
Definition: Utility.h:190
String & Format(const_c_string8 fmt,...)
Definition: String.h:10865
virtual String Description() const
Unicode (UTF-16) string.
Definition: String.h:7911
String Description() const override
constexpr const T & Max(const T &a, const T &b)
Definition: Utility.h:119
64-bit floating point real vector.
constexpr const T & Min(const T &a, const T &b)
Definition: Utility.h:90
virtual CubicFilter * Clone() const
CubicFilter(double B, double C)
CubicFilter * Clone() const override
virtual void Initialize(const T *data, int width, int height)
PCL_HOT_FUNCTION double operator()(double x) const
BicubicFilterInterpolation(int rx, int ry, const CubicFilter &filter)
Mitchell-Netravali parameterized cubic filters.
CubicFilter * Clone() const override
signed long long int64
Definition: Defs.h:610
Mitchell-Netravali cubic filter with B=C=1/3.