diff --git a/CONTRIBUTORS.txt b/CONTRIBUTORS.txt index d582ab63..38e2ea32 100644 --- a/CONTRIBUTORS.txt +++ b/CONTRIBUTORS.txt @@ -132,3 +132,6 @@ - François Boulogne Andres Method for circle perimeter, ellipse perimeter drawing. + +- Thouis Jones + Vectorized operators for arrays of 16-bit ints. diff --git a/skimage/_shared/vectorized_ops.h b/skimage/_shared/vectorized_ops.h new file mode 100644 index 00000000..e426c86e --- /dev/null +++ b/skimage/_shared/vectorized_ops.h @@ -0,0 +1,109 @@ +/* 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/_ctmf.pyx b/skimage/filter/_ctmf.pyx index c133d8d6..5a81cb98 100644 --- a/skimage/filter/_ctmf.pyx +++ b/skimage/filter/_ctmf.pyx @@ -17,6 +17,10 @@ cimport cython from libc.stdlib cimport malloc, free from libc.string cimport memset +cdef extern from "../_shared/vectorized_ops.h": + void add16(np.uint16_t *dest, np.uint16_t *src) + void sub16(np.uint16_t *dest, np.uint16_t *src) + np.import_array() ############################################################################## @@ -58,11 +62,11 @@ cdef struct HistogramPiece: np.uint16_t fine[256] cdef struct Histogram: - HistogramPiece top_left # top-left corner - HistogramPiece top_right # top-right corner - HistogramPiece edge # leading/trailing edge - HistogramPiece bottom_left # bottom-left corner - HistogramPiece bottom_right # bottom-right corner + HistogramPiece top_left # top-left corner + HistogramPiece top_right # top-right corner + HistogramPiece edge # leading/trailing edge + HistogramPiece bottom_left # bottom-left corner + HistogramPiece bottom_right # bottom-right corner # The pixel count has the number of pixels histogrammed in # each of the five compartments for this position. This changes @@ -123,16 +127,16 @@ cdef struct Histograms: # # x --> # - SCoord last_top_left # (-) left side of octagon's top - 1 row - SCoord top_left # (+) -1 row from trailing edge top - SCoord last_top_right # (-) right side of octagon's top - 1 col - 1 row - SCoord top_right # (+) -1 col -1 row from leading edge top - SCoord last_leading_edge # (-) leading edge (right) top stride - 1 row - SCoord leading_edge # (+) leading edge bottom stride - SCoord last_bottom_right # (-) leading edge bottom - 1 col - SCoord bottom_right # (+) right side of octagon's bottom - 1 col - SCoord last_bottom_left # (-) trailing edge bottom - 1 col - SCoord bottom_left # (+) left side of octagon's bottom - 1 col + SCoord last_top_left # (-) left side of octagon's top - 1 row + SCoord top_left # (+) -1 row from trailing edge top + SCoord last_top_right # (-) right side of octagon's top - 1 col - 1 row + SCoord top_right # (+) -1 col -1 row from leading edge top + SCoord last_leading_edge # (-) leading edge (right) top stride - 1 row + SCoord leading_edge # (+) leading edge bottom stride + SCoord last_bottom_right # (-) leading edge bottom - 1 col + SCoord bottom_right # (+) right side of octagon's bottom - 1 col + SCoord last_bottom_left # (-) trailing edge bottom - 1 col + SCoord bottom_left # (+) left side of octagon's bottom - 1 col np.int32_t row_stride # stride between one row and the next np.int32_t col_stride # stride between one column and the next @@ -181,25 +185,25 @@ cdef Histograms *allocate_histograms(np.int32_t rows, SCoord *psc memory_size = (adjusted_stripe_length * - (sizeof(Histogram) + sizeof(PixelCount))+ - sizeof(Histograms)+32) + (sizeof(Histogram) + sizeof(PixelCount)) + + sizeof(Histograms) + 32) ptr = malloc(memory_size) memset(ptr, 0, memory_size) - ph = ptr + ph = ptr if not ptr: return ph ph.memory = ptr - ptr = (ph+1) - ph.pixel_count = ptr - ptr = (ph.pixel_count + adjusted_stripe_length) + ptr = (ph + 1) + ph.pixel_count = ptr + ptr = (ph.pixel_count + adjusted_stripe_length) # # Align histogram memory to a 32-byte boundary # - roundoff = ptr + roundoff = ptr roundoff += 31 roundoff -= roundoff % 32 - ptr = roundoff - ph.histogram = ptr + ptr = roundoff + ph.histogram = ptr # # Fill in the statistical things we keep around # @@ -228,7 +232,7 @@ cdef Histograms *allocate_histograms(np.int32_t rows, # a_2 is the offset from the center to each of the octagon # corners # - a = (radius * 2.0 / 2.414213) + a = ( radius * 2.0 / 2.414213) a_2 = a / 2 if a_2 == 0: a_2 = 1 @@ -323,33 +327,18 @@ cdef void set_stride(Histograms *ph, SCoord *psc): # ############################################################################ cdef inline np.int32_t tl_br_colidx(Histograms *ph, np.int32_t colidx): - return (colidx + 3*ph.radius + ph.current_row) % ph.stripe_length + return (colidx + 3 * ph.radius + ph.current_row) % \ + ph.stripe_length cdef inline np.int32_t tr_bl_colidx(Histograms *ph, np.int32_t colidx): - return (colidx + 3*ph.radius + ph.row_count-ph.current_row) % \ - ph.stripe_length + return (colidx + 3 * ph.radius + ph.row_count - ph.current_row) % \ + ph.stripe_length cdef inline np.int32_t leading_edge_colidx(Histograms *ph, np.int32_t colidx): - return (colidx + 5*ph.radius) % ph.stripe_length + return (colidx + 5 * ph.radius) % ph.stripe_length 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 int i - for i in range(16): - dest[i] += src[i] - -cdef inline void sub16(np.uint16_t *dest, np.uint16_t *src): - cdef int i - for i in range(16): - dest[i] -= src[i] + return (colidx + 3 * ph.radius - 1) % ph.stripe_length ############################################################################ # @@ -540,9 +529,9 @@ cdef inline void update_histogram(Histograms *ph, y = last_coord.y + current_row stride = current_stride+last_coord.stride - if (x >= 0 and x < column_count and - y >= 0 and y < row_count and - ph.mask[stride]): + if (x >= 0 and x < column_count and \ + y >= 0 and y < row_count and \ + ph.mask[stride]): value = ph.data[stride] pixel_count[0] -= 1 hist_piece.fine[value] -= 1 @@ -552,9 +541,9 @@ cdef inline void update_histogram(Histograms *ph, y = coord.y + current_row stride = current_stride + coord.stride - if (x >= 0 and x < column_count and - y >= 0 and y < row_count and - ph.mask[stride]): + if (x >= 0 and x < column_count and \ + y >= 0 and y < row_count and \ + ph.mask[stride]): value = ph.data[stride] pixel_count[0] += 1 hist_piece.fine[value] += 1 @@ -624,8 +613,8 @@ cdef inline np.uint8_t find_median(Histograms *ph): if ph.accumulator_count == 0: return 0 - pixels_below = ((ph.accumulator_count * ph.percent + 50) - / 100) # +50 for roundoff + pixels_below = ((ph.accumulator_count * ph.percent + 50) + / 100) # +50 for roundoff if pixels_below > 0: pixels_below -= 1 @@ -637,7 +626,7 @@ cdef inline np.uint8_t find_median(Histograms *ph): accumulator -= ph.accumulator.coarse[i] update_fine(ph, i) - for j in range(i*16,(i+1)*16): + for j in range(i * 16, (i + 1) * 16): accumulator += ph.accumulator.fine[j] if accumulator > pixels_below: return j @@ -671,9 +660,7 @@ cdef int c_median_filter(np.int32_t rows, cdef: Histograms *ph Histogram *phistogram - int row - int col - int i + int row, col, i np.int32_t top_left_off np.int32_t top_right_off np.int32_t bottom_left_off @@ -682,7 +669,7 @@ cdef int c_median_filter(np.int32_t rows, ph = allocate_histograms(rows, columns, row_stride, col_stride, radius, percent, data, mask, output) if not ph: - return 1 + return 1 for row in range(-radius, rows): # @@ -721,7 +708,7 @@ cdef int c_median_filter(np.int32_t rows, # Update locations and coarse accumulator for the octagon # for points before 0 # - for col in range(-radius, 0 if row >=0 else columns+radius): + for col in range(-radius, 0 if row >= 0 else columns+radius): ph.current_column = col ph.current_stride = row * row_stride + col * col_stride update_current_location(ph) @@ -742,16 +729,15 @@ cdef int c_median_filter(np.int32_t rows, ph.current_stride = row * row_stride + col * col_stride update_current_location(ph) - free_histograms(ph) return 0 -def median_filter( - np.ndarray[dtype=np.uint8_t, ndim=2, negative_indices=False, mode='c'] data, - np.ndarray[dtype=np.uint8_t, ndim=2, negative_indices=False, mode='c'] mask, - np.ndarray[dtype=np.uint8_t, ndim=2, negative_indices=False, mode='c'] output, - int radius, - np.int32_t percent): + +def median_filter(np.ndarray[dtype=np.uint8_t, ndim=2, negative_indices=False, mode='c'] data, + np.ndarray[dtype=np.uint8_t, ndim=2, negative_indices=False, mode='c'] mask, + np.ndarray[dtype=np.uint8_t, ndim=2, negative_indices=False, mode='c'] output, + int radius, + np.int32_t percent): """Median filter with octagon shape and masking. Parameters @@ -773,10 +759,10 @@ def median_filter( """ if percent < 0: - raise ValueError('Median filter percent = %d is less than zero' % \ + raise ValueError('Median filter percent = %d is less than zero' % percent) if percent > 100: - raise ValueError('Median filter percent = %d is greater than 100' % \ + raise ValueError('Median filter percent = %d is greater than 100' % percent) if data.shape[0] != mask.shape[0] or data.shape[1] != mask.shape[1]: raise ValueError('Data shape (%d, %d) is not mask shape (%d, %d)' % @@ -786,10 +772,10 @@ def median_filter( raise ValueError('Data shape (%d, %d) is not output shape (%d, %d)' % (data.shape[0], data.shape[1], output.shape[0], output.shape[1])) - if c_median_filter(data.shape[0], data.shape[1], - data.strides[0], data.strides[1], + if c_median_filter( data.shape[0], data.shape[1], + data.strides[0], data.strides[1], radius, percent, - data.data, - mask.data, - output.data): + data.data, + mask.data, + output.data): raise MemoryError('Failed to allocate scratchpad memory')