Tutorial - Part #4 - Image Subtraction

For image subtraction the package has a module called propersubtract, which implements a main diff function.

[1]:
from copy import copy
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.colors as colors

%matplotlib inline
[2]:
from astropy.visualization import LinearStretch, LogStretch
from astropy.visualization import ZScaleInterval, MinMaxInterval
from astropy.visualization import ImageNormalize
[3]:
palette = copy(plt.cm.gray)
palette.set_bad('r', 1.0)
#palette.set_under('r', 1.0)
#palette.set_over('r', 1.0)
[4]:
import properimage.single_image as si
import properimage.propersubtract as ps
using numpy fft API
[5]:
ref_path = './../../../data/aligned_eso085-030-004.fit'
new_path = './../../../data/aligned_eso085-030-005.fit'

To get the subtraction we need to run this function by using both paths for example:

[6]:
from astropy.io.fits import getdata

plt.figure(figsize=(16, 12))
plt.subplot(121)
ref = getdata(ref_path)
norm = ImageNormalize(ref, interval=ZScaleInterval(),
                      stretch=LinearStretch())

plt.imshow(ref, cmap=plt.cm.gray, norm=norm, interpolation='none')

plt.colorbar(orientation='horizontal')

plt.subplot(122)
ref = getdata(new_path)
norm = ImageNormalize(ref, interval=ZScaleInterval(),
                      stretch=LinearStretch())

plt.imshow(ref, cmap=plt.cm.gray, norm=norm, interpolation='none')

plt.colorbar(orientation='horizontal')
[6]:
<matplotlib.colorbar.Colorbar at 0x7effbf48d070>
../_images/tutorial_Tutorial04_8_1.png
[7]:
result = ps.diff(ref=ref_path, new=new_path, smooth_psf=False, fitted_psf=True,
                 align=False, iterative=False, beta=False, shift=False)
using single psf, gaussian modeled
updating stamp shape to (21,21)
updating stamp shape to (21,21)
S_corr sigma_clipped_stats
mean = -4.210456578948465e-06, median = -4.005563176566472e-06, std = 0.00011940426523510884

Subtraction performed in 17.070908308029175 seconds


The result is a list of numpy arrays.

The arrays are in order: D, P, Scorr, mask

[8]:
D = result[0]
P = result[1]
Scorr = result[2]
mask = result[3]
[9]:
norm = ImageNormalize(D.real, interval=ZScaleInterval(),
                      stretch=LinearStretch())

plt.figure(figsize=(12, 12))
plt.imshow(np.ma.MaskedArray(D.real, mask=mask),
           cmap=palette, norm=norm, interpolation='none')

plt.colorbar(orientation='horizontal')
[9]:
<matplotlib.colorbar.Colorbar at 0x7effbf4e8880>
../_images/tutorial_Tutorial04_12_1.png
[10]:
norm = ImageNormalize(P, interval=MinMaxInterval(),
                      stretch=LogStretch())
xc, yc = np.where(P==P.max())
xc, yc = np.round(xc[0]), np.round(yc[0])
plt.imshow(P[xc-10:xc+10, yc-10:yc+10], norm=norm,
           cmap='viridis', interpolation='none')
plt.colorbar()
[10]:
<matplotlib.colorbar.Colorbar at 0x7effbf381460>
../_images/tutorial_Tutorial04_13_1.png
[11]:
norm = ImageNormalize(Scorr.real, interval=ZScaleInterval(),
                      stretch=LinearStretch())

plt.figure(figsize=(12, 12))
plt.imshow(np.ma.MaskedArray(Scorr.real, mask=mask),
           cmap=palette, norm=norm, interpolation='none')

plt.colorbar(orientation='horizontal')
[11]:
<matplotlib.colorbar.Colorbar at 0x7effbf244700>
../_images/tutorial_Tutorial04_14_1.png
[12]:
plt.hist(np.ma.MaskedArray(D.real, mask=mask).filled(0).flatten(), log=True, bins=500)
#plt.xlim(-10, 10)
plt.show()
../_images/tutorial_Tutorial04_15_0.png
[13]:
simage = np.ma.MaskedArray(Scorr.real, mask=mask).filled(0).flatten()
plt.hist(simage/np.std(simage), log=True, bins=500)
#plt.xlim(-4, 4)
plt.show()
../_images/tutorial_Tutorial04_16_0.png

