algorithm matrix sparse-matrix coordinate-transformation kdtree

Compacting a sparse matrix while preserving outline shape

I am looking for a way to compact a sparse matrix while preserving shape of its outline and (as much as possible) relative distances between non-zero points. So to graphically demonstrate of what I am trying to get:

for a given matrix:

I expect to get the one of the following matrices as a result:

0 1 0    1 1 0    1 1 0
1 1 1    0 1 1    1 1 0
0 1 1    0 1 1    1 0 1

Of course more solutions are possible here and it there is no "better" solution, as long as the algorithm is consistent. The matrices I am working with are 1024x1024 with 15-30k non-zero points. The non-zero points are centres of features I have identified in a 1024x1024 image. I am trying to produce a dense matrix of those features.

I have tried the kD-tree approach where I have calculated n nearest neighbours of each non-zero point. Then, in random order I was visiting each non-zero point and placing its nearest neighbour at appropriate position in a 3x3 matrix with the current point in the centre. This somehow worked but has resulted in still quite sparse matrix and many islands. Taking more than one neighbour has resulted in merging inconsistencies.

Is there some standard way of doing what I am trying to do? Am I missing some pretty standard algorithm out there?

Do you guys have an idea of what would be a good, globally optimisable solution to my problem?

I am implementing this in Python now (with numpy) but I would appreciate any general algorithm suggestions or implementations in any other language...

Thanks in advance!

EDIT 1:

Here is a gist with what I have written so far. I have managed to avoid neighbour conflicts by strict checks of each point's neighbourhood, but I see that some point's are ejected from final matrix. And the whole thing is hell slow...

Input is a CSV in the following format (only cx and cy matter for now):

ParticleID,cx,cy,cz,Volume,Integral 0,Mean 0,Integral 1,Mean 1

EDIT 2: Some comments to Jerry's answer:

1) I guess I did not state my main objectives clearly: the main goal is to make sure that maximum number of points are surrounded by 8 neighbours, without introducing severe distortion to the overall shape of the outline or initial relative placement of points

2) the idea with scaling image seems to be a solution, but how much should I scale it? How to determine optimal scaling factor? I do not want to scale the image to a certain size. I want to remove empty spaces between the non-zero points in my image.

Solution

This is an optimization problem over two objectives: bounding box size and error (in relative distances). As such you've got a Pareto front of possible solutions, each optimizing for a different tradeoff of bounding box size vs. error.

The problem will be much simpler if you put a bound on one of your objectives and just solve for the best you can do on the other objective given that bound. So you could either say "given a particular bounding box size, how can I minimize the error?" or "given a maximum acceptable error, how can I minimize the bounding box size?"

I'm going to assume that you're more likely to have a particular bounding box size in mind than a particular maximum error, so let's start by assuming a bounding box and try to minimize the error of assigning points to unique positions in that bounding box.

You'll want to use an optimization technique that copes well with the discrete flavor of your problem. I'd go with simulated annealing.

Starting with a good seed

Given a bounding box size, you could scale your whole image down to the bounding box size... but some points will end up in the same position. Call this the "naive positioning":

.........................................
.         .         .         .         .              
.         .   2   3 .         .         .              
.    1    .         .         .  4      .              
.........................................
.         .         .         .         .              
.         .         .         .         .              
.         .         .         .         .              
.........................................
.         .         .         .         .              
.         .         .         .         .              
.         .   5     .         .  6      .              
.........................................
.   7     .         .         .         .              
.   8     .         .         .         .              
.         .         .         .         .              
.........................................

Ok, we can deal with this. Define an ordering on your positions. You can do this naively (left-to-right, top-to-bottom) or in a way that better preserves locality (a Hilbert curve). We're going to iterate over those positions and, when we find a position with more than one point, try to distribute the extra points among available empty positions when it makes sense to.

Keep a queue of DisplacedPoints
Keep a stack of EmptyPositions
Iterate over positions in the order you've chosen:
    If this position is empty,
        Push this position onto EmptyPositions
    Else If this position has more than one point,
        Enqueue all but one of these points into DisplacedPoints
    Placement: While there are remaining DisplacedPoints and remaining EmptyPositions,
        Let candidatePoint = DiplacedPoints.Peek
        Let candidatePosition = EmptyPositions.Peek
        Let currentDisplacement = the distance from the current position to candidatePoint's naive position
        Let candidateDisplacement = the distance from candidatePosition to candidatePoint's naive position
        If currentDisplacement >= candidateDisplacement,
            place candidatePoint in candidatePosition
            pop candidatePosition off EmptyPositions
            dequeue candidatePoint off DisplacedPoints
        Else break out of Placement loop
While there are remaining DisplacedPoints and remaining EmptyPositions,
    Dequeue a DisplacedPoint D, pop an EmptyPosition E, and put D into E

If that's too complicated, a simpler approach is to just iterate over the positions; when you encounter a position with too many points, enqueue the extra points, and whenever you encounter an empty position, dequeue a point into it. If you reach the end and there are still points in the queue, iterate backwards over the positions (there shouldn't be any with more than one point anymore) and dequeue points into the empty spaces.

This won't be an optimal solution by any means. It's just a pretty-good first guess to get the simulated annealing off to a good start.

Objective function performance

The naive objective function would be something like the sum (over all pairs of points) of the squared difference between (the distance between points in the original image) and (the distance between points in the scaled down image, times the scaling factor). But with 20k points, that's around 400 million squared differences! We can make this converge to a good solution faster by starting with a dumb but fast objective function to make initial progress, then switching to more expensive (and correct) functions later to get the finer distinctions down.

The fastest, dumbest objective function you could use is O(n): sum over all points their distance from their own naive positioning.

A slower but more correct objective function assigns each point "relevant neighbors": for every point P, we're going to take the points that are in the 3x3 neighborhood of positions surrounding P's naive positioning*. Now you can sum, over all points P, the sum (over all of P's relevant neighbors) the squared difference between the image-distance and the (scaled) assigned-position-distance.

*since the distance from A to B is the same as the distance from B to A, once a point X is one of P's relevant neighbors, we don't have to consider P to be one of X's relevant neighbors. Taking advantage of that gives a 2x speedup.

After those faster objective functions have gotten you closer to optimal, you could use the naive objective function listed above (sum over all pairs squared difference of distances). I'm not sure you actually want this- the naive objective function wants to preserve long-distance structure, but you might only care about preserving distances among neighbors.

EDIT: The overall idea of rescaling should work, but instead of doing a naive rescale and then adjusting it after the fact, another approach is to make the rescaling itself smarter.

Consider using Content-Aware Scaling. Repeatedly alternate between horizontal and vertical rescalings to keep the overall shape as close to the original as possible.

To maintain the spatial relations between nearby points, you can tweak the energy function used in seam generation. Start by assigning the value of 1.0 to positions with a point and 0.0 to positions without points, and apply a gaussian blur. This will cause seam generation to prefer to cut out large empty spaces before it starts to collapse the small spaces between nearby points.

If this is too expensive/ complex, you can try a simpler method:

Take your original example (distinct identities given to points):

Find the mean x and y coordinates of your points. In this example, the mean position is to the lower right of the 3. Logically divide your array up into four quadrants that meet at the mean position.

Each quadrant can now be processed independently. In each quadrant, there are two grid-aligned directions that point towards the center (for example, in the top-left quadrant, down and right point towards the center). We alternate between these two directions, moving points into adjacent empty spaces. If neither direction yields improvement, you're done with that quadrant.

So, taking the top right quadrant as an example, we'll move points down...

Note that we don't push the 4 past the quadrant boundary. Then, we try to push left:

Then we try to push down, but nothing actually moves. Try moving left, and still nothing moves. Top-right quadrant is done. Let's do the other quadrants in the same way:

Now, find the bounding box of points and crop:

2 1
3 4
5 6