algorithms – Area preserving distance optimization

I stumbled across the optimization problem described below and I am looking for how to approach it efficiently. I have the feeling that the problem at hand might be known, but I was not able to find its name or further information on it. So any help to find an efficient algorithm for its solution or a pointer to the name of the problem at hand would be very helpful.

Given $ n times m $ matrix $ M $ of integers from $ (0, w) $ and a set $ P $ of $ w $ points (i.e. tuples of indices in the matrix),
find a permutation $ M’ $ of the values of $ M $, such that the global distance function $ f(M’) = $
$ sum_{i=0}^{n-1}{ sum_{j=0}^{m-1} } { left(i – xleft(P(M'(i, j))right)right)^2 + left(j – yleft(P(M'(i, j))right)right)^2 } $ is minimal.

Explained in a different way, consider an image with $ n times m $ pixels, where each pixel is of one of $ w $ distinct colors. We then want to rearrange the pixels such that all pixels of the same colors from a group around the seed point of the same color. We are allowed the swap the colors of pixels as long as the total number of pixels for each colors stays constant. Below you find an example of the resulting image from a random initial image with a hundred equally likely color as starting point located around randomly pixed seed points. The result was obtained by iteratively picking two points at random and swapping them if the global distance function after swapping is reduced. Doing this for millions of times approaches (at least a local) optimum depicted in the images. However this approach is rather slow for larger values of $ n $ and $ m $. So I am looking for a faster solution to obtain the result. Also a non-discrete solution where polygons would describe the final areas would be appreciated. I was also considering Voronoi tesselations (which can be computed much faster), however they do not work in my case due to the fixed area constraint (i.e. the colored sections should be of equal size in this example).

Area preserving distance optimization