diff --git a/source/MRMesh/MRMesh.vcxproj b/source/MRMesh/MRMesh.vcxproj index 5486fac46b0c..1e9e03f9180c 100644 --- a/source/MRMesh/MRMesh.vcxproj +++ b/source/MRMesh/MRMesh.vcxproj @@ -144,6 +144,7 @@ + @@ -401,6 +402,7 @@ + diff --git a/source/MRMesh/MRMesh.vcxproj.filters b/source/MRMesh/MRMesh.vcxproj.filters index 83aa8d649e75..fce99b06dc49 100644 --- a/source/MRMesh/MRMesh.vcxproj.filters +++ b/source/MRMesh/MRMesh.vcxproj.filters @@ -1272,6 +1272,9 @@ Source Files\AABBTree + + Source Files\PointCloud + @@ -2093,6 +2096,9 @@ Source Files\AABBTree + + Source Files\PointCloud + diff --git a/source/MRMesh/MROutlierPoints.cpp b/source/MRMesh/MROutlierPoints.cpp new file mode 100644 index 000000000000..a2a6f28732df --- /dev/null +++ b/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 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( numVerts ); + if ( calcWeaklyConnectedCached ) + weaklyConnectedStat_ = std::vector( numVerts ); + if ( calcFarSurfaceCached ) + farSurfaceStat_ = std::vector( numVerts ); + if ( calcBadNormalCached ) + badNormalStat_ = std::vector( 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 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 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 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 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 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 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 ) ); +} + +} diff --git a/source/MRMesh/MROutlierPoints.h b/source/MRMesh/MROutlierPoints.h new file mode 100644 index 000000000000..046e9024d2d4 --- /dev/null +++ b/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 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 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& getWeaklyConnectedStat() { return weaklyConnectedStat_; } + +private: + Expected findSmallComponents( ProgressCallback progress = {} ); + Expected findWeaklyConnected( ProgressCallback progress = {} ); + Expected findFarSurface( ProgressCallback progress = {} ); + Expected findAwayNormal( ProgressCallback progress = {} ); + + float radius_ = 1.f; + OutlierParams params_; + + // Cached data + UnionFind unionFindStructure_; // SmallComponents + std::vector weaklyConnectedStat_; // WeaklyConnected + std::vector farSurfaceStat_; // FarSurface + std::vector 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 findOutliers( const PointCloud& pc, const FindOutliersParams& params ); + +}