## Alternatives to Delaunay triangulation

David Winthrop

### David Winthrop (view profile)

さんによって質問されました 2019 年 11 月 4 日

### David Winthrop (view profile)

さんによって コメントされました 2019 年 11 月 7 日
Hello,
I am in need of calculating the volume of a shape that is given as a cloud of points. I thought a first step would be to triangulate using delaunay but the results are not satisfactory. I have attached some images to show what is unsatisfactory. I zoomed in to look at a part of the cloud that is kind of shaped like the cone of an airplane. I have provided two angles so a sense can be gotten. The super elongated triangles do not seem like they will be useful in getting me surface area and volume.
Does anyone know of alternatives that would produce better results?
Thanks.

David Wilson

### David Wilson (view profile)

on 5 Nov 2019
You could try alphashape. Further details on how to extract its volume are given in:
or you could try convhulln.
I'll try to regenerate your cone shape. One could, if they bothered, compute the exact volume for comparison.
%% Make a nose cone
npts = 1e3;
a = 1; % change this to pi/2
l = a*rand(npts, 1);
phi = rand(npts,1)*2*pi;
r = sin(l);
x = cos(l);
y = r.*sin(phi); z = r.*cos(phi);
plot3(x,y,z,'.');
axis equal
grid on
Now compute the convex hull and the volume.
pts = [x,y,z];
[k,vol] = convhulln(pts);
For a sanity check, if you re-run the code above with
a = pi/2 % half a sphere
I get vol = 2.0593 which is not too far off the exact 4/6 pi = 2.0944.
David Winthrop

### David Winthrop (view profile)

on 5 Nov 2019
Richard, it is just the zoom of the MATLAB plotter that cuts off part of the figure when you zoom in. The shape is as closed as a cloud of points can be. I will make an "answer" that has the data in a text file and provide the code I am using to plot it.
David Winthrop

### David Winthrop (view profile)

on 5 Nov 2019
David,
Thanks for the input on how to calculate the volume, but won't convex hull exclude areas of concavity?

サインイン to comment.

## 4 件の回答

### Richard Brown (view profile)

on 7 Nov 2019
Edited by Richard Brown

### Richard Brown (view profile)

on 7 Nov 2019

As David Wilson suggested, use alphashape:
shp = alphashape(U(:, 1), U(:, 2), U(:, 3));
plot(shp)
disp(shp.vol)
I get a volume of 3.59e-8. The picture is here:
Or, I just noticed that the commands are different in R2019a (I did the previous with R2018a by accident):
shp = alphaShape(U(:, 1), U(:, 2), U(:, 3));
plot(shp)
disp(shp.volume)

#### 0 件のコメント

サインイン to comment.

### David Winthrop (view profile)

on 5 Nov 2019
Edited by David Winthrop

### David Winthrop (view profile)

on 5 Nov 2019

I have attached the point cloud as a file Coords.txt. Below is the code I use to plot it. I have found that the boundary code works well but still has some very large and skew triangles in it. It is an improvement. Any further suggestions would be welcome.
fl = fopen('Coords.txt');
T = textscan(fl,'%f %f %f','Delimiter',',');
X = T{1};
Y = T{2};
Z = T{3};
B1 = boundary(X,Y,Z,.8);
figure
plot3(X,Y,Z,'*')
Ax1 = gca;
xlabel('X')
ylabel('Y')
zlabel('Z')
figure
trisurf(B1,X,Y,Z)
Ax2 = gca;
set(findall(Ax2,'type','patch'),'edgecolor',[.2,.2,.2])
xlabel('X')
ylabel('Y')
zlabel('Z')
Thanks

#### 0 件のコメント

サインイン to comment.

### darova (view profile)

on 6 Nov 2019

Volume can be found as summation volume of pyramids. Here is an idea:
clc,clear
x = Coords(:,1);
y = Coords(:,2);
z = Coords(:,3);
% move to origin
x = x-mean(x);
y = y-mean(y);
z = z-mean(z);
% convert to spherical system
[t,p,r] = cart2sph(x,y,z);
tri = delaunay(t,p);
trisurf(tri,x,y,z)
% indices of triangle
i1 = tri(:,1);
i2 = tri(:,2);
i3 = tri(:,3);
% vectors of triangle base
v1 = [x(i1)-x(i2) y(i1)-y(i2) z(i1)-z(i2)];
v2 = [x(i1)-x(i3) y(i1)-y(i3) z(i1)-z(i3)];
A = 1/2*cross(v1,v2,2); % surface of a triangle
V = 1/3*dot(A,[x(i1) y(i1) z(i1)],2); % volume of a triangle
V = sum(abs(V));

David Winthrop

### David Winthrop (view profile)

on 6 Nov 2019
Something is wrong with your volume calculations. If I do this:
fl = fopen('Coords.txt');
T = textscan(fl,'%f %f %f','Delimiter',',');
X = T{1};
Y = T{2};
Z = T{3};
[B1,V1] = boundary(X,Y,Z,.8);
V1
V2 = VolOfTriangulation(B1,X,Y,Z)
The function VolOfTriangulation is simply your code for calculating the volume based on triangulation.
V1 = 3.746596463335299e-08
V2 = 5.908642403710292e-08
darova

### darova (view profile)

on 7 Nov 2019
Maybe it's because of wrong triangulation in this place:
alphashape should be cool in this case then
darova

### darova (view profile)

on 7 Nov 2019
It's because of converting cartesian coordinates to spherical
[t,p,r] = cart2sph(x,y,z);
Geometrically -pi and +pi are the same, but delaunay interprets it as different numbers.
That is why there is a gap between -pi and +pi
Don't know how to resolve it

サインイン to comment.

### David Winthrop (view profile)

on 7 Nov 2019

Thanks everyone!
I found a file exchange function that does a very nice job of displaying the mesh so that it matches the mesh exported from a Finite Element suite. However, I am still in need of calcuating the volume accurately. If I plot my points using that function and use the volume calculated by dorova I do not get the same value that the FEM code gets.
I get 5.901791412991205e-08 whereas the software says the volume is 3.739248319033497e-08. Since I am trying to use MATLAB to track the volumes I need to match the software volume as closely as possible. Does anyone else know of a good method for this?

darova

### darova (view profile)

on 7 Nov 2019
I get 3.7336e-08 using my script. How is this possible?
David Winthrop

### David Winthrop (view profile)

on 7 Nov 2019
Hello Richard. I did see yours and I had tried it before. If you will notice, there are some large triangles near the circular opening as the method you used is covering a convex section. If you feel like it, please download the code I linked to and check it out. It is very nice.
David Winthrop

### David Winthrop (view profile)

on 7 Nov 2019
darova, I am not sure. I will double check .

サインイン to comment.

Translated by