wassname/scikit-image
1
0
Fork 0
mirror of https://github.com/wassname/scikit-image.git synced 2026-07-01 14:36:59 +08:00
Code Issues Packages Projects Releases Wiki Activity

Merge constant time median filter by Lee Kamentsky and Thouis Jones.

Browse Source
This commit is contained in:
Stefan van der Walt
2011-03-13 16:31:02 +02:00
parent 9ff324b0fa 386ad23f55
commit ea3b46d80a
8 changed files with 1056 additions and 1 deletions
Download Patch File Download Diff File Expand all files Collapse all files
CONTRIBUTORS.txt
+3
View File
@@ -30,3 +30,6 @@
- Almar Klein
Binary heap class for graph algorithms
- Lee Kamentsky and Thouis Jones of the CellProfiler team, Broad Institute, MIT
Constant time per pixel median filter
scikits/image/filter/__init__.py
+1 -1
View File
@@ -1,2 +1,2 @@
from lpi_filter import *
from ctmf import median_filter
scikits/image/filter/_ctmf.pyx
+795
View File
@@ -0,0 +1,795 @@
'''_ctmf.pyx - constant time per pixel median filtering
Reference: S. Perreault and P. Hebert, "Median Filtering in Constant Time",
IEEE Transactions on Image Processing, September 2007.
Originally part of CellProfiler, code licensed under both GPL and BSD licenses.
Website: http://www.cellprofiler.org
Copyright (c) 2003-2009 Massachusetts Institute of Technology
Copyright (c) 2009-2011 Broad Institute
All rights reserved.
Original author: Lee Kamentsky
'''
import numpy as np
cimport numpy as np
cimport cython
from libc.stdlib cimport malloc, free
from libc.string cimport memset
np.import_array()
##############################################################################
#
# median_filter - implementation of constant-time median filter with
# octagonal shape. The algorithm is derived from
# Perreault, "Median Filtering in Constant Time",
# IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 16, NO. 9,
# SEPTEMBER 2007.
#
# Inputs:
# a 2d array of uint8's to be median filtered
# a similarly shaped uint8 masking array with "1" indicating a significant
# pixel and "0" indicating a pixel to be masked.
#
# Outputs:
# a 2d median-filtered array
#
##############################################################################
DTYPE_UINT32 = np.uint32
DTYPE_BOOL = np.bool
ctypedef np.uint16_t pixel_count_t
###########
#
# Histograms
#
# There are five separate histograms for the octagonal filter and
# there are two levels (coarse = 16 values, fine = 256 values)
# per histogram. There are four histograms to maintain per position
# representing the four diagonals of the histogram plus one histogram
# for the straight side (which is used for adding and subtracting)
#
###########
cdef struct HistogramPiece:
np.uint16_t coarse[16]
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
# The pixel count has the number of pixels histogrammed in
# each of the five compartments for this position. This changes
# because of the mask
#
cdef struct PixelCount:
pixel_count_t top_left
pixel_count_t top_right
pixel_count_t edge
pixel_count_t bottom_left
pixel_count_t bottom_right
#
# Stride + coordinates: the info we need when computing
# relative offsets from the octagon center
#
cdef struct SCoord:
np.int32_t stride # add the stride to the memory location
np.int32_t x
np.int32_t y
cdef struct Histograms:
void *memory # pointer to the allocated memory
Histogram *histogram # pointer to the histogram memory
PixelCount *pixel_count # pointer to the pixel count memory
np.uint8_t *data # pointer to the image data
np.uint8_t *mask # pointer to the image mask
np.uint8_t *output # pointer to the output array
np.int32_t column_count # number of columns represented by this
# structure
np.int32_t stripe_length # number of columns including "radius" before
# and after
np.int32_t row_count # number of rows available in image
np.int32_t current_column # the column being processed
np.int32_t current_row # the row being processed
np.int32_t current_stride # offset in data and mask to current location
np.int32_t radius # the "radius" of the octagon
np.int32_t a_2 # 1/2 of the length of a side of the octagon
#
#
# The strides are the offsets in the array to the points that need to
# be added or removed from a histogram to shift from the previous row
# to the current row.
# Listed here going clockwise from the trailing edge's top.
# (-) = needs to be removed
# (+) = needs to be added
#
# - -
# 1.=========2
# 1. 2
# +. +- Y
# |. 3 |
# |. 3 |
# -. X \|/
# 5. 4 v
# 5. 4
# +.=========+
#
# 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
np.int32_t row_stride # stride between one row and the next
np.int32_t col_stride # stride between one column and the next
# The accumulator holds the running histogram
#
HistogramPiece accumulator
#
# The running count of pixels in the accumulator
#
np.uint32_t accumulator_count
#
# The percent of pixels within the octagon whose value is
# less than or equal to the median-filtered value (e.g. for
# median, this is 50, for lower quartile it's 25)
#
np.int32_t percent
#
# last_update_column keeps track of the column # of the last update
# to the fine histogram accumulator. Short-term, the median
# stays in one coarse block so only one fine histogram might
# need to be updated
#
np.int32_t last_update_column[16]
############################################################################
#
# allocate_histograms - allocates the Histograms structure for the run
#
############################################################################
cdef Histograms *allocate_histograms(np.int32_t rows,
np.int32_t columns,
np.int32_t row_stride,
np.int32_t col_stride,
np.int32_t radius,
np.int32_t percent,
np.uint8_t *data,
np.uint8_t *mask,
np.uint8_t *output):
cdef:
unsigned int adjusted_stripe_length = columns + 2*radius + 1
unsigned int memory_size
void *ptr
Histograms *ph
size_t roundoff
int a
SCoord *psc
memory_size = (adjusted_stripe_length *
(sizeof(Histogram) + sizeof(PixelCount))+
sizeof(Histograms)+32)
ptr = malloc(memory_size)
memset(ptr, 0, memory_size)
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)
#
# Align histogram memory to a 32-byte boundary
#
roundoff = <size_t>ptr
roundoff += 31
roundoff -= roundoff % 32
ptr = <void *>roundoff
ph.histogram = <Histogram *>ptr
#
# Fill in the statistical things we keep around
#
ph.column_count = columns
ph.row_count = rows
ph.current_column = -radius
ph.stripe_length = adjusted_stripe_length
ph.current_row = 0
ph.radius = radius
ph.percent = percent
ph.row_stride = row_stride
ph.col_stride = col_stride
ph.data = data
ph.mask = mask
ph.output = output
#
# Compute the coordinates of the significant points
# (the SCoords)
#
# First, the length of a side of an octagon, compared
# to what we call the radius is:
# 2*r
# ----------- = a
# (1+sqrt(2))
#
# 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_2 = a / 2
if a_2 == 0:
a_2 = 1
ph.a_2 = a_2
if radius <= a_2:
radius = a_2+1
ph.radius = radius
ph.last_top_left.x = -a_2
ph.last_top_left.y = -radius - 1
ph.top_left.x = -radius
ph.top_left.y = -a_2 - 1
ph.last_top_right.x = a_2 - 1
ph.last_top_right.y = -radius - 1
ph.top_right.x = radius - 1
ph.top_right.y = -a_2 - 1
ph.last_leading_edge.x = radius
ph.last_leading_edge.y = -a_2 - 1
ph.leading_edge.x = radius
ph.leading_edge.y = a_2
ph.last_bottom_right.x = radius
ph.last_bottom_right.y = a_2
ph.bottom_right.x = a_2
ph.bottom_right.y = radius
ph.last_bottom_left.x = -radius-1
ph.last_bottom_left.y = a_2
ph.bottom_left.x = -a_2-1
ph.bottom_left.y = radius
#
# Set the stride of each SCoord based on its x and y
#
set_stride(ph, &ph.last_top_left)
set_stride(ph, &ph.top_left)
set_stride(ph, &ph.last_top_right)
set_stride(ph, &ph.top_right)
set_stride(ph, &ph.last_leading_edge)
set_stride(ph, &ph.leading_edge)
set_stride(ph, &ph.last_bottom_left)
set_stride(ph, &ph.bottom_left)
set_stride(ph, &ph.last_bottom_right)
set_stride(ph, &ph.bottom_right)
return ph
############################################################################
#
# free_histograms - frees the Histograms structure
#
############################################################################
cdef void free_histograms(Histograms *ph):
free(ph.memory)
############################################################################
#
# set_stride - set the stride of a SCoord from its X and Y
#
############################################################################
cdef void set_stride(Histograms *ph, SCoord *psc):
psc.stride = psc.x * ph.col_stride + psc.y * ph.row_stride
############################################################################
#
# <tl,tr,bl,br>_colidx - convert a column index into the histogram
# index for a diagonal
#
# The top-right and bottom left diagonals for one row at one column
# become the diagonals for the next column to the right for the next row.
# Conversely, the top-left and bottom right become the diagonals for the
# previous column.
#
# These functions use the current row number to find the index of
# a particular histogram taking this into account. The indices progress
# forward or backward as you go to successive rows.
#
# The histogram array is, in effect, a circular buffer, so the start
# offset is immaterial - we take advantage of this to make sure that
# the numbers computed before taking the modulus are all positive, including
# those that might be done for columns to the left of 0. We add 3* the radius
# here to account for a row of -radius, a column of -radius and a request for
# a column that is "radius" to the left.
#
############################################################################
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
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
cdef inline np.int32_t leading_edge_colidx(Histograms *ph, np.int32_t colidx):
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]
############################################################################
#
# accumulate_coarse_histogram - accumulate the coarse histogram
# at an index into the accumulator
#
# ph - the Histograms structure that holds the accumulator
# colidx - the index of the column to add
#
############################################################################
cdef inline void accumulate_coarse_histogram(Histograms *ph, np.int32_t colidx):
cdef:
int offset
offset = tr_bl_colidx(ph, colidx)
if ph.pixel_count[offset].top_right > 0:
add16(ph.accumulator.coarse, ph.histogram[offset].top_right.coarse)
ph.accumulator_count += ph.pixel_count[offset].top_right
offset = leading_edge_colidx(ph, colidx)
if ph.pixel_count[offset].edge > 0:
add16(ph.accumulator.coarse, ph.histogram[offset].edge.coarse)
ph.accumulator_count += ph.pixel_count[offset].edge
offset = tl_br_colidx(ph, colidx)
if ph.pixel_count[offset].bottom_right > 0:
add16(ph.accumulator.coarse, ph.histogram[offset].bottom_right.coarse)
ph.accumulator_count += ph.pixel_count[offset].bottom_right
############################################################################
#
# deaccumulate_coarse_histogram - subtract the coarse histogram
# for a given column
#
############################################################################
cdef inline void deaccumulate_coarse_histogram(Histograms *ph, np.int32_t colidx):
cdef:
int offset
#
# The trailing diagonals don't appear until here
#
if colidx <= ph.a_2:
return
offset = tl_br_colidx(ph, colidx)
if ph.pixel_count[offset].top_left > 0:
sub16(ph.accumulator.coarse, ph.histogram[offset].top_left.coarse)
ph.accumulator_count -= ph.pixel_count[offset].top_left
#
# The trailing edge doesn't appear from the border until here
#
if colidx > ph.radius:
offset = trailing_edge_colidx(ph, colidx)
if ph.pixel_count[offset].edge > 0:
sub16(ph.accumulator.coarse, ph.histogram[offset].edge.coarse)
ph.accumulator_count -= ph.pixel_count[offset].edge
offset = tr_bl_colidx(ph, colidx)
if ph.pixel_count[offset].bottom_left > 0:
sub16(ph.accumulator.coarse, ph.histogram[offset].bottom_left.coarse)
ph.accumulator_count -= ph.pixel_count[offset].bottom_left
############################################################################
#
# accumulate_fine_histogram - accumulate one of the 16 fine histograms
#
############################################################################
cdef inline void accumulate_fine_histogram(Histograms *ph,
np.int32_t colidx,
np.uint32_t fineidx):
cdef:
int fineoffset = fineidx * 16
int offset
offset = tr_bl_colidx(ph, colidx)
add16(ph.accumulator.fine + fineoffset,
ph.histogram[offset].top_right.fine + fineoffset)
offset = leading_edge_colidx(ph, colidx)
add16(ph.accumulator.fine + fineoffset,
ph.histogram[offset].edge.fine + fineoffset)
offset = tl_br_colidx(ph, colidx)
add16(ph.accumulator.fine + fineoffset,
ph.histogram[offset].bottom_right.fine + fineoffset)
############################################################################
#
# deaccumulate_fine_histogram - subtract one of the 16 fine histograms
#
############################################################################
cdef inline void deaccumulate_fine_histogram(Histograms *ph,
np.int32_t colidx,
np.uint32_t fineidx):
cdef:
int fineoffset = fineidx * 16
int offset
#
# The trailing diagonals don't appear until here
#
if colidx < ph.a_2:
return
offset = tl_br_colidx(ph, colidx)
sub16(ph.accumulator.fine + fineoffset,
ph.histogram[offset].top_left.fine + fineoffset)
if colidx >= ph.radius:
offset = trailing_edge_colidx(ph, colidx)
sub16(ph.accumulator.fine+fineoffset,
ph.histogram[offset].edge.fine + fineoffset)
offset = tr_bl_colidx(ph, colidx)
sub16(ph.accumulator.fine + fineoffset,
ph.histogram[offset].bottom_left.fine + fineoffset)
############################################################################
#
# accumulate - add the leading edge and subtract the trailing edge
#
############################################################################
cdef inline void accumulate(Histograms *ph):
cdef:
int i
int j
np.int32_t accumulator
accumulate_coarse_histogram(ph, ph.current_column)
deaccumulate_coarse_histogram(ph, ph.current_column)
############################################################################
#
# update_fine - update one of the fine histograms to the current column
#
# The code has two choices:
# redo the fine histogram from scratch - this involves accumulating
# the entire histogram from the top_left.x to the top_right.x,
# the center (edge) histogram from the trailing edge x to the
# top_left.x and then computing a histogram of all points between
# the trailing edge top, the point, (top_left.x,trailing edge top.y)
# and the top_right and the corresponding triangle in the octagon's
# lower half.
#
# accumulate and deaccumulate within the fine histogram from the last
# column computed.
#
# The code below only implements the accumulate; redo and the code
# to choose remains to be done.
############################################################################
cdef inline void update_fine(Histograms *ph, int fineidx):
cdef:
int first_update_column = ph.last_update_column[fineidx]+1
int update_limit = ph.current_column+1
int i
for i in range(first_update_column, update_limit):
accumulate_fine_histogram(ph, i, fineidx)
deaccumulate_fine_histogram(ph, i, fineidx)
ph.last_update_column[fineidx] = ph.current_column
############################################################################
#
# update_histogram - update the coarse and fine levels of a histogram
# based on addition of one value and subtraction of another
#
# ph - Histograms pointer (for access to row_count, column_count)
# hist_piece - coarse and fine histogram to update
# pixel_count- pointer to pixel counter for histogram
# last_coord - coordinate and stride of pixel to remove
# coord - coordinate and stride of pixel to add
#
############################################################################
cdef inline void update_histogram(Histograms *ph,
HistogramPiece *hist_piece,
pixel_count_t *pixel_count,
SCoord *last_coord,
SCoord *coord):
cdef:
np.int32_t current_column = ph.current_column
np.int32_t current_row = ph.current_row
np.int32_t current_stride = ph.current_stride
np.int32_t column_count = ph.column_count
np.int32_t row_count = ph.row_count
np.uint8_t value
np.int32_t stride
np.int32_t x
np.int32_t y
x = last_coord.x + current_column
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]):
value = ph.data[stride]
pixel_count[0] -= 1
hist_piece.fine[value] -= 1
hist_piece.coarse[value / 16] -= 1
x = coord.x + current_column
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]):
value = ph.data[stride]
pixel_count[0] += 1
hist_piece.fine[value] += 1
hist_piece.coarse[value / 16] += 1
############################################################################
#
# update_current_location - update the histograms at the current location
#
############################################################################
cdef inline void update_current_location(Histograms *ph):
cdef:
np.int32_t current_column = ph.current_column
np.int32_t radius = ph.radius
np.int32_t top_left_off = tl_br_colidx(ph, current_column)
np.int32_t top_right_off = tr_bl_colidx(ph, current_column)
np.int32_t bottom_left_off = tr_bl_colidx(ph, current_column)
np.int32_t bottom_right_off = tl_br_colidx(ph, current_column)
np.int32_t leading_edge_off = leading_edge_colidx(ph, current_column)
np.int32_t *coarse_histogram
np.int32_t *fine_histogram
np.int32_t last_xoff
np.int32_t last_yoff
np.int32_t last_stride
np.int32_t xoff
np.int32_t yoff
np.int32_t stride
update_histogram(ph, &ph.histogram[top_left_off].top_left,
&ph.pixel_count[top_left_off].top_left,
&ph.last_top_left,
&ph.top_left)
update_histogram(ph, &ph.histogram[top_right_off].top_right,
&ph.pixel_count[top_right_off].top_right,
&ph.last_top_right,
&ph.top_right)
update_histogram(ph, &ph.histogram[bottom_left_off].bottom_left,
&ph.pixel_count[bottom_left_off].bottom_left,
&ph.last_bottom_left,
&ph.bottom_left)
update_histogram(ph, &ph.histogram[bottom_right_off].bottom_right,
&ph.pixel_count[bottom_right_off].bottom_right,
&ph.last_bottom_right,
&ph.bottom_right)
update_histogram(ph, &ph.histogram[leading_edge_off].edge,
&ph.pixel_count[leading_edge_off].edge,
&ph.last_leading_edge,
&ph.leading_edge)
############################################################################
#
# find_median - search the current accumulator for the median
#
############################################################################
cdef inline np.uint8_t find_median(Histograms *ph):
cdef:
np.uint32_t pixels_below # of pixels below the median
int i
int j
int k
np.uint32_t accumulator
if ph.accumulator_count == 0:
return 0
pixels_below = ((ph.accumulator_count * ph.percent + 50)
/ 100) # +50 for roundoff
if pixels_below > 0:
pixels_below -= 1
accumulator = 0
for i in range(16):
accumulator += ph.accumulator.coarse[i]
if accumulator > pixels_below:
break
accumulator -= ph.accumulator.coarse[i]
update_fine(ph, i)
for j in range(i*16,(i+1)*16):
accumulator += ph.accumulator.fine[j]
if accumulator > pixels_below:
return <np.uint8_t> j
return 0
############################################################################
#
# c_median_filter - median filter algorithm
#
# rows - # of rows in each array
# columns - # of columns in each array
# row_stride - stride from one row to the next in each array
# col_stride - stride from one column to the next in each array
# radius - radius of circle inscribed into octagon
# percent - "median" cutoff: 50 = median, 25 = lower quartile, etc
# data - array of image pixels to be filtered
# mask - mask of significant pixels
# output - array to be filled with filtered pixels
#
############################################################################
cdef int c_median_filter(np.int32_t rows,
np.int32_t columns,
np.int32_t row_stride,
np.int32_t col_stride,
np.int32_t radius,
np.int32_t percent,
np.uint8_t *data,
np.uint8_t *mask,
np.uint8_t *output):
cdef:
Histograms *ph
Histogram *phistogram
int row
int col
int i
np.int32_t top_left_off
np.int32_t top_right_off
np.int32_t bottom_left_off
np.int32_t bottom_right_off
ph = allocate_histograms(rows, columns, row_stride, col_stride,
radius, percent, data, mask, output)
if not ph:
return 1
for row in range(-radius, rows):
#
# Initialize the starting diagonal histograms to zero. The leading
# and trailing histograms descend from above and so are initialized
# when memory is initially set to zero. The diagonals move in
# from the left (top left and bottom right) and right (top right
# and bottom left). One of each needs to be initialized at the
# start of each row
#
tl_br_off = tl_br_colidx(ph, -radius)
tr_bl_off = tr_bl_colidx(ph, columns + radius - 1)
memset(&ph.histogram[tl_br_off].top_left, 0, sizeof(HistogramPiece))
memset(&ph.histogram[tl_br_off].bottom_right, 0, sizeof(HistogramPiece))
memset(&ph.histogram[tr_bl_off].top_right, 0, sizeof(HistogramPiece))
memset(&ph.histogram[tr_bl_off].bottom_left, 0, sizeof(HistogramPiece))
ph.pixel_count[tl_br_off].top_left = 0
ph.pixel_count[tl_br_off].bottom_right = 0
ph.pixel_count[tr_bl_off].top_right = 0
ph.pixel_count[tr_bl_off].bottom_left = 0
#
# Initialize the accumulator (octagon histogram) to zero
#
memset(&ph.accumulator, 0, sizeof(ph.accumulator))
ph.accumulator_count = 0
for i in range(16):
ph.last_update_column[i] = -radius-1
#
# Initialize the current stride to the beginning of the row
#
ph.current_row = row
#
# Update locations and coarse accumulator for the octagon
# for points before 0
#
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)
accumulate(ph)
#
# Update locations and coarse accumulator and compute
# the median for points between 0 and "columns"
#
if row >= 0:
for col in range(0, columns):
ph.current_column = col
ph.current_stride = row * row_stride + col * col_stride
update_current_location(ph)
accumulate(ph)
ph.output[ph.current_stride] = find_median(ph)
for col in range(columns, columns+radius):
ph.current_column = col
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):
"""Median filter with octagon shape and masking.
Parameters
----------
data : (M,N) ndarray, dtype uint8
Input image.
mask : (M,N) array, dtype uint8
A value of 1 indicates a significant pixel, 0
that a pixel is masked.
output : (M,N) array, dtype uint8
Array of same size as the input in which to store
the filtered image.
radius : int
Radius of the inscribed circle to the octagon.
percent : int, optional
The unmasked pixels within the octagon are sorted, and the
value at the `percent`-th index chosen. For example, the
default value of 50 chooses the median pixel.
"""
if percent < 0:
raise ValueError('Median filter percent = %d is less than zero' % \
percent)
if percent > 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)' %
(data.shape[0], data.shape[1],
mask.shape[0], mask.shape[1]))
if data.shape[0] != output.shape[0] or data.shape[1] != output.shape[1]:
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],
radius, percent,
<np.uint8_t *>data.data,
<np.uint8_t *>mask.data,
<np.uint8_t *>output.data):
raise MemoryError('Failed to allocate scratchpad memory')
scikits/image/filter/ctmf.py
+86
View File
@@ -0,0 +1,86 @@
'''ctmf.py - constant time per pixel median filtering with an octagonal shape
Reference: S. Perreault and P. Hebert, "Median Filtering in Constant Time",
IEEE Transactions on Image Processing, September 2007.
Originally part of CellProfiler, code licensed under both GPL and BSD licenses.
Website: http://www.cellprofiler.org
Copyright (c) 2003-2009 Massachusetts Institute of Technology
Copyright (c) 2009-2011 Broad Institute
All rights reserved.
Original author: Lee Kamentsky
'''
import numpy as np
import _ctmf
from rank_order import rank_order
def median_filter(data, mask=None, radius=1, percent=50):
'''Masked median filter with octagon shape.
Parameters
----------
data : (M,N) ndarray, dtype uint8
Input image.
mask : (M,N) ndarray, dtype uint8, optional
A value of 1 indicates a significant pixel, 0
that a pixel is masked. By default, all pixels
are considered.
radius : {int, 1}, optional
The radius of a circle inscribed into the filtering
octagon. Default radius is 1.
percent : {int, 50}, optional
The unmasked pixels within the octagon are sorted, and the
value at the `percent`-th index chosen. For example, the
default value of 50 chooses the median pixel.
Returns
-------
out : (M,N) ndarray, dtype uint8
Filtered array. In areas where the median filter does
not overlap the mask, the filtered result is underfined, but
in practice, it will be the lowest value in the valid area.
'''
if mask is None:
mask = np.ones(data.shape, dtype=np.bool)
mask = np.ascontiguousarray(mask, dtype=np.bool)
if np.all(~ mask):
return data.copy()
#
# Normalize the ranked data to 0-255
#
if (not np.issubdtype(data.dtype, np.int) or
np.min(data) < 0 or np.max(data) > 255):
ranked_data,translation = rank_order(data[mask])
max_ranked_data = np.max(ranked_data)
if max_ranked_data == 0:
return data
if max_ranked_data > 255:
ranked_data = ranked_data * 255 / max_ranked_data
was_ranked = True
else:
ranked_data = data[mask]
was_ranked = False
input = np.zeros(data.shape, np.uint8 )
input[mask] = ranked_data
mask.dtype = np.uint8
output = np.zeros(data.shape, np.uint8)
_ctmf.median_filter(input, mask, output, radius, percent)
if was_ranked:
#
# The translation gives the original value at each ranking.
# We rescale the output to the original ranking and then
# use the translation to look up the original value in the data.
#
if max_ranked_data > 255:
result = translation[output.astype(np.uint32) * max_ranked_data / 255]
else:
result = translation[output]
else:
result = output
return result
scikits/image/filter/rank_order.py
+32
View File
@@ -0,0 +1,32 @@
"""rankorder.py - convert an image of any type to an image of ints whose
pixels have an identical rank order compared to the original image
Originally part of CellProfiler, code licensed under both GPL and BSD licenses.
Website: http://www.cellprofiler.org
Copyright (c) 2003-2009 Massachusetts Institute of Technology
Copyright (c) 2009-2011 Broad Institute
All rights reserved.
Original author: Lee Kamentstky
"""
import numpy
def rank_order(image):
"""Return an image of the same shape where each pixel has the
rank-order value of the corresponding pixel in the image.
The returned image's elements are of type numpy.uint32 which
simplifies processing in C code.
"""
flat_image = image.ravel()
sort_order = flat_image.argsort().astype(numpy.uint32)
flat_image = flat_image[sort_order]
sort_rank = numpy.zeros_like(sort_order)
is_different = flat_image[:-1] != flat_image[1:]
numpy.cumsum(is_different, out=sort_rank[1:])
original_values = numpy.zeros((sort_rank[-1]+1,),image.dtype)
original_values[0] = flat_image[0]
original_values[1:] = flat_image[1:][is_different]
int_image = numpy.zeros_like(sort_order)
int_image[sort_order] = sort_rank
return (int_image.reshape(image.shape), original_values)
scikits/image/filter/setup.py
+33
View File
@@ -0,0 +1,33 @@
#!/usr/bin/env python
import os
import shutil
import hashlib
from scikits.image._build import cython
base_path = os.path.abspath(os.path.dirname(__file__))
def configuration(parent_package='', top_path=None):
from numpy.distutils.misc_util import Configuration, get_numpy_include_dirs
config = Configuration('filter', parent_package, top_path)
config.add_data_dir('tests')
cython(['_ctmf.pyx'], working_path=base_path)
config.add_extension('_ctmf', sources=['_ctmf.c'],
include_dirs=[get_numpy_include_dirs()])
return config
if __name__ == '__main__':
from numpy.distutils.core import setup
setup(maintainer = 'Scikits.Image Developers',
author = 'Scikits.Image Developers',
maintainer_email = 'scikits-image@googlegroups.com',
description = 'Filters',
url = 'http://stefanv.github.com/scikits.image/',
license = 'SciPy License (BSD Style)',
**(configuration(top_path='').todict())
)
scikits/image/filter/tests/test_ctmf.py
+105
View File
@@ -0,0 +1,105 @@
import os.path
import numpy as np
from numpy.testing import *
from scikits.image.filter import median_filter
def test_00_00_zeros():
'''The median filter on an array of all zeros should be zero'''
result = median_filter(np.zeros((10,10)), np.ones((10,10),bool), 3)
assert np.all(result == 0)
def test_00_01_all_masked():
'''Test a completely masked image
Regression test of IMG-1029'''
result = median_filter(np.zeros((10,10)), np.zeros((10,10), bool), 3)
assert (np.all(result == 0))
def test_00_02_all_but_one_masked():
mask = np.zeros((10,10), bool)
mask[5,5] = True
result = median_filter(np.zeros((10,10)), mask, 3)
def test_01_01_mask():
'''The median filter, masking a single value'''
img = np.zeros((10,10))
img[5,5] = 1
mask = np.ones((10,10),bool)
mask[5,5] = False
result = median_filter(img, mask, 3)
assert (np.all(result[mask] == 0))
assert_equal(result[5,5], 1)
def test_02_01_median():
'''A median filter larger than the image = median of image'''
np.random.seed(0)
img = np.random.uniform(size=(9,9))
result = median_filter(img, np.ones((9,9),bool), 20)
assert_equal(result[0,0], np.median(img))
assert (np.all(result == np.median(img)))
def test_02_02_median_bigger():
'''Use an image of more than 255 values to test approximation'''
np.random.seed(0)
img = np.random.uniform(size=(20,20))
result = median_filter(img, np.ones((20,20),bool),40)
sorted = np.ravel(img)
sorted.sort()
min_acceptable = sorted[198]
max_acceptable = sorted[202]
assert (np.all(result >= min_acceptable))
assert (np.all(result <= max_acceptable))
def test_03_01_shape():
'''Make sure the median filter is the expected octagonal shape'''
radius = 5
a_2 = int(radius / 2.414213)
i,j = np.mgrid[-10:11,-10:11]
octagon = np.ones((21,21), bool)
#
# constrain the octagon mask to be the points that are on
# the correct side of the 8 edges
#
octagon[i < -radius] = False
octagon[i > radius] = False
octagon[j < -radius] = False
octagon[j > radius] = False
octagon[i+j < -radius-a_2] = False
octagon[j-i > radius+a_2] = False
octagon[i+j > radius+a_2] = False
octagon[i-j > radius+a_2] = False
np.random.seed(0)
img = np.random.uniform(size=(21,21))
result = median_filter(img, np.ones((21,21),bool), radius)
sorted = img[octagon]
sorted.sort()
min_acceptable = sorted[len(sorted)/2-1]
max_acceptable = sorted[len(sorted)/2+1]
assert (result[10,10] >= min_acceptable)
assert (result[10,10] <= max_acceptable)
def test_04_01_half_masked():
'''Make sure that the median filter can handle large masked areas.'''
img = np.ones((20, 20))
mask = np.ones((20, 20),bool)
mask[10:, :] = False
img[~ mask] = 2
img[1, 1] = 0 # to prevent short circuit for uniform data.
result = median_filter(img, mask, 5)
# in partial coverage areas, the result should be only from the masked pixels
assert (np.all(result[:14, :] == 1))
# in zero coverage areas, the result should be the lowest valud in the valid area
assert (np.all(result[15:, :] == np.min(img[mask])))
def test_default_values():
img = (np.random.random((20, 20)) * 255).astype(np.uint8)
mask = np.ones((20, 20), dtype=np.uint8)
result1 = median_filter(img, mask, radius=1, percent=50)
result2 = median_filter(img)
assert_array_equal(result1, result2)
if __name__ == "__main__":
run_module_suite()
scikits/image/setup.py
+1
View File
@@ -9,6 +9,7 @@ def configuration(parent_package='', top_path=None):
config.add_subpackage('graph')
config.add_subpackage('io')
config.add_subpackage('morphology')
config.add_subpackage('filter')
def add_test_directories(arg, dirname, fnames):
if dirname.split(os.path.sep)[-1] == 'tests':