Astropy.model 2DГаусский выпуск

Question

Я пытаюсь сделать 2D гауссовское изображение, чтобы найти самый яркий объект на изображении.Я пытаюсь следовать примеру 1D здесь http://docs.astropy.org/en/v0.3.2/modeling/index.html. Мой код выглядит как

import numpy as np
import astropy.io.fits as fits
import os
from astropy.modeling import models, fitting
import matplotlib as plt

data = fits.getdata('AzTECC100.fits')
med = np.median(data) 
data = data - med
data = data[0,0,:,:]

w = models.Gaussian2D(data, np.mean(data),np.mean(data), np.std(data),np.std(data),)


fit_w = fitting.LevMarLSQFitter()

max_index = np.where(data >= np.max(data))
x0 = max_index[1] #Middle of X axis
y0 = max_index[0]
sigma = np.std(data)
amp = np.max(data)

x = np.arange(0, data.shape[1], 1)
y = np.arange(0, data.shape[0], 1)


A = 1 / (2*sigma**2)
eq =  amp*np.exp(-A*((x-x0)**2 + (y-y0)**2))
g = fit_w(w, x, y,eq)
plt.figure(figsize=(8,5))
plt.plot(x, y,eq, 'ko')
plt.plot(x, g, label='Gaussian')
circle = Circle((x0, y0), 4, facecolor ='none', edgecolor = 'red', linewidth = 1)
 ax.scatter(x0, y0,s = 10, c = 'red', marker = 'x')

Хотя я получаю ошибку

TypeError                                 Traceback (most recent call 
last)
<ipython-input-33-45f671ba149d> in <module>()
----> 1 g = fit_w(w, x, y,eq)

/Users/samclyne/anaconda/lib/python2.7/site- 
packages/astropy/modeling/fitting.pyc in __call__(self, model, x, y, z, weights, maxiter, acc, epsilon, estimate_jacobian)
    557             self.objective_function, init_values, args=farg, 
Dfun=dfunc,
    558             col_deriv=model_copy.col_fit_deriv, maxfev=maxiter, epsfcn=epsilon,
--> 559             xtol=acc, full_output=True)
    560         _fitter_to_model_params(model_copy, fitparams)
    561         self.fit_info.update(dinfo)

/Users/samclyne/anaconda/lib/python2.7/site- 
packages/scipy/optimize/minpack.pyc in leastsq(func, x0, args, Dfun, 
full_output, col_deriv, ftol, xtol, gtol, maxfev, epsfcn, factor, diag)
    378     m = shape[0]
    379     if n > m:
--> 380         raise TypeError('Improper input: N=%s must not exceed 
M=%s' % (n, m))
    381     if epsfcn is None:
    382         epsfcn = finfo(dtype).eps

TypeError: Improper input: N=65541 must not exceed M=65536

Я понимаю, что числопараметров не может превышать количество точек данных, которые есть, но я не уверен, что вызвало это или где я могу это исправить.

Извините, если я не предоставил достаточно информации в моем вопросе.Я не большой программист, и я новичок в этом сайте

Answer 1

В вашем коде есть несколько ошибок, но одна из причин этого конкретного сбоя заключается в том, что вы не создали модель правильно: первый аргумент w = models.Gaussian2D() должен быть амплитудой, а не data!Данные не имеют ничего общего с «моделью», которая представляет собой теоретическую зависимость / функцию, которую вы планируете подогнать к данным .

Использование данных, сгенерированных так же, как в https://stackoverflow.com/a/50560755/8033585,,

y, x = np.indices((51,51))
data = 3 * np.exp(-0.7 * ((x - 24)**2 + (y - 26)**2))

, а также закомментирование строк, которые не работают, вот модификация вашего кода, которая должна работать (примечание: я увеличил радиус круга в 30 раз, и график интенсивности находится на шкале логарифма):

import astropy.io.fits as fits
import os
from astropy.modeling import models, fitting
import matplotlib.pyplot as plt

#data = fits.getdata('AzTECC100.fits')
med = np.median(data) 
# data = data - med
# data = data[0,0,:,:] # NOT SURE THIS IS NEEDED!

fit_w = fitting.LevMarLSQFitter()

y0, x0 = np.unravel_index(np.argmax(data), data.shape)
sigma = np.std(data)
amp = np.max(data)

w = models.Gaussian2D(amp, x0, y0, sigma, sigma)

yi, xi = np.indices(data.shape)

g = fit_w(w, xi, yi, data)
print(w)
model_data = g(xi, yi)

fig, ax = plt.subplots()
eps = np.min(model_data[model_data > 0]) / 10.0
# use logarithmic scale for sharp Gaussians
ax.imshow(np.log(eps + model_data), label='Gaussian') 
circle = Circle((g.x_mean.value, g.y_mean.value),
                30 * g.x_stddev.value, facecolor ='none',
                edgecolor = 'red', linewidth = 1)

ax.add_patch(circle)
plt.show()

Параметры печатной модели:

Model: Gaussian2D
Inputs: ('x', 'y')
Outputs: ('z',)
Model set size: 1
Parameters:
    amplitude x_mean y_mean     x_stddev        y_stddev    theta
    --------- ------ ------ --------------- --------------- -----
          3.0   24.0   26.0 0.0881184627394 0.0881184627394   0.0

---

ПРИМЕЧАНИЯ:

Вы найдете множество примеров использования astropy.modeling иВажное обсуждение влияния «масштаба» данных изображения (порядка величины) здесь: https://github.com/astropy/astropy/issues/6269 В соответствии с этой проблемой может быть полезно «нормализовать» ваше изображение (изменить масштаб изображения, чтобы привести его в«разумный» диапазон).Также перейдите по ссылкам на записные книжки, которые могут предоставить дополнительные примеры, такие как https://github.com/ibusko/experiments_specviz/blob/master/experiment_normalization.ipynb К сожалению, эти примеры предназначены для данных 1D, но их можно легко расширить до данных 2D.