From 32c0b1ca979fcfd201a435f13851666c01247941 Mon Sep 17 00:00:00 2001 From: "Thouis (Ray) Jones" Date: Tue, 12 Feb 2013 15:40:43 -0500 Subject: [PATCH] Rewrite histogram operations to remove GPL license on Perreault's code --- skimage/filter/_ctmf.pyx | 42 ++----------- skimage/filter/_histogram.h | 119 ++++++++++++++++++++++++++++++++++++ skimage/filter/setup.py | 2 +- 3 files changed, 124 insertions(+), 39 deletions(-) create mode 100644 skimage/filter/_histogram.h diff --git a/skimage/filter/_ctmf.pyx b/skimage/filter/_ctmf.pyx index d24a9693..8133efb0 100644 --- a/skimage/filter/_ctmf.pyx +++ b/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) ############################################################################ # diff --git a/skimage/filter/_histogram.h b/skimage/filter/_histogram.h new file mode 100644 index 00000000..968b664f --- /dev/null +++ b/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 +#elif defined(__MMX__) +#include +#elif defined(__ALTIVEC__) +#include +#endif + +/* Compiler peculiarities */ +#if defined(__GNUC__) +#include +#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 diff --git a/skimage/filter/setup.py b/skimage/filter/setup.py index 650a26a6..cfcb0794 100644 --- a/skimage/filter/setup.py +++ b/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'])