Расчет MSE между numpy массивами

Question

The Scientifi c Вопрос:

У меня много трехмерных томов с цилиндром в них, ориентированным с цилиндром "вертикально" на оси z. Объемы, содержащие цилиндр, невероятно шумные, как супершумный, вы не можете видеть цилиндр в них как человек. Если я усредню вместе тысячи этих объемов, я смогу увидеть цилиндр. Каждый том содержит копию цилиндра, но в некоторых случаях цилиндр может быть не правильно ориентирован, поэтому я хочу выяснить это.

Решение, которое я придумал:

Я взял усредненный объем и спроецировал его на оси z и x (просто проецируя массив numpy), чтобы получить хороший круг в одном направлении и прямоугольник в другом. Затем я беру каждый объем 3D и проецирую каждый вниз по оси Z. SNR все еще настолько плох, что я не могу видеть круг, но если я усредню 2D-срезы, я могу начать видеть круг после усреднения нескольких сотен, и это легко увидеть после усреднения первых 1000. Чтобы вычислить оценку того, как каждый том, который я рассчитывал, вычисляя MSE трехмерных объемов, спроецированных вниз по z, на три других массива, сначала будет среднее значение, спроецированное вниз по Z, затем среднее значение, спроецированное вниз по y или x, и, наконец, массив с нормальное распределение шума в нем.

В настоящее время у меня есть следующее, где RawParticle - это трехмерные данные, а 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, я был бы очень признателен.

Спасибо, что нашли время посмотреть и прочитать.

Answer 1

Если я правильно понял ваш вопрос, вы хотите выяснить, насколько близки ваши средние образцы к среднему. И сравнивая образцы, вы хотите найти выбросы, которые содержат дезориентированный цилиндр. Это очень хорошо соответствует определению L2 norm, поэтому здесь должно работать * 1002. *

Я бы рассчитал среднее 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 и посмотрел, как они ведут себя при разных значениях альфа.

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 (красный). Различные образцы будут вести себя так, как если бы они имели различные альфа-значения, и вы можете сгруппировать их соответствующим образом. И вы должны быть в состоянии обнаружить выбросы, потому что они должны иметь более высокое mse.

Вы сказали, что вам нужно объединить 100-1000 изображений, чтобы увидеть цилиндр, который Признак, что у вас действительно небольшое альфа-значение в вашей проблеме, но в среднее mse должно работать.