Stacking star PSF from image; subpixel center alignment
I have an array (1727,1853) of size (image) in which I have identified stars to simulate a point expansion function. Each array index corresponds to an image coordinate, but the centroid of each star is specified by a subpixel coordinate. I have to do the following
-
Make a 2D slice of each star. I did it using numpy array slicing. However, it is slice by index and I have the coordinates of the subpixel center of gravity, so whatever kind of slice I do will place the star off-center.
-
After I make a 2D slice of each star, I have to stack these arrays on top of each other to create a point-expanding function model. It's easy if the subpixel centers of each star are aligned.
My question is, is this the most efficient (and correct) way of aligning these subpixel coordinates and stacking each 2D slice together?
Hope this was clear. Any help is appreciated. Below is a 2D slice of one of the stars (not very good), however, it is off-center because the numerical fragments by index and stars have subpixel coordinates for their centroids.
source to share
You can express the central coordinates of the pixels in each "slice" relative to the star's center of gravity, and then compute a weighted two-dimensional histogram.
First, some sample data:
import numpy as np
from matplotlib import pyplot as plt
# pixel coordinates (integer)
x, y = np.mgrid[:100, :100]
# centroids (float)
cx, cy = np.random.rand(2, 9) * 100
# a Gaussian kernel to represent the PSF
def gausskern(x, y, cx, cy, sigma):
return np.exp(-((x - cx) ** 2 + (y - cy) ** 2) / (2 * sigma ** 2))
# (nstars, ny, nx)
stars = gausskern(x[None, ...], y[None, ...],
cx[:, None, None], cy[:, None, None], 10)
# add some noise for extra realism
stars += np.random.randn(*stars.shape) * 0.5
fig, ax = plt.subplots(3, 3, figsize=(5, 5))
for ii in xrange(9):
ax.flat[ii].imshow(stars[ii], cmap=plt.cm.hot)
ax.flat[ii].set_axis_off()
fig.tight_layout()
Weighted 2D histogram:
# (nstars, ny, nx) pixel coordinates relative to each centroid
dx = cx[:, None, None] - x[None, ...]
dy = cy[:, None, None] - y[None, ...]
# 2D weighted histogram
bins = np.linspace(-50, 50, 100)
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.pcolormesh(xe, ye, h, cmap=plt.cm.hot)
ax.axhline(0, ls='--', lw=2, c='w')
ax.axvline(0, ls='--', lw=2, c='w')
ax.margins(x=0, y=0)
source to share