diff --git a/requirements.txt b/requirements.txt
index 09de97ab..b1fbbd90 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -3,7 +3,7 @@ scipy
 astropy
 photutils
 scikit-learn
-scikit-image<=0.18.3
+scikit-image>0.17.0, <=0.18.3
 emcee==2.2.1
 nestle
 corner
diff --git a/tests/test_metrics_detection.py b/tests/test_metrics_detection.py
index d328e6c3..fb1d29d4 100644
--- a/tests/test_metrics_detection.py
+++ b/tests/test_metrics_detection.py
@@ -36,44 +36,44 @@ def get_frame_snrmap(example_dataset_adi):
 
 def test_detection_log(get_frame_snrmap):
     res_frame, coord, fwhm = get_frame_snrmap
-    table = detection(res_frame, fwhm, psf=None, mode='log', plot=False)
+    y, x = detection(res_frame, fwhm, psf=None, mode='log', plot=False)
     check = False
-    for i in range(len(table.y)):
-        if np.allclose(table.y[i], coord[0], atol=2) and \
-           np.allclose(table.x[i], coord[1], atol=2):
+    for i in range(len(y)):
+        if np.allclose(y[i], coord[0], atol=2) and \
+           np.allclose(x[i], coord[1], atol=2):
             check = True
     assert check
 
 
 def test_detection_dog(get_frame_snrmap):
     res_frame, coord, fwhm = get_frame_snrmap
-    table = detection(res_frame, fwhm, psf=None, mode='dog', plot=False)
+    y, x = detection(res_frame, fwhm, psf=None, mode='dog', plot=False)
     check = False
-    for i in range(len(table.y)):
-        if np.allclose(table.y[i], coord[0], atol=2) and \
-           np.allclose(table.x[i], coord[1], atol=2):
+    for i in range(len(y)):
+        if np.allclose(y[i], coord[0], atol=2) and \
+           np.allclose(x[i], coord[1], atol=2):
             check = True
     assert check
 
 
 def test_detection_lpeaks(get_frame_snrmap):
     res_frame, coord, fwhm = get_frame_snrmap
-    table = detection(res_frame, fwhm, psf=None, mode='lpeaks', plot=False)
+    y, x = detection(res_frame, fwhm, psf=None, mode='lpeaks', plot=False)
     check = False
-    for i in range(len(table.y)):
-        if np.allclose(table.y[i], coord[0], atol=2) and \
-           np.allclose(table.x[i], coord[1], atol=2):
+    for i in range(len(y)):
+        if np.allclose(y[i], coord[0], atol=2) and \
+           np.allclose(x[i], coord[1], atol=2):
             check = True
     assert check
 
 
 def test_detection_snrmap(get_frame_snrmap):
     res_frame, coord, fwhm = get_frame_snrmap
