<jupyter_start><jupyter_text>Tutorial 2: Feature detectionThis tutorial was adapted from https://github.com/scikit-image/skimage-tutorials/blob/master/lectures/2_feature_detection.ipynb. Feature detection is often the end result of image processing. We'll detect some basic features (edges, points, and circles) in this section, but there are a wealth of feature detectors that are available.  Edge detection Before we start, let's set the default colormap to grayscale and turn off pixel interpolation.<jupyter_code>import matplotlib.pyplot as plt
import numpy as np

plt.rcParams['image.cmap'] = 'gray'
plt.rcParams['image.interpolation'] = 'none'<jupyter_output><empty_output><jupyter_text>We've already discussed edge filtering, using the Sobel filter, in the last section.<jupyter_code>from skimage import data
from skimage import filters

image = data.camera()
pixelated = image[::10, ::10]
gradient = filters.sobel(pixelated)
fig, ax = plt.subplots(ncols=2, figsize=(20, 10))
ax[0].imshow(pixelated)
ax[1].imshow(gradient)<jupyter_output><empty_output><jupyter_text>With the Sobel filter, however, we get back a grayscale image, which essentially tells us the likelihood that a pixel is on the edge of an object.We can apply a bit more logic to *detect* an edge; i.e. we can use that filtered image to make a *decision* whether or not a pixel is on an edge. The simplest way to do that is with thresholding:<jupyter_code>fig, ax = plt.subplots(ncols=2, figsize=(20, 10))
ax[0].imshow(gradient)
ax[1].imshow(gradient > 0.4)<jupyter_output><empty_output><jupyter_text>That approach doesn't do a great job. It's noisy and produces thick edges. Furthermore, it doesn't use our *knowledge* of how edges work: They should be thin and tend to be connected along the direction of the edge.  Canny edge detector The Canny edge detector combines the Sobel filter with a few other steps to give a binary edge image. The steps are as follows:* Gaussian filter* Sobel filter* Non-maximal suppression* Hysteresis thresholding Let's go through these steps one at a time.  Step 1: Gaussian filter As discussed earlier, gradients tend to enhance noise. To combat this effect, we first smooth the image using a gradient filter:<jupyter_code>from skimage import img_as_float

sigma = 1  # Standard-deviation of Gaussian; larger smooths more.
pixelated_float = img_as_float(pixelated)
pixelated_float = pixelated
smooth = filters.gaussian(pixelated_float, sigma)
fig, ax = plt.subplots(ncols=2, figsize=(20, 10))
ax[0].imshow(pixelated_float)
ax[1].imshow(smooth)<jupyter_output><empty_output><jupyter_text>Step 2: Sobel filter Next, we apply our edge filter:<jupyter_code>gradient_magnitude = filters.sobel(smooth)
fig, ax = plt.subplots(ncols=2, figsize=(20, 10))
ax[0].imshow(smooth)
ax[1].imshow(gradient_magnitude)<jupyter_output><empty_output><jupyter_text>Step 3: Non-maximal suppression Goal: Suppress gradients that aren't on an edge Ideally, an edge is thin: In some sense, an edge is infinitely thin, but images are discrete so we want edges to be a single pixel wide. To accomplish this, we thin the image using "non-maximal suppression". This takes the edge-filtered image and thins the response *across* the edge direction; i.e. in the direction of the maximum gradient:<jupyter_code>zoomed_grad = gradient_magnitude[15:25, 5:15]
maximal_mask = np.zeros_like(zoomed_grad)
# This mask is made up for demo purposes
maximal_mask[range(10), (7, 6, 5, 4, 3, 2, 2, 2, 3, 3)] = 1
grad_along_edge = maximal_mask * zoomed_grad
fig, ax = plt.subplots(ncols=2, figsize=(20, 10))
ax[0].imshow(zoomed_grad)
ax[1].imshow(grad_along_edge)<jupyter_output><empty_output><jupyter_text>Obviously, this is a faked version of non-maximal suppression: Pixels are *manually* masked here. The actual algorithm detects the direction of edges, and keeps a pixel only if it has a locally maximal gradient-magnitude in the direction *normal to the edge direction*. It doesn't mask pixels *along* the edge direction since an adjacent edge pixel will be of comparable magnitude.The result of the filter is that an edge is only possible if there are no better edges near it.  Step 4: Hysteresis thresholding Goal: Prefer pixels that are connected to edges The final step is the actual decision-making process. Here, we have two parameters: The low threshold and the high threshold. The high threshold sets the gradient value that you *know* is definitely an edge. The low threshold sets the gradient value that could be an edge, but only if it's connected to a pixel that we know is an edge.These two thresholds are displayed below:<jupyter_code>from skimage import color

low_threshold = 0.2
high_threshold = 0.3
label_image = np.zeros_like(pixelated)
# This uses `gradient_magnitude` which has NOT gone through non-maximal-suppression.
label_image[gradient_magnitude > low_threshold] = 1
label_image[gradient_magnitude > high_threshold] = 2
demo_image = color.label2rgb(label_image, gradient_magnitude,
                             bg_label=0, colors=('yellow', 'red'))
