Added finding outlier points

source/MRMesh/MRMesh.vcxproj

@@ -144,6 +144,7 @@
     <ClInclude Include="MRNormalDenoising.h" />
     <ClInclude Include="MRNormalsToPoints.h" />
     <ClInclude Include="MROffsetContours.h" />
+    <ClInclude Include="MROutlierPoints.h" />
     <ClInclude Include="MRParabola.h" />
     <ClInclude Include="MRBestFitParabola.h" />
     <ClInclude Include="MRParallel.h" />
@@ -401,6 +402,7 @@
     <ClInclude Include="MRScopedValue.h" />
   </ItemGroup>
   <ItemGroup>
+    <ClCompile Include="MROutlierPoints.cpp" />
     <ClCompile Include="MRParallelProgressReporter.cpp" />
     <ClCompile Include="MR2DContoursTriangulation.cpp" />
     <ClCompile Include="MRAABBTreeMaker.cpp" />

source/MRMesh/MRMesh.vcxproj.filters

@@ -1272,6 +1272,9 @@
     <ClInclude Include="MRPointsInBox.h">
       <Filter>Source Files\AABBTree</Filter>
     </ClInclude>
+    <ClInclude Include="MROutlierPoints.h">
+      <Filter>Source Files\PointCloud</Filter>
+    </ClInclude>
   </ItemGroup>
   <ItemGroup>
     <ClCompile Include="MRParallelProgressReporter.cpp">
@@ -2093,6 +2096,9 @@
     <ClCompile Include="MRPointsInBox.cpp">
       <Filter>Source Files\AABBTree</Filter>
     </ClCompile>
+    <ClCompile Include="MROutlierPoints.cpp">
+      <Filter>Source Files\PointCloud</Filter>
+    </ClCompile>
   </ItemGroup>
   <ItemGroup>
     <None Include="..\.editorconfig" />

source/MRMesh/MROutlierPoints.cpp

