Equidistant points across a cube

I need to initialize some three dimensional points, and I want them to be equally spaced throughout a cube. Are there any creative ways to do this?

I am using an iterative Expectation Maximization algorithm and I want my initial vectors to "span" the space evenly.

For example, suppose I have eight points that I want to space equally in a cube sized 1x1x1. I would want the points at the corners of a cube with a side length of 0.333, centered within the larger cube.

A 2D example is below. Notice that the red points are equidistant from eachother and the edges. I want the same for 3D.

In cases where the number of points does not have an integer cube root, I am fine with leaving some "gaps" in the arrangement.

Currently I am taking the cube root of the number of points and using that to calculate the number of points and the desired distance between them. Then I iterate through the points and increment the X, Y and Z coordinates (staggered so that Y doesn't increment until X loops back to 0, same for Z with regard for Y).

If there's an easy way to do this in MATLAB, I'd gladly use it.

You'll have to define the problem in more detail for the cases where the number of points isn't a perfect cube. Hovever, for the cases where the number of points is a cube, you can use:

l=linspace(0,1,n+2);
x=l(2:n+1); y=x; z=x;
[X, Y, Z] = meshgrid(x, y, z);

Then for each position in the matrices, the coordinates of that point are given by the corresponding elements of X, Y, and Z. If you want the points listed in a single matrix, such that each row represents a point, with the three columns for x, y, and z coordinates, then you can say:

points(:,1) = reshape(X, [], 1);
points(:,2) = reshape(Y, [], 1);
points(:,3) = reshape(Z, [], 1);

You now have a list of n^3 points on a grid throughout the unit cube, excluding the boundaries. As others have suggested, you can probably randomly remove some of the points if you want fewer points. This would be easy to do, by using randi([0 n^3], a, 1) to generate a indices of points to remove. (Don't forget to check for duplicates in the matrix returned by randi(), otherwise you might not delete enough points.)