Merge pull request #638 from IainHammond/master

DFT shift fix for new scikit, speed improvement for find_nearest

requirements.txt

@@ -3,7 +3,7 @@ scipy
 astropy
 photutils
 scikit-learn
-scikit-image>0.17.0, <=0.18.3 # unsolved bug in skimage.registration.phase_cross_correlation for versions >0.18.3 up to at least 0.22. Req>0.17 to have the reference mask option implemented in that function.
+scikit-image
 emcee==2.2.1
 nestle
 corner

vip_hci/fm/fakecomp.py

@@ -12,18 +12,13 @@
 import numpy as np
 from scipy import stats
 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
-else:
-    # for photutils version < '0.3' use photutils.centroid_com
-    import photutils.centroid_com as cen_com
+from photutils.centroids import centroid_com
+
 from ..preproc import (cube_crop_frames, frame_shift, frame_crop, cube_shift,
                        frame_rotate)
 from ..var import (frame_center, fit_2dgaussian, fit_2dairydisk, fit_2dmoffat,
@@ -640,7 +635,7 @@ def psf_norm_2d(psf, fwhm, threshold, mask_core, full_output, verbose):
         """Normalize PSF in the 2d case."""
         # we check if the psf is centered and fix it if needed
         cy, cx = frame_center(psf, verbose=False)
-        xcom, ycom = cen_com(psf)
+        xcom, ycom = centroid_com(psf)
         if not (np.allclose(cy, ycom, atol=1e-2) or
                 np.allclose(cx, xcom, atol=1e-2)):
             # first we find the centroid and put it in the center of the array

vip_hci/fm/negfc_fmerit.py

@@ -168,7 +168,6 @@ def chisquare(
     debug: bool, opt
         Whether to debug and plot the post-processed frame after injection of
         the negative fake companion.
-
     Returns
     -------
     out: float
@@ -236,7 +235,7 @@ def chisquare(
         imlib=imlib_sh,
         interpolation=interpolation,
         transmission=transmission,
-        verbose=False,
+        verbose=False
     )
 
     # Perform PCA and extract the zone of interest
@@ -471,7 +470,7 @@ def get_values_optimize(
             collapse=collapse,
             collapse_ifs=collapse_ifs,
             weights=weights,
-            nproc=nproc,
+            nproc=nproc
         )
 
     elif algo == pca_annular or algo == nmf_annular:

vip_hci/fm/negfc_simplex.py

@@ -501,7 +501,7 @@ def firstguess(cube, angs, psfn, planets_xy_coord, ncomp=1, fwhm=4,
                algo_options={}, simplex=True, simplex_options=None, plot=False,
                verbose=True, save=False):
     """Determine a first guess for the position and the flux of a planet using\
-    the negative fake companion techique, as explained in [WER17]_.
+    the negative fake companion technique, as explained in [WER17]_.
 
     This first requires processing the cube without injecting any negative fake
     companion. Once planets or planet candidates are identified, their initial
@@ -532,8 +532,6 @@ def firstguess(cube, angs, psfn, planets_xy_coord, ncomp=1, fwhm=4,
         The number of principal components to use, if the algorithm is PCA. If
         the input cube is 4D, ncomp can be a list of integers, with length
         matching the first dimension of the cube.
-    plsc: float, optional
-        The platescale, in arcsec per pixel.
     fwhm : float, optional
         The FWHM in pixels.
     annulus_width: int, optional

vip_hci/fm/utils_negfc.py

@@ -126,29 +126,23 @@ def find_nearest(array, value, output='index', constraint=None, n=1):
     [output='both']: closest value(s) and index/-ices, respectively.
 
     """
-    if isinstance(array, np.ndarray):
-        pass
-    elif isinstance(array, list):
-        array = np.array(array)
-    else:
-        raise ValueError("Input type for array should be np.ndarray or list.")
-
+    array = np.asarray(array)
     if constraint is None:
-        fm = np.absolute(array-value)
-        idx = fm.argsort()[:n]
+        fm = np.abs(array - value)
+        idx = np.argpartition(fm, n)[:n]
     elif 'floor' in constraint or 'ceil' in constraint:
         indices = np.arange(len(array), dtype=np.int32)
         if 'floor' in constraint:
-            fm = -(array-value)
+            fm = -(array - value)
         else:
-            fm = array-value
+            fm = array - value
         if '=' in constraint:
-            crop_indices = indices[np.where(fm >= 0)]
-            fm = fm[np.where(fm >= 0)]
+            crop_indices = indices[fm >= 0]
+            fm = fm[fm >= 0]
         else:
-            crop_indices = indices[np.where(fm > 0)]
-            fm = fm[np.where(fm > 0)]
-        idx = fm.argsort()[:n]
+            crop_indices = indices[fm > 0]
+            fm = fm[fm > 0]
+        idx = np.argpartition(fm, n)[:n]
         idx = crop_indices[idx]
         if len(idx) == 0:
             msg = "No indices match the constraint ({} w.r.t {:.2f})"

vip_hci/preproc/recentering.py

@@ -31,6 +31,7 @@
            'cube_recenter_via_speckles']
 
 import warnings
+from importlib.metadata import version
 
 import numpy as np
 
@@ -1373,12 +1374,27 @@ def cube_recenter_dft_upsampling(array, upsample_factor=100, subi_size=None,
 
 def _shift_dft(array_rec, array, frnum, upsample_factor, mask, interpolation,
                imlib, border_mode):
-    """Align images using a DFT-based cross-correlation algorithm."""
-    shift_yx = phase_cross_correlation(array_rec[0], array[frnum],
-                                       upsample_factor=upsample_factor,
-                                       reference_mask=mask)
-    y_i, x_i = shift_yx[0]
-
+    """Align images using a DFT-based cross-correlation algorithm, used in
+    cube_recenter_dft_upsampling. See the docstring of
+    skimage.register.phase_cross_correlation for a description of the
+    ``normalization`` parameter which was added in scikit-image 0.19. This
+    should be set to None to maintain the original behaviour of _shift_dft."""
+
+    if version("scikit-image") > "0.18.3":
+        shifts = phase_cross_correlation(array_rec[0], array[frnum],
+                                         upsample_factor=upsample_factor,
+                                         reference_mask=mask,
+                                         normalization=None)
+    else:
+        shifts = phase_cross_correlation(array_rec[0], array[frnum],
+                                         upsample_factor=upsample_factor,
+                                         reference_mask=mask)
+    # from skimage 0.22, phase_cross_correlation returns two more variables
+    # in addition to the array of shifts
+    if len(shifts) == 3:
+        shifts = shifts[0]
+
+    y_i, x_i = shifts
     array_rec_i = frame_shift(array[frnum], shift_y=y_i, shift_x=x_i,
                               imlib=imlib, interpolation=interpolation,
                               border_mode=border_mode)

vip_hci/var/fit_2d.py

@@ -14,13 +14,7 @@
 import numpy as np
 import pandas as pd
 from packaging import version
-import photutils
-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
-else:
-    # for photutils version < '0.3' use photutils.centroid_com
-    import photutils.centroid_com as cen_com
+from photutils.centroids import centroid_com
 from hciplot import plot_frames
 from astropy.modeling import models, fitting
 from astropy.stats import (gaussian_sigma_to_fwhm, gaussian_fwhm_to_sigma,
@@ -232,7 +226,7 @@ def fit_2dgaussian(array, crop=False, cent=None, cropsize=15, fwhmx=4, fwhmy=4,
 
     # Creating the 2D Gaussian model
     init_amplitude = np.ptp(psf_subimage[~bpm_subimage])
-    xcom, ycom = cen_com(psf_subimage)
+    xcom, ycom = centroid_com(psf_subimage)
     gauss = models.Gaussian2D(amplitude=init_amplitude, theta=theta,
                               x_mean=xcom, y_mean=ycom,
                               x_stddev=fwhmx * gaussian_fwhm_to_sigma,
@@ -391,7 +385,7 @@ def fit_2dmoffat(array, crop=False, cent=None, cropsize=15, fwhm=4,
 
     # Creating the 2D Moffat model
     init_amplitude = np.ptp(psf_subimage[~bpm_subimage])
-    xcom, ycom = cen_com(psf_subimage)
+    xcom, ycom = centroid_com(psf_subimage)
     moffat = models.Moffat2D(amplitude=init_amplitude, x_0=xcom, y_0=ycom,
                              gamma=fwhm / 2., alpha=1)
     # Levenberg-Marquardt algorithm
@@ -537,7 +531,7 @@ def fit_2dairydisk(array, crop=False, cent=None, cropsize=15, fwhm=4,
 
     # Creating the 2d Airy disk model
     init_amplitude = np.ptp(psf_subimage[~bpm_subimage])
-    xcom, ycom = cen_com(psf_subimage)
+    xcom, ycom = centroid_com(psf_subimage)
     diam_1st_zero = (fwhm * 2.44) / 1.028
     airy = models.AiryDisk2D(amplitude=init_amplitude, x_0=xcom, y_0=ycom,
                              radius=diam_1st_zero/2.)