@@ -0,0 +1,234 @@
+#include "MROutlierPoints.h"
+#include "MRPointCloud.h"
+#include "MRProgressCallback.h"
+#include "MRBitSetParallelFor.h"
+#include "MRBestFit.h"
+#include "MRPointsInBall.h"
+#include "MRPointsComponents.h"
+
+namespace MR
+{
+
+Expected<void> OutliersDetector::prepare( const PointCloud& pointCloud, float radius, OutlierTypeMask mask, ProgressCallback progress /*= {}*/ )
+{
+    validPoints_ = pointCloud.validPoints;
+
+    if ( !validPoints_.any() )
+        return unexpected( "Empty Point Cloud" );
+
+    if ( !bool( mask & OutlierTypeMask::All ) )
+        return unexpected( "Nothing to look for" );
+
+    maskCached_ = mask;
+    radius_ = radius;
+
+    const bool calcSmallComponentsCached = bool( mask & OutlierTypeMask::SmallComponents );
+    const bool calcWeaklyConnectedCached = bool( mask & OutlierTypeMask::WeaklyConnected );
+    const bool calcFarSurfaceCached = bool( mask & OutlierTypeMask::FarSurface );
+    const bool calcBadNormalCached = bool( mask & OutlierTypeMask::AwayNormal );
+
+
+    const VertBitSet* lastPassVerts = &validPoints_;
+    const auto numVerts = validPoints_.find_last() + 1;
+    const auto numThreads = int( tbb::global_control::active_value( tbb::global_control::max_allowed_parallelism ) );
+
+    if ( calcSmallComponentsCached )
+        unionFindStructure_ = UnionFind<VertId>( numVerts );
+    if ( calcWeaklyConnectedCached )
+        weaklyConnectedStat_ = std::vector<uint8_t>( numVerts );
+    if ( calcFarSurfaceCached )
+        farSurfaceStat_ = std::vector<float>( numVerts );
+    if ( calcBadNormalCached )
+        badNormalStat_ = std::vector<float>( numVerts );
+
+    VertBitSet secondPassVerts;
+    ProgressCallback subProgress = subprogress( progress, 0.f, 0.4f );
+    const auto& points = pointCloud.points;
+    const auto& normals = pointCloud.normals;
+    if ( numThreads > 1 )
+    {
+        secondPassVerts.resize( numVerts );
+        lastPassVerts = &secondPassVerts;
+        const bool continued = BitSetParallelForAllRanged( validPoints_, [&] ( VertId v0, const auto& range )
+        {
+            if ( !contains( validPoints_, v0 ) )
+                return;
+            int count = 0;
+            PointAccumulator plane;
+            Vector3f normalSum;
+            findPointsInBall( pointCloud.getAABBTree(), { pointCloud.points[v0], sqr( radius_ ) },
+                              [&] ( VertId v1, const Vector3f& )
+            {
+                if ( !contains( validPoints_, v1 ) )
+                    return;
+                if ( v1 != v0 )
+                {
+                    ++count;
+                    if ( calcBadNormalCached )
+                        normalSum += normals[v1];
+                    if ( calcFarSurfaceCached )
+                        plane.addPoint( points[v1] );
+                }
+                if ( calcSmallComponentsCached && v0 < v1 )
+                {
+                    if ( v1 >= range.end )
+                        secondPassVerts.set( v0 );
+                    else
+                        unionFindStructure_.unite( v0, v1 );
+                }
+            } );
+            if ( calcWeaklyConnectedCached )
+                weaklyConnectedStat_[int( v0 )] = uint8_t( std::min( count, 255 ) );
+            if ( calcFarSurfaceCached )
+                farSurfaceStat_[int( v0 )] = count >= 3 ? plane.getBestPlanef().distance( points[v0] ) : 0.f;
+            if ( calcBadNormalCached )
+            {
+                if ( normalSum.lengthSq() >= 0.09f )
+                {
+                    normalSum = normalSum.normalized();
+                    badNormalStat_[int( v0 )] = angle( normals[v0], normalSum );
+                }
+                else
+                    badNormalStat_[int( v0 )] = 180.f;
+            }
+        }, subProgress );
+
+        if ( !continued )
+            return unexpectedOperationCanceled();
+        subProgress = subprogress( progress, 0.4f, 1.f );
+    }
+
+    if ( calcSmallComponentsCached )
+    {
+        int counterProcessedVerts = 0;
+        const int lastPassVertsCount = int( lastPassVerts->count() );
+        const float counterMax = float( lastPassVertsCount );
+        const int counterDivider = std::max( lastPassVertsCount / 100, 1 );
+        for ( auto v0 : *lastPassVerts )
+        {
+            findPointsInBall( pointCloud.getAABBTree(), { pointCloud.points[v0], sqr( radius_ ) },
+                              [&] ( VertId v1, const Vector3f& )
+            {
+                if ( v0 < v1 && contains( validPoints_, v1 ) )
+                {
+                    unionFindStructure_.unite( v0, v1 );
+                }
+            } );
+            ++counterProcessedVerts;
+            if ( !reportProgress( subProgress, counterProcessedVerts / counterMax, counterProcessedVerts, counterDivider ) )
+                return unexpectedOperationCanceled();
+        }
+    }
+
+    return {};
+}
+
+void OutliersDetector::setParams( const OutlierParams& params )
+{
+    params_ = params;
+}
+
+Expected<VertBitSet> OutliersDetector::find( OutlierTypeMask mask, ProgressCallback progress /*= {}*/ )
+{
+    mask &= maskCached_;
+    if ( !bool( mask ) )
+        return unexpected( "Nothing to look for" );
+
+    int maxTaskCount = 0;
+    for ( int i = 0; i < 4; ++i )
+        if ( bool( OutlierTypeMask( 1 << i ) & mask ) )
+            ++maxTaskCount;
+
+    VertBitSet result;
+    if ( bool( mask & OutlierTypeMask::SmallComponents ) )
+    {
+        auto res = findSmallComponents( subprogress( progress, 0.f, 1.f / maxTaskCount ) );
+        if ( !res.has_value() )
+            return unexpected( res.error() );
+        result |= *res;
+    }
+    if ( bool( mask & OutlierTypeMask::WeaklyConnected ) )
+    {
+        auto res = findWeaklyConnected();
+        if ( !res.has_value() )
+            return unexpected( res.error() );
+        result |= *res;
+    }
+    if ( bool( mask & OutlierTypeMask::FarSurface ) )
+    {
+        auto res = findFarSurface();
+        if ( !res.has_value() )
+            return unexpected( res.error() );
+        result |= *res;
+    }
+    if ( bool( mask & OutlierTypeMask::AwayNormal ) )
+    {
+        auto res = findAwayNormal();
+        if ( !res.has_value() )
+            return unexpected( res.error() );
+        result |= *res;
+    }
+
+    return result;
+}
+
+Expected<VertBitSet> OutliersDetector::findSmallComponents( ProgressCallback progress /*= {}*/ )
+{
+    auto largeComponentsRes = PointCloudComponents::getLargeComponentsUnion( unionFindStructure_, validPoints_, params_.maxClusterSize + 1, progress );
+
+    if ( !largeComponentsRes.has_value() )
+        return unexpected( largeComponentsRes.error() );
+    return validPoints_ - *largeComponentsRes;
+}
+
+Expected<VertBitSet> OutliersDetector::findWeaklyConnected( ProgressCallback progress /*= {}*/ )
+{
+    VertBitSet result( validPoints_.find_last() + 1 );
+    const bool continued = BitSetParallelFor( validPoints_, [&] ( VertId v )
+    {
+        if ( weaklyConnectedStat_[int( v )] <= params_.maxNeighbors )
+            result.set( v );
+    }, progress );
+    if ( !continued )
+        return unexpectedOperationCanceled();
+    return result;
+}
+
+Expected<VertBitSet> OutliersDetector::findFarSurface( ProgressCallback progress /*= {}*/ )
+{
+    VertBitSet result( validPoints_.find_last() + 1 );
+    const float realMinHeight = params_.minHeight * radius_;
+    const bool continued = BitSetParallelFor( validPoints_, [&] ( VertId v )
+    {
+        if ( farSurfaceStat_[int( v )] > realMinHeight )
+            result.set( v );
+    }, progress );
+    if ( !continued )
+        return unexpectedOperationCanceled();
+    return result;
+}
+
+Expected<VertBitSet> OutliersDetector::findAwayNormal( ProgressCallback progress /*= {}*/ )
+{
+    VertBitSet result( validPoints_.find_last() + 1 );
+    const bool continued = BitSetParallelFor( validPoints_, [&] ( VertId v )
+    {
+        if ( badNormalStat_[int( v )] > params_.minAngle )
+            result.set( v );
+    }, progress );
+    if ( !continued )
+        return unexpectedOperationCanceled();
+    return result;
+}
+
+Expected<VertBitSet> findOutliers( const PointCloud& pc, const FindOutliersParams& params )
+{
+    OutliersDetector finder;
+    auto res = finder.prepare( pc, params.radius, params.mask, subprogress( params.progress, 0.f, 0.8f ) );
+    if ( !res.has_value() )
+        return unexpected( res.error() );
+    finder.setParams( params.finderParams );
+    return finder.find( params.mask, subprogress( params.progress, 0.8f, 1.f ) );
+}
+
+}

