PointCloud sampling

[mariuszhermansdorfer]

Hello!

I'm dealing with relatively large point-clouds (5-10M points). For my workflow, I need to sample these with a regular grid of approx. 10 000 sample points and find the id of the closest points from the point cloud to each sample point.

In Rhino 5M points and 10 000 samples takes approx. 20s to compute. Do you think MeshLib could do this faster?

I need this for a wind simulation study. Here is an example for a much smaller point cloud which runs in real-time:

https://www.youtube.com/watch?v=l33cxovz_3I

[Grantim]

Hello!
Now we have two sampling algorithms in MeshLib

MeshLib/source/MRMesh/MRGridSampling.h

Lines 15 to 18 in 69f5c9f

	/// performs sampling of cloud points;
	/// subdivides point cloud bounding box on voxels of approximately given size and returns at most one point per voxel;
	/// returns std::nullopt if it was terminated by the callback
	MRMESH_API std::optional<VertBitSet> pointGridSampling( const PointCloud& cloud, float voxelSize, const ProgressCallback & cb = {} );

and

MeshLib/source/MRMesh/MRUniformSampling.h

Lines 11 to 14 in 69f5c9f

	/// Sample vertices, removing ones that are too close;
	/// returns std::nullopt if it was terminated by the callback
	/// \ingroup PointCloudGroup
	MRMESH_API std::optional<VertBitSet> pointUniformSampling( const PointCloud& pointCloud, float distance, const ProgressCallback & cb = {} );

it will sampled VertBitSet that can be used to set valid points in PointCloud

[mariuszhermansdorfer]

I'd like to specify the exact location of my sampling points. The above methods don't seem to allow this flexibility, do they?

[Grantim]

Oh, thanks, I think I got it now. Is it correct: you have 10000 sample points and for each of them you need to find closest point from 5M point cloud?

[mariuszhermansdorfer]

Yeah, that's correct. I want to be able to define these sample points wherever I want and for each of them find the closest point from the big point cloud.

[Grantim]

So I can suggest code like this:

std::vector<MR::Vector3f> samplePoints(10000);
std::vector<MR::Vector3f> closestResult(10000);
// fill up with points you need
fillSamplePoints( samplePoints );

// project sample points on bigPointCloud  in parallel
MR::ParallelFor( samplePoints, [&]( size_t i ) 
{
    closestResult[i] = bigPointCloud.points[findProjectionOnPoints( samplePoints[i], bigPointCloud ).vId];
} );

[mariuszhermansdorfer]

Thanks @Grantim, I'll give it a shot and will report back.

[mariuszhermansdorfer]

Hmm, somehow I managed to create a memory leak somewhere and the code breaks without any warnings:

This is how I create the point cloud:

PINVOKE PointCloud* CreatePointCloud( const float coordinates[], const int coordinatesLength )
{
	std::vector<MR::Vector3f> coords( coordinatesLength / 3 );
	for ( int i = 0; i < coordinatesLength; i += 3 )
		coords.emplace_back( Vector3f( coordinates[i], coordinates[i + 1], coordinates[i + 2] ) );

	PointCloud* pc = new PointCloud();
	pc->points.vec_ = coords;

	return pc;
}

And since I am only interested in IDs of points in the point cloud, this is how I modified your suggestion:

PINVOKE int* SamplePointCloud( PointCloud* pc, const float coordinates[], const int coordinatesLength )
{
	std::vector<Vector3f> samplePoints( coordinatesLength / 3 );
	for ( int i = 0; i < coordinatesLength; i += 3 )
		samplePoints.push_back( Vector3f( coordinates[i], coordinates[i + 1], coordinates[i + 2] ) );

	int * pointIds = new int[coordinatesLength / 3];
	MR::ParallelFor( samplePoints, [&] ( size_t i )
	{
		pointIds[i] = ( uint32_t )findProjectionOnPoints( samplePoints[i], *pc ).vId;
	} );

	return pointIds;
}

Any ideas what could be wrong with this code?

[mariuszhermansdorfer]

Ah, sorry rookie mistake with how I created the vectors... This code doesn't crash but only returns -1 as pointIds.

PINVOKE PointCloud* CreatePointCloud( const float coordinates[], const int coordinatesLength )
{
	std::vector<MR::Vector3f> coords( coordinatesLength / 3 );
	for ( size_t i = 0, j = 0; i < coordinatesLength; i += 3, j++ )
		coords[j] = ( Vector3f( coordinates[i], coordinates[i + 1], coordinates[i + 2] ) );

	PointCloud* pc = new PointCloud();
	pc->points.vec_ = coords;

	return pc;
}

PINVOKE int* SamplePointCloud( PointCloud* pc, const float coordinates[], const int coordinatesLength )
{
	std::vector<Vector3f> samplePoints( coordinatesLength / 3 );
	for ( size_t i = 0, j = 0; i < coordinatesLength; i += 3, j++ )
		samplePoints[j] = Vector3f( coordinates[i], coordinates[i + 1], coordinates[i + 2] );

	int* pointIds = new int[coordinatesLength / 3];
	MR::ParallelFor( samplePoints, [&] ( size_t i )
	{
		pointIds[i] = ( uint32_t )findProjectionOnPoints( samplePoints[i], *pc ).vId.get();
	} );

	return pointIds;
}

