Find the peak of each islands in sparse matrix

I have a sparse matrix that contains several islands of unknown size. I would like to find the highest peak of each islands. Consider this matrix as an example:

``````0 0 1 0 0 0 0
0 0 1 2 1 0 0
0 0 0 3 2 1 0
0 1 0 0 0 0 0
0 2 3 4 0 1 0
0 0 1 1 0 1 0
0 0 0 0 0 1 0
``````

In this case I would like to get the position of 3 at (3, 2), the position of 4 at (3, 4) and the position of 1 at (5, 4), (5, 5) or (5, 6). Note: I consider 4-connected neighbours.

Until now I came up with the solution that scans each pixel and if it is not zero it starts a flood fill from that pixel. The flood fill paints the pixels to zero keeping track the maximum value visited. This algorithm works well but I was wondering if there are any faster solution?

In my target platform (iOS) I have access to different vector processing frameworks like BLAS, LAPACK, vDSP that provide some functions for fast finding the maximum value in a vector (or matrix). Maybe one could implement a faster solution with the help of these functions?

• so do you actually store the 0s in your data structure or are you using something else?
– jk.
May 8, 2013 at 11:18
• To continue with @jk. comment - there are several different methods of storing a sparse matrix - other structures may remove the 'for each pixel in image: if pixel is non-zero' loop because you only have non-zero pixels in the stored data.
– user40980
May 8, 2013 at 13:23
• Unfortunately I do not have control over the input format or data structures as I am working on an image processing application and the raw input arrives in the format above. To convert it to any other format I should scan at least once the whole matrix. For the same reason I do not expect faster solution than O(n) but maybe with vector processing instructions the algorithm could be fastened with a constant multiplier.
– MrTJ
May 9, 2013 at 7:52
• how to optimize a specific algorithm using vector instructions is likely to be more a stack overflow question I think
– jk.
May 9, 2013 at 8:05

5 Answers

I once had to solve the exact same problem, as part of commercial image processing software. In the end I opted for a scan-line implementation, because it reads the matrix only once.

Basically, you consider a single row of the matrix and keep track of all 'events' on this scan line. Events are left and right edges of islands. Note that the same island might have multiple left and right edges on the same scan line. Also note that what seems to be 2 seperate islands on your scan line, might turn out to be the same island later.

You also keep a record of each known island with its highest value.

You start with the top row, initialize the events. And then go down one row at a time, each time revealing one more row of the matrix and updating the events. Next to updating the position of the existing events, you need consider also the following topology-changes:

1. Start (top) of new land: a new island record, and new left and right edges are inserted in the scan line.
2. End (bottom) of a piece of land: a left and right event are removed from the scan line.
3. Top of a lake/sea: An island that splits into 2 seperate peninsulae of the same island (2 new events)
4. Bottom of a lake/sea: two peninsulae merge into one. If they didn't already belong to the same island records, then these 2 island records need to be merged.
• This is a very interesting and different approach. I surely will try it and let you know the results! Btw I am also working on an image processing application :)
– MrTJ
May 9, 2013 at 7:50
• If your matrix is somewhat localized (eg. not like a random 0-1 pattern), then this will be significantly faster than fill-based algorithms. May 9, 2013 at 10:21

If I understand correctly, it would seem that you're doing it in the most efficient way. Let me try to summarize your algorithm:

``````Main:
For each pixel in image:
If pixel is non-zero:
Call flood subroutine with pixel.
End
End

Flood subroutine:
Acquire all neighboring pixels of passed pixel of the same height.
Peak = true
For each pixel
If pixel has neighboring pixel with higher height
Peak = false
Call flood subroutine with pixel with higher height
End
End
If Not Peak
For each pixel
Set pixel height to 0
End
Else
Add peak to list
End
``````

This means that although technically you're passing every pixel in the image to perform this flooding subroutine, if a lot of the pixels share the same height, you will later hit a lot of 0 height pixels that were set by the subroutine.

