Конфликтующие ширины для объекта в моем 2d массиве? - PullRequest
Новая
       25

Конфликтующие ширины для объекта в моем 2d массиве?

0 голосов
/ 08 сентября 2018

У меня есть изображение (двумерный массив, форма 600x600) асимметричного шарика, который центрирован и повернут относительно осей массива. Я запустил 2D-гауссовскую подгонку, используя scipy.optimize.curve_fit, и получил приемлемые результаты для параметров подгонки. Ниже приведен мой код и мои результаты для этой части. 

from scipy import optimize
import numpy as np

x0, y0 = np.shape(image)[1]//2, np.shape(image)[0]//2

def gaussian_func(xy, x0, y0, sigma_x, sigma_y, amp, theta, offset):
    x, y = xy
    a = np.cos(theta)**2/(2*sigma_x**2) + np.sin(theta)**2/(2*sigma_y**2)
    b = -np.sin(2*theta)/(4*sigma_x**2) + np.sin(2*theta)/(4*sigma_y**2)
    c = np.sin(theta)**2/(2*sigma_x**2) + np.cos(theta)**2/(2*sigma_y**2)

    inner = a * (x-x0)**2
    inner += 2 * b * (x-x0) * (y-y0)
    inner += c * (y-y0)**2
    return (offset + amp * np.exp(-inner)).ravel()

def Sigma2FWHM(sigma):
    return 2 * np.sqrt(2*np.log(2)) * sigma

def generate(data_set):
    xvec = np.arange(0, np.shape(data_set)[1], 1)
    yvec = np.arange(0, np.shape(data_set)[0], 1)
    X, Y = np.meshgrid(xvec, yvec)
    return X, Y

# Fit subimage of image to 2D Gaussian (subimage helps with computation time)

# Guesses
theta_guess = np.deg2rad(96)
sigma_x, sigma_y = 5, 4
amp = 1 #the image was previously normalized so highest intensity point is 1
subimage = image[y0-14:y0+14, x0-14:x0+14] #the entire blob is contained within this region
offset = np.min(subimage)

guesses = [np.shape(subimage)[1]//2, np.shape(subimage}[0]//2, sigma_x, sigma_y, amp, theta_guess, offset]

# Do fit
xx, yy = generate(subimage)
pred_params, uncert_cov = optimize.curve_fit(gaussian_func, (xx.ravel(), yy.ravel()), subimage.ravel(), p0=guesses)
FWHM_x, FWHM_y = Sigma2FWHM(np.abs(pred_params[2])), Sigma2FWHM(np.abs(pred_params[2]))

# Add back origin to centroid and convert angle back to degrees
x_0, y_0 = pred_params[0]+(x0-14), pred_params[1]+(y0-14)
angle_deg = np.rad2deg(pred_params[5])

print('Predicted FWHM x, y (pixels):', FWHM_x, FWHM_y)

Выход: Predicted FWHM x, y (pixels): 13.503 11.351

У меня сложилось впечатление, что выполнение 2D гауссовой подгонки выберет угол, для которого ширина большой полуоси максимальна (по крайней мере, для меня это имеет смысл), и что выходные данные sigma_x и sigma_y будет действительно соответствовать стандартным отклонениям вдоль большой и малой осей (оси, повернутые на угол подгонки относительно стандартных осей x и y).

Однако из любопытства я вычислил полную ширину на половине максимума вручную на не повернутом изображении по стандартным осям x и y. Я ожидал чего-то меньшего, чем то, что дало мне гауссовское приближение, но результаты, которые я получил, были больше! Ниже приведен мой код для этой части. 

ylim, xlim = np.shape(image)
ypix, xpix = np.where(image==1)
x = np.take(image, ypix[0], axis=0)
y = np.take(image, xpix[0], axis=1)

#Print manually calculated FWHMs (integer estimates)
FWHM_x = len(np.where(x >= 0.5)[0])
FWHM_y = len(np.where(y >= 0.5)[0])
print('Estimated FWHM x, y without rotation:', FWHM_x, FWHM_y)

Выход: Estimated FWHM x, y without rotation: 15 13

Итак, мой вопрос, было ли моё предположение о соответствии Гаусса неверным? Я не могу понять, почему FWHM больше без предполагаемого вращения. Если да, то какую информацию на самом деле дают мне угол и стандартные отклонения, заданные подгонкой? Любое понимание приветствуется, так как это беспокоило меня некоторое время. Спасибо!

...