Приспособление по Гауссу не может вернуть правильные параметры для этого PSF - PullRequest
Новая
       64

Приспособление по Гауссу не может вернуть правильные параметры для этого PSF

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

Чтобы проверить мой 2D-гауссовский код подгонки для изображений ярких объектов, я запускаю его с помощью функции точечного разброса (PSF), которая была разработана и подобрана командой WISE. На своем веб-сайте они перечисляют параметры для центрального PSF для каждой полосы: FWHM вдоль большой и вспомогательной осей и угол положения (это угол от оси y). Вся информация доступна здесь: WISE PSF information , но ниже приведено изображение этих параметров для центрального PSF и соответствующее изображение PSF для полосы 3. WISE PSF parameters

WISE central PSF for W3

Итак, я загрузил соответствующий образ FITS для центрального PSF в полосе 3 (все изображения доступны для загрузки по вышеуказанной ссылке) и попытался запустить на нем свой код. Однако мой код не возвращает ожидаемые параметры, и параметры меняются в зависимости от того, подхожу ли я к подизображению (и в зависимости от его размера), или от целого изображения, что вызывает беспокойство.

Мне интересно, есть ли способ заставить мой код соответствия Гаусса восстановить наиболее точные параметры - или, возможно, другой метод подбора был бы более эффективным. Но меня больше всего беспокоит тот факт, что выходные параметры моей подгонки по Гауссу могут стать совершенно неверными. Ниже мой код. 

from scipy import optimize
import numpy as np
from astropy.io import fits

image = 'wise-w3-psf-wpro-09x09-05x05.fits' #WISE central PSF
stacked_image = fits.getdata(image)

# Center of image (the PSF is centered)
x0 = np.shape(stacked_image)[1]//2
y0 = np.shape(stacked_image)[0]//2

# Normalize image so peak = 1
def normalize(image):
    image *= 1/np.max(image)
    return image

stacked_image = normalize(stacked_image)

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 Sigma2width(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


# METHOD 1: Fit subimage of PSF to Gaussian

# Guesses
theta_guess = np.deg2rad(96) #I believe that the angle in the Gaussian corresponds to CCW from the x-axis (so use WISE position angle + 90 degrees)
sigma_x = 5
sigma_y = 4
amp = 1 #I know this is true since I normalized it
subimage = stacked_image[y0-50:y0+50, x0-50:x0+50]

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

xx, yy = generate(subimage)
pred_params, uncert_cov = optimize.curve_fit(gaussian_func, (xx.ravel(), yy.ravel()), subimage.ravel(), p0=guesses)
width_x, width_y = Sigma2width(np.abs(pred_params[2]))*0.275, Sigma2width(np.abs(pred_params[3]))*0.275 #multiply by pixel scale factor (available on website) to get FWHMs in arcseconds

x_0, y_0 = pred_params[0]+(x0-50), pred_params[1]+(y0-50) #add back origin
theta_deg = np.rad2deg(pred_params[5])
pred_params[5] = theta_deg
pred_params[0] = x_0
pred_params[1] = y_0

if theta_deg < 90:
    pos_angle = theta_deg+90
elif theta_deg >= 90:
    pos_angle = theta_deg-90

print('PREDICTED FWHM x, y in arcsecs:', width_x, width_y)
print('FIT PARAMS [x0, y0, sigma_x, sigma_y, amp, theta, offset]:', pred_params)
print('POSITION ANGLE:', pos_angle)

# Output: PREDICTED FWHM x, y in arcsecs: 6.4917 5.4978
# FIT PARAMS [x0, y0, sigma_x, sigma_y, amp, theta, offset]: [3.195e+02 3.189e+02 1.002e+01 8.489e+00 8.695e-01 8.655e+01 2.613e-02]
# POSITION ANGLE: 176.556


# METHOD 2: Fit whole image to Gaussian

# Guesses
theta_guess = np.deg2rad(96)
sigma_x = 5
sigma_y = 4
amp = 1
offset = np.median(stacked_image)
guesses = [x0, y0, sigma_x, sigma_y, amp, theta_guess, offset]

# Sigmas - manual estimation
ylim, xlim = np.shape(stacked_image)
x, y = np.arange(0, xlim, 1), np.arange(0, ylim, 1)
ypix, xpix = np.where(stacked_image==amp)

y_range = np.take(stacked_image, ypix[0], axis=0)
x_range = np.take(stacked_image, xpix[0], axis=1)

xx, yy = generate(stacked_image)
pred_params, uncert_cov = optimize.curve_fit(gaussian_func, (xx.ravel(), yy.ravel()), stacked_image.ravel(), p0=guesses)
width_x, width_y = Sigma2width(np.abs(pred_params[2]))*0.275, Sigma2width(np.abs(pred_params[3]))*0.275 #in arcsecs
theta = pred_params[5]

print('PREDICTED FWHM x, y in arcsecs:', width_x, width_y)
print('FIT PARAMS [x0, y0, sigma_x, sigma_y, amp, theta, offset]:', pred_params)

# Output: 
# PREDICTED FWHM x, y in arcsecs: 7.088 6.106
# FIT PARAMS [x0, y0, sigma_x, sigma_y, amp, theta, offset]: [3.195e+02 3.190e+02 1.095e+01 9.429e+00 8.378e-01 1.521e+00 7.998e-04]

if theta < 90:
    pos_angle = 90+np.rad2deg(theta)
elif theta >= 90:
    pos_angle = 90-np.rad2deg(theta)

print('POSITION ANGLE:', pos_angle)

# POSITION ANGLE: 177.147

Вы можете видеть в (округленных) выходных данных, что амплитуды, возвращаемые моими подгонками Гаусса, даже не равны 1, а другие параметры (FWHM и углы) также не совпадают с правильными параметрами, показанными в таблице.

Если я подгоняю подизображение, кажется, что амплитуда становится все ближе и ближе (но никогда не достигает) 1, чем меньше я делаю подизображение, но тогда FWHM может стать слишком маленьким по сравнению с реальными значениями. Почему я не получаю правильные результаты, и как я могу сделать подгонку максимально точной?

...