Creating a 3D volume containing an edge interpolated sphere from X,Y,Z coordinates

Hi,

I would like to create a NxNxN volume containing a sphere of radius R, where R can be a fraction of a pixel, so the surface of the sphere is interpolated over multiple voxels in the radial direction.

I wish to embed a sphere of ones (excluding the surface) in a volume of zeros.

How would this be done?

Thanks,

Matt

Here's where I'm starting from.

R = 195;
[X,Y,Z] = sphere(R);
X =round(10*X);
Y =round(10*Y);
Z =round(10*Z);


N = 200;
A = zeros(N,N,N);

P=95;
for i=1:size(X,1)
    for j=1:size(Y,2)
        A(P+X(i,j),P+Y(i,j),P+Z(i,j)) = 1;
    end
end

figure(1);volshow(A);

sizeIn = size(A);
figure(2);slice(double(A),sizeIn(2)/2,sizeIn(1)/2,sizeIn(3)/2);
grid on, colormap gray

Answer 1

[Tim]

This is my best guess as to what you are asking for (updated - misinterpreted your question... hopefully this is closer to what you were looking for).

figure
r = 10;         % Radius
W = r*1.1;      % Volume edge half-length
N = 35;         % Number of voxels along edge

% 3D meshgrid to get radii
[x3, y3, z3] = meshgrid(linspace(-W, W, N), linspace(-W, W, N), linspace(-W, W, N));
rad3 = sqrt(x3.^2 + y3.^2 + z3.^2);

% Voxel of "ones" bounding the surface of the sphere:
volDat = zeros(size(x3));
volDat(rad3 > r - sqrt(3*(W/(N-1)).^2) & rad3 < r + sqrt(3*(W/(N-1)).^2)) = 1;

% Cut-away half to show the results...
volDat(:, ceil(N/2):end, :) = 0;

% Visualization stuff. You will need VOXview from file exchange to
% reproduce.

% The sphere of ones:
pp.edgecolor = 'none';
VOXview(volDat, volDat*0.75, 'colormap', hsv,...
    'patch_props', pp,...
    'CornerXYZ', [x3(1), y3(1), z3(1)],...
    'VoxelSize', 2*W/(N-1),...
    'bounding_box', true);

% The bounding spherical surface:
[XX, YY, ZZ] = sphere(100);
hold on
S = surf(XX*r, YY*r, ZZ*r, 'FaceAlpha', 0.7, 'FaceColor', [0.95, 0.95, 0.95]);
S.EdgeColor = 'none';
hold off
light('position', [0, 1, 5]);

Here is a picture of the result: