偏心加权 numpy histogram2d？

Question

我正在尝试生成 model PSF from a set of observed stars. I'm following the great example provided by ali_m in （下方的 MCVE）

我使用的 5 颗星是这样的：

中心（峰值强度）位于 bins [9, 9]。他们通过 numpy 的 hitsogram2d 组合的结果是这样的：

在 bin [8, 8] 处显示峰值密度。要以 [9, 9] 为中心，我必须获得质心（见下文）：

cx, cy = np.array([1.] * len(stars)), np.array([1.] * len(stars))

相反。这是为什么？

import numpy as np
from matplotlib import pyplot as plt

stars = # Uploaded here: http://pastebin.com/tjLqM9gQ

fig, ax = plt.subplots(2, 3, figsize=(5, 5))
for i in range(5):
    ax.flat[i].imshow(
        stars[i], cmap=plt.cm.viridis, interpolation='nearest',
        origin='lower', vmin=0.)
    ax.flat[i].axhline(9., ls='--', lw=2, c='w')
    ax.flat[i].axvline(9., ls='--', lw=2, c='w')
fig.tight_layout()

# (nstars, ny, nx) pixel coordinates relative to each centroid
# pixel coordinates (integer)
x, y = np.mgrid[:20, :20]
# centroids (float)
cx, cy = np.array([0.] * len(stars)), np.array([0.] * len(stars))
dx = cx[:, None, None] + x[None, ...]
dy = cy[:, None, None] + y[None, ...]

# 2D weighted histogram
bins = np.linspace(0., 20., 20)
h, xe, ye = np.histogram2d(dx.ravel(), dy.ravel(), bins=bins,
                           weights=stars.ravel())
fig, ax = plt.subplots(1, 1, subplot_kw={'aspect': 'equal'})
ax.hold(True)
ax.imshow(h, cmap=plt.cm.viridis, interpolation='nearest',
          origin='lower', vmin=0.)
ax.axhline(8., ls='--', lw=2, c='w')
ax.axvline(8., ls='--', lw=2, c='w')

plt.show()

Answer 1

之所以直方图没有以单星强度分布为中心的点 (9,9) 为中心，是因为生成它的代码围绕直方图的 bin 移动。

正如我在评论中所建议的那样，让事情保持简单。例如。我们不需要绘图来查看问题。另外，我不明白那些 dx dy 是什么，所以让我们避免它们。

然后我们可以通过

计算直方图

import numpy as np

stars = # Uploaded here: http://pastebin.com/tjLqM9gQ

# The argmax of a single star results in (9,9)
single_star_argmax  = np.unravel_index(np.argmax(stars[0]), stars[0].shape)

# Create a meshgrid of coordinates (0,1,...,19) times (0,1,...,19)
y,x = np.mgrid[:len(stars[0,:,0]), :len(stars[0,0,:])]
# duplicating the grids
xcoord, ycoord = np.array([x]*len(stars)),  np.array([y]*len(stars)) 
# compute histogram with coordinates as x,y
# and [20,20] bins    
h, xe, ye = np.histogram2d(xcoord.ravel(), ycoord.ravel(), 
                           bins=[len(stars[0,0,:]), len(stars[0,:,0])],
                           weights=stars.ravel())

# The argmax of the combined stars results in (9,9)
combined_star_argmax =  np.unravel_index(np.argmax(h), h.shape)

print single_star_argmax
print combined_star_argmax
print single_star_argmax == combined_star_argmax
# prints:
# (9, 9)
# (9, 9)
# True

原始代码中唯一的问题确实是在 0 和 20 之间创建 20 个点的行 bins = np.linspace(0., 20., 20)，
0. 1.05263158 2.10526316 ... 18.94736842 20.
这会将 bin 大小缩放到 ~1.05 并让您的 argmax 已经发生 "earlier" 然后预期。
您真正想要的是 0 到 19 之间的 20 个点，np.linspace(0,19,20) 或 np.arange(0,20)
为了避免这样的错误，可以简单地给出原始数组的长度作为参数，bins=20。

偏心加权 numpy histogram2d？

Off-centered weighted numpy histogram2d?

python

numpy

matplotlib

astronomy