diff --git a/.gitignore b/.gitignore index 5cc231d..59b314a 100644 --- a/.gitignore +++ b/.gitignore @@ -1,8 +1,8 @@ *.pyc *.pytest_cache *.vscode -Test_catlogue.fits -test_gpu.png + + *.svg *.fits *.o @@ -13,11 +13,12 @@ test_gpu.png *.output *.pdf *.h5 -bench_script.sh -codeDiagram.drawio -Lockman_facets_process.sh -DRUID.ouput -/Examples_backup + + + + + DRUID.egg-info build -backup \ No newline at end of file +backup +notepad.ipynb diff --git a/DRUID/main.py b/DRUID/main.py index b94709f..1459e38 100644 --- a/DRUID/main.py +++ b/DRUID/main.py @@ -1,19 +1,17 @@ - - -''' +""" File: main.py Author: Rhys Shaw Date: 23/12/2023 Version: 0.0 Description: Main file for DRUID -''' +""" - -version = '0.0.0' +version = "0.0-test" import setproctitle -setproctitle.setproctitle('DRUID') -from .src import utils +setproctitle.setproctitle("DRUID") + +from .src import utils from .src import homology_new as homology from .src import background from matplotlib import colors @@ -33,9 +31,6 @@ import logging - - - DRUID_MESSAGE = """ @@ -59,25 +54,33 @@ For more information see: https://github.com/RhysAlfShaw/DRUID - """.format(version) - - - - - - + """.format( + version +) class sf: - - - - def __init__(self, image : np.ndarray = None, image_path : str = None, mode : str = None, - pb_path : str = None, cutup : bool = False, cutup_size : int = 500, - cutup_buff : int = None, output : bool = True, - area_limit : int = 5, smooth_sigma = 1, nproc : int = 1, GPU : bool = False, - header : astropy.io.fits.header.Header = None, Xoff : int = None, Yoff : int = None,debug_mode=False,remove_edge=True) -> None: + def __init__( + self, + image: np.ndarray = None, + image_path: str = None, + mode: str = None, + pb_path: str = None, + cutup: bool = False, + cutup_size: int = 500, + cutup_buff: int = None, + output: bool = True, + area_limit: int = 5, + smooth_sigma=1, + nproc: int = 1, + GPU: bool = False, + header: astropy.io.fits.header.Header = None, + Xoff: int = None, + Yoff: int = None, + debug_mode=False, + remove_edge=True, + ) -> None: """Initialise the DRUID, here general parameters can be set.. Args: @@ -103,11 +106,11 @@ def __init__(self, image : np.ndarray = None, image_path : str = None, mode : st print(DRUID_MESSAGE) if debug_mode: - logging.basicConfig(format='%(asctime)s - %(message)s', level=logging.DEBUG) - logging.debug('Debug mode enabled') + logging.basicConfig(format="%(asctime)s - %(message)s", level=logging.DEBUG) + logging.debug("Debug mode enabled") else: - logging.basicConfig(format='%(asctime)s - %(message)s', level=logging.INFO) - + logging.basicConfig(format="%(asctime)s - %(message)s", level=logging.INFO) + self.cutup = cutup self.output = output self.image_path = image_path @@ -118,7 +121,7 @@ def __init__(self, image : np.ndarray = None, image_path : str = None, mode : st self.Yoff = Yoff self.cutup_buff = cutup_buff self.remove_edge = remove_edge - + if self.GPU: # Lets try importing the GPU stuff, if it fails then we can just use the CPU. @@ -126,64 +129,65 @@ def __init__(self, image : np.ndarray = None, image_path : str = None, mode : st try: import cupy as cp from cupyx.scipy.ndimage import label as cp_label - + # num_gpus = cp.cuda.runtime.getDeviceCount() - #print(f'Found {num_gpus} GPUs') - + # print(f'Found {num_gpus} GPUs') + # if num_gpus > 0: # #print('GPUs are avalible, GPU functions will now be avalible.') - # GPU_AVALIBLE = True + # GPU_AVALIBLE = True # else: - #print('No GPUs avalible, using CPU') + # print('No GPUs avalible, using CPU') # GPU_AVALIBLE = False except: - - #print('Could not import cupy. DRUID GPU functions will not be avalible') + + # print('Could not import cupy. DRUID GPU functions will not be avalible') GPU_AVALIBLE = False - + self.nproc = nproc - + self.header = header - + if self.image_path is None: self.image = image - + else: self.image, self.header = utils.open_image(self.image_path) - - if mode == 'Radio': - # add 1 pixel padding to all sides of the image this causes issues? - self.image = np.pad(self.image,((1,1),(1,1)),mode='constant',constant_values=0) - - - if self.smooth_sigma !=0: - self.image = utils.smoothing(self.image,self.smooth_sigma) - logging.info('Image smoothed with sigma = {}'.format(self.smooth_sigma)) - + + if mode == "Radio": + pass + # self.image = np.pad( + # self.image, ((1, 1), (1, 1)), mode="constant", constant_values=0 + # ) + + if self.smooth_sigma != 0: + self.image = utils.smoothing(self.image, self.smooth_sigma) + logging.info("Image smoothed with sigma = {}".format(self.smooth_sigma)) + self.mode = mode - # check if the mode is valid - - if self.mode not in ['Radio', 'optical','other']: - raise ValueError('Mode must be either radio, optical or other.') - - + + if self.mode not in ["Radio", "optical", "other"]: + raise ValueError("Mode must be either radio, optical or other.") + self.pb_path = pb_path - + if self.pb_path is not None: - self.pb_image , self.pb_header = utils.open_image(self.pb_path) - - + self.pb_image, self.pb_header = utils.open_image(self.pb_path) + if self.cutup: self.cutup_size = cutup_size - #print('Cutting up image {}x{} into {}x{} cutouts'.format(self.image.shape[0],self.image.shape[1],cutup_size,cutup_size)) if self.cutup_buff is not None: - self.cutouts, self.coords = utils.cut_image_buff(self.image,cutup_size, buffer_size=self.cutup_buff) + self.cutouts, self.coords = utils.cut_image_buff( + self.image, cutup_size, buffer_size=self.cutup_buff + ) else: - self.cutouts, self.coords = utils.cut_image(cutup_size,self.image) - + self.cutouts, self.coords = utils.cut_image(cutup_size, self.image) + if self.pb_path is not None: - self.pb_cutouts, self.pb_coords = utils.cut_image(cutup_size,self.pb_image) - + self.pb_cutouts, self.pb_coords = utils.cut_image( + cutup_size, self.pb_image + ) + else: self.pb_cutouts = None self.pb_coords = None @@ -193,79 +197,94 @@ def __init__(self, image : np.ndarray = None, image_path : str = None, mode : st self.coords = None self.pb_cutouts = None self.pb_coords = None - - - - - - - - - + def phsf(self, lifetime_limit: float = 0, lifetime_limit_fraction: float = 2): + """Performs the persistent homology source finding algorithm. - - - def phsf(self, lifetime_limit : float = 0,lifetime_limit_fraction : float = 2): - - """ Performs the persistent homology source finding algorithm. - Args: - + lifetime_limit (float): The lifetime limit for the persistent homology algorithm. - + Returns: - + None. The catalogue is stored in the self.catalogue attribute. - + """ - - + if self.cutup == True: - + catalogue_list = [] IDoffset = 0 - for i, cutout in tqdm(enumerate(self.cutouts),total=len(self.cutouts),desc='Processing Cutouts',disable=not self.output): - - print("Computing for Cutout number :{}/{}".format(i+1, len(self.cutouts))) - - catalogue = homology.compute_ph_components(cutout,self.local_bg,analysis_threshold_val=self.analysis_threshold_val, - lifetime_limit=lifetime_limit,output=self.output,bg_map=self.bg_map,area_limit=self.area_limit, - GPU=self.GPU,lifetime_limit_fraction=lifetime_limit_fraction,mean_bg=self.mean_bg, - IDoffset=IDoffset,box_size=self.box_size,detection_threshold=self.sigma) + for i, cutout in tqdm( + enumerate(self.cutouts), + total=len(self.cutouts), + desc="Processing Cutouts", + disable=not self.output, + ): + + print( + "Computing for Cutout number :{}/{}".format( + i + 1, len(self.cutouts) + ) + ) + + catalogue = homology.compute_ph_components( + cutout, + self.local_bg, + analysis_threshold_val=self.analysis_threshold_val, + lifetime_limit=lifetime_limit, + output=self.output, + bg_map=self.bg_map, + area_limit=self.area_limit, + GPU=self.GPU, + lifetime_limit_fraction=lifetime_limit_fraction, + mean_bg=self.mean_bg, + IDoffset=IDoffset, + box_size=self.box_size, + detection_threshold=self.sigma, + Cutout_X_offset=self.coords[i][1], + Cutout_Y_offset=self.coords[i][0], + ) if len(catalogue) == 0: continue - + IDoffset += len(catalogue) - - catalogue['Y0_cutout'] = self.coords[i][0] - catalogue['X0_cutout'] = self.coords[i][1] - # corrent the x1 position for the cutout - catalogue['x1'] = catalogue['x1'] + self.coords[i][0] - catalogue['x2'] = catalogue['x2'] + self.coords[i][0] - catalogue['y1'] = catalogue['y1'] + self.coords[i][1] - catalogue['y2'] = catalogue['y2'] + self.coords[i][1] - catalogue['bbox1'] = catalogue['bbox1'] + self.coords[i][0] - catalogue['bbox2'] = catalogue['bbox2'] + self.coords[i][1] - catalogue['bbox3'] = catalogue['bbox3'] + self.coords[i][0] - catalogue['bbox4'] = catalogue['bbox4'] + self.coords[i][1] - catalogue['distance_from_center'] = ((catalogue['x1']-cutout.shape[0]/2)**2 + (catalogue['y1']-cutout.shape[1]/2)**2)**0.5 - catalogue['cutup_number'] = i + + catalogue["Y0_cutout"] = self.coords[i][0] + catalogue["X0_cutout"] = self.coords[i][1] + catalogue["x1"] = catalogue["x1"] + self.coords[i][0] + catalogue["x2"] = catalogue["x2"] + self.coords[i][0] + catalogue["y1"] = catalogue["y1"] + self.coords[i][1] + catalogue["y2"] = catalogue["y2"] + self.coords[i][1] + catalogue["bbox1"] = catalogue["bbox1"] + self.coords[i][0] + catalogue["bbox2"] = catalogue["bbox2"] + self.coords[i][1] + catalogue["bbox3"] = catalogue["bbox3"] + self.coords[i][0] + catalogue["bbox4"] = catalogue["bbox4"] + self.coords[i][1] + catalogue["distance_from_center"] = ( + (catalogue["x1"] - cutout.shape[0] / 2) ** 2 + + (catalogue["y1"] - cutout.shape[1] / 2) ** 2 + ) ** 0.5 + catalogue["cutup_number"] = i catalogue_list.append(catalogue) self.catalogue = pd.concat(catalogue_list) + # print(self.catalogue) # remove duplicates and keep the one closest to its cutout centre. - #print('before duplicated removal :',len(self.catalogue)) - #self.catalogue = utils.remove_duplicates(self.catalogue) + # print('before duplicated removal :',len(self.catalogue)) + # self.catalogue = utils.remove_duplicates(self.catalogue) # drop any with edge_flag == 1 # set edge flag False to 0 - self.catalogue['edge_flag'] = self.catalogue['edge_flag'].astype(int) - + self.catalogue["edge_flag"] = self.catalogue["edge_flag"].astype(int) + if self.remove_edge: self.catalogue = self.catalogue[self.catalogue.edge_flag != 1] - self.catalogue = self.catalogue.sort_values(by=['distance_from_center'],ascending=True) - self.catalogue = self.catalogue.drop_duplicates(subset=['x1','y1','Birth'], keep='first') + self.catalogue = self.catalogue.sort_values( + by=["distance_from_center"], ascending=True + ) + self.catalogue = self.catalogue.drop_duplicates( + subset=["x1", "y1", "Birth"], keep="first" + ) else: - + for i, row in catalogue.iterrows(): if row.edge_flag == 1: if row.bbox1 == 0: @@ -273,27 +292,38 @@ def phsf(self, lifetime_limit : float = 0,lifetime_limit_fraction : float = 2): if row.bbox2 == 0: row.bbox2 = 1 if row.bbox3 == self.image.shape[0]: - row.bbox3 = self.image.shape[0]-1 + row.bbox3 = self.image.shape[0] - 1 if row.bbox4 == self.image.shape[1]: - row.bbox4 = self.image.shape[1]-1 - - self.catalogue.at[i,'bbox1'] = row.bbox1 - self.catalogue.at[i,'bbox2'] = row.bbox2 - self.catalogue.at[i,'bbox3'] = row.bbox3 - self.catalogue.at[i,'bbox4'] = row.bbox4 - + row.bbox4 = self.image.shape[1] - 1 + + self.catalogue.at[i, "bbox1"] = row.bbox1 + self.catalogue.at[i, "bbox2"] = row.bbox2 + self.catalogue.at[i, "bbox3"] = row.bbox3 + self.catalogue.at[i, "bbox4"] = row.bbox4 + else: IDoffset = 0 - catalogue = homology.compute_ph_components(self.image,self.local_bg,analysis_threshold_val=self.analysis_threshold_val, - lifetime_limit=lifetime_limit,output=self.output,bg_map=self.bg_map,area_limit=self.area_limit, - GPU=self.GPU,lifetime_limit_fraction=lifetime_limit_fraction,mean_bg=self.mean_bg, - IDoffset=IDoffset,box_size=self.cutup_size,detection_threshold=self.sigma) + catalogue = homology.compute_ph_components( + self.image, + self.local_bg, + analysis_threshold_val=self.analysis_threshold_val, + lifetime_limit=lifetime_limit, + output=self.output, + bg_map=self.bg_map, + area_limit=self.area_limit, + GPU=self.GPU, + lifetime_limit_fraction=lifetime_limit_fraction, + mean_bg=self.mean_bg, + IDoffset=IDoffset, + box_size=self.cutup_size, + detection_threshold=self.sigma, + ) self.catalogue = catalogue - - self.catalogue['Y0_cutout'] = 0 - self.catalogue['X0_cutout'] = 0 - self.catalogue['edge_flag'] = self.catalogue['edge_flag'].astype(int) - + + self.catalogue["Y0_cutout"] = 0 + self.catalogue["X0_cutout"] = 0 + self.catalogue["edge_flag"] = self.catalogue["edge_flag"].astype(int) + if self.remove_edge: self.catalogue = self.catalogue[self.catalogue.edge_flag != 1] else: @@ -304,21 +334,20 @@ def phsf(self, lifetime_limit : float = 0,lifetime_limit_fraction : float = 2): if row.bbox2 == 0: row.bbox2 = 1 if row.bbox3 == self.image.shape[0]: - row.bbox3 = self.image.shape[0]-1 + row.bbox3 = self.image.shape[0] - 1 if row.bbox4 == self.image.shape[1]: - row.bbox4 = self.image.shape[1]-1 - - self.catalogue.at[i,'bbox1'] = row.bbox1 - self.catalogue.at[i,'bbox2'] = row.bbox2 - self.catalogue.at[i,'bbox3'] = row.bbox3 - self.catalogue.at[i,'bbox4'] = row.bbox4 - - - #print(self.catalogue) - - self.catalogue = self.catalogue.sort_values(by=['lifetime'],ascending=False) - #print('after duplicate removal :',len(self.catalogue)) - + row.bbox4 = self.image.shape[1] - 1 + + self.catalogue.at[i, "bbox1"] = row.bbox1 + self.catalogue.at[i, "bbox2"] = row.bbox2 + self.catalogue.at[i, "bbox3"] = row.bbox3 + self.catalogue.at[i, "bbox4"] = row.bbox4 + + # print(self.catalogue) + + self.catalogue = self.catalogue.sort_values(by=["lifetime"], ascending=False) + # print('after duplicate removal :',len(self.catalogue)) + # do enclosed_i evaluation with the bounding box to ensure we dont use the whole image. # plt.figure(figsize=(20,20)) # plt.imshow(self.image,cmap='gray_r',norm=colors.LogNorm(clip=True)) @@ -334,100 +363,118 @@ def phsf(self, lifetime_limit : float = 0,lifetime_limit_fraction : float = 2): # time.sleep(3) enclosed_i_list = [] t0 = time.time() - for i in tqdm(range(0,len(self.catalogue)),total=len(self.catalogue),desc='Calculating enclosed_i',disable=not self.output): + for i in tqdm( + range(0, len(self.catalogue)), + total=len(self.catalogue), + desc="Calculating enclosed_i", + disable=not self.output, + ): row = self.catalogue.iloc[i] x1 = row.x1 - row.bbox1 + 1 y1 = row.y1 - row.bbox2 + 1 Birth = row.Birth Death = row.Death # is this a new row? - #if row.new_row == 1: - - img = self.image[int(row.bbox1)-1:int(row.bbox3)+1,int(row.bbox2)-1:int(row.bbox4)+1] + # if row.new_row == 1: + + img = self.image[ + int(row.bbox1) - 1 : int(row.bbox3) + 1, + int(row.bbox2) - 1 : int(row.bbox4) + 1, + ] # reduce the cat to just the sources in the bounding box. - - cat = self.catalogue[self.catalogue['x1'] > int(row.bbox1)] - cat = cat[cat['x1'] < int(row.bbox3)] - cat = cat[cat['y1'] > int(row.bbox2)] - cat = cat[cat['y1'] < int(row.bbox4)] - cat['x1'] = cat['x1'] - int(row.bbox1) + 1 - cat['y1'] = cat['y1'] - int(row.bbox2) + 1 - - + + cat = self.catalogue[self.catalogue["x1"] > int(row.bbox1)] + cat = cat[cat["x1"] < int(row.bbox3)] + cat = cat[cat["y1"] > int(row.bbox2)] + cat = cat[cat["y1"] < int(row.bbox4)] + cat["x1"] = cat["x1"] - int(row.bbox1) + 1 + cat["y1"] = cat["y1"] - int(row.bbox2) + 1 + # plt.imshow(img,cmap='gray_r',norm=colors.LogNorm(clip=True,vmin=1E-13,vmax=1E-9)) # plt.scatter(cat.y1,cat.x1,c='r',marker='x') # plt.savefig('test.png') - # plt.close() + # plt.close() # # sleep # time.sleep(2) - - if self.GPU==True: + + if self.GPU == True: import cupy as cp + # this is not the best way to deal with this. We should crop the gpu version of the image. img_gpu = cp.asarray(img, dtype=cp.float64) - enclosed_i = homology.make_point_enclosure_assoc_GPU(0,x1,y1,Birth,Death,cat,img_gpu) + enclosed_i = homology.make_point_enclosure_assoc_GPU( + 0, x1, y1, Birth, Death, cat, img_gpu + ) enclosed_i_list.append(enclosed_i) else: - #print(self.catalogue) - enclosed_i = homology.make_point_enclosure_assoc_CPU(0,x1,y1,Birth,Death,cat,img) + # print(self.catalogue) + enclosed_i = homology.make_point_enclosure_assoc_CPU( + 0, x1, y1, Birth, Death, cat, img + ) enclosed_i_list.append(enclosed_i) - - #print('enclosed_i calculated! t='+str(time.time()-t0)+' s') - self.catalogue['enclosed_i'] = enclosed_i_list - + + # print('enclosed_i calculated! t='+str(time.time()-t0)+' s') + self.catalogue["enclosed_i"] = enclosed_i_list + # print(self.catalogue) + # print("Enclosed_i calculated! TESTESTEST") # correct for first destruction - - #print(len(self.catalogue)) + # t0_parent_tag = time.time() + # self.catalogue = homology.parent_tag_func_vectorized_new(self.catalogue) + # t1_parent_tag = time.time() + # print("Time to assign parent tags: ", t1_parent_tag - t0_parent_tag) + # print(self.catalogue) + # print(len(self.catalogue)) t0_correct_firs = time.time() - #print('BEfore',len(self.catalogue)) - self.catalogue = homology.correct_first_destruction(self.catalogue,output=not self.output) + # print('BEfore',len(self.catalogue)) + self.catalogue = homology.correct_first_destruction( + self.catalogue, output=not self.output + ) t1_correct_firs = time.time() - #print('Time to correct first destruction: ',t1_correct_firs-t0_correct_firs) - #print('After',len(self.catalogue)) + print("Time to correct first destruction: ", t1_correct_firs - t0_correct_firs) + # print('Time to correct first destruction: ',t1_correct_firs-t0_correct_firs) + # print('After',len(self.catalogue)) # parent tag - #print("Assigning parent tags..") + # print("Assigning parent tags..") t0_parent_tag = time.time() - self.catalogue = homology.parent_tag_func_vectorized_new(self.catalogue) + self.catalogue = homology.parent_tag_func_optimized(self.catalogue) t1_parent_tag = time.time() - #print('Time to assign parent tags: ',t1_parent_tag-t0_parent_tag) - #print(self.catalogue['parent_tag'].value_counts()) - #print("Classifying sources in hirearchy..") + + print("Time to assign parent tags: ", t1_parent_tag - t0_parent_tag) + t0_classify = time.time() - self.catalogue['Class'] = self.catalogue.apply(homology.classify_single,axis=1) + self.catalogue["Class"] = self.catalogue.apply(homology.classify_single, axis=1) t1_classify = time.time() - #print('Time to classify sources: ',t1_classify-t0_classify) - # put ID at the front - - - - - - - def set_background(self,detection_threshold,analysis_threshold, - set_bg=None,bg_map_bool=False,box_size=None,mode='mad_std'): - - #print('Setting background..') - - + print("Time to classify sources: ", t1_classify - t0_classify) + + def set_background( + self, + detection_threshold, + analysis_threshold, + set_bg=None, + bg_map_bool=False, + box_size=None, + mode="mad_std", + ): + self.sigma = detection_threshold self.analysis_threshold = analysis_threshold self.bg_map = bg_map_bool # mode should be MAD_Std, RMS or other. - - if mode == 'Radio': - # old verion was called radio. - mode = 'mad_std' + + if mode == "Radio": + # old verion was called radio. + mode = "mad_std" # need to account dor the cutputs if not usinh bg_map. - + # bg_map and cutup are require only the same code. - + # bg_map and no cuput is the same as bg_map and cutup. - + # no cutup and no bg_map is just one estimation for the whole image. - + if self.cutup == True: - + # we want to do the bg_map but for the cutout dims. as the box size is in pixels. self.bg_map = True bg_map_bool = True @@ -439,43 +486,45 @@ def set_background(self,detection_threshold,analysis_threshold, box_size = self.cutup_size else: pass - + self.box_size = box_size + print(self.box_size) if bg_map_bool == True: - #print('Creating a background map. Inputed Box size = ',box_size) + # print('Creating a background map. Inputed Box size = ',box_size) # these will be returned as arrays like a map. - std, mean_bg = background.calculate_background_map(self.image,box_size,mode=mode) - #print('Background map created.') - #print('Mean Background across cutouts: ', np.nanmean(std)) - #print('Median of bg distribution: ', np.nanmean(mean_bg)) - + std, mean_bg = background.calculate_background_map( + self.image, box_size, mode=mode + ) + # print('Background map created.') + # print('Mean Background across cutouts: ', np.nanmean(std)) + # print('Median of bg distribution: ', np.nanmean(mean_bg)) + else: - #print('Not creating a background map.') - std, mean_bg = background.calculate_background(self.image,mode=mode) - #print('Background set to: ',std) - #print('Background mean set to: ',mean_bg) - + # print('Not creating a background map.') + std, mean_bg = background.calculate_background(self.image, mode=mode) + # print('Background set to: ',std) + # print('Background mean set to: ',mean_bg) + if set_bg is not None: # set bg should be a tuple of (std,mean_bg) - #print('User has set the background.') + # print('User has set the background.') std = set_bg[0] mean_bg = set_bg[1] - #print('Background set to: ',std) - #print('Background mean set to: ',mean_bg) - - - self.local_bg = std*self.sigma - self.analysis_threshold_val = std*self.analysis_threshold - self.mean_bg = mean_bg - - - + self.local_bg = std * self.sigma + self.analysis_threshold_val = std * self.analysis_threshold + self.mean_bg = mean_bg - def set_background_old(self,detection_threshold : float,analysis_threshold, - set_bg : float = None, bg_map : bool = None, - box_size : int = 10, mode : str = 'Radio'): - """Sets the background for the source finding algorithm. + def set_background_old( + self, + detection_threshold: float, + analysis_threshold, + set_bg: float = None, + bg_map: bool = None, + box_size: int = 10, + mode: str = "Radio", + ): + """Sets the background for the source finding algorithm. Args: detection_threshold (int): _description_ @@ -485,27 +534,30 @@ def set_background_old(self,detection_threshold : float,analysis_threshold, box_size (int, optional): _description_. Defaults to 10. mode (str, optional): _description_. Defaults to 'Radio'. """ - - + self.bg_map = bg_map self.sigma = detection_threshold self.analysis_threshold = analysis_threshold - - if mode == 'Radio': + + if mode == "Radio": if self.cutup: - + # loop though each cutout and calculate the local background. - + if bg_map is not None: - + # users wants to use background map so lets make it local_bg_list = [] analysis_threshold_list = [] mean_bg_list = [] for i, cutout in enumerate(self.cutouts): - local_bg_map, mean_bg = background.radio_background_map(cutout, box_size) - analysis_threshold_list.append(local_bg_map*self.analysis_threshold) - local_bg_list.append(local_bg_map*self.sigma) + local_bg_map, mean_bg = background.radio_background_map( + cutout, box_size + ) + analysis_threshold_list.append( + local_bg_map * self.analysis_threshold + ) + local_bg_list.append(local_bg_map * self.sigma) mean_bg_list.append(mean_bg) else: local_bg_list = [] @@ -513,110 +565,112 @@ def set_background_old(self,detection_threshold : float,analysis_threshold, mean_bg_list = [] for cutout in self.cutouts: local_bg, mean_bg = background.radio_background(cutout) - analysis_threshold_list.append(local_bg*self.analysis_threshold) - local_bg_list.append(local_bg*self.sigma) + analysis_threshold_list.append( + local_bg * self.analysis_threshold + ) + local_bg_list.append(local_bg * self.sigma) mean_bg_list.append(mean_bg) - local_bg = local_bg_list + local_bg = local_bg_list analysis_threshold = analysis_threshold_list mean_bg = mean_bg_list - + else: # Radio background is calculated using the median absolute deviation of the total image. if bg_map is not None: - local_bg_o, mean_bg = background.radio_background_map(self.image, box_size) - local_bg = local_bg_o*self.sigma - analysis_threshold = local_bg_o*self.analysis_threshold + local_bg_o, mean_bg = background.radio_background_map( + self.image, box_size + ) + local_bg = local_bg_o * self.sigma + analysis_threshold = local_bg_o * self.analysis_threshold else: local_bg_o, mean_bg = background.radio_background(self.image) - local_bg = local_bg_o*self.sigma - analysis_threshold = local_bg_o*self.analysis_threshold - - - if mode == 'Optical': + local_bg = local_bg_o * self.sigma + analysis_threshold = local_bg_o * self.analysis_threshold + + if mode == "Optical": # Optical background is calculated using a random sample of pixels mean_bg, std_bg = background.optical_background(nsamples=1000) - local_bg = mean_bg + self.sigma*std_bg - analysis_threshold = mean_bg + std_bg*self.analysis_threshold + local_bg = mean_bg + self.sigma * std_bg + analysis_threshold = mean_bg + std_bg * self.analysis_threshold - if mode == 'other': + if mode == "other": # If the user has a custom background function, they can pass it in here. - local_bg = set_bg*self.sigma - analysis_threshold = local_bg*self.analysis_threshold - #print('Background set to: ',local_bg) - #print('Analysis threshold set to: ',analysis_threshold) - + local_bg = set_bg * self.sigma + analysis_threshold = local_bg * self.analysis_threshold + # print('Background set to: ',local_bg) + # print('Analysis threshold set to: ',analysis_threshold) + self.analysis_threshold_val = analysis_threshold self.local_bg = local_bg self.mean_bg = mean_bg - #print(self.mean_bg) - - #if bg_map: - #print('Using bg_map for analysis.') - #else: - #if self.cutup: - - # print('Mean Background across cutouts: ', np.nanmean(self.local_bg)) - # print('Median of bg distribution: ', np.nanmean(self.mean_bg)) - - # else: - # print('Background set to: ',self.local_bg) + # print(self.mean_bg) + # if bg_map: + # print('Using bg_map for analysis.') + # else: + # if self.cutup: + # print('Mean Background across cutouts: ', np.nanmean(self.local_bg)) + # print('Median of bg distribution: ', np.nanmean(self.mean_bg)) + # else: + # print('Background set to: ',self.local_bg) - - - - - - - - - def source_characterising(self, use_gpu : bool = False): + def source_characterising(self, use_gpu: bool = False): """Source Characterising function. This function takes the catalogue and the image and calculates the source properties. Args: - use_gpu (bool, optional): Option to use the GPU True to use and False to not, + use_gpu (bool, optional): Option to use the GPU True to use and False to not, requires cupy module and a avalible GPU. Defaults to False. """ - - self.catalogue, self.polygons = source.measure_source_properties(use_gpu=use_gpu,catalogue=self.catalogue, - cutout=self.image,background_map=self.local_bg, - output=self.output,cutupts=self.cutouts,mode=self.mode,header=self.header) - + + self.catalogue, self.polygons = source.measure_source_properties( + use_gpu=use_gpu, + catalogue=self.catalogue, + cutout=self.image, + background_map=self.local_bg, + output=self.output, + cutupts=self.cutouts, + mode=self.mode, + header=self.header, + sigma=self.sigma, + ) + if self.Xoff is not None: # correct for the poistion of the cutout. when using cutout from a larger image. - self.catalogue['Xc'] = self.catalogue['Xc'] + self.Xoff - self.catalogue['bbox1'] = self.catalogue['bbox1'] + self.Xoff - self.catalogue['bbox3'] = self.catalogue['bbox3'] + self.Xoff - + self.catalogue["Xc"] = self.catalogue["Xc"] + self.Xoff + self.catalogue["bbox1"] = self.catalogue["bbox1"] + self.Xoff + self.catalogue["bbox3"] = self.catalogue["bbox3"] + self.Xoff + if self.Yoff is not None: # correct for the poistion of the cutout. when using cutout from a larger image. - self.catalogue['Yc'] = self.catalogue['Yc'] + self.Yoff - self.catalogue['bbox2'] = self.catalogue['bbox2'] + self.Yoff - self.catalogue['bbox4'] = self.catalogue['bbox4'] + self.Yoff + self.catalogue["Yc"] = self.catalogue["Yc"] + self.Yoff + self.catalogue["bbox2"] = self.catalogue["bbox2"] + self.Yoff + self.catalogue["bbox4"] = self.catalogue["bbox4"] + self.Yoff if self.header is not None: - #try: - #print('Converting Xc and Yc to RA and DEC') - Ra, Dec = utils.xy_to_RaDec(self.catalogue['Xc'],self.catalogue['Yc'],self.header,mode=self.mode) - self.catalogue['RA'] = Ra - self.catalogue['DEC'] = Dec - self.catalogue['RA'] = self.catalogue['RA'].astype(float) - self.catalogue['DEC'] = self.catalogue['DEC'].astype(float) - #except: + # try: + # print('Converting Xc and Yc to RA and DEC') + Ra, Dec = utils.xy_to_RaDec( + self.catalogue["Xc"], self.catalogue["Yc"], self.header, mode=self.mode + ) + self.catalogue["RA"] = Ra + self.catalogue["DEC"] = Dec + self.catalogue["RA"] = self.catalogue["RA"].astype(float) + self.catalogue["DEC"] = self.catalogue["DEC"].astype(float) + # except: # pass - + self._set_types_of_dataframe() - - if self.mode == 'optical': - + + if self.mode == "optical": + def ABmag(flux): - return -2.5*np.log10(flux) - - def RONoise(EFFRON,EFFGAIN,EXPTIME,Area): - return np.sqrt(Area)*(EFFRON/EFFGAIN)*EXPTIME + return -2.5 * np.log10(flux) + + def RONoise(EFFRON, EFFGAIN, EXPTIME, Area): + return np.sqrt(Area) * (EFFRON / EFFGAIN) * EXPTIME def SkyNoise(sky): return np.sqrt(sky) @@ -624,127 +678,140 @@ def SkyNoise(sky): def SourceNoise(Flux): return np.sqrt(Flux) - def Flux_err(EFFRON,EFFGAIN,EXPTIME,Area,sky,Flux): - return np.sqrt(RONoise(EFFRON,EFFGAIN,EXPTIME,Area)**2 + SkyNoise(sky) + SourceNoise(Flux)) + def Flux_err(EFFRON, EFFGAIN, EXPTIME, Area, sky, Flux): + return np.sqrt( + RONoise(EFFRON, EFFGAIN, EXPTIME, Area) ** 2 + + SkyNoise(sky) + + SourceNoise(Flux) + ) + + def NOISE(row, local_ng): + return np.sum( + np.random.normal(row["mean_bg"], local_ng, int(row["Area"])) + ) - def NOISE(row,local_ng): - return np.sum(np.random.normal(row['mean_bg'],local_ng,int(row['Area']))) - EFFGAIN = utils.get_EFFGAIN(self.header) EXPTIME = utils.get_EXPTIME(self.header) EFFRON = utils.get_EFFRON(self.header) - #print('EFFGAIN: ',EFFGAIN) - #print('EXPTIME: ',EXPTIME) - #print('EFFRON: ',EFFRON) - #print(self.catalogue['Noise']) - self.catalogue['Flux_total_new'] = (self.catalogue['Flux_total']*EFFGAIN*EXPTIME - self.catalogue['Noise']*EFFGAIN*EXPTIME - RONoise(EFFRON,EFFGAIN,EXPTIME,self.catalogue['Area'])) - #print('Flux_total_new: ',self.catalogue['Flux_total_new']) - self.catalogue['Flux_total_err'] = Flux_err(EFFRON,EFFGAIN,EXPTIME,self.catalogue['Area'],self.catalogue['Noise']*EFFGAIN*EXPTIME,self.catalogue['Flux_total_new']) - #print('Flux_total_err: ',self.catalogue['Flux_total_err']) - self.catalogue['SNR'] = self.catalogue['Flux_total_new']/self.catalogue['Flux_total_err']/self.catalogue['Area'] - #print('SNR: ',self.catalogue['SNR']) - self.catalogue['MAG_err'] = 1/self.catalogue['SNR'] # use the approximate error for the magnitude. - self.catalogue['Flux_total_new'] = self.catalogue['Flux_total_new']/(EFFGAIN*EXPTIME) - self.catalogue['MAG_flux'] = ABmag(self.catalogue['Flux_total_new']) - - - - - - - def create_polygons(self,use_gpu=False): - ''' + # print('EFFGAIN: ',EFFGAIN) + # print('EXPTIME: ',EXPTIME) + # print('EFFRON: ',EFFRON) + # print(self.catalogue['Noise']) + self.catalogue["Flux_total_new"] = ( + self.catalogue["Flux_total"] * EFFGAIN * EXPTIME + - self.catalogue["Noise"] * EFFGAIN * EXPTIME + - RONoise(EFFRON, EFFGAIN, EXPTIME, self.catalogue["Area"]) + ) + # print('Flux_total_new: ',self.catalogue['Flux_total_new']) + self.catalogue["Flux_total_err"] = Flux_err( + EFFRON, + EFFGAIN, + EXPTIME, + self.catalogue["Area"], + self.catalogue["Noise"] * EFFGAIN * EXPTIME, + self.catalogue["Flux_total_new"], + ) + # print('Flux_total_err: ',self.catalogue['Flux_total_err']) + self.catalogue["SNR"] = ( + self.catalogue["Flux_total_new"] + / self.catalogue["Flux_total_err"] + / self.catalogue["Area"] + ) + # print('SNR: ',self.catalogue['SNR']) + self.catalogue["MAG_err"] = ( + 1 / self.catalogue["SNR"] + ) # use the approximate error for the magnitude. + self.catalogue["Flux_total_new"] = self.catalogue["Flux_total_new"] / ( + EFFGAIN * EXPTIME + ) + self.catalogue["MAG_flux"] = ABmag(self.catalogue["Flux_total_new"]) + + def create_polygons(self, use_gpu=False): + """ Creates Polygons/contours best when you just want segmentations and not source charateristics. - ''' - - self.catalogue = source.create_polygons(use_gpu=use_gpu,catalogue=self.catalogue, - cutout=self.image,output=self.output,cutupts=self.cutouts) - - - - - - + """ + self.catalogue = source.create_polygons( + use_gpu=use_gpu, + catalogue=self.catalogue, + cutout=self.image, + output=self.output, + cutupts=self.cutouts, + ) def create_polygons_fast(self): - + # since we have a bounding box, we can just create a polygon in the bounding box. - + polygons = [] - for index, row in tqdm(self.catalogue.iterrows(),total=len(self.catalogue),desc='Creating polygons'): - contour = utils._get_polygons_in_bbox(row.bbox2-2,row.bbox4+2,row.bbox1-2,row.bbox3+2,row.x1,row.y1,row.Birth,row.Death) + for index, row in tqdm( + self.catalogue.iterrows(), + total=len(self.catalogue), + desc="Creating polygons", + ): + contour = utils._get_polygons_in_bbox( + row.bbox2 - 2, + row.bbox4 + 2, + row.bbox1 - 2, + row.bbox3 + 2, + row.x1, + row.y1, + row.Birth, + row.Death, + ) polygons.append(contour) self.polygons = polygons - - - - - - def create_polygons_gpu(self): - ''' + """ Recommended when using GPU acceleration. and not using the bounding box to simplify the polygon creation. - ''' + """ polygons = [] self.image_gpu = cp.asarray(self.image, dtype=cp.float64) - + for index, row in self.catalogue.iterrows(): t0 = time.time() - contour = utils._get_polygons_gpu(row.x1,row.y1,row.Birth,row.Death) + contour = utils._get_polygons_gpu(row.x1, row.y1, row.Birth, row.Death) t1 = time.time() - #print('Time to create polygon: ',t1-t0) + # print('Time to create polygon: ',t1-t0) polygons.append(contour) self.polygons = polygons - - - - - - def _set_types_of_dataframe(self): """ - Sets the catalogue to the correct data types. This is important to allow for writing data. + Sets the catalogue to the correct data types. This is important to allow for writing data. Otherwise the datatypes will remain object which will try to be pickled. - + """ - self.catalogue['ID'] = self.catalogue['ID'].astype(int) - self.catalogue['Birth'] = self.catalogue['Birth'].astype(float) - self.catalogue['Death'] = self.catalogue['Death'].astype(float) - self.catalogue['x1'] = self.catalogue['x1'].astype(float) - self.catalogue['y1'] = self.catalogue['y1'].astype(float) - self.catalogue['x2'] = self.catalogue['x2'].astype(float) - self.catalogue['y2'] = self.catalogue['y2'].astype(float) - self.catalogue['Flux_total'] = self.catalogue['Flux_total'].astype(float) - self.catalogue['Flux_peak'] = self.catalogue['Flux_peak'].astype(float) - self.catalogue['Area'] = self.catalogue['Area'].astype(float) - self.catalogue['Xc'] = self.catalogue['Xc'].astype(float) - self.catalogue['Yc'] = self.catalogue['Yc'].astype(float) - self.catalogue['bbox1'] = self.catalogue['bbox1'].astype(float) - self.catalogue['bbox2'] = self.catalogue['bbox2'].astype(float) - self.catalogue['bbox3'] = self.catalogue['bbox3'].astype(float) - self.catalogue['bbox4'] = self.catalogue['bbox4'].astype(float) - self.catalogue['Maj'] = self.catalogue['Maj'].astype(float) - self.catalogue['Min'] = self.catalogue['Min'].astype(float) - self.catalogue['Pa'] = self.catalogue['Pa'].astype(float) - self.catalogue['parent_tag'] = self.catalogue['parent_tag'].astype(float) - self.catalogue['Class'] = self.catalogue['Class'].astype(float) - if self.cutup: - self.catalogue['Y0_cutout'] = self.catalogue['Y0_cutout'].astype(float) - self.catalogue['X0_cutout'] = self.catalogue['X0_cutout'].astype(float) - self.catalogue['SNR'] = self.catalogue['SNR'].astype(float) - self.catalogue['Noise'] = self.catalogue['Noise'].astype(float) - - - - - - - - def plot_sources(self,cmap,figsize=(10,10),norm='linear',save_path=None): + self.catalogue["ID"] = self.catalogue["ID"].astype(int) + self.catalogue["Birth"] = self.catalogue["Birth"].astype(float) + self.catalogue["Death"] = self.catalogue["Death"].astype(float) + self.catalogue["x1"] = self.catalogue["x1"].astype(float) + self.catalogue["y1"] = self.catalogue["y1"].astype(float) + self.catalogue["x2"] = self.catalogue["x2"].astype(float) + self.catalogue["y2"] = self.catalogue["y2"].astype(float) + self.catalogue["Flux_total"] = self.catalogue["Flux_total"].astype(float) + self.catalogue["Flux_peak"] = self.catalogue["Flux_peak"].astype(float) + self.catalogue["Area"] = self.catalogue["Area"].astype(float) + self.catalogue["Xc"] = self.catalogue["Xc"].astype(float) + self.catalogue["Yc"] = self.catalogue["Yc"].astype(float) + self.catalogue["bbox1"] = self.catalogue["bbox1"].astype(float) + self.catalogue["bbox2"] = self.catalogue["bbox2"].astype(float) + self.catalogue["bbox3"] = self.catalogue["bbox3"].astype(float) + self.catalogue["bbox4"] = self.catalogue["bbox4"].astype(float) + self.catalogue["Maj"] = self.catalogue["Maj"].astype(float) + self.catalogue["Min"] = self.catalogue["Min"].astype(float) + self.catalogue["Pa"] = self.catalogue["Pa"].astype(float) + self.catalogue["parent_tag"] = self.catalogue["parent_tag"].astype(float) + self.catalogue["Class"] = self.catalogue["Class"].astype(float) + if self.cutup: + self.catalogue["Y0_cutout"] = self.catalogue["Y0_cutout"].astype(float) + self.catalogue["X0_cutout"] = self.catalogue["X0_cutout"].astype(float) + self.catalogue["SNR"] = self.catalogue["SNR"].astype(float) + self.catalogue["Noise"] = self.catalogue["Noise"].astype(float) + + def plot_sources(self, cmap, figsize=(10, 10), norm="linear", save_path=None): """Plots the source polygons on the image. Args: @@ -752,111 +819,104 @@ def plot_sources(self,cmap,figsize=(10,10),norm='linear',save_path=None): figsize (tuple, optional): Desired figure size. Defaults to (10,10). norm (str, optional): _description_. Defaults to 'linear'. save_path (str, optional): Save path if you desire to save the figure. Defaults to None. - + """ plt.figure(figsize=figsize) - plt.imshow(self.image,cmap=cmap,origin='lower',norm=norm) - #plt.scatter(self.catalogue['Xc'],self.catalogue['Yc'],s=10,c='r') - + plt.imshow(self.image, cmap=cmap, origin="lower", norm=norm) + # plt.scatter(self.catalogue['Xc'],self.catalogue['Yc'],s=10,c='r') + for i, poly in enumerate(self.polygons): if poly is not None: - plt.plot(poly[:,1],poly[:,0]) + plt.plot(poly[:, 1], poly[:, 0]) if save_path is not None: plt.savefig(save_path) plt.show() - - - def save_catalogue(self,save_path,filetype=None,overwrite=False): + def save_catalogue(self, save_path, filetype=None, overwrite=False): """Save Catalogue to a file. Args: save_path (str): Desired path to save the catalogue. filetype (str, optional): Specify the file type or include the approprate file extention. Defaults to None. overwrite (bool, optional): Overwrite the save file if the name is the same. Defaults to False. - + """ - + # get the extension from the save_path - #print('Saving Catalogue to file: ',save_path) - fileextention = save_path.split('.')[-1] - + # print('Saving Catalogue to file: ',save_path) + fileextention = save_path.split(".")[-1] + if filetype is None: filetype = fileextention - - if filetype == 'csv': - self.catalogue.to_csv(save_path,index=False,overwrite=overwrite) - # print('Catalogue saved to: ',save_path) - - if filetype == 'fits': + + if filetype == "csv": + self.catalogue.to_csv(save_path, index=False, overwrite=overwrite) + # print('Catalogue saved to: ',save_path) + + if filetype == "fits": from astropy.table import Table - # print('Saving to fits with astropy') - enclosed_i = self.catalogue['enclosed_i'] - + + # print('Saving to fits with astropy') + enclosed_i = self.catalogue["enclosed_i"] + for i in range(len(enclosed_i)): for j in range(len(enclosed_i[i])): enclosed_i[i][j] = int(enclosed_i[i][j]) if len(enclosed_i[i]) == 0: enclosed_i[i] = [0] - self.catalogue['enclosed_i'] = enclosed_i - #print(self.catalogue) + self.catalogue["enclosed_i"] = enclosed_i + # print(self.catalogue) t = Table.from_pandas(self.catalogue) - t.write(save_path,overwrite=overwrite) - - - if filetype == 'hdf': - self.catalogue.to_hdf(save_path,key='catalogue',mode='w') - # print('Catalogue saved to: ',save_path) - - if filetype == ('txt' or 'ascii'): - self.catalogue.to_csv(save_path,index=False,overwrite=overwrite) - #print('Catalogue saved to: ',save_path) - - def open_catalogue(self,file_path,filetype=None): + t.write(save_path, overwrite=overwrite) + + if filetype == "hdf": + self.catalogue.to_hdf(save_path, key="catalogue", mode="w") + # print('Catalogue saved to: ',save_path) + + if filetype == ("txt" or "ascii"): + self.catalogue.to_csv(save_path, index=False, overwrite=overwrite) + # print('Catalogue saved to: ',save_path) + + def open_catalogue(self, file_path, filetype=None): from astropy.table import Table - + self.catalogue = Table.read(file_path) - + for i in range(len(self.catalogue)): - self.catalogue['contour'][i] = np.array(self.catalogue['contour'][i]).reshape(-1,2) - - self.catalogue = self.catalogue.to_pandas() - - - def save_polygons_to_ds9(self, filename): + self.catalogue["contour"][i] = np.array( + self.catalogue["contour"][i] + ).reshape(-1, 2) + + self.catalogue = self.catalogue.to_pandas() - ''' + def save_polygons_to_ds9(self, filename): + """ Saves the polygons to a ds9 region file. - ''' + """ - with open(filename, 'w') as f: - f.write('# Region file format: DS9 version 4.1\n') - f.write('global color=green dashlist=8 3 width=1 font="helvetica 10 normal roman" select=1 highlite=1 dash=0 fixed=0 edit=1 move=1 delete=1 include=1 source=1\n') + with open(filename, "w") as f: + f.write("# Region file format: DS9 version 4.1\n") + f.write( + 'global color=green dashlist=8 3 width=1 font="helvetica 10 normal roman" select=1 highlite=1 dash=0 fixed=0 edit=1 move=1 delete=1 include=1 source=1\n' + ) for polygon in self.polygons: - f.write('polygon(') + f.write("polygon(") for i, point in enumerate(polygon): - f.write('{:.2f},{:.2f}'.format(point[1], point[0])) # note this transformation as the index in some CARTA inmages start at -1. + f.write( + "{:.2f},{:.2f}".format(point[1], point[0]) + ) # note this transformation as the index in some CARTA inmages start at -1. if i < len(polygon) - 1: - f.write(',') - f.write(')\n') - - - - def save_polygons_to_hdf5(self, filename): + f.write(",") + f.write(")\n") - ''' + def save_polygons_to_hdf5(self, filename): + """ Saves the polygons to a hdf5 file. - ''' + """ import h5py - hf = h5py.File(filename, 'w') + hf = h5py.File(filename, "w") for i in range(len(self.catalogue)): - key = self.catalogue['ID'][i] - hf.create_dataset(str(key), data=self.catalogue['contour'][i]) + key = self.catalogue["ID"][i] + hf.create_dataset(str(key), data=self.catalogue["contour"][i]) hf.close() - - - - - - \ No newline at end of file diff --git a/DRUID/src/background.py b/DRUID/src/background.py index cb0df8b..c9d1991 100644 --- a/DRUID/src/background.py +++ b/DRUID/src/background.py @@ -9,7 +9,7 @@ """ import numpy as np -from astropy.stats import mad_std, sigma_clip +from astropy.stats import mad_std, sigma_clipped_stats import matplotlib.pyplot as plt from photutils.background import SExtractorBackground @@ -96,7 +96,17 @@ def calculate_background_map(image,box_size,mode='mad_std'): def get_bg_value_from_result_image(original_location_in_full_image, box_size, bg_map): i, j = original_location_in_full_image - result_value = bg_map[i//box_size + 1, j//box_size + 1] + x = i//box_size + y = j//box_size + # print(bg_map.shape) + # print(x,y) + # acount for image size not being a multiple of box size + if x >= bg_map.shape[0]: + x = bg_map.shape[0] - 1 + if y >= bg_map.shape[1]: + y = bg_map.shape[1] - 1 + + result_value = bg_map[x, y] return result_value @@ -146,10 +156,10 @@ def get_optical_background_estimate(image ,mode): _type_: _description_ """ if mode == 'sigma_clip': - median = np.nanmedian(image) - std = sigma_clip(image).std() - return std, median + mean, median, std = sigma_clipped_stats(image,sigma=3,maxiters=5) + return std, median + if mode == 'SEX': bkg = SExtractorBackground() bkg_sigma = bkg.sigma_clip(image).std() diff --git a/DRUID/src/homology_new.py b/DRUID/src/homology_new.py index d5f9f54..d874f91 100644 --- a/DRUID/src/homology_new.py +++ b/DRUID/src/homology_new.py @@ -1,5 +1,3 @@ - - import cripser import numpy as np from ..src import utils @@ -9,25 +7,27 @@ import pandas from tqdm import tqdm -#from multiprocessing import Pool, freeze_support +# from multiprocessing import Pool, freeze_support # used for debugging import pdb -#from collections import deque +# from collections import deque try: import cupy as cp from cupyx.scipy.ndimage import label as cupy_label + GPU_AVAILABLE = True except: GPU_AVAILABLE = False - + # used for debugging import matplotlib.pyplot as plt -def make_point_enclosure_assoc_GPU(id,x1,y1,Birth,Death,pd,img_gpu): + +def make_point_enclosure_assoc_GPU(id, x1, y1, Birth, Death, pd, img_gpu): """Returns a list of the ID of the points that are enclosed by the mask pd point. Uses GPU for computation. @@ -42,22 +42,24 @@ def make_point_enclosure_assoc_GPU(id,x1,y1,Birth,Death,pd,img_gpu): Returns: enclosed_list (list): _description_ """ - # print(pd) - mask = utils.get_mask_GPU(Birth,Death,x1,y1,img_gpu).get() - #pdb.set_trace() - mask_coords = np.column_stack((pd['x1'], pd['y1'])) - points_inside_mask = mask[mask_coords[:, 0].astype(int), mask_coords[:, 1].astype(int)] - encloses_vectorized = pd.iloc[points_inside_mask]['ID'].tolist() + # print(pd) + mask = utils.get_mask_GPU(Birth, Death, x1, y1, img_gpu).get() + # pdb.set_trace() + mask_coords = np.column_stack((pd["x1"], pd["y1"])) + points_inside_mask = mask[ + mask_coords[:, 0].astype(int), mask_coords[:, 1].astype(int) + ] + encloses_vectorized = pd.iloc[points_inside_mask]["ID"].tolist() # remove self from list - #print(row.ID) - #print(encloses_vectorized) - #encloses_vectorized.remove(row.ID) - #print(encloses_vectorized) - #pdb.set_trace() + # print(row.ID) + # print(encloses_vectorized) + # encloses_vectorized.remove(row.ID) + # print(encloses_vectorized) + # pdb.set_trace() return encloses_vectorized -def make_point_enclosure_assoc_CPU(ID,x1,y1,Birth,Death,pd,img): +def make_point_enclosure_assoc_CPU(ID, x1, y1, Birth, Death, pd, img): """Returns a list of the ID of the points that are enclosed by the mask pd point. Uses GPU for computation. @@ -72,26 +74,27 @@ def make_point_enclosure_assoc_CPU(ID,x1,y1,Birth,Death,pd,img): Returns: enclosed_list (list): _description_ """ - - mask = utils.get_mask_CPU(x1,y1,Birth,Death,img) - #print(mask) - #plt.imshow(mask) - #plt.savefig('test.png') - #pdb.set_trace() - mask_coords = np.column_stack((pd['x1'], pd['y1'])) - points_inside_mask = mask[mask_coords[:, 0].astype(int), mask_coords[:, 1].astype(int)] - encloses_vectorized = pd.iloc[points_inside_mask]['ID'].tolist() + + mask = utils.get_mask_CPU(x1, y1, Birth, Death, img) + # print(mask) + # plt.imshow(mask) + # plt.savefig('test.png') + # pdb.set_trace() + mask_coords = np.column_stack((pd["x1"], pd["y1"])) + points_inside_mask = mask[ + mask_coords[:, 0].astype(int), mask_coords[:, 1].astype(int) + ] + encloses_vectorized = pd.iloc[points_inside_mask]["ID"].tolist() # remove self from list - #print(row.ID) - #print(encloses_vectorized) - #encloses_vectorized.remove(ID) - #print(encloses_vectorized) - #pdb.set_trace() + # print(row.ID) + # print(encloses_vectorized) + # encloses_vectorized.remove(ID) + # print(encloses_vectorized) + # pdb.set_trace() return encloses_vectorized def classify_single(row): - """Classifiying the Rows based on orgin. Args: @@ -102,19 +105,19 @@ def classify_single(row): """ if row.new_row == 0: - if len(row.enclosed_i) <=1: # no children - if row.parent_tag==row.ID: # no parent - return 0 # no children, no parent. + if len(row.enclosed_i) <= 1: # no children + if row.parent_tag == row.ID: # no parent + return 0 # no children, no parent. else: - return 2 # no child has parent. + return 2 # no child has parent. else: - if row.parent_tag==row.ID: # no parent - return 4 # has children, no parent. + if row.parent_tag == row.ID: # no parent + return 4 # has children, no parent. else: - return 3 # has children, has parent. + return 3 # has children, has parent. else: - return 1 # new row has children. - + return 1 # new row has children. + def parent_tag_func_vectorized(df): """ @@ -133,22 +136,23 @@ def parent_tag_func_vectorized(df): """ - enclosed_i_dict = {row['ID']: set(row['enclosed_i']) for idx, row in df.iterrows()} - + enclosed_i_dict = {row["ID"]: set(row["enclosed_i"]) for idx, row in df.iterrows()} + def find_parent_tag(row): - #if row.new_row == 0: # only set it if it is not a new row. - + # if row.new_row == 0: # only set it if it is not a new row. + for ID, enclosed_i_set in enclosed_i_dict.items(): if row.ID in enclosed_i_set: return ID - - return np.nan - - #else: # if it is a new row then we already know the parent tag. + + return np.nan + + # else: # if it is a new row then we already know the parent tag. # return row.parent_tags - + return df.apply(find_parent_tag, axis=1) + def parent_tag_func_vectorized_new(df): """ Vectorised implementation of parent tag function. @@ -161,18 +165,137 @@ def parent_tag_func_vectorized_new(df): """ # Create a dictionary that maps each enclosed ID to its parent ID - parent_tag_dict = {enclosed_id: ID for ID, enclosed_i_set in df[['ID', 'enclosed_i']].values for enclosed_id in enclosed_i_set} + def find_parent(row, df): + potential_parents = [] + for idx, other_row in df.iterrows(): + if row["ID"] in other_row["enclosed_i"] and row["ID"] != other_row["ID"]: + potential_parents.append( + (other_row["ID"], len(other_row["enclosed_i"])) + ) - # Define a function that finds the parent tag of a row - def find_parent_tag(row): - return parent_tag_dict.get(row.ID, np.nan) + if potential_parents: + # Sort by length of enclosed_i (descending) and return the ID of the longest + return sorted(potential_parents, key=lambda x: x[1], reverse=True)[0][0] + + # If no parent found, check if this row has the largest enclosed_i + if len(row["enclosed_i"]) == df["enclosed_i"].apply(len).max(): + return row["ID"] # Set its own ID as parent + + return row["ID"] # Apply the function to each row in the DataFrame - df['parent_tag'] = df.apply(find_parent_tag, axis=1) + df["parent_tag"] = df.apply(find_parent, axis=1, args=(df,)) + return df + + +# def parent_tag_func_optimized(df): +# """ +# Optimized version of the parent tag function. + +# Args: +# df: pd.DataFrame - DataFrame for which we calculate the parent tags. + +# Returns: +# df: pd.DataFrame - DataFrame with additional parent_tag column. +# """ +# # Convert enclosed_i lists to sets for faster lookups and map IDs to their indices +# enclosed_sets = df['enclosed_i'].apply(set) +# id_to_index = {id_: idx for idx, id_ in enumerate(df['ID'])} + +# # Create an array of parent indices, initialized to self (i.e., no parent) +# parent_indices = np.arange(len(df)) + +# for row_idx in range(len(df)): +# potential_parents = [id_to_index[id_] for id_ in df['enclosed_i'].iloc[row_idx] if id_ in id_to_index] + +# if potential_parents: +# # Find the parent with the smallest enclosed set +# min_parent = min(potential_parents, key=lambda idx: len(enclosed_sets.iloc[idx])) +# parent_indices[row_idx] = min_parent + +# # Map indices back to IDs +# df['parent_tag'] = df['ID'].iloc[parent_indices].values + +# return df + + +# def parent_tag_func_optimized(df): +# """ +# Corrected and optimized vectorized implementation of parent tag function. + +# Args: +# df: pd.DataFrame - data frame for which we calculate the parent tags. + +# Returns: +# df: pd.DataFrame - Pandas data frame with additional parent tag column. +# """ +# # Create a dictionary mapping IDs to their index in the dataframe +# id_to_index = {id: idx for idx, id in enumerate(df["ID"])} + +# # Convert enclosed_i lists to sets for faster lookup +# enclosed_sets = df["enclosed_i"].apply(set) + +# # Create a matrix of containment relationships +# containment_matrix = enclosed_sets.apply( +# lambda x: [1 if id in x else 0 for id in df["ID"]] +# ) +# containment_matrix = np.array(containment_matrix.tolist()) + +# # Remove self-containment +# np.fill_diagonal(containment_matrix, 0) + +# def find_immediate_parent(row_idx): +# potential_parents = np.where(containment_matrix[:, row_idx] == 1)[0] +# if len(potential_parents) == 0: +# return row_idx # No parent found, return self + +# # Among potential parents, find the one with the smallest enclosed_i set +# parent_sizes = [len(enclosed_sets.iloc[i]) for i in potential_parents] +# immediate_parent_idx = potential_parents[np.argmin(parent_sizes)] +# return immediate_parent_idx + +# # parent_indices = [ +# # find_immediate_parent(i) +# # for i in tqdm(range(len(df)), desc="Finding immediate parents", total=len(df)) +# # ] + +# # rewrite the linline function to a normal function + +# parent_indices = [] +# for i in tqdm(range(len(df)), desc="Finding immediate parents", total=len(df)): +# immediate_parent_idx = find_immediate_parent(i) +# parent_indices.append(immediate_parent_idx) + +# # Map indices back to IDs +# df["parent_tag"] = df["ID"].iloc[parent_indices].values + +# return df + + +def parent_tag_func_optimized(df): + # Create a DataFrame with only IDs that have enclosed_i of length greater than 1 + enclosed_i = df[["ID", "enclosed_i"]] + enclosed_i = enclosed_i[enclosed_i["enclosed_i"].apply(len) > 1] + + # Initialize parent_tag with ID + df["parent_tag"] = df["ID"] + + # Create a mapping from ID to parent ID + mapping = {} + + for idx, row in enclosed_i.iterrows(): + parent_id = row["ID"] + for child_id in row["enclosed_i"]: + if child_id != parent_id: + mapping[child_id] = parent_id + + # Update parent_tag using the mapping + df["parent_tag"] = df["ID"].map(mapping).combine_first(df["parent_tag"]) return df -def correct_first_destruction(pd,output): + +def correct_first_destruction(pd, output): """ Function for correcting for the First destruction of a parent Island. @@ -182,74 +305,77 @@ def correct_first_destruction(pd,output): Returns: pd (pd.DataFrame): The new Catalogue. - + """ - pd['new_row'] = 0 + pd["new_row"] = 0 - for i in tqdm(range(0,len(pd)),total=len(pd),desc='Correcting first destruction',disable=True): + for i in tqdm( + range(0, len(pd)), + total=len(pd), + desc="Correcting first destruction", + disable=output, + ): row = pd.iloc[i] - #print(row) - enlosed_i = row['enclosed_i'] - if len(enlosed_i) > 1: + # print(row) + enlosed_i = row["enclosed_i"] + if len(enlosed_i) > 1: new_row = row.copy() - #print(i) - new_row['Death'] = pd.loc[pd['ID'] == enlosed_i[1]]['Death'] - new_row['parent_tag'] = pd.loc[pd['ID'] == enlosed_i[1]]['ID'] + # print(i) + new_row["Death"] = pd.loc[pd["ID"] == enlosed_i[1]]["Death"] + new_row["parent_tag"] = pd.loc[pd["ID"] == enlosed_i[1]]["ID"] + ## this accounts for a bug were the entire series is placed in the death column. - # not sure on the origin of this but the following corrects for it. It only occationally happends so this is not + # not sure on the origin of this but the following corrects for it. It only occationally happends so this is not # computationally expensive - #print(new_row['Death']) - if type(new_row['Death']) == pandas.core.series.Series: + # print(new_row['Death']) + if type(new_row["Death"]) == pandas.core.series.Series: # get the Death value from the first item in the Series. - # print(new_row['Death']) - new_row['Death'] = new_row['Death'].iloc[0] # - - new_row['new_row'] = 1 - new_row['ID'] = pd['ID'].max() + 1 - new_row['enclosed_i'] = [] - - + # print(new_row['Death']) + new_row["Death"] = new_row["Death"].iloc[0] # + + new_row["new_row"] = 1 + new_row["ID"] = pd["ID"].max() + 1 + new_row["enclosed_i"] = [] + # concat the new row to the dataframe. - - pd = pandas.concat([pd, new_row.to_frame().T], ignore_index=True) - #pd = pd.append(new_row,ignore_index=True) - #pd.loc[len(pd)+1] = new_row - - return pd + pd = pandas.concat([pd, new_row.to_frame().T], ignore_index=True) + # pd = pd.append(new_row,ignore_index=True) + # pd.loc[len(pd)+1] = new_row + return pd -def calculate_area_GPU(Birth,Death,row, img_gpu): +def calculate_area_GPU(Birth, Death, row, img_gpu): """Calcualtes are of source mask (for GPU) Args: - Birth (float): - Death (float): - row (pd.series): - img_gpu (cp.ndarray): + Birth (float): + Death (float): + row (pd.series): + img_gpu (cp.ndarray): Returns: area (float): the calculated area of the source mask. """ - mask = utils.get_mask_GPU(Birth,Death,row.x1,row.y1,img_gpu) - # get bounding box here + mask = utils.get_mask_GPU(Birth, Death, row.x1, row.y1, img_gpu) + # get bounding box here # evalute if mask is True on an edge. bounding_box = utils.bounding_box_gpu(mask) mask = mask.get() edge = utils.check_edge(mask) if edge: edge = True - + area = np.sum(mask) return area, edge, bounding_box -def calculate_area_CPU(Birth,Death,row, img): - - mask = utils.get_mask_CPU(row.x1,row.y1,Birth,Death,img) +def calculate_area_CPU(Birth, Death, row, img): + + mask = utils.get_mask_CPU(row.x1, row.y1, Birth, Death, img) # get bounding box here bounding_box = utils.bounding_box_cpu(mask) edge = utils.check_edge(mask) @@ -258,87 +384,116 @@ def calculate_area_CPU(Birth,Death,row, img): area = np.sum(mask) return area, edge, bounding_box - -def compute_ph_components(img,local_bg,analysis_threshold_val,lifetime_limit, - output=False,bg_map=False,area_limit=3,GPU=False,lifetime_limit_fraction=2, - mean_bg=None, IDoffset=0, box_size=None,detection_threshold=None): - + +def compute_ph_components( + img, + local_bg, + analysis_threshold_val, + lifetime_limit, + output=False, + bg_map=False, + area_limit=3, + GPU=False, + lifetime_limit_fraction=2, + mean_bg=None, + IDoffset=0, + box_size=None, + detection_threshold=None, + Cutout_X_offset=0, + Cutout_Y_offset=0, +): + global GPU_Option GPU_Option = GPU - print('Computing PH components...') + print("Computing PH components...") t0_compute_ph = time.time() - pd = cripser.computePH(-img,maxdim=0) + pd = cripser.computePH(-img, maxdim=0) t1_compute_ph = time.time() - print('Time to compute PH: {}'.format(t1_compute_ph-t0_compute_ph)) - pd = pandas.DataFrame(pd,columns=['dim','Birth','Death','x1','y1','z1','x2','y2','z2'],index=range(1,len(pd)+1)) - pd.drop(columns=['dim','z1','z2'],inplace=True) - pd['lifetime'] = pd['Death'] - pd['Birth'] - pd['Birth'] = -pd['Birth'] - pd['Death'] = -pd['Death'] - - #print("mean_bg: ",mean_bg) + print("Time to compute PH: {}".format(t1_compute_ph - t0_compute_ph)) + pd = pandas.DataFrame( + pd, + columns=["dim", "Birth", "Death", "x1", "y1", "z1", "x2", "y2", "z2"], + index=range(1, len(pd) + 1), + ) + pd.drop(columns=["dim", "z1", "z2"], inplace=True) + pd["lifetime"] = pd["Death"] - pd["Birth"] + pd["Birth"] = -pd["Birth"] + pd["Death"] = -pd["Death"] + + # print("mean_bg: ",mean_bg) # get rid of birth less than 0, helps speed up the code alittle. - mean_bg_temp = np.nanmean(local_bg)/detection_threshold - #print('mean_bg: ',mean_bg_temp) + # mean_bg_temp = np.nanmean(local_bg)/detection_threshold + # print('mean_bg: ',mean_bg_temp) # get rid of alot of defintly not sources. - pd = pd[pd['Birth']>mean_bg_temp] + # pd = pd[pd['Birth']>mean_bg_temp] - pd['bg'] = 0 - pd['edge_flag'] = 0 - pd['mean_bg'] = 0 # this is the mean of the background + pd["bg"] = 0 + pd["edge_flag"] = 0 + pd["mean_bg"] = 0 # this is the mean of the background # assign each row the local bg valuw from the map. if bg_map: - + # this is a slow function. Can it be improved? # for each row we need to assign the local bg value. - pd['bg'] = pd['bg'].astype(float) - pd['mean_bg'] = pd['mean_bg'].astype(float) + pd["bg"] = pd["bg"].astype(float) + pd["mean_bg"] = pd["mean_bg"].astype(float) for index, row in pd.iterrows(): - - pd.loc[index,'bg'] = float(background.get_bg_value_from_result_image((int(row.x1),int(row.y1)), box_size, local_bg)) - #print('Detection_Thresh: ',row['bg']) - #print('Birth :',row['Birth']) - pd.loc[index,'mean_bg'] = float(background.get_bg_value_from_result_image((int(row.x1),int(row.y1)), box_size, mean_bg)) - #print(row['mean_bg']) + + pd.loc[index, "bg"] = float( + background.get_bg_value_from_result_image( + (int(row.x1 + Cutout_Y_offset), int(row.y1 + Cutout_X_offset)), + box_size, + local_bg, + ) + ) + pd.loc[index, "mean_bg"] = float( + background.get_bg_value_from_result_image( + (int(row.x1 + Cutout_Y_offset), int(row.y1 + Cutout_X_offset)), + box_size, + mean_bg, + ) + ) + # evaluate if the death value is below the analysis threshold value at the birth point. # if it is then set the death value to the analysis threshold value. - - Anal_val = background.get_bg_value_from_result_image((int(row.x1),int(row.y1)), box_size, analysis_threshold_val) - - if row['Death'] < Anal_val: - pd.loc[index,'Death'] = Anal_val - - - + + Anal_val = background.get_bg_value_from_result_image( + (int(row.x1 + Cutout_Y_offset), int(row.y1 + Cutout_X_offset)), + box_size, + analysis_threshold_val, + ) + + if row["Death"] < Anal_val: + pd.loc[index, "Death"] = Anal_val + else: # no bg map so just assign the local bg value. asn this should be single value. - pd['bg'] = local_bg - pd['mean_bg'] = mean_bg - + pd["bg"] = local_bg + pd["mean_bg"] = mean_bg + for index, row in pd.iterrows(): - if row['Death'] < analysis_threshold_val: - pd.loc[index,'Death'] = analysis_threshold_val - - #print('Before Cull',len(pd)) - #print('mean_bg: ',np.mean(pd['bg'])) - pd = pd[pd['Birth'] > pd['bg']] # maybe this should be at the beginning. + if row["Death"] < analysis_threshold_val: + pd.loc[index, "Death"] = analysis_threshold_val + + # print('Before Cull',len(pd)) + # print(pd['bg']) + # print('mean_bg: ',np.mean(pd['bg'])) + pd = pd[pd["Birth"] > pd["bg"]] # maybe this should be at the beginning. # also make sure Death is less than Birth - #pd = pd[pd['Death'] < pd['Birth']] - #print('After Cull',len(pd)) + # pd = pd[pd['Death'] < pd['Birth']] + # print('After Cull',len(pd)) # if bg_map: - + # list_of_index_to_drop = [] - - + # for index, row in pd.iterrows(): # # check if local_bg is a map or a value # if row['Birth'] < local_bg[int(row.x1),int(row.y1)]: # list_of_index_to_drop.append(index) - + # pd.drop(list_of_index_to_drop,inplace=True) - - - # # for each row evaluate if death is below analysis thresholdval map value at its birth point. if its below then set Death to bg map value. + + # # for each row evaluate if death is below analysis thresholdval map value at its birth point. if its below then set Death to bg map value. # for index, row in pd.iterrows(): # Analy_val = analysis_threshold_val[int(row.y1),int(row.x1)] # if row['Death'] < Analy_val: @@ -346,75 +501,85 @@ def compute_ph_components(img,local_bg,analysis_threshold_val,lifetime_limit, # # assign each row the local bg value # row['bg'] = local_bg[int(row.y1),int(row.x1)] # row['mean_bg'] = mean_bg[int(row.y1),int(row.x1)] - + # else: # print(local_bg) # pd = pd[pd['Birth']>local_bg] # maybe this should be at the beginning. # pd['Death'] = np.where(pd['Death'] < analysis_threshold_val, analysis_threshold_val, pd['Death']) # pd['bg'] = local_bg # pd['mean_bg'] = mean_bg - - pd['lifetime'] = abs(pd['Death'] - pd['Birth']) - - pd['lifetimeFrac'] = pd['Birth']/pd['Death'] - pd = pd[pd['lifetimeFrac']>lifetime_limit_fraction] - pd = pd[pd['lifetime'] > lifetime_limit] - pd.sort_values(by='lifetime',ascending=False,inplace=True,ignore_index=True) - - pd['ID'] = pd.index + IDoffset - - + + pd["lifetime"] = abs(pd["Death"] - pd["Birth"]) + + pd["lifetimeFrac"] = pd["Birth"] / pd["Death"] + pd = pd[pd["lifetimeFrac"] > lifetime_limit_fraction] + pd = pd[pd["lifetime"] > lifetime_limit] + pd.sort_values(by="lifetime", ascending=False, inplace=True, ignore_index=True) + + pd["ID"] = pd.index + IDoffset + + # print('PD: ',pd) # begins here - + if len(pd) > 0: - - area_list=[] - edge_list=[] + + area_list = [] + edge_list = [] bbox1 = [] bbox2 = [] bbox3 = [] bbox4 = [] - if GPU_Option ==True: + if GPU_Option == True: if GPU_AVAILABLE == True: - img_gpu = cp.asarray(img,dtype=cp.float64) + img_gpu = cp.asarray(img, dtype=cp.float64) # Calculate area and enforce area limit Single Process. - #print('Calculating area with GPU...') + # print('Calculating area with GPU...') t0 = time.time() - - for i in tqdm(range(0,len(pd)),total=len(pd),desc='Calculating area',disable=not output): - + + for i in tqdm( + range(0, len(pd)), + total=len(pd), + desc="Calculating area", + disable=not output, + ): + row = pd.iloc[i] Birth = row.Birth Death = row.Death - area, edge, bbox = calculate_area_GPU(Birth,Death,row,img_gpu) + area, edge, bbox = calculate_area_GPU(Birth, Death, row, img_gpu) area_list.append(area) edge_list.append(edge) bbox1.append(bbox[0].get()) bbox2.append(bbox[1].get()) bbox3.append(bbox[2].get()) bbox4.append(bbox[3].get()) - + t1 = time.time() - #print('Time to calculate area and inital bbox: {}'.format(t1-t0)) - pd['area'] = area_list - pd['edge_flag'] = edge_list - pd['bbox1'] = bbox1 - pd['bbox2'] = bbox2 - pd['bbox3'] = bbox3 - pd['bbox4'] = bbox4 - pd = pd[pd['area']>area_limit] - + # print('Time to calculate area and inital bbox: {}'.format(t1-t0)) + pd["area"] = area_list + pd["edge_flag"] = edge_list + pd["bbox1"] = bbox1 + pd["bbox2"] = bbox2 + pd["bbox3"] = bbox3 + pd["bbox4"] = bbox4 + pd = pd[pd["area"] > area_limit] + else: - - #print('Calculating area with CPU...') + + # print('Calculating area with CPU...') t0 = time.time() - for i in tqdm(range(0,len(pd)),total=len(pd),desc='Calculating area',disable=not output): - + for i in tqdm( + range(0, len(pd)), + total=len(pd), + desc="Calculating area", + disable=not output, + ): + row = pd.iloc[i] Birth = row.Birth Death = row.Death - area, edge, bbox = calculate_area_CPU(Birth,Death,row,img) + area, edge, bbox = calculate_area_CPU(Birth, Death, row, img) area_list.append(area) edge_list.append(edge) bbox1.append(bbox[0]) @@ -423,14 +588,13 @@ def compute_ph_components(img,local_bg,analysis_threshold_val,lifetime_limit, bbox4.append(bbox[3]) t1 = time.time() - #print('Time to calculate area and inital bbox: {}'.format(t1-t0)) - pd['area'] = area_list - pd['edge_flag'] = edge_list - pd['bbox1'] = bbox1 - pd['bbox2'] = bbox2 - pd['bbox3'] = bbox3 - pd['bbox4'] = bbox4 - pd = pd[pd['area']>area_limit] - + # print('Time to calculate area and inital bbox: {}'.format(t1-t0)) + pd["area"] = area_list + pd["edge_flag"] = edge_list + pd["bbox1"] = bbox1 + pd["bbox2"] = bbox2 + pd["bbox3"] = bbox3 + pd["bbox4"] = bbox4 + pd = pd[pd["area"] > area_limit] + return pd - \ No newline at end of file diff --git a/DRUID/src/source.py b/DRUID/src/source.py index 1fce7e9..8b057ee 100644 --- a/DRUID/src/source.py +++ b/DRUID/src/source.py @@ -1,4 +1,3 @@ - """ File: src/source.py Author: Rhys Shaw @@ -15,92 +14,89 @@ from scipy.ndimage import label import pdb from scipy.ndimage import binary_dilation + try: import cupy as cp from cupyx.scipy.ndimage import label as cupy_label except: - - pass - + pass -def create_params_df(cutup : bool, params : list): - - ''' - - Creates a pandas dataframe from the parameters. - - ''' - - params = pd.DataFrame(params,columns=['ID', - 'Birth', - 'Death', - 'x1', - 'y1', - 'x2', - 'y2', - 'Flux_total', - 'Flux_peak', - 'Flux_correction_factor', - 'Area', - 'Xc', - 'Yc', - 'bbox1', - 'bbox2', - 'bbox3', - 'bbox4', - 'Maj', - 'Min', - 'Pa', - 'parent_tag', - 'Class', - 'SNR', - 'Noise', - 'X0_cutout', - 'Y0_cutout', - 'mean_bg', - 'Edge_flag', - 'contour', - 'enclosed_i']) - - return params - - - - - - - - - -def large_mask_red_image_procc_GPU(Birth,Death,x1,y1,image,X0,Y0): - ''' +def create_params_df(cutup: bool, params: list): + """ + + Creates a pandas dataframe from the parameters. + + """ + + params = pd.DataFrame( + params, + columns=[ + "ID", + "Birth", + "Death", + "x1", + "y1", + "x2", + "y2", + "Flux_total", + "Flux_total_err", + "Flux_peak", + "Flux_correction_factor", + "Area", + "Xc", + "Yc", + "bbox1", + "bbox2", + "bbox3", + "bbox4", + "Maj", + "Min", + "Pa", + "parent_tag", + "Class", + "SNR", + "Noise", + "X0_cutout", + "Y0_cutout", + "mean_bg", + "Edge_flag", + "contour", + "enclosed_i", + ], + ) + + return params + + +def large_mask_red_image_procc_GPU(Birth, Death, x1, y1, image, X0, Y0): + """ Does all gpu processing for the large mask. - + return the red_image and red_mask and the bounding box. - - ''' - - mask = cp.zeros(image.shape,dtype=cp.bool_) - mask = cp.logical_or(mask,cp.logical_and(image <= Birth, image > Death)) - + + """ + + mask = cp.zeros(image.shape, dtype=cp.bool_) + mask = cp.logical_or(mask, cp.logical_and(image <= Birth, image > Death)) + # mask_enclosed = self.get_enclosing_mask_gpu(y1,x1,mask) labeled_mask, num_features = cupy_label(mask) - + # Check if the specified pixel is within the mask if 0 <= y1 < mask.shape[1] and 0 <= x1 < mask.shape[0]: label_at_pixel = labeled_mask[x1, y1] - #print(x1) - #print(y1) - #print(label_at_pixel) + # print(x1) + # print(y1) + # print(label_at_pixel) if label_at_pixel != 0: # Extract the connected component containing the specified pixel - component_mask = (labeled_mask == label_at_pixel) - #plt.imshow(component_mask.get()) - #plt.show() - #pdb.set_trace() + component_mask = labeled_mask == label_at_pixel + # plt.imshow(component_mask.get()) + # plt.show() + # pdb.set_trace() non_zero_indices = cp.nonzero(component_mask) # Extract minimum and maximum coordinates @@ -108,52 +104,45 @@ def large_mask_red_image_procc_GPU(Birth,Death,x1,y1,image,X0,Y0): ymin = cp.min(non_zero_indices[0]) xmax = cp.max(non_zero_indices[1]) ymax = cp.max(non_zero_indices[0]) - + # correct the bounding box for the cutout. - xmin = xmin #+ Y0 - xmax = xmax #+ Y0 - ymin = ymin #+ X0 - ymax = ymax #+ X0 - - # images are not being cropped? - - red_image = image[ymin:ymax+1, xmin:xmax+1] - red_mask = component_mask[ymin:ymax+1, xmin:xmax+1] - - return red_image, red_mask, xmin, xmax, ymin, ymax + xmin = xmin # + Y0 + xmax = xmax # + Y0 + ymin = ymin # + X0 + ymax = ymax # + X0 + # images are not being cropped? + red_image = image[ymin : ymax + 1, xmin : xmax + 1] + red_mask = component_mask[ymin : ymax + 1, xmin : xmax + 1] + return red_image, red_mask, xmin, xmax, ymin, ymax +def large_mask_red_image_procc_CPU(Birth, Death, x1, y1, image): + """ + Does all gpu processing for the large mask. + return the red_image and red_mask and the bounding box. + """ -def large_mask_red_image_procc_CPU(Birth,Death,x1,y1,image): - ''' - Does all gpu processing for the large mask. - - return the red_image and red_mask and the bounding box. - - ''' - - mask = np.zeros(image.shape) - mask = np.logical_or(mask,np.logical_and(image <= Birth, image > Death)) + mask = np.logical_or(mask, np.logical_and(image <= Birth, image > Death)) # mask_enclosed = self.get_enclosing_mask_gpu(y1,x1,mask) labeled_mask, num_features = label(mask) - + # Check if the specified pixel is within the mask if 0 <= y1 < mask.shape[1] and 0 <= x1 < mask.shape[0]: - + label_at_pixel = labeled_mask[x1, y1] - + if label_at_pixel != 0: # Extract the connected component containing the specified pixel - component_mask = (labeled_mask == label_at_pixel) - #plt.imshow(component_mask.get()) - #plt.show() - #pdb.set_trace() + component_mask = labeled_mask == label_at_pixel + # plt.imshow(component_mask.get()) + # plt.show() + # pdb.set_trace() non_zero_indices = np.nonzero(component_mask) # Extract minimum and maximum coordinates @@ -161,179 +150,207 @@ def large_mask_red_image_procc_CPU(Birth,Death,x1,y1,image): ymin = np.min(non_zero_indices[0]) xmax = np.max(non_zero_indices[1]) ymax = np.max(non_zero_indices[0]) - + # images are not being cropped? - - red_image = image[ymin:ymax+1, xmin:xmax+1] - red_mask = component_mask[ymin:ymax+1, xmin:xmax+1] - + + red_image = image[ymin : ymax + 1, xmin : xmax + 1] + red_mask = component_mask[ymin : ymax + 1, xmin : xmax + 1] + return red_image, red_mask, xmin, xmax, ymin, ymax -#import cv2 -#import numpy as np -#def dilate_mask_circular(mask, radius): - # Create a circular structuring element using cv2.getStructuringElement +# import cv2 +# import numpy as np + +# def dilate_mask_circular(mask, radius): +# Create a circular structuring element using cv2.getStructuringElement # circular_kernel = cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (2*radius+1, 2*radius+1)) - # Perform dilation using cv2.dilate +# Perform dilation using cv2.dilate # dilated_mask = cv2.dilate(mask.astype(np.uint8), circular_kernel)# # return dilated_mask - - -def curve_of_growth_dilation(mask,image): - converged = False - i = 0 - dilated_mask = np.zeros(image.shape, dtype=bool) - mask_before = mask - while converged == False: - if i == 50: - #print('max iterations reached') - return dilated_mask.astype(int) - # print(i) - i += 1 - - dilated_mask = dilate_mask_circular(mask_before,1) - - # check if the mask has converged - flux = np.sum(image*dilated_mask) - flux_old = np.sum(image*mask_before) - diff = (flux - flux_old) - #print(diff) - if diff<0: - converged = True - dilated_mask = mask_before - - else: - mask_before = dilated_mask - #print('converged {}'.format(i)) - return dilated_mask.astype(int) - - -def measure_source_properties(use_gpu,catalogue=None,cutout=None,background_map=None,output=None,cutupts=None, mode='optical',header=None): - ''' - + + +# def curve_of_growth_dilation(mask,image): +# converged = False +# i = 0 +# dilated_mask = np.zeros(image.shape, dtype=bool) +# mask_before = mask +# while converged == False: +# if i == 50: +# #print('max iterations reached') +# return dilated_mask.astype(int) +# # print(i) +# i += 1 + +# dilated_mask = dilate_mask_circular(mask_before,1) + +# # check if the mask has converged +# flux = np.sum(image*dilated_mask) +# flux_old = np.sum(image*mask_before) +# diff = (flux - flux_old) +# #print(diff) +# if diff<0: +# converged = True +# dilated_mask = mask_before + +# else: +# mask_before = dilated_mask +# #print('converged {}'.format(i)) +# return dilated_mask.astype(int) + + +def measure_source_properties( + use_gpu, + catalogue=None, + cutout=None, + background_map=None, + output=None, + cutupts=None, + mode="optical", + header=None, + sigma=5, +): + """ + Characterising the Source Assuming the input image of of the format of a optical astronomical image. - + # needs to work on cutup images. # work on the whole image. - ''' + """ if header == None: - mode = 'other' - - if mode == 'Radio': - print('Radio mode selected') + mode = "other" + + if mode == "Radio": + print("Radio mode selected") Beam, BMAJ, BMIN, BPA = utils.calculate_beam(header=header) - - if mode == 'optical': - print('Optical mode selected, modeling the PSF as a gaussian.') + + if mode == "optical": + print("Optical mode selected, modeling the PSF as a gaussian.") psf_fwhm_p = utils.get_psf_FWHM(header) - #EFFRON = utils.get_EFFRON(header) - #EFFGAIN = utils.get_EFFGAIN(header) - #EXPTIME = utils.get_EXPTIME(header) - + # EFFRON = utils.get_EFFRON(header) + # EFFGAIN = utils.get_EFFGAIN(header) + # EXPTIME = utils.get_EXPTIME(header) + if use_gpu: - + try: - + import cupy as cp - + except: - - raise ImportError('cupy not installed. GPU acceleration not possible.') - + + raise ImportError("cupy not installed. GPU acceleration not possible.") + image = cutout if use_gpu: image_gpu = cp.asarray(image, dtype=cp.float64) - - - Birth = catalogue['Birth'].to_numpy() - Death = catalogue['Death'].to_numpy() - parent_tag = catalogue['parent_tag'].to_numpy() - Class = catalogue['Class'].to_numpy() - bg = catalogue['bg'].to_numpy() - X0 = catalogue['X0_cutout'].to_numpy() - Y0 = catalogue['Y0_cutout'].to_numpy() - mean_bg = catalogue['mean_bg'].to_numpy() - enclosed_i = catalogue['enclosed_i'].to_numpy() - IDs = catalogue['ID'].to_numpy() - Edge_flags = catalogue['edge_flag'].to_numpy() - bbox1_og = catalogue['bbox1'].to_list() - bbox2_og = catalogue['bbox2'].to_list() - bbox3_og = catalogue['bbox3'].to_list() - bbox4_og = catalogue['bbox4'].to_list() - - x1 = catalogue['x1'].to_numpy() #- 1 #- X0 - y1 = catalogue['y1'].to_numpy() #- 1#- X0 - x2 = catalogue['x2'].to_numpy() #- 1#- Y0 - y2 = catalogue['y2'].to_numpy() #- 1#- X0 - - + + Birth = catalogue["Birth"].to_numpy() + Death = catalogue["Death"].to_numpy() + parent_tag = catalogue["parent_tag"].to_numpy() + Class = catalogue["Class"].to_numpy() + bg = catalogue["bg"].to_numpy() + X0 = catalogue["X0_cutout"].to_numpy() + Y0 = catalogue["Y0_cutout"].to_numpy() + mean_bg = catalogue["mean_bg"].to_numpy() + enclosed_i = catalogue["enclosed_i"].to_numpy() + IDs = catalogue["ID"].to_numpy() + Edge_flags = catalogue["edge_flag"].to_numpy() + bbox1_og = catalogue["bbox1"].to_list() + bbox2_og = catalogue["bbox2"].to_list() + bbox3_og = catalogue["bbox3"].to_list() + bbox4_og = catalogue["bbox4"].to_list() + + x1 = catalogue["x1"].to_numpy() # - 1 #- X0 + y1 = catalogue["y1"].to_numpy() # - 1#- X0 + x2 = catalogue["x2"].to_numpy() # - 1#- Y0 + y2 = catalogue["y2"].to_numpy() # - 1#- X0 + # map X0 and Y0 to the cutout number params = [] polygons = [] - - import matplotlib.pylab as plt - for i, source in tqdm(enumerate(Birth),total=len(Birth),desc='Calculating Source Properties..',disable=not output): - + # import matplotlib.pylab as plt + + for i, source in tqdm( + enumerate(Birth), + total=len(Birth), + desc="Calculating Source Properties..", + disable=not output, + ): + if use_gpu == True: - - cropped_image_gpu = image_gpu[bbox1_og[i]-1:bbox3_og[i]+1,bbox2_og[i]-1:bbox4_og[i]+1] - + + cropped_image_gpu = image_gpu[ + bbox1_og[i] - 1 : bbox3_og[i] + 1, bbox2_og[i] - 1 : bbox4_og[i] + 1 + ] + try: - red_image, red_mask, xmin, xmax, ymin, ymax = large_mask_red_image_procc_GPU(Birth[i],Death[i], - x1[i]-bbox1_og[i]+1, - y1[i]-bbox2_og[i]+1, - cropped_image_gpu, - X0[i],Y0[i]) + red_image, red_mask, xmin, xmax, ymin, ymax = ( + large_mask_red_image_procc_GPU( + Birth[i], + Death[i], + x1[i] - bbox1_og[i] + 1, + y1[i] - bbox2_og[i] + 1, + cropped_image_gpu, + X0[i], + Y0[i], + ) + ) except: - - print('Error in GPU processing!') - print('Source ID: ',IDs[i]) - print('x1:',x1[i]-bbox1_og[i]+1) - print('y1:',y1[i]-bbox2_og[i]+1) - print('Birth:',Birth[i]) - print('Death:',Death[i]) - #print('Cropped_image',cropped_image_gpu.get()) - + + print("Error in GPU processing!") + print("Source ID: ", IDs[i]) + print("x1:", x1[i] - bbox1_og[i] + 1) + print("y1:", y1[i] - bbox2_og[i] + 1) + print("Birth:", Birth[i]) + print("Death:", Death[i]) + # print('Cropped_image',cropped_image_gpu.get()) + red_mask = red_mask.astype(int) - + red_image = red_image.get() red_mask = red_mask.get() xmin = xmin.get() + bbox2_og[i] - xmax = xmax.get() + bbox2_og[i] + xmax = xmax.get() + bbox2_og[i] ymin = ymin.get() + bbox1_og[i] ymax = ymax.get() + bbox1_og[i] - + else: - cropped_image = image[bbox1_og[i]-1:bbox3_og[i]+1,bbox2_og[i]-1:bbox4_og[i]+1] + cropped_image = image[ + bbox1_og[i] - 1 : bbox3_og[i] + 1, bbox2_og[i] - 1 : bbox4_og[i] + 1 + ] try: - red_image, red_mask, xmin, xmax, ymin, ymax = large_mask_red_image_procc_CPU(Birth[i], - Death[i], - int(x1[i])-int(bbox1_og[i])+1, - int(y1[i])-int(bbox2_og[i])+1, - cropped_image) + red_image, red_mask, xmin, xmax, ymin, ymax = ( + large_mask_red_image_procc_CPU( + Birth[i], + Death[i], + int(x1[i]) - int(bbox1_og[i]) + 1, + int(y1[i]) - int(bbox2_og[i]) + 1, + cropped_image, + ) + ) except: print("Error in CPU processing!") - print('Source ID: ',IDs[i]) - print('x1:',x1[i]-bbox1_og[i]+1) - print('y1:',y1[i]-bbox2_og[i]+1) - print('Birth:',Birth[i]) - print('Death:',Death[i]) - + print("Source ID: ", IDs[i]) + print("x1:", x1[i] - bbox1_og[i] + 1) + print("y1:", y1[i] - bbox2_og[i] + 1) + print("Birth:", Birth[i]) + print("Death:", Death[i]) + red_mask = red_mask.astype(int) xmin = xmin + bbox2_og[i] xmax = xmax + bbox2_og[i] ymin = ymin + bbox1_og[i] ymax = ymax + bbox1_og[i] - + # DILATION TEST Not Working. - #print(red_mask) - + # print(red_mask) + # #try: # if Class[i] == 0: # if xmin > 50 and ymin > 50: @@ -354,8 +371,7 @@ def measure_source_properties(use_gpu,catalogue=None,cutout=None,background_map= # plt.imshow(red_image,cmap='gray',origin='lower',vmin=1E-12,vmax=1E-10) # plt.contour(red_mask) # plt.savefig('red_image.png') - - + # #print('red_image shape',red_image.shape) # # dilate the mask until the flux converges. # # check that the shapes are the same @@ -368,9 +384,9 @@ def measure_source_properties(use_gpu,catalogue=None,cutout=None,background_map= # plt.imshow(red_image,cmap='gray',origin='lower',vmin=1E-12,vmax=1E-10) # plt.contour(red_mask) # plt.savefig('red_image_after.png') - + # pdb.set_trace() - + # #print('red_mask shape',red_mask.shape) # # reduce the mask and recalculate the boundig box. # # plt.imshow(red_mask) @@ -382,275 +398,339 @@ def measure_source_properties(use_gpu,catalogue=None,cutout=None,background_map= # ymax = ymaxn - expanx + ymax # red_mask = red_mask[yminn:ymaxn,xminn:xmaxn] # red_image = red_image[yminn:ymaxn,xminn:xmaxn] - - - contour = utils._get_polygons_in_bbox(xmin,xmax,ymin,ymax,x1[i],y1[i],Birth[i],Death[i],red_mask,0,0) - source_props = utils.get_region_props(red_mask,image=red_image) + contour = utils._get_polygons_in_bbox( + xmin, xmax, ymin, ymax, x1[i], y1[i], Birth[i], Death[i], red_mask, 0, 0 + ) + + # TEST # + # from matplotlib.pylab import plt + + # plt.imshow(red_image, cmap="gray", origin="lower") + # plt.imshow(red_mask, alpha=0.5, origin="lower") + # # plot the countour of the red_mask + # plt.plot(contour[:, 1] - xmin, contour[:, 0] - ymin, color="red") + # plt.show() + + # # # # + + source_props = utils.get_region_props(red_mask, image=red_image) source_props = utils.props_to_dict(source_props[0]) - #print(background_map) - background_map = bg[i] - #print('Mean bg',mean_bg) + # print(background_map) + background_map = bg[i] / sigma # background_map + # print('Mean bg',mean_bg) mean_bg_s = mean_bg[i] - - background_map = np.random.normal(mean_bg_s,background_map,red_mask.shape) - - red_background_mask = np.where(red_mask == 0, np.nan, red_mask*background_map) - - Noise = abs(np.nansum(red_background_mask)) - #print('Noise',Noise) - peak_coords = np.where(red_image == source_props['max_intensity']) + + background_map = np.random.normal(mean_bg_s, background_map, red_mask.shape) + + red_background_mask = np.where(red_mask == 0, np.nan, red_mask * background_map) + + Noise = np.nansum(red_background_mask) + # print('Noise',Noise) + peak_coords = np.where(red_image == source_props["max_intensity"]) y_peak_loc = peak_coords[0][0] x_peak_loc = peak_coords[1][0] - #print('x_peak_loc',x_peak_loc) - #print('y_peak_loc',y_peak_loc) + # print('x_peak_loc',x_peak_loc) + # print('y_peak_loc',y_peak_loc) shape = red_image.shape - - Area = source_props['area'] - #print('Class',Class[i]) - - - - if mode == 'Radio' or mode == 'optical': + + Area = source_props["area"] + # print('Class',Class[i]) + SNR = 0 # default value gets overwritten if calculated + if mode == "Radio" or mode == "optical": shape = red_image.shape # if shape is smaller than 100 in any direction then we add padding evenly to each side. if shape[0] < 100: - pad = int((100 - shape[0])/2) - red_image = np.pad(red_image,((pad,pad),(0,0)),mode='constant',constant_values=0) - red_background_mask = np.pad(red_background_mask,((pad,pad),(0,0)),mode='constant',constant_values=0) - red_mask = np.pad(red_mask,((pad,pad),(0,0)),mode='constant',constant_values=0) + pad = int((100 - shape[0]) / 2) + red_image = np.pad( + red_image, ((pad, pad), (0, 0)), mode="constant", constant_values=0 + ) + red_background_mask = np.pad( + red_background_mask, + ((pad, pad), (0, 0)), + mode="constant", + constant_values=0, + ) + red_mask = np.pad( + red_mask, ((pad, pad), (0, 0)), mode="constant", constant_values=0 + ) shape = red_image.shape y = y_peak_loc + pad - + else: x = y_peak_loc y = x_peak_loc - + if shape[1] < 100: - pad = int((100 - shape[1])/2) - red_image = np.pad(red_image,((0,0),(pad,pad)),mode='constant',constant_values=0) - red_background_mask = np.pad(red_background_mask,((0,0),(pad,pad)),mode='constant',constant_values=0) - red_mask = np.pad(red_mask,((0,0),(pad,pad)),mode='constant',constant_values=0) + pad = int((100 - shape[1]) / 2) + red_image = np.pad( + red_image, ((0, 0), (pad, pad)), mode="constant", constant_values=0 + ) + red_background_mask = np.pad( + red_background_mask, + ((0, 0), (pad, pad)), + mode="constant", + constant_values=0, + ) + red_mask = np.pad( + red_mask, ((0, 0), (pad, pad)), mode="constant", constant_values=0 + ) shape = red_image.shape x = x_peak_loc + pad - + else: x = y_peak_loc y = x_peak_loc - - # print(shape) - if mode == 'optical': - - MAJ = psf_fwhm_p/2.355 # maybe use full expression in future. + # print(shape) + + if mode == "optical": + + MAJ = psf_fwhm_p / 2.355 # maybe use full expression in future. MIN = MAJ - BPA = 0 - - Model_Beam = utils.model_beam_func(source_props['max_intensity'],shape, - x,y,MAJ, - MIN,BPA) - #plt.figure(figsize=(10,10)) - #plt.imshow(Model_Beam) - #plt.scatter(x,y,color='red',marker='x',s=10) + BPA = 0 + + Model_Beam = utils.model_beam_func( + source_props["max_intensity"], shape, x, y, MAJ, MIN, BPA + ) + # plt.figure(figsize=(10,10)) + # plt.imshow(Model_Beam) + # plt.scatter(x,y,color='red',marker='x',s=10) # plot the contour of the red_mask - #plt.contour(red_mask) - #plt.savefig('Model_Beam_test.png') - #pdb.set_trace() - + # plt.contour(red_mask) + # plt.savefig('Model_Beam_test.png') + # pdb.set_trace() + else: BMAJ = BMAJ BMIN = BMIN - - Model_Beam = utils.model_beam_func(source_props['max_intensity'],shape, - x,y,BMAJ/2, - BMIN/2,BPA) - - - if mode == 'Radio': - Flux_total = np.nansum(red_mask*red_image - red_background_mask)/Beam + + Model_Beam = utils.model_beam_func( + source_props["max_intensity"], shape, x, y, BMAJ / 2, BMIN / 2, BPA + ) + + if mode == "Radio": + Flux_total = ( + np.nansum(red_mask * red_image - red_background_mask) / Beam + ) else: - - Flux_total = np.nansum(red_mask*red_image - red_background_mask) - - Flux_peak = np.nanmax(red_mask*red_image) - red_background_mask[y_peak_loc,x_peak_loc] + + Flux_total = np.nansum(red_mask * red_image - red_background_mask) + + Flux_peak = ( + np.nanmax(red_mask * red_image) + - red_background_mask[y_peak_loc, x_peak_loc] + ) Flux_correction_factor = utils.flux_correction_factor(red_mask, Model_Beam) - + Flux_total_err = 0 + SNR = 0 if Area < 100: # we assume that the flux cannot be corrected - Flux_total = Flux_total*Flux_correction_factor - - # remove the 1pixel padding - if mode == 'Radio': - padding = 1 - else: - padding = 0 - + Flux_total = Flux_total * Flux_correction_factor + + # Calcualte Radio Flux error + if mode == "Radio": + # adapted from https://github.com/mhardcastle/radioflux/blob/master/radioflux/radioflux.py + rms = bg[i] / sigma + gfactor = 2 * np.sqrt(2 * np.log(2)) + Beam_area = 2 * np.pi * (BMAJ * BMIN) / gfactor + + Flux_total_err = rms * np.sqrt(Area / Beam_area) + SNR = Flux_total / Flux_total_err + + elif mode == "optical": + SNR = 0 + + padding = 0 + else: - - Flux_total = np.nansum(red_mask*red_image - red_background_mask) - Flux_peak = np.nanmax(red_mask*red_image) - red_background_mask[y_peak_loc,x_peak_loc] + + Flux_total = np.nansum(red_mask * red_image - red_background_mask) + Flux_peak = ( + np.nanmax(red_mask * red_image) + - red_background_mask[y_peak_loc, x_peak_loc] + ) Flux_correction_factor = np.nan padding = 0 - - # add some estimation of the PSF. To account for the loss of flux due to the PSF cutting at low SNR. - - - SNR = Flux_total/(Area*bg[i]) - - # Noise = np.std(red_background_mask) - - Xc = source_props['centroid'][1] + xmin - padding - Yc = source_props['centroid'][0] + ymin - padding - - - Xc = Xc #+ X0[i] - Yc = Yc #+ Y0[i] - bbox1 = source_props['bbox'][0] + xmin - padding - bbox2 = source_props['bbox'][1] + ymin - padding - bbox3 = source_props['bbox'][2] + xmin - padding - bbox4 = source_props['bbox'][3] + ymin - padding - - Maj = source_props['major_axis_length'] - Min = source_props['minor_axis_length'] - Pa = source_props['orientation'] - + Flux_total_err = 0 + + # SNR = Flux_total / Flux_total_err + # Noise = np.std(red_background_mask) + + Xc = source_props["centroid"][1] + xmin - padding + Yc = source_props["centroid"][0] + ymin - padding + + Xc = Xc # + X0[i] + Yc = Yc # + Y0[i] + bbox1 = source_props["bbox"][0] + xmin - padding + bbox2 = source_props["bbox"][1] + ymin - padding + bbox3 = source_props["bbox"][2] + xmin - padding + bbox4 = source_props["bbox"][3] + ymin - padding + + Maj = source_props["major_axis_length"] + Min = source_props["minor_axis_length"] + Pa = source_props["orientation"] + + # print(Flux_total_err) if Edge_flags[i] != 1: - - params.append([IDs[i], - Birth[i], - Death[i], - x1[i], - y1[i], - x2[i], - y2[i], - Flux_total, - Flux_peak, # this should be the peak flux. - Flux_correction_factor, - Area, - Xc, - Yc, - bbox1, - bbox2, - bbox3, - bbox4, - Maj, - Min, - Pa, - parent_tag[i], - Class[i], - SNR, - Noise, - X0[i], - Y0[i], - mean_bg_s, - Edge_flags[i], - contour, - enclosed_i[i]]) - + + params.append( + [ + IDs[i], + Birth[i], + Death[i], + x1[i], + y1[i], + x2[i], + y2[i], + Flux_total, + Flux_total_err, + Flux_peak, # this should be the peak flux. + Flux_correction_factor, + Area, + Xc, + Yc, + bbox1, + bbox2, + bbox3, + bbox4, + Maj, + Min, + Pa, + parent_tag[i], + Class[i], + SNR, + Noise, + X0[i], + Y0[i], + mean_bg_s, + Edge_flags[i], + contour, + enclosed_i[i], + ] + ) + polygons.append(contour) - #except: + # except: # print('Error in optical characteristing!') # print('Source ID: ',IDs[i]) # print('Skipping source...') # continue - #print(params) - #print(len(params[0])) - - return create_params_df(False,params), polygons + # print(params) + # print(len(params[0])) + + return create_params_df(False, params), polygons +def create_polygons(use_gpu, catalogue=None, cutout=None, output=None, cutupts=None): -def create_polygons(use_gpu,catalogue=None,cutout=None,output=None,cutupts=None): - if use_gpu: - + try: - + import cupy as cp - + except: - - raise ImportError('cupy not installed. GPU acceleration not possible.') - + + raise ImportError("cupy not installed. GPU acceleration not possible.") + image = cutout if use_gpu: image_gpu = cp.asarray(image, dtype=cp.float64) - - - Birth = catalogue['Birth'].to_numpy() - Death = catalogue['Death'].to_numpy() - - X0 = catalogue['X0_cutout'].to_numpy() - Y0 = catalogue['Y0_cutout'].to_numpy() - IDs = catalogue['ID'].to_numpy() - bbox1_og = catalogue['bbox1'].to_list() - bbox2_og = catalogue['bbox2'].to_list() - bbox3_og = catalogue['bbox3'].to_list() - bbox4_og = catalogue['bbox4'].to_list() - - x1 = catalogue['x1'].to_numpy() #- 1 #- X0 - y1 = catalogue['y1'].to_numpy() #- 1#- X0 + + Birth = catalogue["Birth"].to_numpy() + Death = catalogue["Death"].to_numpy() + + X0 = catalogue["X0_cutout"].to_numpy() + Y0 = catalogue["Y0_cutout"].to_numpy() + IDs = catalogue["ID"].to_numpy() + bbox1_og = catalogue["bbox1"].to_list() + bbox2_og = catalogue["bbox2"].to_list() + bbox3_og = catalogue["bbox3"].to_list() + bbox4_og = catalogue["bbox4"].to_list() + + x1 = catalogue["x1"].to_numpy() # - 1 #- X0 + y1 = catalogue["y1"].to_numpy() # - 1#- X0 polygons = [] - for i, source in tqdm(enumerate(Birth),total=len(Birth),desc='Creating contours..',disable=not output): - + for i, source in tqdm( + enumerate(Birth), + total=len(Birth), + desc="Creating contours..", + disable=not output, + ): + if use_gpu == True: - - cropped_image_gpu = image_gpu[bbox1_og[i]-1:bbox3_og[i]+1,bbox2_og[i]-1:bbox4_og[i]+1] - + + cropped_image_gpu = image_gpu[ + bbox1_og[i] - 1 : bbox3_og[i] + 1, bbox2_og[i] - 1 : bbox4_og[i] + 1 + ] + try: - red_image, red_mask, xmin, xmax, ymin, ymax = large_mask_red_image_procc_GPU(Birth[i],Death[i], - x1[i]-bbox1_og[i]+1, - y1[i]-bbox2_og[i]+1, - cropped_image_gpu, - X0[i],Y0[i]) + red_image, red_mask, xmin, xmax, ymin, ymax = ( + large_mask_red_image_procc_GPU( + Birth[i], + Death[i], + x1[i] - bbox1_og[i] + 1, + y1[i] - bbox2_og[i] + 1, + cropped_image_gpu, + X0[i], + Y0[i], + ) + ) except: - - print('Error in GPU processing!') - print('Source ID: ',IDs[i]) - print('x1:',x1[i]-bbox1_og[i]+1) - print('y1:',y1[i]-bbox2_og[i]+1) - print('Birth:',Birth[i]) - print('Death:',Death[i]) - #print('Cropped_image',cropped_image_gpu.get()) - + + print("Error in GPU processing!") + print("Source ID: ", IDs[i]) + print("x1:", x1[i] - bbox1_og[i] + 1) + print("y1:", y1[i] - bbox2_og[i] + 1) + print("Birth:", Birth[i]) + print("Death:", Death[i]) + # print('Cropped_image',cropped_image_gpu.get()) + red_mask = red_mask.astype(int) - + red_image = red_image.get() red_mask = red_mask.get() xmin = xmin.get() + bbox2_og[i] - xmax = xmax.get() + bbox2_og[i] + xmax = xmax.get() + bbox2_og[i] ymin = ymin.get() + bbox1_og[i] ymax = ymax.get() + bbox1_og[i] - + else: - cropped_image = image[bbox1_og[i]-1:bbox3_og[i]+1,bbox2_og[i]-1:bbox4_og[i]+1] + cropped_image = image[ + bbox1_og[i] - 1 : bbox3_og[i] + 1, bbox2_og[i] - 1 : bbox4_og[i] + 1 + ] try: - red_image, red_mask, xmin, xmax, ymin, ymax = large_mask_red_image_procc_CPU(Birth[i], - Death[i], - int(x1[i])-int(bbox1_og[i])+1, - int(y1[i])-int(bbox2_og[i])+1, - cropped_image) + red_image, red_mask, xmin, xmax, ymin, ymax = ( + large_mask_red_image_procc_CPU( + Birth[i], + Death[i], + int(x1[i]) - int(bbox1_og[i]) + 1, + int(y1[i]) - int(bbox2_og[i]) + 1, + cropped_image, + ) + ) except: print("Error in CPU processing!") - print('Source ID: ',IDs[i]) - print('x1:',x1[i]-bbox1_og[i]+1) - print('y1:',y1[i]-bbox2_og[i]+1) - print('Birth:',Birth[i]) - print('Death:',Death[i]) - + print("Source ID: ", IDs[i]) + print("x1:", x1[i] - bbox1_og[i] + 1) + print("y1:", y1[i] - bbox2_og[i] + 1) + print("Birth:", Birth[i]) + print("Death:", Death[i]) + red_mask = red_mask.astype(int) xmin = xmin + bbox2_og[i] xmax = xmax + bbox2_og[i] ymin = ymin + bbox1_og[i] ymax = ymax + bbox1_og[i] - - - contour = utils._get_polygons_in_bbox(xmin,xmax,ymin,ymax,x1[i],y1[i],Birth[i],Death[i],red_mask,0,0) + contour = utils._get_polygons_in_bbox( + xmin, xmax, ymin, ymax, x1[i], y1[i], Birth[i], Death[i], red_mask, 0, 0 + ) polygons.append(contour) - - catalogue['contour'] = polygons - return catalogue \ No newline at end of file + + catalogue["contour"] = polygons + return catalogue