In [21]:
import matplotlib.pyplot as plt
%matplotlib inline
import matplotlib as mpl
mpl.rcParams['figure.figsize'] = (10, 9)
In [22]:
file = '/Users/schriste/heroespy/data/HEROES_20130910_X-ray_alignment_repeat2/Detector0/position1/Det00_position1_s.evt'
In [23]:
from astropy.io import fits
f = fits.open(file)
In [24]:
f[0].header
Out[24]:
In [25]:
f[1].header['DETECTOR']
Out[25]:
In [26]:
f[1].data[1][2]
Out[26]:
In [27]:
data = f[1].data
In [28]:
len(data)
Out[28]:
In [29]:
data[0]
Out[29]:
In [30]:
f[1].header
Out[30]:
In [31]:
data['RAWX']
Out[31]:
In [32]:
data['RAWY']
Out[32]:
In [33]:
import numpy as np
img, xedges, yedges = np.histogram2d(data['RAWX'], data['RAWY'], bins=200)
In [34]:
plt.imshow(img)
Out[34]:
In [35]:
import skimage
print(skimage.__version__)
from skimage import data, filter, color
from skimage.transform import hough_circle
from skimage.feature import peak_local_max
from skimage.draw import circle_perimeter
from skimage.util import img_as_ubyte
In [36]:
print(img.max())
image = img_as_ubyte(img.astype(np.uint8))
edges = filter.canny(image, sigma=1)
In [37]:
plt.imshow(edges)
Out[37]:
In [38]:
hough_radii = np.arange(15, 30, 2)
hough_res = hough_circle(edges, hough_radii)
In [39]:
centers = []
accums = []
radii = []
for radius, h in zip(hough_radii, hough_res):
# For each radius, extract two circles
peaks = peak_local_max(h, num_peaks=2)
centers.extend(peaks)
accums.extend(h[peaks[:, 0], peaks[:, 1]])
radii.extend([radius, radius])
fig, ax = plt.subplots()
plt.imshow(image, cmap=plt.cm.gray)
collection = []
for idx in np.argsort(accums)[::-1][:5]:
center_x, center_y = centers[idx]
radius = radii[idx]
circle = plt.Circle((center_y, center_x), radius, color='r', fill=False)
ax.add_patch(circle)
plt.show()