[Fedr]

After these lines

	PointCloud* pc = new PointCloud();
	pc->points.vec_ = coords;

please add

	pc->validPoints.resize( pc->points.size(), true );

(otherwise it is considered that all points in the vector are invalid.)

And as to memory leaks, they are because you use raw owning pointers and probably forget to deallocate objects manually. Please replace PointCloud* with std::unique_ptr<PointCloud> or std::shared_ptr<PointCloud> or just PointCloud, e.g.:

PointCloud CreatePointCloud( const float coordinates[], const int coordinatesLength )
{
	PointCloud pc;
	pc.points.reserve( coordinatesLength / 3 );
	for ( size_t i = 0, j = 0; i < coordinatesLength; i += 3, j++ )
		pc.points.push_back( Vector3f( coordinates[i], coordinates[i + 1], coordinates[i + 2] ) );

	pc.validPoints.resize( pc.points.size(), true );
	return pc;
}

[mariuszhermansdorfer]

Thanks @Grantim,

this was the missing piece:

pc->validPoints.resize( pc->points.size(), true );

Now, sampling works but it returns weird results. To make sure we are on the same page, my point cloud is volumetric, this is its side view:

The name findProjectionOnPoints implies that it gets flattened somehow. Is that correct? I'm trying to find the closest point in the 3d space, not projected on a plane.

[mariuszhermansdorfer]

Sorry, my bad. The mistake was elsewhere. Now, everything works very well!

Another question, can I trim the original point cloud with a box while preserving the original point indexing?

[Fedr]

You can set the bits in pc->validPoints to false for the points outside the box, and then call pc->invalidateCaches().

[mariuszhermansdorfer]

Thanks @Fedr, and how do I identify the inside/outside bits given the point cloud and the box?

[Fedr]

In Box3f with have a method bool contains( const Point3f & pt ) const. You can use it to test whether the point is within the box. Then call pc->validPoints.reset( i ) for such points.

[mariuszhermansdorfer]

@Grantim, @Fedr,

I don't know how you guys do it, but you keep delivering amazing solutions!

Remember when I mentioned that Rhino takes 30s to sample the point cloud? In MeshLib this is down to 10ms for 4x the amount of samples!
I'm speechless.

Again, huge thanks for your help! Please close this issue as resolved.

[Grantim]

Thank you!

[mariuszhermansdorfer]

Apologies for reopening. I remember there was a function to create a regular grid of points based on a rectangle and desired resolution. Could you please help me find it?

Inputs: MinCoordinate, MaxCoordinate, Resolution
Output: Vector of points

[Grantim]

I don't think we had anything like this, maybe you mean these functions?

MeshLib/source/MRMesh/MRRegularGridMesh.h

Lines 20 to 27 in 70732d4

	// Creates regular mesh with points in valid grid lattice
	MRMESH_API Mesh makeRegularGridMesh( size_t width, size_t height,
	const RegularGridLatticeValidator& validator,
	const RegularGridLatticePositioner& positioner,
	const RegularGridMeshFaceValidator& faceValidator = {} );
	// Creates regular mesh from monotone (connects point with closed x, y neighbors) points
	MRMESH_API Expected<Mesh, std::string> makeRegularGridMesh( VertCoords pc, ProgressCallback cb = {} );

[mariuszhermansdorfer]

Yes, this one!

How can I use this to make a regular grid with a user-defined resolution based on a rectangle oriented in an arbitrary plane?

[Grantim]

You can call it like this:

// user defined width and height
// RegularGridLatticeValidator - always return true so all points will be present
// RegularGridLatticePositioner - returns position of point with given x and y
auto resMesh = makeRegularGridMesh( width, height, [&]( size_t, size_t ) { return true; },  
[&]( size_t x, size_t y) { return Vector3f(/* defined position of this point based on rectangle orientation*/); } );

this will result as grid mesh with points in positions from RegularGridMeshFaceValidator callback

[mariuszhermansdorfer]

Sorry @Grantim, but I think I need more help here.

Given the following inputs:
rectangle width = 100;
rectangle height = 50;
rectangle minPoint = {0, 0, 0)}:
rectangle normalVector = {0.826724, -0.0802, 0.556862}
grid resolution = 1;

How can I get the desired grid of 100 x 50 points aligned with the input rectangle?

I only understand the width/height part but am not sure how to set the rest.

auto resMesh = makeRegularGridMesh( 100, 50, [&]( size_t, size_t ) { return true; },  
[&]( size_t x, size_t y) { return Vector3f(/* defined position of this point based on rectangle orientation*/); } );

[Grantim]

Hope this code will help

Vector3f minPoint = Vector3f{0, 0, 0)};
Vector3f normalVector = Vector3f{0.826724f, -0.0802f, 0.556862f};
float resolution = 1.0f; // this is step size

auto basis = normalVector.perpendicular(); // find basis vectors of the plane
 
auto resMesh = makeRegularGridMesh( 100, 50, [&]( size_t, size_t ) { return true; },  
[&]( size_t x, size_t y) 
{
    // locate each point with x,y steps along basis 
    return 
        minPoint + // start position
        resolution * float( x ) * basis.first + // do `x` steps with size=`resolution` along first basis vector
        resolution* float( y ) * basis.second; // do `y` steps with size=`resolution` along second basis vector
} );