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.

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.

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)

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.

https://en.wikipedia.org/wiki/Disjoint-set_data_structure

Give it a try. Good luck.