fix matched filter image edge bug

.continuous-integration/appveyor/install-miniconda.ps1

@@ -7,7 +7,7 @@ $MINICONDA_URL = "http://repo.continuum.io/miniconda/"
 
 function DownloadMiniconda ($version, $platform_suffix) {
     $webclient = New-Object System.Net.WebClient
-    $filename = "Miniconda-" + $version + "-Windows-" + $platform_suffix + ".exe"
+    $filename = "Miniconda3-" + $version + "-Windows-" + $platform_suffix + ".exe"
 
     $url = $MINICONDA_URL + $filename
 

CHANGES.md

@@ -1,5 +1,12 @@
 Python module versions
 
+v1.0.1 (10 July 2017)
+=====================
+
+* Fix bug when using masked filter and noise array where objects with member
+  pixels at end of image (maximum y coordinate) were erroneously missed.
+
+
 v1.0.0 (30 September 2016)
 ==========================
 

Makefile

@@ -1,9 +1,9 @@
 
 OS ?= $(shell sh -c 'uname -s | tr "[A-Z]" "[a-z]"')
 
-SOMAJOR = 0
-SOMINOR = 6
-SOBUGFIX = 0
+SOMAJOR = 1
+SOMINOR = 0
+SOBUGFIX = 1
 
 ifeq ($(OS), darwin)
 SONAME = libsep.dylib

sep.pyx

@@ -15,7 +15,7 @@ from cpython.version cimport PY_MAJOR_VERSION
 
 np.import_array()  # To access the numpy C-API.
 
-__version__ = "1.0.0"
+__version__ = "1.0.1"
 
 # -----------------------------------------------------------------------------
 # Definitions from the SEP C library

src/extract.c

@@ -392,6 +392,13 @@ int sep_extract(sep_image *image, float thresh, int thresh_type,
 	    {
 	      free(cdscan);
 	      cdscan = NULL;
+              if (filter_type == SEP_FILTER_MATCHED)
+              {
+                  for (xl=0; xl<stacksize; xl++)
+                  {
+                      sigscan[xl] = -BIG;
+                  }
+              }
 	    }
 	  cdscan = dummyscan;
 	}
@@ -445,7 +452,7 @@ int sep_extract(sep_image *image, float thresh, int thresh_type,
 	  curpixinfo.flag = trunflag;
 
           /* set pixel variance/noise based on noise array */
-          if (isvarnoise) {
+          if (isvarthresh) {
             if (xl == w || yl == h) {
               pixsig = pixvar = 0.0;
             }
@@ -464,7 +471,7 @@ int sep_extract(sep_image *image, float thresh, int thresh_type,
 
             /* set `thresh` (This is needed later, even
              * if filter_type is SEP_FILTER_MATCHED */
-            if (isvarthresh) thresh = relthresh * pixsig;
+            thresh = relthresh * pixsig;
           }
 
           /* luflag: is pixel above thresh (Y/N)? */

test.py

@@ -44,6 +44,16 @@
                   ('flags', np.int64)]
 SUPPORTED_IMAGE_DTYPES = [np.float64, np.float32, np.int32]
 
+# If we have a FITS reader, read in the necessary test images
+if not NO_FITS:
+    image_data = getdata(IMAGE_FNAME)
+    image_refback = getdata(BACKIMAGE_FNAME)
+    image_refrms = getdata(RMSIMAGE_FNAME)
+
+
+# -----------------------------------------------------------------------------
+# Helpers
+
 def assert_allclose_structured(x, y):
     """Assert that two structured arrays are close.
 
@@ -56,11 +66,62 @@ def assert_allclose_structured(x, y):
         else:
             assert_equal(x[name], y[name])
 
-# If we have a FITS reader, read in the necessary test images
-if not NO_FITS:
-    image_data = getdata(IMAGE_FNAME)
-    image_refback = getdata(BACKIMAGE_FNAME)
-    image_refrms = getdata(RMSIMAGE_FNAME)
+
+def matched_filter_snr(data, noise, kernel):
+    """Super slow implementation of matched filter SNR for testing.
+
+    At each pixel, value is
+
+        sum(data[i] * kernel[i] / noise[i]^2)
+        -------------------------------------
+        sqrt(sum(kernel[i]^2 / noise[i]^2))
+    """
+    ctr = kernel.shape[0] // 2, kernel.shape[1] // 2
+    kslice = ((0 - ctr[0], kernel.shape[0] - ctr[0]),  # range in axis 0
+              (0 - ctr[1], kernel.shape[1] - ctr[1]))  # range in axis 1
+    out = np.empty_like(data)
+
+    for y in range(data.shape[0]):
+        jmin = y + kslice[0][0] # min and max indicies to sum over
+        jmax = y + kslice[0][1]
+        kjmin = 0               # min and max kernel indicies to sum over
+        kjmax = kernel.shape[0]
+
+        # if we're over the edge of the image, limit extent
+        if jmin < 0:
+            offset = -jmin
+            jmin += offset
+            kjmin += offset
+        if jmax > data.shape[0]:
+            offset = data.shape[0] - jmax
+            jmax += offset
+            kjmax += offset
+
+        for x in range(data.shape[1]):
+            imin = x + kslice[1][0] # min and max indicies to sum over
+            imax = x + kslice[1][1]
+            kimin = 0               # min and max kernel indicies to sum over
+            kimax = kernel.shape[1]
+
+            # if we're over the edge of the image, limit extent
+            if imin < 0:
+                offset = -imin
+                imin += offset
+                kimin += offset
+            if imax > data.shape[1]:
+                offset = data.shape[1] - imax
+                imax += offset
+                kimax += offset
+
+            d = data[jmin:jmax, imin:imax]
+            n = noise[jmin:jmax, imin:imax]
+            w = 1. / n**2
+            k = kernel[kjmin:kjmax, kimin:kimax]
+            out[y, x] = np.sum(d * k * w) / np.sqrt(np.sum(k**2 * w))
+
+    return out
+
+
 # -----------------------------------------------------------------------------
 # Test versus Source Extractor results
 
@@ -331,6 +392,24 @@ def test_extract_with_noise_convolution():
     assert_approx_equal(objects[1]['y'], 3.)
 
 
+def test_extract_matched_filter_at_edge():
+    """Exercise bug where bright star at end of image not detected
+    with noise array and matched filter on."""
+
+    data = np.zeros((20, 20))
+    err = np.ones_like(data)
+    kernel = np.array([[1., 2., 1.],
+                       [2., 4., 2.],
+                       [1., 2., 1.]])
+
+    data[18:20, 9:12] = kernel[0:2, :]
+
+    objects, pix = sep.extract(data, 2.0, err=err, filter_kernel=kernel,
+                               filter_type="matched", segmentation_map=True)
+    assert len(objects) == 1
+    assert objects["npix"][0] == 6
+
+
 @pytest.mark.skipif(NO_FITS, reason="no FITS reader")
 def test_extract_with_mask():