Surface normals of a sphere

Hello! I am trying to calculate (and plot) the surface normals of a sphere. The sphere is a binary 3D image (see code). I want to use the binary image (not the built-in MATLAB functions for spheres), so I can adapt the script later for processing other shapes as well.

I tried the following (see below), but the "normals" matrix are all NaNs. When I change the meshgrid ([X,Y,Z] = meshgrid(x,y,z)) to ([X,Y,Z] = meshgrid(1:64,1:64,1:64)) I do get real values. However, I want to use the x,y,z values I specify in the meshgrid's input arguments.

Can somebody please tell me how to do this? (N.B.: I am not trying to plot the strictly horizontal and vertical surface normals of the single voxels, but I am looking for something like this http://www.opensg.org/projects/opensg/raw-attachment/wiki/Tutorial/OpenSG1/Geometry/rendered_normals.png)

n = 64;
pixelSpacing = [0.1,0.1,0.1];

x = (-(n/2)+0.5)*pixelSpacing(1):pixelSpacing(1):((n/2)-0.5)*pixelSpacing(1);
y = (-(n/2)+0.5)*pixelSpacing(2):pixelSpacing(2):((n/2)-0.5)*pixelSpacing(2);
z = (-(n/2)+0.5)*pixelSpacing(3):pixelSpacing(3):((n/2)-0.5)*pixelSpacing(3);

% Generate sphere
testImage = zeros([n,n,n]);
for k = 1:1:n
    for j = 1:1:n
        for i = 1:1:n
            if sqrt((x(i)^2)+(y(j)^2)+(z(k)^2)) <= 1
                testImage(i,j,k) = 1;
            end
        end
    end
end

[X,Y,Z] = meshgrid(x,y,z);
surface = isosurface(testImage);
p = patch(surface);
normals = isonormals(X,Y,Z,testImage,p);

Answer 1

[Star Strider]

In your e-mail, you mentioned that you wanted me to post the code I wrote, even though it isn’t what your Question originally requested:

n = 20;                             % Facets On Sphere
[Xs,Ys,Zs] = sphere(n);
[Nx,Ny,Nz] = surfnorm(Xs,Ys,Zs);
% Ns = reshape([Nx,Ny,Nz], n+1, n+1, 3);

figure(1)
surf(Xs,Ys,Zs)
hold on
% plot3(Ns(:,:,1), Ns(:,:,2), Ns(:,:,3))
quiver3(Xs,Ys,Zs, Nx,Ny,Nz)
hold off
grid on
axis equal

produces this plot:

The relevant documentation links for sphere, surfnorm, and quiver3 are for the functions I used here. There are a number of other links on those pages that may provide the information you need to do what you want. You might also want to experiment with the lines I commented-out but kept in the code I posted. They relate to lighting of the surface and the way the normals affect that.