52 #ifndef __PCL_SurfaceSpline_h
53 #define __PCL_SurfaceSpline_h
58 #include <pcl/Diagnostics.h>
71 #ifndef __PCL_BUILDING_PIXINSIGHT_APPLICATION
84 #define __PCL_PSSPLINE_DEFAULT_MAX_POINTS 2100
89 #define __PCL_RSSPLINE_DEFAULT_TREE_BUCKET_CAPACITY 64
95 #define __PCL_RSSPLINE_DEFAULT_SPLINE_MAX_LENGTH 1600
103 #define __PCL_RSSPLINE_DEFAULT_ALLOW_EXTRAPOLATION false
129 namespace RadialBasisFunction
141 Default = ThinPlateSpline
144 inline static bool Validate(
int rbf )
146 return rbf >= 0 && rbf < __number_of_items__;
164 using rbf_type = RadialBasisFunction::value_type;
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 );
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 );
266 template <
typename T>
339 return m_x.Length() >= 3 && m_x.Length() == m_y.Length();
359 for (
int i = 0; i < m_x.Length(); ++i )
360 x[i] = m_x0 + m_x[i]/m_r0;
373 for (
int i = 0; i < m_y.Length(); ++i )
374 y[i] = m_y0 + m_y[i]/m_r0;
433 return Sqrt( m_eps2 );
455 m_eps2 = m_eps*m_eps;
489 PCL_PRECONDITION( order > 1 )
504 return m_havePolynomial;
517 m_havePolynomial = enable;
529 EnablePolynomial( !disable );
560 PCL_PRECONDITION( s >= 0 )
617 PCL_PRECONDITION( x !=
nullptr && y !=
nullptr && z !=
nullptr )
618 PCL_PRECONDITION( n > 2 )
620 if ( x ==
nullptr || y ==
nullptr || z ==
nullptr )
621 throw Error(
"SurfaceSpline::Initialize(): Null vector pointer(s)." );
624 throw Error(
"SurfaceSpline::Initialize(): At least three input nodes must be specified." );
628 if ( m_smoothing <= 0 )
637 for (
int i = 0; i < n; ++i )
649 for (
int i = 0; i < n; ++i )
651 double dx = x[i] - m_x0;
652 double dy = y[i] - m_y0;
653 double r =
Sqrt( dx*dx + dy*dy );
658 throw Error(
"SurfaceSpline::Initialize(): Empty or insignificant interpolation space." );
667 for (
int i = 0; i < n; ++i )
668 for (
int j = i; ++j < n; )
670 double dx = x[i] - x[j];
671 double dy = y[i] - y[j];
684 for (
int i = 0; i < n; ++i )
688 (w !=
nullptr && w[i] > 0) ? w[i] : 1.0F );
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() )
705 int N = n - int( remove.
Length() );
707 throw Error(
"SurfaceSpline::Initialize(): Less than three input nodes left after sanitization." );
714 for (
int j : remove )
716 for ( ; i < j; ++i, ++k )
722 m_weights[k] = P[i].w;
726 for ( ; i < n; ++i, ++k )
732 m_weights[k] = P[i].w;
735 m_spline =
vector(
scalar( 0 ), N + (m_havePolynomial ? ((m_order*(m_order + 1)) >> 1) : 0) );
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() );
770 scalar operator ()(
double x,
double y )
const
772 PCL_PRECONDITION( !m_x.IsEmpty() && !m_y.IsEmpty() )
773 PCL_CHECK( m_order >= 2 )
774 PCL_CHECK( !m_spline.IsEmpty() )
785 int n = m_x.Length();
787 if ( m_havePolynomial )
793 z += m_spline[n+1]*x + m_spline[n+2]*y;
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;
799 for (
int i = 1, j = n+1, i1 = (m_order*(m_order + 1))>>1, ix = 0, iy = 0; i < i1; ++i, ++j )
804 z += m_spline[j] *
PowI( x, ix );
810 z += m_spline[j] *
PowI( x, ix ) *
PowI( y, iy );
821 case RadialBasisFunction::VariableOrder:
822 for (
int i = 0; i < n; ++i )
824 double dx = m_x[i] - x;
825 double dy = m_y[i] - y;
826 double r2 = dx*dx + dy*dy;
830 for (
int j = m_order; --j > 1; )
832 z += m_spline[i] * r2m1 *
pcl::Ln( r2 );
836 case RadialBasisFunction::ThinPlateSpline:
837 for (
int i = 0; i < n; ++i )
839 double dx = m_x[i] - x;
840 double dy = m_y[i] - y;
841 double r2 = dx*dx + dy*dy;
846 case RadialBasisFunction::Gaussian:
847 for (
int i = 0; i < n; ++i )
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) );
854 case RadialBasisFunction::Multiquadric:
855 for (
int i = 0; i < n; ++i )
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) );
862 case RadialBasisFunction::InverseMultiquadric:
863 for (
int i = 0; i < n; ++i )
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) );
870 case RadialBasisFunction::InverseQuadratic:
871 for (
int i = 0; i < n; ++i )
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));
890 template <
typename Tp>
893 return operator ()(
double( p.
x ),
double( p.
y ) );
898 rbf_type m_rbf = RadialBasisFunction::Default;
899 bool m_havePolynomial =
true;
909 float m_smoothing = 0;
910 weight_vector m_weights;
924 : x( a_x ), y( a_y ), z( a_z ), w( a_w )
930 return x == p.x && y == p.y;
935 return (x != p.x) ? x < p.x : y < p.y;
940 friend class DrizzleDataDecoder;
972 using rbf_type = spline::rbf_type;
997 template <
class po
int_list1,
class po
int_list2,
class weight_vector = FVector>
999 float smoothness = 0,
int order = 2,
1000 const weight_vector& W = weight_vector(),
1001 rbf_type rbf = RadialBasisFunction::Default,
1003 bool polynomial =
true )
1005 Initialize( P1, P2, smoothness, W, order, rbf, eps, polynomial );
1017 Initialize( Sx, Sy );
1107 template <
class po
int_list1,
class po
int_list2,
class weight_vector = FVector>
1109 float smoothness = 0,
const weight_vector& W = weight_vector(),
int order = 2,
1110 rbf_type rbf = RadialBasisFunction::Default,
1112 bool polynomial =
true )
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() ) )
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." );
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 );
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 );
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 )
1146 X[i] = double( P1[i].x );
1147 Y[i] = double( P1[i].y );
1149 if ( m_incremental )
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];
1159 ZX[i] = double( P2[i].x );
1160 ZY[i] = double( P2[i].y );
1163 if ( ZX[i] < zxMin )
1165 if ( ZX[i] > zxMax )
1167 if ( ZY[i] < zyMin )
1169 if ( ZY[i] > zyMax )
1173 if ( m_useSimplifiers )
1182 m_epsX = (zxMax - zxMin)/100;
1183 double epsLow = 0, epsHigh = (zxMax - zxMin)/10;
1184 for (
int i = 0; i < 200; ++i )
1187 SS.
Simplify( XXS, YXS, ZXS, X, Y, ZX );
1188 if ( XXS.
Length() > m_maxSplinePoints )
1191 if ( epsHigh - epsLow <= 4*std::numeric_limits<double>::epsilon() )
1196 if ( epsHigh - epsLow <= 2*std::numeric_limits<double>::epsilon() )
1200 m_epsX = (epsLow + epsHigh)/2;
1202 if ( XXS.
Length() > m_maxSplinePoints )
1204 m_truncatedX =
true;
1212 m_epsY = (zyMax - zyMin)/100;
1213 epsLow = 0, epsHigh = (zyMax - zyMin)/10;
1214 for (
int i = 0; i < 200; ++i )
1217 SS.
Simplify( XYS, YYS, ZYS, X, Y, ZY );
1218 if ( XYS.
Length() > m_maxSplinePoints )
1221 if ( epsHigh - epsLow <= 4*std::numeric_limits<double>::epsilon() )
1226 if ( epsHigh - epsLow <= 2*std::numeric_limits<double>::epsilon() )
1230 m_epsY = (epsLow + epsHigh)/2;
1232 if ( XYS.
Length() > m_maxSplinePoints )
1234 m_truncatedY =
true;
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 );
1246 Tx.ValidateOrThrow();
1247 Ty.ValidateOrThrow();
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 );
1258 Tx.ValidateOrThrow();
1259 Ty.ValidateOrThrow();
1282 m_useSimplifiers = m_incremental =
false;
1299 m_epsX = m_epsY = 0;
1300 m_truncatedX = m_truncatedY =
false;
1309 return m_Sx.IsValid() && m_Sy.IsValid() && (!m_incremental || m_H.IsValid());
1336 return m_maxSplinePoints;
1349 PCL_PRECONDITION( n >= 3 )
1350 m_maxSplinePoints =
Max( 3, n );
1361 return m_useSimplifiers;
1373 m_useSimplifiers = enable;
1385 EnableSimplifiers( !disable );
1395 return m_simplifierRejectFraction;
1406 PCL_PRECONDITION( rejectFraction > 0 && rejectFraction < 1 )
1407 m_simplifierRejectFraction =
Range( rejectFraction, 0.0F, 1.0F );
1435 return m_truncatedX;
1445 return m_truncatedY;
1455 return TruncatedX() || TruncatedY();
1497 return m_incremental;
1507 m_incremental = enable;
1517 EnableIncrementalFunction( !disable );
1523 DPoint operator ()(
double x,
double y )
const
1525 DPoint p( m_Sx( x, y ), m_Sy( x, y ) );
1526 if ( m_incremental )
1534 template <
typename T>
1537 return operator ()(
double( p.
x ),
double( p.
y ) );
1546 int m_maxSplinePoints = __PCL_PSSPLINE_DEFAULT_MAX_POINTS;
1551 bool m_incremental =
false;
1557 bool m_useSimplifiers =
false;
1558 float m_simplifierRejectFraction = 0.1F;
1561 bool m_truncatedX =
false;
1562 bool m_truncatedY =
false;
1564 template <
class weight_vector>
1565 class SplineGenerationThread :
public Thread
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 )
1574 PCL_HOT_FUNCTION
void Run()
override
1578 m_S.Initialize( m_X.Begin(), m_Y.Begin(), m_Z.Begin(), m_N, m_W.Begin() );
1581 catch (
const Exception& x )
1587 m_exception = Error(
"Unknown exception" );
1592 void ValidateOrThrow()
const
1594 if ( !m_S.IsValid() )
1606 Exception m_exception;
1609 friend class DrizzleData;
1610 friend class DrizzleDataDecoder;
1611 friend class SplineWorldTransformation;
1676 template <
class po
int_list1,
class po
int_list2,
class weight_vector = FVector>
1678 float smoothness = 0,
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 )
1686 Initialize( P1, P2, smoothness, W, order, allowExtrapolation, maxSplineLength, bucketCapacity, verbose );
1773 template <
class po
int_list1,
class po
int_list2,
class weight_vector = FVector>
1775 float smoothness = 0,
1776 const weight_vector& W = weight_vector(),
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 )
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() )
1790 if ( P1.Length() < 3 || P2.Length() < 3 )
1791 throw Error(
"RecursivePointSurfaceSpline::Initialize(): At least three input nodes must be specified." );
1793 bool weighted = smoothness > 0 && !W.IsEmpty();
1795 if ( P1.Length() > P2.Length() || weighted && P1.Length() >
size_type( W.Length() ) )
1796 throw Error(
"RecursivePointSurfaceSpline::Initialize(): Insufficient data." );
1798 m_extrapolate = allowExtrapolation;
1800 if ( P1.Length() <=
size_type( maxSplineLength ) )
1803 #ifndef __PCL_BUILDING_PIXINSIGHT_APPLICATION
1808 monitor.
Initialize(
"Building surface subsplines", 100 );
1811 for (
const auto& p : P1 )
1813 double x = double( p.x );
1814 double y = double( p.y );
1815 if ( x < m_rect.x0 )
1817 else if ( x > m_rect.x1 )
1819 if ( y < m_rect.y0 )
1821 else if ( y > m_rect.y1 )
1825 m_spline.Initialize( P1, P2, smoothness, W, order,
1826 (order == 2) ? RadialBasisFunction::ThinPlateSpline
1827 : RadialBasisFunction::VariableOrder );
1829 #ifndef __PCL_BUILDING_PIXINSIGHT_APPLICATION
1839 search_rectangle rect = search_coordinate( 0 );
1841 for (
size_type i = 0; i < P1.Length(); ++i )
1843 const auto& p1 = P1[i];
1844 const auto& p2 = P2[i];
1845 data << (weighted ? Node( p1, p2, W[i] ) : Node( p1, p2 ));
1850 else if ( x > rect.x1 )
1854 else if ( y > rect.y1 )
1858 if ( rect.Width() < rect.Height() )
1859 rect.InflateBy( (rect.Height() - rect.Width())/2, search_coordinate( 0 ) );
1861 rect.InflateBy( search_coordinate( 0 ), (rect.Width() - rect.Height())/2 );
1866 m_tree.Build( rect, data, bucketCapacity );
1873 [&](
const search_rectangle& rect,
const node_list& points,
void*& D )
1880 search_rectangle r = rect.InflatedBy( 1.5*
Max( rect.Width(), rect.Height() ) );
1881 node_list N = m_tree.Search( r );
1888 if ( N.Length() >
size_type( maxSplineLength ) )
1890 DPoint c = rect.Center();
1891 double d = 0.61 * (
Max( r.Width(), r.Height() )/2 +
Max( rect.Width(), rect.Height() )/4);
1893 [&](
const auto& a,
const auto& b )
1895 return Abs( a.position.DistanceTo( c ) - d ) <
Abs( b.position.DistanceTo( c ) - d );
1905 for (
int j = 0, l =
Min(
int( N.Length() ), maxSplineLength ); j < l; ++j )
1907 P1 << N[j].position;
1913 subsplineData << SubsplineData( P1, P2, PW, D );
1920 #ifndef __PCL_BUILDING_PIXINSIGHT_APPLICATION
1925 monitor.
Initialize(
"Building recursive surface subsplines", subsplineData.
Length() );
1930 m_parallel ? m_maxProcessors : 1 );
1933 for (
int i = 0, n = 0; i < int( L.Length() ); n += int( L[i++] ) )
1934 threads.
Add(
new SubsplineGenerationThread( threadData,
1942 n +
int( L[i] ) ) );
1955 [](
const search_rectangle&, node_list&,
void*& D )
1961 m_rect = search_coordinate( 0 );
1976 return !m_tree.IsEmpty();
1985 return IsRecursive() || m_spline.IsValid();
1997 DPoint operator ()(
double x,
double y )
const
1999 if ( m_spline.IsValid() )
2001 if ( m_extrapolate || m_rect.IncludesFast( x, y ) )
2002 return m_spline( x, y );
2006 const typename tree::Node* node = m_tree.NodeAt(
search_point( x, y ) );
2007 if ( node ==
nullptr )
2009 if ( !m_extrapolate )
2012 search_rectangle r0 = m_tree.Root()->rect;
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 ) );
2020 node = m_tree.NodeAt(
search_point( r0.x0 + SearchDelta, search_coordinate( y ) ) );
2022 else if ( x >= r0.x1 )
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 ) );
2029 node = m_tree.NodeAt(
search_point( r0.x1 - SearchDelta, search_coordinate( y ) ) );
2031 else if ( y <= r0.y0 )
2032 node = m_tree.NodeAt(
search_point( search_coordinate( x ), r0.y0 + SearchDelta ) );
2034 node = m_tree.NodeAt(
search_point( search_coordinate( x ), r0.y1 - SearchDelta ) );
2036 if ( node ==
nullptr )
2040 if ( node->IsLeaf() )
2042 const typename tree::LeafNode* leaf =
static_cast<const typename tree::LeafNode*
>( node );
2043 if ( leaf->data ==
nullptr )
2050 if ( node->nw !=
nullptr )
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 ) );
2058 leaf = m_tree.LeafNodeAt(
search_point( node->nw->rect.x1 - SearchDelta, node->nw->rect.y1 - SearchDelta ) );
2060 if ( leaf !=
nullptr )
2066 if ( node->ne !=
nullptr )
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 ) );
2074 leaf = m_tree.LeafNodeAt(
search_point( search_coordinate( x ), node->ne->rect.y1 - SearchDelta ) );
2076 if ( leaf !=
nullptr )
2082 if ( node->sw !=
nullptr )
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 ) );
2090 leaf = m_tree.LeafNodeAt(
search_point( node->sw->rect.x1 - SearchDelta, node->sw->rect.y0 + SearchDelta ) );
2092 if ( leaf !=
nullptr )
2098 if ( node->se !=
nullptr )
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 ) );
2106 leaf = m_tree.LeafNodeAt(
search_point( search_coordinate( x ), node->se->rect.y0 + SearchDelta ) );
2108 if ( leaf !=
nullptr )
2116 return p/double( n );
2131 template <
typename T>
2134 return operator ()(
double( p.
x ),
double( p.
y ) );
2149 template <
class po
int1,
class po
int2>
2150 Node(
const point1& p,
const point2& v )
2151 : position( double( p.x ), double( p.y ) )
2152 , value( double( v.x ), double( v.y ) )
2156 template <
class po
int1,
class po
int2>
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 ) )
2164 Node(
const Node& ) =
default;
2166 double& operator [](
int i )
2171 double operator [](
int i )
const
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>;
2184 PointSurfaceSpline m_spline;
2185 search_rectangle m_rect = search_coordinate( 0 );
2186 bool m_extrapolate = __PCL_RSSPLINE_DEFAULT_ALLOW_EXTRAPOLATION;
2188 static constexpr search_coordinate SearchDelta =
2189 2 * std::numeric_limits<search_coordinate>::epsilon();
2194 struct SubsplineData
2196 Array<DPoint> P1, P2;
2198 mutable void** nodeData;
2200 SubsplineData(
const Array<DPoint>& p1,
const Array<DPoint>& p2,
const Array<float>& pw,
void*& nd )
2209 class SubsplineGenerationThread :
public Thread
2213 SubsplineGenerationThread(
const AbstractImage::ThreadData& data,
2214 const Array<SubsplineData>& subsplineData,
2217 bool allowExtrapolation,
2218 int maxSplineLength,
2220 int startIndex,
int endIndex )
2222 , m_subsplineData( subsplineData )
2223 , m_smoothness( smoothness )
2225 , m_allowExtrapolation( allowExtrapolation )
2226 , m_maxSplineLength( maxSplineLength )
2227 , m_bucketCapacity( bucketCapacity )
2228 , m_startIndex( startIndex )
2229 , m_endIndex( endIndex )
2233 PCL_HOT_FUNCTION
void Run()
override
2237 for (
int i = m_startIndex; i < m_endIndex; ++i )
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,
2246 *
reinterpret_cast<RecursivePointSurfaceSpline**
>( d.nodeData ) = s.Release();
2254 operator bool()
const
2261 const AbstractImage::ThreadData& m_data;
2262 const Array<SubsplineData>& m_subsplineData;
2265 bool m_allowExtrapolation;
2266 int m_maxSplineLength;
2267 int m_bucketCapacity;
2268 int m_startIndex, m_endIndex;
2269 bool m_success =
false;
2272 friend class DrizzleData;
2273 friend class DrizzleDataDecoder;
static void RunThreads(ReferenceArray< thread > &threads, ThreadData &data, bool useAffinity=true)
size_type Length() const noexcept
64-bit floating point real vector.
Drizzle integration data parser and generator.
A simple exception with an associated error message.
Generic dynamic matrix of arbitrary dimensions.
A generic point in the two-dimensional space.
component x
Abscissa (horizontal, or X-axis coordinate).
component y
Ordinate (vertical, or Y-axis coordinate).
int Length() const noexcept
Homography geometric transformation.
A process using multiple concurrent execution threads.
Vector surface spline interpolation/approximation in two dimensions.
int MaxSplinePoints() const
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)
PointSurfaceSpline()=default
Vector surface spline interpolation/approximation in two dimensions with recursive subspline generati...
RecursivePointSurfaceSpline(RecursivePointSurfaceSpline &&)=default
RecursivePointSurfaceSpline(const RecursivePointSurfaceSpline &)=delete
~RecursivePointSurfaceSpline() override
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()=default
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)
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
virtual ~SurfaceSplineBase()
SurfaceSplineBase()=default
SurfaceSplineBase(const SurfaceSplineBase &)=default
Two-dimensional interpolating/approximating surface spline.
void DisablePolynomial(bool disable=true)
bool IsPolynomialEnabled() const
~SurfaceSpline() override
void SetShapeParameter(double eps)
void EnablePolynomial(bool enable=true)
SurfaceSpline(const SurfaceSpline &)=default
const weight_vector & Weights() const
void SetRBFType(rbf_type rbf)
SurfaceSpline(SurfaceSpline &&)=default
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.
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
bool operator<(const Array< T, A > &x1, const Array< T, A > &x2) noexcept
Complex< T > Sqrt(const Complex< T > &c) noexcept
Complex< T > Exp(const Complex< T > &c) noexcept
Complex< T > Ln(const Complex< T > &c) noexcept
T Abs(const Complex< T > &c) noexcept
T PowI(T x, int n) noexcept
double Median(const T *__restrict__ i, const T *__restrict__ j)
#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
constexpr const T & Range(const T &x, const T &a, const T &b) noexcept
constexpr const T & Max(const T &a, const T &b) noexcept
Thread synchronization data for status monitoring of parallel image processing tasks.
Auxiliary structure for data sanitization.