Extract cross sections from a plot of multiple spheres.

Hi All,

I want to extract cross sections from a plot of multiple spheres placed in a cube. I currently generate a plot for the packing of spheres using following code:

[X,Y,Z] = sphere;
for SpNum = 1:NumSpheres
    surf( X*Radius(SpNum)+Center(SpNum,1), Y*Radius(SpNum)+Center(SpNum,2), Z*Radius(SpNum)+Center(SpNum,3), ...
        'FaceColor','r' );
    %shading interp;
    hold on;
end
axis tight; daspect([1 1 1]);

In the above code, each sphere could have different radius and they do not overlap (so the centers are also different).

The above code does not however generate cross sections. I want to extract cross sections similar to what we get from say X-ray CT data, a series of images in the Z-direction. I think 'interp2/interp3' and 'slice' functions are the relevant functions, but I am not sure how to use them to generate the cross sections. I would appreciate if anyone could give pointers or provide some sample code for my problem?

-- Thanks in advance.

Update :

As an alternative, I tried using meshgrid (since interp3 requires as its inputs the variables generated using meshgrid) to generate the grid points followed by the function F(X,Y,Z) as follows (where F is the function which represents the spheres in a 3-D space):

    [X,Y,Z] = meshgrid(1:100,1:100,1:100);
    F = zeros(size(X),'uint8');
    for SpNum = 1:NumSpheres
        F( sqrt((X - Center(SpNum,1)).^2 + (Y - Center(SpNum,2)).^2 + (Z -     Center(SpNum,3)).^2) <= Radius(SpNum) ) = 1;
    end
    surf(F);

followed by:

    z = 1;
    I = interp3(X, Y, Z, X*Radius(SpNum)+Center(SpNum,1), Y*Radius(SpNum)+Center(SpNum,2), Z*Radius(SpNum)+Center(SpNum,3), z, 'spline');
    figure, imshow(I);

I know that interp3 is the function to use since it interpolates the values of the function F(X,Y,Z) which represent the spheres at different locations within a bounded box (say 1:100, 1:100, 1:100). The interpolated values at particular 'z' (= 1, 2, 3... 100) should give me each of 100 cross sections (in the form of 2-D images).

The flaw is in the function F itself, since 'surf' throws an error saying that F should be an array - "CData must be an M-by-N matrix or M-by-N-by-3 array".

Can anyone please help.

Answer 1

[Ramaprasad Kulkarni]

I finally figured it. For the benefit of others, here is the code.

    % A 3-D matrix 'F' which has its value at particular coordinate set to 255 if it belongs to any one of the spheres and 0 otherwise.
    [X,Y,Z] = meshgrid(1:100,1:100,1:100);
    F = zeros(size(X));
    for SpNum = 1:NumSpheres
        F( sqrt((X - Center(SpNum,1)).^2 + (Y - Center(SpNum,2)).^2 + (Z - Center(SpNum,3)).^2) <= Radius(SpNum) ) = 255;        
    end
    % Extract cross sections from F using interp3 function along the z-axis.
    I = zeros(size(X));
    for z = 1:100
        I(:,:,z) = interp3(X, Y, Z, F, 1:100, (1:100)', z, 'spline');
    end
    implay(I,4);

You could test and visualize the output by setting Center (a 3-D vector) and Radius of each sphere (some arbitrary NumSpheres) to some random values. The above code will display a window with cross-sections.

Previously, I was trying to use 'surf' to render the spheres which is not right. To render, you have to use the first code snippet. Another mistake I made was using a row vector for the 6th argument instead of column vector.

Hope this helps.

-- Cheers, Ram.