pub fn suppress_non_maximum<I, C>(image: &I, radius: u32) -> Image<Luma<C>>
where
    I: GenericImage<Pixel = Luma<C>>,
    C: Primitive + Ord,
{
    let (width, height) = image.dimensions();
    let mut out: Image<Luma<C>> = Image::new(width, height);
    if width == 0 || height == 0 {
        return out;
    }

    // We divide the image into a grid of blocks of size r * r. We find the maximum
    // value in each block, and then test whether this is in fact the maximum value
    // in the (2r + 1) * (2r + 1) block centered on it. Any pixel that's not maximal
    // within its r * r grid cell can't be a local maximum so we need only perform
    // the (2r + 1) * (2r + 1) search once per r * r grid cell (as opposed to once
    // per pixel in the naive implementation of this algorithm).

    for y in (0..height).step_by(radius as usize + 1) {
        for x in (0..width).step_by(radius as usize + 1) {
            let mut best_x = x;
            let mut best_y = y;
            let mut mi = image.get_pixel(x, y)[0];

            // These mins are necessary for when radius > min(width, height)
            for cy in y..cmp::min(height, y + radius + 1) {
                for cx in x..cmp::min(width, x + radius + 1) {
                    let ci = unsafe { image.unsafe_get_pixel(cx, cy)[0] };
                    if ci < mi {
                        continue;
                    }
                    if ci > mi || (cx, cy) < (best_x, best_y) {
                        best_x = cx;
                        best_y = cy;
                        mi = ci;
                    }
                }
            }

            let x0 = best_x.saturating_sub(radius);
            let x1 = x;
            let x2 = cmp::min(width, x + radius + 1);
            let x3 = cmp::min(width, best_x + radius + 1);

            let y0 = best_y.saturating_sub(radius);
            let y1 = y;
            let y2 = cmp::min(height, y + radius + 1);
            let y3 = cmp::min(height, best_y + radius + 1);

            // Above initial r * r block
            let mut failed = contains_greater_value(image, best_x, best_y, mi, y0, y1, x0, x3);
            // Left of initial r * r block
            failed |= contains_greater_value(image, best_x, best_y, mi, y1, y2, x0, x1);
            // Right of initial r * r block
            failed |= contains_greater_value(image, best_x, best_y, mi, y1, y2, x2, x3);
            // Below initial r * r block
            failed |= contains_greater_value(image, best_x, best_y, mi, y2, y3, x0, x3);

            if !failed {
                unsafe { out.unsafe_put_pixel(best_x, best_y, Luma([mi])) };
            }
        }
    }

    out
}
Returned image has zeroes for all inputs pixels which do not have the greatest
intensity in the (2 * radius + 1) square block centred on them.
Ties are resolved lexicographically.
src/suppress.rs