Rewrite histogram operations to remove GPL license on Perreault's code

skimage/filter/_ctmf.pyx

@@ -17,12 +17,10 @@ cimport cython
 from libc.stdlib cimport malloc, free
 from libc.string cimport memset
 
-cdef extern from "emmintrin.h":
-    ctypedef long __m128i
-    __m128i _mm_load_si128 (__m128i *p)
-    __m128i _mm_adds_epu16 (__m128i a, __m128i b)
-    __m128i _mm_subs_epu16 (__m128i a, __m128i b)
-    void _mm_store_si128 (__m128i *p, __m128i a)
+cdef extern from "_histogram.h":
+    ctypedef unsigned short int uint16_t
+    void add16(uint16_t *dest, uint16_t *src)
+    void sub16(uint16_t *dest, uint16_t *src) 
 
 np.import_array()
 
@@ -341,38 +339,6 @@ cdef inline np.int32_t leading_edge_colidx(Histograms *ph, np.int32_t colidx):
 
 cdef inline np.int32_t trailing_edge_colidx(Histograms *ph, np.int32_t colidx):
     return (colidx + 3*ph.radius - 1) % ph.stripe_length
-#
-# add16 - add 16 consecutive integers
-#
-# Add an array of 16 16-bit integers to an accumulator of 16 16-bit integers
-#
-# TO_DO - optimize using SIMD instructions
-#
-cdef inline void add16(np.uint16_t *dest, np.uint16_t *src):
-    cdef __m128i d, s, *pd, *ps
-    pd = <__m128i *> dest
-    ps = <__m128i *> src
-    d = _mm_load_si128(pd)
-    s = _mm_load_si128(ps)
-    d = _mm_adds_epu16 (d, s)
-    _mm_store_si128(pd, d)
-    d = _mm_load_si128(pd + 1)
-    s = _mm_load_si128(ps + 1)
-    d = _mm_adds_epu16 (d, s)
-    _mm_store_si128(pd + 1, d)
-
-cdef inline void sub16(np.uint16_t *dest, np.uint16_t *src):
-    cdef __m128i d, s, *pd, *ps
-    pd = <__m128i *> dest
-    ps = <__m128i *> src
-    d = _mm_load_si128(pd)
-    s = _mm_load_si128(ps)
-    d = _mm_subs_epu16 (d, s)
-    _mm_store_si128(pd, d)
-    d = _mm_load_si128(pd + 1)
-    s = _mm_load_si128(ps + 1)
-    d = _mm_subs_epu16 (d, s)
-    _mm_store_si128(pd + 1, d)
 
 ############################################################################
 #

skimage/filter/_histogram.h

@@ -0,0 +1,119 @@
+/*
+ * Copyright (C) 2013 Thouis "Ray" Jones
+ * Modified from code
+ * Copyright (C) 2006  Simon Perreault
+ *
+ * Reference: S. Perreault and P. Hébert, "Median Filtering in Constant Time",
+ * IEEE Transactions on Image Processing, September 2007.
+ *
+ */
+
+/* Intrinsic declarations */
+#if defined(__SSE2__)
+#include <emmintrin.h>
+#elif defined(__MMX__)
+#include <mmintrin.h>
+#elif defined(__ALTIVEC__)
+#include <altivec.h>
+#endif
+
+/* Compiler peculiarities */
+#if defined(__GNUC__)
+#include <stdint.h>
+#elif defined(_MSC_VER)
+#define inline __inline
+#endif
+
+/**
+ * Add 16 unsigned 16-bit integers using SSE2, MMX or Altivec, if
+ * available.
+ */
+#if defined(__SSE2__)
+static inline void add16(uint16_t *dest, uint16_t *src)
+{
+    __m128i *d, *s;
+    d = (__m128i *) dest;
+    s = (__m128i *) src;
+    *d = _mm_add_epi16(*d, *s);
+    d++; s++;
+    *d = _mm_add_epi16(*d, *s);
+}
+#elif defined(__MMX__)
+static inline void add16(uint16_t *dest, uint16_t *src)
+{
+    __m64 *d, *s;
+    d = (__m64 *) dest;
+    s = (__m64 *) src;
+    *d = _mm_add_pi16(*d, *s);
+    d++; s++;
+    *d = _mm_add_pi16(*d, *s);
+    d++; s++;
+    *d = _mm_add_pi16(*d, *s);
+    d++; s++;
+    *d = _mm_add_pi16(*d, *s);
+}
+#elif defined(__ALTIVEC__)
+static inline void add16(uint16_t *dest, uint16_t *src)
+{
+    vector unsigned short *d, *s;
+    d = (vector unsigned short *) dest;
+    s = (vector unsigned short *) src;
+    *d = vec_add(*d, *s);
+    d++; s++;
+    *d = vec_add(*d, *s);
+}
+#else
+static inline void add16(uint16_t *dest, uint16_t *src)
+{
+    int i;
+
+    for (i = 0; i < 16; i++) dest[i] += src[i];
+}
+#endif
+
+/**
+ * Subtract 16 unsigned 16-bit integers using SSE2, MMX or Altivec, if
+ * available.
+ */
+#if defined(__SSE2__)
+static inline void sub16(uint16_t *dest, uint16_t *src)
+{
+    __m128i *d, *s;
+    d = (__m128i *) dest;
+    s = (__m128i *) src;
+    *d = _mm_sub_epi16(*d, *s);
+    d++; s++;
+    *d = _mm_sub_epi16(*d, *s);
+}
+#elif defined(__MMX__)
+static inline void sub16(uint16_t *dest, uint16_t *src)
+{
+    __m64 *d, *s;
+    d = (__m64 *) dest;
+    s = (__m64 *) src;
+    *d = _mm_sub_pi16(*d, *s);
+    d++; s++;
+    *d = _mm_sub_pi16(*d, *s);
+    d++; s++;
+    *d = _mm_sub_pi16(*d, *s);
+    d++; s++;
+    *d = _mm_sub_pi16(*d, *s);
+}
+#elif defined(__ALTIVEC__)
+static inline void sub16(uint16_t *dest, uint16_t *src)
+{
+    vector unsigned short *d, *s;
+    d = (vector unsigned short *) dest;
+    s = (vector unsigned short *) src;
+    *d = vec_sub(*d, *s);
+    d++; s++;
+    *d = vec_sub(*d, *s);
+}
+#else
+static inline void sub16(uint16_t *dest, uint16_t *src)
+{
+    int i;
+
+    for (i = 0; i < 16; i++) dest[i] -= src[i];
+}
+#endif

skimage/filter/setup.py

@@ -25,7 +25,7 @@ def configuration(parent_package='', top_path=None):
     cython(['rank/percentile_rank.pyx'], working_path=base_path)
     cython(['rank/bilateral_rank.pyx'], working_path=base_path)
 
-    config.add_extension('_ctmf', sources=['_ctmf.c'],
+    config.add_extension('_ctmf', sources=['_ctmf.c', '_histogram.h'],
         include_dirs=[get_numpy_include_dirs()])
     config.add_extension('_denoise_cy', sources=['_denoise_cy.c'],
         include_dirs=[get_numpy_include_dirs(), '../_shared'])