Cython wrapper to compute sliding window normalized xcorr

.gitignore

@@ -29,3 +29,4 @@ station.dat
 .idea
 .coverage.*
 eqcorrscan/core/subspace_statistic.c
+eqcorrscan/core/sliding_normxcorr.c

eqcorrscan/core/sliding_normxcorr.pyx

@@ -0,0 +1,114 @@
+"""
+Function for computing a sliding window normalised cross-correlation
+function using openCV and numpy.  Written because obspy's xcorr gives different
+results on different systems.  Designed only for small windows.
+
+:copyright:
+    Calum Chamberlain, Chet Hopp.
+
+:license:
+    GNU Lesser General Public License, Version 3
+    (https://www.gnu.org/copyleft/lesser.html)
+"""
+
+from __future__ import absolute_import
+from __future__ import division
+from __future__ import print_function
+from __future__ import unicode_literals
+import numpy as np
+cimport numpy as np
+import cython
+from cv2 import matchTemplate, TM_CCOEFF_NORMED
+
+DTYPE = np.float32
+ctypedef np.float32_t DTYPE_t
+
+@cython.boundscheck(False)
+@cython.wraparound(False)
+def sliding_normxcorr(np.ndarray[DTYPE_t, ndim=1] arr1,
+                      np.ndarray[DTYPE_t, ndim=1] arr2,
+                      int shift_len,
+                      int full_xcorr=False):
+    """
+    Calculate the sliding normalized cross-correlation of a pair of arrays.
+
+    Uses the openCV match_template routine for correlation, and implements a \
+    cythonised sliding-window approach. shift_len must be less than half the \
+    shortest array length to provide reasonable support. Shifts +/- the \
+    shift_len
+
+    :type arr1: np.ndarray
+    :param arr1: First array, shifts will be relative to this array.  Must \
+        be np.float32 datatype
+    :type arr2: np.ndarry
+    :param arr2: Second array to slide through arr 1.  Must be np.float32 \
+        datatype.
+    :type shift_len: int
+    :param shift_len: Length in samples to compute the sliding window for.
+    :type full_xcorr: int
+    :param full_xcorr: Whether to return the full cross-correlation or not. \
+        If False, only the shift for the maximum correlation and the \
+        correlation value will be returned, if True, will also return the \
+        correlation vector. 0=False, 1=True
+
+    :rtype: float
+    :return: Shift (in samples) for maximum correlation
+    :type: float
+    :return: Maximum correlation.
+    :rtype: np.ndarray
+    :return: Full cross-correlation array if full_xcorr=True.
+
+    .. Note:: This is not fast and is a bit silly - should be written in pure \
+        c++ and wrapped.
+    """
+    cdef int i
+    cdef int shift
+    cdef DTYPE_t corr
+    cdef int arr1_len = arr1.shape[0]
+    cdef int arr2_len = arr2.shape[0]
+    cdef int min_len
+    cdef int arr1_start
+    cdef int arr1_end
+    cdef int arr2_start
+    cdef int arr2_end
+    cdef np.ndarray[DTYPE_t, ndim=1] xcorr = np.empty((shift_len * 2) + 1,
+                                                      dtype=DTYPE)
+    if arr1_len <= arr2_len:
+        min_len = arr1_len
+    else:
+        min_len = arr2_len
+    # De-mean
+    arr1 -= arr1.mean()
+    arr2 -= arr2.mean()
+    # Normalize
+    arr1 /= arr1.max()
+    arr2 /= arr2.max()
+    # Check that shift_len is less than half the shortest array length
+    if shift_len * 2 > min_len:
+        raise IndexError('Shift_len is too large, not enough support.')
+    # Do the first half, reduces compute a little
+    arr1_start = 0
+    arr2_end = arr2_len
+    for i in range(shift_len):
+        arr1_end = arr1_len - (shift_len - i)
+        arr2_start = arr2_end - (arr1_end - arr1_start)
+        xcorr[i] = matchTemplate(arr1[arr1_start: arr1_end],
+                                 arr2[arr2_start: arr2_end],
+                                 method=TM_CCOEFF_NORMED)[0][0]
+    arr1_end = arr1_len
+    arr2_start = 0
+    for i in range(shift_len + 1):
+        arr2_end = arr2_len - i
+        arr1_start = arr1_end - (arr2_end - arr2_start)
+        xcorr[i + shift_len] = matchTemplate(arr1[arr1_start: arr1_end],
+                                             arr2[arr2_start: arr2_end],
+                                             method=TM_CCOEFF_NORMED)[0][0]
+
+    shift = xcorr.argmax() - shift_len
+    corr = xcorr.max()
+    if full_xcorr == 1:
+        return shift, corr, xcorr
+    else:
+        return shift, corr
+
+

