Interactive Contour Labels

[mariuszhermansdorfer]

This one is a long shot and not directly related to MeshLib so please bear with me while I explain. In my application, I rely on a GLSL shader to draw contour lines interactively. The majority of the logic takes place in the Geometry Shader, which essentially implements this algorithm.

It works very fast even for complex meshes and I'm happily sharing the code below if you wanted to incorporate it into the MeshInspector.

2023-07-18.16-17-50.mp4

Things slow down, however, when I start displaying labels. The GLSL shader, is only aware of one primitive (triangle) at a time and generates unconnected line segments. I couldn't figure out a way of joining these segments on the GPU, hence I copy them over to the CPU and sort them there. It gets me the results I want, but leads to noticeable slow downs with more complex meshes.

Do you have any suggestions on how to approach it in a more efficient way?

The GLSL shader code:

//vertexShader 
#version 330
layout(location = 0) in vec3 _meshVertex;
layout(location = 1) in vec3 _meshNormal;

out vec3 _normal;

void main() {
_normal = vec3(_meshNormal * 0.01); // offset contours by a small amount to prevent Z-fighting
gl_Position = vec4(_meshVertex + _normal, 1.0);
}

// geometryShader
#version 330 core
layout (triangles) in;
layout (line_strip, max_vertices = 64) out;

uniform mat4 _worldToClip;
uniform vec4 minorColor;
uniform vec4 majorColor;

uniform float majorInterval;
uniform float interval;
uniform int multiplier = 100; // used for float precision while searching for major contours

in vec3 _normal[];
out vec4 color;
out vec3 _positionT;
out vec3 _normalT;


int countAB;
int countBC;
int countCA;

float interpolate(float start, float end, float t){
    return start + t * (end - start);
}
vec3[100] computeIntersections(vec4 start, vec4 end, inout int counter){

    vec3[100] points;
    counter = 0;
    
    float minValue = min(start.z, end.z);
    float maxValue = max(start.z, end.z);
    float startValue = minValue - mod(minValue, interval);
    float delta = maxValue - minValue;
    float nLevels = floor((maxValue - startValue) / interval);

    for(int i = 0; i <= nLevels; i++){
        float value = startValue + (i * interval);
        if (value < minValue || value > maxValue)
            continue;
        float ratio = (value - minValue) / delta;
        if(start.z > end.z) ratio = 1 - ratio;
        vec3 point = vec3(interpolate(start.x, end.x, ratio), interpolate(start.y, end.y, ratio), value);
        points[counter] = point;
        counter++;
    }
    return points;
}
void drawContours(vec3 normalVector, vec3[100] start, vec3[100] end, vec3[100] mid, int startCount, int endCount, int midCount) {
         _normalT = normalVector;

    for (int i = 0; i < 35; i++) {
         if (i >= startCount) break;
         color = mod(round(start[i].z * multiplier), round(majorInterval * multiplier)) != 0 ? minorColor : majorColor;

         _positionT = vec3(start[i].x, start[i].y, start[i].z);
         gl_Position = _worldToClip * vec4(start[i], 1.0);
         EmitVertex();
         _positionT = vec3(end[i].x, end[i].y, end[i].z);
         gl_Position = _worldToClip * vec4(end[i], 1.0);
         EmitVertex();
         EndPrimitive();
     }

     int count = endCount - startCount;
     for (int i = 0; i < 15; i++) {
         if (i >= count) break;
         color = mod(round(end[endCount - 1 - i].z * multiplier), round(majorInterval * multiplier)) != 0 ? minorColor : majorColor;
         _positionT = vec3(end[endCount - 1 - i].x, end[endCount - 1 - i].y, end[endCount - 1 - i].z);
         gl_Position = _worldToClip * vec4(end[endCount - 1 - i], 1.0);
         EmitVertex();
         _positionT = vec3(mid[midCount - 1 - i].x, mid[midCount - 1 - i].y, mid[midCount - 1 - i].z);
         gl_Position = _worldToClip * vec4(mid[midCount - 1 - i], 1.0);
         EmitVertex();
         EndPrimitive();
     }
}
void main() {    
     int counter1;
     int counter2;
     
     vec3[] AB = computeIntersections(gl_in[0].gl_Position, gl_in[1].gl_Position, countAB);
     vec3[] BC = computeIntersections(gl_in[1].gl_Position, gl_in[2].gl_Position, countBC);
     vec3[] CA = computeIntersections(gl_in[2].gl_Position, gl_in[0].gl_Position, countCA);
     
     // A is the lowest point
     if (gl_in[0].gl_Position.z < gl_in[1].gl_Position.z && gl_in[0].gl_Position.z < gl_in[2].gl_Position.z)
     {
         if(countAB < countCA)
            drawContours(_normal[0], AB, CA, BC, countAB, countCA, countBC);
         else
            drawContours(_normal[0], CA, AB, BC, countCA, countAB, countBC);
     } 
     
     // A is the highest point
     else if (gl_in[0].gl_Position.z > gl_in[1].gl_Position.z && gl_in[0].gl_Position.z > gl_in[2].gl_Position.z)
     {
         if(countAB < countCA)
            drawContours(_normal[0], BC, CA, AB, countBC, countCA, countAB);
         else
            drawContours(_normal[0], BC, AB, CA, countBC, countAB, countCA);
     } 
     
     // B lowest & C highest
     else if (gl_in[2].gl_Position.z > gl_in[0].gl_Position.z && gl_in[0].gl_Position.z > gl_in[1].gl_Position.z)
     {
        drawContours(_normal[0], AB, BC, CA, countAB, countBC, countCA);
     }
     // C lowest & B highest
     else
     {
         drawContours(_normal[0], CA, BC, AB, countCA, countBC, countAB);
     }

} 

