PCL
SurfaceSimplifier.h
Go to the documentation of this file.
1 // ____ ______ __
2 // / __ \ / ____// /
3 // / /_/ // / / /
4 // / ____// /___ / /___ PixInsight Class Library
5 // /_/ \____//_____/ PCL 2.4.7
6 // ----------------------------------------------------------------------------
7 // pcl/SurfaceSimplifier.h - Released 2020-12-17T15:46:28Z
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-2020 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  typedef double component;
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  typedef QuadTree<point> tree;
204 
205  typedef typename tree::rectangle rectangle;
206 
207  typedef typename tree::point_list point_list;
208 
209 public:
210 
219  SurfaceSimplifier() = default;
220 
227  SurfaceSimplifier( double tolerance )
228  : m_tolerance( Abs( tolerance ) )
229  {
230  PCL_PRECONDITION( tolerance >= 0 )
231  }
232 
236  SurfaceSimplifier( const SurfaceSimplifier& ) = default;
237 
244  double Tolerance() const
245  {
246  return m_tolerance;
247  }
248 
265  void SetTolerance( double tolerance )
266  {
267  PCL_PRECONDITION( tolerance > 0 )
268  m_tolerance = Abs( tolerance );
269  }
270 
277  bool IsRejectionEnabled() const
278  {
279  return m_enableRejection;
280  }
281 
292  void EnableRejection( bool enabled = true )
293  {
294  m_enableRejection = enabled;
295  }
296 
301  void DisableRejection( bool disable = true )
302  {
303  EnableRejection( !disable );
304  }
305 
313  float RejectFraction() const
314  {
315  return m_rejectFraction;
316  }
317 
336  void SetRejectFraction( float rejectFraction )
337  {
338  PCL_PRECONDITION( rejectFraction > 0 && rejectFraction < 1 )
339  m_rejectFraction = Range( rejectFraction, 0.0F, 1.0F );
340  }
341 
349  {
350  return m_includeCentroids;
351  }
352 
364  void EnableCentroidInclusion( bool enable = true )
365  {
366  m_includeCentroids = enable;
367  }
368 
373  void DisableCentroidInclusion( bool disable = true )
374  {
375  EnableCentroidInclusion( !disable );
376  }
377 
404  template <class C>
405  void Simplify( C& xs, C& ys, C& zs, const C& x, const C& y, const C& z ) const
406  {
407  int n = Min( Min( x.Length(), y.Length() ), z.Length() );
408  if ( n < 4 )
409  {
410  xs = x; ys = y; zs = z;
411  return;
412  }
413 
414  point_list P;
415  for ( int i = 0; i < n; ++i )
416  P << point( x[i], y[i], z[i] );
417 
418  tree R( P, n );
419  P = Simplify( R );
420  if ( int( P.Length() ) >= n )
421  {
422  xs = x; ys = y; zs = z;
423  return;
424  }
425 
426  R.Build( P, n );
427  P = Simplify( R );
428 
429  n = int( P.Length() );
430  xs = C( n );
431  ys = C( n );
432  zs = C( n );
433  for ( int i = 0; i < n; ++i )
434  {
435  xs[i] = typename C::item_type( P[i].x );
436  ys[i] = typename C::item_type( P[i].y );
437  zs[i] = typename C::item_type( P[i].z );
438  }
439  }
440 
441 private:
442 
443  double m_tolerance = __PCL_SURFACE_SIMPLIFIER_DEFAULT_TOLERANCE;
444  float m_rejectFraction = __PCL_SURFACE_SIMPLIFIER_DEFAULT_REJECT_FRACTION;
445  bool m_enableRejection = __PCL_SURFACE_SIMPLIFIER_DEFAULT_REJECTION_ENABLED;
446  bool m_includeCentroids = __PCL_SURFACE_SIMPLIFIER_DEFAULT_INCLUDE_CENTROIDS;
447 
454  point_list Simplify( tree& R ) const
455  {
456  point_list Q;
457  R.Traverse(
458  [&Q,this]( const rectangle& rect, point_list& points, void*& )
459  {
460  int n = int( points.Length() );
461  if ( n < 4 ) // cannot simplify triangles
462  {
463  Q << points;
464  return;
465  }
466 
467  /*
468  * Compute local centroid coordinates.
469  */
470  point p0( 0 );
471  for ( const point& p : points )
472  p0 += p;
473  p0 /= n;
474 
475  /*
476  * Form the covariance matrix.
477  */
478  double xx = 0, xy = 0, xz = 0, yy = 0, yz = 0, zz = 0;
479  for ( const point& p : points )
480  {
481  double dx = p.x - p0.x;
482  double dy = p.y - p0.y;
483  double dz = p.z - p0.z;
484  xx += dx*dx;
485  xy += dx*dy;
486  xz += dx*dz;
487  yy += dy*dy;
488  yz += dy*dz;
489  zz += dz*dz;
490  }
491  int n1 = n - 1;
492  xx /= n1;
493  xy /= n1;
494  xz /= n1;
495  yy /= n1;
496  yz /= n1;
497  zz /= n1;
498  Matrix M( xx, xy, xz,
499  xy, yy, yz,
500  xz, yz, zz );
501 
502  /*
503  * Test all local function values against the plane fitted at the
504  * centroid, with optional outlier rejection. The plane normal
505  * vector is the least eigenvector, from which we can form the
506  * plane equation: ax + by + cz = 0.
507  */
508  ComputeEigenvectors( M );
509  double a = M[0][0];
510  double b = M[1][0];
511  double c = M[2][0];
512  int mr = m_enableRejection ? Max( 1, TruncInt( m_rejectFraction*n ) ) : 1; // max. number of outliers
513  int nr = 0;
514  point_list P = points;
515  for ( point& p : P )
516  {
517  double z = p0.z - (a*(p.x - p0.x) + b*(p.y - p0.y))/c;
518  if ( Abs( p.z - z ) > m_tolerance )
519  {
520  if ( ++nr == mr )
521  {
522  /*
523  * If the current region deviates from the fitted plane
524  * more than allowed by the tolerance parameter after
525  * outlier rejection, try to simplify it further with a
526  * deeper quadtree subdivision.
527  */
528  tree R( points, TruncInt( 0.75*n ) );
529  Q << Simplify( R );
530  return;
531  }
532 
533  // Winsorize outlier function values.
534  p.z = (p.z > z) ? z + m_tolerance : z - m_tolerance;
535  }
536  }
537 
538  /*
539  * If the current region is flat to within the tolerance parameter,
540  * take the convex hull as its simplified point set. This is what
541  * makes our simplification algorithm shape-preserving.
542  */
543  if ( m_includeCentroids )
544  {
545  // If one or more points have been rejected, calculate a new
546  // centroid function value from Winsorized z coordinates.
547  if ( nr > 0 )
548  {
549  p0.z = 0;
550  for ( const point& p : P )
551  p0.z += p.z;
552  p0.z /= n;
553  }
554  Q << p0;
555  }
556  Q << ConvexHull( P );
557  }
558  );
559 
560  return Q;
561  }
562 
567  static void ComputeEigenvectors( Matrix& );
568 
573  static point_list ConvexHull( point_list& );
574 };
575 
576 // ----------------------------------------------------------------------------
577 
578 } // pcl
579 
580 #endif // __PCL_SurfaceSimplifier_h
581 
582 // ----------------------------------------------------------------------------
583 // EOF pcl/SurfaceSimplifier.h - Released 2020-12-17T15:46:28Z
void EnableRejection(bool enabled=true)
void DisableRejection(bool disable=true)
void EnableCentroidInclusion(bool enable=true)
void Simplify(C &xs, C &ys, C &zs, const C &x, const C &y, const C &z) const
PCL root namespace.
Definition: AbstractImage.h:76
void Build(const point_list &points, int bucketCapacity=40)
Definition: QuadTree.h:430
int TruncInt(T x) noexcept
Definition: Math.h:1132
constexpr const T & Max(const T &a, const T &b) noexcept
Definition: Utility.h:119
static void Traverse(const Node *node, F f)
Definition: QuadTree.h:709
64-bit floating-point rectangle in the R^2 space.
void SetTolerance(double tolerance)
64-bit floating point real matrix.
void DisableCentroidInclusion(bool disable=true)
size_type Length() const noexcept
Definition: Array.h:268
constexpr const T & Min(const T &a, const T &b) noexcept
Definition: Utility.h:90
T Abs(const Complex< T > &c) noexcept
Definition: Complex.h:429
Bucket PR quadtree for two-dimensional point data.
Definition: QuadTree.h:111
SurfaceSimplifier(double tolerance)
constexpr const T & Range(const T &x, const T &a, const T &b) noexcept
Definition: Utility.h:190
bool IsCentroidInclusionEnabled() const
void SetRejectFraction(float rejectFraction)