source/MRMesh/MROutlierPoints.h

@@ -0,0 +1,99 @@
+#pragma once
+#include "MRMeshFwd.h"
+#include "MRMesh/MRExpected.h"
+#include "MRUnionFind.h"
+#include "MRBitSet.h"
+#include "MRConstants.h"
+#include "MRFlagOperators.h"
+
+namespace MR
+{
+
+/// Parameters of various criteria for detecting outlier points
+struct OutlierParams
+{
+    /// Maximum points in the outlier component
+    int maxClusterSize = 20;
+    /// Maximum number of adjacent points for an outlier point
+    int maxNeighbors = 7;
+    /// Minimum distance (as proportion of search radius) to the approximate surface from outliers point
+    float minHeight = 0.3f;
+    /// Minimum angle of difference of normal at outlier points
+    /// @note available only if there are normals
+    float minAngle = PI_F / 3.f;
+};
+
+/// Types of outlier points
+enum class OutlierTypeMask
+{
+    SmallComponents = 1 << 0, ///< Small groups of points that are far from the rest
+    WeaklyConnected = 1 << 1, ///< Points that have too few neighbors within the radius
+    FarSurface = 1 << 2, ///< Points far from the surface approximating the nearest points
+    AwayNormal = 1 << 3, ///< Points whose normals differ from the average norm of their nearest neighbors
+    All = SmallComponents | WeaklyConnected | FarSurface | AwayNormal
+};
+MR_MAKE_FLAG_OPERATORS( OutlierTypeMask )
+
+/// A class for searching for outliers of points
+/// @detail The class caches the prepared search results, which allows to speed up the repeat search (while use same radius)
+class MRMESH_CLASS OutliersDetector
+{
+public:
+    OutliersDetector() = default;
+
+    /// Make a preliminary stage of outlier search. Caches the result
+    /// 
+    /// @param pc point cloud
+    /// @param radius radius of the search for neighboring points for analysis
+    /// @param mask mask of the types of outliers that are looking for
+    /// @param progress progress callback function
+    /// @return error text or nothing
+    MRMESH_API Expected<void> prepare( const PointCloud& pc, float radius, OutlierTypeMask mask, ProgressCallback progress = {} ); // calculate caches
+
+    /// Set search parameters
+    MRMESH_API void setParams( const OutlierParams& params );
+    /// Get search parameters
+    MRMESH_API const OutlierParams& getParams() const { return params_; }
+
+    /// Make an outlier search based on preliminary data
+    /// @param mask mask of the types of outliers you are looking for
+    MRMESH_API Expected<VertBitSet> find( OutlierTypeMask mask, ProgressCallback progress = {} ); // unite and calculate actual outliers
+
+    /// Get statistics on the number of neighbors for each point
+    MRMESH_API const std::vector<uint8_t>& getWeaklyConnectedStat() { return weaklyConnectedStat_; }
+
+private:
+    Expected<VertBitSet> findSmallComponents( ProgressCallback progress = {} );
+    Expected<VertBitSet> findWeaklyConnected( ProgressCallback progress = {} );
+    Expected<VertBitSet> findFarSurface( ProgressCallback progress = {} );
+    Expected<VertBitSet> findAwayNormal( ProgressCallback progress = {} );
+
+    float radius_ = 1.f;
+    OutlierParams params_;
+
+    // Cached data
+    UnionFind<VertId> unionFindStructure_; // SmallComponents
+    std::vector<uint8_t> weaklyConnectedStat_; // WeaklyConnected
+    std::vector<float> farSurfaceStat_; // FarSurface
+    std::vector<float> badNormalStat_; // AwayNormal
+
+    OutlierTypeMask maskCached_ = OutlierTypeMask( 0 ); // true means valid cache
+
+    VertBitSet validPoints_;
+};
+
+/// Outlier point search parameters
+struct FindOutliersParams
+{
+    OutlierParams finderParams; ///< Parameters of various criteria for detecting outlier points
+    float radius = 1.f; ///< Radius of the search for neighboring points for analysis
+
+    OutlierTypeMask mask = OutlierTypeMask::All; ///< Mask of the types of outliers that are looking for
+
+    ProgressCallback progress = {}; ///< Progress callback
+};
+
+/// Finding outlier points
+MRMESH_API Expected<VertBitSet> findOutliers( const PointCloud& pc, const FindOutliersParams& params );
+
+}