Intersection volume of two 3d alphashapes

81 views (last 30 days)
Zephyr_
Zephyr_ on 13 Dec 2019
Hi
I have two alphaShapes in 3d that overlap each other and I would like to know the volume of the overlapsing. I already did it in 2d with polygons (but it was the overlaping area, not volume of course) using the functions intersect and polyarea but I can't find something similar for 3d shapes. I use matlab r2018a.

Accepted Answer

Turlough Hughes
Turlough Hughes on 13 Dec 2019
Let's say you have two shapes that were developed as follows using column vectors as inputs:
shp1=alphaShape(x1,y1,z1);
shp2=alphaShape(x2,y2,z2);
You could determine the points of shp1 that are in shp2 and vica versa the points of shp2 that are in shp1 via the inShape function, then use those to get a new shape, shp3, that represents the intersection:
id1=inShape(shp2,x1,y1,z1);
id2=inShape(shp1,x2,y2,z2);
shp3=alphaShape([x1(id1); x2(id2)], [y1(id1); y2(id2)], [z1(id1); z2(id2)]);
You can then get the volume by:
V=volume(shp3);
3 Comments
Turlough Hughes
Turlough Hughes on 7 May 2021
True. In 3D these stars would be the lines of intersection between the faces approximating the two shapes. I don't see a clear way to do this in alphaShapes.
However, I ran some tests to show the discretisation error as a function of computational time and number of points used. It took ~1 sec to generate two spheres and calculate the intersection volume to within 0.4% of the analytical solution. If you can live with that then the above will do. I also calculated the volume a single sphere with the same number of vertices and this came out at about 30-50% more accurate indicating that there is room for improvement.
% Parameters
R = 1;
d = 0.5;
N2 = 30:10:120;
% see: https://mathworld.wolfram.com/Sphere-SphereIntersection.html
fontSize = 14;
lineWidth = 2;
% Analytical solutions
V_intersection_analytical = (1/12)*pi*(4*R + d)*(2*R - d)^2;
V_sphere_analytical = (4/3)*pi*R^3;
%预先配置
[V_int_numerical,V_sphere_numerical,N,t] = deal(zeros(numel(N2),1));
% Numerical
forii = 1:numel(N2)
tic
[x1,y1,z1] = sphere(N2(ii));
P1 = [x1(:) y1(:) z1(:)];% Points - sphere 1
P1 = unique(P1,'rows');
P2 = P1 + [d 0 0];% Points - sphere 2
N(ii) = size(P1,1);% Number of points
shp1 = alphaShape(P1,1.01);
shp2 = alphaShape(P2,1.01);
id1=inShape(shp2,P1);
id2=inShape(shp1,P2);
P3 = unique([P1(id1,:);P2(id2,:)],'rows');
shp3 = alphaShape(P3,1.01);
V_int_numerical(ii)=volume(shp3);% numerical volume of intersection
t(ii) = toc;% time to generate the two spheres and estimate V
V_sphere_numerical(ii) = volume(shp1);% numerical volume of sphere.
end
E_intersection = 100*abs(V_int_numerical - V_intersection_analytical) ./ V_intersection_analytical;
E_sphere = 100*abs(V_sphere_numerical - V_sphere_analytical) ./ V_sphere_analytical;
Then plot the results
hf = figure(); subplot(1,2,1)
ax1 = plot(t,E_intersection,'-.sk','LineWidth',lineWidth);
set(gca,'FontSize',fontSize-2)
xlabel('Computation time, t / (sec)','fontSize',fontSize)
ylabel('Error, E / (%)','fontSize',fontSize)
subplot(1,2,2)
plot(N,E_intersection,'-.sk','LineWidth',lineWidth);
set(gca,'FontSize',fontSize-2)
ylabel('Error, E / (%)','fontSize',fontSize)
xlabel('Number of Points, N / (-)','fontSize',fontSize)
holdon, plot(N,E_sphere,'-or','LineWidth',lineWidth);
set(gca,'FontSize',fontSize-2)
legend('E_{intersection}','E_{sphere}')

Sign in to comment.

更多的Answers (0)

下载188bet金宝搏


Release

R2018a

Community Treasure Hunt

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

Start Hunting!