//fragmentShader
#version 330 core
out vec4 FragColor;
in vec4 color;

void main()
{
  FragColor = color;
}  

The c# code taking an array of float containing vertex coordinates and normal vectors of respective start & end points of each line segment. Here is the layout, where V1 is a vertex, and N1 its corresponding normal vector.

contourPoints = [V1x, V1y, V1z, N1x, N1y, N1z, V2x, V2y, ... ];

        public void UpdateContourLabels()
        {
            var labels = new Dictionary<double, (List<Point3d>, List<Vector3d>)>();

            if (WorkingMesh == null) return;

            int correctForQuadFaces = WorkingMesh.Faces.QuadCount > 10 ? 2 : 1;
            int count = 0;
            if (generatedPointsCount != null)
                count = (int)generatedPointsCount[0] * 6 * correctForQuadFaces;

            for (int i = 0; i < count; i += 6) // Stride is 6 due to two vec3 coming from the geometry shader
            {
                Point3d p = new Point3d
                {
                    X = contourPoints[i],
                    Y = contourPoints[i + 1],
                    Z = Math.Round(contourPoints[i + 2], 2)
                };

                if (!labels.ContainsKey(p.Z))
                {
                    var v = new Vector3d
                    {
                        X = contourPoints[i + 3],
                        Y = contourPoints[i + 4],
                        Z = contourPoints[i + 5]
                    };

                    labels.Add(p.Z, (new List<Point3d> { p }, new List<Vector3d> { v }));
                }
                else
                {
                    bool isClose = false;
                    foreach (var point in labels[p.Z].Item1)
                    {
                        if (point.DistanceToSquared(p) < TerrainData.ContourLabelSpacing.CurrentValue)
                        {
                            isClose = true;
                            break;
                        }
                    }
                    if (!isClose)
                    {
                        var v = new Vector3d
                        {
                            X = contourPoints[i + 3],
                            Y = contourPoints[i + 4],
                            Z = contourPoints[i + 5]
                        };

                        labels[p.Z].Item1.Add(p);
                        labels[p.Z].Item2.Add(v);
                    }
                }
            }
            ContourLabels = labels;
        }

[Grantim]

Hello!

First of all thank you for sharing your code and the article. If I have understood your casea correctly you have following pipeline:

1. Calculate horizontal sections with GPU (and draw it?)
2. Transfering Geometry shader output to CPU (to find label positions)
3. Find label positions
4. Draw labels

Most likely second step is the slowest one, also our experience with geometry shaders shows that using it is quite slow too (compared to rendering without geometry shader)

Suggestion:
You can try using MeshLib function for extracting sections

MeshLib/source/MRMesh/MRExtractIsolines.h

Lines 15 to 16 in 5b62b6d

	/// extracts all plane sections of given mesh
	MRMESH_API PlaneSections extractPlaneSections( const MeshPart & mp, const Plane3f & plane );

It can cover steps 1-2
(we can workl on it to make it faster if need)

---

Newer versions of OpenGl have Storage Buffer, I think it can be used to somehow eliminate 2nd step, but it will require a lot of research and OpenGl coding, and I am not sure that it will work as expected.

