Take slice through surf plot of multiple spheres and find volume of all spheres dissected

I have a 3d sample (ice) which I am filling with spheres (air bubbles) until it is a known density and plotting using bubbleplot3 from the file exchange ( http://www.mathworks.com/matlabcentral/fileexchange/8231-bubbleplot3/content/bubbleplot3.m ). The function I'm using and sample input is below and the plot that is produced is attached.

What I want to be able to do is take a slice through this sample and find the volume of spheres that would be exposed to the surface if this is done. (I need to calculate the porosity of the sample, so I'm interested in how much liquid would be able to seep into the sample if it were cut).

I've tried a few things like populating a large logical 3d matrix with 'air or ice' values but it just gets too big.

Any ideas welcome.

Cheers

sphere_count(1,2,4,0.025,0.885)

function [ s ] = sphere_count( x,y,z,b_d,d )
%sphere_count Summary of this function goes here
%   x,y,z = dimensions of ice sample, b_d is bubble diameter, d is ice
%   sample density

sample_vol = x*y*z;
air_pc = 1-(d./0.917);
air_vol = sample_vol*air_pc;
sphere_vol = 4*pi*((b_d/10)./2)^2;
number_of_spheres = floor(air_vol/sphere_vol);

P = [(x.*rand(1,1)),(y.*rand(1,1)),(z.*rand(1,1))];
k = 1;

while k < number_of_spheres
    P1 = [(x.*rand(1,1)),(y.*rand(1,1)),(z.*rand(1,1))]; % Make new co-ords
    D = pdist2(P1,P,'euclidean'); % calculate distance between
    if D>(b_d/2) % test maximum distance
        P=[P;P1];  %#ok<*AGROW> % concat
        k=k+1;
    end
end

x=P(:,1);y=P(:,2);z=P(:,3); 

a = 0.02; % Variance
b = b_d; % Mean
R = a.*randn(numel(P(:,1)),1) + b;

bubbleplot3(x,y,z,R);

Answer 1

[Prateekshya]

Hello Tom,

To address your problem of determining the volume of spheres (air bubbles) that would be exposed when slicing through a 3D sample, you can follow a structured approach. The challenge is to calculate the intersection of a slicing plane with the spheres and determine the volume of the exposed parts. Here is a step-by-step guide to achieve this:

• Define the Slicing Plane: Decide on the orientation and position of the slicing plane. For simplicity, let's consider a plane parallel to one of the axes (e.g., the XY plane at a specific Z value).
• Calculate Intersections: For each sphere, calculate whether it intersects with the slicing plane. This involves checking if the distance from the sphere's center to the plane is less than or equal to the sphere's radius.
• Calculate Exposed Volume: For spheres that intersect the plane, calculate the volume of the spherical cap that is exposed. The formula for the volume of a spherical cap is:

where h is the height of the cap (distance from the plane to the top of the sphere) and R is the sphere's radius.

• Sum the Volumes: Sum the volumes of all the exposed caps to get the total exposed volume.

Here is a sample MATLAB function to perform the above calculations:

function exposed_volume = calculate_exposed_volume(P, R, slice_z)
% P: Nx3 matrix of sphere centers
% R: Nx1 vector of sphere radii
% slice_z: Z-coordinate of the slicing plane

    exposed_volume = 0; % Initialize total exposed volume
    
    for i = 1:size(P, 1)
        center = P(i, :); % Center of the sphere
        radius = R(i);    % Radius of the sphere
        
        % Distance from the center to the slicing plane
        distance_to_plane = abs(center(3) - slice_z);
        
        if distance_to_plane < radius
            % Height of the spherical cap
            h = radius - distance_to_plane;
            
            % Volume of the spherical cap
            cap_volume = (1/3) * pi * h^2 * (3*radius - h);
            
            % Add to total exposed volume
            exposed_volume = exposed_volume + cap_volume;
        end
    end
end

You can call this function after generating your spheres:

% Example usage
slice_z = 2; % Z-coordinate of the slicing plane
exposed_volume = calculate_exposed_volume(P, R, slice_z);
fprintf('Total exposed volume: %.3f\n', exposed_volume);

I hope this helps!