DEV Community

David Haley
David Haley

Posted on

1

Algorithm edge cases vs notation

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:

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
Enter fullscreen mode Exit fullscreen mode

See below for the input data I was using.

I was using the footprint:

100
000
101
Enter fullscreen mode Exit fullscreen mode

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

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

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
Enter fullscreen mode Exit fullscreen mode

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

In terms of my algorithm, my bug is here:

    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
Enter fullscreen mode Exit fullscreen mode

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,
    )
Enter fullscreen mode Exit fullscreen mode

Image of Timescale

🚀 pgai Vectorizer: SQLAlchemy and LiteLLM Make Vector Search Simple

We built pgai Vectorizer to simplify embedding management for AI applications—without needing a separate database or complex infrastructure. Since launch, developers have created over 3,000 vectorizers on Timescale Cloud, with many more self-hosted.

Read full post →

Top comments (0)

Billboard image

Create up to 10 Postgres Databases on Neon's free plan.

If you're starting a new project, Neon has got your databases covered. No credit cards. No trials. No getting in your way.

Try Neon for Free →

👋 Kindness is contagious

Please leave a ❤️ or a friendly comment on this post if you found it helpful!

Okay