[mariuszhermansdorfer]

Thanks for your prompt response! I'll definitely look into extractiPlaneSections and benchmark its speed against the GLSL shader.

You correctly understand what I'm doing. Step 2 is not the bottleneck here, I am using Transform Feedback and it takes 1-2 ms to complete. The slow part is Step 3 where I'm trying to identify the line segments which belong to the same contour line (have equal Z values) and filter them based on distance to other labels.

[mariuszhermansdorfer]

I have trouble understanding how to work with PlaneSections. Below you can see my attempt at retrieving the contour lines but they mostly result in invalid geometry trying to follow existing mesh edges. Any ideas what I'm doing wrong?

struct RawPolylineArrays {
	float* ContourVertices;
	int* ContourLengths;
	int ArrayLength;
};

PINVOKE RawPolylineArrays CreateContours( Mesh* mesh, float interval)
{
	const MeshPart m = MeshPart( *mesh, nullptr );
	Plane3f plane;
	
	std::vector<float> contourCoordinates;
	std::vector<int> contourLengths;

	Box3f box = mesh->getBoundingBox();
	float bottom = box.min.z + interval;
	float contoursCount = (box.max.z - bottom) / interval;

	for ( size_t i = 0; i < contoursCount; i++ )
	{
		plane = Plane3f( Vector3f::plusZ(), bottom + (i * interval) );
		PlaneSections contours = extractPlaneSections( m, plane );

		for ( const SurfacePath& contour : contours )
		{
			for ( const MeshEdgePoint& point : contour )
			{
				Vector3f& pt = mesh->points[mesh->topology.org( point.e )];
				contourCoordinates.push_back( pt.x );
				contourCoordinates.push_back( pt.y );
				contourCoordinates.push_back( pt.z );
			}
			contourLengths.push_back( static_cast<int>(contour.size() * 3));
		}
	}

	RawPolylineArrays result = RawPolylineArrays();
	result.ArrayLength = static_cast< int >( contourLengths.size());
	result.ContourLengths = contourLengths.data();
	result.ContourVertices = contourCoordinates.data();

	return result;
}

[Grantim]

This block should be changed

			for ( const MeshEdgePoint& point : contour )
			{
---				Vector3f& pt = mesh->points[mesh->topology.org( point.e )];
---				contourCoordinates.push_back( pt.x );
---				contourCoordinates.push_back( pt.y );
---				contourCoordinates.push_back( pt.z );
+++				contourCoordinates.push_back( mesh->edgePoint( point ) )
			}

Mesh edge point is internal point on edge, ant not its origin

using MeshEdgePoint = EdgePoint;

MeshLib/source/MRMesh/MREdgePoint.h

