Find the height of the top surface of a mesh

I have a mesh (triangulation object loaded with stlread) from LiDAR surface data captured from a field. I want to sample the space across an XY grid, where the Z height of each point is the top most surface from the mesh at that point (I am sampling many points and doing RF raytracing). Currently, I am doing this by generating a ray that points directly down from each grid point and calculating the intersection with all the triangles of the mesh (using this).

This approach works, but is quite slow (the mesh is quite large). I imagine there are faster ways to approach this problem, especially considering that I don't need a high degree of accuracy.

Answer 1

[Abhinaya Kennedy]

Hi Rohan,

To speed up the process of sampling the space across an XY grid from the LiDAR surface data captured from a field, you can use an acceleration data structure like a spatial data structure to efficiently find the intersection points. One such data structure is a spatial subdivision structure like a grid, octree, or ‘kd-tree’. These structures can significantly speed up the intersection calculations by reducing the number of triangles that need to be tested for intersection with each ray.

Here's a general approach using a kd-tree:

1. Build a kd-tree: Construct a kd-tree from the triangles of the mesh. This divides the space into smaller regions, allowing for efficient spatial querying.

2. Sampling the space: For each point on the XY grid, instead of casting a ray directly down, you can query the kd-tree to find nearby triangles in the mesh. You can then calculate the intersection of the ray with only these nearby triangles, reducing the number of intersection tests.

3. Adjusting accuracy: You can control the accuracy by adjusting the size of the spatial subdivision. A larger subdivision may reduce accuracy but will speed up the process.

Here's an example of how you can use a kd-tree for this purpose:

% Assuming you have the triangulation object loaded as 'TR'

% Build the kd-tree
kdTree = KDTreeSearcher(TR.Points);

% Define your XY grid (replace with your own grid)
[X, Y] = meshgrid(linspace(min(TR.Points(:,1)), max(TR.Points(:,1)), 100), ...
                  linspace(min(TR.Points(:,2)), max(TR.Points(:,2)), 100));

% Sample the space
Z = zeros(size(X));
for i = 1:numel(X)
    % For each point in the grid, find the nearest triangle
    [idx, ~] = knnsearch(kdTree, [X(i), Y(i)]);
    triangle = TR.ConnectivityList(idx, :);
    % Calculate the Z height based on the intersection with this triangle
    % Perform the intersection calculation here and assign the Z value
    % Z(i) = ...
end 

In this example, `KDTreeSearcher` is used to build a kd-tree from the points of the mesh. Then, for each point in the XY grid, we find the nearest triangle in the kd-tree and calculate the Z height based on the intersection with that triangle.

Using a spatial data structure like a kd-tree can significantly speed up the process of sampling the space across the XY grid, especially when dealing with large meshes. Adjusting the grid resolution and the spatial data structure parameters can help balance speed and accuracy for your specific application.

For further information on the ‘kd-tree’ function, you can refer to the link below.

Create Kd-tree nearest neighbor searcher - MATLAB - MathWorks India

Hope this helps!