Creating grid meshes based on arbitrary shape input

[mariuszhermansdorfer]

Hi @Grantim and @Fedr,

I believe to have seen this before but can't find the corresponding function now. How can I create a gridded mesh based on arbitrary shape input, a mask, and a desired spacing value?

For a rectangle it would look like this:

Input shapes would be closed curves or open meshes (both convex and concave). The resulting grid would be trimmed to the outline, such that no vertex lays outside of it:

The grid would be projected from the World.XY plane, even if the input shape is non-planar. To be clear, the resulting vertices would follow the input shape, kind of like this: Input mesh:

Resulting gridded mesh:

Top view:

I'd appreciate it if you could point me in the right direction!

[Grantim]

Hello!

I believe Distance map functionality is the thing you are looking for https://github.com/MeshInspector/MeshLib/blob/master/source/MRMesh/MRDistanceMap.h

also there is function to create mesh without computing whole distance map https://github.com/MeshInspector/MeshLib/blob/9d16a0e18076e843c6a59e5135705181b9ea7ed7/source/MRMesh/MRRegularGridMesh.h#L20-L24

here validator is function to say if vertex of grid should be present in mesh or not, and positioner is function that returns world position of the vertex

[mariuszhermansdorfer]

Thanks @Grantim!

Could you please show an example on how to call the makeRegularGridMesh function on an existing mesh or a closed curve?

[Grantim]

Sure: for meshes you can use distance map conversion: https://github.com/MeshInspector/MeshLib/blob/9d16a0e18076e843c6a59e5135705181b9ea7ed7/source/MRMesh/MRDistanceMap.h#L170-L173

and then back to mesh: https://github.com/MeshInspector/MeshLib/blob/9d16a0e18076e843c6a59e5135705181b9ea7ed7/source/MRMesh/MRDistanceMap.h#L285-L286

MeshToDistanceMapParams and DistanceMapToWorld structure specifies properties of conversion

for 2d curves there is similar function (each distance map value is distance to the curve in plane): https://github.com/MeshInspector/MeshLib/blob/9d16a0e18076e843c6a59e5135705181b9ea7ed7/source/MRMesh/MRDistanceMap.h#L217-L222

and then back to mesh (same as in first example): https://github.com/MeshInspector/MeshLib/blob/9d16a0e18076e843c6a59e5135705181b9ea7ed7/source/MRMesh/MRDistanceMap.h#L285-L286

---

makeRegularGridMesh - this function can be used to create mesh without creating DistanceMap:

    makeRegularGridMesh( resolutionX, resolutionY, [&]( size_t x, size_t y )
    {
        return !isBlueVertex(x,y); // blue vertices are invalid and should not be present in result mesh
    },
                                            [&]( size_t x, size_t y )
    {
        return originPoint + Vector3f( x * pixelSizeX, y * pixelSizeY, calcVertexHeight( x, y ) );
    } );

please note that it is pseudocode and it may rise some warnings or errors, its purpose to demonstrate the idea of makeRegularGridMesh function

---

Also, we had some discussion about distance maps here https://github.com/MeshInspector/MeshLib/issues/1382

[mariuszhermansdorfer]

Super! It was #1382 that I was looking for and DistanceMap is the way to go about it.

One more question, is there a way to limit the output of a DistanceMap to a mask?

My input mesh is a circle, the mask is a closed curve in an arbitrary plane. I'd like the resulting mesh to be only limited to the vertices whose projection to the World.XY plane falls within the outline of the mask.

My current approach would be to start with the mask, create a DistanceMap out of it, and raycast from its vertices to the input mesh. But maybe there is already a better approach within MeshLib?

[Grantim]

I think algorithm here should look like this:

1. Project input curve to 2d space:

Polyline2 polyline( inputCurve );

https://github.com/MeshInspector/MeshLib/blob/0f80451a175b17ba459ceab8846e2395df6ecb3d/source/MRMesh/MRPolyline.h#L28-L29

2. Then build MeshToDistanceMapParams based on Curve like this:

