In this notebook we show how to use the main methods of astroalign
.
Full documentation available at https://astroalign.readthedocs.io/.
This notebook was created for astroalign version 2.1
import astroalign as aa
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
np.random.seed(seed=12)
aa.__version__
We will create two images of the same portion of the sky. Our source image will end up being (400, 400) pixels, have 10 stars randomly distributed. It will have a Gaussian PSF with a sigma of 2 pixels and Gaussian background noise.
The target image will be (200, 200) pixels, will be rotated 30 degrees clockwise with respect to source image. It will have a Gaussian PSF with a sigma of 1.5 pixels and Gaussian background noise.
Notice that source image is made so to have better resolution (and seeing) than the target.
h, w = img_shape = (200, 200)
n_stars = 10
pos_x = np.random.randint(10, w - 10, n_stars)
pos_y = np.random.randint(10, h - 10, n_stars)
fluxes = 200.0 + np.random.rand(n_stars) * 300.0
img_posxy = np.array([[x, y] for x, y in zip(pos_x, pos_y)], dtype="float64")
img = np.zeros(img_shape)
for x, y, f in zip(pos_x, pos_y, fluxes):
img[y, x] = f
# Let's rotate and expand the original image
from scipy.ndimage import rotate, zoom
img_rotated = rotate(img, angle=30.0, reshape=False)
img_rotated = zoom(img_rotated, 1.5, order=2)
# Let's add a Gaussian PSF response with different seeing for both images
from scipy.ndimage.filters import gaussian_filter
img = gaussian_filter(img, sigma=2.0, mode='constant')
img_rotated = gaussian_filter(img_rotated, sigma=1.5, mode='constant')
# Let's add some noise to the images
noise_dc = 5.0
noise_std = np.sqrt(noise_dc)
img += np.random.normal(loc=noise_dc, scale=noise_std, size=img.shape)
img_rotated += np.random.normal(loc=noise_dc, scale=noise_std, size=img_rotated.shape)
img_aligned, footprint = aa.register(img, img_rotated, detection_sigma=3.0)
# Try this too!! Solving using only lists of (x, y) positions
# t, __ = aa.find_transform(img_posxy, img_rotated)
# img_aligned, footprint = aa.apply_transform(t, img, img_rotated)
fig, axes = plt.subplots(2, 2, figsize=(10, 10))
axes[0, 0].imshow(img, cmap='gray', interpolation='none', origin='lower')
axes[0, 0].axis('off')
axes[0, 0].set_title("Source Image")
axes[0, 1].imshow(img_rotated, cmap='gray', interpolation='none', origin='lower')
axes[0, 1].axis('off')
axes[0, 1].set_title("Target Image")
axes[1, 0].imshow(img_aligned, cmap='gray', interpolation='none', origin='lower')
axes[1, 0].axis('off')
axes[1, 0].set_title("Source Image aligned with Target")
axes[1, 1].imshow(footprint, cmap='gray', interpolation='none', origin='lower')
axes[1, 1].axis('off')
axes[1, 1].set_title("Footprint of the transformation")
axes[1, 0].axis('off')
plt.tight_layout()
plt.show()
In some occasions we need to check for the parameters of the transformation before applying it. Or maybe we just need to know the parameters for further analysis.
Another situation is where we don't have the image arrays, but we have lists of the (x, y) positions of the stars.
In these cases, we can call astroalign.find_transform
.
p, (pos_img, pos_img_rot) = aa.find_transform(img_posxy, img_rotated)
print("Rotation: {:.2f} degrees".format(p.rotation * 180.0 / np.pi))
print("\nScale factor: {:.2f}".format(p.scale))
print("\nTranslation: (x, y) = ({:.2f}, {:.2f})".format(*p.translation))
print("\nTranformation matrix:\n{}".format(p.params))
print("\nPoint correspondence:")
for (x1, y1), (x2, y2) in zip(pos_img, pos_img_rot):
print("({:.2f}, {:.2f}) is source --> ({:.2f}, {:.2f}) in target"
.format(x1, y1, x2, y2))
pos_img
and pos_img_rot
¶fig, axes = plt.subplots(2, 2, figsize=(10, 10))
colors = ['r', 'g', 'b', 'y', 'cyan', 'w', 'm']
axes[0, 0].imshow(img, cmap='gray', interpolation='none', origin='lower')
axes[0, 0].axis('off')
axes[0, 0].set_title("Source Image")
for (xp, yp), c in zip(pos_img[:len(colors)], colors):
circ = plt.Circle((xp, yp), 4, fill=False, edgecolor=c, linewidth=2)
axes[0, 0].add_patch(circ)
axes[0, 1].imshow(img_rotated, cmap='gray', interpolation='none', origin='lower')
axes[0, 1].axis('off')
axes[0, 1].set_title("Target Image")
for (xp, yp), c in zip(pos_img_rot[:len(colors)], colors):
circ = plt.Circle((xp, yp), 4 * p.scale, fill=False, edgecolor=c, linewidth=2)
axes[0, 1].add_patch(circ)
axes[1, 1].imshow(img_aligned, cmap='gray', interpolation='none', origin='lower')
axes[1, 1].axis('off')
axes[1, 1].set_title("Source Image aligned with Target")
for (xp, yp), c in zip(pos_img_rot[:len(colors)], colors):
circ = plt.Circle((xp, yp), 4 * p.scale, fill=False, edgecolor=c, linewidth=2)
axes[1, 1].add_patch(circ)
axes[1, 0].axis('off')
plt.tight_layout()
plt.show()