PCL
SurfaceSimplifier.h
Go to the documentation of this file.
1 // ____ ______ __
2 // / __ \ / ____// /
3 // / /_/ // / / /
4 // / ____// /___ / /___ PixInsight Class Library
5 // /_/ \____//_____/ PCL 2.7.0
6 // ----------------------------------------------------------------------------
7 // pcl/SurfaceSimplifier.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_SurfaceSimplifier_h
53 #define __PCL_SurfaceSimplifier_h
54 
56 
57 #include <pcl/Defs.h>
58 #include <pcl/Diagnostics.h>
59 
60 #include <pcl/Matrix.h>
61 #include <pcl/QuadTree.h>
62 
63 namespace pcl
64 {
65 
66 // ----------------------------------------------------------------------------
67 
68 /*
69  * Default tolerance for surface simplification.
70  *
71  * The tolerance parameter determines the maximum absolute deviation of
72  * function values from a plane allowed to simplify the subset of points within
73  * a subregion of the input point space. This parameter is specified in
74  * bivariate function values, i.e. in Z-axis units.
75  */
76 #define __PCL_SURFACE_SIMPLIFIER_DEFAULT_TOLERANCE 0.01
77 
78 /*
79  * Default state of the surface simplification outlier rejection feature.
80  *
81  * When enabled, a prescribed fraction of outlier points will be rejected on
82  * each subregion for estimation of local curvature.
83  */
84 #define __PCL_SURFACE_SIMPLIFIER_DEFAULT_REJECTION_ENABLED true
85 
86 /*
87  * Default rejection fraction for surface simplification.
88  *
89  * This parameter specifies the maximum fraction of outliers allowed for
90  * simplification of a subset of points in a region of the input point space.
91  * Point rejection makes the surface simplification algorithm robust to outlier
92  * function values.
93  */
94 #define __PCL_SURFACE_SIMPLIFIER_DEFAULT_REJECT_FRACTION 0.2F
95 
96 /*
97  * Whether to include subregion centroids in simplified point sets.
98  *
99  * When enabled, the centroid of each simplified subregion will also be
100  * included in the corresponding list of simplified points. This can improve
101  * the shape preservation behavior of the surface simplification algorithm, at
102  * the cost of a reduced amount of additional points.
103  */
104 #define __PCL_SURFACE_SIMPLIFIER_DEFAULT_INCLUDE_CENTROIDS false
105 
106 // ----------------------------------------------------------------------------
107 
145 class PCL_CLASS SurfaceSimplifier
146 {
147 private:
148 
155  struct point
156  {
157  using component = double;
158 
159  component x, y, z;
160 
161  point() = default;
162  point( const point& ) = default;
163 
164  template <typename T>
165  point( T a_x, T a_y, T a_z )
166  : x( double( a_x ) )
167  , y( double( a_y ) )
168  , z( double( a_z ) )
169  {
170  }
171 
172  template <typename T>
173  point( T k )
174  {
175  x = y = z = double( k );
176  }
177 
178  component operator[]( int i ) const
179  {
180  return i ? y : x;
181  }
182 
183  point& operator +=( const point& p )
184  {
185  x += p.x; y += p.y; z += p.z;
186  return *this;
187  }
188 
189  template <typename T>
190  point& operator /=( T k )
191  {
192  x /= k; y /= k; z /= k;
193  return *this;
194  }
195 
196  double SquaredDistanceTo( const point& p ) const
197  {
198  double dx = p.x - x, dy = p.y - y;
199  return dx*dx + dy*dy;
200  }
201  };
202 
203  using tree = QuadTree<point>;
204  using rectangle = typename tree::rectangle;
205  using point_list = typename tree::point_list;
206 
207 public:
208 
217  SurfaceSimplifier() = default;
218 
223  SurfaceSimplifier( double tolerance )
224  : m_tolerance( Abs( tolerance ) )
225  {
226  PCL_PRECONDITION( tolerance >= 0 )
227  }
228 
233 
240  double Tolerance() const
241  {
242  return m_tolerance;
243  }
244 
261  void SetTolerance( double tolerance )
262  {
263  PCL_PRECONDITION( tolerance > 0 )
264  m_tolerance = Abs( tolerance );
265  }
266 
273  bool IsRejectionEnabled() const
274  {
275  return m_enableRejection;
276  }
277 
288  void EnableRejection( bool enabled = true )
289  {
290  m_enableRejection = enabled;
291  }
292 
297  void DisableRejection( bool disable = true )
298  {
299  EnableRejection( !disable );
300  }
301 
309  float RejectFraction() const
310  {
311  return m_rejectFraction;
312  }
313 
332  void SetRejectFraction( float rejectFraction )
333  {
334  PCL_PRECONDITION( rejectFraction > 0 && rejectFraction < 1 )
335  m_rejectFraction = Range( rejectFraction, 0.0F, 1.0F );
336  }
337 
345  {
346  return m_includeCentroids;
347  }
348 
360  void EnableCentroidInclusion( bool enable = true )
361  {
362  m_includeCentroids = enable;
363  }
364 
369  void DisableCentroidInclusion( bool disable = true )
370  {
371  EnableCentroidInclusion( !disable );
372  }
373 
400  template <class C>
401  void Simplify( C& xs, C& ys, C& zs, const C& x, const C& y, const C& z ) const
402  {
403  int n = int( Min( Min( x.Length(), y.Length() ), z.Length() ) );
404  if ( n < 4 )
405  {
406  xs = x; ys = y; zs = z;
407  return;
408  }
409 
410  point_list P;
411  for ( int i = 0; i < n; ++i )
412  P << point( x[i], y[i], z[i] );
413 
414  tree R( P, n );
415  P = Simplify( R );
416  if ( int( P.Length() ) >= n )
417  {
418  xs = x; ys = y; zs = z;
419  return;
420  }
421 
422  R.Build( P, n );
423  P = Simplify( R );
424 
425  n = int( P.Length() );
426  xs = C( n );
427  ys = C( n );
428  zs = C( n );
429  for ( int i = 0; i < n; ++i )
430  {
431  xs[i] = typename C::item_type( P[i].x );
432  ys[i] = typename C::item_type( P[i].y );
433  zs[i] = typename C::item_type( P[i].z );
434  }
435  }
436 
437 private:
438 
439  double m_tolerance = __PCL_SURFACE_SIMPLIFIER_DEFAULT_TOLERANCE;
440  float m_rejectFraction = __PCL_SURFACE_SIMPLIFIER_DEFAULT_REJECT_FRACTION;
441  bool m_enableRejection = __PCL_SURFACE_SIMPLIFIER_DEFAULT_REJECTION_ENABLED;
442  bool m_includeCentroids = __PCL_SURFACE_SIMPLIFIER_DEFAULT_INCLUDE_CENTROIDS;
443 
450  point_list Simplify( tree& R ) const
451  {
452  point_list Q;
453  R.Traverse(
454  [&Q,this]( const rectangle& rect, point_list& points, void*& )
455  {
456  int n = int( points.Length() );
457  if ( n < 4 ) // cannot simplify triangles
458  {
459  Q << points;
460  return;
461  }
462 
463  /*
464  * Compute local centroid coordinates.
465  */
466  point p0( 0 );
467  for ( const point& p : points )
468  p0 += p;
469  p0 /= n;
470 
471  /*
472  * Form the covariance matrix.
473  */
474  double xx = 0, xy = 0, xz = 0, yy = 0, yz = 0, zz = 0;
475  for ( const point& p : points )
476  {
477  double dx = p.x - p0.x;
478  double dy = p.y - p0.y;
479  double dz = p.z - p0.z;
480  xx += dx*dx;
481  xy += dx*dy;
482  xz += dx*dz;
483  yy += dy*dy;
484  yz += dy*dz;
485  zz += dz*dz;
486  }
487  int n1 = n - 1;
488  xx /= n1;
489  xy /= n1;
490  xz /= n1;
491  yy /= n1;
492  yz /= n1;
493  zz /= n1;
494  Matrix M( xx, xy, xz,
495  xy, yy, yz,
496  xz, yz, zz );
497 
498  /*
499  * Test all local function values against the plane fitted at the
500  * centroid, with optional outlier rejection. The plane normal
501  * vector is the least eigenvector, from which we can form the
502  * plane equation: ax + by + cz = 0.
503  */
504  ComputeEigenvectors( M );
505  double a = M[0][0];
506  double b = M[1][0];
507  double c = M[2][0];
508  int mr = m_enableRejection ? Max( 1, TruncInt( m_rejectFraction*n ) ) : 1; // max. number of outliers
509  int nr = 0;
510  point_list P = points;
511  for ( point& p : P )
512  {
513  double z = p0.z - (a*(p.x - p0.x) + b*(p.y - p0.y))/c;
514  if ( Abs( p.z - z ) > m_tolerance )
515  {
516  if ( ++nr == mr )
517  {
518  /*
519  * If the current region deviates from the fitted plane
520  * more than allowed by the tolerance parameter after
521  * outlier rejection, try to simplify it further with a
522  * deeper quadtree subdivision.
523  */
524  tree R( points, TruncInt( 0.75*n ) );
525  Q << Simplify( R );
526  return;
527  }
528 
529  // Winsorize outlier function values.
530  p.z = (p.z > z) ? z + m_tolerance : z - m_tolerance;
531  }
532  }
533 
534  /*
535  * If the current region is flat to within the tolerance parameter,
536  * take the convex hull as its simplified point set. This is what
537  * makes our simplification algorithm shape-preserving.
538  */
539  if ( m_includeCentroids )
540  {
541  // If one or more points have been rejected, calculate a new
542  // centroid function value from Winsorized z coordinates.
543  if ( nr > 0 )
544  {
545  p0.z = 0;
546  for ( const point& p : P )
547  p0.z += p.z;
548  p0.z /= n;
549  }
550  Q << p0;
551  }
552  Q << ConvexHull( P );
553  }
554  );
555 
556  return Q;
557  }
558 
563  static void ComputeEigenvectors( Matrix& );
564 
569  static point_list ConvexHull( point_list& );
570 };
571 
572 // ----------------------------------------------------------------------------
573 
574 } // pcl
575 
576 #endif // __PCL_SurfaceSimplifier_h
577 
578 // ----------------------------------------------------------------------------
579 // EOF pcl/SurfaceSimplifier.h - Released 2024-06-18T15:48:54Z
A generic rectangle in the two-dimensional space.
Definition: Rectangle.h:314
64-bit floating point real matrix.
Bucket PR quadtree for two-dimensional point data.
Definition: QuadTree.h:124
Shape-preserving simplification of 2-D surfaces.
bool IsCentroidInclusionEnabled() const
void Simplify(C &xs, C &ys, C &zs, const C &x, const C &y, const C &z) const
void EnableCentroidInclusion(bool enable=true)
void DisableRejection(bool disable=true)
SurfaceSimplifier(const SurfaceSimplifier &)=default
SurfaceSimplifier(double tolerance)
void EnableRejection(bool enabled=true)
void DisableCentroidInclusion(bool disable=true)
void SetRejectFraction(float rejectFraction)
void SetTolerance(double tolerance)
T Abs(const Complex< T > &c) noexcept
Definition: Complex.h:429
int TruncInt(T x) noexcept
Definition: Math.h:1132
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