File Exchange

## point2trimesh( ) — Distance Between Point and Triangulated Surface

version 1.0.0.0 (38.2 KB) by Daniel Frisch

### Daniel Frisch (view profile)

Distance between a point and a triangulated surface in 3D. Can insert the nearest point as vertex.

Updated 25 Sep 2016

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')

### Cite As

Daniel Frisch (2020). point2trimesh( ) — Distance Between Point and Triangulated Surface (https://www.mathworks.com/matlabcentral/fileexchange/52882-point2trimesh-distance-between-point-and-triangulated-surface), MATLAB Central File Exchange. Retrieved .

Daniel Frisch

### Daniel Frisch (view profile)

Hello KAE. point2trimesh() requires that it is given a valid 2.5D surface dataset in the 'Faces' and 'Vertices' parameters. There should be no isolated points in 'vertices' that are not part of any triangle in 'Faces'. If the algorithm comes by such an isolated, unconnected vertex, it gets stuck. This would be a good item for improvement, but until then just remove such unnecessary points before passing the vertices to point2trimesh().

KAE

### KAE (view profile)

Very useful. Can you help me understand why a cloud of points may result in vertices that are not connected to any face? I saw your code in the comments to delete these, but it would be helpful to understand what causes them.

Chaochao

MoHa

### MoHa (view profile)

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)
Regards

Daniel Frisch

### Daniel Frisch (view profile)

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

### MoHa (view profile)

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

### JS (view profile)

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.

ByeongHoon Kim

### ByeongHoon Kim (view profile)

Amazing function!!

Nataliya Perevoshchikova

### Nataliya Perevoshchikova (view profile)

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.

Daniel Frisch

### Daniel Frisch (view profile)

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);

Nataliya Perevoshchikova

### Nataliya Perevoshchikova (view profile)

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

Nataliya Perevoshchikova

### Nataliya Perevoshchikova (view profile)

Thank you Daniel.

Daniel Frisch

### Daniel Frisch (view profile)

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)

Nataliya Perevoshchikova

### Nataliya Perevoshchikova (view profile)

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,:)?

Nataliya Perevoshchikova

### Nataliya Perevoshchikova (view profile)

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

Daniel Frisch

### Daniel Frisch (view profile)

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.

Robin Pillar

### Robin Pillar (view profile)

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

Daniel Frisch

### Daniel Frisch (view profile)

@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.

Daniel Frisch

### Daniel Frisch (view profile)

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);

Matyas Varga

### Matyas Varga (view profile)

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

Shengmin Jin

### Shengmin Jin (view profile)

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

Daniel Frisch

### Daniel Frisch (view profile)

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

Shengmin Jin

### Shengmin Jin (view profile)

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

John DeBusscher

### John DeBusscher (view profile)

and was able to get beautiful plots in no time. Thanks!

Chantel Charlebois

### Chantel Charlebois (view profile)

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.

Chris Fischer

StellaBou

Ghazal

Julian Erskine

dang vantuan

### dang vantuan (view profile)

I have to change the 'linear' to 'parallel'
parser = inputParser;
parser.FunctionName = mfilename;
But the computational time doesn't change.

Daniel Frisch

### Daniel Frisch (view profile)

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.

dang vantuan

### dang vantuan (view profile)

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.

Daniel Frisch

### Daniel Frisch (view profile)

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.

Daniel Frisch

### Daniel Frisch (view profile)

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

### zhenyu (view profile)

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.

Daniel Frisch

### Daniel Frisch (view profile)

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

### zhenyu (view profile)

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

Peter Zacharzewski

### Peter Zacharzewski (view profile)

thanks for sharing

Giovanni Calzone

### Giovanni Calzone (view profile)

great function.... perfect

Julian Kreimeier

Florian Bernard

Michael Stritt

### Michael Stritt (view profile)

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

Nate_52

### Nate_52 (view profile)

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.

Daniel Frisch

### Daniel Frisch (view profile)

@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

### Nate_52 (view profile)

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!

Daniel Frisch

### Daniel Frisch (view profile)

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));

Ernst Schwartz

### Ernst Schwartz (view profile)

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

Bhuvendhraa Rudrusamy

### Bhuvendhraa Rudrusamy (view profile)

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.

Daniel de Malmazet

### Daniel de Malmazet (view profile)

Great work. Thanks for sharing!

William Morris

godfreap

### godfreap (view profile)

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!

Daniel Frisch

### Daniel Frisch (view profile)

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

godfreap

### godfreap (view profile)

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.

Daniel Frisch

### Daniel Frisch (view profile)

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.

Johannes Korsawe

### Johannes Korsawe (view profile)

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.

Johannes Korsawe

### Johannes Korsawe (view profile)

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.