Intersection volume of two 3d alphashapes

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.

Answer 1

[Turlough Hughes]

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

[Zephyr_]

Thank you !

[Sterling Baird]

Are the vertices of the intersecting volume necessarily part of the group of original points? It seems like the "stars" would be missed in the following 2D example:

[Turlough Hughes]

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;

% Preallocations
[V_int_numerical,V_sphere_numerical,N,t] = deal(zeros(numel(N2),1));

% Numerical
for ii = 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)

hold on, plot(N,E_sphere,'-or','LineWidth',lineWidth);
set(gca,'FontSize',fontSize-2)
legend('E_{intersection}','E_{sphere}')