{
"cells": [
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"from __future__ import division, print_function\n",
"%matplotlib inline"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Color and exposure"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt\n",
"import numpy as np"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Exercise:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Create three images; each should display a red, green, or blue channel of the original image."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"import skdemo\n",
"from skimage import data\n",
"color_image = data.chelsea()\n",
"\n",
"red_image = np.zeros_like(color_image)\n",
"green_image = np.zeros_like(color_image)\n",
"blue_image = np.zeros_like(color_image)\n",
"\n",
"red_image[:, :, 0] = color_image[:, :, 0]\n",
"green_image[:, :, 1] = color_image[:, :, 1]\n",
"blue_image[:, :, 2] = color_image[:, :, 2]\n",
"\n",
"skdemo.imshow_all(red_image, green_image, blue_image)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Alternatively, we can \"unpack\" our RGB image into red, green, and blue channels using `np.rollaxis`"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"np.rollaxis(color_image, axis=-1).shape"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The `axis` argument tells `rollaxis` which axis to... ahem... roll to the front."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"red_image = np.zeros_like(color_image)\n",
"green_image = np.zeros_like(color_image)\n",
"blue_image = np.zeros_like(color_image)\n",
"\n",
"r, g, b = np.rollaxis(color_image, axis=-1)\n",
"red_image[:, :, 0] = r\n",
"green_image[:, :, 1] = g\n",
"blue_image[:, :, 2] = b\n",
"\n",
"skdemo.imshow_all(red_image, green_image, blue_image)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Histograms"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Question:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The following attempt to plot as histogram does not work:\n",
"\n",
"```python\n",
"plt.hist(color_image);\n",
"```\n",
"\n",
"The code above raises the following error:\n",
"\n",
"> ValueError: x must be 1D or 2D\n",
"\n",
"As the error suggests, a 3D, color image will not work with a histogram"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"\n",
"plt.hist(color_image.flatten());"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"These bins are really coarse, so lets use the `bins` argument to `plt.hist` to get a better result:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"plt.hist(color_image.flatten(), bins=100);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Why does that look so strange? (Think about it before continuing.)\n",
"\n",
"The pixel values are integers, and we're dividing up the data into 100 bins. If the number of integers isn't a divisor of the range of integers, then some bins will have more integers than others. Let's take a look at the number of unique bins:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"np.ptp(color_image)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"`np.ptp` gives the difference between the min and max values. If we used 115 bins, every bin would have approximately 2 integers, and we should get a good-looking histogram. If we use a bit fewer bins, then some bins will have 3 integers instead of 2; if go a bit higher, then some bins will have a bit fewer integers."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"fig, axes = plt.subplots(ncols=3, figsize=(18, 5))\n",
"\n",
"for ax, bins in zip(axes, (100, 115, 130)):\n",
" ax.hist(color_image.flatten(), bins=bins);\n",
" ax.set_ylim(ymax=8000)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Exercise:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Create a function to **plot a histogram for each color channel** in a single plot."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"for color, channel in zip('rgb', np.rollaxis(color_image, axis=-1)):\n",
" plt.hist(channel.flatten(), color=color, alpha=0.5)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Bonus: Instead of using `plt.hist`, use `plt.plot` or `plt.fill_between` in combination with `histogram` from `skimage.exposure`."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"from skimage import exposure\n",
"exposure.histogram?"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"for color, channel in zip('rgb', np.rollaxis(color_image, axis=-1)):\n",
" hist, bin_centers = exposure.histogram(channel)\n",
" plt.fill_between(bin_centers, hist, color=color, alpha=0.3)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Color spaces"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Exercise:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Use the CIELAB color space to **isolate the eyes** in the `chelsea` image. Plot the L, a, b components to get a feel for what they do, and see the [wikipedia page](http://en.wikipedia.org/wiki/Lab_colorspace) for more info. Bonus: **Isolate the nose** too."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"from skimage import color\n",
"\n",
"lab_image = color.rgb2lab(color_image)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"luminance, a, b = np.rollaxis(lab_image, axis=-1)\n",
"titles = ['luminance', 'a-component', 'b-component']\n",
"skdemo.imshow_all(luminance, a, b, titles=titles)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Notice the eyes are really dark in the middle image, i.e. the \"A\"-component of the LAB color space. We should use that for thresholding. We don't want to just guess values for thresholding; instead, let's look at a histogram:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"skdemo.imshow_with_histogram(lab_image);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"I've cheated here and reused our RGB histogram for LAB space. The image is wrong because the image isn't RGB, but the histograms are still valid. Note that we're plotting in RGB order, so the \"A\" commponent here is green. Since the eyes are very dark in the \"A\" image above, we'd like to select the left tail of the green histogram---roughly at 0:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"plt.imshow(a < 0)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To get the nose, notice that it's quite bright in the \"A\" component, so lets grab the right tail of the green histogram (> 30)."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"plt.imshow(a > 30)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Combining those two, plus an extra bit of thresholding to blacken the pupil, gives:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"eyes = (a < 0) & (b > 5)\n",
"nose = a > 30\n",
"plt.imshow(eyes | nose)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Note that we want to pixels that are either the nose or the eyes; you might be tempted to write `eyes & nose`, but that would only keep pixels that are both part of the eyes and the nose, which is an empty set."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Image filtering"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Replicate image from demo\n",
"image = data.camera()\n",
"pixelated = image[::10, ::10]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In the image filtering section, our first attempt at applying a difference filter didn't produce good results"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"horizontal_edges = pixelated[1:, :] - pixelated[:-1, :]\n",
"vertical_edges = pixelated[:, 1:] - pixelated[:, :-1]\n",
"skdemo.imshow_all(horizontal_edges, vertical_edges)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The reason for the noise is that these images are unsigned, 8-bit integers:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"print pixelated.dtype, horizontal_edges.dtype, vertical_edges.dtype"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"If you subtract two values, and the result is negative, then that value wraps around (underflow); thus, -1 becomes 255. Similarly, values exceeding 255, will wrap around (overflow) so that 256 becomes -1.\n",
"\n",
"To fix this, we'll convert to floating-point images:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"from skimage import img_as_float\n",
"pixelated_float = img_as_float(pixelated)\n",
"horizontal_edges = pixelated_float[1:, :] - pixelated_float[:-1, :]\n",
"vertical_edges = pixelated_float[:, 1:] - pixelated_float[:, :-1]\n",
"skdemo.imshow_all(horizontal_edges, vertical_edges)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Exercise:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Create a simple difference filter to **find the horizontal or vertical edges** of an image. Try to ensure that the filtering operation doesn't shift the edge position preferentially. (Don't use slicing to produce the difference image; use convolution.)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"def imshow_edges(image, horizontal_edge_kernel, vertical_edge_kernel):\n",
" horizontal_edges = convolve(image, horizontal_edge_kernel)\n",
" vertical_edges = convolve(image, vertical_edge_kernel)\n",
" skdemo.imshow_all(horizontal_edges, vertical_edges)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can emulate the simple differencing operation that we did earlier with the following edge kernels:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"from skimage import img_as_float\n",
"from scipy.ndimage import convolve\n",
"\n",
"horizontal_edge_kernel = np.array([[1], [-1]])\n",
"vertical_edge_kernel = np.array([[1, -1]])\n",
"\n",
"# As discussed earlier, use a floating-point image to prevent overflow\n",
"image = img_as_float(pixelated)\n",
"imshow_edges(image, horizontal_edge_kernel, vertical_edge_kernel)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"But, as we said earlier, this tends to skew the edges toward one corner of the image. To preserve the position of edge points, we can create a kernel that's centered on each pixel:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"horizontal_edge_kernel = np.array([[1], [0], [-1]])\n",
"vertical_edge_kernel = np.array([[1, 0, -1]])\n",
"\n",
"imshow_edges(image, horizontal_edge_kernel, vertical_edge_kernel)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can verify that this doesn't skew the edges by looking at the bright square from earlier"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Replicate image from demo\n",
"bright_square = np.zeros((7, 7), dtype=float)\n",
"bright_square[2:5, 2:5] = 1\n",
"\n",
"image = bright_square\n",
"horizontal_edges = convolve(image, horizontal_edge_kernel)\n",
"vertical_edges = convolve(image, vertical_edge_kernel)\n",
"skdemo.imshow_all(horizontal_edges, vertical_edges)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Exercise:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Using a couple of the filters in the `filter` module, **find the direction of the maximum gradient** in an image."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**NOTE:** When developing this exercise, it was assumed that `vsobel` and `hsobel` from `skimage.filter` returned the response from the vertical-edge and horizontal-edge Sobel kernels, pictured in the links below:\n",
"\n",
"* http://scikit-image.org/docs/dev/api/skimage.filter.html#vsobel\n",
"* http://scikit-image.org/docs/dev/api/skimage.filter.html#hsobel\n",
"\n",
"As described in those docs, however, the **absolute value** of the response is returned. While this result is typically preferred, it's not in this case, since the sign of the response contributes to the gradient angle.\n",
"\n",
"To get around this oversight, we'll copy the edge kernels from the documentation:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"from skimage import filter\n",
"\n",
"# Kernel copied from `vsobel` docstring.\n",
"# Vertical edge-reponse is the *horizontal* gradient.\n",
"dx_kernel = np.array([\n",
" [1, 0, -1],\n",
" [2, 0, -2],\n",
" [1, 0, -1],\n",
"])\n",
"# Rotate array by 90 degrees to get y-gradient kernel\n",
"dy_kernel = np.rot90(dx_kernel)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"image_45 = np.tril(-np.ones([7, 7]))\n",
"image_135 = np.rot90(image_45)\n",
"skdemo.imshow_all(image_45, image_135)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"def angle_image(image):\n",
" image = img_as_float(image)\n",
" dy = convolve(image, dy_kernel)\n",
" dx = convolve(image, dx_kernel)\n",
"\n",
" angle = np.arctan2(dy, dx)\n",
" return np.degrees(angle)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"%precision 0\n",
"angle_image(image_45)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"angle_image(image_135)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Feature detection"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Exercise:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Use the response from the Hough transform to **detect the position and radius of each coin**."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"First, let's replicate the data from the tutorial"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"from skimage.transform import hough_circle\n",
"\n",
"image = data.coins()[0:95, 180:370]\n",
"edges = filter.canny(image, sigma=3, low_threshold=10, high_threshold=60)\n",
"hough_radii = np.arange(15, 30, 2)\n",
"hough_response = hough_circle(edges, hough_radii)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"from skimage.feature import peak_local_max\n",
"\n",
"centers = []\n",
"likelihood = []\n",
"\n",
"for h in hough_response:\n",
" peaks = peak_local_max(h)\n",
" centers.extend(peaks)\n",
" likelihood.extend(h[peaks[:, 0], peaks[:, 1]])\n",
"\n",
"for i in np.argsort(likelihood)[-3:]:\n",
" row, column = centers[i]\n",
" plt.plot(column, row, 'ro')\n",
" plt.imshow(image)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"coin_centers"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"a = np.arange(10)\n",
"a[np.arange(3)]"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"centers"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"from skimage.feature import peak_local_max\n",
"\n",
"radii = []\n",
"centers = []\n",
"likelihood = []\n",
"for h, r in zip(hough_response, hough_radii):\n",
" peaks = peak_local_max(h)\n",
" radii.extend([r] * len(peaks))\n",
" centers.extend(peaks)\n",
" likelihood.extend(h[peaks[:, 0], peaks[:, 1]])"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"from skimage.draw import circle_perimeter\n",
"from skimage.color import gray2rgb\n",
"\n",
"drawn_image = gray2rgb(image)\n",
"plt.imshow(drawn_image)\n",
"for i in np.argsort(likelihood)[-3:]:\n",
" row, column = centers[i]\n",
" rr, cc = circle_perimeter(row, column, radii[i])\n",
" drawn_image[rr, cc, 1:] = 0\n",
" plt.plot(column, row, 'ro')"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"circle_perimeter?"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.4.3"
}
},
"nbformat": 4,
"nbformat_minor": 0
}