This is a primitive track finder in CCD images. It thresholds the image to find "blobs" above threshold and then loops through all found blobs to group them according to their separation distance.
If one really wants to look for line-like features to ID muon tracks, a better approach for grouping blobs could be something like a Hough transform.
A simple calculation of the principal moments of the blob groups is available in the classes provided but is not used in this notebook. Those calculations are based on the paper "Cosmic rays and other nonsense in astronomical CCD imagers" by Don Groom, Exp. Astron. 14:55, 2002. The code was written for the DECO high-school outreach program at UW-Madison in 2013.
In [2]:
from astropy.io import fits
import numpy as np
import matplotlib.pyplot as plt
from skimage import measure
from astropy.visualization import astropy_mpl_style
plt.style.use(astropy_mpl_style)
In [3]:
class Blob:
"""Class that defines a 'blob' in an image: the contour of a set of pixels
with values above a given threshold.
"""
def __init__(self, x, y):
"""Define a counter by its contour lines (an list of points in the xy
plane), the contour centroid, and its enclosed area.
Parameters
----------
x : list or array_like
x-values of blob contour.
y : list or array_like
y-values of blob contour.
"""
self.x = x
self.y = y
self.xc = np.mean(x)
self.yc = np.mean(y)
# Find the area inside the contour
self.area = 0.
n = len(x)
for i in range(0, n):
self.area += 0.5*(y[i]+y[i-1])*(x[i]-x[i-1])
def distance(self, blob):
"""Calculate the distance between the centroid of this blob contour and
another one in the xy plane.
Parameters
----------
blob : Blob
A second blob.
Returns
-------
dist : float
Euclidean distance between two blob centroids.
"""
return np.sqrt((self.xc - blob.xc)**2 + (self.yc-blob.yc)**2)
In [4]:
class BlobGroup:
"""A list of blobs that is grouped or associated in some way, i.e., if
their contour centroids are relatively close together."""
def __init__(self):
"""Initialize a list of stored blobs and the bounding rectangle which
defines the group.
"""
self.blobs = []
self.xmin = 1e10
self.xmax = -1e10
self.ymin = 1e10
self.ymax = -1e10
def addBlob(self, blob):
"""Add a blob to the group and enlarge the bounding rectangle of the
group.
"""
self.blobs.append(blob)
self.xmin = min(self.xmin, blob.x.min())
self.xmax = max(self.xmax, blob.x.max())
self.ymin = min(self.ymin, blob.y.min())
self.ymax = max(self.ymax, blob.y.max())
self.cov = None
def getBoundingBox(self):
"""Get the bounding rectangle of the group."""
return (self.xmin, self.xmax, self.ymin, self.ymax)
def getSquareBoundingBox(self):
"""Get the bounding rectangle, redefined to give it a square aspect
ratio."""
xmin, xmax, ymin, ymax = (self.xmin, self.xmax, self.ymin, self.ymax)
xL = np.abs(xmax - xmin)
yL = np.abs(ymax - ymin)
if xL > yL:
ymin -= 0.5*(xL-yL)
ymax += 0.5*(xL-yL)
else:
xmin -= 0.5*(yL-xL)
xmax += 0.5*(yL-xL)
return (xmin, xmax, ymin, ymax)
def getSubImage(self, image):
"""Given an image, extract the section of the image corresponding to
the bounding box of the blob group.
"""
ny,nx = image.shape
x0,x1,y0,y1 = self.getBoundingBox()
# Account for all the weird row/column magic in the image table...
i0,i1 = [ny - int(t) for t in (y1,y0)]
j0,j1 = [int(t) for t in (x0,x1)]
# Add a pixel buffer around the bounds, and check the ranges
buf = 1
i0 = 0 if i0-buf < 0 else i0-buf
i1 = ny-1 if i1 > ny-1 else i1+buf
j0 = 0 if j0-buf < 0 else j0-buf
j1 = nx-1 if j1 > nx-1 else j1+buf
return image[i0:i1, j0:j1]
def getRawMoment(self, image, p, q):
"""Calculate the image moment given by
M_{ij}=\sum_x\sum_y x^p y^q I(x,y)
where I(x,y) is the image intensity at location x,y.
"""
nx,ny = image.shape
Mpq = 0.
if p == 0 and q == 0:
Mpq = np.sum(image)
else:
for i in range(0,nx):
x = 0.5 + i
for j in range(0,ny):
y = 0.5 + j
Mpq += x**p * y**q * image[i,j]
return Mpq
def getCovariance(self, image):
"""Get the raw moments of the image region inside the bounding box
defined by this blob group and calculate the image covariance
matrix."""
if self.cov is None:
subImage = self.getSubImage(image).transpose()
M00 = self.getRawMoment(subImage, 0, 0)
M10 = self.getRawMoment(subImage, 1, 0)
M01 = self.getRawMoment(subImage, 0, 1)
M11 = self.getRawMoment(subImage, 1, 1)
M20 = self.getRawMoment(subImage, 2, 0)
M02 = self.getRawMoment(subImage, 0, 2)
xbar = M10/M00
ybar = M01/M00
self.cov = np.vstack([[M20/M00 - xbar*xbar, M11/M00 - xbar*ybar],
[M11/M00 - xbar*ybar, M02/M00 - ybar*ybar]])
return self.cov
def getPrincipalMoments(self, image):
"""Return the maximum and minimum eigenvalues of the covariance matrix,
as well as the angle theta between the maximum eigenvector and the
x-axis."""
cov = self.getCovariance(image)
u20 = cov[0,0]
u11 = cov[0,1]
u02 = cov[1,1]
theta = 0.5 * np.arctan2(2*u11, u20-u02)
l1 = 0.5*(u20+u02) + 0.5*np.sqrt(4*u11**2 + (u20-u02)**2)
l2 = 0.5*(u20+u02) - 0.5*np.sqrt(4*u11**2 + (u20-u02)**2)
return l1, l2, theta
Two functions to find blobs using a threshhold.
findBlobs uses the marching squares algorithm to identify a blob's outer contours.
groupBlobs is called after findBlobs and just loops over all blobs ID-ed in the image and agglomerates anything within a given distance threshold of another blob.
In [5]:
def findBlobs(image, threshold, minArea=2.):
"""Pass through an image and find a set of blobs/contours above a set
threshold value. The minArea parameter is used to exclude blobs with an
area below this value.
"""
blobs = []
ny, nx = image.shape
# Find contours using the Marching Squares algorithm in the scikit package.
contours = measure.find_contours(image, threshold)
for contour in contours:
x = contour[:,1]
y = ny - contour[:,0]
blob = Blob(x, y)
if blob.area >= minArea:
blobs.append(blob)
return blobs
In [6]:
def groupBlobs(blobs, maxDist):
"""Given a list of blobs, group them by distance between the centroids of
any two blobs. If the centroids are more distant than maxDist, create a
new blob group.
"""
n = len(blobs)
groups = []
if n >= 1:
# Single-pass clustering algorithm: make the first blob the nucleus of
# a blob group. Then loop through each blob and add either add it to
# this group (depending on the distance measure) or make it the
# nucleus of a new blob group
bg = BlobGroup()
bg.addBlob(blobs[0])
groups.append(bg)
for i in range(1, n):
bi = blobs[i]
isGrouped = False
for group in groups:
# Calculate distance measure for a blob and a blob group:
# blob just has to be < maxDist from any other blob in the group
for bj in group.blobs:
if bi.distance(bj) < maxDist:
group.addBlob(bi)
isGrouped = True
break
if not isGrouped:
bg = BlobGroup()
bg.addBlob(bi)
groups.append(bg)
return groups
In [7]:
hdus = fits.open('/global/project/projectdirs/desi/spectro/redux/daily/preproc/20191208/00031136/preproc-z3-00031136.fits')
for hdu in hdus:
print(hdu.header['EXTNAME'])
In [8]:
img = hdus['IMAGE'].data
mask = hdus['MASK'].data
readnoise = hdus['READNOISE'].data
In [9]:
plt.subplots(1,1, figsize=(12,9), tight_layout=True)
plt.imshow(img, cmap='gray', vmin=0, vmax=2000)
plt.colorbar()
Out[9]:
In [10]:
plt.subplots(1,1, figsize=(12,9), tight_layout=True)
plt.imshow(mask, cmap='gray', vmin=0, vmax=1)
plt.colorbar()
Out[10]:
ID blobs as regions above a threshold, and choose a minimum area to avoid identifying noisy individual pixels. The threshold used was tuned (barely) by visual inspection of the result.
The group blobs according to the Euclidean distance between them, up to some maximum distance in pixel space. Again, this was tuned in a very crude way by looking at the output.
In [11]:
blobs = findBlobs(mask, threshold=0.5, minArea=2)
groups = groupBlobs(blobs, maxDist=30.)
In [12]:
plt.subplots(1,1, figsize=(12,9), tight_layout=True)
plt.imshow(mask, cmap='gray', vmin=0, vmax=1)
for blob in blobs:
plt.plot(blob.x, blob.y, linewidth=2, color='#00dd00')
plt.colorbar()
Out[12]:
In [13]:
plt.subplots(1,1, figsize=(12,9), tight_layout=True)
plt.imshow(mask, cmap='gray', vmin=0, vmax=1.)
for i, group in enumerate(groups):
if len(group.blobs) > 5:
for blob in group.blobs:
plt.plot(blob.x, blob.y, linewidth=2, color='#00dd00')
plt.colorbar()
Out[13]:
In [ ]: