diff --git a/CONTRIBUTORS.txt b/CONTRIBUTORS.txt index 590da154..d5496943 100644 --- a/CONTRIBUTORS.txt +++ b/CONTRIBUTORS.txt @@ -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 diff --git a/scikits/image/filter/__init__.py b/scikits/image/filter/__init__.py index 4bdff030..afcc9db5 100644 --- a/scikits/image/filter/__init__.py +++ b/scikits/image/filter/__init__.py @@ -1,2 +1,2 @@ from lpi_filter import * - +from ctmf import median_filter diff --git a/scikits/image/filter/_ctmf.pyx b/scikits/image/filter/_ctmf.pyx new file mode 100644 index 00000000..c133d8d6 --- /dev/null +++ b/scikits/image/filter/_ctmf.pyx @@ -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 = ptr + if not ptr: + return ph + ph.memory = ptr + 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 += 31 + roundoff -= roundoff % 32 + ptr = roundoff + ph.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 = (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 + +############################################################################ +# +# _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 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, + data.data, + mask.data, + output.data): + raise MemoryError('Failed to allocate scratchpad memory') diff --git a/scikits/image/filter/ctmf.py b/scikits/image/filter/ctmf.py new file mode 100644 index 00000000..952fbc62 --- /dev/null +++ b/scikits/image/filter/ctmf.py @@ -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 diff --git a/scikits/image/filter/rank_order.py b/scikits/image/filter/rank_order.py new file mode 100644 index 00000000..399641f8 --- /dev/null +++ b/scikits/image/filter/rank_order.py @@ -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) + + diff --git a/scikits/image/filter/setup.py b/scikits/image/filter/setup.py new file mode 100644 index 00000000..d211d7f4 --- /dev/null +++ b/scikits/image/filter/setup.py @@ -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()) + ) diff --git a/scikits/image/filter/tests/test_ctmf.py b/scikits/image/filter/tests/test_ctmf.py new file mode 100644 index 00000000..02795b3c --- /dev/null +++ b/scikits/image/filter/tests/test_ctmf.py @@ -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() diff --git a/scikits/image/setup.py b/scikits/image/setup.py index 7d1b8a9f..d0f08285 100644 --- a/scikits/image/setup.py +++ b/scikits/image/setup.py @@ -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':