Additional features and minor bug fixes

src/vip_hci/fm/fakecomp.py

@@ -824,7 +824,7 @@ def psf_norm_2d(psf, fwhm, threshold, mask_core, full_output, verbose):
 
 def cube_planet_free(planet_parameter, cube, angs, psfn, imlib='vip-fft',
                      interpolation='lanczos4', transmission=None,
-                     radial_gradient=False):
+                     radial_gradient=False, weights=None):
     """Return a cube in which we have injected negative fake companion at the\
     position/flux given by planet_parameter.
 
@@ -863,6 +863,10 @@ def cube_planet_free(planet_parameter, cube, angs, psfn, imlib='vip-fft',
         at the very edge of a physical mask (e.g. ALC) or the effect on the
         light distribution of a marginally extended source near the IWA of the
         coronagraph.
+    weights : 1d array, optional
+        If provided, the negative fake companion fluxes will be scaled according
+        to these weights before injection in the cube. Can reflect changes in
+        the observing conditions throughout the sequence.
 
     Returns
     -------
@@ -906,8 +910,12 @@ def cube_planet_free(planet_parameter, cube, angs, psfn, imlib='vip-fft',
                                                 transmission=transmission,
                                                 radial_gradient=radial_gradient)
         else:
+            if weights is None:
+                flevel = -planet_parameter[i, 2]
+            else:
+                flevel = -planet_parameter[i, 2]*weights
             cpf = cube_inject_companions(cube_temp, psfn, angs, n_branches=1,
-                                         flevel=-planet_parameter[i, 2],
+                                         flevel=flevel,
                                          rad_dists=[planet_parameter[i, 0]],
                                          theta=planet_parameter[i, 1],
                                          imlib=imlib, verbose=False,

src/vip_hci/metrics/snr_source.py

@@ -66,6 +66,9 @@ def snrmap(array, fwhm, approximated=False, plot=False, known_sources=None,
     use2alone: bool, optional
         Whether to use array2 alone to estimate the noise (might be useful to
         estimate the snr of extended disk features).
+    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.
     verbose: bool, optional
         Whether to print timing or not.
     **kwargs : dictionary, optional

src/vip_hci/preproc/badframes.py

@@ -80,7 +80,8 @@ def cube_detect_badfr_pxstats(array, mode='annulus', in_radius=10, width=10,
     n = array.shape[0]
 
     res = cube_basic_stats(array, mode, radius=in_radius,
-                           inner_radius=in_radius, size=width, full_output=True)
+                           inner_radius=in_radius, size=width,
+                           full_output=True)
     if method == 'mean':
         mean_values = res[0]
     else:
@@ -97,12 +98,7 @@ def cube_detect_badfr_pxstats(array, mode='annulus', in_radius=10, width=10,
     top_boundary = np.empty([n])
     bot_boundary = np.empty([n])
     for i in range(n):
-        if mode == 'annulus':
-            i_mean_value = get_annulus_segments(array[i], width=width,
-                                                inner_radius=in_radius,
-                                                mode="val")[0].mean()
-        elif mode == 'circle':
-            i_mean_value = mean_values[i]
+        i_mean_value = mean_values[i]
         top_boundary[i] = mean_smooth[i] + top_sigma*sigma
         bot_boundary[i] = mean_smooth[i] - low_sigma*sigma
         if i_mean_value > top_boundary[i] or i_mean_value < bot_boundary[i]:

src/vip_hci/preproc/recentering.py

@@ -1251,6 +1251,8 @@ def cube_recenter_dft_upsampling(array, upsample_factor=100, subi_size=None,
         start_time = time_ini()
 
     check_array(array, dim=3)
+    if array.shape[0] < 2:
+        raise TypeError("Input array cube should have more than 1 frame.")
     if mask is not None:
         if mask.shape != array.shape[-2:]:
             msg = "If provided, mask should have same shape as frames"
@@ -1296,7 +1298,7 @@ def cube_recenter_dft_upsampling(array, upsample_factor=100, subi_size=None,
     if nproc is None:
         nproc = cpu_count() // 2  # Hyper-threading doubles the # of cores
 
-    if nproc == 1:
+    if nproc == 1 or n_frames < nproc:
         for i in Progressbar(range(1, n_frames),
                              desc="frames", verbose=verbose):
             y[i], x[i], array_rec[i] = _shift_dft(array_rec, array_rec, i,

src/vip_hci/preproc/rescaling.py

@@ -807,6 +807,7 @@ def find_scal_vector(
     interpolation="lanczos4",
     hpf=False,
     fwhm_max=5,
+    check_vec=True,
     **kwargs
 ):
     """
