diff --git a/doc/examples/plot_windowed_histogram.py b/doc/examples/plot_windowed_histogram.py new file mode 100644 index 00000000..7b0bcb3a --- /dev/null +++ b/doc/examples/plot_windowed_histogram.py @@ -0,0 +1,129 @@ +""" +======================== +Sliding window histogram +======================== + +This example extracts a single coin from the coins image and generates a +histogram of its greyscale values. + +It then computes a sliding window histogram of the complete image using +rank.windowed_histogram. The local histogram for the region surrounding +each pixel in the image is compared to that of the single coin, with +a similarity measure being computed and displayed. + +To demonstrate the rotational invariance of the technique, the same +test is performed on a version of the coins image rotated by 45 degrees. +""" +import numpy as np +import matplotlib +import matplotlib.pyplot as plt + +from skimage import data +from skimage.util.dtype import dtype_range +from skimage.util import img_as_ubyte +from skimage import exposure +from skimage.morphology import disk +from skimage.filter import rank +from skimage import transform + + +matplotlib.rcParams['font.size'] = 9 + + +def windowed_histogram_similarity(image, selem, reference_hist, n_bins): + # Compute normalized windowed histogram feature vector for each pixel + px_histograms = rank.windowed_histogram(image, selem, n_bins=n_bins) + + # Reshape coin histogram to (1,1,N) for broadcast when we want to use it in + # arithmetic operations with the windowed histograms fro the image + reference_hist = reference_hist.reshape((1,1) + reference_hist.shape) + + # Compute Chi squared distance metric: sum((X-Y)**2 / (X+Y); + # a measure of distance between histograms + X = px_histograms + Y = reference_hist + num = (X-Y)*(X-Y) + denom = X+Y + frac = num / denom + frac[denom==0] = 0 + chi_sqr = np.sum(frac, axis=2) * 0.5 + + # Generate a similarity measure. It needs to be low when distance is high. + # and high when distance is low; taking the reciprocal will do this. + # Chi squared will always be >= 0. Add small value to prevent divide by 0. + # Square the denominator to push low values toward 0; this makes the + # high similarity regions stand out in the figure created below; this + # us just done for aesthetics. + similarity = 1 / (chi_sqr + 1.0e-6)**2 + + return similarity + + +# Load the coins image +img = img_as_ubyte(data.coins()) +# img = img_as_ubyte(plt.imread('../../skimage/data/coins.png')) + +# Quantize to 16 levels of grayscale; this way the output image will have a +# 16-dimensional feature vector per pixel +quantized_img = img/16 + +# Select the coin from the 4th column, second row. +# Co-ordinate ordering: [x1,y1,x2,y2] +coin_coords = [184,100,228,148] # 44 x 44 region +coin = quantized_img[coin_coords[1]:coin_coords[3], coin_coords[0]:coin_coords[2]] + +# Compute coin histogram and normalize +coin_hist, _ = np.histogram(coin.flatten(), bins=16, range=(0,16)) +coin_hist = coin_hist.astype(float) / np.sum(coin_hist) + + +# Compute a disk shaped mask that will define the shape of our sliding window +# Example coin is ~44px across, so make a disk 61px wide (2*rad+1) to be big +# enough for other coins too. +selem = disk(30) + + +# Compute the similarity across the complete image +similarity = windowed_histogram_similarity(quantized_img, selem, coin_hist, + coin_hist.shape[0]) + +# Now try a rotated image +rotated_img = img_as_ubyte(transform.rotate(img, 45.0, resize=True)) +# Quantize to 16 levels as before +quantized_rotated_image = rotated_img/16 +# Similarity on rotated image +rotated_similarity = windowed_histogram_similarity(quantized_rotated_image, + selem, coin_hist, + coin_hist.shape[0]) + + + +# Plot it all +fig, axes = plt.subplots(nrows=5, figsize=(6, 18)) +ax0, ax1, ax2, ax3, ax4 = axes + +ax0.imshow(img, cmap='gray') +ax0.set_title('Original image') +ax0.axis('off') + +ax1.imshow(quantized_img, cmap='gray') +ax1.set_title('Quantized image') +ax1.axis('off') + +ax2.imshow(coin, cmap='gray') +ax2.set_title('Coin from 2nd row, 4th column') +ax2.axis('off') + +ax3.imshow(img, cmap='gray') +# While jet is not a great colormap, it makes the high similarity areas +# stand out +ax3.imshow(similarity, cmap='jet', alpha=0.5) +ax3.set_title('Original image with overlayed similarity') +ax3.axis('off') + +ax4.imshow(rotated_img, cmap='gray') +ax4.imshow(rotated_similarity, cmap='jet', alpha=0.5) +ax4.set_title('Rotated image with overlayed similarity') +ax4.axis('off') + +plt.show() diff --git a/skimage/filter/rank/generic.py b/skimage/filter/rank/generic.py index 487d34b5..2af2428e 100644 --- a/skimage/filter/rank/generic.py +++ b/skimage/filter/rank/generic.py @@ -868,8 +868,8 @@ def otsu(image, selem, out=None, mask=None, shift_x=False, shift_y=False): mask=mask, shift_x=shift_x, shift_y=shift_y) -def windowed_histogram(image, selem, out=None, mask=None, shift_x=False, shift_y=False): - """Sliding window histogram +def windowed_histogram(image, selem, out=None, mask=None, shift_x=False, shift_y=False, n_bins=None): + """Normalized sliding window histogram Parameters ---------- @@ -886,15 +886,19 @@ def windowed_histogram(image, selem, out=None, mask=None, shift_x=False, shift_y Offset added to the structuring element center point. Shift is bounded to the structuring element sizes (center must be inside the given structuring element). + n_bins : int or None + The number of histogram bins. Will default to image.max() + 1 if None + is passed. Returns ------- - out : 3-D array (same dtype as input image) whose extra dimension - Output image. - - References - ---------- - .. [otsu] http://en.wikipedia.org/wiki/Otsu's_method + out : 3-D array with float dtype of dimensions (H,W,N), where (H,W) are + the dimensions of the input image and N is n_bins or image.max()+1 + if no value is provided as a parameter. Effectively, each pixel + is an N-dimensional feature vector that is the histogram. + The sum of the elements in the feature vector will be 1, unless + no pixels in the window were covered by both selem and mask, in which + case all elements will be 0. Examples -------- @@ -903,9 +907,14 @@ def windowed_histogram(image, selem, out=None, mask=None, shift_x=False, shift_y >>> from skimage.morphology import disk >>> img = data.camera() >>> hist_img = windowed_histogram(img, disk(5)) - >>> thresh_image = img >= local_otsu """ - return _apply_vector_per_pixel(generic_cy._windowed_hist, image, selem, out=out, - mask=mask, shift_x=shift_x, shift_y=shift_y, pixel_size=image.max()+1) + if n_bins is None: + n_bins = image.max() + 1 + + return _apply_vector_per_pixel(generic_cy._windowed_hist, image, selem, + out=out, mask=mask, + shift_x=shift_x, shift_y=shift_y, + out_dtype=np.double, + pixel_size=n_bins) diff --git a/skimage/filter/rank/generic_cy.pyx b/skimage/filter/rank/generic_cy.pyx index 88766f44..bafd2432 100644 --- a/skimage/filter/rank/generic_cy.pyx +++ b/skimage/filter/rank/generic_cy.pyx @@ -378,8 +378,14 @@ cdef inline void _kernel_win_hist(dtype_t_out[:] out, Py_ssize_t* histo, Py_ssize_t s0, Py_ssize_t s1): cdef Py_ssize_t i cdef Py_ssize_t max_i - for i in xrange(out.shape[0]): - out[i] = histo[i] + cdef double scale + if pop: + scale = 1.0 / pop + for i in xrange(out.shape[0]): + out[i] = (histo[i] * scale) + else: + for i in xrange(out.shape[0]): + out[i] = 0 diff --git a/skimage/filter/rank/tests/test_rank.py b/skimage/filter/rank/tests/test_rank.py index 390a7790..48d02462 100644 --- a/skimage/filter/rank/tests/test_rank.py +++ b/skimage/filter/rank/tests/test_rank.py @@ -551,22 +551,35 @@ def test_windowed_histogram(): [0, 1, 1, 1, 0], [0, 0, 0, 0, 0]], dtype=np.uint8) elem = np.ones((3, 3), dtype=np.uint8) - out8 = np.empty(image8.shape+(2,), dtype=np.uint8) + outf = np.empty(image8.shape+(2,), dtype=float) mask = np.ones(image8.shape, dtype=np.uint8) + # Population so we can normalize the expected output while maintaining + # code readability + pop = np.array([[4, 6, 6, 6, 4], + [6, 9, 9, 9, 6], + [6, 9, 9, 9, 6], + [6, 9, 9, 9, 6], + [4, 6, 6, 6, 4]], dtype=float) + r0 = np.array([[3, 4, 3, 4, 3], - [4, 5, 3, 5, 4], - [3, 3, 0, 3, 3], - [4, 5, 3, 5, 4], - [3, 4, 3, 4, 3]], dtype=np.uint8) + [4, 5, 3, 5, 4], + [3, 3, 0, 3, 3], + [4, 5, 3, 5, 4], + [3, 4, 3, 4, 3]], dtype=float) / pop r1 = np.array([[1, 2, 3, 2, 1], - [2, 4, 6, 4, 2], - [3, 6, 9, 6, 3], - [2, 4, 6, 4, 2], - [1, 2, 3, 2, 1]], dtype=np.uint8) - rank.windowed_histogram(image=image8, selem=elem, out=out8, mask=mask) - assert_array_equal(r0, out8[:,:,0]) - assert_array_equal(r1, out8[:,:,1]) + [2, 4, 6, 4, 2], + [3, 6, 9, 6, 3], + [2, 4, 6, 4, 2], + [1, 2, 3, 2, 1]], dtype=float) / pop + rank.windowed_histogram(image=image8, selem=elem, out=outf, mask=mask) + assert_array_equal(r0, outf[:,:,0]) + assert_array_equal(r1, outf[:,:,1]) + + # Test n_bins parameter + larger_output = rank.windowed_histogram(image=image8, selem=elem, + mask=mask, n_bins=5) + assert larger_output.shape[2] == 5 if __name__ == "__main__":