SurfaceSimplifier.h

Go to the documentation of this file.

//     ____   ______ __
//    / __ \ / ____// /
//   / /_/ // /    / /
//  / ____// /___ / /___   PixInsight Class Library
// /_/     \____//_____/   PCL 2.4.7
// ----------------------------------------------------------------------------
// pcl/SurfaceSimplifier.h - Released 2020-12-17T15:46:28Z
// ----------------------------------------------------------------------------
// This file is part of the PixInsight Class Library (PCL).
// PCL is a multiplatform C++ framework for development of PixInsight modules.
//
// Copyright (c) 2003-2020 Pleiades Astrophoto S.L. All Rights Reserved.
//
// Redistribution and use in both source and binary forms, with or without
// modification, is permitted provided that the following conditions are met:
//
// 1. All redistributions of source code must retain the above copyright
//    notice, this list of conditions and the following disclaimer.
//
// 2. All redistributions in binary form must reproduce the above copyright
//    notice, this list of conditions and the following disclaimer in the
//    documentation and/or other materials provided with the distribution.
//
// 3. Neither the names "PixInsight" and "Pleiades Astrophoto", nor the names
//    of their contributors, may be used to endorse or promote products derived
//    from this software without specific prior written permission. For written
//    permission, please contact info@pixinsight.com.
//
// 4. All products derived from this software, in any form whatsoever, must
//    reproduce the following acknowledgment in the end-user documentation
//    and/or other materials provided with the product:
//
//    "This product is based on software from the PixInsight project, developed
//    by Pleiades Astrophoto and its contributors (https://pixinsight.com/)."
//
//    Alternatively, if that is where third-party acknowledgments normally
//    appear, this acknowledgment must be reproduced in the product itself.
//
// THIS SOFTWARE IS PROVIDED BY PLEIADES ASTROPHOTO AND ITS CONTRIBUTORS
// "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
// TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL PLEIADES ASTROPHOTO OR ITS
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, BUSINESS
// INTERRUPTION; PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; AND LOSS OF USE,
// DATA OR PROFITS) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
// POSSIBILITY OF SUCH DAMAGE.
// ----------------------------------------------------------------------------

#ifndef __PCL_SurfaceSimplifier_h
#define __PCL_SurfaceSimplifier_h


#include <pcl/Defs.h>
#include <pcl/Diagnostics.h>

#include <pcl/Matrix.h>
#include <pcl/QuadTree.h>

namespace pcl
{

// ----------------------------------------------------------------------------

/*
 * Default tolerance for surface simplification.
 *
 * The tolerance parameter determines the maximum absolute deviation of
 * function values from a plane allowed to simplify the subset of points within
 * a subregion of the input point space. This parameter is specified in
 * bivariate function values, i.e. in Z-axis units.
 */
#define __PCL_SURFACE_SIMPLIFIER_DEFAULT_TOLERANCE          0.01

/*
 * Default state of the surface simplification outlier rejection feature.
 *
 * When enabled, a prescribed fraction of outlier points will be rejected on
 * each subregion for estimation of local curvature.
 */
#define __PCL_SURFACE_SIMPLIFIER_DEFAULT_REJECTION_ENABLED  true

/*
 * Default rejection fraction for surface simplification.
 *
 * This parameter specifies the maximum fraction of outliers allowed for
 * simplification of a subset of points in a region of the input point space.
 * Point rejection makes the surface simplification algorithm robust to outlier
 * function values.
 */
#define __PCL_SURFACE_SIMPLIFIER_DEFAULT_REJECT_FRACTION    0.2F

/*
 * Whether to include subregion centroids in simplified point sets.
 *
 * When enabled, the centroid of each simplified subregion will also be
 * included in the corresponding list of simplified points. This can improve
 * the shape preservation behavior of the surface simplification algorithm, at
 * the cost of a reduced amount of additional points.
 */
#define __PCL_SURFACE_SIMPLIFIER_DEFAULT_INCLUDE_CENTROIDS  false

// ----------------------------------------------------------------------------

class PCL_CLASS SurfaceSimplifier
{
private:

   struct point
   {
      typedef double component;

      component x, y, z;

      point() = default;
      point( const point& ) = default;

      template <typename T>
      point( T a_x, T a_y, T a_z )
         : x( double( a_x ) )
         , y( double( a_y ) )
         , z( double( a_z ) )
      {
      }

      template <typename T>
      point( T k )
      {
         x = y = z = double( k );
      }

      component operator[]( int i ) const
      {
         return i ? y : x;
      }

      point& operator +=( const point& p )
      {
         x += p.x; y += p.y; z += p.z;
         return *this;
      }

      template <typename T>
      point& operator /=( T k )
      {
         x /= k; y /= k; z /= k;
         return *this;
      }

      double SquaredDistanceTo( const point& p ) const
      {
         double dx = p.x - x, dy = p.y - y;
         return dx*dx + dy*dy;
      }
   };

   typedef QuadTree<point>             tree;

   typedef typename tree::rectangle    rectangle;

   typedef typename tree::point_list   point_list;

public:

   SurfaceSimplifier() = default;

   SurfaceSimplifier( double tolerance )
      : m_tolerance( Abs( tolerance ) )
   {
      PCL_PRECONDITION( tolerance >= 0 )
   }

   SurfaceSimplifier( const SurfaceSimplifier& ) = default;

   double Tolerance() const
   {
      return m_tolerance;
   }

   void SetTolerance( double tolerance )
   {
      PCL_PRECONDITION( tolerance > 0 )
      m_tolerance = Abs( tolerance );
   }

   bool IsRejectionEnabled() const
   {
      return m_enableRejection;
   }

