Merge remote-tracking branch 'refs/remotes/origin/master'

README.md

@@ -63,10 +63,21 @@ To install the DRP along with its dependencies, you need to run the following st
     cd lvmdrp
     ```
 
+3. Optional, but recommended: compile c-code for better performance
 
-3. Install the DRP package in the current python environment (see [contributing](#contributing-to-lvm-drp-development) section below for a replacement of this step):
+    ```bash
+    pushd lvmdrp/python/cextern/fast_median/src
+    make
+    popd
+    ```
+    This should have created a fast_median.so (Linux) or fast_median.dylib (MacOS) in the src directory.
+    Set environment variable `LVMDRP_LIB_DIR` if you move the shared library. If you leave it in src/, no need to set this.
+    We will automate this step eventually ;-)
+
+4. Install the DRP package in the current python environment (see [contributing](#contributing-to-lvm-drp-development) section below for a replacement of this step). Make sure you're back in the lvmdrp directory.
 
     ```bash
+    cd lvmdrp
     pip install .
     ```
 

python/cextern/README.txt

@@ -0,0 +1,8 @@
+C code to be compiled on setup goes here.
+
+fast_median.cpp is a very fast median filter for 1d and 2d ndarrays
+with float or double data. The glue code is in fast_median.py
+
+compile using the Makefile in src/
+
+

python/cextern/fast_median/fast_median.py

@@ -0,0 +1,86 @@
+import ctypes
+import numpy
+from numpy.ctypeslib import ndpointer
+from scipy.ndimage import median_filter
+
+HAVE_MEDIAN_SO = False
+try:
+    import platform
+    from os import path, environ
+    s = platform.system()
+    resources_dir = environ.get('LVMDRP_LIB_DIR') or path.join(path.dirname(__file__), 'src')
+    if s=='Linux':
+        resources_dir = path.join(resources_dir, 'fast_median.so')
+        #print(resources_dir)
+        lib = ctypes.cdll.LoadLibrary(resources_dir)
+    elif s=='Darwin':
+        resources_dir = path.join(resources_dir, 'fast_median.dylib')
+        #print(resources_dir)
+        lib = ctypes.cdll.LoadLibrary(resources_dir)
+    else:
+        raise Exception('Unknown platform: '+s)
+
+    #template <typename T>
+    #void median_filter_2d(int x, int y, int hx, int hy, int blockhint, const T* in, T* out);
+    lib.median_filter_2d_float.argtypes = (ctypes.c_int, ctypes.c_int,
+                                        ctypes.c_int, ctypes.c_int,
+                                        ctypes.c_int,
+                                        ndpointer(ctypes.c_float, flags="C_CONTIGUOUS"),
+                                        ndpointer(ctypes.c_float, flags="C_CONTIGUOUS"))
+    lib.median_filter_2d_float.restype = None
+    lib.median_filter_2d_double.argtypes = (ctypes.c_int, ctypes.c_int,
+                                        ctypes.c_int, ctypes.c_int,
+                                        ctypes.c_int,
+                                        ndpointer(ctypes.c_double, flags="C_CONTIGUOUS"),
+                                        ndpointer(ctypes.c_double, flags="C_CONTIGUOUS"))
+    lib.median_filter_2d_double.restype = None
+
+    lib.median_filter_1d_float.argtypes = (ctypes.c_int, 
+                                        ctypes.c_int,
+                                        ctypes.c_int,
+                                        ndpointer(ctypes.c_float, flags="C_CONTIGUOUS"),
+                                        ndpointer(ctypes.c_float, flags="C_CONTIGUOUS"))
+    lib.median_filter_1d_float.restype = None
+    lib.median_filter_1d_double.argtypes = (ctypes.c_int,
+                                        ctypes.c_int,
+                                        ctypes.c_int,
+                                        ndpointer(ctypes.c_double, flags="C_CONTIGUOUS"),
+                                        ndpointer(ctypes.c_double, flags="C_CONTIGUOUS"))
+    lib.median_filter_1d_double.restype = None
+
+    HAVE_MEDIAN_SO = True
+except Exception:
+    print(Exception)
+    pass
+
+print('HAVE_MEDIAN_SO ',HAVE_MEDIAN_SO)
+
+if HAVE_MEDIAN_SO is True:
+    def fast_median_filter_2d(input, size=None):
+        assert len(input.shape) == 2, 'fast_median_filter_2d requires 2-dimensional input'
+        assert len(size) == 2, 'fast_median_filter_2d requires two element size tuple, list or array'
+        out = numpy.empty(input.shape, dtype=input.dtype)
+        # Note we have to flip the shape and box x/y to match python's indexing
+        if input.dtype==numpy.float32:
+            lib.median_filter_2d_float(input.shape[1], input.shape[0], size[1], size[0], 0, input, out)
+        elif input.dtype==numpy.float64:
+            lib.median_filter_2d_double(input.shape[1], input.shape[0], size[1], size[0], 0, input, out)
+        else:
+            raise TypeError('median_filter_2d not implemented for type '+str(input.dtype))
+        return out
+
+    def fast_median_filter_1d(input, size=None):
+        assert len(input.shape) == 1, 'fast_median_filter_1d requires 1-dimensional input'
+        out = numpy.empty(input.shape, dtype=input.dtype)
+        if input.dtype==numpy.float32:
+            lib.median_filter_1d_float(input.shape[0], size, 0, input, out)
+        elif input.dtype==numpy.float64:
+            lib.median_filter_1d_double(input.shape[0], size, 0, input, out)
+        else:
+            raise TypeError('median_filter_1d not implemented for type '+str(input.dtype))
+        return out
+else:
+    def fast_median_filter_2d(input, size=None):
+        return median_filter(input, size=size)
+    def fast_median_filter_1d(input, size=None):
+        return median_filter(input, size=size)

python/cextern/fast_median/src/Makefile

@@ -0,0 +1,14 @@
+UNAME_S := $(shell uname -s)
+ifeq ($(UNAME_S),Linux)
+	FLAGS = -fPIC -shared -O3 -march=native
+	TRG = fast_median.so
+	COMPILER = g++
+endif
+ifeq ($(UNAME_S),Darwin)
+        FLAGS = -dynamiclib -current_version 1.0 -compatibility_version 1.0 -O3 -march=native
+	TRG = fast_median.dylib
+	COMPILER = clang++
+endif
+
+all: fast_median.cpp fast_median.hpp
+	$(COMPILER) $(FLAGS) -o $(TRG) fast_median.cpp