Import stuff



In [1]:

    
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import fitsio
import deblend

Read in big file. Image is first extension, catalog is second extension. Also, fill in the uninitialized 'visible' column.



In [2]:

    
fits = fitsio.FITS("/Users/josh/src/WeakLensingDeblending/LSST_i.fits")
img = fits[0].read()
data = fits[1].read()
rad = 2048*0.2
data['visible'] = ((abs(data['dx']) < rad) & (abs(data['dy']) < rad))

We'll access items by group, so make a uniqified group array. Also, get the group sizes so we can ignore the insanely large blends.



In [3]:

    
grp_ids, grp_id_idx = np.unique(data['grp_id'][:], return_index=True)
grp_sizes = data['grp_size'][grp_id_idx]

Some functions to grab input parameters for the deblender given a group id.



In [4]:

    
def get_obj_ids(grp_id, data):
    return np.nonzero(data['grp_id']==grp_id)

def get_deblend_input(grp_id, data, img):
    w = get_obj_ids(grp_id, data)
    xmin = np.min(data['xmin'][w])
    xmax = np.max(data['xmax'][w])
    ymin = np.min(data['ymin'][w])
    ymax = np.max(data['ymax'][w])
    stamp = img[ymin:ymax, xmin:xmax]
    centroids = [(a['dx']/0.2+2048-xmin, a['dy']/0.2+2048-ymin) for a in data[w]]
    return stamp, centroids

Plot a particularly blendy blend.



In [5]:

    
stamp, centroids = get_deblend_input(grp_ids[7], data, img)
fig=plt.figure(figsize=(15,9))
plt.imshow(np.log10(stamp))
plt.colorbar()
plt.scatter(*map(list, zip(*centroids)), marker='+', color='k', s=50)
plt.show()
fig=plt.figure(figsize=(15,9))
plt.imshow(stamp, vmax=10000)
plt.colorbar()
plt.scatter(*map(list, zip(*centroids)), marker='+', color='k', s=50)









    



-c:3: RuntimeWarning: divide by zero encountered in log10






    












    Out[5]:





<matplotlib.collections.PathCollection at 0x1146c1ad0>

Deblend by group.



In [6]:

    
def deblend_group(grp_id, data, img, plot=False, log=True):
    stamp, centroids = get_deblend_input(grp_id, data, img)
    s = stamp.shape
    db_centroids = [(c[0]-s[0]*0.5, c[1]-s[1]*0.5) for c in centroids]
    _, _, children = deblend.deblend(stamp, db_centroids)

    if plot:
        nplots = len(children)+1
        nx = 5
        ny = np.ceil(1.*nplots/nx)
        fig = plt.figure(figsize=(18,18))
        ax = fig.add_subplot(ny, nx, 1)
        if log:
            vmax = np.log10(stamp.max())
            vmin = vmax-4
            im = ax.imshow(np.log10(stamp), vmin=vmin, vmax=vmax)
        else:
            vmax = stamp.max()
            vmin = 1e-4*vmax
            im = ax.imshow(stamp, vmin=vmin, vmax=vmax)
        ax.scatter(*map(list, zip(*centroids)), marker='+', color='k', s=300)

        for i, child in enumerate(children):
            ax = fig.add_subplot(ny, nx, i+2)
            if log:
                ax.imshow(np.log10(child), vmin=vmin, vmax=vmax)
            else:
                ax.imshow(child, vmin=vmin, vmax=vmax)
            ax.scatter(*centroids[i], marker='+', color='k', s=300)
        
    return children

good = np.nonzero((grp_sizes>=2) & (grp_sizes<=4))[0]
for g in good[0:10]:
    w = get_obj_ids(grp_ids[g], data)
    if not data[w]['visible'].all():
        continue
    deblend_group(grp_ids[g], data, img, plot=True, log=True)









    



deblend.py:33: RuntimeWarning: invalid value encountered in divide
  template_fractions = [template/template_sum * (template_sum != 0) for template in templates]
-c:16: RuntimeWarning: divide by zero encountered in log10
-c:26: RuntimeWarning: divide by zero encountered in log10



In [ ]: