Create a sphere using surface equation

The surface equation of a sphere is (x-a)^2+(y-b)^2+(z-c)^2-r^2=0.

What I pretend is to create a sphere surface using the equation above. And then get an output matrix with the x,y,z values of the sphere surface nodes...

I don't want to use the sphere function because I intend to create multiple spheres and get the x,y,z values to each one of them. With sphere function I dont know if it would be possible.

Can somebody help me please?

Answer 1

[Adam Danz]

"I don't want to use the sphere function because I intend to create multiple spheres and get the x,y,z values to each one of them."

The solution below contains a function produceSphereCoord() that is based on Matlab's sphere() function but doesn't produce the surface plot. At the top the script below, you can set the number of nodes each sphere will have. Then you can set the number of spheres to create. The 3D coordinates of one sphere is created with [-1,1] normalized units. It's then replicated as many times as requested and the coordinates are stored in a [n-by-3] cell array "xyz" where the columns are values of [x,y,z] coordinates and there's one row for each sphere.

n = 20;           %number of nodes per sphere
nSpheres = 6;     %number of normalized spheres (-1:1)
xyz = cell(nSpheres, 3); %stores you [x,y,z] coordinates for each sphere

% Create normalized sphere
[x, y, z] = produceSphereCoord(n);

% replicate and store in cell array
xyz(:,1) = repmat({x}, nSpheres,1); 
xyz(:,2) = repmat({y}, nSpheres,1); 
xyz(:,3) = repmat({z}, nSpheres,1); 

function [x, y, z] = produceSphereCoord(n)
% n is the number of nodes
theta = (-n:2:n)/n*pi;
phi = (-n:2:n)'/n*pi/2;
cosphi = cos(phi); cosphi(1) = 0; cosphi(n+1) = 0;
sintheta = sin(theta); sintheta(1) = 0; sintheta(n+1) = 0;
x = cosphi*cos(theta);
y = cosphi*sintheta;
z = sin(phi)*ones(1,n+1);
end

You can then scale them as needed with the help of cellfun(). In this example, I scale the 6 spheres to 6 different sizes.

scales = {10 8 6 4 2 .5}'; 
scalesMat = repmat(scales, 1,3); 
xyzScaled = cellfun(@(n,s)n.*s, xyz, scalesMat, 'UniformOutput', false);

To plot the results,

cla()
hold on
grid on
box on
colors = jet(size(xyzScaled,1)); 
arrayfun(@(row)plot3(xyzScaled{row,1},xyzScaled{row,2},xyzScaled{row,3}, 'o', 'color', colors(row,:)), 1:size(xyzScaled,1))
axis equal
view(3)

[Cristiana Abilheira]

Nice, thanks. But I'm kinda out of spherical coordinate world. How can I inform the program where I want the sphere to be centered?

[Adam Danz]

The spheres are created in normalized units. That means they are centered at (0,0) and extend from -1 to +1 in all directions.

To change their center, add the amount of offset to the x,y,z values. For example, if the sphere is centered at (0,0,0) and you want it to be centered at (2,-3,0),

[x, y, z] = produceSphereCoord(n);
x = x+2; 
y = y-3;

To change their size, multiply. For example, if the sphere extends from -1 to +1 (diameter of 2) and you want it to extend from -5 to +5,

[x, y, z] = produceSphereCoord(n);
x = x*5; 
y = y*5; 
z = z*5; 

[Jenny Lee]

Hi,

Thanks for this thread it's really helped me!

Can I ask, how do you decide on the number of nodes per sphere (n)?

For example if I want to make a sphere with a diameter of 4, should i set the number of nodes to be 33, because that's the volume of the sphere? (V=3/4πr^3)

Apologies if the answer is obvious - I'm also new to spherical coordinate calculations.

[Adam Danz]

The number of nodes determines how many points will be used to define the surface of the sphere. If the number of nodes is 20, then X, Y and Z will all be 21x21 points (20 unique points since the first and last coordinate will be the same).

The "scales" variable determines the radius of each sphere. For diameter 4, set "scales" to 2.

[Xiangjie Wang]

Your solutions are really helpful!!!

I wanted to plan some parts of the sphere but I didn't know the meaning of output of the function sphere. Now I am sure that they are:

x = cosphi*cos(theta);

y = cosphi*sintheta;

z = sin(phi)*ones(1,n+1);

By setting some elements in x y z as NaN, I can plot some parts of the sphere now!!!

[Adam Danz]

Thanks Xiangjie Wang. Here's a small demo showing how to plot part of a sphere within a range of theta and phi values. You can adapt it to your needs if it's close to what you're looking for.

% Inputs
thetaRange = [0,pi/2];  % range of theta [min, max]
phiRange = [0,pi/2];    % range of phi [min, max]
center = [5,15,-5];     % center of sphere
radius = 12;            % sphere radius
n = 25;                 % number of verticies per line

theta = linspace(thetaRange(1),thetaRange(2),n);
phi = linspace(phiRange(1),phiRange(2),n)';
cosphi = cos(phi); 
sintheta = sin(theta); 
x = cosphi*cos(theta)*radius+center(1);
y = cosphi*sintheta*radius+center(2);
z = sin(phi)*radius+center(3)*ones(1,n);

surf(x,y,z)
axis equal
grid on

[Xiangjie Wang]

I am not sure whether your script can work for what I am plotting. The partial sphere I am working on is not that continuous, so I just creat x, y, z as three n by n matrices and set the part I don't want to plot out as NaN. It could increase the storage demand because there are too many NaN elements in the matrix. I am wondering whether you have some good idea for what I am doing.

Thank you Adam!

[Adam Danz]

Your approach sounds right for this plot.

[Xiangjie Wang]

Thank you Adam!

Since I am calculating 256000 such surfaces with x,y,z as 2000 by 2000 matrices, it will lead to a really large storage burden. But with your confirmation, I will not feel sorry for such 500GB .mat files ^_^

[Adam Danz]

In that case, it sounds like it could be improved. Replacing unwanted vertices with NaNs is OK in some cases but with lots of objects, it becomes suboptimal since it's storing a lot of unwanted data.

Could you be more specific about the rings you're creating? For example, are you always creating 2 rings around the surface that are orthogonal to each other? Are the widths and orientations of the rings always the same?

[Adam Danz]

Here's a simple example that creates two bands around a circle.

% Parameters
width = 1/10;   % Normalized width of band; ie, 1/10 = 1/10th of the surface
n = 50;         % number of points along each band
m = 5;          % number of points along the width
radius = 5;     % radius of sphere
center = [-2 5 4]; % center of sphere

% Create band along the XY plane
theta = linspace(0,2*pi,n);             % row vector
phi = linspace(-width, width, m)';      % col vector

cosphi = cos(phi); 
sintheta = sin(theta); 
x = cosphi*cos(theta)*radius+center(1);
y = cosphi*sintheta*radius+center(2);
z = sin(phi)*radius*ones(1,n)+center(3);

% Create band along the YZ plane
x2 = sin(phi)*radius*ones(1,n)+center(1); % switched x and z
y2 = cosphi*sintheta*radius+center(2);
z2 = cosphi*cos(theta)*radius+center(3);  % switched x and z

% Plot the first band
surf(x,y,z, 'FaceColor', 'k')

% Plot a second band
hold on
surf(x2,y2,z2, 'FaceColor', 'k')
axis equal
grid on

% Plot center
plot3(center(1), center(2), center(3), 'r*')

[Xiangjie Wang]

Thank you Adam for your reply!

The shape of the partial sphere is not fixed and depends on the calculation results. In the figure I posted, it looks like two rings(with varied width) but in some other situation, it looks more like a partial sphere(the area increases).

Do you have any good ideas for this?

[Adam Danz]

You may want to ask a new question since this is off-topic. If you do so, include a description of how the partial spheres are defined. It may be helpful to provide code for your current approach as well. That may be helpful to create a solution more optimal than your current nan-fill approach.