From d8ebea6a1836df9fee8f14f0ffcc5227040af5af Mon Sep 17 00:00:00 2001 From: VChristiaens Date: Wed, 18 Jan 2023 20:23:26 +0100 Subject: [PATCH 01/14] Clarified verbose printed text for bad pixel clump algorithm --- vip_hci/preproc/badpixremoval.py | 17 ++++++++++++----- 1 file changed, 12 insertions(+), 5 deletions(-) diff --git a/vip_hci/preproc/badpixremoval.py b/vip_hci/preproc/badpixremoval.py index 837227d3..d7932819 100644 --- a/vip_hci/preproc/badpixremoval.py +++ b/vip_hci/preproc/badpixremoval.py @@ -740,11 +740,18 @@ def bp_removal_2d(array_corr, cy, cx, fwhm, sig, protect_mask, bpm_mask_ori, while nbpix_tbc > 0 and nit < max_nit: nit = nit+1 if verbose: - print("Iteration {}: {} bpix in total, {} to be " - "corrected".format(nit, nbpix_tot, nbpix_tbc)) - array_corr = sigma_filter(array_corr, bpix_map, neighbor_box=neighbor_box, - min_neighbors=nneig, half_res_y=half_res_y, - verbose=verbose) + msg = "Iteration {}: {} bad pixels identified".format(nit, + nbpix_tot) + if bpm_mask_ori is not None: + nbpix_ori = np.sum(bpm_mask_ori) + msg += " ({} new ones)".format(nbpix_tot-nbpix_ori) + if protect_mask: + msg += ", {} to be corrected".format(nbpix_tbc) + print(msg) + array_corr = sigma_filter(array_corr, bpix_map, + neighbor_box=neighbor_box, + min_neighbors=nneig, half_res_y=half_res_y, + verbose=verbose) bp = clip_array(array_corr, sig, sig, bpm_mask_ori, out_good=False, neighbor=True, num_neighbor=neighbor_box, mad=mad, half_res_y=half_res_y) From db70e99e0a4eed6a1aa25c75a483bfb0f8b44e96 Mon Sep 17 00:00:00 2001 From: VChristiaens Date: Wed, 18 Jan 2023 20:24:26 +0100 Subject: [PATCH 02/14] Added option to fully pad annulus with apertures, in function identifying independent apertures at a given radius --- vip_hci/metrics/snr_source.py | 36 ++++++++++++++++++++++++++++++++++- 1 file changed, 35 insertions(+), 1 deletion(-) diff --git a/vip_hci/metrics/snr_source.py b/vip_hci/metrics/snr_source.py index a9372b41..59c97e39 100644 --- a/vip_hci/metrics/snr_source.py +++ b/vip_hci/metrics/snr_source.py @@ -218,8 +218,38 @@ def _snr_approx(array, source_xy, fwhm, centery, centerx): return sourcey, sourcex, snr_value def indep_ap_centers(array, source_xy, fwhm, exclude_negative_lobes=False, - exclude_theta_range=None): + exclude_theta_range=None, no_gap=False): + """ + Define independent aperture centers at a given radial separation, starting + from a test location provided with source_xy. + + Parameters + ---------- + array : numpy ndarray, 2d + Frame in which the apertures will be defined (its dimensions are used). + source_xy : tuple of floats + X and Y coordinates of the planet or test speckle. + fwhm : float + Size in pixels of the FWHM, corresponding to the diameter of the + non-overlapping apertures. + exclude_negative_lobes : bool, opt + Whether to include the adjacent aperture lobes to the tested location + or not. Can be set to True if the image shows significant neg lobes. + exclude_theta_range : tuple of 2 floats or None, opt + If provided, range of trigonometric angles in deg (measured from + positive x axis), to be avoided for apertures used for noise estimation. + WARNING: this is to be used wisely, e.g. only if a known authentic + circumstellar signal is biasing the SNR estimate. + no_gap: bool, opt + Whether an overlapping aperture is defined between the first and last + non-overlapping aperture (at the end of making a full circle), in order + to leave no gap. False by default. + Returns + ------- + (yy, xx) : tuple of 2 numpy ndarray + Tuple containing y and x coordinates of the apertures + """ sourcex, sourcey = source_xy centery, centerx = frame_center(array) sep = dist(centery, centerx, float(sourcey), float(sourcex)) @@ -244,6 +274,10 @@ def indep_ap_centers(array, source_xy, fwhm, exclude_negative_lobes=False, angle = np.arcsin(fwhm / 2. / sep) * 2 number_apertures = int(np.floor(2 * np.pi / angle)) + if no_gap: + # if requested, add an (overlapping) aperture to avoid a gap + number_apertures += 1 + yy = [] xx = [] yy_all = np.zeros(number_apertures) From 8d23f6b618514586c758f65eaa8382752c65247a Mon Sep 17 00:00:00 2001 From: VChristiaens Date: Tue, 31 Jan 2023 16:03:25 +0100 Subject: [PATCH 03/14] Adapted try/except statements re: shared_memory module import to avoid error in Python 3.7 automated tests --- vip_hci/preproc/badpixremoval.py | 13 +++++++------ vip_hci/preproc/cosmetics.py | 9 +++++---- 2 files changed, 12 insertions(+), 10 deletions(-) diff --git a/vip_hci/preproc/badpixremoval.py b/vip_hci/preproc/badpixremoval.py index 82b28756..cc8fdbc1 100644 --- a/vip_hci/preproc/badpixremoval.py +++ b/vip_hci/preproc/badpixremoval.py @@ -37,6 +37,7 @@ from multiprocessing import Process import multiprocessing from multiprocessing import set_start_method +shared_mem = True try: from multiprocessing import shared_memory except ImportError: @@ -45,9 +46,9 @@ print('Trying to import shared_memory directly(for python 3.7)') import shared_memory except ModuleNotFoundError: - print('Use shared_memory on python 3.7 to activate') - print('multiprocessing on badpixels using..') - print('pip install shared-memory38') + shared_mem = False + print("WARNING: multiprocessing unavailable for bad pixel correction.") + print('Either pip install shared-memory38, or upgrade to python>=3.8') import warnings try: @@ -283,7 +284,7 @@ def cube_fix_badpix_isolated(array, bpm_mask=None, correct_only=False, if bpm_mask.ndim == 2: bpm_mask = [bpm_mask]*n_frames bpm_mask = np.array(bpm_mask) - if nproc==1: + if nproc==1 or not shared_mem: for i in Progressbar(range(n_frames), desc="processing frames"): if bpm_mask is not None: bpm_mask_tmp = bpm_mask[i] @@ -886,7 +887,7 @@ def bp_removal_2d(array_corr, cy, cx, fwhm, sig, protect_mask, bpm_mask_ori, cx = [cx]*n_z if isinstance(fwhm, (float, int)): fwhm = [fwhm]*n_z - if nproc==1: + if nproc==1 or not shared_mem: bpix_map_cumul = np.zeros_like(array_corr) for i in range(n_z): if verbose: @@ -942,7 +943,7 @@ def _mp_clump_slow(args): neighbor_box = max(3, fwhm_round) # to not replace a companion nneig = sum(np.arange(3, neighbor_box+2, 2)) - if nproc==1: + if nproc==1 or not shared_mem: for i in range(n_z): if verbose: print('Using serial approach') diff --git a/vip_hci/preproc/cosmetics.py b/vip_hci/preproc/cosmetics.py index 7bd0d82c..791419ed 100644 --- a/vip_hci/preproc/cosmetics.py +++ b/vip_hci/preproc/cosmetics.py @@ -24,6 +24,7 @@ from multiprocessing import Process import multiprocessing from multiprocessing import set_start_method +shared_mem = True try: from multiprocessing import shared_memory except ImportError: @@ -32,9 +33,9 @@ print('Trying to import shared_memory directly(for python 3.7)') import shared_memory except ModuleNotFoundError: - print('Use shared_memory on python 3.7 to activate') - print('multiprocessing on badpixels using..') - print('pip install shared-memory38') + shared_mem = False + print("WARNING: multiprocessing unavailable for bad pixel correction.") + print('Either pip install shared-memory38, or upgrade to python>=3.8') def cube_crop_frames(array, size, xy=None, force=False, verbose=True, @@ -356,7 +357,7 @@ def cube_correct_nan(cube, neighbor_box=3, min_neighbors=3, verbose=False, if nproc is None: nproc = cpu_count()//2 n_z = obj_tmp.shape[0] - if nproc == 1: + if nproc == 1 or not shared_mem: for zz in range(n_z): obj_tmp[zz], nnanpix = nan_corr_2d(obj_tmp[zz], neighbor_box, min_neighbors, half_res_y, From 7f37f2744a0718cb13ddab5cc465f7291dbda9d1 Mon Sep 17 00:00:00 2001 From: VChristiaens Date: Wed, 1 Feb 2023 20:31:56 +0100 Subject: [PATCH 04/14] Added option for exclusion mask to bad pixel identification/correction routines --- vip_hci/preproc/badpixremoval.py | 432 +++++++++++++++++++++---------- 1 file changed, 293 insertions(+), 139 deletions(-) diff --git a/vip_hci/preproc/badpixremoval.py b/vip_hci/preproc/badpixremoval.py index cc8fdbc1..2e5d464e 100644 --- a/vip_hci/preproc/badpixremoval.py +++ b/vip_hci/preproc/badpixremoval.py @@ -14,6 +14,7 @@ """ +import warnings __author__ = 'V. Christiaens, Carlos Alberto Gomez Gonzalez' __all__ = ['frame_fix_badpix_isolated', 'cube_fix_badpix_isolated', @@ -39,18 +40,17 @@ from multiprocessing import set_start_method shared_mem = True try: - from multiprocessing import shared_memory + from multiprocessing import shared_memory except ImportError: - print('Failed to import shared_memory from multiprocessing') - try: - print('Trying to import shared_memory directly(for python 3.7)') - import shared_memory - except ModuleNotFoundError: - shared_mem = False - print("WARNING: multiprocessing unavailable for bad pixel correction.") - print('Either pip install shared-memory38, or upgrade to python>=3.8') + print('Failed to import shared_memory from multiprocessing') + try: + print('Trying to import shared_memory directly(for python 3.7)') + import shared_memory + except ModuleNotFoundError: + shared_mem = False + print("WARNING: multiprocessing unavailable for bad pixel correction.") + print('Either pip install shared-memory38, or upgrade to python>=3.8') -import warnings try: from numba import njit no_numba = False @@ -62,8 +62,8 @@ def frame_fix_badpix_isolated(array, bpm_mask=None, correct_only=False, sigma_clip=3, num_neig=5, size=5, protect_mask=0, - cxy=None, mad=False, ignore_nan=True, - verbose=True, full_output=False): + excl_mask=None, cxy=None, mad=False, + ignore_nan=True, verbose=True, full_output=False): """ Corrects the bad pixels, marked in the bad pixel mask. The bad pixel is replaced by the median of the adjacent pixels. This function is very fast but works only with isolated (sparse) pixels. @@ -73,9 +73,9 @@ def frame_fix_badpix_isolated(array, bpm_mask=None, correct_only=False, array : numpy ndarray Input 2d array. bpm_mask : numpy ndarray, optional - Input bad pixel map. Zeros frame where the bad pixels have a value of - 1. If None is provided a bad pixel map will be created using - sigma clip statistics. + Input bad pixel map. Should have same size as array. Binary map where + the bad pixels have a value of 1. If None is provided a bad pixel map + will be created using sigma clip statistics. correct_only : bool, opt If True and bpix_map is provided, will only correct for provided bad pixels. Else, the algorithm will determine (more) bad pixels. @@ -93,6 +93,11 @@ def frame_fix_badpix_isolated(array, bpm_mask=None, correct_only=False, If larger than 0, radius of a circular aperture (at the center of the frames) in which no bad pixels will be identified. This can be useful to protect the star and vicinity. + excl_mask : numpy ndarray, optional + Binary mask with 1 in areas that should not be considered as good + neighbouring pixels during the identification of bad pixels. These + should not be considered as bad pixels to be corrected neither (i.e. + different to bpm_mask). cxy: None or tuple If protect_mask is True, this is the location of the star centroid in the images. If None, assumes the star is already centered. If a tuple, @@ -126,7 +131,15 @@ def frame_fix_badpix_isolated(array, bpm_mask=None, correct_only=False, raise ValueError(msg) if bpm_mask is not None: + msg = "Input bad pixel mask should have same shape as array\n" + assert bpm_mask.shape == array.shape, msg bpm_mask = bpm_mask.astype('bool') + if excl_mask is None: + excl_mask = np.zeros(array.shape, dtype=bool) + else: + msg = "Input exclusion mask should have same shape as array\n" + assert excl_mask.shape == array.shape, msg + ind_excl = np.where(excl_mask) if verbose: start = time_ini() @@ -143,6 +156,10 @@ def frame_fix_badpix_isolated(array, bpm_mask=None, correct_only=False, cx, cy = cxy if bpm_mask is None or not correct_only: + if bpm_mask is None: + bpm_mask = np.zeros(array.shape, dtype=bool) + bpm_mask = bpm_mask+excl_mask + bpm_mask[np.where(bpm_mask > 1)] = 1 ori_nan_mask = np.where(np.isnan(frame)) ind = clip_array(frame, sigma_clip, sigma_clip, bpm_mask, neighbor=neigh, num_neighbor=num_neig, mad=mad) @@ -153,6 +170,7 @@ def frame_fix_badpix_isolated(array, bpm_mask=None, correct_only=False, if protect_mask: cir = disk((cy, cx), protect_mask, shape=bpm_mask.shape) bpm_mask[cir] = 0 + bpm_mask[ind_excl] = 0 bpm_mask = bpm_mask.astype('bool') smoothed = median_filter(frame, size, mode='mirror') @@ -173,9 +191,10 @@ def frame_fix_badpix_isolated(array, bpm_mask=None, correct_only=False, def cube_fix_badpix_isolated(array, bpm_mask=None, correct_only=False, sigma_clip=3, num_neig=5, size=5, - frame_by_frame=False, protect_mask=0, cxy=None, - mad=False, ignore_nan=True, verbose=True, - full_output=False, nproc=1): + frame_by_frame=False, protect_mask=0, + excl_mask=None, cxy=None, mad=False, + ignore_nan=True, verbose=True, full_output=False, + nproc=1): """ Corrects the bad pixels, marked in the bad pixel mask. The bad pixel is replaced by the median of the adjacent pixels. This function is very fast but works only with isolated (sparse) pixels. @@ -208,6 +227,11 @@ def cube_fix_badpix_isolated(array, bpm_mask=None, correct_only=False, If larger than 0, radius of a circular aperture (at the center of the frames) in which no bad pixels will be identified. This can be useful to protect the star and vicinity. + excl_mask : numpy ndarray, optional + Binary mask with 1 in areas that should not be considered as good + neighbouring pixels during the identification of bad pixels. These + should not be considered as bad pixels to be corrected neither (i.e. + different to bpm_mask). cxy: None, tuple or 2d numpy ndarray If protect_mask is True, this is the location of the star centroid in the images. If None, assumes the star is already centered. If a tuple, @@ -226,17 +250,19 @@ def cube_fix_badpix_isolated(array, bpm_mask=None, correct_only=False, Whether to return as well the cube of bad pixel maps and the cube of defined annuli. nproc: int, optional - This feature is added following ADACS update. Refers to the number of processors - available for calculations. Choosing a number >1 enables multiprocessing for the - correction of frames. This happens only when frame_by_frame=True. + This feature is added following ADACS update. Refers to the number of + processors available for calculations. Choosing a number >1 enables + multiprocessing for the correction of frames. This happens only when + ``frame_by_frame=True''. Return ------ array_out : numpy ndarray Cube with bad pixels corrected. bpm_mask: 2d or 3d array [if full_output is True] - The bad pixel map or the cube of bpix maps + The bad pixel map or the cube of bad pixel maps """ + if array.ndim != 3: raise TypeError('Array is not a 3d array or cube') if size % 2 == 0: @@ -246,6 +272,8 @@ def cube_fix_badpix_isolated(array, bpm_mask=None, correct_only=False, raise ValueError(msg) if bpm_mask is not None: + msg = "Input bad pixel mask should have same last 2 dims as array\n" + assert bpm_mask.shape[-2:] == array.shape[-2:], msg bpm_mask = bpm_mask.astype('bool') if verbose: @@ -284,76 +312,129 @@ def cube_fix_badpix_isolated(array, bpm_mask=None, correct_only=False, if bpm_mask.ndim == 2: bpm_mask = [bpm_mask]*n_frames bpm_mask = np.array(bpm_mask) - if nproc==1 or not shared_mem: + if nproc == 1 or not shared_mem: for i in Progressbar(range(n_frames), desc="processing frames"): if bpm_mask is not None: bpm_mask_tmp = bpm_mask[i] else: bpm_mask_tmp = None + if excl_mask is not None: + excl_mask_tmp = excl_mask[i] + else: + excl_mask_tmp = None res = frame_fix_badpix_isolated(array[i], bpm_mask=bpm_mask_tmp, sigma_clip=sigma_clip, num_neig=num_neig, size=size, protect_mask=protect_mask, - verbose=False, cxy=(cx[i], cy[i]), + excl_mask=excl_mask_tmp, + verbose=False, + cxy=(cx[i], cy[i]), ignore_nan=ignore_nan, full_output=True) array_out[i] = res[0] final_bpm[i] = res[1] else: print("Processing using ADACS' multiprocessing approach...") - #dummy calling the function to create cached version of the code prior to forking + # dummy calling the function to create cached version of the code prior to forking if bpm_mask is not None: bpm_mask_dum = bpm_mask[0] else: bpm_mask_dum = None - #point of dummy call - frame_fix_badpix_isolated(array[0], bpm_mask=bpm_mask_dum, sigma_clip=sigma_clip, num_neig=num_neig, size=size, protect_mask=protect_mask, verbose=False, cxy=(cx[0], cy[0]), ignore_nan=ignore_nan, full_output=False) - # multiprocessing included only in the frame-by-frame branch of the if statement above. + if excl_mask is not None: + excl_mask_dum = excl_mask[0] + else: + excl_mask_dum = None + # point of dummy call + frame_fix_badpix_isolated(array[0], bpm_mask=bpm_mask_dum, + sigma_clip=sigma_clip, num_neig=num_neig, + size=size, protect_mask=protect_mask, + excl_mask=excl_mask_dum, verbose=False, + cxy=(cx[0], cy[0]), ignore_nan=ignore_nan, + full_output=False) + # multiprocessing included only in the frame-by-frame branch of the + # if statement above. # creating shared memory buffer for the cube (array) - shm_array_out = shared_memory.SharedMemory(create=True, size=array.nbytes) - # creating a shared array_out version that is the shm_array_out buffer above. - shared_array_out = np.ndarray(array.shape, dtype=array.dtype, buffer=shm_array_out.buf) + shm_arr = shared_memory.SharedMemory(create=True, size=array.nbytes) + # creating a shared array_out version that is the shm_array_out + # buffer above. + sh_arr = np.ndarray(array.shape, dtype=array.dtype, + buffer=shm_arr.buf) # creating shared memory buffer for the final bad pixel mask cube. - shm_final_bpm = shared_memory.SharedMemory(create=True, size=final_bpm.nbytes) - # creating a shared final_bpm version that is in the shm_final_bpm buffer above. - shared_final_bpm = np.ndarray(final_bpm.shape, dtype=final_bpm.dtype, buffer=shm_final_bpm.buf) - - #function that calls frame_fix_badpix_isolated using the similar arguments as in if nproc==1 branch above. - def mp_clean_isolated (j,frame, bpm_mask=None, sigma_clip=3, num_neig=5, size=5, protect_mask=0, verbose=False, cxy=None, ignore_nan=True, full_output=True): - shared_array_out[j], shared_final_bpm[j] = frame_fix_badpix_isolated(frame, bpm_mask=bpm_mask, sigma_clip=sigma_clip, num_neig=num_neig, size=size, protect_mask=protect_mask, verbose=verbose, cxy=cxy, ignore_nan=ignore_nan, full_output=full_output) - #function that unwraps the arguments and passes them to mp_clean_isolated. + shm_fbpm = shared_memory.SharedMemory(create=True, + size=final_bpm.nbytes) + # creating a shared final_bpm version that is in the shm_final_bpm + # buffer above. + sh_fbpm = np.ndarray(final_bpm.shape, dtype=final_bpm.dtype, + buffer=shm_fbpm.buf) + + # function that calls frame_fix_badpix_isolated using the similar + # arguments as in if nproc==1 branch above. + def mp_clean_isolated(j, frame, bpm_mask=None, sigma_clip=3, + num_neig=5, size=5, protect_mask=0, + excl_mask=None, verbose=False, cxy=None, + ignore_nan=True, full_output=True): + sh_res = frame_fix_badpix_isolated(frame, bpm_mask, + sigma_clip=sigma_clip, + num_neig=num_neig, size=size, + protect_mask=protect_mask, + excl_mask=excl_mask, + verbose=verbose, cxy=cxy, + ignore_nan=ignore_nan, + full_output=full_output) + sh_arr[j], sh_fbpm[j] = sh_res + # function that unwraps the arguments and passes them to + # mp_clean_isolated. global _mp_clean_isolated - def _mp_clean_isolated (args): - pargs=args[0:2] - kwargs=args[2] + + def _mp_clean_isolated(args): + pargs = args[0:2] + kwargs = args[2] mp_clean_isolated(*pargs, **kwargs) - - context=multiprocessing.get_context('fork') - pool=context.Pool(processes=nproc, maxtasksperchild=1) - args=[] + context = multiprocessing.get_context('fork') + pool = context.Pool(processes=nproc, maxtasksperchild=1) + + args = [] for j in range(n_frames): if bpm_mask is not None: bpm_mask_tmp = bpm_mask[j] else: bpm_mask_tmp = None - dict_kwargs={'bpm_mask' : bpm_mask_tmp, 'sigma_clip': sigma_clip, 'num_neig': num_neig, 'size' : size, 'protect_mask': protect_mask, 'cxy' : (cx[j], cy[j]), 'ignore_nan': ignore_nan} + if excl_mask is not None: + excl_mask_tmp = excl_mask[j] + else: + excl_mask_tmp = None + dict_kwargs = {'bpm_mask': bpm_mask_tmp, + 'sigma_clip': sigma_clip, 'num_neig': num_neig, + 'size': size, 'protect_mask': protect_mask, + 'excl_mask': excl_mask_tmp, 'cxy': (cx[j], cy[j]), + 'ignore_nan': ignore_nan} args.append([j, array[j], dict_kwargs]) - + try: - pool.map_async(_mp_clean_isolated, args, chunksize=1 ).get(timeout=10_000_000) + pool.map_async(_mp_clean_isolated, args, + chunksize=1).get(timeout=10_000_000) finally: pool.close() pool.join() - array_out[:]=shared_array_out[:] - final_bpm[:]=shared_final_bpm[:] - shm_array_out.close() - shm_array_out.unlink() - shm_final_bpm.close() - shm_final_bpm.unlink() + array_out[:] = sh_arr[:] + final_bpm[:] = sh_fbpm[:] + shm_arr.close() + shm_arr.unlink() + shm_fbpm.close() + shm_fbpm.unlink() count_bp = np.sum(final_bpm) else: + if excl_mask is None: + excl_mask = np.zeros(array.shape, dtype=bool) + elif excl_mask is None: + msg = "Input exclusion mask should have same shape as array\n" + assert excl_mask.shape == array.shape, msg + ind_excl = np.where(excl_mask) if bpm_mask is None or not correct_only: + if bpm_mask is None: + bpm_mask = np.zeros(array.shape, dtype=bool) + bpm_mask = bpm_mask+excl_mask ori_nan_mask = np.where(np.isnan(np.nanmean(array, axis=0))) ind = clip_array(np.nanmean(array, axis=0), sigma_clip, sigma_clip, bpm_mask, neighbor=neigh, num_neighbor=num_neig, @@ -365,6 +446,7 @@ def _mp_clean_isolated (args): if protect_mask: cir = disk((cy, cx), protect_mask, shape=final_bpm.shape) final_bpm[cir] = 0 + final_bpm[ind_excl] = 0 final_bpm = final_bpm.astype('bool') else: if bpm_mask.ndim == 3: @@ -631,10 +713,10 @@ def bp_removal_2d(array, cy, cx, fwhm, sig, protect_mask, r_in_std, cy, cx = frame_center(array) if ndims == 2: array_corr, bpix_map, ann_frame_cumul = bp_removal_2d(array, cy, cx, - fwhm, sig, - protect_mask, - r_in_std, r_out_std, - verbose) + fwhm, sig, + protect_mask, + r_in_std, r_out_std, + verbose) if ndims == 3: array_corr = array.copy() n_z = array.shape[0] @@ -661,8 +743,9 @@ def bp_removal_2d(array, cy, cx, fwhm, sig, protect_mask, r_in_std, def cube_fix_badpix_clump(array, bpm_mask=None, correct_only=False, cy=None, cx=None, fwhm=4., sig=4., protect_mask=0, - half_res_y=False, min_thr=None, max_nit=15, mad=True, - verbose=True, full_output=False, nproc=1): + excl_mask=None, half_res_y=False, min_thr=None, + max_nit=15, mad=True, verbose=True, full_output=False, + nproc=1): """ Function to identify and correct clumps of bad pixels. Very fast when a bad pixel map is provided. If a bad pixel map is not provided, the bad pixel @@ -676,9 +759,10 @@ def cube_fix_badpix_clump(array, bpm_mask=None, correct_only=False, cy=None, array : 3D or 2D array Input 3d cube or 2d image. bpm_mask: 3D or 2D array, opt - Input bad pixel array. Should have same dimensions as array. If not - provided, the algorithm will attempt to identify bad pixel clumps - automatically. + Input bad pixel array. If 2D and array is 3D: should have same last 2 + dimensions as array. If 3D, should have exact same dimensions as input + array. If not provided, the algorithm will attempt to identify bad pixel + clumps automatically. correct_only : bool, opt If True and bpix_map is provided, will only correct for provided bad pixels. Else, the algorithm will determine (more) bad pixels. @@ -699,6 +783,11 @@ def cube_fix_badpix_clump(array, bpm_mask=None, correct_only=False, cy=None, If larger than 0, radius of a circular aperture (at the center of the frames) in which no bad pixels will be identified. This can be useful to protect the star and vicinity. + excl_mask : numpy ndarray, optional + Binary mask with same dimensions as array, with 1 in areas that should + not be considered as good neighbouring pixels during the identification + of bad pixels. These should not be considered as bad pixels to be + corrected neither (i.e. different to bpm_mask). half_res_y: bool, {True,False}, optional Whether the input data has only half the angular resolution vertically compared to horizontally (e.g. the case of SINFONI data); in other words @@ -744,15 +833,21 @@ def cube_fix_badpix_clump(array, bpm_mask=None, correct_only=False, cy=None, assert ndims == 2 or ndims == 3, "Object is not two or three dimensional.\n" if bpm_mask is not None: - if bpm_mask.shape[-2:] != array.shape[-2:]: - raise TypeError("Bad pixel map has wrong y/x dimensions.") + msg = "Input bad pixel mask should have same last 2 dims as array\n" + assert bpm_mask.shape[-2:] == array.shape[-2:], msg + bpm_mask = bpm_mask.astype('bool') if correct_only and bpm_mask is None: msg = "Bad pixel map should be provided if correct_only is True." raise ValueError(msg) def bp_removal_2d(array_corr, cy, cx, fwhm, sig, protect_mask, bpm_mask_ori, - min_thr, half_res_y, mad, verbose): + excl_mask, min_thr, half_res_y, mad, verbose): + + msg = "Input exclusion mask should have same shape as array\n" + assert excl_mask.shape == array.shape, msg + ind_excl = np.where(excl_mask) + n_x = array_corr.shape[1] n_y = array_corr.shape[0] @@ -788,7 +883,9 @@ def bp_removal_2d(array_corr, cy, cx, fwhm, sig, protect_mask, bpm_mask_ori, circl_new = [] # 3/ Create a bad pixel map, by detecting them with clip_array - bp = clip_array(array_corr, sig, sig, bpm_mask_ori, out_good=False, + bpm_mask = bpm_mask_ori+excl_mask + bpm_mask[np.where(bpm_mask > 1)] = 1 + bp = clip_array(array_corr, sig, sig, bpm_mask, out_good=False, neighbor=True, num_neighbor=neighbor_box, mad=mad, half_res_y=half_res_y) bpix_map = np.zeros_like(array_corr) @@ -808,6 +905,7 @@ def bp_removal_2d(array_corr, cy, cx, fwhm, sig, protect_mask, bpm_mask_ori, bpix_map[np.where(cond1 & cond2)] = 0 nbpix_tot = int(np.sum(bpix_map)) bpix_map[circl_new] = 0 + bpix_map[ind_excl] = 0 nbpix_tbc = int(np.sum(bpix_map)) bpix_map_cumul = np.zeros_like(bpix_map) bpix_map_cumul[:] = bpix_map[:] @@ -829,7 +927,9 @@ def bp_removal_2d(array_corr, cy, cx, fwhm, sig, protect_mask, bpm_mask_ori, neighbor_box=neighbor_box, min_neighbors=nneig, half_res_y=half_res_y, verbose=verbose) - bp = clip_array(array_corr, sig, sig, bpm_mask_ori, out_good=False, + bpm_mask += bpix_map_cumul + bpm_mask[np.where(bpm_mask > 1)] = 1 + bp = clip_array(array_corr, sig, sig, bpm_mask, out_good=False, neighbor=True, num_neighbor=neighbor_box, mad=mad, half_res_y=half_res_y) bpix_map = np.zeros_like(array_corr) @@ -840,6 +940,7 @@ def bp_removal_2d(array_corr, cy, cx, fwhm, sig, protect_mask, bpm_mask_ori, bpix_map[np.where(cond1 & cond2)] = 0 nbpix_tot = int(np.sum(bpix_map)) bpix_map[circl_new] = 0 + bpix_map[ind_excl] = 0 nbpix_tbc = int(np.sum(bpix_map)) bpix_map_cumul = bpix_map_cumul+bpix_map @@ -864,20 +965,29 @@ def bp_removal_2d(array_corr, cy, cx, fwhm, sig, protect_mask, bpm_mask_ori, cy, cx = frame_center(array) array_corr, bpix_map_cumul = bp_removal_2d(array_corr, cy, cx, fwhm, sig, protect_mask, - bpm_mask, min_thr, - half_res_y, mad, verbose) + bpm_mask, excl_mask, + min_thr, half_res_y, mad, + verbose) else: fwhm_round = int(round(fwhm)) fwhm_round = fwhm_round+1-(fwhm_round % 2) # make it odd neighbor_box = max(3, fwhm_round) # to not replace a companion nneig = sum(np.arange(3, neighbor_box+2, 2)) array_corr = sigma_filter(array_corr, bpm_mask, neighbor_box, nneig, - half_res_y, verbose) + half_res_y, verbose) bpix_map_cumul = bpm_mask if ndims == 3: n_z = array_corr.shape[0] if bpm_mask is None or not correct_only: + if bpm_mask is None: + bpm_mask = np.zeros(array_corr.shape, dtype=bool) + elif bpm_mask.ndim == 2: + bpm_mask = np.array([bpm_mask]*n_z, dtype=bool) + if excl_mask is None: + excl_mask = np.zeros(array_corr.shape, dtype=bool) + elif excl_mask.ndim == 2: + excl_mask = np.array([excl_mask]*n_z, dtype=bool) if cy is None or cx is None: cy, cx = frame_center(array) cy = [cy]*n_z @@ -887,48 +997,62 @@ def bp_removal_2d(array_corr, cy, cx, fwhm, sig, protect_mask, bpm_mask_ori, cx = [cx]*n_z if isinstance(fwhm, (float, int)): fwhm = [fwhm]*n_z - if nproc==1 or not shared_mem: + if nproc == 1 or not shared_mem: bpix_map_cumul = np.zeros_like(array_corr) for i in range(n_z): if verbose: print('************Frame # ', i, ' *************') - array_corr[i], bpix_map_cumul[i] = bp_removal_2d(array_corr[i], cy[i], - cx[i], fwhm[i], - sig, protect_mask, - bpm_mask, min_thr, - half_res_y, mad, - verbose) + res = bp_removal_2d(array_corr[i], cy[i], cx[i], fwhm[i], + sig, protect_mask, bpm_mask[i], + excl_mask[i], min_thr, half_res_y, mad, + verbose) + array_corr[i], bpix_map_cumul[i] = res else: - msg="Cleaning frames using ADACS' multiprocessing appraoch" - print(msg) - #creating shared memory buffer space for the image cube. - shm_clump= shared_memory.SharedMemory(create=True, size=array_corr.nbytes) - obj_tmp_shared_clump = np.ndarray(array_corr.shape, dtype=array_corr.dtype, buffer=shm_clump.buf) - #creating shared memory buffer space for the bad pixel cube. - shm_clump_bpix= shared_memory.SharedMemory(create=True, size=array_corr.nbytes) - #works with dtype=obj_tmp.dtype but not dtype=int - bpix_map_cumul_shared= np.ndarray(array_corr.shape, dtype=array_corr.dtype, buffer=shm_clump_bpix.buf) - def mp_clump_slow(j, array_corr, cy, cx, fwhm, sig, protect_mask, bpm_mask, min_thr, half_res_y, mad, verbose): - obj_tmp_shared_clump[j], bpix_map_cumul_shared[j] = bp_removal_2d(array_corr, cy, cx, fwhm, sig, protect_mask, bpm_mask, min_thr, half_res_y, mad, verbose) - + msg = "Cleaning frames using ADACS' multiprocessing approach" + print(msg) + # creating shared memory buffer space for the image cube. + shm_clump = shared_memory.SharedMemory(create=True, + size=array_corr.nbytes) + obj_tmp_shared_clump = np.ndarray(array_corr.shape, + dtype=array_corr.dtype, + buffer=shm_clump.buf) + # creating shared memory buffer space for the bad pixel cube. + shm_clump_bpix = shared_memory.SharedMemory(create=True, + size=array_corr.nbytes) + # works with dtype=obj_tmp.dtype but not dtype=int + bpix_map_cumul_shared = np.ndarray(array_corr.shape, + dtype=array_corr.dtype, + buffer=shm_clump_bpix.buf) + + def mp_clump_slow(j, array_corr, cy, cx, fwhm, sig, protect_mask, + bpm_mask, excl_mask, min_thr, half_res_y, mad, + verbose): + res = bp_removal_2d(array_corr, cy, cx, fwhm, sig, + protect_mask, bpm_mask, excl_mask, + min_thr, half_res_y, mad, verbose) + obj_tmp_shared_clump[j], bpix_map_cumul_shared[j] = res global _mp_clump_slow - + def _mp_clump_slow(args): mp_clump_slow(*args) - - context=multiprocessing.get_context('fork') - pool=context.Pool(processes=nproc, maxtasksperchild=1) - args=[] + + context = multiprocessing.get_context('fork') + pool = context.Pool(processes=nproc, maxtasksperchild=1) + args = [] for i in range(n_z): - args.append([i,array_corr[i], cy[i], cx[i], fwhm[i], sig, protect_mask, bpm_mask, min_thr, half_res_y, mad, verbose ]) + args.append([i, array_corr[i], cy[i], cx[i], fwhm[i], sig, + protect_mask, bpm_mask[i], excl_mask[i], min_thr, + half_res_y, mad, verbose]) try: - pool.map_async(_mp_clump_slow, args, chunksize=1 ).get(timeout=10_000_000) - finally: + pool.map_async(_mp_clump_slow, args, chunksize=1).get( + timeout=10_000_000) + finally: pool.close() pool.join() - bpix_map_cumul = np.zeros_like(array_corr, dtype=array_corr.dtype) - bpix_map_cumul[:]=bpix_map_cumul_shared[:] - array_corr[:]=obj_tmp_shared_clump[:] + bpix_map_cumul = np.zeros_like(array_corr, + dtype=array_corr.dtype) + bpix_map_cumul[:] = bpix_map_cumul_shared[:] + array_corr[:] = obj_tmp_shared_clump[:] shm_clump.close() shm_clump.unlink() shm_clump_bpix.close() @@ -943,7 +1067,7 @@ def _mp_clump_slow(args): neighbor_box = max(3, fwhm_round) # to not replace a companion nneig = sum(np.arange(3, neighbor_box+2, 2)) - if nproc==1 or not shared_mem: + if nproc == 1 or not shared_mem: for i in range(n_z): if verbose: print('Using serial approach') @@ -952,51 +1076,62 @@ def _mp_clump_slow(args): bpm = bpm_mask[i] else: bpm = bpm_mask - array_corr[i] = sigma_filter(array_corr[i], bpm, neighbor_box, - nneig, half_res_y, verbose) + array_corr[i] = sigma_filter(array_corr[i], bpm, + neighbor_box, nneig, half_res_y, + verbose) else: - msg="Cleaning frames using ADACS' multiprocessing appraoch" + msg = "Cleaning frames using ADACS' multiprocessing approach" print(msg) - #dummy calling sigma_filter function to create a cached version of the numba function + # dummy calling sigma_filter function to create a cached version of the numba function if bpm_mask.ndim == 3: dummy_bpm = bpm_mask[0] - else: + else: dummy_bpm = bpm_mask - #Actual dummy call is here. - sigma_filter(array_corr[0], dummy_bpm, neighbor_box, nneig, half_res_y, verbose) - #creating shared memory that each process writes into. - shm_clump = shared_memory.SharedMemory(create=True, size=array_corr.nbytes) - #creating an array that uses shared memory buffer and has the properties of array_corr. - obj_tmp_shared_clump = np.ndarray(array_corr.shape, dtype=array_corr.dtype, buffer=shm_clump.buf) - #function that is called repeatedly by each process. - def mp_clean_clump(j, array_corr, bpm, neighbor_box, nneig, half_res_y, verbose): - obj_tmp_shared_clump[j] = sigma_filter(array_corr, bpm, neighbor_box, nneig, half_res_y, verbose) - + # Actual dummy call is here. + sigma_filter(array_corr[0], dummy_bpm, + neighbor_box, nneig, half_res_y, verbose) + # creating shared memory that each process writes into. + shm_clump = shared_memory.SharedMemory( + create=True, size=array_corr.nbytes) + # creating an array that uses shared memory buffer and has the properties of array_corr. + obj_tmp_shared_clump = np.ndarray( + array_corr.shape, dtype=array_corr.dtype, buffer=shm_clump.buf) + # function that is called repeatedly by each process. + + def mp_clean_clump(j, array_corr, bpm, neighbor_box, nneig, + half_res_y, verbose): + obj_tmp_shared_clump[j] = sigma_filter(array_corr, bpm, + neighbor_box, nneig, + half_res_y, verbose) + global _mp_clean_clump - #function that converts the args into bite-sized pieces for mp_clean_clump. + # function that converts the args into bite-sized pieces for mp_clean_clump. + def _mp_clean_clump(args): mp_clean_clump(*args) - context=multiprocessing.get_context('fork') - pool=context.Pool(processes=nproc, maxtasksperchild=1) - args=[] + context = multiprocessing.get_context('fork') + pool = context.Pool(processes=nproc, maxtasksperchild=1) + args = [] for j in range(n_z): if bpm_mask.ndim == 3: bpm = bpm_mask[j] else: bpm = bpm_mask - args.append([j,array_corr[j], bpm, neighbor_box, nneig, half_res_y, verbose]) + args.append([j, array_corr[j], bpm, neighbor_box, + nneig, half_res_y, verbose]) try: - pool.map_async(_mp_clean_clump, args, chunksize=1 ).get(timeout=10_000_000) - finally: + pool.map_async(_mp_clean_clump, args, chunksize=1).get( + timeout=10_000_000) + finally: pool.close() pool.join() - array_corr[:]=obj_tmp_shared_clump[:] + array_corr[:] = obj_tmp_shared_clump[:] shm_clump.close() - shm_clump.unlink() + shm_clump.unlink() bpix_map_cumul = bpm_mask # make it a binary map - bpix_map_cumul[np.where(bpix_map_cumul>1)] = 1 + bpix_map_cumul[np.where(bpix_map_cumul > 1)] = 1 if full_output: return array_corr, bpix_map_cumul @@ -1218,14 +1353,15 @@ def _res_scaled_images(array, lbdas, fluxes, mask, cy, cx): return array_out -def cube_fix_badpix_interp(array, bpm_mask, mode='fft', fwhm=4., kernel_sz=None, - psf=None, half_res_y=False, nit=500, tol=1, nproc=1, - full_output=False, **kwargs): +def cube_fix_badpix_interp(array, bpm_mask, mode='fft', excl_mask=None, fwhm=4., + kernel_sz=None, psf=None, half_res_y=False, nit=500, + tol=1, nproc=1, full_output=False, **kwargs): """ Function to correct clumps of bad pixels by interpolation with either a user-defined kernel (through astropy.convolution) or through the FFT-based algorithm described in [AAC01]_. A bad pixel map must be provided (e.g. found with function `cube_fix_badpix_clump`). + Parameters ---------- array : 3D or 2D array @@ -1237,6 +1373,12 @@ def cube_fix_badpix_interp(array, bpm_mask, mode='fft', fwhm=4., kernel_sz=None, mode: str, optional {'fft', 'gauss', 'psf'} Can be either a 2D Gaussian ('gauss') or an input normalized PSF ('psf'). + excl_mask: 3D or 2D array, optional + [Only used if mode != 'fft'] Input exclusion mask array. Pixels in the + exclusion mask will neither be used for interpolation, nor replaced as + bad pixels. excl_mask should have same x,y dimensions as array. If 2D, + but input array is 3D, the same exclusion mask will be assumed for all + frames. fwhm: float, 1D array or tuple of 2 floats, opt If mode is 'gauss', the fwhm of the Gaussian. kernel_sz: int or None, optional @@ -1289,6 +1431,14 @@ def cube_fix_badpix_interp(array, bpm_mask, mode='fft', fwhm=4., kernel_sz=None, if bpm_mask.shape[-2:] != array.shape[-2:]: raise TypeError("Bad pixel map has wrong y/x dimensions.") + if excl_mask is None: + excl_mask = np.zeros(array.shape, dtype=bool) + elif excl_mask.ndim == 2 and array.ndim == 3: + nz = array.shape[0] + excl_mask = np.array([excl_mask]*nz, dtype=bool) + msg = "Input exclusion mask should have same shape as array\n" + assert excl_mask.shape[-2:] == array.shape[-2:], msg + if np.sum(bpm_mask) == 0: msg = "Warning: no bad pixel found in bad pixel map. " msg += "Returning input array as is." @@ -1320,18 +1470,22 @@ def squash_v(array): if ndims == 2: array_corr = squash_v(array) bpm_mask = squash_v(bpm_mask) + bpm_mask = squash_v(excl_mask) else: new_array_corr = [] new_bpm_mask = [] + new_excl_mask = [] for z in range(nz): new_array_corr.append(squash_v(array_corr[z])) new_bpm_mask.append(squash_v(bpm_mask[z])) + new_excl_mask.append(squash_v(excl_mask[z])) array_corr = np.array(new_array_corr) bpm_mask = np.array(new_bpm_mask) + excl_mask = np.array(new_excl_mask) if mode != 'fft': # first replace all bad pixels with NaNs - they will be interpolated - array_corr[np.where(bpm_mask)] = np.nan + array_corr[np.where(bpm_mask+excl_mask)] = np.nan if ndims == 2: array_corr_filt = frame_filter_lowpass(array_corr, mode=mode, fwhm_size=fwhm, @@ -1351,7 +1505,8 @@ def squash_v(array): elif psf.ndim == 2: psf = [psf]*nz elif psf.shape[0] != nz: - raise ValueError("input psf must have same z dimension as array") + raise ValueError( + "input psf must have same z dimension as array") for z in range(nz): array_corr_filt[z] = frame_filter_lowpass(array_corr[z], mode=mode, @@ -1375,7 +1530,7 @@ def squash_v(array): else: array_corr, recon_cube = res else: - if bpm_mask.ndim==2: + if bpm_mask.ndim == 2: bpm_mask = [bpm_mask]*nz if nproc is None: nproc = cpu_count()//2 @@ -1780,7 +1935,7 @@ def frame_fix_badpix_fft(array, bpm_mask, nit=500, tol=1, pad_fac=2, verbose=True, full_output=False): """ Function to interpolate bad pixels with the FFT-based algorithm in [AAC01]_. - + Parameters ---------- array : 2D ndarray @@ -1801,7 +1956,7 @@ def frame_fix_badpix_fft(array, bpm_mask, nit=500, tol=1, pad_fac=2, full_output: bool Whether to also return the reconstructed estimate f_hat of the input array. - + Returns ------- array_corr: 2D ndarray or list of 2D ndarray @@ -1823,7 +1978,6 @@ def frame_fix_badpix_fft(array, bpm_mask, nit=500, tol=1, pad_fac=2, nit_max = nit return_list = False - final_array_corr = [] final_f_est = [] @@ -1858,7 +2012,7 @@ def frame_fix_badpix_fft(array, bpm_mask, nit=500, tol=1, pad_fac=2, # 2. compute the new F_i - ## handle cases with no conjugate: + # handle cases with no conjugate: cond1 = (ind[0] == 0) and (ind[1] == 0) cond2 = (ind[0] == ny / 2) and (ind[1] == 0) cond3 = (ind[0] == 0) and (ind[1] == nx / 2) @@ -1898,7 +2052,7 @@ def frame_fix_badpix_fft(array, bpm_mask, nit=500, tol=1, pad_fac=2, # 5. Calculate new error - to check if still larger than tolerance Eg = np.sum(np.power(np.abs(G_i.flatten()), 2))/npix - if (return_list and it in nit) or (it==nit_max-1) or (Eg < tol): + if (return_list and it in nit) or (it == nit_max-1) or (Eg < tol): array_corr = g + np.fft.ifft2(F_est).real * (1 - w) # crop zeros to return to initial size From 176cd44532af312e188167b4894b75909255ee9b Mon Sep 17 00:00:00 2001 From: VChristiaens Date: Thu, 2 Feb 2023 00:02:35 +0100 Subject: [PATCH 05/14] Added option for exclusion mask to bad pixel identification/correction routines --- vip_hci/preproc/badpixremoval.py | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/vip_hci/preproc/badpixremoval.py b/vip_hci/preproc/badpixremoval.py index 2e5d464e..b6caff23 100644 --- a/vip_hci/preproc/badpixremoval.py +++ b/vip_hci/preproc/badpixremoval.py @@ -426,14 +426,18 @@ def _mp_clean_isolated(args): count_bp = np.sum(final_bpm) else: if excl_mask is None: - excl_mask = np.zeros(array.shape, dtype=bool) + excl_mask = np.zeros(array.shape[-2:], dtype=bool) + elif excl_mask.ndim == 3: + excl_mask = np.median(excl_mask, axis=0) elif excl_mask is None: - msg = "Input exclusion mask should have same shape as array\n" - assert excl_mask.shape == array.shape, msg + msg = "Input exclusion mask should have same last 2 dims as array\n" + assert excl_mask.shape == array.shape[-2:], msg ind_excl = np.where(excl_mask) if bpm_mask is None or not correct_only: if bpm_mask is None: - bpm_mask = np.zeros(array.shape, dtype=bool) + bpm_mask = np.zeros(array.shape[-2:], dtype=bool) + elif bpm_mask.ndim == 3: + bpm_mask = np.median(bpm_mask, axis=0) bpm_mask = bpm_mask+excl_mask ori_nan_mask = np.where(np.isnan(np.nanmean(array, axis=0))) ind = clip_array(np.nanmean(array, axis=0), sigma_clip, sigma_clip, @@ -845,7 +849,7 @@ def bp_removal_2d(array_corr, cy, cx, fwhm, sig, protect_mask, bpm_mask_ori, excl_mask, min_thr, half_res_y, mad, verbose): msg = "Input exclusion mask should have same shape as array\n" - assert excl_mask.shape == array.shape, msg + assert excl_mask.shape == array_corr.shape, msg ind_excl = np.where(excl_mask) n_x = array_corr.shape[1] From 39e6ec4d623daddb2fd1cf3daf2a792bb2f06d7b Mon Sep 17 00:00:00 2001 From: VChristiaens Date: Thu, 2 Feb 2023 13:09:34 +0100 Subject: [PATCH 06/14] Added option for exclusion mask to bad pixel identification/correction routines --- vip_hci/preproc/badpixremoval.py | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) diff --git a/vip_hci/preproc/badpixremoval.py b/vip_hci/preproc/badpixremoval.py index b6caff23..14f7cf48 100644 --- a/vip_hci/preproc/badpixremoval.py +++ b/vip_hci/preproc/badpixremoval.py @@ -489,8 +489,8 @@ def cube_fix_badpix_annuli(array, fwhm, cy=None, cx=None, sig=5., provided location of the star), in an input frame or cube. This function is faster than ``cube_fix_badpix_clump``; hence to be preferred in all cases where there is only one bright source with circularly - symmetric PSF. The bad pixel values are replaced by:\ - ann_median + random_poisson;\ + symmetric PSF. The bad pixel values are replaced by: + ann_median + random_poisson; where ann_median is the median of the annulus, and random_poisson is random noise picked from a Poisson distribution centered on ann_median. @@ -887,8 +887,7 @@ def bp_removal_2d(array_corr, cy, cx, fwhm, sig, protect_mask, bpm_mask_ori, circl_new = [] # 3/ Create a bad pixel map, by detecting them with clip_array - bpm_mask = bpm_mask_ori+excl_mask - bpm_mask[np.where(bpm_mask > 1)] = 1 + bpm_mask = bpm_mask_ori.astype(bool)+excl_mask.astype(bool) bp = clip_array(array_corr, sig, sig, bpm_mask, out_good=False, neighbor=True, num_neighbor=neighbor_box, mad=mad, half_res_y=half_res_y) @@ -911,7 +910,7 @@ def bp_removal_2d(array_corr, cy, cx, fwhm, sig, protect_mask, bpm_mask_ori, bpix_map[circl_new] = 0 bpix_map[ind_excl] = 0 nbpix_tbc = int(np.sum(bpix_map)) - bpix_map_cumul = np.zeros_like(bpix_map) + bpix_map_cumul = np.zeros(bpix_map.shape, dtype=bool) bpix_map_cumul[:] = bpix_map[:] nit = 0 @@ -932,11 +931,10 @@ def bp_removal_2d(array_corr, cy, cx, fwhm, sig, protect_mask, bpm_mask_ori, min_neighbors=nneig, half_res_y=half_res_y, verbose=verbose) bpm_mask += bpix_map_cumul - bpm_mask[np.where(bpm_mask > 1)] = 1 bp = clip_array(array_corr, sig, sig, bpm_mask, out_good=False, neighbor=True, num_neighbor=neighbor_box, mad=mad, half_res_y=half_res_y) - bpix_map = np.zeros_like(array_corr) + bpix_map = np.zeros(array_corr, dtype=bool) bpix_map[bp] = 1 if min_thr is not None: cond1 = array_corr > min_thr[0] From edc60103461920023f71890e7ff7db64505ba7f4 Mon Sep 17 00:00:00 2001 From: VChristiaens Date: Thu, 2 Feb 2023 13:21:30 +0100 Subject: [PATCH 07/14] Solved photutils deprecation warnings --- vip_hci/fm/fakecomp.py | 69 ++++++++++++++++++----------------- vip_hci/metrics/contrcurve.py | 37 ++++++++++--------- vip_hci/metrics/snr_source.py | 28 +++++++------- vip_hci/preproc/badframes.py | 7 +++- 4 files changed, 76 insertions(+), 65 deletions(-) diff --git a/vip_hci/fm/fakecomp.py b/vip_hci/fm/fakecomp.py index c141bb5e..43e186fc 100644 --- a/vip_hci/fm/fakecomp.py +++ b/vip_hci/fm/fakecomp.py @@ -17,6 +17,10 @@ from scipy.interpolate import interp1d from packaging import version import photutils +try: + from photutils.aperture import aperture_photometry, CircularAperture +except: + from photutils import aperture_photometry, CircularAperture if version.parse(photutils.__version__) >= version.parse('0.3'): # for photutils version >= '0.3' use photutils.centroids.centroid_com from photutils.centroids import centroid_com as cen_com @@ -32,8 +36,8 @@ def cube_inject_companions(array, psf_template, angle_list, flevel, rad_dists, plsc=None, n_branches=1, theta=0, imlib='vip-fft', - interpolation='lanczos4', transmission=None, - radial_gradient=False, full_output=False, + interpolation='lanczos4', transmission=None, + radial_gradient=False, full_output=False, verbose=False, nproc=1): """ Injects fake companions in branches, at given radial distances. @@ -112,14 +116,14 @@ def cube_inject_companions(array, psf_template, angle_list, flevel, rad_dists, """ if nproc is None: nproc = cpu_count()//2 - + def _cube_inject_adi(array, psf_template, angle_list, flevel, plsc, rad_dists, n_branches=1, theta=0, imlib='vip-fft', interpolation='lanczos4', transmission=None, radial_gradient=False, verbose=False): if np.isscalar(flevel): flevel = np.ones_like(angle_list)*flevel - + if transmission is not None: # last radial separation should be beyond the edge of frame interp_trans = interp1d(transmission[0], transmission[1]) @@ -159,9 +163,9 @@ def _cube_inject_adi(array, psf_template, angle_list, flevel, plsc, # check the effect of transmission on a single PSF tmp psf_trans = frame_rotate(fc_fr_rad[0], - -(ang*180/np.pi-angle_list[0]), - imlib=imlib_rot, - interpolation=interpolation) + -(ang*180/np.pi-angle_list[0]), + imlib=imlib_rot, + interpolation=interpolation) # shift_y = rad * np.sin(ang - np.deg2rad(angle_list[0])) # shift_x = rad * np.cos(ang - np.deg2rad(angle_list[0])) @@ -174,19 +178,19 @@ def _cube_inject_adi(array, psf_template, angle_list, flevel, plsc, fc_fr_rad = interp_trans(rad)*fc_fr if nproc == 1: for fr in range(nframes): - array_out[fr] += _frame_shift_fcp(fc_fr_rad[fr], - array[fr], rad, ang, - angle_list[fr], - flevel[fr], size_fc, - imlib_sh, imlib_rot, + array_out[fr] += _frame_shift_fcp(fc_fr_rad[fr], + array[fr], rad, ang, + angle_list[fr], + flevel[fr], size_fc, + imlib_sh, imlib_rot, interpolation, - transmission, + transmission, radial_gradient) else: res = pool_map(nproc, _frame_shift_fcp, iterable(fc_fr_rad), - iterable(array), rad, ang, - iterable(angle_list), iterable(flevel), - size_fc, imlib_sh, imlib_rot, interpolation, + iterable(array), rad, ang, + iterable(angle_list), iterable(flevel), + size_fc, imlib_sh, imlib_rot, interpolation, transmission, radial_gradient) array_out += np.array(res) @@ -318,34 +322,33 @@ def _cube_inject_adi(array, psf_template, angle_list, flevel, plsc, return array_out -def _frame_shift_fcp(fc_fr_rad, array, rad, ang, derot_ang, flevel, size_fc, - imlib_sh, imlib_rot, interpolation, transmission, +def _frame_shift_fcp(fc_fr_rad, array, rad, ang, derot_ang, flevel, size_fc, + imlib_sh, imlib_rot, interpolation, transmission, radial_gradient): """ Specific cube shift algorithm for injection of fake companions """ - + ceny, cenx = frame_center(array) sizey = array.shape[-2] sizex = array.shape[-1] - + array_sh = np.zeros_like(array) - + w = int(np.ceil(size_fc/2)) if size_fc % 2: # new convention w -= 1 sty = int(ceny) - w stx = int(cenx) - w - + shift_y = rad * np.sin(ang - np.deg2rad(derot_ang)) shift_x = rad * np.cos(ang - np.deg2rad(derot_ang)) if transmission is not None and radial_gradient: - fc_fr_ang = frame_rotate(fc_fr_rad, -(ang*180/np.pi -derot_ang), + fc_fr_ang = frame_rotate(fc_fr_rad, -(ang*180/np.pi - derot_ang), imlib=imlib_rot, interpolation=interpolation) else: fc_fr_ang = fc_fr_rad.copy() - # sub-px shift (within PSF template frame) dsy = shift_y-int(shift_y) dsx = shift_x-int(shift_x) @@ -372,9 +375,9 @@ def _frame_shift_fcp(fc_fr_rad, array, rad, ang, derot_ang, flevel, size_fc, if xN > sizex: p_xN -= xN-sizex xN = sizex - + array_sh[y0:yN, x0:xN] = flevel*fc_fr_ang[p_y0:p_yN, p_x0:p_xN] - + return array_sh @@ -580,7 +583,7 @@ def collapse_psf_cube(array, size, fwhm=4, verbose=True, collapse='mean'): def normalize_psf(array, fwhm='fit', size=None, threshold=None, mask_core=None, model='gauss', imlib='vip-fft', interpolation='lanczos4', - force_odd=True, correct_outliers=True, full_output=False, + force_odd=True, correct_outliers=True, full_output=False, verbose=True, debug=False): """ Normalizes a PSF (2d or 3d array), to have the flux in a 1xFWHM aperture equal to one. It also allows to crop the array and center the PSF @@ -656,7 +659,7 @@ def psf_norm_2d(psf, fwhm, threshold, mask_core, full_output, verbose): shiftx, shifty = centrx - cx, centry - cy psf = frame_shift(psf, -shifty, -shiftx, imlib=imlib, interpolation=interpolation) - + for _ in range(2): centry, centrx = fit_2d(psf, full_output=False, debug=False) if np.isnan(centry) or np.isnan(centrx): @@ -667,9 +670,9 @@ def psf_norm_2d(psf, fwhm, threshold, mask_core, full_output, verbose): interpolation=interpolation) # we check whether the flux is normalized and fix it if needed - fwhm_aper = photutils.CircularAperture((cx, cy), fwhm/2) - fwhm_aper_phot = photutils.aperture_photometry(psf, fwhm_aper, - method='exact') + fwhm_aper = CircularAperture((cx, cy), fwhm/2) + fwhm_aper_phot = aperture_photometry(psf, fwhm_aper, + method='exact') fwhm_flux = np.array(fwhm_aper_phot['aperture_sum']) if fwhm_flux > 1.1 or fwhm_flux < 0.9: @@ -779,9 +782,9 @@ def psf_norm_2d(psf, fwhm, threshold, mask_core, full_output, verbose): print_precision(fwhm) # Replace outliers if needed if correct_outliers: - if np.sum(np.isnan(fwhm))>0: + if np.sum(np.isnan(fwhm)) > 0: for f in range(n): - if np.isnan(fwhm[f]) and f!=0 and f!=n-1: + if np.isnan(fwhm[f]) and f != 0 and f != n-1: fwhm[f] = np.nanmean(np.array([fwhm[f-1], fwhm[f+1]])) elif np.isnan(fwhm[f]): diff --git a/vip_hci/metrics/contrcurve.py b/vip_hci/metrics/contrcurve.py index 5de27773..8b92440e 100644 --- a/vip_hci/metrics/contrcurve.py +++ b/vip_hci/metrics/contrcurve.py @@ -13,7 +13,10 @@ import numpy as np import pandas as pd -import photutils +try: + from photutils.aperture import aperture_photometry, CircularAperture +except: + from photutils import aperture_photometry, CircularAperture from inspect import getfullargspec from scipy.interpolate import InterpolatedUnivariateSpline from scipy import stats @@ -23,16 +26,16 @@ from ..fm import (cube_inject_companions, frame_inject_companion, normalize_psf) from ..config import time_ini, timing -from ..config.utils_conf import sep, vip_figsize, vip_figdpi +from ..config.utils_conf import vip_figsize, vip_figdpi from ..var import frame_center, dist def contrast_curve(cube, angle_list, psf_template, fwhm, pxscale, starphot, algo, sigma=5, nbranch=1, theta=0, inner_rad=1, fc_rad_sep=3, - noise_sep=1, wedge=(0, 360), fc_snr=100, student=True, - transmission=None, smooth=True, interp_order=2, plot=True, - dpi=vip_figdpi, debug=False, verbose=True, full_output=False, - save_plot=None, object_name=None, frame_size=None, + noise_sep=1, wedge=(0, 360), fc_snr=100, student=True, + transmission=None, smooth=True, interp_order=2, plot=True, + dpi=vip_figdpi, debug=False, verbose=True, full_output=False, + save_plot=None, object_name=None, frame_size=None, fix_y_lim=(), figsize=vip_figsize, **algo_dict): """ Computes the contrast curve at a given confidence (``sigma``) level for an ADI cube or ADI+IFS cube. The contrast is calculated as @@ -174,7 +177,7 @@ def contrast_curve(cube, angle_list, psf_template, fwhm, pxscale, starphot, if transmission is not None: if len(transmission) != 2 and len(transmission) != cube.shape[0]+1: msg = 'Wrong shape for transmission should be 2xn_rad or (nch+1) ' - msg +='x n_rad, instead of {}'.format(transmission.shape) + msg += 'x n_rad, instead of {}'.format(transmission.shape) raise TypeError(msg) if isinstance(fwhm, (np.ndarray, list)): @@ -204,7 +207,7 @@ def contrast_curve(cube, angle_list, psf_template, fwhm, pxscale, starphot, verbose_thru = True res_throug = throughput(cube, angle_list, psf_template, fwhm, algo=algo, nbranch=nbranch, theta=theta, inner_rad=inner_rad, - fc_rad_sep=fc_rad_sep, wedge=wedge, fc_snr=fc_snr, + fc_rad_sep=fc_rad_sep, wedge=wedge, fc_snr=fc_snr, full_output=True, verbose=verbose_thru, **algo_dict) vector_radd = res_throug[3] if res_throug[0].shape[0] > 1: @@ -247,7 +250,7 @@ def contrast_curve(cube, angle_list, psf_template, fwhm, pxscale, starphot, ntransmission[0] = trans_rad_list ntransmission[j+1] = trans_list transmission = ntransmission.copy() - if t_nz>2: #take the mean transmission over all wavelengths + if t_nz > 2: # take the mean transmission over all wavelengths ntransmission = np.zeros([2, len(trans_rad_list)]) ntransmission[0] = transmission[0] ntransmission[1] = np.mean(transmission[1:], axis=0) @@ -657,11 +660,11 @@ def throughput(cube, angle_list, psf_template, fwhm, algo, nbranch=1, theta=0, wedge=wedge) if scaling is not None: noise_noscal, _, _ = noise_per_annulus(frame_nofc_noscal, - separation=fwhm_med, + separation=fwhm_med, fwhm=fwhm_med, wedge=wedge) else: noise_noscal = noise.copy() - + vector_radd = vector_radd[inner_rad-1:] noise = noise[inner_rad-1:] res_level = res_level[inner_rad-1:] @@ -707,7 +710,7 @@ def throughput(cube, angle_list, psf_template, fwhm, algo, nbranch=1, theta=0, rad_dists=[radvec[i]], theta=br*angle_branch + theta, - nproc=nproc, imlib=imlib, + nproc=nproc, imlib=imlib, interpolation=interpolation, verbose=False) y = cy + radvec[i] * np.sin(np.deg2rad(br * angle_branch + @@ -936,8 +939,8 @@ def find_coords(rad, sep, init_angle, fin_angle): yy += centery xx += centerx - apertures = photutils.CircularAperture(np.array((xx, yy)).T, fwhm/2) - fluxes = photutils.aperture_photometry(array, apertures) + apertures = CircularAperture(np.array((xx, yy)).T, fwhm/2) + fluxes = aperture_photometry(array, apertures) fluxes = np.array(fluxes['aperture_sum']) noise_ann = np.std(fluxes) @@ -1003,9 +1006,9 @@ def aperture_flux(array, yc, xc, fwhm, ap_factor=1, mean=False, verbose=False): values = array[ind] obj_flux = np.mean(values) else: - aper = photutils.CircularAperture((x, y), (ap_factor*fwhm)/2) - obj_flux = photutils.aperture_photometry(array, aper, - method='exact') + aper = CircularAperture((x, y), (ap_factor*fwhm)/2) + obj_flux = aperture_photometry(array, aper, + method='exact') obj_flux = np.array(obj_flux['aperture_sum']) flux[i] = obj_flux diff --git a/vip_hci/metrics/snr_source.py b/vip_hci/metrics/snr_source.py index 1206e5d1..f4c37030 100644 --- a/vip_hci/metrics/snr_source.py +++ b/vip_hci/metrics/snr_source.py @@ -14,7 +14,10 @@ 'frame_report'] import numpy as np -import photutils +try: + from photutils.aperture import aperture_photometry, CircularAperture +except: + from photutils import aperture_photometry, CircularAperture from scipy.stats import norm, t from hciplot import plot_frames from skimage.draw import disk, circle_perimeter @@ -217,6 +220,7 @@ def _snr_approx(array, source_xy, fwhm, centery, centerx): snr_value = signal / noise return sourcey, sourcex, snr_value + def indep_ap_centers(array, source_xy, fwhm, exclude_negative_lobes=False, exclude_theta_range=None, no_gap=False): """ @@ -311,8 +315,9 @@ def indep_ap_centers(array, source_xy, fwhm, exclude_negative_lobes=False, return yy, xx + def snr(array, source_xy, fwhm, full_output=False, array2=None, use2alone=False, - exclude_negative_lobes=False, exclude_theta_range=None, plot=False, + exclude_negative_lobes=False, exclude_theta_range=None, plot=False, verbose=False): """ Calculate the S/N (signal to noise ratio) of a test resolution element @@ -395,14 +400,12 @@ def snr(array, source_xy, fwhm, full_output=False, array2=None, use2alone=False, rad = fwhm/2. - apertures = photutils.CircularAperture( - zip(xx, yy), r=rad) # Coordinates (X,Y) - fluxes = photutils.aperture_photometry(array, apertures, method='exact') + apertures = CircularAperture(zip(xx, yy), r=rad) # Coordinates (X,Y) + fluxes = aperture_photometry(array, apertures, method='exact') fluxes = np.array(fluxes['aperture_sum']) if array2 is not None: - fluxes2 = photutils.aperture_photometry(array2, apertures, - method='exact') + fluxes2 = aperture_photometry(array2, apertures, method='exact') fluxes2 = np.array(fluxes2['aperture_sum']) if use2alone: fluxes = np.concatenate(([fluxes[0]], fluxes2[:])) @@ -538,9 +541,8 @@ def frame_report(array, fwhm, source_xy=None, verbose=True): print('Coords of chosen px (X,Y) = {:.1f}, {:.1f}'.format(x, y)) # we get integrated flux on aperture with diameter=1FWHM - aper = photutils.CircularAperture((x, y), r=fwhm / 2.) - obj_flux_i = photutils.aperture_photometry(array, aper, - method='exact') + aper = CircularAperture((x, y), r=fwhm / 2.) + obj_flux_i = aperture_photometry(array, aper, method='exact') obj_flux_i = obj_flux_i['aperture_sum'][0] # we get the mean and stddev of SNRs on aperture @@ -578,8 +580,8 @@ def frame_report(array, fwhm, source_xy=None, verbose=True): print('Coords of Max px (X,Y) = {:.1f}, {:.1f}'.format(x, y)) # we get integrated flux on aperture with diameter=1FWHM - aper = photutils.CircularAperture((x, y), r=fwhm / 2.) - obj_flux_i = photutils.aperture_photometry(array, aper, method='exact') + aper = CircularAperture((x, y), r=fwhm / 2.) + obj_flux_i = aperture_photometry(array, aper, method='exact') obj_flux_i = obj_flux_i['aperture_sum'][0] # we get the mean and stddev of SNRs on aperture @@ -607,4 +609,4 @@ def frame_report(array, fwhm, source_xy=None, verbose=True): print(msg3.format(stdsnr_i)) print(sep) - return source_xy, obj_flux, snr_centpx, meansnr_pixels \ No newline at end of file + return source_xy, obj_flux, snr_centpx, meansnr_pixels diff --git a/vip_hci/preproc/badframes.py b/vip_hci/preproc/badframes.py index 633e660d..78971dc5 100644 --- a/vip_hci/preproc/badframes.py +++ b/vip_hci/preproc/badframes.py @@ -12,7 +12,10 @@ import numpy as np import pandas as pn from matplotlib import pyplot as plt -from photutils import DAOStarFinder +try: + from photutils.detection import DAOStarFinder +except: + from photutils import DAOStarFinder from astropy.stats import sigma_clip from ..var import get_annulus_segments from ..config import time_ini, timing, check_array @@ -372,4 +375,4 @@ def cube_detect_badfr_correlation(array, frame_ref, crop_size=30, if full_output: return good_index_list, bad_index_list, distances else: - return good_index_list, bad_index_list \ No newline at end of file + return good_index_list, bad_index_list From dabf0eb52cdb1fd3d7e42a1a14b22d728de95a9a Mon Sep 17 00:00:00 2001 From: VChristiaens Date: Thu, 2 Feb 2023 19:42:49 +0100 Subject: [PATCH 08/14] Added option for exclusion mask to bad pixel identification/correction routines --- vip_hci/preproc/badpixremoval.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/vip_hci/preproc/badpixremoval.py b/vip_hci/preproc/badpixremoval.py index 14f7cf48..a91f5604 100644 --- a/vip_hci/preproc/badpixremoval.py +++ b/vip_hci/preproc/badpixremoval.py @@ -934,7 +934,7 @@ def bp_removal_2d(array_corr, cy, cx, fwhm, sig, protect_mask, bpm_mask_ori, bp = clip_array(array_corr, sig, sig, bpm_mask, out_good=False, neighbor=True, num_neighbor=neighbor_box, mad=mad, half_res_y=half_res_y) - bpix_map = np.zeros(array_corr, dtype=bool) + bpix_map = np.zeros(array_corr.shape, dtype=bool) bpix_map[bp] = 1 if min_thr is not None: cond1 = array_corr > min_thr[0] From 6967907bb91743d17d3922862116355968b28241 Mon Sep 17 00:00:00 2001 From: VChristiaens Date: Thu, 2 Feb 2023 23:38:35 +0100 Subject: [PATCH 09/14] Added option for exclusion mask to bad pixel identification/correction routines --- vip_hci/preproc/badpixremoval.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/vip_hci/preproc/badpixremoval.py b/vip_hci/preproc/badpixremoval.py index a91f5604..1148f0a4 100644 --- a/vip_hci/preproc/badpixremoval.py +++ b/vip_hci/preproc/badpixremoval.py @@ -930,7 +930,7 @@ def bp_removal_2d(array_corr, cy, cx, fwhm, sig, protect_mask, bpm_mask_ori, neighbor_box=neighbor_box, min_neighbors=nneig, half_res_y=half_res_y, verbose=verbose) - bpm_mask += bpix_map_cumul + bpm_mask = None # known bad ones are corrected above +=bpix_map_cumul bp = clip_array(array_corr, sig, sig, bpm_mask, out_good=False, neighbor=True, num_neighbor=neighbor_box, mad=mad, half_res_y=half_res_y) @@ -965,6 +965,8 @@ def bp_removal_2d(array_corr, cy, cx, fwhm, sig, protect_mask, bpm_mask_ori, if bpm_mask is None or not correct_only: if (cy is None or cx is None) and protect_mask: cy, cx = frame_center(array) + if excl_mask is None: + excl_mask = np.zeros(array_corr.shape, dtype=bool) array_corr, bpix_map_cumul = bp_removal_2d(array_corr, cy, cx, fwhm, sig, protect_mask, bpm_mask, excl_mask, From 8e778c8bb706a85ae7fd2d2cfe5a619639c49b6b Mon Sep 17 00:00:00 2001 From: VChristiaens Date: Fri, 3 Feb 2023 12:51:08 +0100 Subject: [PATCH 10/14] Added option for exclusion mask to bad pixel identification/correction routines --- vip_hci/preproc/badpixremoval.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/vip_hci/preproc/badpixremoval.py b/vip_hci/preproc/badpixremoval.py index 1148f0a4..ecb15602 100644 --- a/vip_hci/preproc/badpixremoval.py +++ b/vip_hci/preproc/badpixremoval.py @@ -887,7 +887,9 @@ def bp_removal_2d(array_corr, cy, cx, fwhm, sig, protect_mask, bpm_mask_ori, circl_new = [] # 3/ Create a bad pixel map, by detecting them with clip_array - bpm_mask = bpm_mask_ori.astype(bool)+excl_mask.astype(bool) + bpm_mask = excl_mask.copy() + if bpm_mask_ori is not None: + bpm_mask += bpm_mask_ori.astype(bool) bp = clip_array(array_corr, sig, sig, bpm_mask, out_good=False, neighbor=True, num_neighbor=neighbor_box, mad=mad, half_res_y=half_res_y) From b4720c2b53d022b9ff155a1689cf76c60fd99b05 Mon Sep 17 00:00:00 2001 From: VChristiaens Date: Fri, 3 Feb 2023 15:46:38 +0100 Subject: [PATCH 11/14] Optimized NEGFC for nmf-annular (in addition to existing optimization for pca-annular) --- vip_hci/fm/negfc_fmerit.py | 50 +++++++++++++++++++++++--------------- 1 file changed, 30 insertions(+), 20 deletions(-) diff --git a/vip_hci/fm/negfc_fmerit.py b/vip_hci/fm/negfc_fmerit.py index 1adf3ffb..b0e3a50c 100644 --- a/vip_hci/fm/negfc_fmerit.py +++ b/vip_hci/fm/negfc_fmerit.py @@ -13,7 +13,7 @@ from ..fm import cube_inject_companions from ..var import (frame_center, get_annular_wedge, cube_filter_highpass, get_annulus_segments) -from ..psfsub import pca_annulus, pca_annular, pca +from ..psfsub import pca_annulus, pca_annular, nmf_annular, pca from ..preproc import cube_crop_frames @@ -313,11 +313,11 @@ def get_values_optimize(cube, angs, ncomp, annulus_width, aperture_radius, res = pca_annulus(cube, angs, ncomp, annulus_width, r_guess, cube_ref, svd_mode, scaling, imlib=imlib, interpolation=interpolation, collapse=collapse, - collapse_ifs=collapse_ifs, weights=weights, + collapse_ifs=collapse_ifs, weights=weights, nproc=nproc) - elif algo == pca_annular: - + elif algo == pca_annular or algo == nmf_annular: + tol = algo_options.get('tol', 1e-1) min_frames_lib = algo_options.get('min_frames_lib', 2) max_frames_lib = algo_options.get('max_frames_lib', 200) @@ -333,15 +333,25 @@ def get_values_optimize(cube, angs, ncomp, annulus_width, aperture_radius, else: crop_cube = cube pad = 0 - res_tmp = pca_annular(crop_cube, angs, radius_int=radius_int, fwhm=fwhm, - asize=annulus_width, delta_rot=delta_rot, - ncomp=ncomp, svd_mode=svd_mode, scaling=scaling, - imlib=imlib, interpolation=interpolation, - collapse=collapse, collapse_ifs=collapse_ifs, - weights=weights, tol=tol, nproc=nproc, - min_frames_lib=min_frames_lib, - max_frames_lib=max_frames_lib, full_output=False, - verbose=False) + if algo == pca_annular: + res_tmp = algo(crop_cube, angs, radius_int=radius_int, fwhm=fwhm, + asize=annulus_width, delta_rot=delta_rot, + ncomp=ncomp, svd_mode=svd_mode, scaling=scaling, + imlib=imlib, interpolation=interpolation, + collapse=collapse, collapse_ifs=collapse_ifs, + weights=weights, tol=tol, nproc=nproc, + min_frames_lib=min_frames_lib, + max_frames_lib=max_frames_lib, full_output=False, + verbose=False) + else: + res_tmp = algo(crop_cube, angs, radius_int=radius_int, fwhm=fwhm, + asize=annulus_width, delta_rot=delta_rot, + ncomp=ncomp, scaling=scaling, imlib=imlib, + interpolation=interpolation, collapse=collapse, + weights=weights, nproc=nproc, + min_frames_lib=min_frames_lib, + max_frames_lib=max_frames_lib, full_output=False, + verbose=False) # pad again now res = np.pad(res_tmp, pad, mode='constant', constant_values=0) @@ -358,25 +368,25 @@ def get_values_optimize(cube, angs, ncomp, annulus_width, aperture_radius, indices = disk((posy, posx), radius=aperture_radius*fwhm) yy, xx = indices - + # also consider indices of the annulus for pca_annulus if algo == pca_annulus: fr_size = res.shape[-1] inner_rad = r_guess-annulus_width/2 yy_a, xx_a = get_annulus_segments((fr_size, fr_size), inner_rad, annulus_width, nsegm=1)[0] - ## only consider overlapping indices + # only consider overlapping indices yy_f = [] xx_f = [] for i in range(len(yy)): - ind_y = np.where(yy_a==yy[i]) + ind_y = np.where(yy_a == yy[i]) for j in ind_y[0]: - if xx[i]==xx_a[j]: + if xx[i] == xx_a[j]: yy_f.append(yy[i]) xx_f.append(xx[i]) yy = np.array(yy_f, dtype=int) xx = np.array(xx_f, dtype=int) - + if collapse is None: values = res[:, yy, xx].ravel() else: @@ -527,7 +537,7 @@ def get_mu_and_sigma(cube, angs, ncomp, annulus_width, aperture_radius, fwhm, pad = int((cube.shape[1]-crop_sz)/2) crop_cube = cube_crop_frames(cube, crop_sz, verbose=False) else: - pad=0 + pad = 0 crop_cube = cube pca_res_tmp = pca_annular(crop_cube, angs, radius_int=radius_int, @@ -606,4 +616,4 @@ def get_mu_and_sigma(cube, angs, ncomp, annulus_width, aperture_radius, fwhm, ddof = min(int(npx*(1.-(1./area)))+1, npx-1) sigma = np.std(all_res, ddof=ddof) - return mu, sigma \ No newline at end of file + return mu, sigma From fb2d57e43699de5aa5a96fa9a48b4d114e055816 Mon Sep 17 00:00:00 2001 From: VChristiaens Date: Fri, 3 Feb 2023 17:32:19 +0100 Subject: [PATCH 12/14] Added option for exclusion mask to bad pixel identification/correction routines --- vip_hci/preproc/badpixremoval.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/vip_hci/preproc/badpixremoval.py b/vip_hci/preproc/badpixremoval.py index ecb15602..e8a0a01b 100644 --- a/vip_hci/preproc/badpixremoval.py +++ b/vip_hci/preproc/badpixremoval.py @@ -863,8 +863,16 @@ def bp_removal_2d(array_corr, cy, cx, fwhm, sig, protect_mask, bpm_mask_ori, n_y = int(n_y/2) frame = array_corr.copy() array_corr = np.zeros([n_y, n_x]) + excl_mask_corr = np.zeros([n_y, n_x]) for yy in range(n_y): array_corr[yy] = frame[2*yy] + excl_mask_corr[yy] = excl_mask[2*yy] + excl_mask = excl_mask_corr + if bpm_mask_ori is not None: + bpm_mask_tmp = np.zeros([n_y, n_x]) + for yy in range(n_y): + bpm_mask_tmp[yy] = bpm_mask_ori[2*yy] + bpm_mask_ori = bpm_mask_tmp fwhm_round = int(round(fwhm)) # This should reduce the chance to accidently correct a bright planet: From 09362ca63d71142ef9eee289057c9fa52dca9874 Mon Sep 17 00:00:00 2001 From: VChristiaens Date: Mon, 6 Feb 2023 19:16:12 +0100 Subject: [PATCH 13/14] Bug fix in function to correct bad pixels through interpolation in the case of half y resolution --- vip_hci/preproc/badpixremoval.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/vip_hci/preproc/badpixremoval.py b/vip_hci/preproc/badpixremoval.py index e8a0a01b..d5740d0d 100644 --- a/vip_hci/preproc/badpixremoval.py +++ b/vip_hci/preproc/badpixremoval.py @@ -1476,10 +1476,10 @@ def squash_v(array): if ny % 2: raise ValueError("Input array y dimension should be even") nny = ny//2 - new_array_corr = np.zeros([nny, nx]) + new_array = np.zeros([nny, nx]) for y in range(nny): - new_array_corr[y] = array_corr[int(y*2)] - return new_array_corr + new_array[y] = array[int(y*2)] + return new_array if ndims == 2: array_corr = squash_v(array) From 8d3ba2f7f5f3b97a1fa757b9eb63369398427a84 Mon Sep 17 00:00:00 2001 From: VChristiaens Date: Thu, 9 Feb 2023 00:50:45 +0100 Subject: [PATCH 14/14] Added new bpm_mask, excl_mask and bad_values options to cube_fix_badpix_annuli --- vip_hci/preproc/badpixremoval.py | 117 ++++++++++++++---- vip_hci/preproc/recentering.py | 204 +++++++++++++++---------------- 2 files changed, 196 insertions(+), 125 deletions(-) diff --git a/vip_hci/preproc/badpixremoval.py b/vip_hci/preproc/badpixremoval.py index d5740d0d..7e52d046 100644 --- a/vip_hci/preproc/badpixremoval.py +++ b/vip_hci/preproc/badpixremoval.py @@ -480,10 +480,11 @@ def _mp_clean_isolated(args): return array_out -def cube_fix_badpix_annuli(array, fwhm, cy=None, cx=None, sig=5., - protect_mask=0, r_in_std=50, r_out_std=None, - verbose=True, half_res_y=False, min_thr=None, - max_thr=None, full_output=False): +def cube_fix_badpix_annuli(array, fwhm, cy=None, cx=None, sig=5., bpm_mask=None, + protect_mask=0, excl_mask=None, r_in_std=50, + r_out_std=None, verbose=True, half_res_y=False, + min_thr=None, max_thr=None, bad_values=None, + full_output=False): """ Function to correct the bad pixels annulus per annulus (centered on the provided location of the star), in an input frame or cube. @@ -511,10 +512,20 @@ def cube_fix_badpix_annuli(array, fwhm, cy=None, cx=None, sig=5., sig: Float scalar, optional Number of stddev above or below the median of the pixels in the same annulus, to consider a pixel as bad. + bpm_mask: 3D or 2D array, opt + Input bad pixel array. If 2D and array is 3D: should have same last 2 + dimensions as array. If 3D, should have exact same dimensions as input + array. If not provided, the algorithm will attempt to identify bad pixel + clumps automatically. protect_mask : int or float, optional If larger than 0, radius of a circular aperture (at the center of the frames) in which no bad pixels will be identified. This can be useful to protect the star and vicinity. + excl_mask : numpy ndarray, optional + Binary mask with 1 in areas that should not be considered as good + neighbouring pixels during the identification of bad pixels. These + should not be considered as bad pixels to be corrected neither (i.e. + different to bpm_mask). r_in_std: float or None, optional Inner radius (in pixels) of the annulus used for the calculation of the standard deviation of the background noise - used as a min threshold @@ -535,6 +546,9 @@ def cube_fix_badpix_annuli(array, fwhm, cy=None, cx=None, sig=5., Any pixel whose value is lower (resp. larger) than this threshold will be automatically considered bad and hence sigma_filtered. If None, it is not used. + bad_values: list or None, optional + If not None, should correspond to a list of known bad values (e.g. 0). + These pixels will be added to the input bad pixel map. full_output: bool, {False,True}, optional Whether to return as well the cube of bad pixel maps and the cube of defined annuli. @@ -558,9 +572,25 @@ def cube_fix_badpix_annuli(array, fwhm, cy=None, cx=None, sig=5., if max_thr is None: max_thr = np.amax(array)-1 - def bp_removal_2d(array, cy, cx, fwhm, sig, protect_mask, r_in_std, - r_out_std, verbose): + if bpm_mask is not None: + msg = "Input bad pixel mask should have same last 2 dims as array\n" + assert bpm_mask.shape[-2:] == array.shape[-2:], msg + bpm_mask = bpm_mask.astype('bool') + + if bad_values is not None: + if bpm_mask is None: + bpm_mask = np.zeros(array.shape, dtype=bool) + for bad in bad_values: + bpm_mask[np.where(array == bad)] = 1 + + def bp_removal_2d(array, cy, cx, fwhm, sig, protect_mask, bpm_mask_ori, + excl_mask, r_in_std, r_out_std, verbose): + msg = "Input exclusion mask should have same shape as array\n" + assert excl_mask.shape == array.shape, msg + ind_excl = np.where(excl_mask) + + frame = array.copy() n_x = array.shape[1] n_y = array.shape[0] @@ -572,10 +602,17 @@ def bp_removal_2d(array, cy, cx, fwhm, sig, protect_mask, r_in_std, raise ValueError(msg+msg2) n_y = int(n_y/2) cy = int(cy/2) - frame = array.copy() array = np.zeros([n_y, n_x]) + excl_mask_corr = np.zeros([n_y, n_x]) for yy in range(n_y): array[yy] = frame[2*yy] + excl_mask_corr[yy] = excl_mask[2*yy] + excl_mask = excl_mask_corr + if bpm_mask_ori is not None: + bpm_mask_tmp = np.zeros([n_y, n_x]) + for yy in range(n_y): + bpm_mask_tmp[yy] = bpm_mask_ori[2*yy] + bpm_mask_ori = bpm_mask_tmp # 1/ Stddev of background if r_in_std or r_out_std: @@ -614,6 +651,12 @@ def bp_removal_2d(array, cy, cx, fwhm, sig, protect_mask, r_in_std, std_neig = np.zeros(nrad) neighbours = np.zeros([nrad, n_y*n_x]) + bpm_mask = excl_mask.copy() + if bpm_mask_ori is not None: + bpm_mask += bpm_mask_ori.astype(bool) + + ind_bad = np.where(bpm_mask) + for rr in range(nrad): if rr > int(d_bord_max/ann_width): # just to merge farthest annuli with very few elements @@ -641,6 +684,7 @@ def bp_removal_2d(array, cy, cx, fwhm, sig, protect_mask, r_in_std, big_ell_frame[big_ell_idx] = 1 if rr != 0: sma_ell_frame[small_ell_idx] = 1 + sma_ell_frame[ind_bad] = 1 ann_frame = big_ell_frame - sma_ell_frame n_neig[rr] = ann_frame[np.where(ann_frame)].shape[0] neighbours[rr, :n_neig[rr]] = array[np.where(ann_frame)] @@ -711,32 +755,50 @@ def bp_removal_2d(array, cy, cx, fwhm, sig, protect_mask, r_in_std, bpix_map[yy] = frame_bpix[int(yy/2)] ann_frame_cumul[yy] = ann_frame[int(yy/2)] + # Excluded pixels should take back their original values + array_corr[ind_excl] = frame[ind_excl] + bpix_map[ind_excl] = 0 + return array_corr, bpix_map, ann_frame_cumul if cy is None or cx is None: cy, cx = frame_center(array) if ndims == 2: + if excl_mask is None: + excl_mask = np.zeros(array.shape, dtype=bool) array_corr, bpix_map, ann_frame_cumul = bp_removal_2d(array, cy, cx, fwhm, sig, protect_mask, - r_in_std, r_out_std, + bpm_mask, + excl_mask, + r_in_std, + r_out_std, verbose) if ndims == 3: array_corr = array.copy() n_z = array.shape[0] bpix_map = np.zeros_like(array) ann_frame_cumul = np.zeros_like(array) - if isinstance(fwhm, (int, float)): + if np.isscalar(fwhm): fwhm = [fwhm]*n_z - if isinstance(cy, (float, int)) and isinstance(cx, (float, int)): + if np.isscalar(cx) and np.isscalar(cy): cy = [cy]*n_z cx = [cx]*n_z + if bpm_mask is None: + bpm_mask = np.zeros(array_corr.shape, dtype=bool) + elif bpm_mask.ndim == 2: + bpm_mask = np.array([bpm_mask]*n_z, dtype=bool) + if excl_mask is None: + excl_mask = np.zeros(array_corr.shape, dtype=bool) + elif excl_mask.ndim == 2: + excl_mask = np.array([excl_mask]*n_z, dtype=bool) for i in range(n_z): if verbose: print('************Frame # ', i, ' *************') print('centroid assumed at coords:', cx[i], cy[i]) res_i = bp_removal_2d(array[i], cy[i], cx[i], fwhm[i], sig, - protect_mask, r_in_std, r_out_std, verbose) + protect_mask, bpm_mask[i], excl_mask[i], + r_in_std, r_out_std, verbose) array_corr[i], bpix_map[i], ann_frame_cumul[i] = res_i if full_output: @@ -748,8 +810,8 @@ def bp_removal_2d(array, cy, cx, fwhm, sig, protect_mask, r_in_std, def cube_fix_badpix_clump(array, bpm_mask=None, correct_only=False, cy=None, cx=None, fwhm=4., sig=4., protect_mask=0, excl_mask=None, half_res_y=False, min_thr=None, - max_nit=15, mad=True, verbose=True, full_output=False, - nproc=1): + max_nit=15, mad=True, bad_values=None, verbose=True, + full_output=False, nproc=1): """ Function to identify and correct clumps of bad pixels. Very fast when a bad pixel map is provided. If a bad pixel map is not provided, the bad pixel @@ -813,6 +875,9 @@ def cube_fix_badpix_clump(array, bpm_mask=None, correct_only=False, cy=None, mad : {False, True}, bool optional If True, the median absolute deviation will be used instead of the standard deviation. + bad_values: list or None, optional + If not None, should correspond to a list of known bad values (e.g. 0). + These pixels will be added to the input bad pixel map. verbose: bool, {False,True}, optional Whether to print the number of bad pixels and number of iterations required for each frame. @@ -820,9 +885,9 @@ def cube_fix_badpix_clump(array, bpm_mask=None, correct_only=False, cy=None, Whether to return as well the cube of bad pixel maps and the cube of defined annuli. nproc: int, optional - This feature is added following ADACS update. Refers to the number of processors - available for calculations. Choosing a number >1 enables multiprocessing for the - correction of frames. + This feature is added following ADACS update. Refers to the number of + processors available for calculations. Choosing a number >1 enables + multiprocessing for the correction of frames. Returns ------- @@ -841,6 +906,12 @@ def cube_fix_badpix_clump(array, bpm_mask=None, correct_only=False, cy=None, assert bpm_mask.shape[-2:] == array.shape[-2:], msg bpm_mask = bpm_mask.astype('bool') + if bad_values is not None: + if bpm_mask is None: + bpm_mask = np.zeros(array_corr.shape, dtype=bool) + for bad in bad_values: + bpm_mask[np.where(array_corr == bad)] = 1 + if correct_only and bpm_mask is None: msg = "Bad pixel map should be provided if correct_only is True." raise ValueError(msg) @@ -930,9 +1001,9 @@ def bp_removal_2d(array_corr, cy, cx, fwhm, sig, protect_mask, bpm_mask_ori, if verbose: msg = "Iteration {}: {} bad pixels identified".format(nit, nbpix_tot) - if bpm_mask_ori is not None: - nbpix_ori = np.sum(bpm_mask_ori) - msg += " ({} new ones)".format(nbpix_tot-nbpix_ori) + # if bpm_mask_ori is not None: + # nbpix_ori = np.sum(bpm_mask_ori) + # msg += " ({} new ones)".format(nbpix_tot-nbpix_ori) if protect_mask: msg += ", {} to be corrected".format(nbpix_tbc) print(msg) @@ -1006,10 +1077,10 @@ def bp_removal_2d(array_corr, cy, cx, fwhm, sig, protect_mask, bpm_mask_ori, cy, cx = frame_center(array) cy = [cy]*n_z cx = [cx]*n_z - elif isinstance(cy, (float, int)) and isinstance(cx, (float, int)): + elif np.isscalar(cy) and np.isscalar(cx): cy = [cy]*n_z cx = [cx]*n_z - if isinstance(fwhm, (float, int)): + if np.isscalar(fwhm): fwhm = [fwhm]*n_z if nproc == 1 or not shared_mem: bpix_map_cumul = np.zeros_like(array_corr) @@ -1073,7 +1144,7 @@ def _mp_clump_slow(args): shm_clump_bpix.unlink() else: - if isinstance(fwhm, (float, int)): + if np.isscalar(fwhm): fwhm_round = int(round(fwhm)) else: fwhm_round = int(np.median(fwhm)) @@ -1314,7 +1385,7 @@ def _res_scaled_images(array, lbdas, fluxes, mask, cy, cx): array_out = np.zeros_like(cube) array_res = np.zeros_like(cube) final_bpm = np.zeros_like(cube) - if isinstance(cy, (float, int)) and isinstance(cx, (float, int)): + if np.isscalar(cy) and np.isscalar(cx): cy = [cy]*n_z cx = [cx]*n_z for i in range(n_z): diff --git a/vip_hci/preproc/recentering.py b/vip_hci/preproc/recentering.py index 214e85cd..3f1cc000 100644 --- a/vip_hci/preproc/recentering.py +++ b/vip_hci/preproc/recentering.py @@ -256,7 +256,7 @@ def frame_shift(array, shift_y, shift_x, imlib='vip-fft', return array_shifted -def cube_shift(cube, shift_y, shift_x, imlib='vip-fft', +def cube_shift(cube, shift_y, shift_x, imlib='vip-fft', interpolation='lanczos4', border_mode='reflect', nproc=None): """ Shifts the X-Y coordinates of a cube or 3D array by x and y values. @@ -293,7 +293,7 @@ def cube_shift(cube, shift_y, shift_x, imlib='vip-fft', if nproc is None: nproc = cpu_count()//2 - + if nproc == 1: cube_out = np.zeros_like(cube) for i in range(cube.shape[0]): @@ -366,7 +366,7 @@ def frame_center_satspots(array, xy, subi_size=19, sigfactor=6, shift=False, ceny, cenx : floats Center Y,X coordinates of the true center. *Only returned if ``shift=True``.* - + Note ---- We are solving a linear system: @@ -647,11 +647,11 @@ def cube_recenter_satspots(array, xy, subi_size=19, sigfactor=6, plot=True, return array_rec -def frame_center_radon(array, cropsize=None, hsize_ini=1., step_ini=0.1, - n_iter=5, tol=0.05, mask_center=None, nproc=None, - satspots_cfg=None, theta_0=0, delta_theta=5, - gauss_fit=True, hpf=True, filter_fwhm=8, imlib='vip-fft', - interpolation='lanczos4', full_output=False, +def frame_center_radon(array, cropsize=None, hsize_ini=1., step_ini=0.1, + n_iter=5, tol=0.05, mask_center=None, nproc=None, + satspots_cfg=None, theta_0=0, delta_theta=5, + gauss_fit=True, hpf=True, filter_fwhm=8, imlib='vip-fft', + interpolation='lanczos4', full_output=False, verbose=True, plot=True, debug=False): """ Finding the center of a broadband (co-added) frame with speckles and satellite spots elongated towards the star (center). We use the Radon @@ -733,17 +733,17 @@ def frame_center_radon(array, cropsize=None, hsize_ini=1., step_ini=0.1, """ from .cosmetics import frame_crop - + if array.ndim != 2: raise TypeError('Input array is not a frame or 2d array') if verbose: start_time = time_ini() - - def _center_radon(array, cropsize=None, hsize=1., step=0.1, - mask_center=None, nproc=None, satspots_cfg=None, - theta_0=0, delta_theta=5, gauss_fit=False, - imlib='vip-fft', interpolation='lanczos4', + + def _center_radon(array, cropsize=None, hsize=1., step=0.1, + mask_center=None, nproc=None, satspots_cfg=None, + theta_0=0, delta_theta=5, gauss_fit=False, + imlib='vip-fft', interpolation='lanczos4', verbose=True, plot=True, debug=False): frame = array.copy() @@ -760,62 +760,62 @@ def _center_radon(array, cropsize=None, hsize=1., step=0.1, if not isinstance(mask_center, int): raise TypeError radint = mask_center - + coords = [(y, x) for y in listyx for x in listyx] cent, _ = frame_center(frame) - + frame = get_annulus_segments(frame, radint, cent-radint, mode="mask")[0] - + if debug: if satspots_cfg is not None: samples = 10 if satspots_cfg == 'x': - theta = np.hstack((np.linspace(start=45-delta_theta, - stop=45+delta_theta, + theta = np.hstack((np.linspace(start=45-delta_theta, + stop=45+delta_theta, num=samples, endpoint=False), - np.linspace(start=135-delta_theta, - stop=135+delta_theta, + np.linspace(start=135-delta_theta, + stop=135+delta_theta, num=samples, endpoint=False), - np.linspace(start=225-delta_theta, - stop=225+delta_theta, + np.linspace(start=225-delta_theta, + stop=225+delta_theta, num=samples, endpoint=False), - np.linspace(start=315-delta_theta, - stop=315+delta_theta, + np.linspace(start=315-delta_theta, + stop=315+delta_theta, num=samples, endpoint=False))) elif satspots_cfg == '+': - theta = np.hstack((np.linspace(start=-delta_theta, - stop=delta_theta, + theta = np.hstack((np.linspace(start=-delta_theta, + stop=delta_theta, num=samples, endpoint=False), - np.linspace(start=90-delta_theta, - stop=90+delta_theta, + np.linspace(start=90-delta_theta, + stop=90+delta_theta, num=samples, endpoint=False), - np.linspace(start=180-delta_theta, - stop=180+delta_theta, + np.linspace(start=180-delta_theta, + stop=180+delta_theta, num=samples, endpoint=False), - np.linspace(start=270-delta_theta, - stop=270+delta_theta, + np.linspace(start=270-delta_theta, + stop=270+delta_theta, num=samples, endpoint=False))) elif satspots_cfg == 'custom': - theta = np.hstack((np.linspace(start=90-theta_0-delta_theta, - stop=90-theta_0+delta_theta, + theta = np.hstack((np.linspace(start=90-theta_0-delta_theta, + stop=90-theta_0+delta_theta, num=samples, endpoint=False), - np.linspace(start=180-theta_0-delta_theta, - stop=180-theta_0+delta_theta, + np.linspace(start=180-theta_0-delta_theta, + stop=180-theta_0+delta_theta, num=samples, endpoint=False), - np.linspace(start=270-theta_0-delta_theta, - stop=270-theta_0+delta_theta, + np.linspace(start=270-theta_0-delta_theta, + stop=270-theta_0+delta_theta, num=samples, endpoint=False), - np.linspace(start=360-theta_0-delta_theta, - stop=360-theta_0+delta_theta, - num=samples, endpoint=False))) + np.linspace(start=360-theta_0-delta_theta, + stop=360-theta_0+delta_theta, + num=samples, endpoint=False))) else: msg = "If not None, satspots_cfg can only be 'x' or '+'." raise ValueError(msg) @@ -828,27 +828,27 @@ def _center_radon(array, cropsize=None, hsize=1., step=0.1, sinogram = radon(frame, theta=theta, circle=True) plot_frames((frame, sinogram)) print(np.sum(np.abs(sinogram[int(cent), :]))) - + if nproc is None: nproc = cpu_count() // 2 # Hyper-threading doubles the # of cores - + if nproc == 1: costf = [] for coord in coords: - res = _radon_costf(frame, cent, radint, coord, satspots_cfg, + res = _radon_costf(frame, cent, radint, coord, satspots_cfg, theta_0, delta_theta, imlib, interpolation) costf.append(res) costf = np.array(costf) elif nproc > 1: - res = pool_map(nproc, _radon_costf, frame, cent, radint, - iterable(coords), satspots_cfg, theta_0, delta_theta, + res = pool_map(nproc, _radon_costf, frame, cent, radint, + iterable(coords), satspots_cfg, theta_0, delta_theta, imlib, interpolation) costf = np.array(res) - + if verbose: msg = 'Done {} radon transform calls distributed in {} processes' print(msg.format(len(coords), nproc)) - + cost_bound = costf.reshape(listyx.shape[0], listyx.shape[0]) if plot: plt.contour(cost_bound, cmap='CMRmap', origin='lower') @@ -857,57 +857,57 @@ def _center_radon(array, cropsize=None, hsize=1., step=0.1, plt.colorbar() plt.grid('off') plt.show() - + if gauss_fit or full_output: # fit a 2d gaussian to the surface - fit_res = fit_2dgaussian(cost_bound-np.amin(cost_bound), crop=False, - threshold=False, sigfactor=1, debug=debug, + fit_res = fit_2dgaussian(cost_bound-np.amin(cost_bound), crop=False, + threshold=False, sigfactor=1, debug=debug, full_output=True) - ## optimal shift -> optimal position + # optimal shift -> optimal position opt_yind = float(fit_res['centroid_y']) opt_xind = float(fit_res['centroid_x']) opt_yshift = -hsize + opt_yind*step opt_xshift = -hsize + opt_xind*step optimy = ori_cent_y - opt_yshift optimx = ori_cent_x - opt_xshift - - ## find uncertainty on centering + + # find uncertainty on centering unc_y = float(fit_res['fwhm_y'])*step unc_x = float(fit_res['fwhm_x'])*step - dyx = (unc_y, unc_x) #np.sqrt(unc_y**2 + unc_x**2) - + dyx = (unc_y, unc_x) # np.sqrt(unc_y**2 + unc_x**2) + # Replace the position found by Gaussian fit if not gauss_fit: # OLD CODE: # argm = np.argmax(costf) # index of 1st max in 1d cost function 'surface' # optimy, optimx = coords[argm] - + # maxima in the 2d cost function surface num_max = np.where(cost_bound == cost_bound.max())[0].shape[0] ind_maximay, ind_maximax = np.where(cost_bound == cost_bound.max()) argmy = ind_maximay[int(np.ceil(num_max/2)) - 1] argmx = ind_maximax[int(np.ceil(num_max/2)) - 1] - y_grid = np.array(coords)[:, 0].reshape(listyx.shape[0], + y_grid = np.array(coords)[:, 0].reshape(listyx.shape[0], listyx.shape[0]) - x_grid = np.array(coords)[:, 1].reshape(listyx.shape[0], + x_grid = np.array(coords)[:, 1].reshape(listyx.shape[0], listyx.shape[0]) optimy = ori_cent_y-y_grid[argmy, 0] # subtract optimal shift optimx = ori_cent_x-x_grid[0, argmx] # subtract optimal shift - + if verbose: print('Cost function max: {}'.format(costf.max())) - #print('Cost function # maxima: {}'.format(num_max)) + # print('Cost function # maxima: {}'.format(num_max)) msg = 'Finished grid search radon optimization. dy={:.3f}, dx={:.3f}' print(msg.format(opt_yshift, opt_xshift)) timing(start_time) - + return optimy, optimx, opt_yshift, opt_xshift, dyx, cost_bound - + # high-pass filtering if requested if hpf: - array = frame_filter_highpass(array, mode='gauss-subt', + array = frame_filter_highpass(array, mode='gauss-subt', fwhm_size=filter_fwhm) - + ori_cent_y, ori_cent_x = frame_center(array) hsize = hsize_ini step = step_ini @@ -916,18 +916,18 @@ def _center_radon(array, cropsize=None, hsize=1., step=0.1, for i in range(n_iter): if verbose: print("*** Iteration {}/{} ***".format(i+1, n_iter)) - res = _center_radon(array, cropsize=cropsize, hsize=hsize, step=step, - mask_center=mask_center, nproc=nproc, - satspots_cfg=satspots_cfg, theta_0=theta_0, - delta_theta=delta_theta, gauss_fit=gauss_fit, - imlib=imlib, interpolation=interpolation, + res = _center_radon(array, cropsize=cropsize, hsize=hsize, step=step, + mask_center=mask_center, nproc=nproc, + satspots_cfg=satspots_cfg, theta_0=theta_0, + delta_theta=delta_theta, gauss_fit=gauss_fit, + imlib=imlib, interpolation=interpolation, verbose=verbose, plot=plot, debug=debug) _, _, y_shift, x_shift, dyx, cost_bound = res - array = frame_shift(array, y_shift, x_shift, imlib=imlib, + array = frame_shift(array, y_shift, x_shift, imlib=imlib, interpolation=interpolation) opt_yshift += y_shift opt_xshift += x_shift - + abs_shift = np.sqrt(y_shift**2 + x_shift**2) if abs_shift < tol: if i == 0: @@ -939,25 +939,25 @@ def _center_radon(array, cropsize=None, hsize=1., step=0.1, print(msg.format(i+1, step)) break # refine box - max_sh = np.amax(np.abs(np.array([y_shift,x_shift]))) + max_sh = np.amax(np.abs(np.array([y_shift, x_shift]))) hsize = 2*max_sh step = hsize/10. optimy = ori_cent_y-opt_yshift optimx = ori_cent_x-opt_xshift - + if verbose: print("Star (x,y) location: {:.2f}, {:.2f}".format(optimx, optimy)) - print("Final cumulative (x,y) shifts: {:.2f}, {:.2f}".format(opt_xshift, + print("Final cumulative (x,y) shifts: {:.2f}, {:.2f}".format(opt_xshift, opt_yshift)) if full_output: return optimy, optimx, dyx, cost_bound else: return optimy, optimx - -def _radon_costf(frame, cent, radint, coords, satspots_cfg=None, theta_0=0, + +def _radon_costf(frame, cent, radint, coords, satspots_cfg=None, theta_0=0, delta_theta=5, imlib='vip-fft', interpolation='lanczos4'): """ Radon cost function used in frame_center_radon(). """ @@ -971,45 +971,45 @@ def _radon_costf(frame, cent, radint, coords, satspots_cfg=None, theta_0=0, endpoint=False) elif satspots_cfg == '+': samples = 10 - theta = np.hstack((np.linspace(start=-delta_theta, stop=delta_theta, + theta = np.hstack((np.linspace(start=-delta_theta, stop=delta_theta, num=samples, endpoint=False), - np.linspace(start=90-delta_theta, + np.linspace(start=90-delta_theta, stop=90+delta_theta, num=samples, endpoint=False), - np.linspace(start=180-delta_theta, + np.linspace(start=180-delta_theta, stop=180+delta_theta, num=samples, endpoint=False), - np.linspace(start=270-delta_theta, + np.linspace(start=270-delta_theta, stop=270+delta_theta, num=samples, endpoint=False))) elif satspots_cfg == 'x': samples = 10 - theta = np.hstack((np.linspace(start=45-delta_theta, + theta = np.hstack((np.linspace(start=45-delta_theta, stop=45+delta_theta, num=samples, endpoint=False), - np.linspace(start=135-delta_theta, + np.linspace(start=135-delta_theta, stop=135+delta_theta, num=samples, endpoint=False), - np.linspace(start=225-delta_theta, + np.linspace(start=225-delta_theta, stop=225+delta_theta, num=samples, endpoint=False), - np.linspace(start=315-delta_theta, + np.linspace(start=315-delta_theta, stop=315+delta_theta, num=samples, endpoint=False))) elif satspots_cfg == 'custom': samples = 10 - theta = np.hstack((np.linspace(start=theta_0-delta_theta, - stop=theta_0+delta_theta, + theta = np.hstack((np.linspace(start=theta_0-delta_theta, + stop=theta_0+delta_theta, num=samples, endpoint=False), - np.linspace(start=theta_0+90-delta_theta, - stop=theta_0+90+delta_theta, + np.linspace(start=theta_0+90-delta_theta, + stop=theta_0+90+delta_theta, num=samples, endpoint=False), - np.linspace(start=theta_0+180-delta_theta, - stop=theta_0+180+delta_theta, + np.linspace(start=theta_0+180-delta_theta, + stop=theta_0+180+delta_theta, num=samples, endpoint=False), - np.linspace(start=theta_0+270-delta_theta, - stop=theta_0+270+delta_theta, - num=samples, endpoint=False))) + np.linspace(start=theta_0+270-delta_theta, + stop=theta_0+270+delta_theta, + num=samples, endpoint=False))) sinogram = radon(frame_shifted_ann, theta=theta, circle=True) costf = np.sum(np.abs(sinogram[int(cent), :])) return costf @@ -1065,13 +1065,13 @@ def cube_recenter_radon(array, full_output=False, verbose=True, imlib='vip-fft', n_frames = array.shape[0] x = np.zeros((n_frames)) y = np.zeros((n_frames)) - dyx = np.zeros((n_frames,2)) + dyx = np.zeros((n_frames, 2)) cy, cx = frame_center(array[0]) array_rec = array.copy() - for i in Progressbar(range(n_frames), desc="Recentering frames...", + for i in Progressbar(range(n_frames), desc="Recentering frames...", verbose=verbose): - res = frame_center_radon(array[i], verbose=False, plot=False, + res = frame_center_radon(array[i], verbose=False, plot=False, imlib=imlib, interpolation=interpolation, full_output=True, **kwargs) y[i] = res[0] @@ -1373,14 +1373,14 @@ def cube_recenter_2dfit(array, xy=None, fwhm=4, subi_size=5, model='gauss', guess parameters for the double gaussian. E.g.: params_2g = {'fwhm_neg': 3.5, 'fwhm_pos': (3.5,4.2), 'theta_neg': 48., 'theta_pos':145., 'neg_amp': 0.5} - + - fwhm_neg: float or tuple with fwhm of neg gaussian - fwhm_pos: can be a tuple for x and y axes of pos gaussian (replaces fwhm) - theta_neg: trigonometric angle of the x axis of the neg gaussian (deg) - theta_pos: trigonometric angle of the x axis of the pos gaussian (deg) - neg_amp: amplitude of the neg gaussian wrt the amp of the positive one - + Note: it is always recommended to provide theta_pos and theta_neg for a better fit. threshold : bool, optional @@ -1833,7 +1833,7 @@ def cube_recenter_via_speckles(cube_sci, cube_ref=None, alignment_iter=5, crop_sz = int(6*fwhm) if not crop_sz % 2: # size should be odd and small, but at least 7 for 2D fit - if crop_sz>7: + if crop_sz > 7: crop_sz -= 1 else: crop_sz += 1 @@ -2028,4 +2028,4 @@ def _fit_2dannulus(array, fwhm=4, crop=False, cent=None, cropsize=15, return mean_y, mean_x else: final_hole_rad = best_rad[i_max, j_max]/fwhm - return mean_y, mean_x, final_hole_rad \ No newline at end of file + return mean_y, mean_x, final_hole_rad