This is related to the quantities derived in Zackay et al. works. \(S_{corr} = P_D \otimes D\)

[14]:
D, P, Scorr, mask = ps.diff(ref=ref_path, new=new_path,
                            align=False, iterative=False, beta=True)
using single psf, gaussian modeled
updating stamp shape to (21,21)
updating stamp shape to (21,21)
Found that beta = [ 0.98135014  0.04777034 -0.02548551]
Took only 13.118748426437378 awesome seconds
The solution was with cost 0.3210708100934259
S_corr sigma_clipped_stats
mean = -2.786345837641917e-06, median = -2.857150869955863e-06, std = 0.00011675875233117702

Subtraction performed in 32.61811280250549 seconds


The result is a list of numpy arrays.

The arrays are in order: D, P, Scorr, mask

[15]:
norm = ImageNormalize(D.real, interval=ZScaleInterval(),
                      stretch=LinearStretch())

plt.figure(figsize=(12, 12))
plt.imshow(np.ma.MaskedArray(D.real, mask=mask),
           cmap=palette, norm=norm, interpolation='none')
plt.colorbar(orientation='horizontal')
[15]:
<matplotlib.colorbar.Colorbar at 0x7effbc13cd90>
../_images/tutorial_Tutorial04_20_1.png
[16]:
norm = ImageNormalize(P, interval=MinMaxInterval(),
                      stretch=LogStretch())
xc, yc = np.where(P==P.max())
xc, yc = np.round(xc[0]), np.round(yc[0])
plt.imshow(P[xc-10:xc+10, yc-10:yc+10], norm=norm,
           cmap='viridis', interpolation='none')
plt.colorbar()
[16]:
<matplotlib.colorbar.Colorbar at 0x7effbc405ca0>
../_images/tutorial_Tutorial04_21_1.png
[17]:
norm = ImageNormalize(Scorr.real, interval=ZScaleInterval(),
                      stretch=LinearStretch())

plt.figure(figsize=(12, 12))
plt.imshow(np.ma.MaskedArray(Scorr.real, mask=mask),
           cmap=palette, norm=norm, interpolation='none')

plt.colorbar(orientation='horizontal')
[17]:
<matplotlib.colorbar.Colorbar at 0x7effb74757f0>
../_images/tutorial_Tutorial04_22_1.png
[18]:
plt.hist(np.ma.MaskedArray(D.real, mask=mask).filled(0).flatten(), log=True, bins=500)
#plt.xlim(-10, 10)
plt.show()
../_images/tutorial_Tutorial04_23_0.png
[19]:
simage = np.ma.MaskedArray(Scorr.real, mask=mask).filled(0).flatten()
plt.hist(simage/np.std(simage), log=True, bins=500)
#plt.xlim(-4, 4)
plt.show()
../_images/tutorial_Tutorial04_24_0.png

We have the option of using the iterative methods without beta

[20]:
D, P, Scorr, mask = ps.diff(ref=ref_path, new=new_path, align=False, iterative=True, beta=False)
using single psf, gaussian modeled
updating stamp shape to (21,21)
updating stamp shape to (21,21)
Found that shift = [ 0.04702745 -0.02644272]
Took only 7.601429224014282 awesome seconds
The solution was with cost 0.32143964639703293
S_corr sigma_clipped_stats
mean = -4.25172462215914e-06, median = -4.014162635757632e-06, std = 0.00011881913562472867

Subtraction performed in 23.768754482269287 seconds


The result is a list of numpy arrays.

The arrays are in order: D, P, Scorr, mask

[21]:
norm = ImageNormalize(D.real, interval=ZScaleInterval(),
                      stretch=LinearStretch())

plt.figure(figsize=(12, 12))
plt.imshow(np.ma.MaskedArray(D.real, mask=mask),
           cmap=palette, norm=norm, interpolation='none')
plt.colorbar(orientation='horizontal')
[21]:
<matplotlib.colorbar.Colorbar at 0x7effbc0e4550>
../_images/tutorial_Tutorial04_28_1.png
[22]:
norm = ImageNormalize(P, interval=MinMaxInterval(),
                      stretch=LogStretch())