plt.imshow(demo_image)<jupyter_output><empty_output><jupyter_text>The **red points** here are above `high_threshold` and are seed points for edges. The **yellow points** are edges if connected (possibly by other yellow points) to seed points; i.e. isolated groups of yellow points will not be detected as edges.Note that the demo above is on the edge image *before* non-maximal suppression, but in reality, this would be done on the image *after* non-maximal suppression. There isn't currently an easy way to get at the intermediate result.  All-in-one: The Canny edge detector Now we're ready to look at the actual result:<jupyter_code>canny_edge = feature.canny(pixelated)
fig, ax = plt.subplots(ncols=2, figsize=(20, 10))
ax[0].imshow(pixelated)
ax[1].imshow(canny_edge)<jupyter_output><empty_output><jupyter_text>Hough transforms Hough transforms are a general class of operations that make up a step in feature detection. Just like we saw with edge detection, Hough transforms produce a result that we can use to detect a feature. (The distinction between the "filter" that we used for edge detection and the "transform" that we use here is a bit arbitrary.)  Circle detection To explore the Hough transform, let's take the *circular* Hough transform as our example. Let's grab an image with some circles:<jupyter_code>image = data.coins()[0:95, 180:370]
plt.imshow(image);<jupyter_output><empty_output><jupyter_text>We can use the Canny edge filter to get a pretty good representation of these circles:<jupyter_code>edges = feature.canny(image, sigma=3, low_threshold=10, high_threshold=60)
plt.imshow(edges);<jupyter_output><empty_output><jupyter_text>While it looks like we've extracted circles, this doesn't give us much if what we want to do is *measure* these circles. For example, what are the size and position of the circles in the above image? The edge image doesn't really tell us much about that.We'll use the Hough transform to extract circle positions and radii:<jupyter_code>from skimage.transform import hough_circle
 
hough_radii = np.arange(15, 30, 2)
hough_response = hough_circle(edges, hough_radii)<jupyter_output><empty_output><jupyter_text>Here, the circular Hough transform actually uses the edge image from before. We also have to define the radii that we want to search for in our image.So... what's the actual result that we get back?<jupyter_code>print(edges.shape)
print(hough_response.shape)<jupyter_output><empty_output><jupyter_text>We can see that the last two dimensions of the response are exactly the same as the original image, so the response is image-like. Then what does the first dimension correspond to?As always, you can get a better feel for the result by plotting it:<jupyter_code># Use max value to intelligently rescale the data for plotting.
h_max = hough_response.max()
fig, axes = plt.subplots(ncols=4, nrows=2, figsize=(20, 6))
for h, r, ax in zip(hough_response, hough_radii, axes.ravel()):
    ax.imshow(h, vmax=0.5*h_max)
    ax.axis('off')
    ax.set_title('radius %.0f' % r)
fig.tight_layout()<jupyter_output><empty_output><jupyter_text>Exercise: Use the response from the Hough transform to **detect the position of each coin**.<jupyter_code>from skimage.feature import peak_local_max

# A list of (row, column) tuples indicating 
# the coordinates of local peaks:
centers = []
# A list of likelihood values for each of the
# local peaks above:
likelihood = []

# Hint
# ----
# To add a single element to a list, use `append`:
# In [0]: mylist = [2]
# ...     mylist.append(3)
# Out[0]: [2, 3]
#
# To add multiple elements to a list, use `extend`:
# In [1]: mylist = [2]
# ...     mylist.extend([3, 4])
# Out[1]: [2, 3, 4]

for h in hough_response:
    # Your code here...
    pass<jupyter_output><empty_output><jupyter_text>Advanced Exercise:Now repeat the above procedure for a more realistic image; e.g. detecting cell nuclei:<jupyter_code>from skimage import io
img = io.imread('../img/tissue.jpg', as_grey=True)
plt.imshow(img);
edges = feature.canny(img)
plt.imshow(edges);<jupyter_output><empty_output><jupyter_text>Need to tweak Canny's optional arguments... have a look at `low_threshold` and `high_threshold`, too...<jupyter_code>edges = feature.canny(img, sigma=3)
plt.imshow(edges);
# Your code here<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.08806848582129481, 'fraction_numerical': 0.01337613697164259, 'mean_word_length': 5.521981856245638, 'pattern_counts': {'":': 0, '<': 65, '<?xml version=': 0, '>': 68, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '16', 'dir': 'jupyter-structured-clean-dedup', 'id': '474512', 'license': 'MIT', 'n_tokens_mistral': 2996, 'n_tokens_neox': 2809, 'n_words': 1126, 'path': 'tutorials/2-skimage_feature_detection.ipynb', 'repo_name': 'mbeyeler/2018-neurohack-cv', 'size': '9345'}