MeshToDistanceMapParams params;
params.xRange = Vector3f::plusX() * pixelSize;
params.yRange = Vector3f::plusY() * pixelSize;
auto box = polyline.getBoundingBox();
const int offsetPixels = 2;
params.orgPoint = to3dim( box.min ) - offsetPixels * ( params.xRange + params.yRange ); // 2 pixels offset of minimum point
params.resolution = Vector3i( box.size() / pixelSize ) + ( 2 * offsetPixels + 1 ) * Vector3i::diagonal( 1 ); // + 1 to `ceil` division result
params.direction = Vector3f::plusZ(); // may be Vector3f::minusZ(); as well

note that orgPoint.z is responsible for raycast origin Z coordinate

3. Then just computeDistanceMap

auto dm = computeDistanceMap( inputMesh, params );

At this point distance map has all values around input curve

4. Invalidate distance map values outside of curve:

DistanceMapToWorld dmToWorld(params); // params to backproject pixels to world
ParallelFor( 0, int( dm.size() ), [&]( int i ) 
{
    auto intPos = dm.toPos( i );
    auto coord = dm.unproject( intPos.x, intPos.y, dmToWorld );
    if ( !coord )
        return; // pixel value is already invalid
    if ( !isPointInsidePolyline( polyline, to2dim( *coord ) )
        dm.unset( i );
} );

5. Now distance map only has valid values inside curve so simple distanceMapToMesh will work:

auto mesh = distanceMapToMesh( dm, dmToWorld );

[mariuszhermansdorfer]

Hey @Grantim,

I'm testing the above code and have a few questions still. Could you please help me resolve these?

1. DistanceMap with input mesh works super good!
2. DistanceMap with input curve is not really what I need. If I understand it correctly, it creates a grid with values measuring distance to the curve itself:

Top view:

Perspective view:

I'd like to have a gridded mesh conforming to the curve. Similar to how this triangle mesh works, but with a regular grid:

I could of course, first create the triangle mesh, and then grid it using the distance map as above but maybe there is a better way. My current code:

PINVOKE Mesh* CurveToGridMesh( const float coordinates[], const int coordinatesLength, const float resolution )
{
	float res = resolution < 0.01 ? 0.01f : resolution;

	std::vector<Vector2f> coordinateVector;
	coordinateVector.reserve( coordinatesLength / 3 );

	for ( size_t i = 0; i < coordinatesLength; i += 3 )
	{
		coordinateVector.emplace_back( coordinates[i], coordinates[i + 1] );
	}
	DecimatePolylineSettings2 dps{ .maxError = MAXERROR };
	decimateContour( coordinateVector, dps );

	coordinateVector.emplace_back( coordinateVector[0] ); // repeat first vertex since addFromPoints removes it
	Polyline2 polyline;
	polyline.addFromPoints( coordinateVector.data(), coordinateVector.size() );

	auto box = polyline.getBoundingBox();

	ContourToDistanceMapParams params;
	params.orgPoint = box.min;
	params.resolution = Vector2i( static_cast< int >( std::ceil( box.size().x / res) ), static_cast< int >( std::ceil( box.size().y / res) ));


	DistanceMap proposedMap = distanceMapFromContours( polyline, params );

	Mesh* gridMesh = new Mesh( distanceMapToMesh( proposedMap, DistanceMapToWorld( params ) ) );

	return gridMesh;
}

3. DistanceMap with input mesh and curve, where the resulting grid follows the mesh, but is masked to the curve's projection onto the mesh The code you've provided unfortunately results in invalid meshes:

PINVOKE Mesh* ProjectGridToMesh( const Mesh* mesh, const float resolution, const float coordinates[], const int coordinatesLength )
{
	float res = resolution < 0.01 ? 0.01f : resolution;

	std::vector<Vector2f> coordinateVector;
	coordinateVector.reserve( coordinatesLength / 3 );

	for ( size_t i = 0; i < coordinatesLength; i += 3 )
	{
		coordinateVector.emplace_back( coordinates[i], coordinates[i + 1] );
	}
	DecimatePolylineSettings2 dps{ .maxError = MAXERROR };
	decimateContour( coordinateVector, dps );

	coordinateVector.emplace_back( coordinateVector[0] ); // repeat first vertex since addFromPoints removes it
	Polyline2 polyline;
	polyline.addFromPoints( coordinateVector.data(), coordinateVector.size() );

	MeshToDistanceMapParams params;
	params.xRange = Vector3f::plusX() * res;
	params.yRange = Vector3f::plusY() * res;
	auto box = polyline.getBoundingBox();
	const float offsetPixels = 2;
	
	params.orgPoint = to3dim( box.min ) - offsetPixels * ( params.xRange + params.yRange );
	Vector2f r = box.size() / res  + ( 2.0f * offsetPixels + 1.0f ) * Vector2f::diagonal( 1 );
	params.resolution = Vector2i( static_cast< int > (r.x), static_cast< int > (r.y));
	params.direction = Vector3f::plusZ(); 

	auto dm = computeDistanceMap( *mesh, params );

	DistanceMapToWorld dmToWorld( params ); // params to backproject pixels to world
	ParallelFor( 0, int( dm.size() ), [&] ( int i )
	{
		auto intPos = dm.toPos( i );
		auto coord = dm.unproject( intPos.x, intPos.y, dmToWorld );
		if ( !coord )
			return; // pixel value is already invalid
		if ( !isPointInsidePolyline( polyline, to2dim( *coord ) ))
			dm.unset( i );
	} );

	Mesh* result = new Mesh(distanceMapToMesh( dm, dmToWorld ));

	return result;
}

[Grantim]

Hello!

1. Good to hear!
2. Distance map built from polyline really have distances to polyline measured in xy plane. To get mesh that you need you need to manually call makeRegularGridMesh, or triangulate contour and create Distance map out of this mesh. (If contour is always planar there is an easy way to call makeRegularGridMesh otherwise you need to specify surface function for the contour)
3. Code looks fine, can you please attach data to reproduce issue. Only small update for the code

PINVOKE Mesh* ProjectGridToMesh( const Mesh* mesh, const float resolution, const float coordinates[], const int coordinatesLength )
{
	float res = resolution < 0.01 ? 0.01f : resolution;

	std::vector<Vector2f> coordinateVector;
	coordinateVector.reserve( coordinatesLength / 3 );

	for ( size_t i = 0; i < coordinatesLength; i += 3 )
	{
		coordinateVector.emplace_back( coordinates[i], coordinates[i + 1] );
	}
	coordinateVector.emplace_back( coordinateVector[0] ); // this is needed to understand that contour is closed
        // better create polyline before decimating, because decimateContour will create it too if you pass simple contours
        Polyline2 polyline( coordinateVector );
	DecimatePolylineSettings2 dps{ .maxError = MAXERROR };
	decimateContour( polyline, dps );

	MeshToDistanceMapParams params;
	params.xRange = Vector3f::plusX() * res;
	params.yRange = Vector3f::plusY() * res;
	auto box = polyline.getBoundingBox();
	const float offsetPixels = 2;
	
	params.orgPoint = to3dim( box.min ) - offsetPixels * ( params.xRange + params.yRange );
	Vector2f r = box.size() / res  + ( 2.0f * offsetPixels + 1.0f ) * Vector2f::diagonal( 1 );
	params.resolution = Vector2i( static_cast< int > (r.x), static_cast< int > (r.y));
	params.direction = Vector3f::minusZ(); // mb try minusZ
        params.orgPoint.z = mesh.getBoundingBox().max.z + offsetPixels; // make sure that origin point is higher than mesh

	auto dm = computeDistanceMap( *mesh, params );

	DistanceMapToWorld dmToWorld( params ); // params to backproject pixels to world
	ParallelFor( 0, int( dm.size() ), [&] ( int i )
	{
		auto intPos = dm.toPos( i );
		auto coord = dm.unproject( intPos.x, intPos.y, dmToWorld );
		if ( !coord )
			return; // pixel value is already invalid
		if ( !isPointInsidePolyline( polyline, to2dim( *coord ) ))
			dm.unset( i );
	} );

	Mesh* result = new Mesh(distanceMapToMesh( dm, dmToWorld ));

	return result;
}

[mariuszhermansdorfer]

Thanks @Grantim!

For 2. I'll triangulate the input contour first and run a distance map on this mesh. A few comments to 3:

Polyline2 doesn't have a constructor which would take a std::vector<Vector2f>:

Polyline2 polyline( coordinateVector );

So I replaced it with:

	Polyline2 polyline;
	polyline.addFromPoints( coordinateVector.data(), coordinateVector.size() );
	DecimatePolylineSettings2 dps{ .maxError = MAXERROR };
	decimatePolyline( polyline, dps );

But it still results in an invalid mesh. Attached the sample data I'm working with Projection.zip

[Grantim]

Oh really, it should look this

Polyline2 polyline( { coordinateVector } );

---

I'll take a look at the data and come back to you

[Grantim]

There is only mesh in the file attached, could you please also provide curve and parameters

---

EDIT: I found bspline object in obj file, unfortunately we now does not have import functionality for it

[mariuszhermansdorfer]

OK, this is a super simple polyline with the following vertices:

v -18.360707888498766 -0.32774690489587788 14.422012291528848
v -7.6229552286595457 -1.083248162243623 14.422012291528848
v -4.9839390508489556 7.5979974866691968 14.422012291528848
v -18.873149284948191 7.8239787801343388 14.422012291528848

But you can draw any closed polyline and a resolution of say 1.2 meters. It always results in an invalid mesh, independently of the inputs.

[Grantim]

Thanks @mariuszhermansdorfer ! There was an error in setting up parameters for distance map, it should be like this:

    MeshToDistanceMapParams params;
    auto box = polyline.getBoundingBox();
    const float offsetPixels = 2;
    box.include( box.min - res * offsetPixels * Vector2f::diagonal( 1 ) );
    box.include( box.max + res * offsetPixels * Vector2f::diagonal( 1 ) );

    params.orgPoint = to3dim( box.min );
    Vector2f r = box.size() / res + Vector2f::diagonal( 1 );
    params.xRange = Vector3f::plusX() * box.size().x;
    params.yRange = Vector3f::plusY() * box.size().y;
    params.resolution = Vector2i( static_cast< int > ( r.x ), static_cast< int > ( r.y ) );
    params.direction = Vector3f::minusZ(); // mb try minusZ
    params.orgPoint.z = mesh.getBoundingBox().max.z + offsetPixels; // make sure that origin point is higher than mesh

xRange - is full size of distance map X side (previously I said that it was pixel size, it it was mistake) yRange - is full size of distance map Y side (previously I said that it was pixel size, it it was mistake)

So full code should look like this

PINVOKE Mesh* ProjectGridToMesh( const Mesh* mesh, const float resolution, const float coordinates[], const int coordinatesLength )
{
	float res = resolution < 0.01 ? 0.01f : resolution;

	std::vector<Vector2f> coordinateVector;
	coordinateVector.reserve( coordinatesLength / 3 );

	for ( size_t i = 0; i < coordinatesLength; i += 3 )
	{
		coordinateVector.emplace_back( coordinates[i], coordinates[i + 1] );
	}
	coordinateVector.emplace_back( coordinateVector[0] ); // this is needed to understand that contour is closed
        // better create polyline before decimating, because decimateContour will create it too if you pass simple contours
        Polyline2 polyline( { coordinateVector } ); // added {} here
	DecimatePolylineSettings2 dps{ .maxError = MAXERROR };
	decimateContour( polyline, dps );

	MeshToDistanceMapParams params;
	auto box = polyline.getBoundingBox();
	const float offsetPixels = 2;
	box.include( box.min - res * offsetPixels * Vector2f::diagonal( 1 ) );
	box.include( box.max + res * offsetPixels * Vector2f::diagonal( 1 ) );

	params.orgPoint = to3dim( box.min );
	Vector2f r = box.size() / res + Vector2f::diagonal( 1 );
	params.xRange = Vector3f::plusX() * box.size().x;
	params.yRange = Vector3f::plusY() * box.size().y;
	params.resolution = Vector2i( static_cast< int > ( r.x ), static_cast< int > ( r.y ) );
	params.direction = Vector3f::minusZ(); // mb try minusZ
	params.orgPoint.z = mesh.getBoundingBox().max.z + offsetPixels; // make sure that origin point is higher than mesh

	auto dm = computeDistanceMap( *mesh, params );

	DistanceMapToWorld dmToWorld( params ); // params to backproject pixels to world
	ParallelFor( 0, int( dm.size() ), [&] ( int i )
	{
		auto intPos = dm.toPos( i );
		auto coord = dm.unproject( intPos.x, intPos.y, dmToWorld );
		if ( !coord )
			return; // pixel value is already invalid
		if ( !isPointInsidePolyline( polyline, to2dim( *coord ) ))
			dm.unset( i );
	} );

	Mesh* result = new Mesh(distanceMapToMesh( dm, dmToWorld ));

	return result;
}

[mariuszhermansdorfer]

Thanks @Grantim! Confirmed :)

