mirror of
https://github.com/wassname/scikit-image.git
synced 2026-07-09 04:48:03 +08:00
Fixed some issues with the new windowed_histogram function in filter.rank.
It used to be able to output uint8 histogram, whose max pixel counts of 255 could easily overflow. Addressed by limiting the output type to float. Normalized histograms are now generated, otherwise behaviour is unpredictable at boundaries or when pixels are not permitted by a mask. The new optional n_bins parameter allows the caller to specify the size of the histogram generated. Having it fixed to image.max()+1 could result in feature vectors being shorter than desired just due to a value not being used in an image. An example has been added to the docs, that demonstrates the application of windows histograms in object matching; a single coin is extracted and found by chi squared histogram matching.
This commit is contained in:
Geoffrey French
@@ -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()
|
||||
@@ -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)
|
||||
|
||||
@@ -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] = <dtype_t_out>histo[i]
|
||||
cdef double scale
|
||||
if pop:
|
||||
scale = 1.0 / pop
|
||||
for i in xrange(out.shape[0]):
|
||||
out[i] = <dtype_t_out>(histo[i] * scale)
|
||||
else:
|
||||
for i in xrange(out.shape[0]):
|
||||
out[i] = <dtype_t_out>0
|
||||
|
||||
|
||||
|
||||
|
||||
@@ -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__":
|
||||
|
||||
Reference in New Issue
Block a user