What is wrong with my cython implementation of erosion operation of mathematical morphology

I have produced a naive implementation of "erosion". The performance is not relevant since I just trying to understand the algorithm. However, the output of my implementation does not match the one I get from scipy.ndimage. What is wrong with my implementation ?

Here is my implementation with a small test case:

import numpy as np
from PIL import Image

# a small image to play with a cross structuring element
imgmat = np.array([
    [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
    [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
    [0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0],
    [0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,0,0,0,0,0,0,0,0],
    [0,0,0,0,0,1,1,0,0,1,1,1,1,0,0,0,1,1,0,0,1,0,0,0,0,0,0,0,0,0],
    [0,0,0,0,0,1,1,0,0,1,1,1,1,0,0,0,1,1,1,0,1,0,0,0,0,0,0,0,0,0],
    [0,0,0,0,0,1,1,0,0,1,1,1,1,0,0,0,1,1,1,1,1,0,0,0,0,0,0,0,0,0],
    [0,0,0,0,0,1,1,0,0,1,1,1,1,0,0,0,1,1,1,1,1,0,0,0,0,0,0,0,0,0],
    [0,0,0,0,0,1,1,0,0,1,1,1,1,0,0,0,0,1,1,1,1,0,0,0,0,0,0,0,0,0],
    [0,0,0,0,0,1,1,0,0,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
    [0,0,0,0,0,1,1,0,0,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
    [0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,1,1,1,1,1,0,0,1,0,1,1,0,0,0,0],
    [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,0,1,0,1,1,0,0,0,0],
    [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,1,1,0,0,0,0],
    [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,1,1,0,0,0,0],
    [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,0,0,0,0],
])
imgmat2 = np.where(imgmat == 0, 0, 255).astype(np.uint8)
imarr = Image.fromarray(imgmat2).resize((100, 200))
imarr = np.array(imgrrr)
imarr = np.where(imarr == 0, 0, 1)

se_mat3 = np.array([
    [0,1,0],
    [1,1,1],
    [0,1,0]
])
se_mat31 = np.where(se_mat3 == 1, 0, 1)

The imarr is .

My implementation of erosion:

%%cython -a
import numpy as np
cimport numpy as cnp

cdef erosionC(cnp.ndarray[cnp.int_t, ndim=2] img, 
              cnp.ndarray[cnp.int_t, ndim=2] B, cnp.ndarray[cnp.int_t, ndim=2] X):
    """
    X: image coordinates
    struct_element_mat: black and white image, black region is considered as the shape 
                    of structuring element
                    
    This operation checks whether (B *includes* X) = $B subset X$ 
    as per defined in
    Serra (Jean), « Introduction to mathematical morphology », 
    Computer Vision, Graphics, and Image Processing,
    vol. 35, nᵒ 3 (septembre 1986). 
    URL : https://linkinghub.elsevier.com/retrieve/pii/0734189X86900022.. 
    doi: 10.1016/0734-189X(86)90002-2
    Consulted le 6 août 2020, p. 283‑305.
    """
    cdef cnp.ndarray[cnp.int_t, ndim=1] a, x, bx
    cdef cnp.ndarray[cnp.int_t, ndim=2] Bx, B_frame, Xcp, b
    cdef bint check
    a = B[0] # get an anchor point from the structuring element coordinates
    B_frame = B - a # express the se element coordinates in with respect to anchor point
    Xcp = X.copy() 
    b = img.copy()
    for x in X: # X contains the foreground coordinates in the image
        Bx = B_frame + x # translate relative coordinates with respect to foreground coordinates considering it as the anchor point
        check = True # this is erosion so if any of the se coordinates is not in foreground coordinates we consider it a miss
        for bx in Bx: # Bx contains all the translated coordinates of se
            if bx not in Xcp:
                check = False
        if check:
            b[x[0], x[1]] = 1 # if there is a hit
        else:
            b[x[0], x[1]] = 0 # if there is no hit
    return b
        
def erosion(img: np.ndarray, struct_el_mat: np.ndarray, foregroundValue = 0):
    B = np.argwhere(struct_el_mat == 0)
    X = np.argwhere(img == foregroundValue)
    nimg = erosionC(img, B, X)
    return np.where(nimg == 1, 255, 0)

The calling code for both is:

from scipy import ndimage as nd

err = nd.binary_erosion(imarr, se_mat3)

imerrCustom = erosion(imarr, se_mat31, foregroundValue=1)

err produces

imerrCustom produces

%%cython -a
import numpy as np
cimport numpy as cnp

cdef dilation_c(cnp.ndarray[cnp.uint8_t, ndim=2] X, 
                cnp.ndarray[cnp.uint8_t, ndim=2] SE):
    """
    X: boolean image
    SE: structuring element matrix
    origin: coordinate of the origin of the structuring element
                    
    This operation checks whether (B *hits* X) = $B cap X not = emptyset$ 
    as per defined in
    Serra (Jean), « Introduction to mathematical morphology », 
    Computer Vision, Graphics, and Image Processing,
    vol. 35, nᵒ 3 (septembre 1986). 
    URL : https://linkinghub.elsevier.com/retrieve/pii/0734189X86900022.. 
    doi: 10.1016/0734-189X(86)90002-2
    Consulted le 6 août 2020, p. 283‑305.
    
    The algorithm adapts DILDIRECT of
    Najman (Laurent) et Talbot (Hugues), 
    Mathematical morphology: from theory to applications, 
    2013. ISBN : 9781118600788, p. 329
    to the formula given in 
    Jähne (Bernd), 
    Digital image processing, 
    6th rev. and ext. ed, Berlin ; New York, 
    2005. TA1637 .J34 2005. 
    ISBN : 978-3-540-24035-8.

    """
    cdef cnp.ndarray[cnp.uint8_t, ndim=2] O
    cdef list elst
    cdef int r, c, X_rows, X_cols, SE_rows, SE_cols, se_r, se_c
    cdef cnp.ndarray[cnp.int_t, ndim=1] bp
    cdef list conds
    cdef bint check, b, p, cond
    O = np.zeros_like(X)
    X_rows, X_cols = X.shape[:2]
    SE_rows, SE_cols = SE.shape[:2]
    # a boolean convolution
    for r in range(0, X_rows-SE_rows):
        for c in range(0, X_cols - SE_cols):
            conds = []
            for se_r in range(SE_rows):
                for se_c in range(SE_cols):
                    b = <bint>SE[se_r, se_c]
                    p = <bint>X[se_r+r, se_c+c]
                    conds.append(b and p)
            O[r,c] = <cnp.uint8_t>any(conds)
    return O
    
        
def dilation_erosion(
    img: np.ndarray, 
    struct_el_mat: np.ndarray, 
    foregroundValue: int = 1,
    isErosion: bool = False):
    """
    img: image matrix
    struct_el: NxN mesh grid of the structuring element whose center is SE's origin
              structuring element is encoded as 1 
    foregroundValue: value to be considered as foreground in the image
    """
    B = struct_el_mat.astype(np.uint8)
    if isErosion:
        X = np.where(img == foregroundValue, 0, 1).astype(np.uint8)
    else:
        X = np.where(img == foregroundValue, 1, 0).astype(np.uint8)

    nimg = dilation_c(X, B)
    foreground, background = (255, 0) if foregroundValue == 1 else (0, 1)
    if isErosion:
        return np.where(nimg == 1, background, foreground).astype(np.uint8)
    else:
        return np.where(nimg == 1, foreground, background).astype(np.uint8)
    # return nimg