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

조회 수: 5 (최근 30일)
Matt Southwick
Matt Southwick 2020년 7월 28일
편집: Tim 2020년 10월 29일
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