Consider getting 3 random normal variables x1,x2,x3 and calculating the vector
That's points uniformly distributed on a sphere. To distribute them uniformly in the ball (the sphere interior), the standard way is to take an U(0,1) distributed variable r, then multiply by cbrt(r). Because the volume element dV is proportional to r^(2) dr (plus one integration), this causes the points to have uniform 3D density f(r)=1.
What we want is the 2D equivalent, but still with 3 variables. It might seem that the only change necessary is to use sqrt(r) instead, then abuse the fact that x3 is independent from the others to ignore it, and the shadow will be right - but that doesn't take into account the x3^(2) in the square root.
Finding E[sqrt(x1^(2)+x2^(2)+x3^(2))/sqrt(x1^(2)+x2^(2))] and multiplying by it won't work, it will just scale the points. Calculating the average when withholding each of x1,x2,x3 also doesn't work. But what we can do is find a _conditional_ expectation, conditioned on the only spherically symmetric value (so sqrt(x1^(2)+x2^(2)+x3^(2))), then multiply by that based on known x1,x2,x3. Pretty sure you can do that. You'd have to integrate r/sqrt(x1^(2)+x2^(2)) over the surface sqrt(x1^(2)+x2^(2)+x3^(2))=r, weighted with the right probability density. I'm sure there's some chi distribution magic you can do there, but I'm not sure how.
Ignoring all that, if you have a spherically symmetric density f(r), you want to integrate f(sqrt(k^(2)+z^(2)))dz from 0 to sqrt(1-k^(2)), where k^(2) = x^(2) + y^(2). This should be equal to a constant for all k. I don't know if that's even possible for any f, but that's another way forward.
Given that you're only interested in the result being roughly ok, I'd try a variant of this algorithm:
1. Pick a random projection angle.
2. Generate a new 3D point in the ball, uniformly.
3. Calculate the minimum distance to other points in the projected plane, and also the minimum distance in 3D.
4. Reject if it's too close to other points based on a weighed combination of these factors, weighing the 2D distance more; repeat from 2. Else add point and repeat from 1. Decrease minimum distance by 5% or so if many attempts in a row fail.
Perhaps do something more after it's finished adding points, e.g. move points slightly in a preferred direction. Or perhaps keep a long list of random projections (at least 50, I'd say), and decide based on those. I'm not sure if what you want is even possible - if it isn't, you'll notice that you can't make any projection better/flatter/more uniform without making others worse.
EDIT: I think this should work: 1. generate a point in the 4-ball, divide coordinates by sqrt(x1^(2)+x2^(2)+x3^(2)+x4^(2)) 2. discard one coordinate.