[mariuszhermansdorfer]

Hey @Grantim, I'm revisiting this topic with two, hopefully minor, feature requests:

1. As far as I can understand DistanceMap internally performs RayMeshIntersections: https://github.com/MeshInspector/MeshLib/blob/9fe3308bc26b71ea3c886a6ae5013bfa7e6c9062/source/MRMesh/MRDistanceMap.cpp#L318-L339

Could there be an option to use MRCuda-based raycasting? My use cases quite often require 4M+ grid points to be created and DistanceMap creation slows down noticeably.

2. DistanceMap needs to have a normal direction to perform the raycasting. How would you convert the vertical parts of a cylinder to a gridded mesh? In this case, each face has its own normal.

[EDIT] Would it make sense to voxelize this kind of geometry maybe? Would this be faster or slower than raycasting?

[Grantim]

Hello @mariuszhermansdorfer !

1. Could there be an option to use MRCuda-based raycasting?

Sure, it can be done, but now we don't have CUDA implementation, we will add it.

2. Distance map does not really require face normals, it just raycasts in given direction (all rays are parallel) to first intersection, so in this case result will be simple circle (if raycast direction is parallel to cylinder walls)

Would it make sense to voxelize this kind of geometry maybe? Would this be faster or slower than raycasting?

It depends on what result you want to achieve, distance map is 2d grid, while voxelization is 3d grid, so in common it is N times slower

