Algorithm edge cases vs notation
Tags: algorithms, debugging, python, programming
Published: 2024-02-11 07:57:04

I just spent 2.5 hours debugging an infinite loop in this algorithm, caused by an edge case (center point of structuring element is excluded). I'd missed a single but crucial requirement…

Here is the [looping code](https://github.com/dchaley/deepcell-imaging/blob/515b5040e01d5005d2122dffccf77aa5078c35dc/benchmarking/h_maxima/benchmark_utils.py#L17):

```
def opencv_reconstruct(
    marker: np.ndarray, mask: np.ndarray, kernel: np.ndarray, anchor: tuple = (-1, -1)
):
    # .item() converts the numpy scalar to a python scalar
    pad_value = np.min(marker).item()

    while True:
        expanded = cv2.dilate(
            src=marker,
            kernel=kernel,
            borderType=cv2.BORDER_CONSTANT,
            borderValue=pad_value,
            anchor=anchor,
        )
        expanded = np.fmin(expanded, mask)

        # Termination criterion: Expansion didn't change the image at all
        if (marker == expanded).all():
            return expanded

        marker = expanded
```

See below for the input data I was using.

I was using the footprint:
```
100
000
101
```

Crucially, the center point is false– which I gather is a bit unusual in these operations, every example I've ever seen has a true in the center.

Anyhow, I was confused, because Vincent'93 says the algorithm terminates but mine never converged. Let's say the above is an implementation of the following:

![Screenshot of Vincent'93 parallel reconstruction algorithm](https://dev-to-uploads.s3.amazonaws.com/uploads/articles/s5kgg1nw8nrj2gvw76sa.jpeg)

opencv for its part takes the footprint very literally. If you didn't ask for the center point, it won't give you the center point.

But Vincent'93 actually *defines* reconstruction as including the center point:

![Vincent '93 dilation step of parallel reconstruction](https://dev-to-uploads.s3.amazonaws.com/uploads/articles/g5ca33dryz7bb694ugj1.png)

In English this is saying,
- for each pixel in the image,
- the dilated pixel value is equal to the maximum of:
  - the pixel value of each point q, where q is in the neighborhood of p, or p itself

It's that last little `∪` that does the trick…

If we change the algorithm to select the center point (hackity hack)
```
kernel[1, 1] = True
```

then … voilà, no more infinite loop (and we get the correct answer to boot).

In terms of my algorithm, my bug is [here](https://github.com/dchaley/deepcell-imaging/blob/515b5040e01d5005d2122dffccf77aa5078c35dc/benchmarking/h_maxima/fasthybridreconstruct.pyx#L105):

```
    for footprint_row_offset in range(-footprint_center_row, footprint_center_row + 1):
        for footprint_col_offset in range(
            -footprint_center_col, footprint_center_col + 1
        ):
            # Skip this point if not in the footprint.
            footprint_x = footprint_center_row + footprint_row_offset
            footprint_y = footprint_center_col + footprint_col_offset
            if not footprint[footprint_x, footprint_y]:
                continue
```

I skip if its not in the footprint– I should also ensure the center is always kept.

That union symbol really mattered 😎  all because of an odd footprint that doesn't include the center.


Input data:
```
    seed = np.array(
        [
            [
                49,
                -1,
                -54,
                -16,
                55,
                -20,
                56,
                41,
                -5,
                -28,
            ],
            [
                -18,
                56,
                -16,
                56,
                -26,
                56,
                41,
                -6,
                7,
                7,
            ],
            [
                -17,
                37,
                56,
                -16,
                56,
                41,
                -35,
                -69,
                -15,
                95,
            ],
            [
                -31,
                -17,
                -18,
                -2,
                -19,
                18,
                41,
                45,
                95,
                -0,
            ],
            [
                31,
                7,
                42,
                -22,
                10,
                -52,
                45,
                95,
                4,
                28,
            ],
            [
                31,
                42,
                31,
                10,
                31,
                -50,
                -21,
                45,
                95,
                2,
            ],
            [
                42,
                31,
                31,
                44,
                86,
                -54,
                107,
                14,
                45,
                17,
            ],
            [
                -103,
                45,
                31,
                66,
                -63,
                66,
                40,
                -64,
                40,
                49,
            ],
            [
                -46,
                41,
                42,
                -28,
                107,
                45,
                -36,
                60,
                -26,
                18,
            ],
            [
                41,
                -29,
                -37,
                45,
                23,
                48,
                40,
                -26,
                4,
                -40,
            ],
        ],
        dtype=np.int16,
    )
    mask = np.array(
        [
            [
                49,
                -1,
                -54,
                35,
                55,
                -20,
                122,
                46,
                17,
                -28,
            ],
            [
                37,
                57,
                -16,
                77,
                -26,
                109,
                52,
                -6,
                7,
                7,
            ],
            [
                36,
                91,
                56,
                27,
                56,
                94,
                -35,
                -69,
                -15,
                126,
            ],
            [
                -31,
                -17,
                -18,
                -2,
                -19,
                18,
                41,
                99,
                112,
                -0,
            ],
            [
                31,
                7,
                44,
                -22,
                92,
                -52,
                102,
                123,
                4,
                28,
            ],
            [
                50,
                82,
                81,
                10,
                61,
                -50,
                -21,
                45,
                95,
                2,
            ],
            [
                42,
                105,
                31,
                52,
                86,
                -54,
                124,
                14,
                72,
                17,
            ],
            [
                -103,
                97,
                31,
                107,
                -63,
                126,
                60,
                -64,
                124,
                93,
            ],
            [
                -46,
                48,
                45,
                -28,
                107,
                86,
                -36,
                113,
                -0,
                18,
            ],
            [
                51,
                125,
                -37,
                66,
                23,
                48,
                45,
                -26,
                4,
                -40,
            ],
        ],
        dtype=np.int16,
    )
```