在二进制图像的 Numpy 矩阵上通过矩阵算术计算均方误差答案

【问题标题】：Calculating Mean Squared Error through Matrix Arithmetic on Numpy Matrices of Binary Images在二进制图像的 Numpy 矩阵上通过矩阵算术计算均方误差
【发布时间】：2019-04-04 02:41:14
【问题描述】：

我有 2 张二值图像，一张是地面实况，一张是我制作的图像分割。

我正在尝试计算均方距离...

Let G = {g1, g2, . . . , gN} be the points in the ground truth boundary.
Let B = {b1, b2, . . . , bM} be the points in the segmented boundary.
Define d(p, p0) be a measure of distance between points p and p0 (e.g. Euclidean, city block, etc.)

使用以下算法在两个图像之间。

def MSD(A,G):
    '''
    Takes a thresholded binary image, and a ground truth img(binary), and computes the mean squared absolute difference
    :param A: The thresholded binary image
    :param G: The ground truth img
    :return:
    '''
    sim = np.bitwise_xor(A,G)
    sum = 0
    for i in range(0,sim.shape[0]):
        for j in range(0,sim.shape[1]):
            if (sim[i,j] == True):
                min = 9999999
                for k in range(0,sim.shape[0]):
                    for l in range(0,sim.shape[1]):
                        if (sim[k, l] == True):
                            e = abs(i-k) + abs(j-l)
                            if e < min:
                                min = e
                                mink = k
                                minl = l
                sum += min
    return sum/(sim.shape[0]*sim.shape[1])

这个算法太慢了，而且永远不会完成。

example 和 this example（答案 3）可能会展示如何使用矩阵算术获得均方误差的方法，但我不明白这些示例有何意义或它们为何起作用。

【问题讨论】：

您的代码与方程式不匹配。通过计算对称差 (xor)，您的代码是对称的。如果我交换 A 和 G 我会得到相同的结果。这不适用于方程，这不是对称距离测量。 B (A) 中不在 G 中的像素会增加平均距离。 G 中的像素不在 B 中。

标签： python arrays numpy image-processing image-segmentation

【解决方案1】：

因此，如果我正确理解了您的公式和代码，您将拥有一个（二进制）图像B 和一个（基本事实）图像G。 “点”由任一图像具有True（或至少非零）值的像素位置定义。从您的bitwise_xor 我推断出两个图像具有相同的形状(M,N)。

所以d^2(b,g) 的数量最多是一个(M*N, M*N) 大小的数组，将B 的每个像素与G 的每个像素相关联。更好的是：如果B 中存在m 非零值，G 中存在n 非零值，我们只需要一个形状(m,n)。除非您的图像很大，否则我们可以跟踪这么大的数量。这将消耗内存，但我们将通过矢量化赢得大量 CPU 时间。因此，对于每个m，我们只需要针对每个n 可能值找到这个距离的最小值。然后总结每个最小值。请注意，下面的解决方案使用了极端矢量化，如果图像很大，它很容易占用您的内存。

假设曼哈顿距离（您的代码中似乎缺少d^2 中的正方形）：

import numpy as np

# generate dummy data
M,N = 100,100
B = np.random.rand(M,N) > 0.5
G = np.random.rand(M,N) > 0.5

def MSD(B, G):
    # get indices of nonzero pixels
    nnz_B = B.nonzero() # (x_inds, y_inds) tuple, x_inds and y_inds are shape (m,)
    nnz_G = G.nonzero() # (x_inds', y_inds') each with shape (n,)

    # np.array(nnz_B) has shape (2,m)
    # compute squared Manhattan distance
    dist2 = abs(np.array(nnz_B)[...,None] - np.array(nnz_G)[:,None,:]).sum(axis=0)**2 # shape (m,n)
    # alternatively: Euclidean for comparison:
    #dist2 = ((np.array(nnz_B)[...,None] - np.array(nnz_G)[:,None,:])**2).sum(axis=0)

    mindist2 = dist2.min(axis=-1) # shape (m,) of minimum square distances

    return mindist2.mean() # sum divided by m, i.e. the MSD itself

print(MSD(B, G))

如果上面使用了太多内存，我们可以在nnz_B 的元素上引入一个循环，并且只对nnz_G 的元素进行向量化。这将占用更多的 CPU 资源和更少的内存。这种权衡是矢量化的典型特征。

【讨论】：

@CrisLuengo 它可能看起来是对称的，但事实并非如此。这正是公式所说的。 dist2 成对距离是“对称的”（就像 d^2(b,g) 一样），但最小化发生在它的第二维和第一维的总和。如果这还不够令人信服：MSD(B, G) != MSD(G, B) 如果我将度量标准更改为欧几里得，我会得到与你完全相同的结果；）但是你让我去那里一秒钟！
是的，我让自己去那里一分钟。 :) 抱歉，我还在习惯 Numpy 的索引方式。您确实沿二维距离矩阵的一个维度取最小值，然后沿另一个维度求和。完全有道理。

【解决方案2】：

计算此距离的一种有效方法是使用距离变换。 SciPy 在 ndimage 包中有一个实现：scipy.ndimage.morphology.distance_transform_edt。

这个想法是为地面真实图像G的背景计算距离变换。这将生成一个新图像D，对于G 中的每个非零像素，该图像为0，对于G 中的每个零像素，将有到最近的非零像素的距离。

接下来，对于B（或您发布的代码中的A）中的每个非零像素，您查看D 中的相应像素。这是该像素到G 的距离。因此，只需对 D 中 B 不为零的所有值进行平均即可获得结果。

import numpy as np
import scipy.ndimage as nd
import matplotlib.pyplot as pp

# Create some test data
img = pp.imread('erika.tif')         # a random image
G = img > 120                        # the ground truth
img = img + np.random.normal(0, 20, img.shape)
B = img > 120                        # the other image

D = nd.morphology.distance_transform_edt(~G)
msd = np.mean(D[B]**2)

【讨论】：