[mariuszhermansdorfer]

Thanks @Grantim!

What I'm trying to achieve is that the points would be regularly distributed on the vertical surface of the cylinder (or any arbitrarily shaped mesh). Something like this but in a regular grid:

[Grantim]

Thanks for explanation!

Regular DistanceMap approach will not work here, but you still can create custom DistanceMap, for example in radial coordinates:

• lets say that X axis of the distance map is central axis of cylinder
• Y axis of distance map is angle of the ray in plane parallel to cylinder basis
• Raycast each ray to the cylinder and fill Distance map with results
• call makeRegularGridMesh function with custom callbacks to follow radial coordinate system

[mariuszhermansdorfer]

Thanks, but the cylinder is just a very specific case. A more general case could be this one:

Or this one:

Essentially, there is no obvious way to determine the direction of raycasting.

[Grantim]

I see, so in this case I think there is no good way to achieve regular mesh, (because in common case there could be undercuts or some topology issues). Maybe there is another way to solve your problem without creating regular grid mesh?

You can try voxelization but I it will output different mesh than regular grid

auto volume = meshToDistanceVolume(mesh).volume(); // need to setup parameters here especially `maxDistSq = sqr(2*voxelSize)`
auto mesh = marchingCubes(volume ); // also need to setup parameters here