eqcorrscan/core/sliding_xcorr.py

@@ -0,0 +1,76 @@
+"""
+Function for computing a sliding window normalised cross-correlation
+function using openCV and numpy.  Written because obspy's xcorr gives different
+results on different systems.  Designed only for small windows.
+
+:copyright:
+    Calum Chamberlain, Chet Hopp.
+
+:license:
+    GNU Lesser General Public License, Version 3
+    (https://www.gnu.org/copyleft/lesser.html)
+"""
+
+from __future__ import absolute_import
+from __future__ import division
+from __future__ import print_function
+from __future__ import unicode_literals
+import numpy as np
+from eqcorrscan.core import sliding_xcorrlib
+
+
+def sliding_xcorr(arr1, arr2, shift_len, full_xcorr=False):
+    """
+    Calculate the sliding normalized cross-correlation of a pair of arrays.
+
+    Uses the openCV match_template routine for correlation, and implements a \
+    cythonised sliding-window approach. shift_len must be less than half the \
+    shortest array length to provide reasonable support. Shifts +/- the \
+    shift_len
+
+    :type arr1: np.ndarray
+    :param arr1: First array, shifts will be relative to this array.  Must \
+        be np.float32 datatype
+    :type arr2: np.ndarry
+    :param arr2: Second array to slide through arr 1.  Must be np.float32 \
+        datatype.
+    :type shift_len: int
+    :param shift_len: Length in samples to compute the sliding window for.
+    :type full_xcorr: int
+    :param full_xcorr: Whether to return the full cross-correlation or not. \
+        If False, only the shift for the maximum correlation and the \
+        correlation value will be returned, if True, will also return the \
+        correlation vector. 0=False, 1=True
+
+    :rtype: float
+    :return: Shift (in samples) for maximum correlation
+    :type: float
+    :return: Maximum correlation.
+    :rtype: np.ndarray
+    :return: Full cross-correlation array if full_xcorr=True.
+    """
+    # De-mean
+    arr1 -= arr1.mean()
+    arr2 -= arr2.mean()
+    # Normalize
+    arr1 /= arr1.max()
+    arr2 /= arr2.max()
+    # Make contiguous
+    arr1 = np.ascontiguousarray(arr1, np.float32)
+    arr2 = np.ascontiguousarray(arr2, np.float32)
+    # Check that shift_len is less than half the shortest array length
+    if shift_len * 2 > np.min([len(arr1), len(arr2)]):
+        raise IndexError('Shift_len is too large, not enough support.')
+    xcorr = np.empty((shift_len * 2) + 1, dtype=np.float32, order='C')
+    print(len(arr1))
+    print(len(arr2))
+    res = sliding_xcorrlib.xcorr(arr1, arr2, xcorr, shift_len, len(arr1),
+                                 len(arr2))
+    if res:
+        raise MemoryError('Internal C error %s' % res)
+    shift = xcorr.argmax() - shift_len
+    corr = xcorr.max()
+    if full_xcorr:
+        return shift, corr, xcorr
+    else:
+        return shift, corr

eqcorrscan/core/sliding_xcorrlib.c

