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