[mariuszhermansdorfer]

The intended use case is for sunlight analysis. The surfaces above represent walls of buildings. For each point on the surface, I would perform raycasting with a series of vectors representing sun directions throughout a given day. For each point, I can count how many hours it received direct sunlight and visualize this as a mesh.

Creating a regular grid assures a uniform distribution of these sample points. But maybe we could approach this in a different way? Maybe poisson disk sampling, or any other uniform sampling method could work? I tried remeshing with a set edge length, but it was approx. 10x slower than DistanceMap.

[Edit] Here is an example using only the ground plane for analysis, I'd like to extend this to building walls as well:

https://github.com/MeshInspector/MeshLib/assets/49192999/13fa81ee-ff46-4db6-85e8-10150e82f9e3

[Grantim]

As I understand You raycast from the plane to the sun, maybe it can be done vice versa? Raycast from the sun to the ground, so you will both get information about sunlight on plane and about sunlight on building, with no need in any regular mesh, will it work?

[mariuszhermansdorfer]

Sure, but I still need the targets on the building and then a way of color coding the results (vertex colors of a mesh).

[Grantim]

In this way re-mesh looks like best option here. I think can call it in parallel for all buildings, so it may improve performance

[mariuszhermansdorfer]

Super, I'll resort to re-meshing for nonplanar surfaces. Would be cool if you could add the CUDA-based DistanceMap version for maximum performance.

