Difficulties implementing the Hysteresis step of Canny Algorithm in Halide without define_extern func

Question

• The problem is that when a pixel marked as weak-edge (between the two thresholds) changes to strong-edge (accepted, as described here) it is required to apply the same logic to your connected neighbors recursively (tracking the edges).
• In an imperative language, when changing from weak to strong edge, a Stack could be used to store the position (x,y). Then, at the end, process the neighbors while the stack is not empty, updating the Stack as needed. But, how to implement anything similar in pure Halide, without define_extern func?

I have used this code to hysteresis, but lacks the dynamic recursion and/or stack to do hysteresis on the neighbors when needed, which is what I can't find how to implement:

magnitude = input(x, y);
// mask receives only 0, 1, or 2.
mask(x, y) = select(magnitude > high_threshold, 2/*strong-edge*/, magnitude < low_threshold, 0/*no-edge*/, 1/*weak-edge*/);

// when mask(x,y) == 1 checks the neighbors to decide.
hysteresis(x, y) = select(mask(x, y) == 0, 0,   mask(x, y) == 2, 255,
        mask(x-1, y-1) == 2 || mask(x, y-1) == 2 || mask(x+1, y-1) == 2 ||
        mask(x-1, y) == 2 || mask(x+1, y) == 2 ||
        mask(x-1, y+1) == 2 || mask(x, y+1) == 2 || mask(x+1, y+1) == 2, 255/*weak-to-strong edge*/, 0);

The doubt, is there a way to, with recursion, stack, or anything else, do something like this:

if (hysteresis(x, y) above changes from weak to strong edge, do) {
  hysteresis(x-1, y-1); hysteresis(x, y-1); hysteresis(x+1, y-1);
  hysteresis(x-1, y); hysteresis(x+1, y);
  hysteresis(x-1, y+1); hysteresis(x, y+1); hysteresis(x+1, y+1);
}

Answer 1

Short answer: No.

There's no way to use non-image data structures (like a stack), and no way to do dynamic recursion. It's not clear that Halide would really add much value here, because that algorithm doesn't seem to be tileable, parallelizable or vectorizable as it is written.

You may however be able to re-write the algorithm as one that makes iterative sweeps over the image flipping edges from weak to strong. It can be thought of as a cellular automata on three states (weak, strong, not an edge) run to completion, and we could vectorize/parallelize each pass. See test/correctness/gameoflife.cpp in the Halide repo for an example. I think this way of doing it would have a lousy computational complexity though. You'd be doing work at every pixel, not just at the live edge of pixels that are flipping.

You could also run it as a cellular automata that does in-place updates along some wavefront, e.g. do a sweep from top to bottom, bottom to top, left to right, and right to left. You can then vectorize along the wavefront. The schedule would be similar to an IIR (see https://github.com/halide/CVPR2015/blob/master/RecursiveFilter/IirBlur.cpp). That approach would handle a linear edge along any direction, but any fixed numbers of sweeps would miss a spiral moving from weak to strong.

But rather than contort your code in these ways, I would just use a different algorithm, or use define_extern.