You could probably even optimize this a bit to both select all pixels of the same height as well as determine if there are neighboring pixels with higher height in the same algorithm to avoid having to perform that check afterwards.

The only problem with this algorithm is that it may be rather memory heavy, since it uses recursion to hold the pixels with the same color. You could probably improve on this a bit by selecting the neighboring pixels with higher height, setting currently selected pixels to 0 if found, then freeing the resources prior to calling the flood subroutine for each pixel with higher height.

I hope that helps.

• Also check out non-recursive flood fill versions such as this and this, which should more efficient. You may be able to use GCD for this as well to parallelize it. May 8, 2013 at 11:32
• Thank you very much for your answer, it made me think further. As far as I understand your version supposes that there is only one "mountain" per island, i.e. only one local maxima. I did not use this information in my implementation (the "If pixel has neighboring pixel with higher height" condition) but I called the flood subroutine for each neighbouring pixel that has a non-zero height. But in fact it seems my input data contains only one local maxima per island (I have to analyse it carefully) and in this case I can use your optimisation.
– MrTJ
May 8, 2013 at 12:07
• However if you have any idea how could I utilise vector processor features for accelerating the code, please let me know. Vector processor has very fast (apparently atomic) instructions for getting for example the maximum value and position of the maxima in un scalar vector or matrix.
– MrTJ
May 8, 2013 at 12:14
• @MrTJ I think that's not possible. It would be good for determining the max and min, but it couldn't tell you things like local peaks. If anything, perhaps it'd be possible to get good approximations by analyzing sections of the map at a time for the min and max. It really depends on how efficient / precise it has to be. Also, the subroutine gets called for every higher pixel, and since the subroutine sets pixels to 0, all higher pixels not yet set to 0 will still be called (which logically could contain many such mountains per island).
– Neil
May 8, 2013 at 12:39

If I understand your current algorithm correctly it is:

``````loop through all n values in your matrix looking for non 0 ones in O(n)
for each of island, flood fill over the m non zero elements and set them to 0 in O(m)
``````

given you matrix is sparse i.e. n >> m then this algorithm is going to be O(n) complexity.

We can get better complexity by using a data structure that doesn't store the 0s. for instance if you added your non 0 data to a dictionary with the key being the index of the data (i,j) then we can still get any data in O(1) and rely on a key not being present meaning the value is 0.

``````we can then for each non zero element we haven't looked at yet, run the flood fill
and remove every key we look at from the dictionary in a complexity of O(m).
``````

Now your question doesn't ask for a lower algorithmic complexity solution but a faster one, is this new algorithm going to help?

Well, it depends:

• If your matrix is big enough and sparse enough then this should outperform your current implementation. For a smaller matrix you original algorithm may well end up being faster as it could well have smaller constant factors.
• Can you even change your data structures? (answer here appears to be no)
• Thank you for your hints. Unfortunately I do not have control over the input data structure as my algorithm is used in a video processing application and the input matrix is a preprocessed video signal. So in all case I have to seek through the whole matrix. I hope some speed up in the fact the I could use a vector processor, so for example MAX(vector) operations could be atomic. If I could use these operations somehow, maybe I could come up a solution that still has O(n) complexity but might be constant times faster then the current one.
– MrTJ
May 9, 2013 at 7:37
• yes constructing a dictionary from your data then would be O(n) anyway.
– jk.
May 9, 2013 at 7:58

Sounds like a variation of blob coloring. (Google is your FRIEND!)

You're going to make a list of blobs, initially empty. Scan for a nonzero value. When you find one, do 4-way-connectivity to find the blob pixels, and call the "color" the highest value. When you run out of 4-way-connected pixels, restart your scan to find a new blob.

I believe this is an example of Disjoints Sets where
Island: Set
Peak: Great Number of the set (or "Set's parent")

I would go with "Disjoint Sets Forest" where each sets keeps a pointer to the "peak" of the island.

Give it a try. Good luck.