PCL
SurfaceSpline.h
Go to the documentation of this file.
1 // ____ ______ __
2 // / __ \ / ____// /
3 // / /_/ // / / /
4 // / ____// /___ / /___ PixInsight Class Library
5 // /_/ \____//_____/ PCL 2.7.0
6 // ----------------------------------------------------------------------------
7 // pcl/SurfaceSpline.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_SurfaceSpline_h
53 #define __PCL_SurfaceSpline_h
54 
56 
57 #include <pcl/Defs.h>
58 #include <pcl/Diagnostics.h>
59 
60 #include <pcl/AbstractImage.h>
61 #include <pcl/Array.h>
62 #include <pcl/AutoPointer.h>
63 #include <pcl/Homography.h>
64 #include <pcl/ParallelProcess.h>
65 #include <pcl/Point.h>
66 #include <pcl/QuadTree.h>
67 #include <pcl/SurfaceSimplifier.h>
68 #include <pcl/Thread.h>
69 #include <pcl/Vector.h>
70 
71 #ifndef __PCL_BUILDING_PIXINSIGHT_APPLICATION
72 # include <pcl/Console.h>
73 # include <pcl/StandardStatus.h>
74 #endif
75 
76 namespace pcl
77 {
78 
79 // ----------------------------------------------------------------------------
80 
81 /*
82  * The maximum allowed number of points in a point surface spline by default.
83  */
84 #define __PCL_PSSPLINE_DEFAULT_MAX_POINTS 2100
85 
86 /*
87  * Default quadtree bucket capacity for recursive surface subspline generation.
88  */
89 #define __PCL_RSSPLINE_DEFAULT_TREE_BUCKET_CAPACITY 64
90 
91 /*
92  * Default maximum spline length for a non-recursive quadtree node in a
93  * recursive surface spline.
94  */
95 #define __PCL_RSSPLINE_DEFAULT_SPLINE_MAX_LENGTH 1600
96 
97 /*
98  * Whether to allow extrapolation outside the interpolation region for
99  * recursive surface splines. Extrapolation is disabled by default because
100  * recursively defined subsplines are slightly more prone to oscillation than
101  * normal surface splines.
102  */
103 #define __PCL_RSSPLINE_DEFAULT_ALLOW_EXTRAPOLATION false
104 
105 // ----------------------------------------------------------------------------
106 
129 namespace RadialBasisFunction
130 {
131  enum value_type
132  {
133  Unknown = -1,
134  VariableOrder = 0, // phi(r) = (r^2)^(m-1) * Ln( r^2 )
135  ThinPlateSpline, // phi(r) = r^2 * Ln( r )
136  Gaussian, // phi(r) = Exp( -(eps*r)^2 )
137  Multiquadric, // phi(r) = Sqrt( 1 + (eps*r)^2 )
138  InverseMultiquadric, // phi(r) = 1/Sqrt( 1 + (eps*r)^2 )
139  InverseQuadratic, // phi(r) = 1/( 1 + (eps*r)^2 )
140  __number_of_items__,
141  Default = ThinPlateSpline
142  };
143 
144  inline static bool Validate( int rbf )
145  {
146  return rbf >= 0 && rbf < __number_of_items__;
147  }
148 }
149 
150 // ----------------------------------------------------------------------------
151 
156 class PCL_CLASS SurfaceSplineBase
157 {
158 public:
159 
164  using rbf_type = RadialBasisFunction::value_type;
165 
166 protected:
167 
171  SurfaceSplineBase() = default;
172 
177 
182 
187  {
188  }
189 
193  SurfaceSplineBase& operator =( const SurfaceSplineBase& ) = default;
194 
198  SurfaceSplineBase& operator =( SurfaceSplineBase&& ) = default;
199 
203  static void Generate( float* __restrict__ c,
204  rbf_type, double e2, bool polynomial,
205  const float* __restrict__ x, const float* __restrict__ y, const float* __restrict__ z,
206  int n, int m, float r, const float* __restrict__ w );
207 
211  static void Generate( double* __restrict__ c,
212  rbf_type, double e2, bool polynomial,
213  const double* __restrict__ x, const double* __restrict__ y, const double* __restrict__ z,
214  int n, int m, float r, const float* __restrict__ w );
215 };
216 
217 // ----------------------------------------------------------------------------
218 
266 template <typename T>
267 class PCL_CLASS SurfaceSpline : private SurfaceSplineBase
268 {
269 public:
270 
275  using scalar = T;
276 
282 
287  using weight = float;
288 
293 
298  using rbf_type = SurfaceSplineBase::rbf_type;
299 
304  SurfaceSpline() = default;
305 
309  SurfaceSpline( const SurfaceSpline& ) = default;
310 
314  SurfaceSpline( SurfaceSpline&& ) = default;
315 
319  ~SurfaceSpline() override
320  {
321  }
322 
326  SurfaceSpline& operator =( const SurfaceSpline& ) = default;
327 
331  SurfaceSpline& operator =( SurfaceSpline&& ) = default;
332 
337  bool IsValid() const
338  {
339  return m_x.Length() >= 3 && m_x.Length() == m_y.Length();
340  }
341 
345  int Length() const
346  {
347  return m_x.Length();
348  }
349 
355  vector X() const
356  {
357  vector x( m_x.Length() );
358  if ( IsValid() )
359  for ( int i = 0; i < m_x.Length(); ++i )
360  x[i] = m_x0 + m_x[i]/m_r0;
361  return x;
362  }
363 
369  vector Y() const
370  {
371  vector y( m_y.Length() );
372  if ( IsValid() )
373  for ( int i = 0; i < m_y.Length(); ++i )
374  y[i] = m_y0 + m_y[i]/m_r0;
375  return y;
376  }
377 
383  const vector& Z() const
384  {
385  return m_z;
386  }
387 
394  const weight_vector& Weights() const
395  {
396  return m_weights;
397  }
398 
404  rbf_type RBFType() const
405  {
406  return m_rbf;
407  }
408 
417  void SetRBFType( rbf_type rbf )
418  {
419  Clear();
420  m_rbf = rbf;
421  }
422 
431  double ShapeParameter() const
432  {
433  return Sqrt( m_eps2 );
434  }
435 
451  void SetShapeParameter( double eps )
452  {
453  Clear();
454  m_eps = Abs( eps );
455  m_eps2 = m_eps*m_eps;
456  }
457 
461  int Order() const
462  {
463  return m_order;
464  }
465 
487  void SetOrder( int order )
488  {
489  PCL_PRECONDITION( order > 1 )
490  Clear();
491  m_order = pcl::Max( 2, order );
492  }
493 
502  bool IsPolynomialEnabled() const
503  {
504  return m_havePolynomial;
505  }
506 
514  void EnablePolynomial( bool enable = true )
515  {
516  Clear();
517  m_havePolynomial = enable;
518  }
519 
527  void DisablePolynomial( bool disable = true )
528  {
529  EnablePolynomial( !disable );
530  }
531 
536  float Smoothing() const
537  {
538  return m_smoothing;
539  }
540 
558  void SetSmoothing( float s )
559  {
560  PCL_PRECONDITION( s >= 0 )
561  Clear();
562  m_smoothing = pcl::Max( 0.0F, s );
563  }
564 
615  void Initialize( const scalar* x, const scalar* y, const scalar* z, int n, const weight* w = nullptr )
616  {
617  PCL_PRECONDITION( x != nullptr && y != nullptr && z != nullptr )
618  PCL_PRECONDITION( n > 2 )
619 
620  if ( x == nullptr || y == nullptr || z == nullptr )
621  throw Error( "SurfaceSpline::Initialize(): Null vector pointer(s)." );
622 
623  if ( n < 3 )
624  throw Error( "SurfaceSpline::Initialize(): At least three input nodes must be specified." );
625 
626  Clear();
627 
628  if ( m_smoothing <= 0 )
629  w = nullptr;
630 
631  try
632  {
633  /*
634  * Find mean coordinates.
635  */
636  m_x0 = m_y0 = 0;
637  for ( int i = 0; i < n; ++i )
638  {
639  m_x0 += x[i];
640  m_y0 += y[i];
641  }
642  m_x0 /= n;
643  m_y0 /= n;
644 
645  /*
646  * Find radius of largest containing circle.
647  */
648  m_r0 = 0;
649  for ( int i = 0; i < n; ++i )
650  {
651  double dx = x[i] - m_x0;
652  double dy = y[i] - m_y0;
653  double r = Sqrt( dx*dx + dy*dy );
654  if ( r > m_r0 )
655  m_r0 = r;
656  }
657  if ( 1 + m_r0 == 1 )
658  throw Error( "SurfaceSpline::Initialize(): Empty or insignificant interpolation space." );
659  m_r0 = 1/m_r0;
660 
661  /*
662  * If requested, compute an optimal shape parameter.
663  */
664  if ( m_eps == 0 )
665  {
666  Array<double> R2;
667  for ( int i = 0; i < n; ++i )
668  for ( int j = i; ++j < n; )
669  {
670  double dx = x[i] - x[j];
671  double dy = y[i] - y[j];
672  R2 << dx*dx + dy*dy;
673  }
674  m_eps2 = 1/(m_r0*4*Sqrt( Median( R2.Begin(), R2.End() ) ));
675  }
676  else
677  m_eps2 = m_r0*m_eps;
678  m_eps2 *= m_eps2;
679 
680  /*
681  * Build point list with normalized node coordinates.
682  */
683  Array<NodeData> P;
684  for ( int i = 0; i < n; ++i )
685  P << NodeData( m_r0*(x[i] - m_x0),
686  m_r0*(y[i] - m_y0),
687  z[i],
688  (w != nullptr && w[i] > 0) ? w[i] : 1.0F );
689 
690  /*
691  * Find duplicate input nodes. Two nodes are considered equal if their
692  * (normalized) coordinates don't differ more than the machine epsilon
693  * for the floating point type T.
694  */
695  P.Sort();
696  Array<int> remove;
697  for ( int i = 0, j = 1; j < n; ++i, ++j )
698  if ( Abs( P[i].x - P[j].x ) <= std::numeric_limits<scalar>::epsilon() )
699  if ( Abs( P[i].y - P[j].y ) <= std::numeric_limits<scalar>::epsilon() )
700  remove << i;
701 
702  /*
703  * Build working vectors, excluding duplicate input nodes.
704  */
705  int N = n - int( remove.Length() );
706  if ( N < 3 )
707  throw Error( "SurfaceSpline::Initialize(): Less than three input nodes left after sanitization." );
708  m_x = vector( N );
709  m_y = vector( N );
710  m_z = vector( N );
711  if ( w != nullptr )
712  m_weights = weight_vector( N );
713  int i = 0, k = 0;
714  for ( int j : remove )
715  {
716  for ( ; i < j; ++i, ++k )
717  {
718  m_x[k] = P[i].x;
719  m_y[k] = P[i].y;
720  m_z[k] = P[i].z;
721  if ( w != nullptr )
722  m_weights[k] = P[i].w;
723  }
724  ++i;
725  }
726  for ( ; i < n; ++i, ++k )
727  {
728  m_x[k] = P[i].x;
729  m_y[k] = P[i].y;
730  m_z[k] = P[i].z;
731  if ( w != nullptr )
732  m_weights[k] = P[i].w;
733  }
734 
735  m_spline = vector( scalar( 0 ), N + (m_havePolynomial ? ((m_order*(m_order + 1)) >> 1) : 0) );
736 
737  Generate( m_spline.Begin(), m_rbf, m_eps2, m_havePolynomial,
738  m_x.Begin(), m_y.Begin(), m_z.Begin(), N, m_order,
739  m_smoothing, m_weights.Begin() );
740  }
741  catch ( ... )
742  {
743  Clear();
744  throw;
745  }
746  }
747 
752  void Clear()
753  {
754  m_x.Clear();
755  m_y.Clear();
756  m_z.Clear();
757  m_weights.Clear();
758  m_spline.Clear();
759  }
760 
770  scalar operator ()( double x, double y ) const
771  {
772  PCL_PRECONDITION( !m_x.IsEmpty() && !m_y.IsEmpty() )
773  PCL_CHECK( m_order >= 2 )
774  PCL_CHECK( !m_spline.IsEmpty() )
775 
776  /*
777  * Normalized interpolation coordinates.
778  */
779  x = m_r0*(x - m_x0);
780  y = m_r0*(y - m_y0);
781 
782  /*
783  * Add polynomial part of the surface spline.
784  */
785  int n = m_x.Length();
786  double z = 0;
787  if ( m_havePolynomial )
788  {
789  z += m_spline[n];
790  switch ( m_order )
791  {
792  case 2:
793  z += m_spline[n+1]*x + m_spline[n+2]*y;
794  break;
795  case 3:
796  z += (m_spline[n+1] + m_spline[n+3]*x + m_spline[n+4]*y)*x + (m_spline[n+2] + m_spline[n+5]*y)*y;
797  break;
798  default:
799  for ( int i = 1, j = n+1, i1 = (m_order*(m_order + 1))>>1, ix = 0, iy = 0; i < i1; ++i, ++j )
800  if ( ix == 0 )
801  {
802  ix = iy + 1;
803  iy = 0;
804  z += m_spline[j] * PowI( x, ix );
805  }
806  else
807  {
808  --ix;
809  ++iy;
810  z += m_spline[j] * PowI( x, ix ) * PowI( y, iy );
811  }
812  break;
813  }
814  }
815 
816  /*
817  * Add radial basis functions.
818  */
819  switch ( m_rbf )
820  {
821  case RadialBasisFunction::VariableOrder:
822  for ( int i = 0; i < n; ++i )
823  {
824  double dx = m_x[i] - x;
825  double dy = m_y[i] - y;
826  double r2 = dx*dx + dy*dy;
827  if ( r2 != 0 )
828  {
829  double r2m1 = r2;
830  for ( int j = m_order; --j > 1; )
831  r2m1 *= r2;
832  z += m_spline[i] * r2m1 * pcl::Ln( r2 );
833  }
834  }
835  break;
836  case RadialBasisFunction::ThinPlateSpline:
837  for ( int i = 0; i < n; ++i )
838  {
839  double dx = m_x[i] - x;
840  double dy = m_y[i] - y;
841  double r2 = dx*dx + dy*dy;
842  if ( r2 != 0 )
843  z += m_spline[i] * r2 * pcl::Ln( Sqrt( r2 ) );
844  }
845  break;
846  case RadialBasisFunction::Gaussian:
847  for ( int i = 0; i < n; ++i )
848  {
849  double dx = m_x[i] - x;
850  double dy = m_y[i] - y;
851  z += m_spline[i] * pcl::Exp( -m_eps2 * (dx*dx + dy*dy) );
852  }
853  break;
854  case RadialBasisFunction::Multiquadric:
855  for ( int i = 0; i < n; ++i )
856  {
857  double dx = m_x[i] - x;
858  double dy = m_y[i] - y;
859  z += m_spline[i] * pcl::Sqrt( 1 + m_eps2 * (dx*dx + dy*dy) );
860  }
861  break;
862  case RadialBasisFunction::InverseMultiquadric:
863  for ( int i = 0; i < n; ++i )
864  {
865  double dx = m_x[i] - x;
866  double dy = m_y[i] - y;
867  z += m_spline[i] / pcl::Sqrt( 1 + m_eps2 * (dx*dx + dy*dy) );
868  }
869  break;
870  case RadialBasisFunction::InverseQuadratic:
871  for ( int i = 0; i < n; ++i )
872  {
873  double dx = m_x[i] - x;
874  double dy = m_y[i] - y;
875  z += m_spline[i] / (1 + m_eps2 * (dx*dx + dy*dy));
876  }
877  break;
878  default:
879  break;
880  }
881 
882  return scalar( z );
883  }
884 
890  template <typename Tp>
891  scalar operator ()( const GenericPoint<Tp>& p ) const
892  {
893  return operator ()( double( p.x ), double( p.y ) );
894  }
895 
896 protected:
897 
898  rbf_type m_rbf = RadialBasisFunction::Default;
899  bool m_havePolynomial = true; // use 1st order polynomial for stabilization
900  double m_eps = 0; // shape parameter, 0 -> find optimal value automatically
901  double m_eps2 = 0; // squared shape parameter
902  vector m_x; // vector of normalized X node coordinates
903  vector m_y; // vector of normalized Y node coordinates
904  vector m_z; // vector of function values - for reference only
905  double m_r0 = 1; // scaling factor for normalization of node coordinates
906  double m_x0 = 0; // zero offset for normalization of X node coordinates
907  double m_y0 = 0; // zero offset for normalization of Y node coordinates
908  int m_order = 2; // derivative order >= 2
909  float m_smoothing = 0; // <= 0 -> interpolating spline, > 0 -> smoothing factor of approximating spline
910  weight_vector m_weights; // optional node weights for approximating spline
911  vector m_spline; // coefficients of the 2-D surface spline, polynomial coeffs. at tail
912 
918  struct NodeData
919  {
920  scalar x, y, z;
921  weight w;
922 
923  NodeData( scalar a_x, scalar a_y, scalar a_z, weight a_w )
924  : x( a_x ), y( a_y ), z( a_z ), w( a_w )
925  {
926  }
927 
928  bool operator ==( const NodeData& p ) const
929  {
930  return x == p.x && y == p.y;
931  }
932 
933  bool operator <( const NodeData& p ) const
934  {
935  return (x != p.x) ? x < p.x : y < p.y;
936  }
937  };
938 
939  friend class DrizzleData;
940  friend class DrizzleDataDecoder;
941  friend class SplineWorldTransformation;
942 };
943 
944 // ----------------------------------------------------------------------------
945 
959 class PCL_CLASS PointSurfaceSpline
960 {
961 public:
962 
967 
972  using rbf_type = spline::rbf_type;
973 
978  PointSurfaceSpline() = default;
979 
984 
989 
997  template <class point_list1, class point_list2, class weight_vector = FVector>
998  PointSurfaceSpline( const point_list1& P1, const point_list2& P2,
999  float smoothness = 0, int order = 2,
1000  const weight_vector& W = weight_vector(),
1001  rbf_type rbf = RadialBasisFunction::Default,
1002  double eps = 0,
1003  bool polynomial = true )
1004  {
1005  Initialize( P1, P2, smoothness, W, order, rbf, eps, polynomial );
1006  }
1007 
1015  PointSurfaceSpline( const spline& Sx, const spline& Sy )
1016  {
1017  Initialize( Sx, Sy );
1018  }
1019 
1023  PointSurfaceSpline& operator =( const PointSurfaceSpline& other ) = default;
1024 
1028  PointSurfaceSpline& operator =( PointSurfaceSpline&& ) = default;
1029 
1107  template <class point_list1, class point_list2, class weight_vector = FVector>
1108  void Initialize( const point_list1& P1, const point_list2& P2,
1109  float smoothness = 0, const weight_vector& W = weight_vector(), int order = 2,
1110  rbf_type rbf = RadialBasisFunction::Default,
1111  double eps = 0,
1112  bool polynomial = true )
1113  {
1114  PCL_PRECONDITION( P1.Length() >= 3 )
1115  PCL_PRECONDITION( P1.Length() <= P2.Length() )
1116  PCL_PRECONDITION( order >= 2 )
1117  PCL_PRECONDITION( W.IsEmpty() || P1.Length() <= size_type( W.Length() ) )
1118 
1119  Clear();
1120 
1121  int N = int( P1.Length() );
1122  if ( N < 3 || P2.Length() < 3 )
1123  throw Error( "PointSurfaceSpline::Initialize(): At least three input nodes must be specified." );
1124  if ( int( P2.Length() ) < N || !W.IsEmpty() && int( W.Length() ) < N )
1125  throw Error( "PointSurfaceSpline::Initialize(): Insufficient data." );
1126 
1127  m_Sx.SetRBFType( rbf );
1128  m_Sx.SetOrder( order );
1129  m_Sx.SetShapeParameter( eps );
1130  m_Sx.EnablePolynomial( polynomial );
1131  m_Sx.SetSmoothing( smoothness );
1132 
1133  m_Sy.SetRBFType( rbf );
1134  m_Sy.SetOrder( order );
1135  m_Sy.SetShapeParameter( eps );
1136  m_Sy.EnablePolynomial( polynomial );
1137  m_Sy.SetSmoothing( smoothness );
1138 
1139  DVector X( N ), Y( N ), ZX( N ), ZY( N );
1140  double zxMin = std::numeric_limits<double>::max();
1141  double zxMax = -std::numeric_limits<double>::max();
1142  double zyMin = std::numeric_limits<double>::max();
1143  double zyMax = -std::numeric_limits<double>::max();
1144  for ( int i = 0; i < N; ++i )
1145  {
1146  X[i] = double( P1[i].x );
1147  Y[i] = double( P1[i].y );
1148 
1149  if ( m_incremental )
1150  {
1151  ZX[i] = X[i];
1152  ZY[i] = Y[i];
1153  m_H.Apply( ZX[i], ZY[i] );
1154  ZX[i] = double( P2[i].x ) - ZX[i];
1155  ZY[i] = double( P2[i].y ) - ZY[i];
1156  }
1157  else
1158  {
1159  ZX[i] = double( P2[i].x );
1160  ZY[i] = double( P2[i].y );
1161  }
1162 
1163  if ( ZX[i] < zxMin )
1164  zxMin = ZX[i];
1165  if ( ZX[i] > zxMax )
1166  zxMax = ZX[i];
1167  if ( ZY[i] < zyMin )
1168  zyMin = ZY[i];
1169  if ( ZY[i] > zyMax )
1170  zyMax = ZY[i];
1171  }
1172 
1173  if ( m_useSimplifiers )
1174  {
1175  SurfaceSimplifier SS;
1176  SS.EnableRejection();
1177  SS.SetRejectFraction( m_simplifierRejectFraction );
1179 
1180  // Simplified surface, X axis.
1181  DVector XXS, YXS, ZXS;
1182  m_epsX = (zxMax - zxMin)/100;
1183  double epsLow = 0, epsHigh = (zxMax - zxMin)/10;
1184  for ( int i = 0; i < 200; ++i )
1185  {
1186  SS.SetTolerance( m_epsX );
1187  SS.Simplify( XXS, YXS, ZXS, X, Y, ZX );
1188  if ( XXS.Length() > m_maxSplinePoints )
1189  {
1190  epsLow = m_epsX;
1191  if ( epsHigh - epsLow <= 4*std::numeric_limits<double>::epsilon() )
1192  epsHigh *= 2;
1193  }
1194  else
1195  {
1196  if ( epsHigh - epsLow <= 2*std::numeric_limits<double>::epsilon() )
1197  break;
1198  epsHigh = m_epsX;
1199  }
1200  m_epsX = (epsLow + epsHigh)/2;
1201  }
1202  if ( XXS.Length() > m_maxSplinePoints )
1203  {
1204  m_truncatedX = true;
1205  XXS = DVector( XXS.Begin(), m_maxSplinePoints );
1206  YXS = DVector( YXS.Begin(), m_maxSplinePoints );
1207  ZXS = DVector( ZXS.Begin(), m_maxSplinePoints );
1208  }
1209 
1210  // Simplified surface, Y axis.
1211  DVector XYS, YYS, ZYS;
1212  m_epsY = (zyMax - zyMin)/100;
1213  epsLow = 0, epsHigh = (zyMax - zyMin)/10;
1214  for ( int i = 0; i < 200; ++i )
1215  {
1216  SS.SetTolerance( m_epsY );
1217  SS.Simplify( XYS, YYS, ZYS, X, Y, ZY );
1218  if ( XYS.Length() > m_maxSplinePoints )
1219  {
1220  epsLow = m_epsY;
1221  if ( epsHigh - epsLow <= 4*std::numeric_limits<double>::epsilon() )
1222  epsHigh *= 2;
1223  }
1224  else
1225  {
1226  if ( epsHigh - epsLow <= 2*std::numeric_limits<double>::epsilon() )
1227  break;
1228  epsHigh = m_epsY;
1229  }
1230  m_epsY = (epsLow + epsHigh)/2;
1231  }
1232  if ( XYS.Length() > m_maxSplinePoints )
1233  {
1234  m_truncatedY = true;
1235  XYS = DVector( XYS.Begin(), m_maxSplinePoints );
1236  YYS = DVector( YYS.Begin(), m_maxSplinePoints );
1237  ZYS = DVector( ZYS.Begin(), m_maxSplinePoints );
1238  }
1239 
1240  SplineGenerationThread<weight_vector> Tx( m_Sx, XXS, YXS, ZXS, XXS.Length(), weight_vector() );
1241  SplineGenerationThread<weight_vector> Ty( m_Sy, XYS, YYS, ZYS, XYS.Length(), weight_vector() );
1242  Tx.Start( ThreadPriority::DefaultMax );
1243  Ty.Start( ThreadPriority::DefaultMax );
1244  Tx.Wait();
1245  Ty.Wait();
1246  Tx.ValidateOrThrow();
1247  Ty.ValidateOrThrow();
1248  }
1249  else
1250  {
1251  m_truncatedX = m_truncatedY = N > m_maxSplinePoints;
1252  SplineGenerationThread<weight_vector> Tx( m_Sx, X, Y, ZX, Min( N, m_maxSplinePoints ), W );
1253  SplineGenerationThread<weight_vector> Ty( m_Sy, X, Y, ZY, Min( N, m_maxSplinePoints ), W );
1254  Tx.Start( ThreadPriority::DefaultMax );
1255  Ty.Start( ThreadPriority::DefaultMax );
1256  Tx.Wait();
1257  Ty.Wait();
1258  Tx.ValidateOrThrow();
1259  Ty.ValidateOrThrow();
1260  }
1261  }
1262 
1279  void Initialize( const spline& Sx, const spline& Sy )
1280  {
1281  Clear();
1282  m_useSimplifiers = m_incremental = false;
1283  m_H = Homography();
1284  if ( Sx.IsValid() && Sy.IsValid() )
1285  {
1286  m_Sx = Sx;
1287  m_Sy = Sy;
1288  }
1289  }
1290 
1295  void Clear()
1296  {
1297  m_Sx.Clear();
1298  m_Sy.Clear();
1299  m_epsX = m_epsY = 0;
1300  m_truncatedX = m_truncatedY = false;
1301  }
1302 
1307  bool IsValid() const
1308  {
1309  return m_Sx.IsValid() && m_Sy.IsValid() && (!m_incremental || m_H.IsValid());
1310  }
1311 
1316  const spline& SplineX() const
1317  {
1318  return m_Sx;
1319  }
1320 
1325  const spline& SplineY() const
1326  {
1327  return m_Sy;
1328  }
1329 
1334  int MaxSplinePoints() const
1335  {
1336  return m_maxSplinePoints;
1337  }
1338 
1347  void SetMaxSplinePoints( int n )
1348  {
1349  PCL_PRECONDITION( n >= 3 )
1350  m_maxSplinePoints = Max( 3, n );
1351  }
1352 
1359  bool SimplifiersEnabled() const
1360  {
1361  return m_useSimplifiers;
1362  }
1363 
1371  void EnableSimplifiers( bool enable = true )
1372  {
1373  m_useSimplifiers = enable;
1374  }
1375 
1383  void DisableSimplifiers( bool disable = true )
1384  {
1385  EnableSimplifiers( !disable );
1386  }
1387 
1394  {
1395  return m_simplifierRejectFraction;
1396  }
1397 
1404  void SetSimplifierRejectFraction( float rejectFraction )
1405  {
1406  PCL_PRECONDITION( rejectFraction > 0 && rejectFraction < 1 )
1407  m_simplifierRejectFraction = Range( rejectFraction, 0.0F, 1.0F );
1408  }
1409 
1414  double ErrorX() const
1415  {
1416  return m_epsX;
1417  }
1418 
1423  double ErrorY() const
1424  {
1425  return m_epsY;
1426  }
1427 
1433  bool TruncatedX() const
1434  {
1435  return m_truncatedX;
1436  }
1437 
1443  bool TruncatedY() const
1444  {
1445  return m_truncatedY;
1446  }
1447 
1453  bool Truncated() const
1454  {
1455  return TruncatedX() || TruncatedY();
1456  }
1457 
1465  const Matrix& LinearFunction() const
1466  {
1467  return m_H;
1468  }
1469 
1477  void SetLinearFunction( const Matrix& H )
1478  {
1479  m_H = H;
1480  }
1481 
1496  {
1497  return m_incremental;
1498  }
1499 
1505  void EnableIncrementalFunction( bool enable = true )
1506  {
1507  m_incremental = enable;
1508  }
1509 
1515  void DisableIncrementalFunction( bool disable = true )
1516  {
1517  EnableIncrementalFunction( !disable );
1518  }
1519 
1523  DPoint operator ()( double x, double y ) const
1524  {
1525  DPoint p( m_Sx( x, y ), m_Sy( x, y ) );
1526  if ( m_incremental )
1527  p += m_H( x, y );
1528  return p;
1529  }
1530 
1534  template <typename T>
1535  DPoint operator ()( const GenericPoint<T>& p ) const
1536  {
1537  return operator ()( double( p.x ), double( p.y ) );
1538  }
1539 
1540 private:
1541 
1542  /*
1543  * The surface splines in the X and Y plane directions.
1544  */
1545  spline m_Sx, m_Sy;
1546  int m_maxSplinePoints = __PCL_PSSPLINE_DEFAULT_MAX_POINTS;
1547 
1548  /*
1549  * Incremental surface splines.
1550  */
1551  bool m_incremental = false; // true => fit differences w.r.t a projective transformation
1552  Homography m_H; // base projective transformation when m_incremental = true
1553 
1554  /*
1555  * Surface simplification.
1556  */
1557  bool m_useSimplifiers = false;
1558  float m_simplifierRejectFraction = 0.1F;
1559  double m_epsX = 0; // simplification error on the X axis
1560  double m_epsY = 0; // simplification error on the Y axis
1561  bool m_truncatedX = false; // true => truncated vectors in the X axis
1562  bool m_truncatedY = false; // true => truncated vectors in the Y axis
1563 
1564  template <class weight_vector>
1565  class SplineGenerationThread : public Thread
1566  {
1567  public:
1568 
1569  SplineGenerationThread( spline& S, const DVector& X, const DVector& Y, const DVector& Z, int N, const weight_vector& W )
1570  : m_S( S ), m_X( X ), m_Y( Y ), m_Z( Z ), m_N( N ), m_W( W )
1571  {
1572  }
1573 
1574  PCL_HOT_FUNCTION void Run() override
1575  {
1576  try
1577  {
1578  m_S.Initialize( m_X.Begin(), m_Y.Begin(), m_Z.Begin(), m_N, m_W.Begin() );
1579  return;
1580  }
1581  catch ( const Exception& x )
1582  {
1583  m_exception = x;
1584  }
1585  catch ( ... )
1586  {
1587  m_exception = Error( "Unknown exception" );
1588  }
1589  m_S.Clear();
1590  }
1591 
1592  void ValidateOrThrow() const
1593  {
1594  if ( !m_S.IsValid() )
1595  throw m_exception;
1596  }
1597 
1598  private:
1599 
1600  spline& m_S;
1601  const DVector& m_X;
1602  const DVector& m_Y;
1603  const DVector& m_Z;
1604  int m_N;
1605  weight_vector m_W; // ### N.B. do not use a reference: the W ctor. argument can be a temporary object.
1606  Exception m_exception;
1607  };
1608 
1609  friend class DrizzleData;
1610  friend class DrizzleDataDecoder;
1611  friend class SplineWorldTransformation;
1612 };
1613 
1614 // ----------------------------------------------------------------------------
1615 
1634 {
1635 public:
1636 
1642 
1650 
1658 
1663 
1668 
1676  template <class point_list1, class point_list2, class weight_vector = FVector>
1677  RecursivePointSurfaceSpline( const point_list1& P1, const point_list2& P2,
1678  float smoothness = 0,
1679  int order = 2,
1680  const weight_vector& W = weight_vector(),
1681  bool allowExtrapolation = __PCL_RSSPLINE_DEFAULT_ALLOW_EXTRAPOLATION,
1682  int maxSplineLength = __PCL_RSSPLINE_DEFAULT_SPLINE_MAX_LENGTH,
1683  int bucketCapacity = __PCL_RSSPLINE_DEFAULT_TREE_BUCKET_CAPACITY,
1684  bool verbose = true )
1685  {
1686  Initialize( P1, P2, smoothness, W, order, allowExtrapolation, maxSplineLength, bucketCapacity, verbose );
1687  }
1688 
1693  {
1694  Clear();
1695  }
1696 
1773  template <class point_list1, class point_list2, class weight_vector = FVector>
1774  void Initialize( const point_list1& P1, const point_list2& P2,
1775  float smoothness = 0,
1776  const weight_vector& W = weight_vector(),
1777  int order = 2,
1778  bool allowExtrapolation = __PCL_RSSPLINE_DEFAULT_ALLOW_EXTRAPOLATION,
1779  int maxSplineLength = __PCL_RSSPLINE_DEFAULT_SPLINE_MAX_LENGTH,
1780  int bucketCapacity = __PCL_RSSPLINE_DEFAULT_TREE_BUCKET_CAPACITY,
1781  bool verbose = true )
1782  {
1783  PCL_PRECONDITION( P1.Length() >= 3 )
1784  PCL_PRECONDITION( P1.Length() <= P2.Length() )
1785  PCL_PRECONDITION( order >= 2 )
1786  PCL_PRECONDITION( W.IsEmpty() || P1.Length() <= W.Length() )
1787 
1788  Clear();
1789 
1790  if ( P1.Length() < 3 || P2.Length() < 3 )
1791  throw Error( "RecursivePointSurfaceSpline::Initialize(): At least three input nodes must be specified." );
1792 
1793  bool weighted = smoothness > 0 && !W.IsEmpty();
1794 
1795  if ( P1.Length() > P2.Length() || weighted && P1.Length() > size_type( W.Length() ) )
1796  throw Error( "RecursivePointSurfaceSpline::Initialize(): Insufficient data." );
1797 
1798  m_extrapolate = allowExtrapolation;
1799 
1800  if ( P1.Length() <= size_type( maxSplineLength ) )
1801  {
1802  StatusMonitor monitor;
1803 #ifndef __PCL_BUILDING_PIXINSIGHT_APPLICATION
1804  StandardStatus status;
1805  if ( verbose )
1806  {
1807  monitor.SetCallback( &status );
1808  monitor.Initialize( "Building surface subsplines", 100 );
1809  }
1810 #endif
1811  for ( const auto& p : P1 )
1812  {
1813  double x = double( p.x );
1814  double y = double( p.y );
1815  if ( x < m_rect.x0 )
1816  m_rect.x0 = x;
1817  else if ( x > m_rect.x1 )
1818  m_rect.x1 = x;
1819  if ( y < m_rect.y0 )
1820  m_rect.y0 = y;
1821  else if ( y > m_rect.y1 )
1822  m_rect.y1 = y;
1823  }
1824 
1825  m_spline.Initialize( P1, P2, smoothness, W, order,
1826  (order == 2) ? RadialBasisFunction::ThinPlateSpline
1827  : RadialBasisFunction::VariableOrder );
1828 
1829 #ifndef __PCL_BUILDING_PIXINSIGHT_APPLICATION
1830  if ( verbose )
1831  monitor.Complete();
1832 #endif
1833  }
1834  else
1835  {
1836  /*
1837  * Find the total interpolation region.
1838  */
1839  search_rectangle rect = search_coordinate( 0 );
1840  node_list data;
1841  for ( size_type i = 0; i < P1.Length(); ++i )
1842  {
1843  const auto& p1 = P1[i];
1844  const auto& p2 = P2[i];
1845  data << (weighted ? Node( p1, p2, W[i] ) : Node( p1, p2 ));
1846  double x = p1.x;
1847  double y = p1.y;
1848  if ( x < rect.x0 )
1849  rect.x0 = x;
1850  else if ( x > rect.x1 )
1851  rect.x1 = x;
1852  if ( y < rect.y0 )
1853  rect.y0 = y;
1854  else if ( y > rect.y1 )
1855  rect.y1 = y;
1856  }
1857  // Force a square interpolation region for optimal quadtree behavior.
1858  if ( rect.Width() < rect.Height() )
1859  rect.InflateBy( (rect.Height() - rect.Width())/2, search_coordinate( 0 ) );
1860  else
1861  rect.InflateBy( search_coordinate( 0 ), (rect.Width() - rect.Height())/2 );
1862 
1863  /*
1864  * Build the interpolation quadtree.
1865  */
1866  m_tree.Build( rect, data, bucketCapacity );
1867 
1868  /*
1869  * Build subspline interpolation data.
1870  */
1871  Array<SubsplineData> subsplineData;
1872  m_tree.Traverse(
1873  [&]( const search_rectangle& rect, const node_list& points, void*& D )
1874  {
1875  /*
1876  * Ensure redundancy for all subsplines by gathering a set of
1877  * interpolation points in an extended rectangular region around
1878  * each quadtree node.
1879  */
1880  search_rectangle r = rect.InflatedBy( 1.5*Max( rect.Width(), rect.Height() ) );
1881  node_list N = m_tree.Search( r );
1882 
1883  /*
1884  * Sort the cloud of interpolation points by proximity to a
1885  * circle centered at the argument of the minimum RBF value
1886  * (approximately 0.61 in normalized coordinates).
1887  */
1888  if ( N.Length() > size_type( maxSplineLength ) )
1889  {
1890  DPoint c = rect.Center();
1891  double d = 0.61 * (Max( r.Width(), r.Height() )/2 + Max( rect.Width(), rect.Height() )/4);
1892  N.Sort(
1893  [&]( const auto& a, const auto& b )
1894  {
1895  return Abs( a.position.DistanceTo( c ) - d ) < Abs( b.position.DistanceTo( c ) - d );
1896  } );
1897  }
1898 
1899  /*
1900  * Get the optimal subset of at most maxSplineLength redundant
1901  * interpolation points for this quadtree node.
1902  */
1903  Array<DPoint> P1, P2;
1904  Array<float> PW;
1905  for ( int j = 0, l = Min( int( N.Length() ), maxSplineLength ); j < l; ++j )
1906  {
1907  P1 << N[j].position;
1908  P2 << N[j].value;
1909  if ( weighted )
1910  PW << N[j].weight;
1911  }
1912 
1913  subsplineData << SubsplineData( P1, P2, PW, D );
1914  } );
1915 
1916  /*
1917  * Generate all subsplines.
1918  */
1919  StatusMonitor monitor;
1920 #ifndef __PCL_BUILDING_PIXINSIGHT_APPLICATION
1921  StandardStatus status;
1922  if ( verbose )
1923  {
1924  monitor.SetCallback( &status );
1925  monitor.Initialize( "Building recursive surface subsplines", subsplineData.Length() );
1926  }
1927 #endif
1928  Array<size_type> L = Thread::OptimalThreadLoads( subsplineData.Length(),
1929  1/*overheadLimit*/,
1930  m_parallel ? m_maxProcessors : 1 );
1931  AbstractImage::ThreadData threadData( monitor, subsplineData.Length() );
1933  for ( int i = 0, n = 0; i < int( L.Length() ); n += int( L[i++] ) )
1934  threads.Add( new SubsplineGenerationThread( threadData,
1935  subsplineData,
1936  smoothness,
1937  order,
1938  allowExtrapolation,
1939  maxSplineLength,
1940  bucketCapacity,
1941  n,
1942  n + int( L[i] ) ) );
1943  AbstractImage::RunThreads( threads, threadData );
1944  threads.Destroy();
1945  }
1946  }
1947 
1952  void Clear()
1953  {
1954  m_tree.Traverse(
1955  []( const search_rectangle&, node_list&, void*& D )
1956  {
1957  delete reinterpret_cast<RecursivePointSurfaceSpline*&>( D );
1958  } );
1959  m_tree.Clear();
1960  m_spline.Clear();
1961  m_rect = search_coordinate( 0 );
1962  }
1963 
1974  bool IsRecursive() const
1975  {
1976  return !m_tree.IsEmpty();
1977  }
1978 
1983  bool IsValid() const
1984  {
1985  return IsRecursive() || m_spline.IsValid();
1986  }
1987 
1997  DPoint operator ()( double x, double y ) const
1998  {
1999  if ( m_spline.IsValid() )
2000  {
2001  if ( m_extrapolate || m_rect.IncludesFast( x, y ) )
2002  return m_spline( x, y );
2003  }
2004  else
2005  {
2006  const typename tree::Node* node = m_tree.NodeAt( search_point( x, y ) );
2007  if ( node == nullptr )
2008  {
2009  if ( !m_extrapolate )
2010  return 0.0;
2011 
2012  search_rectangle r0 = m_tree.Root()->rect;
2013  if ( x <= r0.x0 )
2014  {
2015  if ( y <= r0.y0 )
2016  node = m_tree.NodeAt( search_point( r0.x0 + SearchDelta, r0.y0 + SearchDelta ) );
2017  else if ( y >= r0.y1 )
2018  node = m_tree.NodeAt( search_point( r0.x0 + SearchDelta, r0.y1 - SearchDelta ) );
2019  else
2020  node = m_tree.NodeAt( search_point( r0.x0 + SearchDelta, search_coordinate( y ) ) );
2021  }
2022  else if ( x >= r0.x1 )
2023  {
2024  if ( y <= r0.y0 )
2025  node = m_tree.NodeAt( search_point( r0.x1 - SearchDelta, r0.y0 + SearchDelta ) );
2026  else if ( y >= r0.y1 )
2027  node = m_tree.NodeAt( search_point( r0.x1 - SearchDelta, r0.y1 - SearchDelta ) );
2028  else
2029  node = m_tree.NodeAt( search_point( r0.x1 - SearchDelta, search_coordinate( y ) ) );
2030  }
2031  else if ( y <= r0.y0 )
2032  node = m_tree.NodeAt( search_point( search_coordinate( x ), r0.y0 + SearchDelta ) );
2033  else // y >= r0.y1
2034  node = m_tree.NodeAt( search_point( search_coordinate( x ), r0.y1 - SearchDelta ) );
2035 
2036  if ( node == nullptr ) // ?!
2037  return 0.0;
2038  }
2039 
2040  if ( node->IsLeaf() )
2041  {
2042  const typename tree::LeafNode* leaf = static_cast<const typename tree::LeafNode*>( node );
2043  if ( leaf->data == nullptr ) // ?!
2044  return 0.0;
2045  return reinterpret_cast<RecursivePointSurfaceSpline*>( leaf->data )->operator()( x, y );
2046  }
2047 
2048  DPoint p = 0.0;
2049  int n = 0;
2050  if ( node->nw != nullptr )
2051  {
2052  const typename tree::LeafNode* leaf;
2053  if ( y <= node->nw->rect.y1 )
2054  leaf = m_tree.LeafNodeAt( search_point( node->nw->rect.x1 - SearchDelta, search_coordinate( y ) ) );
2055  else if ( x <= node->nw->rect.x1 )
2056  leaf = m_tree.LeafNodeAt( search_point( search_coordinate( x ), node->nw->rect.y1 - SearchDelta ) );
2057  else
2058  leaf = m_tree.LeafNodeAt( search_point( node->nw->rect.x1 - SearchDelta, node->nw->rect.y1 - SearchDelta ) );
2059 
2060  if ( leaf != nullptr )
2061  {
2062  p += reinterpret_cast<RecursivePointSurfaceSpline*>( leaf->data )->operator()( x, y );
2063  ++n;
2064  }
2065  }
2066  if ( node->ne != nullptr )
2067  {
2068  const typename tree::LeafNode* leaf;
2069  if ( y <= node->ne->rect.y1 )
2070  leaf = m_tree.LeafNodeAt( search_point( node->ne->rect.x0 + SearchDelta, search_coordinate( y ) ) );
2071  else if ( x <= node->ne->rect.x0 )
2072  leaf = m_tree.LeafNodeAt( search_point( node->ne->rect.x0 + SearchDelta, node->ne->rect.y1 - SearchDelta ) );
2073  else
2074  leaf = m_tree.LeafNodeAt( search_point( search_coordinate( x ), node->ne->rect.y1 - SearchDelta ) );
2075 
2076  if ( leaf != nullptr )
2077  {
2078  p += reinterpret_cast<RecursivePointSurfaceSpline*>( leaf->data )->operator()( x, y );
2079  ++n;
2080  }
2081  }
2082  if ( node->sw != nullptr )
2083  {
2084  const typename tree::LeafNode* leaf;
2085  if ( y >= node->sw->rect.y0 )
2086  leaf = m_tree.LeafNodeAt( search_point( node->sw->rect.x1 - SearchDelta, search_coordinate( y ) ) );
2087  else if ( x <= node->sw->rect.x1 )
2088  leaf = m_tree.LeafNodeAt( search_point( search_coordinate( x ), node->sw->rect.y0 + SearchDelta ) );
2089  else
2090  leaf = m_tree.LeafNodeAt( search_point( node->sw->rect.x1 - SearchDelta, node->sw->rect.y0 + SearchDelta ) );
2091 
2092  if ( leaf != nullptr )
2093  {
2094  p += reinterpret_cast<RecursivePointSurfaceSpline*>( leaf->data )->operator()( x, y );
2095  ++n;
2096  }
2097  }
2098  if ( node->se != nullptr )
2099  {
2100  const typename tree::LeafNode* leaf;
2101  if ( y >= node->se->rect.y0 )
2102  leaf = m_tree.LeafNodeAt( search_point( node->se->rect.x0 + SearchDelta, search_coordinate( y ) ) );
2103  else if ( x <= node->se->rect.x0 )
2104  leaf = m_tree.LeafNodeAt( search_point( node->se->rect.x0 + SearchDelta, node->se->rect.y0 + SearchDelta ) );
2105  else
2106  leaf = m_tree.LeafNodeAt( search_point( search_coordinate( x ), node->se->rect.y0 + SearchDelta ) );
2107 
2108  if ( leaf != nullptr )
2109  {
2110  p += reinterpret_cast<RecursivePointSurfaceSpline*>( leaf->data )->operator()( x, y );
2111  ++n;
2112  }
2113  }
2114 
2115  if ( n > 0 )
2116  return p/double( n );
2117  }
2118 
2119  return 0.0;
2120  }
2121 
2131  template <typename T>
2132  DPoint operator ()( const GenericPoint<T>& p ) const
2133  {
2134  return operator ()( double( p.x ), double( p.y ) );
2135  }
2136 
2137 private:
2138 
2139  /*
2140  * Interpolation data point, QuadTree-compatible.
2141  */
2142  struct Node
2143  {
2144  DPoint position, value;
2145  float weight;
2146 
2147  using component = DPoint::component;
2148 
2149  template <class point1, class point2>
2150  Node( const point1& p, const point2& v )
2151  : position( double( p.x ), double( p.y ) )
2152  , value( double( v.x ), double( v.y ) )
2153  {
2154  }
2155 
2156  template <class point1, class point2>
2157  Node( const point1& p, const point2& v, float w )
2158  : position( double( p.x ), double( p.y ) )
2159  , value( double( v.x ), double( v.y ) )
2160  , weight( w )
2161  {
2162  }
2163 
2164  Node( const Node& ) = default;
2165 
2166  double& operator []( int i )
2167  {
2168  return position[i];
2169  }
2170 
2171  double operator []( int i ) const
2172  {
2173  return position[i];
2174  }
2175  };
2176 
2177  using tree = QuadTree<Node>;
2178  using node_list = typename tree::point_list;
2179  using search_rectangle = typename tree::rectangle;
2180  using search_coordinate = typename tree::coordinate;
2181  using search_point = GenericPoint<search_coordinate>;
2182 
2183  tree m_tree; // the tree of subsplines
2184  PointSurfaceSpline m_spline; // final point spline if there is no further recursion
2185  search_rectangle m_rect = search_coordinate( 0 ); // the interpolation region for this subspline
2186  bool m_extrapolate = __PCL_RSSPLINE_DEFAULT_ALLOW_EXTRAPOLATION;
2187 
2188  static constexpr search_coordinate SearchDelta =
2189  2 * std::numeric_limits<search_coordinate>::epsilon();
2190 
2191  /*
2192  * Parallel subspline generation.
2193  */
2194  struct SubsplineData
2195  {
2196  Array<DPoint> P1, P2;
2197  Array<float> PW;
2198  mutable void** nodeData;
2199 
2200  SubsplineData( const Array<DPoint>& p1, const Array<DPoint>& p2, const Array<float>& pw, void*& nd )
2201  : P1( p1 )
2202  , P2( p2 )
2203  , PW( pw )
2204  , nodeData( &nd )
2205  {
2206  }
2207  };
2208 
2209  class SubsplineGenerationThread : public Thread
2210  {
2211  public:
2212 
2213  SubsplineGenerationThread( const AbstractImage::ThreadData& data,
2214  const Array<SubsplineData>& subsplineData,
2215  float smoothness,
2216  int order,
2217  bool allowExtrapolation,
2218  int maxSplineLength,
2219  int bucketCapacity,
2220  int startIndex, int endIndex )
2221  : m_data( data )
2222  , m_subsplineData( subsplineData )
2223  , m_smoothness( smoothness )
2224  , m_order( order )
2225  , m_allowExtrapolation( allowExtrapolation )
2226  , m_maxSplineLength( maxSplineLength )
2227  , m_bucketCapacity( bucketCapacity )
2228  , m_startIndex( startIndex )
2229  , m_endIndex( endIndex )
2230  {
2231  }
2232 
2233  PCL_HOT_FUNCTION void Run() override
2234  {
2236 
2237  for ( int i = m_startIndex; i < m_endIndex; ++i )
2238  {
2239  const SubsplineData& d = m_subsplineData[i];
2240  AutoPointer<RecursivePointSurfaceSpline> s(
2241  new RecursivePointSurfaceSpline( d.P1, d.P2, m_smoothness, m_order, d.PW,
2242  m_allowExtrapolation,
2243  m_maxSplineLength,
2244  m_bucketCapacity,
2245  false/*verbose*/ ) );
2246  *reinterpret_cast<RecursivePointSurfaceSpline**>( d.nodeData ) = s.Release();
2247 
2249  }
2250 
2251  m_success = true;
2252  }
2253 
2254  operator bool() const
2255  {
2256  return m_success;
2257  }
2258 
2259  private:
2260 
2261  const AbstractImage::ThreadData& m_data;
2262  const Array<SubsplineData>& m_subsplineData;
2263  float m_smoothness;
2264  int m_order;
2265  bool m_allowExtrapolation;
2266  int m_maxSplineLength;
2267  int m_bucketCapacity;
2268  int m_startIndex, m_endIndex;
2269  bool m_success = false;
2270  };
2271 
2272  friend class DrizzleData;
2273  friend class DrizzleDataDecoder;
2274 };
2275 
2276 // ----------------------------------------------------------------------------
2277 
2278 } // pcl
2279 
2280 #endif // __PCL_SurfaceSpline_h
2281 
2282 // ----------------------------------------------------------------------------
2283 // EOF pcl/SurfaceSpline.h - Released 2024-06-18T15:48:54Z
static void RunThreads(ReferenceArray< thread > &threads, ThreadData &data, bool useAffinity=true)
iterator Begin()
Definition: Array.h:426
size_type Length() const noexcept
Definition: Array.h:266
iterator End()
Definition: Array.h:451
64-bit floating point real vector.
Drizzle integration data parser and generator.
Definition: DrizzleData.h:143
A simple exception with an associated error message.
Definition: Exception.h:239
Generic dynamic matrix of arbitrary dimensions.
Definition: Matrix.h:123
A generic point in the two-dimensional space.
Definition: Point.h:100
component x
Abscissa (horizontal, or X-axis coordinate).
Definition: Point.h:111
component y
Ordinate (vertical, or Y-axis coordinate).
Definition: Point.h:112
iterator Begin()
Definition: Vector.h:1900
int Length() const noexcept
Definition: Vector.h:1784
Homography geometric transformation.
Definition: Homography.h:85
A process using multiple concurrent execution threads.
Vector surface spline interpolation/approximation in two dimensions.
const Matrix & LinearFunction() const
void SetLinearFunction(const Matrix &H)
void EnableSimplifiers(bool enable=true)
PointSurfaceSpline(PointSurfaceSpline &&)=default
void Initialize(const point_list1 &P1, const point_list2 &P2, float smoothness=0, const weight_vector &W=weight_vector(), int order=2, rbf_type rbf=RadialBasisFunction::Default, double eps=0, bool polynomial=true)
bool SimplifiersEnabled() const
float SimplifierRejectFraction() const
const spline & SplineY() const
void DisableIncrementalFunction(bool disable=true)
const spline & SplineX() const
PointSurfaceSpline(const PointSurfaceSpline &)=default
void DisableSimplifiers(bool disable=true)
PointSurfaceSpline(const spline &Sx, const spline &Sy)
void EnableIncrementalFunction(bool enable=true)
void SetSimplifierRejectFraction(float rejectFraction)
PointSurfaceSpline(const point_list1 &P1, const point_list2 &P2, float smoothness=0, int order=2, const weight_vector &W=weight_vector(), rbf_type rbf=RadialBasisFunction::Default, double eps=0, bool polynomial=true)
bool IncrementalFunctionEnabled() const
void SetMaxSplinePoints(int n)
void Initialize(const spline &Sx, const spline &Sy)
Vector surface spline interpolation/approximation in two dimensions with recursive subspline generati...
RecursivePointSurfaceSpline(RecursivePointSurfaceSpline &&)=default
RecursivePointSurfaceSpline(const RecursivePointSurfaceSpline &)=delete
void Initialize(const point_list1 &P1, const point_list2 &P2, float smoothness=0, const weight_vector &W=weight_vector(), int order=2, bool allowExtrapolation=__PCL_RSSPLINE_DEFAULT_ALLOW_EXTRAPOLATION, int maxSplineLength=__PCL_RSSPLINE_DEFAULT_SPLINE_MAX_LENGTH, int bucketCapacity=__PCL_RSSPLINE_DEFAULT_TREE_BUCKET_CAPACITY, bool verbose=true)
RecursivePointSurfaceSpline(const point_list1 &P1, const point_list2 &P2, float smoothness=0, int order=2, const weight_vector &W=weight_vector(), bool allowExtrapolation=__PCL_RSSPLINE_DEFAULT_ALLOW_EXTRAPOLATION, int maxSplineLength=__PCL_RSSPLINE_DEFAULT_SPLINE_MAX_LENGTH, int bucketCapacity=__PCL_RSSPLINE_DEFAULT_TREE_BUCKET_CAPACITY, bool verbose=true)
Dynamic array of pointers to objects providing direct iteration and element access by reference.
void Destroy(iterator i, size_type n=1)
void Add(const ReferenceArray &x)
Surface spline world coordinate transformation.
A status monitoring callback that sends progress information to the process console.
An asynchronous status monitoring system.
void SetCallback(StatusCallback *callback)
void Initialize(const String &info, size_type count=0)
Shape-preserving simplification of 2-D surfaces.
void Simplify(C &xs, C &ys, C &zs, const C &x, const C &y, const C &z) const
void EnableCentroidInclusion(bool enable=true)
void EnableRejection(bool enabled=true)
void SetRejectFraction(float rejectFraction)
void SetTolerance(double tolerance)
Base class of two-dimensional surface splines.
SurfaceSplineBase(SurfaceSplineBase &&)=default
static void Generate(float *__restrict__ c, rbf_type, double e2, bool polynomial, const float *__restrict__ x, const float *__restrict__ y, const float *__restrict__ z, int n, int m, float r, const float *__restrict__ w)
static void Generate(double *__restrict__ c, rbf_type, double e2, bool polynomial, const double *__restrict__ x, const double *__restrict__ y, const double *__restrict__ z, int n, int m, float r, const float *__restrict__ w)
RadialBasisFunction::value_type rbf_type
SurfaceSplineBase(const SurfaceSplineBase &)=default
Two-dimensional interpolating/approximating surface spline.
SurfaceSpline()=default
void DisablePolynomial(bool disable=true)
int Length() const
vector X() const
bool IsPolynomialEnabled() const
~SurfaceSpline() override
rbf_type RBFType() const
void SetShapeParameter(double eps)
void EnablePolynomial(bool enable=true)
SurfaceSpline(const SurfaceSpline &)=default
const weight_vector & Weights() const
void SetRBFType(rbf_type rbf)
vector Y() const
const vector & Z() const
float Smoothing() const
SurfaceSpline(SurfaceSpline &&)=default
void SetOrder(int order)
bool IsValid() const
void Initialize(const scalar *x, const scalar *y, const scalar *z, int n, const weight *w=nullptr)
void SetSmoothing(float s)
double ShapeParameter() const
Client-side interface to a PixInsight thread.
Definition: Thread.h:130
static Array< size_type > OptimalThreadLoads(size_type count, size_type overheadLimit=1u, int maxThreads=PCL_MAX_PROCESSORS)
bool operator==(const Array< T, A > &x1, const Array< T, A > &x2) noexcept
Definition: Array.h:2267
bool operator<(const Array< T, A > &x1, const Array< T, A > &x2) noexcept
Definition: Array.h:2278
Complex< T > Sqrt(const Complex< T > &c) noexcept
Definition: Complex.h:674
Complex< T > Exp(const Complex< T > &c) noexcept
Definition: Complex.h:714
Complex< T > Ln(const Complex< T > &c) noexcept
Definition: Complex.h:725
T Abs(const Complex< T > &c) noexcept
Definition: Complex.h:429
T PowI(T x, int n) noexcept
Definition: Math.h:1755
size_t size_type
Definition: Defs.h:609
double Median(const T *__restrict__ i, const T *__restrict__ j)
Definition: Math.h:2878
#define INIT_THREAD_MONITOR()
Declares and initializes local variables used for synchronization of thread status monitoring.
#define UPDATE_THREAD_MONITOR(N)
Synchronized status monitoring of a set of image processing threads.
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
Thread synchronization data for status monitoring of parallel image processing tasks.
Auxiliary structure for data sanitization.