python - 计算numpy数组之间的MSE

Question

科学问题:

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

我想出的解决方案:

我取了平均体积并将其沿 z 轴和 x 轴投影(仅投影 numpy 数组)，以便在一个方向上得到一个漂亮的圆圈，在另一个方向上得到一个矩形。然后我将每个 3D 体积都投影到 Z 轴上。 SNR 仍然很糟糕，我看不到一个圆，但是如果我对 2D 切片求平均，我可以在平均几百个之后开始看到一个圆，并且在前 1000 个平均后很容易看到。为了计算每个体积我如何计算 3D 体积的 MSE 相对于其他三个阵列向下投影 z 的分数，第一个是向下投影 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 重复此操作，然后再次将原始粒子与噪声阵列进行比较。

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



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

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

感谢您花时间查看和阅读。

Answer 1

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

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

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

我用参数 alpha 用 2D 圆做了一些小测试这表明图片中有多少噪音和多少圆圈。
( 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应该管用。