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

6 次查看(过去 30 天)
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

回答(1 个)

Tim
Tim 2020-10-29
编辑:Tim 2020-10-29
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:

类别

Help CenterFile Exchange 中查找有关 Surface and Mesh Plots 的更多信息

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by