Merge pull request #568 from VChristiaens/master

Final JOSS paper edits + additional minor features

.github/workflows/ci.yml

@@ -17,7 +17,7 @@ jobs:
     steps:
     - uses: actions/checkout@v3
     - name: Set up Python ${{ matrix.python-version }}
-      uses: actions/setup-python@v2
+      uses: actions/setup-python@v4
       with:
         python-version: ${{ matrix.python-version }}
     - name: Install dependencies

paper.md

@@ -120,10 +120,10 @@ options (e.g. PCA-based sky subtraction, image centering or bad frame trimming).
 processing pipeline including similar options as ``VIP`` regarding preprocessing 
 [@pynpoint:2019]. Nonetheless, the PCA implementation in ``VIP`` offers a much
 wider diversity of options, such as the possibility to carry it out in 
-concentric annuli, and considering a parallactic angle threshold when creating 
+concentric annuli, and to consider a parallactic angle threshold when creating 
 the PCA library. Depending on the high-contrast imaging dataset at hand, 
 different post-processing methods and reduction parameters can lead to better 
-speckle suppression, hence help with the detection of fainter companions 
+speckle suppression, helping with the detection of fainter companions 
 [@Dahlqvist:2021]. In that regard, ``VIP`` is thus better equipped than other 
 existing toolkits. It is also worth mentioning that FFT-based methods are 
 implemented in ``VIP`` (default option) for all image operations (rotation, 
@@ -165,11 +165,11 @@ of debris disks [@Milli:2017b; @Milli:2019], or the development of new
 high-contrast imaging algorithms 
 [@Gomez:2018;  @Dahlqvist:2020;  @Pairet:2021;  @Dahlqvist:2021]. 
 
-Given the rapid expansion of ``VIP``, we summarize here all novelties that were 
+Given the rapid expansion of ``VIP``, we summarize here all new features that were 
 brought to the package over the past five years. Specifically, the rest of this 
 manuscript summarizes all major changes since v0.7.0 [@Gomez:2017], that are 
-included in the latest release of ``VIP`` (v1.3.4). At a structural level, 
-``VIP`` underwent a major change since version v1.1.0, which aimed to migrate 
+included in the latest release of ``VIP`` (v1.3.5). At a structural level, 
+``VIP`` underwent a major change since version v1.1.0, which migrated it
 towards a more streamlined and easy-to-use architecture. The package now 
 revolves around five major modules (`fm`, `invprob`, `metrics`, `preproc` and 
 `psfsub`, as described above) complemented by four additional modules containing 

vip_hci/__init__.py

@@ -1,4 +1,4 @@
-__version__ = "1.3.4"
+__version__ = "1.3.5"
 
 
 from . import preproc

vip_hci/metrics/snr_source.py

@@ -9,6 +9,7 @@
 __author__ = 'Carlos Alberto Gomez Gonzalez, O. Absil @ ULg, V. Christiaens'
 __all__ = ['snr',
            'snrmap',
+           'indep_ap_centers',
            'significance',
            'frame_report']
 
@@ -216,6 +217,65 @@ 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):
+
+    sourcex, sourcey = source_xy
+    centery, centerx = frame_center(array)
+    sep = dist(centery, centerx, float(sourcey), float(sourcex))
+    theta_0 = np.rad2deg(np.arctan2(sourcey - centery, sourcex - centerx))
+
+    if exclude_theta_range is not None:
+        exc_theta_range = list(exclude_theta_range)
+
+    if not sep > (fwhm / 2) + 1:
+        raise RuntimeError('`source_xy` is too close to the frame center')
+
+    # sens = 'clock'  # counterclock
+    # assumes clockwise rotation when building test apertures
+    # change sign and conditions if counterclockwise
+    sign = -1
+    if exclude_theta_range is not None:
+        if theta_0 > exc_theta_range[0] and theta_0 < exc_theta_range[1]:
+            exc_theta_range[0] += 360
+        while theta_0 < exc_theta_range[1]:
+            theta_0 += 360
+    theta = theta_0
+
+    angle = np.arcsin(fwhm / 2. / sep) * 2
+    number_apertures = int(np.floor(2 * np.pi / angle))
+    yy = []
+    xx = []
+    yy_all = np.zeros(number_apertures)
+    xx_all = np.zeros(number_apertures)
+    cosangle = np.cos(angle)
+    sinangle = np.sin(angle)
+    xx.append(sourcex - centerx)
+    yy.append(sourcey - centery)
+    xx_all[0] = sourcex - centerx
+    yy_all[0] = sourcey - centery
+
+    for i in range(number_apertures - 1):
+        xx_all[i + 1] = cosangle * xx_all[i] - sign * sinangle * yy_all[i]
+        yy_all[i + 1] = cosangle * yy_all[i] + sign * sinangle * xx_all[i]
+        theta += sign * np.rad2deg(angle)
+        if exclude_negative_lobes and (i == 0 or i == number_apertures - 2):
+            continue
+        if exclude_theta_range is None:
+            xx.append(cosangle * xx_all[i] - sign * sinangle * yy_all[i])
+            yy.append(cosangle * yy_all[i] + sign * sinangle * xx_all[i])
+        else:
+            if theta < exc_theta_range[0] or theta > exc_theta_range[1]:
+                xx.append(cosangle * xx_all[i] - sign * sinangle * yy_all[i])
+                yy.append(cosangle * yy_all[i] + sign * sinangle * xx_all[i])
+
+    xx = np.array(xx)
+    yy = np.array(yy)
+
+    xx += centerx
+    yy += centery
+
+    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, 
@@ -295,59 +355,10 @@ def snr(array, source_xy, fwhm, full_output=False, array2=None, use2alone=False,
             raise TypeError('`array2` has not the same shape as input array')
 
     sourcex, sourcey = source_xy
-    centery, centerx = frame_center(array)
-    sep = dist(centery, centerx, float(sourcey), float(sourcex))
-    theta_0 = np.rad2deg(np.arctan2(sourcey-centery,sourcex-centerx))
-
-    if exclude_theta_range is not None:
-        exc_theta_range = list(exclude_theta_range)   
-
-    if not sep > (fwhm/2)+1:
-        raise RuntimeError('`source_xy` is too close to the frame center')
 
-    #sens = 'clock'  # counterclock
-    # assumes clockwise rotation when building test apertures
-    # change sign and conditions if counterclockwise
-    sign = -1
-    if exclude_theta_range is not None:
-        if theta_0 > exc_theta_range[0] and theta_0 < exc_theta_range[1]:
-            exc_theta_range[0] += 360
-        while theta_0<exc_theta_range[1]:
-            theta_0 += 360
-    theta = theta_0
+    yy, xx = indep_ap_centers(array, source_xy, fwhm, exclude_negative_lobes,
+                              exclude_theta_range)
 
-    angle = np.arcsin(fwhm/2./sep)*2
-    number_apertures = int(np.floor(2*np.pi/angle))
-    yy = []
-    xx = []
-    yy_all = np.zeros(number_apertures)
-    xx_all = np.zeros(number_apertures)
-    cosangle = np.cos(angle)
-    sinangle = np.sin(angle)
-    xx.append(sourcex - centerx)
-    yy.append(sourcey - centery)
-    xx_all[0] = sourcex - centerx
-    yy_all[0] = sourcey - centery
-
-    for i in range(number_apertures-1):
-        xx_all[i+1]=cosangle*xx_all[i] - sign*sinangle*yy_all[i]
-        yy_all[i+1]=cosangle*yy_all[i] + sign*sinangle*xx_all[i]
-        theta += sign*np.rad2deg(angle)
-        if exclude_negative_lobes and (i==0 or i==number_apertures-2):
-            continue
-        if exclude_theta_range is None:
-            xx.append(cosangle*xx_all[i] - sign*sinangle*yy_all[i])
-            yy.append(cosangle*yy_all[i] + sign*sinangle*xx_all[i])
-        else:
-            if theta<exc_theta_range[0] or theta>exc_theta_range[1]:
-                xx.append(cosangle*xx_all[i] - sign*sinangle*yy_all[i])
-                yy.append(cosangle*yy_all[i] + sign*sinangle*xx_all[i])
-
-    xx=np.array(xx)
-    yy=np.array(yy)
-
-    xx += centerx
-    yy += centery
     rad = fwhm/2.
 
     apertures = photutils.CircularAperture(