Repository: lightkurve/lightkurve

Issue:
`centroid_quadratic()` incorrectly failed in a boundary case, when `mask` is specified and fluxes are all negative
<!-- Fill in the information below before opening an issue. -->

#### Problem description
<!-- Provide a clear and concise description of the issue. -->

`lk.utils.centroid_quadratic()` would fail and return `(nan, nan)` in a specific boundary case:
- `mask` argument is specified, and
- the flux in `data` is all negative.

While all negative flux should not happen in normal flux data, they could possibly arise when the data is from a difference flux.


#### Example

See: https://gist.github.com/orionlee/fb9a855c3f4fe935ef653d72910768ac


#### Expected behavior
<!-- Describe the behavior you expected and how it differs from the behavior observed in the example. -->
It should still return a reasonable centroid.


#### Environment

-  lightkurve version (e.g. 1.0b6): 2.4.2


Discussion:

Patch:
diff --git a/src/lightkurve/periodogram.py b/src/lightkurve/periodogram.py
index 2c7c9932..bb3a4211 100644
--- a/src/lightkurve/periodogram.py
+++ b/src/lightkurve/periodogram.py
@@ -650,7 +650,8 @@ class LombScarglePeriodogram(Periodogram):
         ls_method="fast",
         **kwargs
     ):
-        """Creates a `Periodogram` from a LightCurve using the Lomb-Scargle method.
+        """Creates a `Periodogram` from a LightCurve using the Lomb-Scargle method in
+        `astropy`'s `~astropy.timeseries.LombScargle`.
 
         By default, the periodogram will be created for a regular grid of
         frequencies from one frequency separation to the Nyquist frequency,
@@ -768,9 +769,11 @@ class LombScarglePeriodogram(Periodogram):
             (`'amplitude'`).
         ls_method : str
             Default: `'fast'`. Passed to the `method` keyword of
-            `astropy.stats.LombScargle()`.
+            `astropy.timeseries.LombScargle()`.
         kwargs : dict
-            Keyword arguments passed to `astropy.stats.LombScargle()`
+            Keyword arguments passed to
+            `LombScargle() <astropy.timeseries.LombScargle>`
+
 
         Returns
         -------
@@ -1019,7 +1022,8 @@ class BoxLeastSquaresPeriodogram(Periodogram):
 
     @staticmethod
     def from_lightcurve(lc, **kwargs):
-        """Creates a `Periodogram` from a LightCurve using the Box Least Squares (BLS) method.
+        """Creates a `Periodogram` from a LightCurve using the Box Least Squares (BLS) method
+        in `astropy`'s `~astropy.timeseries.BoxLeastSquares`.
 
         Parameters
         ----------
