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 pointsurface 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')
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 
Create scripts with code, output, and formatted text in a single executable document.
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 (view profile)
Amazing function!!
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 (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', 1e6);
Nataliya Perevoshchikova (view profile)
what if
isequal(surface_points, Vl(corresponding_vertices_ID,:)) gives ans=0 as the result assertion failed?
Nataliya Perevoshchikova (view profile)
Thank you Daniel.
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 (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 (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 renumbered 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 (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 (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 (view profile)
@Matyas, I didn't really use references for this, as it just deals with wellknown 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 (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 (view profile)
Hello, nice code. Do you have any references you used to create it? Thanks
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 (view profile)
Hello Shengmin, please send me some example code so I can reproduce the error. (daniel.frisch@kit.edu)
Shengmin Jin (view profile)
Sometimes I got error like 'Vertex 157 is not connected to any face.' What is going on?
John DeBusscher (view profile)
Used this tool alongside https://www.mathworks.com/matlabcentral/fileexchange/22409stlfilereader
and was able to get beautiful plots in no time. Thanks!
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 (view profile)
StellaBou (view profile)
Ghazal (view profile)
Julian Erskine (view profile)
dang vantuan (view profile)
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.
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 (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 (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 (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 (view profile)
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.
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 (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 (view profile)
thanks for sharing
Giovanni Calzone (view profile)
great function.... perfect
Julian Kreimeier (view profile)
Florian Bernard (view profile)
Michael Stritt (view profile)
Great function, I use it a lot for my work and it helps a lot
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 (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 (view profile)
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 (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 (view profile)
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!
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 (view profile)
Great work. Thanks for sharing!
William Morris (view profile)
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 (view profile)
Please send an example to daniel.frisch at student.kit.edu.
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 (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 Kdtree 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 (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 (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.