python - 计算 numpy 数组之间的 MSE
问题描述
科学问题:
我有很多 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 相加重复此操作,然后再次将原始粒子与噪声数组进行比较。
但是,当我比较应该都是圆形的投影时,我的输出给出了最高的 MSE,如下所示:
我对 MSE 的理解是,其他两个人群应该有高 MSE,而我同意的人群应该有低 MSE。也许我的数据对于这种类型的分析来说太嘈杂了?但如果这是真的,那么我真的不知道该怎么做我正在做的事情。
如果有人可以浏览我的代码或启发我对 MSE 的理解,我将非常感激。
感谢您花时间浏览和阅读。
解决方案
如果我正确理解了您的问题,您想弄清楚您的不同样本与平均值的接近程度。通过比较样本,您可以找到包含迷失方向圆柱体的异常值。这非常符合 的定义L2 norm
,所以MSE
应该在这里工作。
我将计算所有样本的平均 3D 图像,然后计算每个样本与该平均值的距离。然后我会比较这些值。
将样本与人工噪声图像进行比较的想法还不错,但我不确定正态分布和您的归一化是否按您的计划进行。我可以是苹果和橘子。而且我认为沿不同轴查看投影不是一个好主意,只需比较 3D 图像即可。
我用 2D 圆圈做了一些小测试,参数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
应该可以。
推荐阅读
- javascript - Javascript:for循环检索和打印一定数量的API数据
- java - 如何使用 RESTAPI 在 Java 中检索特定的 JIRA 字段
- ssl - RabbitMQ - Tls 连接 - Golang
- c# - 如何在 C# 中使用搜索逻辑?
- php - 如何根据数据库结果将多个 PHP 复选框设置为选中?
- visual-studio - 在 Visual Studio 2017 中,如何/在哪里配置解决方案资源管理器的解决方案和文件夹视图?
- java - Spark 无法按包含文件夹路径的字符串进行过滤
- cucumber - Gherkin 场景、可重用步骤方法或特定方法
- android - 搜索 Android RecyclerView
- python-3.x - 删除文件夹中除最新的特定文件之外的所有文件