wassname/scikit-image
1
0
Fork 0
mirror of https://github.com/wassname/scikit-image.git synced 2026-07-21 12:50:27 +08:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #432 from thouis/sse_ctmf

Browse Source 
Use SSE/MMX/Altivec for faster CTMF
This commit is contained in:
Johannes Schönberger
2013-02-22 07:50:06 -08:00
parent c25ec7f036 d2b98a93e9
commit 9f14edaa57
3 changed files with 172 additions and 74 deletions
Download Patch File Download Diff File Expand all files Collapse all files
CONTRIBUTORS.txt
+3
View File
@@ -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.
skimage/_shared/vectorized_ops.h
+109
View File
@@ -0,0 +1,109 @@
/* 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/_ctmf.pyx
+60 -74
View File
@@ -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 = <Histograms *>ptr
ph = <Histograms *> ptr
if not ptr:
return ph
ph.memory = ptr
ptr = <void *>(ph+1)
ph.pixel_count = <PixelCount *>ptr
ptr = <void *>(ph.pixel_count + adjusted_stripe_length)
ptr = <void *> (ph + 1)
ph.pixel_count = <PixelCount *> ptr
ptr = <void *> (ph.pixel_count + adjusted_stripe_length)
#
# Align histogram memory to a 32-byte boundary
#
roundoff = <size_t>ptr
roundoff = <size_t> ptr
roundoff += 31
roundoff -= roundoff % 32
ptr = <void *>roundoff
ph.histogram = <Histogram *>ptr
ptr = <void *> roundoff
ph.histogram = <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 = <int>(<np.float64_t>radius * 2.0 / 2.414213)
a = <int> (<np.float64_t> 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 <np.int32_t> (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 <np.int32_t> (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 <np.int32_t> (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 <np.int32_t> (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 = <np.uint32_t> ((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 <np.uint8_t> 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(<np.int32_t> data.shape[0], <np.int32_t> data.shape[1],
<np.int32_t> data.strides[0], <np.int32_t> data.strides[1],
radius, percent,
<np.uint8_t *>data.data,
<np.uint8_t *>mask.data,
<np.uint8_t *>output.data):
<np.uint8_t *> data.data,
<np.uint8_t *> mask.data,
<np.uint8_t *> output.data):
raise MemoryError('Failed to allocate scratchpad memory')