-    table = detection(res_frame, fwhm, psf=None, mode='snrmapf',
+    y, x = detection(res_frame, fwhm, psf=None, mode='snrmapf',
                       plot=False, snr_thresh=5, nproc=2)
     check = False
-    for i in range(len(table.y)):
-        if np.allclose(table.y[i], coord[0], atol=2) and \
-           np.allclose(table.x[i], coord[1], atol=2):
+    for i in range(len(y)):
+        if np.allclose(y[i], coord[0], atol=2) and \
+           np.allclose(x[i], coord[1], atol=2):
             check = True
     assert check
diff --git a/vip_hci/__init__.py b/vip_hci/__init__.py
index 10a06fe9..b02c6f08 100644
--- a/vip_hci/__init__.py
+++ b/vip_hci/__init__.py
@@ -1,4 +1,4 @@
-__version__ = "1.3.2"
+__version__ = "1.3.3"
 
 
 from . import preproc
diff --git a/vip_hci/fits/fits.py b/vip_hci/fits/fits.py
index b52b7ab3..da0745c5 100644
--- a/vip_hci/fits/fits.py
+++ b/vip_hci/fits/fits.py
@@ -22,7 +22,7 @@ def open_fits(fitsfilename, n=0, header=False, ignore_missing_end=False,
               precision=np.float32, return_memmap=False, verbose=True,
               **kwargs):
     """
-    Load a fits file into a memory as numpy array.
+    Load a fits file into memory as numpy array.
 
     Parameters
     ----------
diff --git a/vip_hci/fm/negfc_fmerit.py b/vip_hci/fm/negfc_fmerit.py
index f04951e6..1adf3ffb 100644
--- a/vip_hci/fm/negfc_fmerit.py
+++ b/vip_hci/fm/negfc_fmerit.py
@@ -36,8 +36,8 @@ def chisquare(modelParameters, cube, angs, psfs_norm, fwhm, annulus_width,
     Parameters
     ----------
     modelParameters: tuple
-        The model parameters, typically (r, theta, flux). Where flux can be a
-        1d array if cube is 4d.
+        The model parameters: (r, theta, flux) for a 3D input cube, or
+        (r, theta, f1, ..., fN) for a 4D cube with N spectral channels.
     cube: 3d or 4d numpy ndarray
         Input ADI or ADI+IFS cube.
     angs: numpy.array
@@ -317,7 +317,7 @@ def get_values_optimize(cube, angs, ncomp, annulus_width, aperture_radius,
                           nproc=nproc)
 
     elif algo == pca_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)
@@ -327,12 +327,12 @@ def get_values_optimize(cube, angs, ncomp, annulus_width, aperture_radius,
         crop_sz = int(2*np.ceil(radius_int+annulus_width+1))
         if not crop_sz % 2:
             crop_sz += 1
-        if crop_sz < cube.shape[1] and crop_sz < cube.shape[2]:
-            pad = int((cube.shape[1]-crop_sz)/2)
+        if crop_sz < cube.shape[-2] and crop_sz < cube.shape[-1]:
+            pad = int((cube.shape[-2]-crop_sz)/2)
             crop_cube = cube_crop_frames(cube, crop_sz, verbose=False)
         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,
@@ -527,6 +527,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
             crop_cube = cube
 
         pca_res_tmp = pca_annular(crop_cube, angs, radius_int=radius_int,
diff --git a/vip_hci/fm/negfc_simplex.py b/vip_hci/fm/negfc_simplex.py
index 04c0f6de..c16e35e2 100644
--- a/vip_hci/fm/negfc_simplex.py
+++ b/vip_hci/fm/negfc_simplex.py
@@ -609,6 +609,7 @@ def firstguess(cube, angs, psfn, ncomp, planets_xy_coord, fwhm=4,
                                          transmission=transmission,
                                          mu_sigma=mu_sigma, weights=weights,
                                          plot=plot, verbose=verbose, save=save)
+
         r_pre = res_init[0]
         theta_pre = res_init[1]
         f_pre = res_init[2:]
diff --git a/vip_hci/metrics/detection.py b/vip_hci/metrics/detection.py
index ddb5de36..dc2813b7 100644
--- a/vip_hci/metrics/detection.py
+++ b/vip_hci/metrics/detection.py
@@ -80,7 +80,7 @@ def detection(array, fwhm=4, psf=None, mode='lpeaks', bkg_sigma=5,
         Two vectors with the y and x coordinates of the centers of the sources
         (potential planets).
     If full_output is True then a table with all the candidates that passed the
-    2d Gaussian fit constrains and their S/N is returned.
+    2d Gaussian fit constrains and the S/N threshold is returned.
 
     Note
     ----
@@ -377,9 +377,9 @@ def print_abort():
         print(table_full)
 
     if full_output:
-        return table_full
-    else:
         return table
+    else:
+        return yy_final, xx_final
 
 
 def peak_coordinates(obj_tmp, fwhm, approx_peak=None, search_box=None,
diff --git a/vip_hci/preproc/badpixremoval.py b/vip_hci/preproc/badpixremoval.py
index 151a1a8d..e3d00880 100644
--- a/vip_hci/preproc/badpixremoval.py
+++ b/vip_hci/preproc/badpixremoval.py
@@ -29,7 +29,8 @@
 from astropy.stats import sigma_clipped_stats
 from multiprocessing import cpu_count
 from ..stats import clip_array, sigma_filter
-from ..var import frame_center, get_annulus_segments, frame_filter_lowpass
+from ..var import (frame_center, get_annulus_segments, frame_filter_lowpass, 
+                   frame_filter_highpass)
 from ..config import timing, time_ini, Progressbar
 from ..config.utils_conf import pool_map, iterable
 from .rescaling import find_scal_vector, frame_rescaling
@@ -45,9 +46,10 @@
     no_numba = True
 
 
-def frame_fix_badpix_isolated(array, bpm_mask=None, sigma_clip=3, num_neig=5,
-                              size=5, protect_mask=0, cxy=None, mad=False,
-                              ignore_nan=True, verbose=True, full_output=False):
+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):
     """ 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.
@@ -60,6 +62,9 @@ def frame_fix_badpix_isolated(array, bpm_mask=None, sigma_clip=3, num_neig=5,
         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.
+    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.
     sigma_clip : int, optional
         In case no bad pixel mask is provided all the pixels above and below
         sigma_clip*STDDEV will be marked as bad.
@@ -102,6 +107,9 @@ def frame_fix_badpix_isolated(array, bpm_mask=None, sigma_clip=3, num_neig=5,
         raise TypeError('Array is not a 2d array or single frame')
     if size % 2 == 0:
         raise TypeError('Size of the median blur kernel must be an odd integer')
+    if correct_only and bpm_mask is None:
+        msg = "Bad pixel map should be provided if correct_only is True."
+        raise ValueError(msg)
 
     if bpm_mask is not None:
         bpm_mask = bpm_mask.astype('bool')
@@ -120,10 +128,10 @@ def frame_fix_badpix_isolated(array, bpm_mask=None, sigma_clip=3, num_neig=5,
     else:
         cx, cy = cxy
 
-    if bpm_mask is None:
+    if bpm_mask is None or not correct_only:
         ori_nan_mask = np.where(np.isnan(frame))
-        ind = clip_array(frame, sigma_clip, sigma_clip, neighbor=neigh,
-                         num_neighbor=num_neig, mad=mad)
+        ind = clip_array(frame, sigma_clip, sigma_clip, bpm_mask, 
+                         neighbor=neigh, num_neighbor=num_neig, mad=mad)
         bpm_mask = np.zeros_like(frame)
         bpm_mask[ind] = 1
         if ignore_nan:
@@ -149,9 +157,10 @@ def frame_fix_badpix_isolated(array, bpm_mask=None, sigma_clip=3, num_neig=5,
         return array_out
 
 
-def cube_fix_badpix_isolated(array, bpm_mask=None, sigma_clip=3, num_neig=5,
-                             size=5, frame_by_frame=False, protect_mask=0,
-                             cxy=None, mad=False, ignore_nan=True, verbose=True,
+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):
     """ 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
@@ -165,6 +174,9 @@ def cube_fix_badpix_isolated(array, bpm_mask=None, sigma_clip=3, num_neig=5,
         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 per frame 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.
     sigma_clip : int, optional
         In case no bad pixel mask is provided all the pixels above and below
         sigma_clip*STDDEV will be marked as bad.
@@ -211,9 +223,13 @@ def cube_fix_badpix_isolated(array, bpm_mask=None, sigma_clip=3, num_neig=5,
         raise TypeError('Array is not a 3d array or cube')
     if size % 2 == 0:
         raise TypeError('Size of the median blur kernel must be an odd integer')
-
+    if correct_only and bpm_mask is None:
+        msg = "Bad pixel map should be provided if correct_only is True."
+        raise ValueError(msg)
+        
     if bpm_mask is not None:
         bpm_mask = bpm_mask.astype('bool')
+        
 
     if verbose:
         start = time_ini()
@@ -267,10 +283,11 @@ def cube_fix_badpix_isolated(array, bpm_mask=None, sigma_clip=3, num_neig=5,
             final_bpm[i] = res[1]
         count_bp = np.sum(final_bpm)
     else:
-        if bpm_mask is None:
+        if bpm_mask is None or not correct_only:
             ori_nan_mask = np.where(np.isnan(np.nanmean(array, axis=0)))
             ind = clip_array(np.nanmean(array, axis=0), sigma_clip, sigma_clip,
-                             neighbor=neigh, num_neighbor=num_neig, mad=mad)
+                             bpm_mask, neighbor=neigh, num_neighbor=num_neig, 
+                             mad=mad)
             final_bpm = np.zeros_like(array[0], dtype=bool)
             final_bpm[ind] = 1
             if ignore_nan:
@@ -574,10 +591,10 @@ def bp_removal_2d(obj_tmp, cy, cx, fwhm, sig, protect_mask, r_in_std,
         return obj_tmp
 
 
-def cube_fix_badpix_clump(array, bpm_mask=None, cy=None, cx=None, fwhm=4.,
-                          sig=4., protect_mask=0, verbose=True,
-                          half_res_y=False, min_thr=None, max_nit=15, mad=True,
-                          full_output=False):
+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):
     """
     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
@@ -586,8 +603,6 @@ def cube_fix_badpix_clump(array, bpm_mask=None, cy=None, cx=None, fwhm=4.,
     the box is set by the closest odd integer larger than fwhm (to avoid
     accidentally replacing point sources).
 
-
-
     Parameters
     ----------
     array : 3D or 2D array
@@ -596,6 +611,9 @@ def cube_fix_badpix_clump(array, bpm_mask=None, cy=None, cx=None, fwhm=4.,
         Input bad pixel array. Should have same dimenstions as 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.
     cy,cx : float or 1D array. opt
         Vector with approximate y and x coordinates of the star for each channel
         (cube_like), or single 2-elements vector (frame_like). These will be
@@ -613,9 +631,6 @@ def cube_fix_badpix_clump(array, bpm_mask=None, cy=None, cx=None, fwhm=4.,
         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.
-    verbose: bool, {False,True}, optional
-        Whether to print the number of bad pixels and number of iterations
-        required for each frame.
     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
@@ -637,6 +652,9 @@ def cube_fix_badpix_clump(array, bpm_mask=None, cy=None, cx=None, fwhm=4.,
     mad : {False, True}, bool optional
         If True, the median absolute deviation will be used instead of the
         standard deviation.
+    verbose: bool, {False,True}, optional
+        Whether to print the number of bad pixels and number of iterations
+        required for each frame.
     full_output: bool, {False,True}, optional
         Whether to return as well the cube of bad pixel maps and the cube of
         defined annuli.
@@ -656,9 +674,14 @@ def cube_fix_badpix_clump(array, bpm_mask=None, cy=None, cx=None, fwhm=4.,
     if bpm_mask is not None:
         if bpm_mask.shape[-2:] != array.shape[-2:]:
             raise TypeError("Bad pixel map has wrong y/x dimensions.")
+            
+    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(obj_tmp, cy, cx, fwhm, sig, protect_mask, min_thr,
-                      half_res_y, mad, verbose):
+    def bp_removal_2d(obj_tmp, cy, cx, fwhm, sig, protect_mask, bpm_mask_ori,
+                      min_thr, half_res_y, mad, verbose):
         n_x = obj_tmp.shape[1]
         n_y = obj_tmp.shape[0]
 
@@ -694,8 +717,8 @@ def bp_removal_2d(obj_tmp, cy, cx, fwhm, sig, protect_mask, min_thr,
             circl_new = []
 
         # 3/ Create a bad pixel map, by detecting them with clip_array
-        bp = clip_array(obj_tmp, sig, sig, out_good=False, neighbor=True,
-                        num_neighbor=neighbor_box, mad=mad,
+        bp = clip_array(obj_tmp, sig, sig, bpm_mask_ori, out_good=False, 
+                        neighbor=True, num_neighbor=neighbor_box, mad=mad,
                         half_res_y=half_res_y)
         bpix_map = np.zeros_like(obj_tmp)
         bpix_map[bp] = 1
@@ -712,9 +735,9 @@ def bp_removal_2d(obj_tmp, cy, cx, fwhm, sig, protect_mask, min_thr,
             cond1 = obj_tmp > min_thr[0]
             cond2 = obj_tmp < min_thr[1]
             bpix_map[np.where(cond1 & cond2)] = 0
-        nbpix_tot = np.sum(bpix_map)
+        nbpix_tot = int(np.sum(bpix_map))
         bpix_map[circl_new] = 0
-        nbpix_tbc = np.sum(bpix_map)
+        nbpix_tbc = int(np.sum(bpix_map))
         bpix_map_cumul = np.zeros_like(bpix_map)
         bpix_map_cumul[:] = bpix_map[:]
         nit = 0
@@ -728,8 +751,8 @@ def bp_removal_2d(obj_tmp, cy, cx, fwhm, sig, protect_mask, min_thr,
             obj_tmp = sigma_filter(obj_tmp, bpix_map, neighbor_box=neighbor_box,
                                    min_neighbors=nneig, half_res_y=half_res_y,
                                    verbose=verbose)
-            bp = clip_array(obj_tmp, sig, sig, out_good=False, neighbor=True,
-                            num_neighbor=neighbor_box, mad=mad,
+            bp = clip_array(obj_tmp, sig, sig, bpm_mask_ori, out_good=False, 
+                            neighbor=True, num_neighbor=neighbor_box, mad=mad,
                             half_res_y=half_res_y)
             bpix_map = np.zeros_like(obj_tmp)
             bpix_map[bp] = 1
@@ -737,9 +760,9 @@ def bp_removal_2d(obj_tmp, cy, cx, fwhm, sig, protect_mask, min_thr,
                 cond1 = obj_tmp > min_thr[0]
                 cond2 = obj_tmp < min_thr[1]
                 bpix_map[np.where(cond1 & cond2)] = 0
-            nbpix_tot = np.sum(bpix_map)
+            nbpix_tot = int(np.sum(bpix_map))
             bpix_map[circl_new] = 0
-            nbpix_tbc = np.sum(bpix_map)
+            nbpix_tbc = int(np.sum(bpix_map))
             bpix_map_cumul = bpix_map_cumul+bpix_map
 
         if verbose:
@@ -758,12 +781,13 @@ def bp_removal_2d(obj_tmp, cy, cx, fwhm, sig, protect_mask, min_thr,
         return obj_tmp, bpix_map_cumul
 
     if ndims == 2:
-        if bpm_mask is None:
+        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)
             obj_tmp, bpix_map_cumul = bp_removal_2d(obj_tmp, cy, cx, fwhm, sig,
-                                                    protect_mask, min_thr,
-                                                    half_res_y, mad, verbose)
+                                                    protect_mask, bpm_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
@@ -775,7 +799,7 @@ def bp_removal_2d(obj_tmp, cy, cx, fwhm, sig, protect_mask, min_thr,
 
     if ndims == 3:
         n_z = obj_tmp.shape[0]
-        if bpm_mask is None:
+        if bpm_mask is None or not correct_only:
             if cy is None or cx is None:
                 cy, cx = frame_center(array)
                 cy = [cy]*n_z
@@ -792,7 +816,7 @@ def bp_removal_2d(obj_tmp, cy, cx, fwhm, sig, protect_mask, min_thr,
                 obj_tmp[i], bpix_map_cumul[i] = bp_removal_2d(obj_tmp[i], cy[i],
                                                               cx[i], fwhm[i],
                                                               sig, protect_mask,
-                                                              min_thr,
+                                                              bpm_mask, min_thr, 
                                                               half_res_y, mad,
                                                               verbose)
         else:
@@ -1054,7 +1078,7 @@ 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').
-    fwhm: float or 1D array, opt
+    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
         Size of the kernel in pixels for 2D Gaussian and Moffat convolutions.
@@ -1089,7 +1113,6 @@ def cube_fix_badpix_interp(array, bpm_mask, mode='fft', fwhm=4., kernel_sz=None,
     -------
     obj_tmp: 2d or 3d array; the bad pixel corrected frame/cube.
     """
-
     obj_tmp = array.copy()
     ndims = obj_tmp.ndim
     assert ndims == 2 or ndims == 3, "Object is not two or three dimensional.\n"
@@ -1101,11 +1124,11 @@ def cube_fix_badpix_interp(array, bpm_mask, mode='fft', fwhm=4., kernel_sz=None,
         msg = "Warning: no bad pixel found in bad pixel map. "
         msg += "Returning input array as is."
         print(msg)
-        return obj_tmp
+        return array
 
-    ny, nx = obj_tmp.shape[-2:]
+    ny, nx = array.shape[-2:]
     if ndims == 3:
-        nz = obj_tmp.shape[0]
+        nz = array.shape[0]
 
     if bpm_mask.ndim == 2 and ndims == 3:
         master_bpm = np.zeros([nz, ny, nx])
@@ -1113,10 +1136,6 @@ def cube_fix_badpix_interp(array, bpm_mask, mode='fft', fwhm=4., kernel_sz=None,
             master_bpm[z] = bpm_mask[np.newaxis, :, :]
         bpm_mask = master_bpm.copy()
 
-    # below: superseded by mask option?
-    # first replace all bad pixels with NaN values - they will be interpolated
-    # obj_tmp[np.where(bpm_mask)] = np.nan
-
     if half_res_y:
         # squash vertically
         def squash_v(array):
@@ -1126,11 +1145,11 @@ def squash_v(array):
             nny = ny//2
             new_obj_tmp = np.zeros([nny, nx])
             for y in range(nny):
-                new_obj_tmp[y] = array[int(y*2)]
+                new_obj_tmp[y] = obj_tmp[int(y*2)]
             return new_obj_tmp
 
         if ndims == 2:
-            obj_tmp = squash_v(obj_tmp)
+            obj_tmp = squash_v(array)
             bpm_mask = squash_v(bpm_mask)
         else:
             new_obj_tmp = []
@@ -1141,18 +1160,22 @@ def squash_v(array):
             obj_tmp = np.array(new_obj_tmp)
             bpm_mask = np.array(new_bpm_mask)
 
+
     if mode != 'fft':
+        # first replace all bad pixels with NaN values - they will be interpolated
+        obj_tmp[np.where(bpm_mask)] = np.nan
         if ndims == 2:
-            obj_tmp_corr = frame_filter_lowpass(obj_tmp, mode=mode, 
-                                                fwhm_size=fwhm,
-                                                conv_mode='convfft',
+            obj_tmp_corr = frame_filter_lowpass(obj_tmp.copy(), mode=mode, 
+                                                fwhm_size=fwhm, conv_mode='conv',
                                                 kernel_sz=kernel_sz, psf=psf,
-                                                mask=bpm_mask, iterate=True,
+                                                iterate=True,
                                                 half_res_y=half_res_y, **kwargs)
         else:
             obj_tmp_corr = obj_tmp.copy()
             if np.isscalar(fwhm):
                 fwhm = [fwhm]*nz
+            elif len(fwhm) == 2 and len(fwhm) != nz:
+                fwhm = [fwhm]*nz
             if psf is None:
                 psf = [psf]*nz
             elif psf.ndim == 2:
@@ -1161,14 +1184,13 @@ def squash_v(array):
                 raise ValueError("input psf must have same z dimension as array")
             for z in range(nz):
                 obj_tmp_corr[z] = frame_filter_lowpass(obj_tmp[z], mode=mode,
-                                                       fwhm_size=fwhm[z],
-                                                       conv_mode='convfft',
-                                                       kernel_sz=kernel_sz,
-                                                       psf=psf[z], 
-                                                       mask=bpm_mask[z],
-                                                       iterate=True,
-                                                       half_res_y=half_res_y,
-                                                       **kwargs)
+                                                        fwhm_size=fwhm[z],
+                                                        conv_mode='conv',
+                                                        kernel_sz=kernel_sz,
+                                                        psf=psf[z], 
+                                                        iterate=True,
+                                                        half_res_y=half_res_y,
+                                                        **kwargs)
     
         # replace only the bad pixels (obj_tmp_corr is low-pass filtered)
         obj_tmp[np.where(bpm_mask)] = obj_tmp_corr[np.where(bpm_mask)]
@@ -1563,8 +1585,8 @@ def _correct_ann_outliers(obj_tmp, ann_width, sig, med_neig, std_neig,
                                  half_res_y=False)
 
 
-def frame_fix_badpix_fft(array, bpm_mask, nit=500, tol=1, verbose=True, 
-                         full_output=False):
+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]_.
     
@@ -1574,13 +1596,16 @@ def frame_fix_badpix_fft(array, bpm_mask, nit=500, tol=1, verbose=True,
         Input image.
     bpm_mask : 2D ndarray
         Bad pixel map.
-    nit : int
+    nit : int, opt
         Number of iterations.
-    tol: float
+    tol: float, opt
         Tolerance in terms of E_g (see [AAC01]_). The iterative process is 
         stopped if the error E_g gets lower than this tolerance.
+    pad_fac: int or float, opt
+        Padding factor before calculating 2D-FFT.
     verbose: bool
-        Whether to print additional information during processing
+        Whether to print additional information during processing, incl. 
+        progress bar.
     full_output: bool
         Whether to also return the reconstructed estimate f_hat of the input 
         array. 
@@ -1601,15 +1626,17 @@ def frame_fix_badpix_fft(array, bpm_mask, nit=500, tol=1, verbose=True,
 
     # Pad zeros for better results
     ini_y, ini_x = array.shape
-    pad_fac = (int(2*ini_x/ini_y), 2)
+    pad_fac = (int(pad_fac*ini_x/ini_y), pad_fac)
     g = frame_pad(array, pad_fac, keep_parity=False, fillwith=0)
     w = frame_pad(1-bpm_mask, pad_fac, keep_parity=False, fillwith=0)
 
     # Following AAC01 notations:
     g *= w
+    if verbose:
+        start=time_ini()
     G_i = np.fft.fft2(g)
     W = np.fft.fft2(w)
-
+        
     # Initialisation
     dims = g.shape
     ny, nx = dims
@@ -1618,7 +1645,7 @@ def frame_fix_badpix_fft(array, bpm_mask, nit=500, tol=1, verbose=True,
     it = 0
     Eg = tol + 1
 
-    while it < nit and Eg > tol:
+    for it in Progressbar(range(nit), desc="iterative bad pixel correction"):
         # 1. select line as max(G_i) and infer conjugate coordinates
         ind = np.unravel_index(np.argmax(np.abs(G_i.real[:, 0: nx // 2])),
                                (ny, nx // 2))
@@ -1668,12 +1695,15 @@ def frame_fix_badpix_fft(array, bpm_mask, nit=500, tol=1, verbose=True,
         # 5. Calculate new error - to check if still larger than tolerance
         Eg = np.sum(np.power(np.abs(G_i.flatten()), 2))/npix
 
-        it += 1
+        if Eg < tol:
+            break
         
     if verbose:
         msg = "FFT-interpolation terminated after {} iterations (Eg={})"
         print(msg.format(it, Eg))
         
+    if verbose:  
+        timing(start)
     bpix_corr = g + np.fft.ifft2(F_est).real * (1-w)
 
     # crop zeros to return to initial size
diff --git a/vip_hci/preproc/recentering.py b/vip_hci/preproc/recentering.py
index e08f296e..9fdb9d5b 100644
--- a/vip_hci/preproc/recentering.py
+++ b/vip_hci/preproc/recentering.py
@@ -1188,8 +1188,8 @@ def cube_recenter_dft_upsampling(array, center_fr1=None, negative=False,
     n_frames, sizey, sizex = array.shape
     if subi_size is not None:
         if center_fr1 is None:
-            print('`cx_1` or `cy_1` not be provided')
-            print('Using the coordinated of the 1st frame center for '
+            print('`cx_1` or `cy_1` not provided')
+            print('Using the coordinates of the 1st frame center for '
                   'the Gaussian 2d fit')
             cy_1, cx_1 = frame_center(array[0])
         else:
diff --git a/vip_hci/psfsub/pca_local.py b/vip_hci/psfsub/pca_local.py
index b66a25e7..a2c33dc7 100644
--- a/vip_hci/psfsub/pca_local.py
+++ b/vip_hci/psfsub/pca_local.py
@@ -181,7 +181,7 @@ def pca_annular(cube, angle_list, cube_ref=None, scale_list=None, radius_int=0,
         See the documentation of the ``vip_hci.preproc.frame_rotate`` function.
     collapse : {'median', 'mean', 'sum', 'trimmean'}, str optional
         Sets the way of collapsing the frames for producing a final image.
-    collapse_ifs : {'median', 'mean', 'sum', 'trimmean'}, str optional
+    collapse_ifs : {'median', 'mean', 'sum', 'trimmean', 'absmean'}, str opt
         Sets how spectral residual frames should be combined to produce an
         mSDI image.
     ifs_collapse_range: str 'all' or tuple of 2 int
@@ -265,7 +265,7 @@ def pca_annular(cube, angle_list, cube_ref=None, scale_list=None, radius_int=0,
                                    svd_mode, nproc, min_frames_lib,
                                    max_frames_lib, tol, scaling, imlib,
                                    interpolation, collapse, True, verbose,
-                                   cube_ref_tmp, weights, cube_sig,
+                                   cube_ref_tmp, theta_init, weights, cube_sig,
                                    **rot_options)
             cube_out.append(res_pca[0])
             cube_der.append(res_pca[1])
@@ -290,13 +290,10 @@ def pca_annular(cube, angle_list, cube_ref=None, scale_list=None, radius_int=0,
         fwhm = int(np.round(np.mean(fwhm)))
         n_annuli = int((y_in / 2 - radius_int) / asize)
 
-        if scale_list is None:
-            raise ValueError('Scaling factors vector must be provided')
-        else:
-            if np.array(scale_list).ndim > 1:
-                raise ValueError('Scaling factors vector is not 1d')
-            if not scale_list.shape[0] == z:
-                raise ValueError('Scaling factors vector has wrong length')
+        if np.array(scale_list).ndim > 1:
+            raise ValueError('Scaling factors vector is not 1d')
+        if not scale_list.shape[0] == z:
+            raise ValueError('Scaling factors vector has wrong length')
 
         if not isinstance(ncomp, tuple):
             raise TypeError("`ncomp` must be a tuple of two integers when "
@@ -440,7 +437,6 @@ def _pca_sdi_fr(fr, scal, radius_int, fwhm, asize, n_segments, delta_sep, ncomp,
                         interpolation=interpolation, collapse=collapse)
     return frame_desc
 
-
 def _pca_adi_rdi(cube, angle_list, radius_int=0, fwhm=4, asize=2, n_segments=1,
                  delta_rot=1, ncomp=1, svd_mode='lapack', nproc=None,
                  min_frames_lib=2, max_frames_lib=200, tol=1e-1, scaling=None,
@@ -462,7 +458,7 @@ def _pca_adi_rdi(cube, angle_list, radius_int=0, fwhm=4, asize=2, n_segments=1,
 
     if isinstance(delta_rot, tuple):
         delta_rot = np.linspace(delta_rot[0], delta_rot[1], num=n_annuli)
-    elif isinstance(delta_rot, (int, float)):
+    elif np.isscalar(delta_rot):
         delta_rot = [delta_rot] * n_annuli
 
     if isinstance(n_segments, int):
diff --git a/vip_hci/stats/clip_sigma.py b/vip_hci/stats/clip_sigma.py
index 59baaa6b..b9744f9f 100644
--- a/vip_hci/stats/clip_sigma.py
+++ b/vip_hci/stats/clip_sigma.py
@@ -191,8 +191,9 @@ def _sigma_filter(frame_tmp, bpix_map, neighbor_box=3, min_neighbors=3,
                                    min_neighbors=3, verbose=False)
 
 
-def clip_array(array, lower_sigma, upper_sigma, out_good=False, neighbor=False,
-               num_neighbor=3, mad=False, half_res_y=False):
+def clip_array(array, lower_sigma, upper_sigma, bpm_mask_ori=None, 
+               out_good=False, neighbor=False, num_neighbor=3, mad=False, 
+               half_res_y=False):
     """Sigma clipping for detecting outlying values in 2d array. If the
     parameter 'neighbor' is True the clipping can be performed in a local patch
     around each pixel, whose size depends on 'neighbor' parameter.
@@ -205,8 +206,13 @@ def clip_array(array, lower_sigma, upper_sigma, out_good=False, neighbor=False,
         Value for sigma, lower boundary.
     upper_sigma : float
         Value for sigma, upper boundary.
+    bpm_mask_ori : 2d numpy ndarray or None
+        Known (static) bad pixels. Additional bad pixels will be identified,
+        on top of those ones. They are not considered as good neighbours during
+        sigma-clipping.
     out_good : bool, optional
-        For choosing different outputs.
+        For choosing different outputs. True to return indices of good pixels, 
+        False for indices of bad pixels.
     neighbor : bool, optional
         For clipping over the median of the contiguous pixels.
     num_neighbor : int, optional
@@ -233,12 +239,12 @@ def clip_array(array, lower_sigma, upper_sigma, out_good=False, neighbor=False,
 
     """
 
-    def _clip_array(array, lower_sigma, upper_sigma, out_good=False,
-                    neighbor=False, num_neighbor=3, mad=False,
+    def _clip_array(array, lower_sigma, upper_sigma, bpm_mask_ori, 
+                    out_good=False, neighbor=False, num_neighbor=3, mad=False,
                     half_res_y=False):
         """Sigma clipping for detecting outlying values in 2d array. If the
-        parameter 'neighbor' is True the clipping can be performed in a local patch
-        around each pixel, whose size depends on 'neighbor' parameter.
+        parameter 'neighbor' is True the clipping can be performed in a local 
+        patch around each pixel, whose size depends on 'neighbor' parameter.
 
         Parameters
         ----------
@@ -248,6 +254,10 @@ def _clip_array(array, lower_sigma, upper_sigma, out_good=False,
             Value for sigma, lower boundary.
         upper_sigma : float
             Value for sigma, upper boundary.
+        bpm_mask_ori : 2d numpy ndarray or None
+            Known (e.g. static) bad pixel map. Additional bad pixels will be 
+            identified, on top of these ones. They are not considered as good 
+            neighbours during sigma-clipping.
         out_good : bool, optional
             For choosing different outputs.
         neighbor : bool optional
@@ -274,14 +284,20 @@ def _clip_array(array, lower_sigma, upper_sigma, out_good=False,
         """
         if array.ndim != 2:
             raise TypeError("Input array is not two dimensional (frame)\n")
-
+        if bpm_mask_ori is None:
+            gpm_ori = np.ones(array.shape)
+        else:
+            gpm_ori = np.ones(array.shape)-bpm_mask_ori
+        
         ny, nx = array.shape
-        bpm = array.copy()
-        gpm = array.copy()
+        bpm = np.ones(array.shape)
+        gpm = np.zeros(array.shape)
 
         if neighbor and num_neighbor:
             for y in range(ny):
                 for x in range(nx):
+                    if not gpm_ori[y,x]:
+                        continue
                     # 0/ Determine the box around each pixel
                     half_box_x = int(np.floor(num_neighbor/2.))
                     if half_res_y:
@@ -313,7 +329,9 @@ def _clip_array(array, lower_sigma, upper_sigma, out_good=False,
 
                     sub_arr = array[y-hbox_b:y+hbox_t+1,
                                     x-hbox_l:x+hbox_r+1]
-                    neighbours = sub_arr.flatten()
+                    gp_arr = gpm_ori[y-hbox_b:y+hbox_t+1,
+                                     x-hbox_l:x+hbox_r+1]
+                    neighbours = sub_arr[np.where(gp_arr)].flatten()
                     neigh_list = []
                     remove_itself = True
                     for i in range(neighbours.shape[0]):
@@ -351,14 +369,15 @@ def _clip_array(array, lower_sigma, upper_sigma, out_good=False,
         else:
             bad = np.where(bpm)
             return bad
-
+    
     if no_numba:
-        return _clip_array(array, lower_sigma, upper_sigma, out_good=out_good,
-                           neighbor=neighbor, num_neighbor=num_neighbor,
-                           mad=mad, half_res_y=half_res_y)
+        return _clip_array(array, lower_sigma, upper_sigma, bpm_mask_ori, 
+                           out_good=out_good, neighbor=neighbor, 
+                           num_neighbor=num_neighbor, mad=mad, 
+                           half_res_y=half_res_y)
     else:
         _clip_array_numba = njit(_clip_array)
-        return _clip_array_numba(array, lower_sigma, upper_sigma,
+        return _clip_array_numba(array, lower_sigma, upper_sigma, bpm_mask_ori,
                                  out_good=out_good, neighbor=neighbor,
                                  num_neighbor=num_neighbor, mad=mad,
                                  half_res_y=half_res_y)
\ No newline at end of file
diff --git a/vip_hci/var/filters.py b/vip_hci/var/filters.py
index a5c1024b..7b854219 100644
--- a/vip_hci/var/filters.py
+++ b/vip_hci/var/filters.py
@@ -417,8 +417,9 @@ def frame_filter_lowpass(array, mode='gauss', median_size=5, fwhm_size=5,
         Type of low-pass filtering.
     median_size : int, optional
         Size of the median box for filtering the low-pass median filter.
-    fwhm_size : float, optional
-        Size of the Gaussian kernel for the low-pass Gaussian filter.
+    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 correspon
     conv_mode : {'conv', 'convfft'}, str optional
         'conv' uses the multidimensional gaussian filter from scipy.ndimage and
         'convfft' uses the fft convolution with a 2d Gaussian kernel.
@@ -432,7 +433,7 @@ def frame_filter_lowpass(array, mode='gauss', median_size=5, fwhm_size=5,
     mask: numpy ndarray, optional
         Binary mask indicating where the low-pass filtered image should be
         interpolated with astropy.convolution. This option can be useful if the
-        low-pass filtered image is aimed to capture low-spatial frequency sky
+        low-pass filtered image is used to capture low-spatial frequency sky
         signal, while avoiding a stellar halo (set to one in the binary mask).
         Note: only works with Gaussian kernel or PSF convolution.
     iterate: bool, opt
@@ -471,16 +472,27 @@ def frame_filter_lowpass(array, mode='gauss', median_size=5, fwhm_size=5,
         filtered = median_filter(array, median_size, mode='nearest')
     elif mode == 'gauss':
         # 2d Gaussian filter
-        sigma = fwhm_size * gaussian_fwhm_to_sigma
         kernel_sz_y = kernel_sz
-        if half_res_y:
-            sigma_y = max(1, sigma//2)
+        if np.isscalar(fwhm_size):
+            sigma = fwhm_size * gaussian_fwhm_to_sigma
+            sigma_y = sigma
+        else:
+            if len(fwhm_size) != 2:
+                msg = "If not a scalar, fwhm_size must be of length 2"
+                raise TypeError(msg)
+            sigma_y = fwhm_size[0] * gaussian_fwhm_to_sigma
+            sigma = fwhm_size[1] * gaussian_fwhm_to_sigma
             if kernel_sz is not None:
-                kernel_sz_y = kernel_sz//2
+                kernel_sz_y = int(kernel_sz*fwhm_size[0]/fwhm_size[1])
+                if kernel_sz_y % 2 != kernel_sz % 2:
+                    kernel_sz_y += 1
+            
+        if half_res_y:
+            sigma_y = max(1, sigma_y//2)
+            if kernel_sz_y is not None:
+                kernel_sz_y = kernel_sz_y//2
                 if kernel_sz_y % 2 != kernel_sz % 2:
                     kernel_sz_y += 1
-        else:
-            sigma_y = sigma
 
         if conv_mode == 'conv':
             filtered = convolve(array, Gaussian2DKernel(x_stddev=sigma,
diff --git a/vip_hci/var/fit_2d.py b/vip_hci/var/fit_2d.py
index 367e9eeb..48000b9d 100644
--- a/vip_hci/var/fit_2d.py
+++ b/vip_hci/var/fit_2d.py
@@ -211,7 +211,8 @@ def fit_2dgaussian(array, crop=False, cent=None, cropsize=15, fwhmx=4, fwhmy=4,
 
         imside = array.shape[0]
         psf_subimage, suby, subx = get_square(array, min(cropsize, imside),
-                                              ceny, cenx, position=True)
+                                              ceny, cenx, position=True, 
+                                              verbose=False)
     else:
         psf_subimage = array.copy()