I assume you want the samples to be generated from a uniform distribution in the range from 0 to 20.
You could use the reduction method to generate data samples of the desired constraint. That is, generate 10 samples and check whether the samples comply with the desired condition. If not, generate additional samples until the conditions are met.
numberOfSamples = 10;
randomSamples = 20*rand(numberOfSamples,3);
validPoints = reductionMethod(randomSamples);
while(size(validPoints,1) < numberOfSamples)
newPointsNeeded = numberOfSamples - size(validPoints,1);
newPoints = 20*rand(newPointsNeeded,3);
newSampleSpace = [validPoints;newPoints];
validPoints = reductionMethod(newSampleSpace);
function validPoints = reductionMethod(randomSamples)
distanceAmongPts = squareform(pdist(randomSamples));
distanceAmongPts(distanceAmongPts == 0) = 20*sqrt(3);
wrongDist = (distanceAmongPts < 1);
wrongPoint = sum(wrongDist,2) > 0;
randomSamples(wrongPoint,:) = ;
validPoints = randomSamples;
The above algorithm samples points from uniform distribution but I do not think it is scalable to 6400 points. The maximum I saw is around 1000 points. I think some improvements can be done in the above algorithm in the removing section of points less than unit distance apart.
You can think of the 20x20x20 surface equally divided into 8000 cubes. And from the center of those cubes, select 6400 cube centers.
x = 0.5:19.5;
y = 0.5:19.5;
z = 0.5:19.5;
centers = combvec(x,y,z);
centerPoints = randperm(8000,6400);
randomCenters = centers(:,centerPoints);