@@ -845,6 +846,10 @@ def find_scal_vector(
     fwhm_max: float, optional
         Maximum FWHM of the PSF across all wavelengths, in pixels. Only used if
         hpf is set to True.
+    check_vec: True, optional
+        Whether to check (and force) the output vector to be >=1. This will set
+        the smallest scaling factor to 1, and scale accordingly the other
+        factors.
     **kwargs: optional
         Optional arguments to the scipy.optimize.minimize function
 
@@ -908,7 +913,8 @@ def find_scal_vector(
         scal_vec[z] = scal_fac
         flux_vec[z] = flux_fac
 
-    scal_vec = check_scal_vector(scal_vec)
+    if check_vec:
+        scal_vec = check_scal_vector(scal_vec)
 
     return scal_vec, flux_vec
 

src/vip_hci/psfsub/pca_fullfr.py

@@ -1032,6 +1032,22 @@ def _adi_rdi_pca(
                 interpolation=interpolation,
                 **rot_options,
             )
+            if smooth is not None:
+                nout = len(gridre)
+                pca_res = cube_filter_lowpass(gridre[0], mode='gauss',
+                                              fwhm_size=smooth, verbose=False)
+                if full_output==True and source_xy is not None:
+                    frame = frame_filter_lowpass(gridre[1], mode='gauss',
+                                                 fwhm_size=smooth)
+                else:
+                    frame = gridre[1]
+                gridre_new = [pca_res, frame]
+                if nout>2:
+                    diff = nout-2
+                    for d in range(diff):
+                        gridre_new.append(gridre[2+d])
+                gridre = tuple(gridre_new)
+
             return gridre
 
 

src/vip_hci/var/filters.py

@@ -52,6 +52,8 @@
     no_opencv = True
 import numpy as np
 from scipy.ndimage import median_filter
+from skimage.filters import median
+from skimage.morphology import disk
 from skimage.restoration import richardson_lucy
 from astropy.convolution import (convolve_fft, convolve, Gaussian2DKernel)
 from astropy.convolution import interpolate_replace_nans as interp_nan
@@ -62,7 +64,7 @@
 
 def cube_filter_iuwt(cube, coeff=5, rel_coeff=1, full_output=False):
     """
-    Isotropic Undecimated Wavelet Transform filtering, as implemented in
+    Isotropic Undecimated Wavelet Transform filtering, as implemented in\
     [KEN15]_ and detailed in [DAB15]_.
 
     Parameters
@@ -203,15 +205,18 @@ def frame_filter_highpass(array, mode, median_size=5, kernel_size=5,
             ``kernel_size``) and using the ``convolve_fft`` Astropy function.
         ``median-subt``
             subtracts a median low-pass filtered version of the image.
+        ``cmedian-subt``
+            subtracts a circular median low-pass filtered version of the image.
         ``gauss-subt``
             subtracts a Gaussian low-pass filtered version of the image.
         ``fourier-butter``
             applies a high-pass 2D Butterworth filter in Fourier domain.
         ``hann``
             uses a Hann window.
 
-    median_size : int, optional
-        Size of the median box for the ``median-subt`` filter.
+    median_size : int or float, optional
+        Size of the median box for the ``median-subt`` and ``cmedian-subt``
+        filters. In the former case, input median_size should be int.
     kernel_size : int, optional
         Size of the Laplacian kernel used in ``laplacian`` mode. It must be an
         positive odd integer value.
@@ -357,6 +362,12 @@ def round_away(x):
                                         median_size=median_size)
         filtered = array - medianed
 
+    elif mode == 'cmedian-subt':
+        # Subtracting the low_pass filtered (median) image from the image itself
+        medianed = frame_filter_lowpass(array, 'cmedian',
+                                        median_size=median_size)
+        filtered = array - medianed
+
     elif mode == 'gauss-subt':
         # Subtracting the low_pass filtered (median) image from the image itself
         gaussed = frame_filter_lowpass(array, 'gauss', fwhm_size=fwhm_size,
@@ -412,10 +423,11 @@ def frame_filter_lowpass(array, mode='gauss', median_size=5, fwhm_size=5,
     ----------
     array : numpy ndarray
         Input array, 2d frame.
-    mode : {'median', 'gauss', 'psf'}, str optional
-        Type of low-pass filtering.
-    median_size : int, optional
-        Size of the median box for filtering the low-pass median filter.
+    mode : {'median', 'cmedian', 'gauss', 'psf'}, str optional
+        Type of low-pass filtering. 'cmedian' stands for circular median.
+    median_size : int or float, optional
+        Size of the median box for the low-pass median filter ('median' and
+        'cmedian' modes). Should be an int for the 'median' mode.
     fwhm_size : float or tuple of 2 floats, optional
         Size of the Gaussian kernel for the low-pass Gaussian filter. If a
         tuple is provided, it should correspond to y and x kernel sizes,
@@ -457,8 +469,6 @@ def frame_filter_lowpass(array, mode='gauss', median_size=5, fwhm_size=5,
     """
     if array.ndim != 2:
         raise TypeError('Input array is not a frame or 2d array.')
-    if not isinstance(median_size, int):
-        raise ValueError('`Median_size` must be integer')
 
     if mask is not None:
         if mode == 'median':
@@ -468,8 +478,14 @@ def frame_filter_lowpass(array, mode='gauss', median_size=5, fwhm_size=5,
             raise TypeError(msg)
 
     if mode == 'median':
+        if not isinstance(median_size, int):
+            raise ValueError('`Median_size` must be integer')
         # creating the low_pass filtered (median) image
         filtered = median_filter(array, median_size, mode='nearest')
+    elif mode == 'cmedian':
+        circle = disk(median_size)
+        # apply median filter with given footprint = structuring element
+        filtered = median(array, footprint=circle)
     elif mode == 'gauss':
         # 2d Gaussian filter
         kernel_sz_y = kernel_sz
@@ -615,7 +631,7 @@ def cube_filter_lowpass(array, mode='gauss', median_size=5, fwhm_size=5,
 
 def frame_deconvolution(array, psf, n_it=30):
     """
-    Iterative image deconvolution following the scikit-image implementation
+    Deconvolve image iteratively following the scikit-image implementation\
     of the Richardson-Lucy algorithm, described in [RIC72]_ and [LUC74]_.
 
     Considering an image that has been convolved by the point spread function