# adds parent dir to python path
import sys
sys.path.insert(0, '..')
import os
import numpy as np
from astropy.io import fits
import matplotlib.pyplot as plt
%matplotlib inline
import ImageTools as it
The Petrosian Radius is the R at which 1η(R)=0.2, where
1η(R)=πr2∑RrI(r)I(R)
def one_over_eta(rs, fs, R):
return ((np.pi*R**2) / (np.sum(fs[rs<=R]))) * fs[rs==R]
def petrosian_radius(rs, fs):
r_candidates = np.array([one_over_eta(rs, fs, R) for R in rs])
return rs[np.square(r_candidates-0.2).argmin()]
The Petrosian Flux is defined as: Fp=2Rp∑rI(r) Where Rp is the Petrosian Radius
def petrosian_flux(rs, fs, R_p):
return np.sum(fs[rs<=2*R_p])
Defining the effective radius using the Petrosian Flux
Re≈R501−P3(R90/R50)P4
Where:
Rx=x100=∑RxrI(r)∑2RprI(r),P3=8∗10−6,P4=8.47
def R_x(rs, fs, x):
x /= 100
sum_ratio = np.cumsum(fs)/np.sum(fs)
return rs[np.square(sum_ratio-x).argmin()]
def petrosian_Re(R_50, R_90):
P_3 = 8e-6
P_4 = 8.47
return R_50 / (1 - P_3*(R_90/R_50)**P_4)
home = os.getenv('HOME')
img_dir = os.path.join(home, 'Documents/astro_data/orig_images')
def get_img(_id, band):
file_name = 'GDS_{}_{}.fits'.format(_id, band)
return fits.getdata(os.path.join(img_dir, file_name))
img = get_img('deep2_10042', 'h')
segmap = get_img('deep2_10042', 'segmap')
old_re = it.effective_radius(img, segmap==10042)
plt.imshow(img, cmap='gray')
cx, cy = it.img_center(img, segmap==10042)
xs, ys = np.meshgrid(np.arange(img.shape[0]), np.arange(img.shape[1]).T)
rs = np.sqrt(np.square(xs-cx)+np.square(ys-cy))
rs = rs.flatten()
fs = img.flatten()
sorted_rs = np.argsort(rs)
rs = rs[sorted_rs]
fs = fs[sorted_rs]
R_p = petrosian_radius(rs, fs)
F_p = petrosian_flux(rs, fs, R_p)
p_mask = rs <= 2*R_p
rs = rs[p_mask]
fs = fs[p_mask]
R_50 = R_x(rs, fs, 50)
R_90 = R_x(rs, fs, 90)
R_e = petrosian_Re(R_50, R_90)
print(f'Petrosian Radius:{R_p}')
print(f'Petrosian Flux:{fp}')
print(f'R_50:{R_50}')
print(f'R_90:{R_90}')
print(f'Effective Radius:{R_e}')
make_cir = lambda r,c: plt.Circle((cx, cy), r, color=c, fill=False)
cir_p_radius = make_cir(R_p, 'b')
cir_R_50 = make_cir(R_50, 'y')
cir_R_90 = make_cir(R_90, 'r')
cir_R_e = make_cir(R_e, 'g')
cir_old_re = make_cir(old_re, 'w')
f, a = plt.subplots(figsize=(15,15))
a.imshow(img, cmap='gray')
a.add_artist(cir_p_radius)
a.add_artist(cir_R_50)
a.add_artist(cir_R_90)
a.add_artist(cir_R_e)
a.add_artist(cir_old_re)
a.legend([cir_p_radius, cir_R_50, cir_R_90, cir_R_e, cir_old_re],
['Petrosian Radius', 'R_50', 'R_90', 'R_e', 'Old Re'],
fontsize=20)
plt.show()