Is there an established algorithm for moving a set of randomly placed points into an evenly spaced set of points? I have one million points that are distributed randomly across the plane (see diagram below).
I want to move the million points, so they form a 1000×1000 grid, while moving them as little as possible. What mathematical methods can I use to achieve this?

I want to get a grid with distances $h_x$ and $h_y$ as shown here:

 A: I don't know the answer, so will start with a simplification, then attempt to move upwards from that...
Consider first the 1 dimensional case, where all the $y$-coordinates are $0$.
If there is only one point at $(x_0,0)$ then it must be trivially moved to $(0,0)$.
If there are two points at $(x_0,0)$ and $(x_1,0)$, then the first must move to $(0,0)$. The second can stay where it is, giving us $h_x=x_1$.
If there are three points at $(x_0,0)$, $(x_1,0)$ and $(x_2,0)$, then it becomes less straightforward. The first must move to $(0,0)$. The second will move to $(h_x,0)$ and the third will move to $(2h_x,0)$.
The distance moved by the first point is unimportant, because it must be moved to $(0,0)$ regardless.
The distance moved by the second point is $|x_1-h_x|$ and the distance moved by the third point is $|x_2-2h_x|$. Probably the best approach is to minimise the sum of the squares of the distances moved.
This will be given by $S=\left(x_1-h_x\right)^2+\left(x_2-2h_x\right)^2=x_1^2-2h_xx_1+h_x^2+x_2^2-4h_xx_2+4h_x^2$
We want to choose the value $h_x$ that minimises this, so solve $\frac {\partial S}{\partial x_h}=0$
$\frac {\partial S}{\partial x_h}=-2x_1+2h_x-4x_2+8h_x$
$-2x_1+2h_x-4x_2+8h_x=0 \Rightarrow 10h_x=2x_1+4x_2 \Rightarrow h_x=\frac{x_1+2x_2}{5}$
If there are $n$ points at $(x_0,0)$, $(x_1,0)$, ... and $(x_{n-1},0)$, then we have the following generalisation:
The first must move to $(0,0)$. The second will move to $(h_x,0)$, the third will move to $(2h_x,0)$ and the $k$th will move to $((k-1)h_x,0)$.
The sum of the squares of the distances moved is $S=\left(x_1-h_x\right)^2+\left(x_2-2h_x\right)^2+...+\left(x_{n-1}-(n-1)h_x\right)^2$
$S=\Sigma_1^{n-1} (x_k-kh_x)^2=\Sigma_1^{n-1} (x_k^2-2kh_xx_k+k^2h_x^2)$
$\frac {\partial S}{\partial x_h}=-2\Sigma_1^{n-1} kx_k +2h_x\Sigma_1^{n-1} (k^2)$
Solving $\frac {\partial S}{\partial x_h}=0$ gives $h_x=\frac{\Sigma_1^{n-1} kx_k}{\Sigma_1^{n-1} (k^2)}$
I think that this solves the one-dimensional problem.
Moving to two dimensions has one main difficulty. In one dimension it is possible to identify the position of each point in relation to the points immediately before and after it. It doesn't even matter if there are two points with the same value. I can see that ordering the points will be the first complication in moving to two dimensions. Once you achieve that, my method can be easily adapted to two dimensions. 
A: Here is a greedy approach I just thought up.  It may work reasonably.  This feels to me like on of the NP problems that finding the optimum solution is hard but finding a pretty good one is easy.
Run the leftmost part of your grid through the $500^{th}$ leftmost point and the rightmost part of your grid through the $500^{th}$ rightmost point.  Similarly for the top an bottom.  That gives you (maybe a trial) $h_x,h_y$ and locates each grid point.  Now start from the corners and assign each grid point the closest point not yet used.  Work your way into the center.  Having made the assignments you can perturb $h_x,h_y$ and the lower left corner of the grid slightly to see if you can make an improvement in the sum of the distances.  
You could also try some annealing.  Find two points that are assigned to nearby grid points.  See if flipping their assignments improves things.
