1"""
2========================
3Sliding window histogram
4========================
5
6Histogram matching can be used for object detection in images [1]_. This
7example extracts a single coin from the ``skimage.data.coins`` image and uses
8histogram matching to attempt to locate it within the original image.
9
10First, a box-shaped region of the image containing the target coin is
11extracted and a histogram of its grayscale values is computed.
12
13Next, for each pixel in the test image, a histogram of the grayscale values in
14a region of the image surrounding the pixel is computed.
15``skimage.filters.rank.windowed_histogram`` is used for this task, as it employs
16an efficient sliding window based algorithm that is able to compute these
17histograms quickly [2]_. The local histogram for the region surrounding each
18pixel in the image is compared to that of the single coin, with a similarity
19measure being computed and displayed.
20
21The histogram of the single coin is computed using ``numpy.histogram`` on a box
22shaped region surrounding the coin, while the sliding window histograms are
23computed using a disc shaped structural element of a slightly different size.
24This is done in aid of demonstrating that the technique still finds similarity
25in spite of these differences.
26
27To demonstrate the rotational invariance of the technique, the same test is
28performed on a version of the coins image rotated by 45 degrees.
29
30References
31----------
32.. [1] Porikli, F. "Integral Histogram: A Fast Way to Extract Histograms
33       in Cartesian Spaces" CVPR, 2005. Vol. 1. IEEE, 2005
34.. [2] S.Perreault and P.Hebert. Median filtering in constant time.
35       Trans. Image Processing, 16(9):2389-2394, 2007.
36"""
37
38import numpy as np
39import matplotlib
40import matplotlib.pyplot as plt
41
42from skimage import data, transform
43from skimage.util import img_as_ubyte
44from skimage.morphology import disk
45from skimage.filters import rank
46
47
48matplotlib.rcParams['font.size'] = 9
49
50
51def windowed_histogram_similarity(image, footprint, reference_hist, n_bins):
52    # Compute normalized windowed histogram feature vector for each pixel
53    px_histograms = rank.windowed_histogram(image, footprint, n_bins=n_bins)
54
55    # Reshape coin histogram to (1,1,N) for broadcast when we want to use it in
56    # arithmetic operations with the windowed histograms from the image
57    reference_hist = reference_hist.reshape((1, 1) + reference_hist.shape)
58
59    # Compute Chi squared distance metric: sum((X-Y)^2 / (X+Y));
60    # a measure of distance between histograms
61    X = px_histograms
62    Y = reference_hist
63
64    num = (X - Y) ** 2
65    denom = X + Y
66    denom[denom == 0] = np.infty
67    frac = num / denom
68
69    chi_sqr = 0.5 * np.sum(frac, axis=2)
70
71    # Generate a similarity measure. It needs to be low when distance is high
72    # and high when distance is low; taking the reciprocal will do this.
73    # Chi squared will always be >= 0, add small value to prevent divide by 0.
74    similarity = 1 / (chi_sqr + 1.0e-4)
75
76    return similarity
77
78
79# Load the `skimage.data.coins` image
80img = img_as_ubyte(data.coins())
81
82# Quantize to 16 levels of grayscale; this way the output image will have a
83# 16-dimensional feature vector per pixel
84quantized_img = img // 16
85
86# Select the coin from the 4th column, second row.
87# Co-ordinate ordering: [x1,y1,x2,y2]
88coin_coords = [184, 100, 228, 148]   # 44 x 44 region
89coin = quantized_img[coin_coords[1]:coin_coords[3],
90                     coin_coords[0]:coin_coords[2]]
91
92# Compute coin histogram and normalize
93coin_hist, _ = np.histogram(coin.flatten(), bins=16, range=(0, 16))
94coin_hist = coin_hist.astype(float) / np.sum(coin_hist)
95
96# Compute a disk shaped mask that will define the shape of our sliding window
97# Example coin is ~44px across, so make a disk 61px wide (2 * rad + 1) to be
98# big enough for other coins too.
99footprint = disk(30)
100
101# Compute the similarity across the complete image
102similarity = windowed_histogram_similarity(quantized_img, footprint, coin_hist,
103                                           coin_hist.shape[0])
104
105# Now try a rotated image
106rotated_img = img_as_ubyte(transform.rotate(img, 45.0, resize=True))
107# Quantize to 16 levels as before
108quantized_rotated_image = rotated_img // 16
109# Similarity on rotated image
110rotated_similarity = windowed_histogram_similarity(quantized_rotated_image,
111                                                   footprint, coin_hist,
112                                                   coin_hist.shape[0])
113
114fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(10, 10))
115
116axes[0, 0].imshow(quantized_img, cmap='gray')
117axes[0, 0].set_title('Quantized image')
118axes[0, 0].axis('off')
119
120axes[0, 1].imshow(coin, cmap='gray')
121axes[0, 1].set_title('Coin from 2nd row, 4th column')
122axes[0, 1].axis('off')
123
124axes[1, 0].imshow(img, cmap='gray')
125axes[1, 0].imshow(similarity, cmap='hot', alpha=0.5)
126axes[1, 0].set_title('Original image with overlaid similarity')
127axes[1, 0].axis('off')
128
129axes[1, 1].imshow(rotated_img, cmap='gray')
130axes[1, 1].imshow(rotated_similarity, cmap='hot', alpha=0.5)
131axes[1, 1].set_title('Rotated image with overlaid similarity')
132axes[1, 1].axis('off')
133
134plt.tight_layout()
135plt.show()
136