@@ -0,0 +1,73 @@
+/*-------------------------------------------------------------------------
+# Filename: sliding_xcorr.c
+#  Purpose: Sliding window, normalized cross-correlation using openCV
+#   Author: Calum John Chamberlain
+#  Licence: LGPL, Verison 3
+#--------------------------------------------------------------------------*/
+#include <Python.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <memory.h>
+#include <float.h>
+
+static PyObject *
+xcorr(PyObject *self, PyObject *args)
+//float *arr1, float *arr2, float *ccs, const int shift_len, const int arr1_len, const int arr2_len)
+{
+    int i;
+    int shift;
+    int min_len;
+    int arr1_start;
+    int arr1_end;
+    int arr2_start;
+    int arr2_end;
+    float *arr1
+    float *arr2
+    float *ccs
+    int shift_len
+    int arr1_len
+    int arr2_len
+
+    // Extract the arguments and convert to useful C things
+    if (!PyArg_ParseTuple(args, "s", &command))
+        return NULL;
+
+
+    printf("Array 1 length: %d \n", arr1_len);
+    printf("Array 2 length: %d \n", arr2_len);
+    if (arr1_len <= arr2_len)
+    {
+        min_len = arr1_len;
+    }
+    else
+    {
+        min_len = arr2_len;
+    }
+    printf("Minimum_length = %d \n", min_len);
+    arr1_start = 0;
+    arr2_end = arr2_len;
+//    for (i=0;i<min_len;i++)
+//    {
+//        arr1_end = arr1_len - (shift_len - i);
+//        arr2_start = arr2_end - (arr1_end - arr1_start);
+//        printf("arr1_start: %d \n", arr1_start);
+//        printf("arr1_end: %d \n", arr1_end);
+//        printf("arr2_start: %d \n", arr2_start);
+//        printf("arr2_end: %d \n", arr2_end);
+//    }
+    return 0;
+}
+
+static PyMethodDef XcorrMethods[] = {
+    {"xcorr",  xcorr, METH_VARARGS,
+     "Compute the sliding window cross correlation"},
+    {NULL, NULL, 0, NULL}        /* Sentinel */
+};
+
+PyMODINIT_FUNC
+initsliding_xcorrlib(void)
+{
+    (void) Py_InitModule("sliding_xcorrlib", XcorrMethods);
+}

eqcorrscan/tests/sliding_normxcorr_test.py

@@ -0,0 +1,48 @@
+"""
+Functions for testing the utils.stacking functions
+"""
+from __future__ import absolute_import
+from __future__ import division
+from __future__ import print_function
+from __future__ import unicode_literals
+from eqcorrscan.core import sliding_normxcorr
+import unittest
+import numpy as np
+
+
+class TestSliding(unittest.TestCase):
+    """Test the sliding window correlation."""
+    def test_sine(self):
+        arr1 = np.sin(np.arange(10000)/200.0)
+        arr2 = np.sin(np.arange(10000)/200.0+10)
+        shift, ccs = sliding_normxcorr.\
+            sliding_normxcorr(arr1.astype(np.float32),
+                              arr2.astype(np.float32),
+                              600)
+        self.assertEqual(round(ccs, 5), 1.0)
+        self.assertEqual(shift, -513)  # Interpolation rounding
+
+    def test_fail(self):
+        arr1 = np.sin(np.arange(10000)/200.0)
+        arr2 = np.sin(np.arange(10000)/200.0+10)
+        with self.assertRaises(IndexError):
+            shift, ccs = sliding_normxcorr.\
+                sliding_normxcorr(arr1.astype(np.float32),
+                                  arr2.astype(np.float32),
+                                  100000)
+
+    def test_full(self):
+        arr1 = np.sin(np.arange(10000)/200.0)
+        arr2 = np.sin(np.arange(10000)/200.0+10)
+        shift, ccs, xcorr = sliding_normxcorr.\
+            sliding_normxcorr(arr1.astype(np.float32),
+                              arr2.astype(np.float32),
+                              10, full_xcorr=True)
+        self.assertEqual(len(xcorr), 21)
+
+
+if __name__ == '__main__':
+    """
+    Run stacking tests
+    """
+    unittest.main()

eqcorrscan/tests/stacking_test.py

@@ -61,7 +61,7 @@ def test_phase_weighted_stack(self):
         self.assertEqual(len(synth[0].data), len(stack[0].data))
 
     def test_align_traces(self):
-        """Test the utils.stacking.align_traces fucntion."""
+        """Test the utils.stacking.align_traces function."""
         # Generate synth data
         import numpy as np
         from obspy import Trace

eqcorrscan/tests/subspace_test.py

@@ -190,7 +190,8 @@ def test_create_multiplexed_aligned(self):
         detector = subspace.Detector()
         detector.construct(streams=templates, lowcut=2, highcut=9,
                            filt_order=4, sampling_rate=20, multiplex=True,
-                           name=str('Tester'), align=True, shift_len=0.2)
+                           name=str('Tester'), align=True, shift_len=3.0,
+                           reject=0.2)
         for u in detector.data:
             identity = np.dot(u.T, u).astype(np.float16)
             self.assertTrue(np.allclose(identity,