Source code for asterism.analysis_tasks.source_detection.deblending.denclue_deblending

"""
Overview
--------
This module provides the implementation of the :class:`Table`

Classes and Inheritance Structure
----------------------------------------------
.. inheritance-diagram::
    DoDENCLUEDeblendingTask



mary::
    Table


Summary
---------
.. autosummary::
    :nosignatures:

    DoDENCLUEDeblendingTask
    do_denclue_deblending
    do_denclue_of_cluster
    gauss_laplace_down_sampling

Module API
-----------
"""
from __future__ import division, absolute_import, print_function

import  numpy as np
from .deblending import  DeblendedProducts
from ....core.clustering.density_based.denclue import  DENCLUE
from ....pipeline_manager.analysis_tasks import AnalysisTask
from ....core.image_processing.features import BlobDetectionGaussLaplace
from ....core.geometry.distance import dist_eval
__author__ = 'andrea tramacere'

[docs]def denclue_plot(denclue_debl,parent_cluster,target_coordinates=None,target_coordinates_w_array=None): from asterism.plotting.plot_tools import show_image,plot_contour import pylab as plt if target_coordinates is not None: x_t=target_coordinates[:,0] y_t=target_coordinates[:,1] image,off_set_x,off_set_y,masked_pixels=parent_cluster.get_cluster_Image(border=1,bkg=0) im,ax=show_image(image.array) if denclue_debl.attractors_position_array is not None: x=denclue_debl.attractors_position_array[:,0] y=denclue_debl.attractors_position_array[:,1] dx=off_set_x dy=off_set_y ax.scatter(x-dx,y-dy,marker='o', edgecolor='black', linewidth='0',facecolor='white',s=100) plot_contour(ax,parent_cluster.contour_shape[:,0]-dx,parent_cluster.contour_shape[:,1]-dy,'black') if target_coordinates_w_array is not None: ax.scatter(x_t-dx,y_t-dy,marker='+',color='black',s=10) for cl in denclue_debl.clusters_list: cl.add_contour_points() if cl.contour_shape is not None: plot_contour(ax,cl.contour_shape[:,0]-dx,cl.contour_shape[:,1]-dy,'w') plt.plot(cl.x_c-dx,cl.y_c-dy,'+',ms=10) print("|dbs cl=%d, x_c=%f, y_c=%f"%(cl.ID,cl.x_c-dx,cl.y_c-dy)) plt.show()
[docs]def plot_final_clusters(parent_cluster,deblended_clusters_list): from asterism.plotting.plot_tools import plot_contour,show_image import pylab as plt image,off_set_x,off_set_y,masked_pixels=parent_cluster.get_cluster_Image(border=1,bkg=0) dx=off_set_x dy=off_set_y im,ax=show_image(image.array) for cl in deblended_clusters_list: if cl.contour_shape is not None: plot_contour(ax,cl.contour_shape[:,0]-dx,cl.contour_shape[:,1]-dy,'w') plt.show()
[docs]def get_pixel_matching_gauss_laplace_lm(parent_cluster,h,local_maxima_xy,off_x,off_y): lm_pix_coords=np.copy(local_maxima_xy) for ID,lm in enumerate(local_maxima_xy): #print ('x,y ',lm[0],lm[1],lm[0]+off_x,lm[1]+off_y) d=dist_eval(parent_cluster.cartesian.coords,x_c=lm[0]+off_x,y_c=lm[1]+off_y,metric='euclidean') msk1=d<h if parent_cluster.flux[msk1].size>0: lm_id=np.argmax(parent_cluster.flux[msk1]) #move back to cluster image coords, shift to cluster coords is done later #print ('f1',parent_cluster.flux[msk1].max()) lm_pix_coords[ID]=parent_cluster.cartesian.coords[msk1][lm_id]-[off_x,off_y] return lm_pix_coords
[docs]def gauss_laplace_down_sampling(parent_cluster,h,verbose=False,gl_th_rel=0.1,h_grid=3): """ Performs the down sampling of the target coordinates to use in the Denclue. * the cluster image is extracted * local maxima are extracted using the :class:`BlobDetectionGaussLaplace` * if number of local maxima is leq 1 then * target_coordinates are set to None * target_coordinates_w_array are set to None * if number of local maxima is geq 2 then * target_coordinates are obtained from a grid with size = max(1,h*0.5) * target_coordinates_w_array are obtained from the image flux at grid coordinates * coordinates and w_array corresponding to local maxima are added to the target arrays Parameters ---------- parent_cluster : h : float sets the side of the grid and the max_sigma=h, and min_sigma=h*0.5 used in the :class:`BlobDetectionGaussLaplace` Returns ------- target_coordinates : target_coordinates_w_array : """ bkg=parent_cluster.flux.min() image,off_x,off_y,masked=parent_cluster.get_cluster_image_array(bkg=bkg) h2=2.0*h h1=max(1.0,h2*0.5) #h_grid=max(1,np.int(h2*0.5)) if verbose == True: print("|downsampling start") print('|h=',h,) print('|h1,h2=',h1, h2) print('|h_grid=', h_grid) blob_log_det=BlobDetectionGaussLaplace(min_sigma=h1, max_sigma=h2, num_sigma=2, threshold_rel=gl_th_rel, overlap=0.5) local_maxima_xy,local_maxima_flux=blob_log_det.apply(image,coords_in_xy=True,verbose=verbose) #print ("local maxima",local_maxima_xy,local_maxima_flux) #print ("pixel matched local maxima",local_maxima_xy,local_maxima_flux) # Local maxima added to target point final list if local_maxima_flux is not None: local_maxima_xy = get_pixel_matching_gauss_laplace_lm(parent_cluster, h, local_maxima_xy, off_x, off_y) if verbose == True: print('|local maxima', local_maxima_flux.size) r=local_maxima_xy[:,1].astype(np.int16) c=local_maxima_xy[:,0].astype(np.int16) m=masked[r,c]==False r=r[m] c=c[m] x=c y=r else: x=None y=None #grid to extract target points with step given by kernel width step=max(1.0,np.int(h_grid)) row_sel=np.arange(0,image.shape[0],step,dtype=np.int16) col_sel=np.arange(0,image.shape[1],step,dtype=np.int16) mgrid=np.zeros(image.shape, dtype=np.bool) for row in row_sel: mgrid[row][col_sel]=True #target points are extracted from grid vertices #and appended to targets from local maxima y_grid,x_grid=np.where(np.logical_not(masked)*mgrid) if local_maxima_flux is not None: x_t=np.append(x,x_grid).astype(np.int16) y_t=np.append(y,y_grid).astype(np.int16) else: x_t=x_grid y_t=y_grid #fluxes are of target points are extracted target_coordinates_w_array=image[np.int16(y_t),np.int16(x_t)] #cluster image offset is added to target coordinates x_t=x_t+off_x y_t=y_t+off_y target_coordinates=np.column_stack((x_t.astype(np.int16),y_t.astype(np.int16))) if verbose == True: print("|downsampling stop") return target_coordinates,target_coordinates_w_array,local_maxima_xy,local_maxima_flux
[docs]def do_denclue_of_cluster(parent_cluster, min_size, kernel, eps, h_frac, h_min, attr_dbs_K, attr_dbs_eps, mask_unchanged, k_table_size, R_max_frac, R_max_kern_th, gl_th_rel=0.1, h_max=None, plot=False, digitize_attractors=False, gl_downsampling=False, h_grid=3, verbose=False): """ This function sets-up the parameters for the :class:`.DENCLUE` partitioning of the parent clusters. The following features are Parameters ---------- parent_cluster min_size kernel eps h_frac attr_dbs_K attr_dbs_eps mask_unchanged k_table_size R_max_frac R_max_kern_th plot digitize_attractors gl_downsampling Returns ------- """ print ("|debelnding for cluster",parent_cluster.ID, "r_cluter",parent_cluster.r_cluster,"r_max",parent_cluster.r_max,'size',parent_cluster.cl_len,min_size) if parent_cluster.cl_len>min_size: h=h_frac*np.sqrt(parent_cluster.r_max**2+parent_cluster.r_cluster**2) if h_min is not None: h=max(h,h_min) if h_max is not None: h=min(h,h_max) if R_max_frac is not None: R_max_method='fixed' R_max=R_max_frac*parent_cluster.r_cluster elif R_max_kern_th is not None: R_max_method='kern_rel_th' R_max=None else: R_max_method='auto' R_max=None if verbose==True: print ("|eps",eps) print ("|h",h) denclue_debl=DENCLUE(eps=eps,h=h,kernel=kernel,k_table_size=k_table_size,R_max_method=R_max_method,R_max_th_rel=R_max_kern_th,R_max=R_max) if gl_downsampling == True: target_coordinates,target_coordinates_w_array,local_maxima_xy,local_maxima_flux=gauss_laplace_down_sampling(parent_cluster,h,gl_th_rel,h_grid=h_grid) else: target_coordinates=None target_coordinates_w_array=None denclue_debl.run(parent_cluster.cartesian.coords, dbscan_eps=attr_dbs_eps, dbscan_K=attr_dbs_K, weight_array=parent_cluster.pts_weight, mask_unchanged=mask_unchanged, R_max=R_max, make_convolved_image=False, digitize_attractors=digitize_attractors, target_coordinates=target_coordinates, target_coordinates_w_array=target_coordinates_w_array, verbose=verbose) print("|denclue done") if plot is True: denclue_plot(denclue_debl,parent_cluster,target_coordinates=target_coordinates,target_coordinates_w_array=target_coordinates_w_array) children_clusters_list=denclue_debl.clusters_list if len(children_clusters_list)==1: children_clusters_list=[] else: children_clusters_list=[] return children_clusters_list
[docs]def do_denclue_deblending(clusters_list, eps=0.1, h_frac=0.20, h_min=1.0, h_max=None, kernel='gauss', gl_downsampling=False, h_grid=3, gl_th_rel=0.1, attr_dbs_eps=1.0, attr_dbs_K=4.0, min_size=9, mask_unchanged=0.5, k_table_size=None, R_max_frac=1.0, R_max_kern_th=None, plot=False, verbose=False, digitize_attractors=False): """ This function implements the top-level algorithm for the Denclue-based deblending: * each parent cluster in the `cluster_list` is partitioned by the :func:`do_denclue_of_cluster` function. * the parent cluster with his children clusters are used to build :class:`DeblendedProducts` object * a list of class:`DeblendedProducts` object is returned Parameters ---------- clusters_list eps h_frac kernel gl_downsampling attr_dbs_eps attr_dbs_K children_min_frac_integ_flux children_min_frac_peak_flux children_min_frac_size children_bright_frac_peak_flux min_size mask_unchanged k_table_size R_max_frac R_max_kern_th plot denclue_catalog verbose do_denclue_deblending out_catalog digitize_attractors Returns ------- """ deblended_prod_list=[] for parent_cluster in clusters_list: print("|-") denclue_clusters=do_denclue_of_cluster(parent_cluster, min_size, kernel, eps, h_frac, h_min, attr_dbs_K, attr_dbs_eps, mask_unchanged, k_table_size, gl_th_rel=gl_th_rel, h_max=h_max, R_max_frac=R_max_frac, R_max_kern_th=R_max_kern_th, plot=plot, digitize_attractors=digitize_attractors, gl_downsampling=gl_downsampling, h_grid=h_grid, verbose=verbose) print("|-") #if plot==True: # plot_final_clusters(parent_cluster,denclue_clusters) deblended_prod_list.append(DeblendedProducts(parent=parent_cluster, children_list=denclue_clusters,)) return deblended_prod_list
[docs]class DoDENCLUEDeblendingTask(AnalysisTask): """ """ def __init__(self, name='denclue_deblending', func=do_denclue_deblending, parser=None, process=None): super(DoDENCLUEDeblendingTask,self).__init__(name=name,func=func,parser=parser,process=process) self.parameters_list.add_par('-gl_downsampling',help='if True, a GuassLaplace downsampling of the target coordiantes is performed ',action='store_true') self.parameters_list.add_par('-h_grid', type=np.int, help=' grid mesh in pixels for the downsampling', default=3) self.parameters_list.add_par('-gl_th_rel', type=np.float, help='relative threshold for the GaussLaplace local maxima detection', default=0.1) self.parameters_list.add_par('-eps', type=np.float ,help='sets the stop threshold for the recursive attractor update rule',default=0.01) self.parameters_list.add_par('-digitize_attractors', help='if True, will approximate the final attractor position to the coordinate of the closest image pixel ',action='store_true') self.parameters_list.add_par('-kernel',type=str, help='smoothing kernel',default='gauss') self.parameters_list.add_par('-k_table_size',type=np.int, help='if not None, the kernel is precomputed in a look-up table over a grid with size equal to k_table_size', default=None) self.parameters_list.add_par('-h_frac', type=np.float ,help='sets the width of the kernel as h=h_frac*sqrt(r_max^2+r_cluster^2)',default=0.1 ) self.parameters_list.add_par('-h_min', type=np.float ,help='kernel width obtained by h_frac can not be lower than h_min',default=1.0 ) self.parameters_list.add_par('-h_max', type=np.float ,help='kernel width obtained by h_frac can not be lerger than h_max',default=None ) self.parameters_list.add_par('-mask_unchanged', type=np.float ,help='th value to flag and attractor as unchanged',default=None ) self.parameters_list.add_par('-min_size',type=np.int, help=' sets the minimum size in pixels to perform a deblending', default=9) self.parameters_list.add_par('-attr_dbs_eps', type=np.float ,help='eps for dbscan of attractors',default=1.0 ) self.parameters_list.add_par('-attr_dbs_K', type=np.float ,help='K for dbscan of attractors',default=4.0 ) self.parameters_list.add_par('-verbose',help='set verbosity on',action='store_true') self.parameters_list.add_par('-plot',help='plot ',action='store_true') group_method=self.parameters_list.add_mex_group('R_max_method') group_method.add_par('-R_max_frac', type=np.float ,help='max kernel influence radius as fraction of parent cluster r_max',default=None ) group_method.add_par('-R_max_kern_th', type=np.float, help='max kernel influence radius when kernel<=th with th relative to kenerl max value',default=None) self.func=func