计算 numpy 数组之间的 MSE

Question

科学问题：

我有很多 3D 体积，其中都有一个圆柱体，圆柱体在 z 轴上“直立”。包含圆柱体的体积非常嘈杂，就像超级嘈杂一样，您无法像人类一样看到其中的圆柱体。如果我将这些体积的 1000 个平均放在一起，我可以看到圆柱体。每个卷都包含一个圆柱体的副本，但在少数情况下，圆柱体的方向可能不正确，所以我想要一种方法来解决这个问题。

我想出的解决方案：

我获取了平均体积并将其投影到 z 和 x 轴上（仅投影 numpy 数组），以便在一个方向上得到一个漂亮的圆圈，在另一个方向上得到一个矩形。然后我获取每个 3D 体积并沿 Z 轴投影每个单独的体积。 SNR 仍然很差，我看不到一个圆圈，但如果我平均 2D 切片，我可以在平均几百个之后开始看到一个圆圈，并且在前 1000 个平均后很容易看到。为了计算每个体积的分数，我计算了 3D 体积相对于其他三个阵列沿 z 投影的 MSE，第一个是沿 Z 投影的平均值，然后是沿 y 或 x 投影的平均值，最后是一个带有噪声在其中的正态分布。

目前我有以下内容，其中 RawParticle 是 3D 数据，Ave 是平均值：

def normalise(array):
    min = np.amin(array)
    max = np.amax(array)
    normarray = (array - min) / (max - min)

    return normarray

def Noise(mag):
    NoiseArray = np.random.normal(0, mag, size=(200,200,200))
    return NoiseArray

#3D volume (normally use a for loop to iterate through al particles but for this example just showing one)
RawParticleProjected = np.sum(RawParticle, 0)
RawParticleProjectedNorm = normalise(RawParticleProjected)
#Average
AveProjected = np.sum(Ave, 0)
AveProjectedNorm = normalise(AveProjected)
#Noise Array
NoiseArray = Noise(0.5)
NoiseNorm = normalise(NoiseArray)


#Mean squared error
MSE = (np.square(np.subtract(RawParticleProjectedNorm, AveProjectedNorm))).mean()

然后我将 Ave 与轴 1 相加重复此操作，然后再次将 Raw 粒子与 Noise 数组进行比较。

但是，当我比较应该都是圆形的投影时，我的输出给出了最高的 MSE，如下所示：

我对 MSE 的理解是，其他两个人群应该有高 MSE，而我同意的人群应该有低 MSE。也许我的数据对于这种类型的分析来说太嘈杂了？但如果这是真的，那么我真的不知道该怎么做。

如果有人可以浏览我的代码或启发我对 MSE 的理解，我将非常感激。

感谢您花时间浏览和阅读。

Answer 1

如果我正确理解了您的问题，您想知道不同样本与平均值的接近程度。 通过比较样本，您可以找到包含迷失方向圆柱体的异常值。 这非常符合L2 norm 的定义，所以MSE 应该在这里工作。

我会计算所有样本的平均 3D 图像，然后计算每个样本与该平均值的距离。然后我会比较这些值。

将样本与人工噪声图像进行比较的想法不错，但我不确定正态分布和归一化是否按您的计划进行。我可以是苹果和橘子。 而且我认为查看沿不同轴的投影不是一个好主意， 只需比较 3D 图像。

我用参数alpha 对二维圆进行了一些小测试，该参数指示图片中有多少噪声和多少圆。 （alpha=0 表示只有噪音，alpha=1 表示只有圆圈`）

import numpy as np
import matplotlib.pyplot as plt

grid_size = 20
radius = 5
mag = 1

def get_circle_stencil(radius):
    xx, yy = np.meshgrid(np.linspace(-grid_size/2+1/2, grid_size/2-1/2, grid_size),
                         np.linspace(-grid_size/2+1/2, grid_size/2-1/2, grid_size))
    dist = np.sqrt(xx**2 + yy**2)
    inner = dist < (radius - 1/2)
    return inner.astype(float)

def create_noise(mag, n_dim=2):
    # return np.random.normal(0, mag, size=(grid_size,)*n_dim)
    return np.random.uniform(0, mag, size=(grid_size,)*n_dim)

def create_noisy_sample(alpha, n_dim=2):
    return (np.random.uniform(0, 1-alpha, size=(grid_size,)*n_dim) + 
            alpha*get_circle_stencil(radius))


fig = plt.figure()
ax = fig.subplots(nrows=3, ncols=3)
np.unravel_index(3, shape=(3, 3))
alpha_list = np.arange(9) / 10
for i, alpha in enumerate(alpha_list):
    r, c = np.unravel_index(i, shape=(3, 3))
    ax[r][c].imshow(*norm(create_noisy_sample(alpha=alpha)), cmap='Greys')
    ax[r][c].set_title(f"alpha={alpha}")
    ax[r][c].xaxis.set_ticklabels([])
    ax[r][c].yaxis.set_ticklabels([])

我尝试了一些指标（mse、cosine similarity 和 binary cross entropy 并查看它们在不同 alpha 值下的表现。

def normalize(*args):
    return [a / np.linalg.norm(a) for a in args]

def cosim(a, b):
    return np.sum(a * b)

def mse(a, b):
    return np.sqrt(np.sum((a-b)**2))

def bce(a, b):
    # binary cross entropy implemented from tensorflow / keras
    eps = 1e-7
    res = a * np.log(b + eps)
    res += (1 - a) * np.log(1 - b + eps)
    return np.mean(-res)

我比较了NoiseA-NoiseB、Circle-Circle、Circle-Noise、Noise-Sample、Circle-Sample

alpha = 0.1
noise = create_noise(mag=1, grid_size=grid_size)
noise_b = create_noise(mag=1, grid_size=grid_size)
circle_reference = get_circle_stencil(radius=radius, grid_size=grid_size)
sample = create_noise(mag=1, grid_size=grid_size) + alpha * circle_reference

print('NoiseA-NoiseB:', mse(*norm(noise, noise_b)))    # 0.718
print('Circle-Circle:', mse(*norm(circle, circle)))    # 0.000
print('Circle-Noise:', mse(*norm(circle, noise)))      # 1.168
print('Noise-Sample:', mse(*norm(noise, sample)))      # 0.697
print('Circle-Sample:', mse(*norm(circle, sample)))    # 1.100

print('NoiseA-NoiseB:', cosim(*norm(noise, noise_b)))  # 0.741
print('Circle-Circle:', cosim(*norm(circle, circle)))  # 1.000
print('Circle-Noise:', cosim(*norm(circle, noise)))    # 0.317
print('Noise-Sample:', cosim(*norm(noise, sample)))    # 0.757
print('Circle-Sample:', cosim(*norm(circle, sample)))  # 0.393

print('NoiseA-NoiseB:', bce(*norm(noise, noise_b)))    # 0.194
print('Circle-Circle:', bce(*norm(circle, circle)))    # 0.057
print('Circle-Noise:', bce(*norm(circle, noise)))      # 0.111
print('Noise-Circle:', bce(*norm(noise, circle)))      # 0.636
print('Noise-Sample:', bce(*norm(noise, sample)))      # 0.192
print('Circle-Sample:', bce(*norm(circle, sample)))    # 0.104

n = 1000
ns = np.zeros(n)
cs = np.zeros(n)
for i, alpha in enumerate(np.linspace(0, 1, n)):
    sample = create_noisy_sample(alpha=alpha)
    ns[i] = mse(*norm(noise, sample))
    cs[i] = mse(*norm(circle, sample))

fig, ax = plt.subplots()
ax.plot(np.linspace(0, 1, n), ns, c='b', label='noise-sample')
ax.plot(np.linspace(0, 1, n), cs, c='r', label='circle-sample')
ax.set_xlabel('alpha')
ax.set_ylabel('mse')
ax.legend()

对于您的问题，我只想查看比较 circle-sample（红色）。 不同的样本会表现得好像它们具有不同的 alpha 值，您可以对它们进行相应的分组。而且您应该能够检测到异常值，因为它们应该具有更高的mse。

你说你必须结合 100-1000 张图片才能看到圆柱体，这表明你的问题中的 alpha 值非常小，但 平均 mse 应该工作。