Fix neighbour lookup to use correct 4-connectivity in fillminima.c

[sptzmllr]

The previous implementation incorrectly added only diagonal neighbours to the queue, ignoring directly adjacent pixels (up, down, left, right). This could lead to checkerboard-like artifacts. The fix replaces the nested loop with a simple 4-directional lookup using static offset arrays.

In the screenshot you can see the checkboard-like artifact. Sometimes they can be a couple of DN's in a sink, which I believe depends on the starting element in the queue.

Full disclosure: I only tested this in my C++ Implementation because i'am currently not able to run python-fmask on my machine.

[neilflood]

Thanks @sptzmllr
This looks worthwhile, and must have been quite an effort from you to track it down. Well done.

Do you know if this problem was a part of the original Soille & Gratin paper? My recollection is that I got most of the C code directly from the paper, but I no longer have access to that journal, so I can't check.

Most likely it is some sort of transcription error on my part, but it would be nice to know.

I will do some investigating with this implementation, thank you for the full disclosure :-).

Have you done a full C++ implementation of the Fmask algorithm?

[neilflood]

The real error is in line 194, in the test which is supposed to exclude the central pixel (a pixel is never its own neighbour). I messed up the boolean logic, and so also excluded all horizontal and vertical neighbours, leaving only the diagonal neighbours. It should be

if (!((ii == 0) && (jj == 0))) {

I think that the correct solution is to fix that logic. However, you are suggesting we go to 4-way connectivity instead. Is there a strong reason for this? My feeling is that it would be better to leave it at 8-way connectivity, but just fix the central pixel test. I am open to your thoughts, though (as well as being very grateful to you for investigating it).

[sptzmllr]

Do you know if this problem was a part of the original Soille & Gratin paper? My recollection is that I got most of the C code directly from the paper, but I no longer have access to that journal, so I can't check.

The code in the original paper is "preudo C" and only describes/uses the operations you implemented like fifo_add, fifo_first and so on. The neighbor relation is formulated quite openly and also mentions four-way connectivity. Here is the part from Soille & Gratin 1994:

I think that the correct solution is to fix that logic. However, you are suggesting we go to 4-way connectivity instead. Is there a strong reason for this? My feeling is that it would be better to leave it at 8-way connectivity, but just fix the central pixel test. I am open to your thoughts, though (as well as being very grateful to you for investigating it).

The original paper mentions 4-way connectivity, and that is perfectly fine for the job, since you can reach every pixel that way. I also think this should be faster because the code checks to see if it can use the neighboring pixel after the neighbous() function in line 285.

if (img2val == hMax) 

I think with a 4-way connection you have less operations overall and less discarded neighbors. But feel free to benchmark this :)

Have you done a full C++ implementation of the fmask algorithm?

I'm currently writing one for SAGA-GIS, but it's not finished yet. I think it will be done in a few weeks.
What I mean by "my C++ implementation" is the Fill Minima tool, which is based on your C-code. SAGA has a design philosophy where you can call tools from within tools, so I made the Fill Minima routine a tool. And since the licenses were compatible, I ported your code for use in SAGA.

You can find the tool description here I also credited you and cited this project. However: This tool got a major update since the last release and is now quite fast and versatile. I also added a different algorithm for floating-point data (Barnes et al. 2014)

Overall, i believe the 4-way connectivity is faster and produces correct results. The array-based solution is also more maintainable and easier to understand, but that’s just my personal opinion.

[neilflood]

Thank you for your thoughts, that all makes sense.

I managed to get hold of a copy of the paper, and could remind myself of the rest of it. You are right, that paragraph you quoted does indicate that either 4, 6 or 8-way connectedness is fine, and the caption for Fig 4 does suggest 4-way.

I did some timings and tests. Firstly, the speed is much better, at least for the fillMinima operation. Using 4-way saves about 25% on that step. For the whole of the Fmask algorithm, that is only about 2% of the total, so if 8-way were better, I would still use it.

However, I did tests of the whole thing on a couple of Landsat images, and found only tiny differences, whether I used the incorrect diagonal 4-way, your correct 4-way, or full 8-way, so there does not seem to be any value in using 8-way. They are not identical, but you are probably right about the arbitrary queue order being a factor.

I also agree that your hard-coded arrays are a simpler and neater way to represent the set of neighbours.

I have tested your fork directly, too, and it works fine, so I will merge this. Many thanks for your efforts, I really appreciate it.