[Grantim]

Hello again! In new version of meshlib we've added cuda version of computeDistanceMap function, please give it a try

[mariuszhermansdorfer]

Thanks @Grantim! It works as expected, but unfortunately it seems that the bottleneck is elsewhere and it's caused by the distanceMapToMesh() call. Switching between the CPU and GPU distanceMaps results only in a minor speedup.

DistanceMap computeDistanceMap(const Mesh& mesh, const MeshToDistanceMapParams& params, bool useGPU) {
	if (MR::Cuda::isCudaAvailable() && useGPU)
		return MR::Cuda::computeDistanceMap(mesh, params);
	else
		return computeDistanceMap(MeshPart(mesh), params);
}

PINVOKE Mesh* GridRemesh( const Mesh* mesh, const float resolution, bool useGPU)
{
	float res = std::max(resolution, MIN_RESOLUTION);

	const auto& bbox = mesh->getBoundingBox();
	auto boxSize = bbox.size();

	MeshToDistanceMapParams params;
	params.direction = Vector3f( 0, 0, 1 );
	params.allowNegativeValues = false;
	params.orgPoint = bbox.min;
	params.xRange = Vector3f( boxSize.x, 0.f, 0.f );
	params.yRange = Vector3f( 0.f, boxSize.y, 0.f );
	params.resolution = Vector2i( static_cast< int >( std::ceil( boxSize.x / res ) ), static_cast< int >( std::ceil( boxSize.y / res ) ) );
	
	auto proposedMap = computeDistanceMap(*mesh, params, useGPU);

	Mesh* gridMesh = new Mesh( distanceMapToMesh( proposedMap, params ).value() );

	return gridMesh;
}

[Grantim]

Oh, I see, in fact we have already optimized this algorithm and I don't think that it can be speed up much more, but we'll have a look at this

[mariuszhermansdorfer]

Having looked at the code here: https://github.com/MeshInspector/MeshLib/blob/a0fa86fad845c58c07930ab441b23ee3dfa01d1d/source/MRMesh/MRRegularGridMesh.cpp#L15

I have a suggestion which could work in this specific case. Since we know all dimensions of the grid, we can safely assume a regular triangulation pattern (n, n + 1, n * stride). Couldn’t we leverage this knowledge to use this algorithm:

1. Add all grid points as vertices to the mesh
2. Create all faces based on the grid pattern
3. Loop over all grid points, if invalid, remove all corresponding faces
4. Pack mesh

I’m not sure about its speed, but it seems that it could avoid a lot of checks going on in the current logic.

[Grantim]

Old implementation of the function looked like this: https://github.com/MeshInspector/MeshLib/blob/50a76c45f9b71b7803cf94a1f5f10ab95f0fbd63/source/MRMesh/MRRegularGridMesh.cpp New version is faster and have better packing (optimal for parallel decimation).