diff --git a/src/lightkurve/utils.py b/src/lightkurve/utils.py
index f63897ba..fe965d3f 100644
--- a/src/lightkurve/utils.py
+++ b/src/lightkurve/utils.py
@@ -627,10 +627,20 @@ def centroid_quadratic(data, mask=None):
     """
     if isinstance(data, u.Quantity):
         data = data.value
+
+    if np.issubdtype(data.dtype, int):
+        # multiple code paths below require data be of float type
+        # proactively convert int to float once and for all.
+        data = data.astype(float)
+
     # Step 1: identify the patch of 3x3 pixels (z_)
     # that is centered on the brightest pixel (xx, yy)
     if mask is not None:
-        data = data * mask
+        # mask handling.
+        # Issue 1401 demonstrates that using 'data' to find the max will break when all flux is negative
+        # set masked pixels NaN (instead of 0) to resolve it.
+        data = data.copy()
+        data[~mask] = np.nan
     arg_data_max = np.nanargmax(data)
     yy, xx = np.unravel_index(arg_data_max, data.shape)
     # Make sure the 3x3 patch does not leave the TPF bounds
@@ -644,6 +654,13 @@ def centroid_quadratic(data, mask=None):
         xx = data.shape[1] - 2
 
     z_ = data[yy - 1 : yy + 2, xx - 1 : xx + 2]
+    if np.any(np.isnan(z_)):
+        # handle edge case the 3X3 patch has NaN
+        # Need some finite value for NaN pixels for the
+        # quadratic fit below: use the mean of the 3x3 patch
+        # to reduce the skew
+        z_ = z_.copy()
+        z_[np.isnan(z_)] = np.nanmean(z_)
 
     # Next, we will fit the coefficients of the bivariate quadratic with the
     # help of a design matrix (A) as defined by Eqn 20 in Vakili & Hogg


Test patch:
diff --git a/tests/test_utils.py b/tests/test_utils.py
index d58da191..601b36e5 100644
--- a/tests/test_utils.py
+++ b/tests/test_utils.py
@@ -165,25 +165,66 @@ def test_centroid_quadratic():
     assert np.isclose(row, 5) & np.isclose(col, 1.5)
 
 
-def test_centroid_quadratic_robustness():
-    """Test quadratic centroids in edge cases; regression test for #610."""
+# create a mask that'd mask out first 2 rows
+a_mask = np.full((5, 5), True, dtype=bool)
+a_mask[0:2, :] = False
+
+@pytest.mark.parametrize("data_dtype, mask", [
+    (float, None),
+    (float, a_mask),
+    (int, None),
+    (int, a_mask),
+    ])
+def test_centroid_quadratic_robustness(data_dtype, mask):
+    """Test quadratic centroids in edge cases; regression test for #610, etc."""
     # Brightest pixel in upper left
-    data = np.zeros((5, 5))
+    data = np.zeros((5, 5), dtype=data_dtype)
     data[0, 0] = 1
-    centroid_quadratic(data)
+    centroid_quadratic(data, mask=mask)
+    col, row = centroid_quadratic(data, mask=mask)
+    if mask is None:
+        assert np.isfinite(col) & np.isfinite(row)
+    else:
+        # the mask (of top 2 rows) would make the eligible pixels to be uniformly at 0
+        # no centroid can be calculated in such case.
+        assert np.isnan(col) & np.isnan(row)
 
     # Brightest pixel in bottom right
-    data = np.zeros((5, 5))
+    data = np.zeros((5, 5), dtype=data_dtype)
     data[-1, -1] = 1
-    centroid_quadratic(data)
+    col, row = centroid_quadratic(data, mask=mask)
+    assert np.isfinite(col) & np.isfinite(row)
 
     # Data contains a NaN
-    data = np.zeros((5, 5))
-    data[0, 0] = np.nan
-    data[-1, -1] = 10
-    col, row = centroid_quadratic(data)
+    if data_dtype is float:  # NaN only applicable to float
+        data = np.zeros((5, 5), dtype=data_dtype)
+        data[0, 0] = np.nan
+        data[-1, -1] = 10
+        col, row = centroid_quadratic(data, mask=mask)
+        assert np.isfinite(col) & np.isfinite(row)
+
+    # has NaN in the identified 3x3 patch
+    if data_dtype is float:  # NaN only applicable to float
+        data = np.zeros((5, 5), dtype=data_dtype)
+        data[3, 2] = 10
+        data[3, 3] = np.nan
+        col, row = centroid_quadratic(data, mask=mask)
+        assert np.isfinite(col) & np.isfinite(row)
+
+    # Data contains are all negative (#1401 case `mask` is specified)
+    # e.g., in a difference image
+    data = np.full((5, 5), -9, dtype=data_dtype)
+    data[3, 2] = -5
+    col, row = centroid_quadratic(data, mask=mask)
     assert np.isfinite(col) & np.isfinite(row)
 
+    # case the 3x3 patch contains masked pixels
+    if mask is not None:
+        data = np.zeros((5, 5), dtype=data_dtype)
+        data[2, 1] = 10  # the row above is masked
+        col, row = centroid_quadratic(data, mask=mask)
+        assert np.isfinite(col) & np.isfinite(row)
+
 
 def test_show_citation_instructions():
     show_citation_instructions()