xc, yc = np.where(P==P.max())
xc, yc = np.round(xc[0]), np.round(yc[0])
plt.imshow(P[xc-10:xc+10, yc-10:yc+10], norm=norm,
           cmap='viridis', interpolation='none')
plt.colorbar()
[22]:
<matplotlib.colorbar.Colorbar at 0x7effbc170df0>
../_images/tutorial_Tutorial04_29_1.png
[23]:
norm = ImageNormalize(Scorr.real, interval=ZScaleInterval(),
                      stretch=LinearStretch())

plt.figure(figsize=(12, 12))
plt.imshow(np.ma.MaskedArray(Scorr.real, mask=mask),
           cmap=palette, norm=norm, interpolation='none')

plt.colorbar(orientation='horizontal')
[23]:
<matplotlib.colorbar.Colorbar at 0x7effbc7070d0>
../_images/tutorial_Tutorial04_30_1.png
[24]:
plt.hist(np.ma.MaskedArray(D.real, mask=mask).filled(0).flatten(), log=True, bins=500)
#plt.xlim(-10, 10)
plt.show()
../_images/tutorial_Tutorial04_31_0.png
[25]:
simage = np.ma.MaskedArray(Scorr.real, mask=mask).filled(0).flatten()
plt.hist(simage/np.std(simage), log=True, bins=500)
#plt.xlim(-4, 4)
plt.show()
../_images/tutorial_Tutorial04_32_0.png

We have the option of using the iterative methods without beta

[26]:
D, P, Scorr, mask = ps.diff(ref=ref_path, new=new_path, smooth_psf=False, fitted_psf=True,
                            align=False, iterative=True, beta=True, shift=True)
using single psf, gaussian modeled
updating stamp shape to (21,21)
updating stamp shape to (21,21)
Found that beta = [ 0.98135014  0.04777034 -0.02548551]
Took only 9.965375185012817 awesome seconds
The solution was with cost 0.3210708100934259
S_corr sigma_clipped_stats
mean = -2.786345837641917e-06, median = -2.857150869955863e-06, std = 0.00011675875233117702

Subtraction performed in 24.975581407546997 seconds


The result is a list of numpy arrays.

The arrays are in order: D, P, Scorr, mask

[27]:
norm = ImageNormalize(D.real, interval=ZScaleInterval(),
                      stretch=LinearStretch())

plt.figure(figsize=(22, 22))
plt.imshow(np.ma.MaskedArray(D.real, mask=mask),
           cmap=palette, norm=norm, interpolation='none')
plt.colorbar(orientation='horizontal')
[27]:
<matplotlib.colorbar.Colorbar at 0x7effbf2f07c0>
../_images/tutorial_Tutorial04_36_1.png
[28]:
norm = ImageNormalize(P, interval=MinMaxInterval(),
                      stretch=LogStretch())
xc, yc = np.where(P==P.max())
xc, yc = np.round(xc[0]), np.round(yc[0])
plt.imshow(P[xc-10:xc+10, yc-10:yc+10], norm=norm,
           cmap='viridis', interpolation='none')
plt.colorbar()
[28]:
<matplotlib.colorbar.Colorbar at 0x7effbf43dc10>
../_images/tutorial_Tutorial04_37_1.png
[29]:
norm = ImageNormalize(Scorr.real, interval=ZScaleInterval(),
                      stretch=LinearStretch())

plt.figure(figsize=(12, 12))
plt.imshow(np.ma.MaskedArray(Scorr.real, mask=mask),
           cmap=palette, norm=norm, interpolation='none')

plt.colorbar(orientation='horizontal')
[29]:
<matplotlib.colorbar.Colorbar at 0x7effb7654700>
../_images/tutorial_Tutorial04_38_1.png
[30]:
plt.hist(np.ma.MaskedArray(D.real, mask=mask).filled(0).flatten(), log=True, bins=500)
#plt.xlim(-10, 10)
plt.show()
../_images/tutorial_Tutorial04_39_0.png
[31]:
simage = np.ma.MaskedArray(Scorr.real, mask=mask).filled(0).flatten()
plt.hist(simage/np.std(simage), log=True, bins=500)
#plt.xlim(-10, 4)
plt.show()
../_images/tutorial_Tutorial04_40_0.png