Astropy.model 2DGaussian问题 [英] Astropy.model 2DGaussian issue
问题描述
我正在尝试使2D高斯拟合图像以找到图像中最亮的对象.我正在尝试遵循此处的一维示例 http://docs.astropy .org/en/v0.3.2/modeling/index.html .我的代码看起来像
I'm trying to a 2D Gaussian to a fits image in order to find the brightest object in the image. I am trying to follow the 1D example here http://docs.astropy.org/en/v0.3.2/modeling/index.html. My code looks like
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
我知道参数的数量不能超过数据点的数量，但是我不确定是什么原因或在哪里可以修复它.
I understand that the number of parameters cannot exceed the amount of data points there are but I'm uncertain whats causing it or where I can fix it.
很抱歉，如果我在提问中没有提供足够的信息.我不是一个优秀的编码人员，而且是这个网站的新手
I'm sorry if I don't give enough info in my question. I'm not a great coder and I'm new to this website
推荐答案
您的代码中有几个错误，但造成此特定崩溃的一个错误是您没有正确创建模型:w = models.Gaussian2D()的第一个参数应为是幅度而不是data！数据与模型"无关，模型"是您计划适合数据的理论依赖/功能.
There are several mistakes in your code but the one responsible for this particular crash is that you do not create the model correctly: the first argument to w = models.Gaussian2D() should be the amplitude and not data! Data have nothing to do with a "model" which is a theoretical dependence/function that you plan to fit to data.
使用类似 https://stackoverflow.com/a/50560755/8033585 中生成的数据，具体来说，
Using data similarly generated as in https://stackoverflow.com/a/50560755/8033585, specifically,
y, x = np.indices((51,51))
data = 3 * np.exp(-0.7 * ((x - 24)**2 + (y - 26)**2))
并注释掉不起作用的行，这是对代码的修改，该行应该起作用(注意:我将圆的半径增加了30倍，并且强度图在对数刻度上):
and also commenting out lines that do not work, here is a modification of your code that should work (note: I increased the radius of the circle by 30x and intensity plot is on the log scale):
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/ibusko/experiments_specviz/blob/master/experiment_normalization.ipynb 不幸的是，这些示例适用于一维数据，但可以轻松扩展为二维数据.
NOTES:
You will find a great deal of examples of using astropy.modeling and an important discussion about the effect of the "scale" of the image data (order of magnitude) here: https://github.com/astropy/astropy/issues/6269 According to that issue, it may be helpful to "normalize" your image (rescale intensity to bring it into a "reasonable" range). Also follow links to notebooks that may provide further examples, such as https://github.com/ibusko/experiments_specviz/blob/master/experiment_normalization.ipynb Unfortunately these examples are for 1D data but can be easily expanded to 2D data.
这篇关于Astropy.model 2DGaussian问题的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
