File Exchange

image thumbnail

point2trimesh( ) — Distance Between Point and Triangulated Surface

version 1.0.0.0 (38.2 KB) by Daniel Frisch
Distance between a point and a triangulated surface in 3D. Can insert the nearest point as vertex.

32 Downloads

Updated 25 Sep 2016

View License

point2trimesh - Distance between a point and a triangulated surface in 3D

The shortest line connecting a point and a triangulation in 3D is
computed. The nearest point on the surface as well as the distance
is returned. The distance is signed according to face normals to
identify on which side of the surface the query point resides.
The implementation is optimized for speed, and depending on your
application you can use linear or parallel computation.

Point insertion functionality
(this feature is experimental and not optimized for speed):
If the function is called with more than two output arguments,
the surface points are included into the given triangulation
and Delaunay conditions are restored locally. If triangles with small
angles occur, additional vertices are inserted to eliminate them
if possible.

Algorithm: From every query point,
- the nearest vertex
- the nearest point on the edges and
- the nearest point on the triangle's surfaces
is calculated and the minimum distance out of these three is returned.

------------------------------------------------------

[ distances, surface_points, normal_vectors,
faces2, vertices2, corresponding_vertices_ID, new_faces_ID ]
= point2trimesh( varargin )

Inputs:
Pairs of parameter names and corresponding values.
Structure arrays are expanded into separate inputs,
where each field name corresponds to an input parameter name.
Parameters:
- 'Faces' (#faces x 3)
Triangulation connectivity list. Each row defines a face, elements are vertex IDs.
- 'Vertices' (#vertices x 3)
Point matrix. Columns are x,y,z coordinates, row numbers are vertex IDs.
- 'QueryPoints' (#qPoints x 3)
Columns are x,y,z coordinates; each row defines a query point. Can be empty.
- 'UseSubSurface' (1 x 1)
Logical. If true (default), the distance to edges and surfaces is only calculated
for faces that are connected to the vertex nearest to the query point.
This speeds up the calculation but if the distance between two opposite parts
of the surface is less than the spacing of the vertices, wrong results are produced.
In the vectorized algorithm, 'SubSurface' is always false.
- 'MaxDistance' (1 x 1)
If the distance between a surface_point and its nearest vertex is within this range,
no new vertex is inserted into the mesh. This helps avoiding
triangles with small angles. (default: 1/10 the smallest inradius)

Outputs:
- distances (#qPoints x 1)
Vector with the point-surface distances; sign depends on normal vectors.
- surface_points (#qPoints x 3)
Matrix with the corresponding nearest points on the surface.
- normal_vectors (#qPoints x 3)
Normal vectors of the surface_points
- faces2 (#faces2 x 3)
Connectivity matrix of the triangulation including the surface_points as dedicated vertices
- vertices2 (#vertices2 x 3)
Point/Vertex matrix of the triangulation including the surface_points as dedicated vertices
- corresponding_vertices_ID (#qPoints x 1)
Vector with the IDs of the vertices corresponding to the query points
- new_faces_ID
Vector with the IDs of the new or modified faces (to give them a different color, for example)

Usage example:
FV.faces = [5 3 1; 3 2 1; 3 4 2; 4 6 2];
FV.vertices = [2.5 8.0 1; 6.5 8.0 2; 2.5 5.0 1; 6.5 5.0 0; 1.0 6.5 1; 8.0 6.5 1.5];
points = [2 4 2; 4 6 2; 5 6 2];
[distances,surface_points] = point2trimesh(FV, 'QueryPoints', points);
patch(FV,'FaceAlpha',.5); xlabel('x'); ylabel('y'); zlabel('z'); axis equal; hold on
plot3M = @(XYZ,varargin) plot3(XYZ(:,1),XYZ(:,2),XYZ(:,3),varargin{:});
plot3M(points,'*r')
plot3M(surface_points,'*k')
plot3M(reshape([shiftdim(points,-1);shiftdim(surface_points,-1);shiftdim(points,-1)*NaN],[],3),'k')

Comments and Ratings (55)

Chaochao

MoHa

Thanks for the feedback. yes correct i want to color each point in the point cloud by its distance to the surface.
I have following three questions:
1- I used the function of Luigi Giaccari for the triangulation, because it enabled me to create my best surface. the function returned only one list (points contained in triangles nx3 array). Which Method returned Vertices and Faces?
2- the distances of each point to all triangles are calculated then the empty distance each points is stored. the function works so well?
so the calculation time becomes very long (for milion points and triangles).
3- Do you have any plan to use classifications such as Supervoxell, Kdtree or Octree in your function? so the calculation time is very much reduced. like the function M2C in Cloudcompare. (https://github.com/CloudCompare/CloudCompare/blob/master/CC/include/DistanceComputationTools.h)
Thank you very much for your answers.
Regards

If I understand you correctly, you want to color each point in the point cloud by its distance to the surface. This can be easily done. 1) Calculate the distance. Therefore pass your surface with the arguments 'Faces' and 'Vertices'; also pass your point cloud as QP=[300000 x 3] matrix after the argument 'QueryPoints'. 2) With the scatter() function, generate a plot with colored points: scatter3(x,y,z,10,c); where x=QP(:,1) etc and c (color) are the distances returned by point2trimesh().

MoHa

Hello, very nice Code, thanks for sharing
I have a meshed surface and a point cloud consists of 300000 points that have only x, y, z coordinates. How can I use this function to represent the distance between surface and the point cloud as a colored figure? Thank you very much for your help and proposal.

JS

Very nice code - Thanks a log.

However, watch out with the default settings. In special cases the closesd point might be a wrong one. Use "'UseSubSurface',false" as additional input to always have a correct output.

Cheers.

Amazing function!!

Thanks Daniel,

Some times it works and in times when it does not work I need to translate the surface that needs to be projected.

OK in this case there was a surface_point on an edge or a face in FV very close to an existing vertex in FV.vertices. Inserting a vertex at this place (for faces2, vertices2) might result in unvaroably small and spiky triangles. So if for one QueryPoint points(k,:) the closest point on the surface surface_points(k,:) is very close to an already existing vertex, then that existing vertex is mentioned in corresponding_vertices_ID(k). The distance that is considered as "very close" here can be configured with the property 'MaxDistance': point2trimesh(FV, 'QueryPoints', points, 'MaxDistance', 1e-6);

what if
isequal(surface_points, Vl(corresponding_vertices_ID,:)) gives ans=0 as the result assertion failed?

Thank you Daniel.

Hello Nataliya,
the first part of your Vl is the same as Vo; the new vertices are just appended. You should make the patch from Fl and Vl. The vector new_faces_ID was just intended for visualization: you can give them a different color and see where something was changed. Here is an extended version of the example provided in the function header that uses point insertion functionality:

FV.faces = [5 3 1; 3 2 1; 3 4 2; 4 6 2];
FV.vertices = [2.5 8.0 1; 6.5 8.0 2; 2.5 5.0 1; 6.5 5.0 0; 1.0 6.5 1; 8.0 6.5 1.5];
points = [2 4 2; 4 6 2; 5 6 2];

[ distances, surface_points, normal_vectors, faces2, vertices2, corresponding_vertices_ID, new_faces_ID ] = point2trimesh(FV, 'QueryPoints', points);
assert( isequal(surface_points, vertices2(corresponding_vertices_ID,:)) )

FV2.faces = faces2;
FV2.vertices = vertices2;
col = zeros( size(faces2,1), 3);
col(new_faces_ID,3) = .5;

pt = patch(FV2, 'FaceAlpha',.5, 'FaceVertexCData',col, 'FaceColor','flat');
xlabel('x'); ylabel('y'); zlabel('z'); axis equal; hold on

plot3M = @(XYZ,varargin) plot3(XYZ(:,1),XYZ(:,2),XYZ(:,3),varargin{:});
plot3M(points,'*r')
plot3M(vertices2(corresponding_vertices_ID,:),'*k') % equivalent to: plot3M(surface_points,'*k')
plot3M(reshape([shiftdim(points,-1);shiftdim(surface_points,-1);shiftdim(points,-1)*NaN],[],3),'k')

view(3)

Daniel,

I did read it wrongly. new_faces_ID corresponds to new created faces. Would it be possible to generate the corresponding_faces_ID for faces in order to patch VI(corresponding_vertices_ID,:) and Fl(corresponding_faces_ID,:)?

Hello Daniel,

Thank you for the code. I am using it to project the vertices1 to the surface (vertices2, faces2).
Surface.vertices=[Vo];
Surface.faces=[Fo];
points= [V(:,1) V(:,2) V(:,3);];
[distances,surface_points,Fl,Vl,corresponding_vertices_ID,new_faces_ID] = point2trimesh(Surface,'QueryPoints', points);

I can see that Fl and Vl were recreated (new vertices and faces were created). Do their ID were re-numbered as well? What is their order (from 1 to n)?
In my understanding corresponding_vertices_ID and new_faces_ID correspond to vertices and faces of Fl and Vl!? If this is so, using them I should be able to patch the projected points into surface
[V]=Vl(corresponding_vertices_ID,:);
[F]=Fl(new_faces_ID,:);
patch('Faces',F,'Vertices',V,'r','FaceAlpha',faceAlpha);

But I cannot do this,
Error using patch
Vectors must be the same length.

Where can be the problem?

Best Regards,
Nataliya

Hello Robin, this is not directly possible as the nearest thing on the triangulated surface can be either a vertex, an edge between two vertices, or a face with three vertices. But you can call point2trimesh with all its output parameters, thereby using its point insertion functionality: for each query point, the nearest point on the surface will become a new vertex (if it wasn't a vertex before). In the argument corresponding_vertices_ID, you will then get the corresponding vertex to each input parameter.

Hi, I need the ID of the nearest surface to each point, is it possible to do this with your code?

@Matyas, I didn't really use references for this, as it just deals with well-known operations of high scool level: find the distance between points, lines, and surfaces. The only "advanced" thing are the barycentric coordinates, which you can find explained on Wikipedia. point2trimesh is just a pragmatic function for prototyping. Surely you can develop a more efficient software for your specific application.

Hello Shengmin, your dataset contains vertices that are not connected to any face. Since your surface represents the convex hull of the all the vertices, you can just remove those unconnected vertices like this:

% Find vertices without face
NVert = size(FV.vertices,1);
has_faces = ismember(1:NVert, FV.faces);
% Remove vertices without face
FV2.vertices = FV.vertices(has_faces,:);
new_ind = cumsum(has_faces);
FV2.faces = FV.faces;
for iVert = 1:NVert
FV2.faces(FV2.faces==iVert) = new_ind(iVert);
end
[distances,surface_points] = point2trimesh(FV2, 'QueryPoints', points);

Hello, nice code. Do you have any references you used to create it? Thanks

Hi Daniel, I've sent you the code and data, please help to take a look. I' really appreciate it. Thanks!

Hello Shengmin, please send me some example code so I can reproduce the error. (daniel.frisch@kit.edu)

Sometimes I got error like 'Vertex 157 is not connected to any face.' What is going on?

Used this tool alongside https://www.mathworks.com/matlabcentral/fileexchange/22409-stl-file-reader
and was able to get beautiful plots in no time. Thanks!

Never mind what I said below, I think the issue was the order of outputs? I called the function [~,~,~,faces2,vertices2]=point2trimesh(input); hoping to get the faces and vertices and I think I ended up getting the vertices and corresponding_vertices_ID.

StellaBou

Ghazal

I have to change the 'linear' to 'parallel'
parser = inputParser;
parser.FunctionName = mfilename;
parser.addParameter('Faces' ,[] ,faceChk);
parser.addParameter('Vertices' ,[] ,vertChk);
parser.addParameter('QueryPoints' ,[] ,pointChk);
parser.addParameter('MaxDistance' ,[] ,distChk);
parser.addParameter('UseSubSurface' ,true ,logicChk);
parser.addParameter('Algorithm' ,'linear',charChk);
But the computational time doesn't change.

For fastest computation, call the function with not more than two output arguments, and set 'UseSubSurface',1, and try out the various algorithms: 'Algorithm','linear' / 'Algorithm','parallel' / 'Algorithm','vectorized' to find out which one is fastest in your specific application.

Thank you for your great job! I have a question. Suppose I calculate the distance of 3000 points to surface, the times calculated is to 3 seconds, How to improve the accelerate of computing or not? for example, I wanna the times of computing with 3000 points decrease to 0.3 seconds. Is it possible?
Thank you very much.

I'm sorry the mesh is delaunized iterativley after the shortest distance to one point has been calculated and the nearest point on the surface has been added as vertex. This behavior might not be optimal and should be changed in a future release.

Everybody who just want the correct distance of any point to your trimesh, set "UseSubSurface"=false (slow calculation without simplifications), and call the function with only one or two output arguments! If your geometry is nearly convex, you can try "UseSubSurface"=true (faster for big trimeshes), and see if it yields the same results.

After the shortest distance points are inserted as additional vertices, the mesh is "delaunized" to remove triangles with unwanted properties (very small and very big angles). With this, the shape of the mesh can change a bit. But shortest distance and nearest point should stay the same. If not, please send me an example.

Generally you should call point2trimesh with more than two output arguments only if you really need the points of shortest distance as dedicated vertices in your mesh.

zhenyu

Thank you for your quick reply! and I know I made a mistake yesterday.
But I'm sorry to bother you to ask another question.I find the results "distances, surface_points" if the output is only this two is different from I have all the outputs.I want to know that why there is the differences and if I only what to know the "distances, surface_points" ,which way should I choose.
when you add the surface_points to the original mesh, I see you do some optimization other than connect the points to the vertices simply.and I really want to know if I want to add the points into the mesh, can I just replace the surfaces_points with the points.
I am very grateful for your answer.

Hi zhenyu! When you use SubSurface, the shortest distance to only a subset of the entire surface is calculated, to save computation time. So the distance with SubSurface can be longer, that's all right.

But it was really a good idea from you to try both. Now you have seen that you can get wrong results with SubSurface - your geometry is probably highly concave. So in order to really get the shortest distance to the entire surface, disable SubSurface.

zhenyu

First, thank you very much for your code. But I have a problem when use your code.
the distances when using SubSurface should be smaller than not using it, isn't it? but I have a bigger one at some points. I don't know why and I want you can help me with this. thank you so much

thanks for sharing

great function.... perfect

Great function, I use it a lot for my work and it helps a lot

Nate_52

Done! I did as you said and collected f.

Thanks for the quick response.

For the case of points on the edge or vertex just the index of one of the faces is fine for what I need, which is just what your code returns.

@Nate_52: Maybe I should have defined struct output arguments, then anything could have been added easily. But consider if the nearest point is on an edge or even a vertex, there is more than one corresponding face.

If you still want one of these faces you can use
[d,p,f] = processPoint(...);
everywhere processPoint() is called and collect the face f in the same way the distance d is collected. Then you can return the nearest faces. (In vectorized algorithm, it's similar but a little more work.)

But doing it by post processing is not a bad idea here.

Nate_52

Nate_52

Nice work! This is super helpful.

One question: I am looking to get the barycentric coordinates of the closest point on the surface. I can get what I need through post processing, but is there an easy way to get the index of the closest face to each point from your function?

Thanks again!

Hello Bhuvendhraa Rudrusamy,
delaunayTriangulation() with 3D points creates a volume of tetrahedrons and no surface of triangles. So you first need to get the convex hull of that volume. Then, some vertices inside the volume may not be part of the hull and have to be removed from the surface net.

% Volume Triangulation with Tetrahedrons
model = rand(50,3);
DT = delaunayTriangulation(model(:,1),model(:,2),model(:,3));

% Construct Surface Net; Remove Unconnected Inner Points
K = DT.convexHull;
vert_ind = 1:size(model,1);
Lia = ismember(vert_ind,K);
model_surf = model(Lia,:);
K_surf = changem(K,1:nnz(Lia),find(Lia));
distances = point2trimesh('Faces',K_surf, 'Vertices',model_surf, 'QueryPoints',model)

% Visualization
trisurf(K,DT.Points(:,1),DT.Points(:,2),DT.Points(:,3), 'FaceAlpha',.2)
axis square; hold on
scatter3Mat = @(XYZ,varargin) scatter3(XYZ(:,1),XYZ(:,2),XYZ(:,3), varargin{:});
h = scatter3Mat(model,'SizeData',50, 'LineWidth',5, 'CData',double(Lia));

I had to add
I=I(:); D=D(:); on lines 575, 576, 634 and 636 to ensure the vectors are column vectors. otherwise, great work!

Given "model" and "probe" made of 3D point clouds. Transforming model to surface via Delaunay Triangulation
DT = delaunayTriangulation(x,y,z);

Looking at the DT properties:
DT =

delaunayTriangulation with properties:

Points: [6940x3 double]
ConnectivityList: [45469x4 double]
Constraints: []
While given input for the function is:
facesx3
verticesx3

Now I have additional column for faces? Can i know how model 3D point cloud can be transformed to surface to suite this script? Thanks for any help.

Great work. Thanks for sharing!

godfreap

Ignore my earlier comment - I didn't have the correct properties called. Very nice, does exactly what it advertises. If you're having problems, just make sure to read the documentation very carefully!

Please send an example to daniel.frisch at student.kit.edu.

godfreap

This seems like a good script, but I think it might have problems with some shapes. I've had a few instances where it seems to call the incorrect surface the closest.

Thanks for the hint! KNNSEARCH with KDTreeSearcher really looks interesting and could be used here. It divides the search space into many buckets. The query point also falls into one bucket, and only points in this bucket and the neighbor buckets have to be considered for distance calculation.

However this requires Statistics and Machine Learning Toolbox, and the Kd-tree for a specific surface has to be grown in advance and must always be passed to the function.

When I need point2trimesh the next time and the profiler tells me it is the bottleneck, I will include this.

Perhaps an idea for improvement: There are new functions for the triangulation class, and also for point clouds (such as "rangesearch" etc.). Maybe one could partially make use of them, as they are fully optimized and very fast. E.g. for reading a distance out of a former computed table instead of calculating it in some sort of loop for each vertex.
Just my 2c.

This just came in very handy!
For long times i have thought about writing such a function by myself, but here, Daniel already did an excellent job, so no further need.

Thanks a lot. Works like a charm.

Updates

1.0.0.0

Full documentation on this page.

1.0.0.0

t

1.0.0.0

Description formatting

1.0.0.0

Description formatting

1.0.0.0

Photo upload

MATLAB Release Compatibility
Created with R2015a
Compatible with any release
Platform Compatibility
Windows macOS Linux