我对数学形态学侵蚀操作的 cython 实现有什么问题答案

【问题标题】：What is wrong with my cython implementation of erosion operation of mathematical morphology我对数学形态学侵蚀操作的 cython 实现有什么问题
【发布时间】：2020-08-06 12:19:47
【问题描述】：

我已经产生了“侵蚀”的幼稚实现。性能无关紧要，因为我只是想了解算法。但是，我的实现输出与我从scipy.ndimage 得到的输出不匹配。我的实现有什么问题？

这是我的一个小测试用例的实现：

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)

imarr 是。

我的侵蚀实现：


%%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)

两者的调用代码是：

from scipy import ndimage as nd

err = nd.binary_erosion(imarr, se_mat3)

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

err 产生

imerrCustom 产生

【问题讨论】：

好的，我正在删除另一个问题
输出是什么，它与预期的有什么不同？我正在尝试阅读您的代码，但无法理解它。 if bx not in Xcp我只是无法理解逻辑，对我来说它看起来根本不像扩张。我建议您使用图像和 SE 作为矩阵，而不是坐标列表。这将大大简化您的代码。至少从某种意义上说，我能够理解发生了什么......这是一个类似的问题，我能够大大简化过于复杂的代码：stackoverflow.com/q/63098284/7328782
按你说的添加了图片。 bx not in Xcp 检查翻译后的结构元素的坐标是否包含在图像前景的坐标中。我正在添加 cmets 以使其更容易消化
所以imarr不是二元矩阵，它有不同的层次。您应该在应用二元运算符之前对其进行阈值处理，以便您知道每个函数将什么视为前景和背景。 ndimage 可能认为任何非零值都是前景，但我不确定。强制它二进制，你肯定会知道。另外，您似乎认为零是前景？这与任何其他实现的做法相反。将您的腐蚀与 ndimage 的膨胀进行比较。
imarr 是一个二进制矩阵，只是论文没有区分0或1作为前景，所以我想我会让用户指定，因此`X = np.argwhere (img == foregroundValue)`，这不是 SE 的情况，我假设 0 作为 SE 坐标的指标

标签： numpy image-processing cython mathematical-morphology

【解决方案1】：

最后，我仍然不确定，但是在阅读了几篇论文之后，我认为我将X 解释为前景坐标是错误的。它可能应该是正在迭代的整个图像。正如我所说，我不确定这种解释是否也正确。但是我做了一个迭代图像的新实现，它给出了一个更合理的结果。我在这里分享它，希望它可以帮助某人：

%%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

【讨论】：