MATLAB Answers

Alternatives to Delaunay triangulation

21 views (last 30 days)
David Winthrop
David Winthrop on 4 Nov 2019
Commented: M.S. Khan on 31 May 2020
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?
dots_s.png
dots2_s.png
triangs_s.png
triangs2_s.png
Thanks.
  4 Comments
David Winthrop
David Winthrop on 5 Nov 2019
David,
Thanks for the input on how to calculate the volume, but won't convex hull exclude areas of concavity?

Sign in to comment.

Answers (4)

Richard Brown
Richard Brown on 7 Nov 2019
Edited: Richard Brown on 7 Nov 2019
As David Wilson suggested, use alphashape:
U = csvread('Coords.txt');
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):
U = csvread('Coords.txt');
shp = alphaShape(U(:, 1), U(:, 2), U(:, 3));
plot(shp)
disp(shp.volume)

David Winthrop
David Winthrop on 5 Nov 2019
Edited: David Winthrop 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])
hlink = linkprop([Ax1;Ax2],{'xlim','ylim','zlim','cameraposition'});
xlabel('X')
ylabel('Y')
zlabel('Z')
Thanks

darova
darova on 6 Nov 2019
Volume can be found as summation volume of pyramids. Here is an idea:
img1.png
clc,clear
load Coords.txt
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));
  3 Comments
darova
darova 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

Sign in to comment.


David Winthrop
David Winthrop 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?
  5 Comments
M.S. Khan
M.S. Khan on 31 May 2020
David, could you plz explain how you applied trgingulations.
regards,

Sign in to comment.

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by