PCL
SurfaceSimplifier.h
Go to the documentation of this file.
1 // ____ ______ __
2 // / __ \ / ____// /
3 // / /_/ // / / /
4 // / ____// /___ / /___ PixInsight Class Library
5 // /_/ \____//_____/ PCL 2.1.19
6 // ----------------------------------------------------------------------------
7 // pcl/SurfaceSimplifier.h - Released 2019-11-07T10:59:34Z
8 // ----------------------------------------------------------------------------
9 // This file is part of the PixInsight Class Library (PCL).
10 // PCL is a multiplatform C++ framework for development of PixInsight modules.
11 //
12 // Copyright (c) 2003-2019 Pleiades Astrophoto S.L. All Rights Reserved.
13 //
14 // Redistribution and use in both source and binary forms, with or without
15 // modification, is permitted provided that the following conditions are met:
16 //
17 // 1. All redistributions of source code must retain the above copyright
18 // notice, this list of conditions and the following disclaimer.
19 //
20 // 2. All redistributions in binary form must reproduce the above copyright
21 // notice, this list of conditions and the following disclaimer in the
22 // documentation and/or other materials provided with the distribution.
23 //
24 // 3. Neither the names "PixInsight" and "Pleiades Astrophoto", nor the names
25 // of their contributors, may be used to endorse or promote products derived
26 // from this software without specific prior written permission. For written
27 // permission, please contact info@pixinsight.com.
28 //
29 // 4. All products derived from this software, in any form whatsoever, must
30 // reproduce the following acknowledgment in the end-user documentation
31 // and/or other materials provided with the product:
32 //
33 // "This product is based on software from the PixInsight project, developed
34 // by Pleiades Astrophoto and its contributors (http://pixinsight.com/)."
35 //
36 // Alternatively, if that is where third-party acknowledgments normally
37 // appear, this acknowledgment must be reproduced in the product itself.
38 //
39 // THIS SOFTWARE IS PROVIDED BY PLEIADES ASTROPHOTO AND ITS CONTRIBUTORS
40 // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
41 // TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
42 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL PLEIADES ASTROPHOTO OR ITS
43 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
44 // EXEMPLARY OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, BUSINESS
45 // INTERRUPTION; PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; AND LOSS OF USE,
46 // DATA OR PROFITS) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
47 // CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
48 // ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
49 // POSSIBILITY OF SUCH DAMAGE.
50 // ----------------------------------------------------------------------------
51 
52 #ifndef __PCL_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 ) ), y( double( a_y ) ), z( double( a_z ) )
167  {
168  }
169 
170  template <typename T>
171  point( T k )
172  {
173  x = y = z = double( k );
174  }
175 
176  component operator[]( int i ) const
177  {
178  return i ? y : x;
179  }
180 
181  point& operator +=( const point& p )
182  {
183  x += p.x; y += p.y; z += p.z;
184  return *this;
185  }
186 
187  template <typename T>
188  point& operator /=( T k )
189  {
190  x /= k; y /= k; z /= k;
191  return *this;
192  }
193 
194  double SquaredDistanceTo( const point& p ) const
195  {
196  double dx = p.x - x, dy = p.y - y;
197  return dx*dx + dy*dy;
198  }
199  };
200 
201  typedef QuadTree<point> tree;
202 
203  typedef typename tree::rectangle rectangle;
204 
205  typedef typename tree::point_list point_list;
206 
207 public:
208 
217  SurfaceSimplifier() = default;
218 
225  SurfaceSimplifier( double tolerance ) :
226  m_tolerance( Abs( tolerance ) )
227  {
228  PCL_PRECONDITION( tolerance >= 0 )
229  }
230 
234  SurfaceSimplifier( const SurfaceSimplifier& ) = default;
235 
242  double Tolerance() const
243  {
244  return m_tolerance;
245  }
246 
263  void SetTolerance( double tolerance )
264  {
265  PCL_PRECONDITION( tolerance > 0 )
266  m_tolerance = Abs( tolerance );
267  }
268 
275  bool IsRejectionEnabled() const
276  {
277  return m_enableRejection;
278  }
279 
290  void EnableRejection( bool enabled = true )
291  {
292  m_enableRejection = enabled;
293  }
294 
299  void DisableRejection( bool disable = true )
300  {
301  EnableRejection( !disable );
302  }
303 
311  float RejectFraction() const
312  {
313  return m_rejectFraction;
314  }
315 
334  void SetRejectFraction( float rejectFraction )
335  {
336  PCL_PRECONDITION( rejectFraction > 0 && rejectFraction < 1 )
337  m_rejectFraction = Range( rejectFraction, 0.0F, 1.0F );
338  }
339 
347  {
348  return m_includeCentroids;
349  }
350 
362  void EnableCentroidInclusion( bool enable = true )
363  {
364  m_includeCentroids = enable;
365  }
366 
371  void DisableCentroidInclusion( bool disable = true )
372  {
373  EnableCentroidInclusion( !disable );
374  }
375 
402  template <class C>
403  void Simplify( C& xs, C& ys, C& zs, const C& x, const C& y, const C& z ) const
404  {
405  int n = Min( Min( x.Length(), y.Length() ), z.Length() );
406  if ( n < 4 )
407  {
408  xs = x; ys = y; zs = z;
409  return;
410  }
411 
412  point_list P;
413  for ( int i = 0; i < n; ++i )
414  P << point( x[i], y[i], z[i] );
415 
416  tree R( P, n );
417  P = Simplify( R );
418  if ( int( P.Length() ) >= n )
419  {
420  xs = x; ys = y; zs = z;
421  return;
422  }
423 
424  R.Build( P, n );
425  P = Simplify( R );
426 
427  n = int( P.Length() );
428  xs = C( n );
429  ys = C( n );
430  zs = C( n );
431  for ( int i = 0; i < n; ++i )
432  {
433  xs[i] = typename C::item_type( P[i].x );
434  ys[i] = typename C::item_type( P[i].y );
435  zs[i] = typename C::item_type( P[i].z );
436  }
437  }
438 
439 private:
440 
441  double m_tolerance = __PCL_SURFACE_SIMPLIFIER_DEFAULT_TOLERANCE;
442  float m_rejectFraction = __PCL_SURFACE_SIMPLIFIER_DEFAULT_REJECT_FRACTION;
443  bool m_enableRejection = __PCL_SURFACE_SIMPLIFIER_DEFAULT_REJECTION_ENABLED;
444  bool m_includeCentroids = __PCL_SURFACE_SIMPLIFIER_DEFAULT_INCLUDE_CENTROIDS;
445 
452  point_list Simplify( tree& R ) const
453  {
454  point_list Q;
455  R.Traverse(
456  [&Q,this]( const rectangle& rect, point_list& points, void*& )
457  {
458  int n = int( points.Length() );
459  if ( n < 4 ) // cannot simplify triangles
460  {
461  Q << points;
462  return;
463  }
464 
465  /*
466  * Compute local centroid coordinates.
467  */
468  point p0( 0 );
469  for ( const point& p : points )
470  p0 += p;
471  p0 /= n;
472 
473  /*
474  * Form the covariance matrix.
475  */
476  double xx = 0, xy = 0, xz = 0, yy = 0, yz = 0, zz = 0;
477  for ( const point& p : points )
478  {
479  double dx = p.x - p0.x;
480  double dy = p.y - p0.y;
481  double dz = p.z - p0.z;
482  xx += dx*dx;
483  xy += dx*dy;
484  xz += dx*dz;
485  yy += dy*dy;
486  yz += dy*dz;
487  zz += dz*dz;
488  }
489  int n1 = n - 1;
490  xx /= n1;
491  xy /= n1;
492  xz /= n1;
493  yy /= n1;
494  yz /= n1;
495  zz /= n1;
496  Matrix M( xx, xy, xz,
497  xy, yy, yz,
498  xz, yz, zz );
499 
500  /*
501  * Test all local function values against the plane fitted at the
502  * centroid, with optional outlier rejection. The plane normal
503  * vector is the least eigenvector, from which we can form the
504  * plane equation: ax + by + cz = 0.
505  */
506  ComputeEigenvectors( M );
507  double a = M[0][0];
508  double b = M[1][0];
509  double c = M[2][0];
510  int mr = m_enableRejection ? Max( 1, TruncInt( m_rejectFraction*n ) ) : 1; // max. number of outliers
511  int nr = 0;
512  point_list P = points;
513  for ( point& p : P )
514  {
515  double z = p0.z - (a*(p.x - p0.x) + b*(p.y - p0.y))/c;
516  if ( Abs( p.z - z ) > m_tolerance )
517  {
518  if ( ++nr == mr )
519  {
520  /*
521  * If the current region deviates from the fitted plane
522  * more than allowed by the tolerance parameter after
523  * outlier rejection, try to simplify it further with a
524  * deeper quadtree subdivision.
525  */
526  tree R( points, TruncInt( 0.75*n ) );
527  Q << Simplify( R );
528  return;
529  }
530 
531  // Winsorize outlier function values.
532  p.z = (p.z > z) ? z + m_tolerance : z - m_tolerance;
533  }
534  }
535 
536  /*
537  * If the current region is flat to within the tolerance parameter,
538  * take the convex hull as its simplified point set. This is what
539  * makes our simplification algorithm shape-preserving.
540  */
541  if ( m_includeCentroids )
542  {
543  // If one or more points have been rejected, calculate a new
544  // centroid function value from Winsorized z coordinates.
545  if ( nr > 0 )
546  {
547  p0.z = 0;
548  for ( const point& p : P )
549  p0.z += p.z;
550  p0.z /= n;
551  }
552  Q << p0;
553  }
554  Q << ConvexHull( P );
555  }
556  );
557 
558  return Q;
559  }
560 
565  static void ComputeEigenvectors( Matrix& );
566 
571  static point_list ConvexHull( point_list& );
572 };
573 
574 // ----------------------------------------------------------------------------
575 
576 } // pcl
577 
578 #endif // __PCL_SurfaceSimplifier_h
579 
580 // ----------------------------------------------------------------------------
581 // EOF pcl/SurfaceSimplifier.h - Released 2019-11-07T10:59:34Z
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:402
static void Traverse(const Node *node, F f)
Definition: QuadTree.h:634
constexpr const T & Range(const T &x, const T &a, const T &b)
Definition: Utility.h:190
64-bit floating-point rectangle in the R^2 space.
void SetTolerance(double tolerance)
64-bit floating point real matrix.
T Abs(const Complex< T > &c)
Definition: Complex.h:420
void DisableCentroidInclusion(bool disable=true)
constexpr const T & Max(const T &a, const T &b)
Definition: Utility.h:119
int TruncInt(T x)
Definition: Math.h:1035
constexpr const T & Min(const T &a, const T &b)
Definition: Utility.h:90
size_type Length() const
Definition: Array.h:265
Bucket PR quadtree for two-dimensional point data.
Definition: QuadTree.h:111
SurfaceSimplifier(double tolerance)
bool IsCentroidInclusionEnabled() const
void SetRejectFraction(float rejectFraction)