   void EnableRejection( bool enabled = true )
   {
      m_enableRejection = enabled;
   }

   void DisableRejection( bool disable = true )
   {
      EnableRejection( !disable );
   }

   float RejectFraction() const
   {
      return m_rejectFraction;
   }

   void SetRejectFraction( float rejectFraction )
   {
      PCL_PRECONDITION( rejectFraction > 0 && rejectFraction < 1 )
      m_rejectFraction = Range( rejectFraction, 0.0F, 1.0F );
   }

   bool IsCentroidInclusionEnabled() const
   {
      return m_includeCentroids;
   }

   void EnableCentroidInclusion( bool enable = true )
   {
      m_includeCentroids = enable;
   }

   void DisableCentroidInclusion( bool disable = true )
   {
      EnableCentroidInclusion( !disable );
   }

   template <class C>
   void Simplify( C& xs, C& ys, C& zs, const C& x, const C& y, const C& z ) const
   {
      int n = Min( Min( x.Length(), y.Length() ), z.Length() );
      if ( n < 4 )
      {
         xs = x; ys = y; zs = z;
         return;
      }

      point_list P;
      for ( int i = 0; i < n; ++i )
         P << point( x[i], y[i], z[i] );

      tree R( P, n );
      P = Simplify( R );
      if ( int( P.Length() ) >= n )
      {
         xs = x; ys = y; zs = z;
         return;
      }

      R.Build( P, n );
      P = Simplify( R );

      n = int( P.Length() );
      xs = C( n );
      ys = C( n );
      zs = C( n );
      for ( int i = 0; i < n; ++i )
      {
         xs[i] = typename C::item_type( P[i].x );
         ys[i] = typename C::item_type( P[i].y );
         zs[i] = typename C::item_type( P[i].z );
      }
   }

private:

   double m_tolerance = __PCL_SURFACE_SIMPLIFIER_DEFAULT_TOLERANCE;
   float  m_rejectFraction = __PCL_SURFACE_SIMPLIFIER_DEFAULT_REJECT_FRACTION;
   bool   m_enableRejection = __PCL_SURFACE_SIMPLIFIER_DEFAULT_REJECTION_ENABLED;
   bool   m_includeCentroids = __PCL_SURFACE_SIMPLIFIER_DEFAULT_INCLUDE_CENTROIDS;

   point_list Simplify( tree& R ) const
   {
      point_list Q;
      R.Traverse(
         [&Q,this]( const rectangle& rect, point_list& points, void*& )
         {
            int n = int( points.Length() );
            if ( n < 4 ) // cannot simplify triangles
            {
               Q << points;
               return;
            }

            /*
             * Compute local centroid coordinates.
             */
            point p0( 0 );
            for ( const point& p : points )
               p0 += p;
            p0 /= n;

            /*
             * Form the covariance matrix.
             */
            double xx = 0, xy = 0, xz = 0, yy = 0, yz = 0, zz = 0;
            for ( const point& p : points )
            {
               double dx = p.x - p0.x;
               double dy = p.y - p0.y;
               double dz = p.z - p0.z;
               xx += dx*dx;
               xy += dx*dy;
               xz += dx*dz;
               yy += dy*dy;
               yz += dy*dz;
               zz += dz*dz;
            }
            int n1 = n - 1;
            xx /= n1;
            xy /= n1;
            xz /= n1;
            yy /= n1;
            yz /= n1;
            zz /= n1;
            Matrix M( xx, xy, xz,
                      xy, yy, yz,
                      xz, yz, zz );

            /*
             * Test all local function values against the plane fitted at the
             * centroid, with optional outlier rejection. The plane normal
             * vector is the least eigenvector, from which we can form the
             * plane equation: ax + by + cz = 0.
             */
            ComputeEigenvectors( M );
            double a = M[0][0];
            double b = M[1][0];
            double c = M[2][0];
            int mr = m_enableRejection ? Max( 1, TruncInt( m_rejectFraction*n ) ) : 1; // max. number of outliers
            int nr = 0;
            point_list P = points;
            for ( point& p : P )
            {
               double z = p0.z - (a*(p.x - p0.x) + b*(p.y - p0.y))/c;
               if ( Abs( p.z - z ) > m_tolerance )
               {
                  if ( ++nr == mr )
                  {
                     /*
                      * If the current region deviates from the fitted plane
                      * more than allowed by the tolerance parameter after
                      * outlier rejection, try to simplify it further with a
                      * deeper quadtree subdivision.
                      */
                     tree R( points, TruncInt( 0.75*n ) );
                     Q << Simplify( R );
                     return;
                  }

                  // Winsorize outlier function values.
                  p.z = (p.z > z) ? z + m_tolerance : z - m_tolerance;
               }
            }

            /*
             * If the current region is flat to within the tolerance parameter,
             * take the convex hull as its simplified point set. This is what
             * makes our simplification algorithm shape-preserving.
             */
            if ( m_includeCentroids )
            {
               // If one or more points have been rejected, calculate a new
               // centroid function value from Winsorized z coordinates.
               if ( nr > 0 )
               {
                  p0.z = 0;
                  for ( const point& p : P )
                     p0.z += p.z;
                  p0.z /= n;
               }
               Q << p0;
            }
            Q << ConvexHull( P );
         }
      );

      return Q;
   }

   static void ComputeEigenvectors( Matrix& );

   static point_list ConvexHull( point_list& );
};

// ----------------------------------------------------------------------------

}  // pcl

#endif   // __PCL_SurfaceSimplifier_h

// ----------------------------------------------------------------------------
// EOF pcl/SurfaceSimplifier.h - Released 2020-12-17T15:46:28Z