Lines 8 to 12 in 5b62b6d

	/// encodes a point on an edge of mesh or of polyline
	struct EdgePoint
	{
	EdgeId e;
	float a = 0; ///< a in [0,1], a=0 => point is in org( e ), a=1 => point is in dest( e )

[mariuszhermansdorfer]

Thanks, this helped solve the issues! It works well now but can't beat the GLSL code on speed. In the top-left corner you can see the measured time in ms. My GLSL routine takes 2-3 ms even for very low interval values. You mentioned there is potential to speed this logic up, how far do you think this could be pushed?

2023-07-19.12-16-13.mp4

[Grantim]

First I would recommend to find sections in parallel

	Box3f box = mesh->getBoundingBox();
	float bottom = box.min.z + interval;
	float contoursCount = (box.max.z - bottom) / interval;

	std::vector<PlaneSections> sections(contoursCount);
	ParalleFor( sections, [&]( size_t i ) 
	{
		plane = Plane3f( Vector3f::plusZ(), bottom + (i * interval) );
		sections[i] = extractPlaneSections( m, plane );
	} );

	size_t coordsSize = 0;
	size_t lengthSize = 0;

	for ( const auto& s : sections )
	{
		lengthSize += s.size();
		for ( const auto& path : s )
			coordsSize += path.size() * 3;
	}

	std::vector<float> contourCoordinates;
	contourCoordinates.reserve( coordsSize );
	std::vector<int> contourLengths;
	contourLengths.reserve( lengthSize );

	for ( size_t i = 0; i < contoursCount; i++ )
	{
		for ( const SurfacePath& contour : contours )
		{
			for ( const MeshEdgePoint& point : contour )
			{
				Vector3f& pt = mesh->points[mesh->topology.org( point.e )];
				contourCoordinates.push_back( pt.x );
				contourCoordinates.push_back( pt.y );
				contourCoordinates.push_back( pt.z );
			}
			contourLengths.push_back( static_cast<int>(contour.size() * 3));
		}
	}

I think it should inprove performace, also we can improve performance of single extractPlaneSections call (we will have a look), but I am not sure that it can beat shader. Anyway I think that whole pipeline with calculating label positions can work faster with extractPlaneSections rather than passing data back and forth and sorting it.

[mariuszhermansdorfer]

Something is wrong with the above code. I commented everything out except for this part and it still results in a memory leak and a crash:

Box3f box = mesh->getBoundingBox();
	float bottom = box.min.z + interval;
	float contoursCount = (box.max.z - bottom) / interval;

	std::vector<PlaneSections> sections(contoursCount);
	ParallelFor( sections, [&]( size_t i ) 
	{
		plane = Plane3f( Vector3f::plusZ(), bottom + (i * interval) );
		sections[i] = extractPlaneSections( m, plane );
	} );

[Grantim]

Oh sure, there is race for plane arg

        Box3f box = mesh->getBoundingBox();
	float bottom = box.min.z + interval;
	float contoursCount = (box.max.z - bottom) / interval;

	std::vector<PlaneSections> sections(contoursCount);
	ParallelFor( sections, [&]( size_t i ) 
	{
		auto localPlane = Plane3f( Vector3f::plusZ(), bottom + (i * interval) );
		sections[i] = extractPlaneSections( m, localPlane );
	} );

[mariuszhermansdorfer]

Thanks @Grantim, switching to parallel execution speeds things up by 5-8x on my machine:

2023-07-19.14-20-24.mp4

I totally agree that this is a more robust way of working with labels but ideally it would be even faster than this. For my application, contour lines are 'just' an analysis mode running on top of other mesh processing operations. The overhead should be as minimal as possible.

For reference, there were some minor bugs in the code, I'm posting a working version below. I also removed the size allocation to vectors, there is no measurable difference in performance but the code looks cleaner:

	const MeshPart m = MeshPart( *mesh, nullptr );

	Box3f box = mesh->getBoundingBox();
	float bottom = std::ceil( box.min.z / interval ) * interval;
	float contoursCount = (box.max.z - bottom) / interval;

	std::vector<PlaneSections> sections( (size_t)contoursCount);
	ParallelFor( sections, [&] ( size_t i )
	{
		Plane3f	plane = Plane3f( Vector3f::plusZ(), bottom + ( i * interval ) );
	sections[i] = extractPlaneSections( m, plane );
	} );

	std::vector<float> contourCoordinates;
	std::vector<int> contourLengths;
	
	for ( const PlaneSections& contours : sections )
	{
		for ( const SurfacePath& contour : contours )
		{
			for ( const MeshEdgePoint& point : contour )
			{
				Vector3f pt = mesh->edgePoint( point );
				contourCoordinates.push_back( pt.x );
				contourCoordinates.push_back( pt.y );
				contourCoordinates.push_back( pt.z );
			}
			contourLengths.push_back( ( uint32_t )( contour.size() * 3 ) );
		}
	}

[Grantim]

Thanks, we will work on extractPlaneSections speedup!

[mariuszhermansdorfer]

Independently of speeding things up, do you have any suggestions on how to generate labels for each contour line?

There will be a user provided spacing value (float). For each contour, I'd like to get points which are: spacing * 0.5, spacing * 1.5, spacing * 2.5, ... away from the contour start. For each of these points, I'd like to have the corresponding mesh normal vector such that I can orient the label correctly.

[Grantim]

I think it also should be done in parallel for each section:

• find all points label positions (int the same ParallelFor loop mb). I would suggest finding it as index (index of SurfacePath, index of MeshEdgePoint) of sections (not to write geometry code to find intermediate positions). You can use tbb::enumerable_thread_specific<T> for acumulating data in ParallelFor (you can find examples in MeshLib code)

• normal can be found by two functions
  

  MeshLib/source/MRMesh/MRMesh.h

  Lines 128 to 130 in 5b62b6d

  	// returns pseudonormal in corresponding face/edge/vertex;
  	// unlike normal( const MeshTriPoint & p ), this is not a smooth function
  	[[nodiscard]] MRMESH_API Vector3f pseudonormal( const MeshTriPoint & p, const FaceBitSet * region = nullptr ) const;

  
  

  MeshLib/source/MRMesh/MRMesh.h

  Lines 119 to 120 in 5b62b6d

  	// computes normal in three vertices of p's triangle, then interpolates them using barycentric coordinates
  	[[nodiscard]] MRMESH_API Vector3f normal( const MeshTriPoint & p ) const;

  
  

  MeshLib/source/MRMesh/MRMeshTriPoint.h

  Line 29 in 5b62b6d

  [[nodiscard]] MeshTriPoint( const MeshEdgePoint & ep ) : e( ep.e ), bary( ep.a, 0 ) { }

[mariuszhermansdorfer]

find all points label positions (int the same ParallelFor loop mb). I would suggest finding it as index (index of SurfacePath, index of MeshEdgePoint) of sections (not to write geometry code to find intermediate positions)

I'm not sure I understand this. How should I accommodate for the spacing value? Should I keep adding individual segments' lengths or is there a better way of finding points located at a specified distance from the beginning of a SurfacePath?

[Grantim]

I think simple adding segments should be fine, unfortunately I don't have any other idea now. There is decimation algorithm for polylines, but casting it to polyline and decimation will most likely take more time.

[Grantim]

We have improved extractPlaneSections a bit. Hope it will speed up things

[mariuszhermansdorfer]

I'll investigate and will keep you posted. On another note, I switched to writing the results directly to the output arrays and parallelized this as well. This shaved off approx. 10-20%.

PINVOKE RawPolylineArrays CreateContours( Mesh* mesh, float interval )
{
	const MeshPart m = MeshPart( *mesh, nullptr );

	const Box3f box = mesh->getBoundingBox();
	const float bottom = std::ceil( box.min.z / interval ) * interval;
	const float contoursCount = ( box.max.z - bottom ) / interval;

	std::vector<PlaneSections> sections( ( size_t )contoursCount );
	ParallelFor( sections, [&] ( size_t i )
	{
	Plane3f	plane = Plane3f( Vector3f::plusZ(), bottom + ( i * interval ) );
	sections[i] = extractPlaneSections( m, plane );
	} );

	size_t lengthSize = 0;
	size_t coordsSize = 0;
	std::vector<size_t> lengthLookup;
	std::vector<size_t> coordsLookup;

	for ( const auto& s : sections )
	{
		lengthLookup.emplace_back( lengthSize );
		coordsLookup.emplace_back( coordsSize );

		lengthSize += s.size();
		for ( const auto& path : s )
			coordsSize += path.size() * 3;
	}

	RawPolylineArrays result = RawPolylineArrays();
	result.ArrayLength = ( uint32_t )lengthSize;
	result.ContourLengths = new int[result.ArrayLength];
	result.ContourVertices = new float[coordsSize];

	ParallelFor( sections, [&] ( size_t i )
	{
	size_t localCoordsStart = coordsLookup[i];
	size_t localLengthStart = lengthLookup[i];
	size_t j = 0;
	size_t k = 0;
	for ( const SurfacePath& contour : sections[i] )
	{
		for ( const MeshEdgePoint& point : contour )
		{
			Vector3f pt = mesh->edgePoint( point );
			result.ContourVertices[localCoordsStart + j++] = pt.x;
			result.ContourVertices[localCoordsStart + j++] = pt.y;
			result.ContourVertices[localCoordsStart + j++] = pt.z;
		}
		result.ContourLengths[localLengthStart + k++] = ( uint32_t )( contour.size() * 3 );
	}

	} );

	return result;
}

[mariuszhermansdorfer]

We have improved extractPlaneSections a bit. Hope it will speed up things

It's roughly 2.5 times faster now! Super cool! For my test mesh and 0.1 interval I'm down to 4ms (screen capture slows it down a bit):

2023-07-19.17-54-56.mp4

Next step is to add the labels.

[mariuszhermansdorfer]

Following your suggestion @Grantim, I have implemented a first prototype of contour label display. It has a very minimal overhead and works as expected:

2023-07-19.22-47-46.mp4

Code below. I'd appreciate if you could have a look and suggest any optimizations.

struct RawPolylineArrays {
	float* ContourVertices;
	int* ContourVerticesLengths;
	int ContourCount;
	float* LabelVertices;
	float* LabelNormals;
	int LabelCount;
};

PINVOKE RawPolylineArrays CreateContours( Mesh* mesh, float interval, bool showLabels, float spacing)
{
	const MeshPart m = MeshPart( *mesh, nullptr );

	const Box3f box = mesh->getBoundingBox();
	const float bottom = std::ceil( box.min.z / interval ) * interval;
	const float contoursCount = ( box.max.z - bottom ) / interval;

	std::vector<PlaneSections> sections( ( size_t )contoursCount );
	ParallelFor( sections, [&] ( size_t i )
	{
	Plane3f	plane = Plane3f( Vector3f::plusZ(), bottom + ( i * interval ) );
	sections[i] = extractPlaneSections( m, plane );
	} );

	size_t lengthSize = 0;
	size_t coordsSize = 0;
	std::vector<size_t> lengthLookup;
	std::vector<size_t> coordsLookup;

	for ( const auto& s : sections )
	{
		lengthLookup.emplace_back( lengthSize );
		coordsLookup.emplace_back( coordsSize );

		lengthSize += s.size();
		for ( const auto& path : s )
			coordsSize += path.size() * 3;
	}

	RawPolylineArrays result = RawPolylineArrays();
	result.ContourCount = ( uint32_t )lengthSize;
	result.ContourVerticesLengths = new int[result.ContourCount];
	result.ContourVertices = new float[coordsSize];

	std::vector<std::vector<float>> labelPositions( sections.size() );
	std::vector<std::vector<float>> labelNormals( sections.size() );

	ParallelFor( sections, [&] ( size_t i )
	{
	size_t localCoordsStart = coordsLookup[i];
	size_t localLengthStart = lengthLookup[i];
	size_t j = 0;
	size_t k = 0;
	std::vector<float> contourLabels;
	std::vector<float> contourLabelNormals;
	for ( const SurfacePath& contour : sections[i] )
	{
		float cumulativeDist = 0;
		// Add first point manually to avoid if statements in the loop
		Vector3f pt = mesh->edgePoint( contour[0] );
		result.ContourVertices[localCoordsStart + j++] = pt.x;
		result.ContourVertices[localCoordsStart + j++] = pt.y;
		result.ContourVertices[localCoordsStart + j++] = pt.z;

		for ( size_t l = 1; l < contour.size(); l++ )
		{
			pt = mesh->edgePoint( contour[l] );
			result.ContourVertices[localCoordsStart + j++] = pt.x;
			result.ContourVertices[localCoordsStart + j++] = pt.y;
			result.ContourVertices[localCoordsStart + j++] = pt.z;

			if ( !showLabels )
				continue;

			Vector3f diff = pt - mesh->edgePoint( contour[l - 1] );

			if ( cumulativeDist + diff.length() >= spacing )
			{
				contourLabels.emplace_back( pt.x );
				contourLabels.emplace_back( pt.y );
				contourLabels.emplace_back( pt.z );

				MeshTriPoint mtp = MeshTriPoint( contour[l] );
				Vector3f normal = mesh->normal( mtp );
				contourLabelNormals.emplace_back( normal.x );
				contourLabelNormals.emplace_back( normal.y );
				contourLabelNormals.emplace_back( normal.z );

				cumulativeDist = 0;
			}
			else
			{
				cumulativeDist += diff.length();
			}
		}
		result.ContourVerticesLengths[localLengthStart + k++] = ( uint32_t )( contour.size() * 3 );
	}
	if ( showLabels )
	{
		labelPositions[i] = contourLabels;
		labelNormals[i] = contourLabelNormals;
	}
	} );

	if ( !showLabels )
		return result;


	size_t labelCount = 0;
	for ( const auto& l : labelPositions )
	{
		labelCount += l.size();
	}

	result.LabelCount = ( int32_t )labelCount;
	result.LabelVertices = new float[result.LabelCount];
	result.LabelNormals = new float[result.LabelCount];

	size_t k = 0;
	for ( size_t i = 0; i < labelPositions.size(); i++ )
	{
		for ( size_t j = 0; j < labelPositions[i].size(); j++ )
		{
			result.LabelVertices[k] = labelPositions[i][j];
			result.LabelNormals[k++] = labelNormals[i][j];
		}	
	}


	return result;
}

[Grantim]

Hello, I have a look at your sample and have ideas for minor changes:

1. Instead of const Box3f box = mesh->getBoundingBox(); better use const Box3f box = mesh->computeBoundingBox();
getBoundingBox takes cached box from mesh's AABBTree, but if there is no one this function will build it that can be slower than just calculating box directly.

if input mesh does not lose caches and you have some operations that require AABBTree getBoundingBox is OK here (extractPlaneSections does not require the tree)

2. Looks like you calculate diff.length() twice: once in if condition and once in else block, length has sqrt inside, so it is prefered to avoid excessive calls, better call it once before if statement.

3. I think you can use std::move here:

	if ( showLabels )
	{
---		labelPositions[i] = contourLabels;
---		labelNormals[i] = contourLabelNormals;
+++		labelPositions[i] = std::move( contourLabels );
+++		labelNormals[i] = std::move ( contourLabelNormals );
	}

or even that way

---	std::vector<float> contourLabels;
---	std::vector<float> contourLabelNormals;
+++	auto& contourLabels = labelPositions[i];
+++	auto& contourLabelNormals = labelNormals[i];

---	if ( showLabels )
---	{
---		labelPositions[i] = contourLabels;
---		labelNormals[i] = contourLabelNormals;
---	}

4. I think std::move can be used in the last copy block:

	size_t k = 0;
	for ( size_t i = 0; i < labelPositions.size(); i++ )
	{
		auto addition = labelPositions[i].size();
		std::move( labelPositions[i].begin(), labelPositions[i].end(), result.LabelVertices[k] );
		std::move( labelNormals[i].begin(), labelNormals[i].end(), result.LabelNormals[k] );
		k += addition;
	}

P.S. I am not sure if this changes will have any significant effort on performance, but I think they can.

[Fedr]

I think it is important to understand which part of your code takes most of the time (e.g. profile the code).

If the time of section extraction is still dominating, then we can think about more accelerations there. For example, you can make use of the fact that all planes are parallel to OXY. And instead of extractPlaneSections call newly added function with custom lambda:

    return extractIsolines( topology, [height, &points=mesh.points] ( VertId v ) { return points[v] - height; }, region );

I guess it can save some 15% time more.

More optimizations inside extractIsolines are possible, but it is important to know whether it is still the most time consuming operation, or we need to switch our focus on other functions.

[mariuszhermansdorfer]

Thanks @Grantim, I ported all of your suggestions except for the last move statement. It doesn't compile for some reason.

@Fedr, as you can see below, extractPlaneSections is still the biggest offender here:

extractPlaneSections: 4.4159 ms
rest: 0.2735 ms
labels: 0.0247 ms

Could you please provide a sample on how to use the extractIsolines function in this context?

[Fedr]

Instead of

Plane3f	plane = Plane3f( Vector3f::plusZ(), bottom + ( i * interval ) );
sections[i] = extractPlaneSections( m, plane );

please try the following code:

sections[i] = extractIsolines( m.topology,
    [height=bottom + ( i * interval ), &points=m.points] ( VertId v ) { return points[v] - height; } );

[mariuszhermansdorfer]

I get the following errors:

[Fedr]

We have added new overload of extractIsolines function only today. Please pull the latest MeshLib.

[mariuszhermansdorfer]

I pulled the latest source but still get these error messages:

[Fedr]

Oh, sorry, please add .z:

sections[i] = extractIsolines( m.topology,
    [height=bottom + ( i * interval ), &points=m.points] ( VertId v ) { return points[v].z - height; } );

[mariuszhermansdorfer]

Now it worked. Here are the results:

extractPlaneSections: 4.249 ms
rest: 0.5207 ms
labels: 0.1264 ms

There is no noticeable difference. If anything it seems to be slightly slower than the previous method.

[Fedr]

I see, ok, we will think what can be optimized more. But it will take a while.

[Fedr]

Today we added one more function in MeshLib:

/// extracts all sections of given mesh with the plane z=zLevel;
/// this function works faster than general extractPlaneSections(...) for the same plane
/// if the sections cross relatively small number of mesh triangles and AABB tree has already been constructed
[[nodiscard]] MRMESH_API PlaneSections extractXYPlaneSections( const MeshPart & mp, float zLevel );

Please give it a try, it might accelerate contour lines extraction from your terrain, if the conditions in the comment are met ( 1. most plane sections are relatively small and 2. cached AABB tree is typically available).

[mariuszhermansdorfer]

Thanks @Fedr. You guys keep delivering awesome quality code! Switching to the new function brings the total execution time down to 1-2ms for my test mesh. This is including finding label positions and copying data to managed code. This is perfect for my application.

2023-07-24.09-00-52.mp4