{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Networks of Image Regions with ndimage\n", "\n", "> Tyger Tyger, burning bright,
\n", "> In the forests of the night;
\n", "> What immortal hand or eye,
\n", "> Could frame thy fearful symmetry?\n", ">\n", "> — William Blake, *The Tyger*\n", "\n", "You probably know that digital images are made up of *pixels*. Generally, you\n", "should not think of these as little squares, but as *point samples* of the light\n", "signal *measured on a regular grid* [^alvyraysmith].\n", "\n", "Further, when processing\n", "images, we often deal with objects much larger than individual pixels.\n", "In a landscape, the sky, earth, trees, and rocks each span many\n", "pixels. A common structure to represent these is the Region Adjacency Graph,\n", "or RAG. Its *nodes* hold properties of each region in the image, and its\n", "*links* hold the spatial relationships between the regions. Two nodes are\n", "linked whenever their corresponding regions touch each other in the input\n", "image.\n", "\n", "Building such a structure could be a complicated\n", "affair, and even more difficult\n", "when images are not two-dimensional but 3D and even 4D, as is\n", "common in microscopy, materials science, and climatology, among others. But\n", "here we will show you how to produce a RAG in a few lines of code using\n", "NetworkX (a Python library to analyze graphs and networks), and\n", "a filter from SciPy's N-dimensional image processing submodule, `ndimage`.\n", "\n", "> **The origins of Elegant SciPy {.callout}**\n", ">\n", "> *(A note from Juan.)*\n", ">\n", "> This chapter gets a special mention because it inspired the whole book.\n", "> Vighnesh Birodkar wrote this code snippet as an undergraduate while\n", "> participating in Google Summer of Code (GSoC) 2014. When I saw this bit of\n", "> code, it blew me away. For the purposes of this book, it touches on many\n", "> aspects of scientific Python. By the time you're done with this chapter, you\n", "> should be able to process arrays of *any* dimension, rather than thinking of\n", "> them only as 1D lists or 2D tables. More than that, you'll understand the\n", "> basics of image filtering and network processing." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import networkx as nx\n", "import numpy as np\n", "from scipy import ndimage as ndi\n", "\n", "def add_edge_filter(values, graph):\n", " center = values[len(values) // 2]\n", " for neighbor in values:\n", " if neighbor != center and not graph.has_edge(center, neighbor):\n", " graph.add_edge(center, neighbor)\n", " return 0.0\n", "\n", "\n", "def build_rag(labels, image):\n", " g = nx.Graph()\n", " footprint = ndi.generate_binary_structure(labels.ndim, connectivity=1)\n", " _ = ndi.generic_filter(labels, add_edge_filter, footprint=footprint,\n", " mode='nearest', extra_arguments=(g,))\n", " return g" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "There are a few things going on here: images being represented as numpy arrays,\n", "*filtering* of these images using `scipy.ndimage`, and building of the image\n", "regions into a graph (network) using the NetworkX library. We'll go over these\n", "in turn.\n", "\n", "## Images are just numpy arrays\n", "\n", "In the previous chapter, we saw that numpy arrays can efficiently represent\n", "tabular data, and are a convenient way to perform computations on it.\n", "It turns out that arrays are equally adept at representing images.\n", "\n", "Here's how to create an image of white noise using just numpy, and display it\n", "with matplotlib. First, we import the necessary packages, and use the\n", "`matplotlib inline` IPython magic to make our images appear below the code:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Make plots appear inline, set custom plotting style\n", "%matplotlib inline\n", "import matplotlib.pyplot as plt\n", "plt.style.use('style/elegant.mplstyle')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Then, we \"make some noise\" and display it as an image:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "random_image = np.random.rand(500, 500)\n", "plt.imshow(random_image);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "This `imshow` function displays a numpy array as an image. The converse is also true: an image\n", "can be considered as a numpy array. For this example we use the scikit-image\n", "library, a collection of image processing tools built on top of NumPy and SciPy.\n", "\n", "Here is a PNG image from the scikit-image repository. It is a black and white\n", "(sometimes called \"grayscale\") picture of some ancient Roman coins from\n", "Pompeii, obtained from the Brooklyn Museum [^coins-source]:\n", "\n", "![Coins image from the Brooklyn Museum](https://raw.githubusercontent.com/scikit-image/scikit-image/v0.10.1/skimage/data/coins.png)\n", "\n", "Here is the coin image loaded with scikit-image:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from skimage import io\n", "url_coins = ('https://raw.githubusercontent.com/scikit-image/scikit-image/'\n", " 'v0.10.1/skimage/data/coins.png')\n", "coins = io.imread(url_coins)\n", "print(\"Type:\", type(coins), \"Shape:\", coins.shape, \"Data type:\", coins.dtype)\n", "plt.imshow(coins);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "A grayscale image can be represented as a *2-dimensional array*, with each array\n", "element containing the grayscale intensity at that position. So, **an image is\n", "just a numpy array**.\n", "\n", "Color images are a *3-dimensional* array, where the first two dimensions\n", "represent the spatial positions of the image, while the final dimension represents\n", "color channels, typically the three primary additive colors of red, green, and blue.\n", "To show what we can do with these dimensions, let's play with this photo of astronaut Eileen Collins:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "url_astronaut = ('https://raw.githubusercontent.com/scikit-image/scikit-image/'\n", " 'master/skimage/data/astronaut.png')\n", "astro = io.imread(url_astronaut)\n", "print(\"Type:\", type(astro), \"Shape:\", astro.shape, \"Data type:\", astro.dtype)\n", "plt.imshow(astro);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "This image is *just numpy arrays*. Adding a green square to the image is easy\n", "once you realize this, using simple numpy slicing:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "astro_sq = np.copy(astro)\n", "astro_sq[50:100, 50:100] = [0, 255, 0] # red, green, blue\n", "plt.imshow(astro_sq);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "You can also use a boolean *mask*, an array of `True` or `False` values.\n", "We saw these in Chapter 2 as a way to select rows of a table. In this case, we\n", "can use an array of the same shape as the image to select pixels:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "astro_sq = np.copy(astro)\n", "sq_mask = np.zeros(astro.shape[:2], bool)\n", "sq_mask[50:100, 50:100] = True\n", "astro_sq[sq_mask] = [0, 255, 0]\n", "plt.imshow(astro_sq);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "\n", "\n", "**Exercise:** We just saw how to select a square and paint it green. Can you\n", "extend that to other shapes and colors? Create a function to draw a blue grid\n", "onto a color image, and apply it to the `astronaut` image of Eileen Collins\n", "(above). Your function should take\n", "two parameters: the input image, and the grid spacing.\n", "Use the following template to help you get started." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def overlay_grid(image, spacing=128):\n", " \"\"\"Return an image with a grid overlay, using the provided spacing.\n", "\n", " Parameters\n", " ----------\n", " image : array, shape (M, N, 3)\n", " The input image.\n", " spacing : int\n", " The spacing between the grid lines.\n", "\n", " Returns\n", " -------\n", " image_gridded : array, shape (M, N, 3)\n", " The original image with a blue grid superimposed.\n", " \"\"\"\n", " image_gridded = image.copy()\n", " pass # replace this line with your code...\n", " return image_gridded\n", "\n", "# plt.imshow(overlay_grid(astro, 128)); # uncomment this line to test your function" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "**Solution:** We can use numpy slicing to select the rows of the grid, set them\n", "to blue, and then select the columns, and set them to blue as well:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def overlay_grid(image, spacing=128):\n", " \"\"\"Return an image with a grid overlay, using the provided spacing.\n", "\n", " Parameters\n", " ----------\n", " image : array, shape (M, N, 3)\n", " The input image.\n", " spacing : int\n", " The spacing between the grid lines.\n", "\n", " Returns\n", " -------\n", " image_gridded : array, shape (M, N, 3)\n", " The original image with a blue grid superimposed.\n", " \"\"\"\n", " image_gridded = image.copy()\n", " image_gridded[spacing:-1:spacing, :] = [0, 0, 255]\n", " image_gridded[:, spacing:-1:spacing] = [0, 0, 255]\n", " return image_gridded\n", "\n", "plt.imshow(overlay_grid(astro, 128));" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "Note that we used `-1` to mean the last value of the axis, as is standard in\n", "Python indexing. You can omit this value, but the meaning is slightly\n", "different. Without it (i.e. `spacing::spacing`), you go all the way to the end\n", "of the array, including the final row/column. When you use it as the stop\n", "index, you prevent the final row from being selected. In the case of a grid\n", "overlay, this is probably the desired behavior.\n", "\n", "\n", "\n", "\n", "\n", "## Filters in signal processing\n", "\n", "Filtering is one of the most fundamental and common operations in image\n", "processing. You can filter an image to remove noise, to enhance features, or to\n", "detect edges between objects in the image.\n", "\n", "To understand filters, it's easiest to start with a 1D signal, instead of an image. For\n", "example, you might measure the light arriving at your end of a fiber-optic cable.\n", "If you *sample* the signal every millisecond (ms) for 100ms, you end up with an\n", "array of length 100. Suppose that after 30ms the light signal is turned on, and\n", "30ms later, it is switched off. You end up with a signal like this:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "sig = np.zeros(100, np.float) #\n", "sig[30:60] = 1 # signal = 1 during the period 30-60ms because light is observed\n", "fig, ax = plt.subplots()\n", "ax.plot(sig);\n", "ax.set_ylim(-0.1, 1.1);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "To find *when* the light is turned on, you can *delay* it by 1ms, then\n", "*subtract* the original from the delayed signal. This way, when the signal is\n", "unchanged from one millisecond to the next, the subtraction will give zero,\n", "but when the signal *increases*, you will get a positive signal.\n", "\n", "When the signal *decreases*, we will get a negative signal. If we are\n", "only interested in pinpointing the time when the light was turned on, we can\n", "*clip* the difference signal, so that any negative values are converted to 0." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "sigdelta = sig[1:] # sigdelta[0] equals sig[1], and so on\n", "sigdiff = sigdelta - sig[:-1]\n", "sigon = np.clip(sigdiff, 0, np.inf)\n", "fig, ax = plt.subplots()\n", "ax.plot(sigon)\n", "ax.set_ylim(-0.1, 1.1)\n", "print('Signal on at:', 1 + np.flatnonzero(sigon)[0], 'ms')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "(Here we have used NumPy's `flatnonzero` function to get the first index where\n", "the `sigon` array is not equal to 0.)\n", "\n", "It turns out that this can be accomplished by an signal processing operation\n", "called *convolution*. At every point of the signal, we compute the dot-product\n", "between the values surrounding that point and a *kernel* or *filter*, which is a\n", "predetermined vector of values. Depending on the kernel, then, the convolution\n", "shows a different feature of the signal.\n", "\n", "Now, think of what happens when the kernel is (1, 0, -1), the difference\n", "filter, for a signal `s`. At any position `i`, the convolution result is\n", "`1*s[i+1] + 0*s[i] - 1*s[i-1]`, that is, `s[i+1] - s[i-1]`.\n", "Thus, when the values adjacent to `s[i]` are identical, the convolution gives 0, but when\n", "`s[i+1] > s[i-1]` (the signal is increasing), it gives a positive value, and,\n", "conversely, when `s[i+1] < s[i-1]`, it gives a negative value. You can think\n", "of this as an estimate of the derivative of the input function.\n", "\n", "In general, the formula for convolution is:\n", "$s'(t) = \\sum_{j=t-\\tau}^{t}{s(j)f(t-j)}$\n", "where $s$ is the signal, $s'$ is the filtered signal, $f$ is the filter, and\n", "$\\tau$ is the length of the filter.\n", "\n", "In scipy, you can use the `scipy.ndimage.convolve` to work on this." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "diff = np.array([1, 0, -1])\n", "from scipy import ndimage as ndi\n", "dsig = ndi.convolve(sig, diff)\n", "plt.plot(dsig);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "Signals are usually *noisy* though, not perfect as above:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "np.random.seed(0)\n", "sig = sig + np.random.normal(0, 0.3, size=sig.shape)\n", "plt.plot(sig);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "The plain difference filter can amplify that noise:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.plot(ndi.convolve(sig, diff));" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "In such cases, you can add smoothing to the filter. The most common form of\n", "smoothing is *Gaussian* smoothing, which takes the weighted average of\n", "neighboring points in the signal using the\n", "[Gaussian function](https://en.wikipedia.org/wiki/Gaussian_function). We can\n", "write a function to make a Gaussian smoothing kernel as follows:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def gaussian_kernel(size, sigma):\n", " \"\"\"Make a 1D Gaussian kernel of the specified size and standard deviation.\n", "\n", " The size should be an odd number and at least ~6 times greater than sigma\n", " to ensure sufficient coverage.\n", " \"\"\"\n", " positions = np.arange(size) - size // 2\n", " kernel_raw = np.exp(-positions**2 / (2 * sigma**2))\n", " kernel_normalized = kernel_raw / np.sum(kernel_raw)\n", " return kernel_normalized" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A really nice feature feature of convolution is that it's *associative*,\n", "meaning if you want to find the derivative of the smoothed signal, you can\n", "equivalently convolve the signal with the smoothed difference filter! This can\n", "save a lot of computation time, because you can smooth just the filter, which\n", "is usually much smaller than the data." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "smooth_diff = ndi.convolve(gaussian_kernel(25, 3), diff)\n", "plt.plot(smooth_diff);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "This smoothed difference filter looks for an edge in the central position,\n", "but also for that difference to continue. This continuation happens in the case\n", "of a true\n", "edge, but not in \"spurious\" edges caused by noise. Check out the result:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "sdsig = ndi.convolve(sig, smooth_diff)\n", "plt.plot(sdsig);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "Although it still looks wobbly, the *signal-to-noise ratio* (SNR),\n", "is much greater in this version than when using the simple difference filter.\n", "\n", "> **A note about filtering** {.callout}\n", ">\n", "> This operation is called filtering because, in physical electrical circuits,\n", "> many of these operations are implemented by hardware that allows certain\n", "> kinds of current through, while blocking others; these hardware components\n", "> are called filters. For example, a common filter that removes high-frequency\n", "> voltage fluctuations from a current is called a *low-pass filter*.\n", "\n", "## Filtering images (2D filters)\n", "\n", "Now that you've seen filtering in 1D, we hope you'll find it straightforward to\n", "extend these concepts to 2D signals, such as images. Here's a 2D difference\n", "filter for finding the edges in the coins image:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "coins = coins.astype(float) / 255 # prevents overflow errors\n", "diff2d = np.array([[0, 1, 0], [1, 0, -1], [0, -1, 0]])\n", "coins_edges = ndi.convolve(coins, diff2d)\n", "io.imshow(coins_edges);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "The principle is the same as the 1D filter: at every point in the image, place the\n", "filter, compute the dot-product of the filter's values with the image values, and\n", "place the result at the same location in the output image. And, as with the 1D\n", "difference filter, when the filter is placed on a location with little variation, the\n", "dot-product cancels out to zero, whereas, placed on a location where the\n", "image brightness is changing, the values multiplied by 1 will be different from\n", "those multiplied by -1, and the filtered output will be a positive or negative\n", "value (depending on whether the image is brighter towards the bottom-right\n", "or top-left at that point).\n", "\n", "Just as with the 1D filter, you can get more sophisticated and smooth out\n", "noise right within the filter. The *Sobel* filter is designed to do just that.\n", "It comes in horizontal and vertical varieties, to find edges with that\n", "orientation in the data.\n", "Let's start with the horizontal filter first.\n", "To find a horizontal edge in a picture, you might try the following filter:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# column vector (vertical) to find horizontal edges\n", "hdiff = np.array([[1], [0], [-1]])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "However, as we saw with 1D filters, this will result in a noisy estimate of the\n", "edges in the image. But rather than using Gaussian smoothing, which can cause\n", "blurry edges, the Sobel filter uses the property that edges in images tend to\n", "be continuous: a picture of the ocean, for example, will contain a horizontal\n", "edge along an entire line, not just at specific points of the image. So the\n", "Sobel filter smooths the vertical filter horizontally: it looks for a strong\n", "edge at the central position that is corroborated by the adjacent positions:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "hsobel = np.array([[ 1, 2, 1],\n", " [ 0, 0, 0],\n", " [-1, -2, -1]])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The vertical Sobel filter is simply the transpose of the horizontal:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "vsobel = hsobel.T" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can then find the horizontal and vertical edges in the coins image:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Some custom x-axis labelling to make our plots easier to read\n", "def reduce_xaxis_labels(ax, factor):\n", " \"\"\"Show only every ith label to prevent crowding on x-axis\n", " e.g. factor = 2 would plot every second x-axis label,\n", " starting at the first.\n", "\n", " Parameters\n", " ----------\n", " ax : matplotlib plot axis to be adjusted\n", " factor : int, factor to reduce the number of x-axis labels by\n", " \"\"\"\n", " plt.setp(ax.xaxis.get_ticklabels(), visible=False)\n", " for label in ax.xaxis.get_ticklabels()[::factor]:\n", " label.set_visible(True)\n", "\n", "\n", "coins_h = ndi.convolve(coins, hsobel)\n", "coins_v = ndi.convolve(coins, vsobel)\n", "\n", "fig, axes = plt.subplots(nrows=1, ncols=2)\n", "axes[0].imshow(coins_h, cmap=plt.cm.RdBu)\n", "axes[1].imshow(coins_v, cmap=plt.cm.RdBu)\n", "for ax in axes:\n", " reduce_xaxis_labels(ax, 2)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "And finally, just like the Pythagorean theorem, you can argue that the edge\n", "magnitude in *any* direction is equal to the square root of the sum of squares\n", "of the horizontal and vertical components:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "coins_sobel = np.sqrt(coins_h**2 + coins_v**2)\n", "plt.imshow(coins_sobel, cmap='viridis');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "## Generic filters: arbitrary functions of neighborhood values\n", "\n", "In addition to dot-products, implemented by `ndi.convolve`, SciPy lets you\n", "define a filter that is an *arbitrary function* of the points in a neighborhood,\n", "implemented in `ndi.generic_filter`. This can let you express arbitrarily\n", "complex filters.\n", "\n", "For example, suppose an image represents median house values in a county,\n", "with a 100m x 100m resolution. The local council decides to tax house sales as\n", "$10,000 plus 5% of the 90th percentile of house prices in a 1km radius. (So,\n", "selling a house in an expensive neighborhood costs more.) With\n", "`generic_filter`, we can produce the map of the tax rate everywhere in the map:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from skimage import morphology\n", "def tax(prices):\n", " return 10000 + 0.05 * np.percentile(prices, 90)\n", "house_price_map = (0.5 + np.random.rand(100, 100)) * 1e6\n", "footprint = morphology.disk(radius=10)\n", "tax_rate_map = ndi.generic_filter(house_price_map, tax, footprint=footprint)\n", "plt.imshow(tax_rate_map)\n", "plt.colorbar();" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "\n", "\n", "**Exercise:** Conway's Game of Life.\n", "\n", "Suggested by Nicolas Rougier.\n", "\n", "Conway's [Game of Life](https://en.wikipedia.org/wiki/Conway%27s_Game_of_Life) is a\n", "seemingly simple construct in which \"cells\" on a regular square grid live or die\n", "according to the cells in their immediate surroundings. At every timestep, we\n", "determine the state of position (i, j) according to its previous state and that\n", "of its 8 neighbors (above, below, left, right, and diagonals):\n", "\n", "- a live cell with only one live neighbor or none dies.\n", "- a live cell with two or three live neighbors lives on for another generation.\n", "- a live cell with four or more live neighbors dies, as if from overpopulation.\n", "- a dead cell with exactly three live neighbors becomes alive, as if by\n", " reproduction.\n", "\n", "Although the rules sound like a contrived math problem, they in fact give rise\n", "to incredible patterns, starting with gliders (small patterns of live cells\n", "that slowly move in each generation) and glider guns (stationary patterns that\n", "sprout off gliders), all the way up to prime number generator machines (see,\n", "for example,\n", "[this page](http://www.njohnston.ca/2009/08/generating-sequences-of-primes-in-conways-game-of-life/)),\n", "and even\n", "[simulating Game of Life itself](https://www.youtube.com/watch?v=xP5-iIeKXE8)!\n", "\n", "Can you implement the Game of Life using `ndi.generic_filter`?\n", "\n", "\n", "\n", "**Solution:**\n", "\n", "Nicolas Rougier (@rougier) provides a NumPy-only solution on his 100 NumPy\n", "Exercises page (Exercise #79):" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def next_generation(Z):\n", " N = (Z[0:-2,0:-2] + Z[0:-2,1:-1] + Z[0:-2,2:] +\n", " Z[1:-1,0:-2] + Z[1:-1,2:] +\n", " Z[2: ,0:-2] + Z[2: ,1:-1] + Z[2: ,2:])\n", "\n", " # Apply rules\n", " birth = (N==3) & (Z[1:-1,1:-1]==0)\n", " survive = ((N==2) | (N==3)) & (Z[1:-1,1:-1]==1)\n", " Z[...] = 0\n", " Z[1:-1,1:-1][birth | survive] = 1\n", " return Z" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Then we can start a board with:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "random_board = np.random.randint(0, 2, size=(50, 50))\n", "n_generations = 100\n", "for generation in range(n_generations):\n", " random_board = next_generation(random_board)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Using generic filter makes it even easier:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def nextgen_filter(values):\n", " center = values[len(values) // 2]\n", " neighbors_count = np.sum(values) - center\n", " if neighbors_count == 3 or (center and neighbors_count == 2):\n", " return 1.\n", " else:\n", " return 0.\n", "\n", "def next_generation(board):\n", " return ndi.generic_filter(board, nextgen_filter,\n", " size=3, mode='constant')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The nice thing is that some formulations of the Game of Life use what's known\n", "as a *toroidal board*, which means that the left and right ends \"wrap around\"\n", "and connect to each other, as well as the top and bottom ends. With\n", "`generic_filter`, it's trivial to modify our solution to incorporate this:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def next_generation_toroidal(board):\n", " return ndi.generic_filter(board, nextgen_filter,\n", " size=3, mode='wrap')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can now simulate this toroidal board for a few generations:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "random_board = np.random.randint(0, 2, size=(50, 50))\n", "n_generations = 100\n", "for generation in range(n_generations):\n", " random_board = next_generation_toroidal(random_board)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "\n", "\n", "\n", "\n", "**Exercise:** Sobel gradient magnitude.\n", "\n", "Above, we saw how we can combine the output of two different filters, the\n", "horizontal Sobel filter, and the vertical one. Can you write a function that\n", "does this in a single pass using `ndi.generic_filter`?\n", "\n", "" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "hsobel = np.array([[ 1, 2, 1],\n", " [ 0, 0, 0],\n", " [-1, -2, -1]])\n", "\n", "vsobel = hsobel.T\n", "\n", "hsobel_r = np.ravel(hsobel)\n", "vsobel_r = np.ravel(vsobel)\n", "\n", "def sobel_magnitude_filter(values):\n", " h_edge = values @ hsobel_r\n", " v_edge = values @ vsobel_r\n", " return np.hypot(h_edge, v_edge)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we can try it out on the coins image:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "sobel_mag = ndi.generic_filter(coins, sobel_magnitude_filter, size=3)\n", "plt.imshow(sobel_mag, cmap='viridis');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "\n", "\n", "\n", "\n", "## Graphs and the NetworkX library\n", "\n", "Graphs are a natural representation for an astonishing variety of data. Pages\n", "on the world wide web, for example, can comprise nodes, while links between\n", "those pages can be, well, links. Or, in biology, so-called *transcription\n", "networks* have nodes represent genes and edges connect genes that have a direct\n", "influence on each other's expression.\n", "\n", "> **Note: graphs and networks {.callout}**\n", ">\n", "> In this context, the term \"graph\" is synonymous with \"network\", not with\n", "> \"plot\". Mathematicians and computer scientists invented slightly different\n", "> words to discuss these: graph = network, vertex = node, edge = link = arc. As\n", "> most people do, we will be using these terms interchangeably.\n", ">\n", "> You might be slightly more familiar with the network terminology: a network\n", "> consists of *nodes* and *links* between the nodes. Equivalently, a graph\n", "> consists of *vertices* and *edges* between the vertices. In NetworkX, you\n", "> have `Graph` objects consisting of `nodes` and `edges` between the nodes, and\n", "> this is probably the most common usage.\n", "\n", "To introduce you to graphs, we will reproduce some results from the paper\n", "[\"Structural properties of the *Caenorhabditis elegans* neuronal network\"](http://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1001066), by Lav Varshney *et al*, 2011.\n", "\n", "In our example, we will represent neurons in the nematode worm's nervous system as\n", "nodes, and place an edge between two nodes when a neuron makes a synapse with\n", "another. (*Synapses* are the chemical connections through which neurons\n", "communicate.) The worm is an awesome example of neural connectivity analysis\n", "because every worm (of this species) has the same number of neurons (302), and the\n", "connections between them are all known. This has resulted in the fantastic Openworm\n", "project [^openworm], which we encourage you to read more about.\n", "\n", "You can download the neuronal dataset in Excel format from the WormAtlas\n", "database at [http://www.wormatlas.org/neuronalwiring.html#Connectivitydata](http://www.wormatlas.org/neuronalwiring.html#Connectivitydata).\n", "The `pandas` library allows one to read an Excel table over the web, so we will\n", "use it here to read in the data, then feed that into NetworkX." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import pandas as pd\n", "connectome_url = 'http://www.wormatlas.org/images/NeuronConnect.xls'\n", "conn = pd.read_excel(connectome_url)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "`conn` now contains a pandas DataFrame, with rows of the form:\n", "\n", "[Neuron1, Neuron2, connection type, strength]\n", "\n", "We are only going to examine the connectome of chemical synapses, so we filter\n", "out other synapse types as follows:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "conn_edges = [(n1, n2, {'weight': s})\n", " for n1, n2, t, s in conn.itertuples(index=False, name=None)\n", " if t.startswith('S')]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "(Look at the WormAtlas page for a description of the different connection types.)\n", "\n", "We use `weight` in a dictionary above because it is a special keyword for\n", "edge properties in NetworkX. We then build the graph using NetworkX's\n", "`DiGraph` class:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import networkx as nx\n", "wormbrain = nx.DiGraph()\n", "wormbrain.add_edges_from(conn_edges)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can now examine some of the properties of this network. One of the\n", "first things researchers ask about directed networks is which nodes are\n", "the most critical to information flow within it. Nodes with high\n", "*betweenness centrality* are those that belong to the shortest path between\n", "many different pairs of nodes. Think of a rail network: certain stations will\n", "connect to many lines, so that you will be forced to change lines there\n", "for many different trips. They are the ones with high betweenness\n", "centrality.\n", "\n", "With NetworkX, we can find similarly important neurons with ease. In the\n", "NetworkX API documentation [^nxdoc], under \"centrality\", the docstring\n", "for `betweenness_centrality` [^bwcdoc] specifies a function that takes a\n", "graph as input and returns a dictionary mapping node IDs to betweenness\n", "centrality values (floating point values)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "centrality = nx.betweenness_centrality(wormbrain)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we can find the neurons with highest centrality using the Python built-in\n", "function `sorted`:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "central = sorted(centrality, key=centrality.get, reverse=True)\n", "print(central[:5])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This returns the neurons AVAR, AVAL, PVCR, PVT, and PVCL, which have been\n", "implicated in how the worm responds to prodding: the AVA neurons link\n", "the worm's front touch receptors (among others) to neurons responsible\n", "for backward motion, while the PVC neurons link the rear touch receptors to\n", "forward motion.\n", "\n", "Varshney *et al* study the properties of a *strongly connected component*\n", "of 237 neurons, out of a total of 279. In graphs, a\n", "*connected component* is a set of nodes that are reachable by some path\n", "through all the links. The connectome is a *directed* graph, meaning the\n", "edges *point* from one node to the other, rather than merely connecting\n", "them. In this case, a strongly connected component is one where all nodes\n", "are reachable from each other by traversing links *in the correct direction*.\n", "So A $\\rightarrow$ B $\\rightarrow$ C is not strongly connected, because there is no way to get to\n", "A from B or C. However, A $\\rightarrow$ B $\\rightarrow$ C $\\rightarrow$ A *is* strongly connected.\n", "\n", "In a neuronal circuit, you can think of the strongly connected component\n", "as the \"brain\" of the circuit, where the processing happens, while nodes\n", "upstream of it are inputs, and nodes downstream are outputs.\n", "\n", "> **Cycles in neuronal networks {.callout}**\n", ">\n", "> The idea of cyclical neuronal circuits dates back to the 1950s. Here's a\n", "> lovely paragraph about this idea from an article in *Nautilus*,\n", "> \"The Man Who Tried to Redeem the World With Logic\", by Amanda Gefter:\n", ">\n", "> > If one were to see a lightning bolt flash on the sky, the eyes would send a signal to the brain, shuffling it through a chain of neurons. Starting with any given neuron in the chain, you could retrace the signal's steps and figure out just how long ago the lightning struck. Unless, that is, the chain is a loop. In that case, the information encoding the lightning bolt just spins in circles, endlessly. It bears no connection to the time at which the lightning actually occurred. It becomes, as McCulloch put it, \"an idea wrenched out of time.\" In other words, a memory.\n", "\n", "NetworkX makes straightforward work out of getting the largest strongly\n", "connected component from our `wormbrain` network:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "sccs = (wormbrain.subgraph(c) for c in nx.strongly_connected_components(wormbrain))\n", "giantscc = max(sccs, key=len)\n", "print(f'The largest strongly connected component has '\n", " f'{giantscc.number_of_nodes()} nodes, out of '\n", " f'{wormbrain.number_of_nodes()} total.')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As noted in the paper, the size of this component is *smaller* than\n", "expected by chance, demonstrating that the network is segregated into\n", "input, central, and output layers.\n", "\n", "Now we reproduce figure 6B from the paper, the survival function of the\n", "in-degree distribution. First, compute the relevant quantities:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "in_degrees = list(dict(wormbrain.in_degree()).values())\n", "in_deg_distrib = np.bincount(in_degrees)\n", "avg_in_degree = np.mean(in_degrees)\n", "cumfreq = np.cumsum(in_deg_distrib) / np.sum(in_deg_distrib)\n", "survival = 1 - cumfreq" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Then, plot using Matplotlib:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots()\n", "ax.loglog(np.arange(1, len(survival) + 1), survival)\n", "ax.set_xlabel('in-degree distribution')\n", "ax.set_ylabel('fraction of neurons with higher in-degree distribution')\n", "ax.scatter(avg_in_degree, 0.0022, marker='v')\n", "ax.text(avg_in_degree - 0.5, 0.003, 'mean=%.2f' % avg_in_degree)\n", "ax.set_ylim(0.002, 1.0);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "There you have it: a reproduction of a scientific analysis, using SciPy. We are\n", "missing the line fit.... But that's what exercises are for.\n", "\n", "\n", "\n", "This exercise is a bit of a preview for chapter 7 (optimization):\n", "use `scipy.optimize.curve_fit` to fit the tail of the\n", "in-degree survival function to a power-law,\n", "$f(d) \\sim d^{-\\gamma}, d > d_0$,\n", "for $d_0 = 10$ (the red line in Figure 6B of the paper), and modify the plot\n", "to include that line.\n", "\n", "\n", "\n", "**Solution:** Let's look at the start of the docstring for `curve_fit`:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "```\n", "Use non-linear least squares to fit a function, f, to data.\n", "\n", "Assumes ``ydata = f(xdata, *params) + eps``\n", "\n", "Parameters\n", "----------\n", "f : callable\n", " The model function, f(x, ...). It must take the independent\n", " variable as the first argument and the parameters to fit as\n", " separate remaining arguments.\n", "xdata : An M-length sequence or an (k,M)-shaped array\n", " for functions with k predictors.\n", " The independent variable where the data is measured.\n", "ydata : M-length sequence\n", " The dependent data --- nominally f(xdata, ...)\n", "```\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It looks like we just need to provide a function that takes in a data point,\n", "and some parameters, and returns the predicted value. In our case, we want the\n", "cumulative remaining frequency, $f(d)$ to be proportional to $d^{-\\gamma}$.\n", "That means we need $f(d) = \\alpha d^{-gamma}$:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def fraction_higher(degree, alpha, gamma):\n", " return alpha * degree ** (-gamma)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Then, we need our x and y data to fit, *for $d > 10$*:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "x = 1 + np.arange(len(survival))\n", "valid = x > 10\n", "x = x[valid]\n", "y = survival[valid]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can now use `curve_fit` to obtain fit parameters:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from scipy.optimize import curve_fit\n", "\n", "alpha_fit, gamma_fit = curve_fit(fraction_higher, x, y)[0]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's plot the results to see how we did:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "y_fit = fraction_higher(x, alpha_fit, gamma_fit)\n", "\n", "fig, ax = plt.subplots()\n", "ax.loglog(np.arange(1, len(survival) + 1), survival)\n", "ax.set_xlabel('in-degree distribution')\n", "ax.set_ylabel('fraction of neurons with higher in-degree distribution')\n", "ax.scatter(avg_in_degree, 0.0022, marker='v')\n", "ax.text(avg_in_degree - 0.5, 0.003, 'mean=%.2f' % avg_in_degree)\n", "ax.set_ylim(0.002, 1.0)\n", "ax.loglog(x, y_fit, c='red');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "Voilà! A full Figure 6B, fit and all!\n", "\n", "\n", "\n", "\n", "\n", "You now should have a fundamental understanding of graphs as a scientific\n", "abstraction, and how to easily manipulate and analyse them using Python and\n", "NetworkX. Now, we move on to a particular kind of graph used in image\n", "processing and computer vision.\n", "\n", "## Region adjacency graphs\n", "\n", "A Region Adjacency Graph (RAG) is a representation of an image that is useful\n", "for *segmentation*: the division of images into meaningful regions (or\n", "*segments*). If you've seen Terminator 2, you've seen segmentation:\n", "\n", "![Terminator vision](../images/terminator-vision.png)\n", "\n", "Segmentation is one of those problems that humans do trivially, all the time,\n", "without thinking, whereas computers have a hard time of it. To\n", "understand this difficulty, look at this image:\n", "\n", "![Face (Eileen Collins)](http://i.imgur.com/ky5qwIS.png)\n", "\n", "While you see a face, a computer only sees a bunch of numbers:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "```\n", " 58688888888888899998898988888666532121\n", " 66888886888998999999899998888888865421\n", " 66665566566689999999999998888888888653\n", " 66668899998655688999899988888668665554\n", " 66888899998888888889988888665666666543\n", " 66888888886868868889998888666688888865\n", " 66666443334556688889988866666666668866\n", " 66884235221446588889988665644644444666\n", " 86864486233664666889886655464321242345\n", " 86666658333685588888866655659381366324\n", " 88866686688666866888886658588422485434\n", " 88888888888688688888866566686666565444\n", " 88888888868666888888866556688666686555\n", " 88888988888888888888886656888688886666\n", " 88889999989998888888886666888888868886\n", " 88889998888888888888886566888888888866\n", " 88888998888888688888666566868868888888\n", " 68888999888888888868886656888888888866\n", " 68888999998888888688888655688888888866\n", " 68888999886686668886888656566888888886\n", " 88888888886668888888888656558888888886\n", " 68888886665668888889888555555688888886\n", " 86868868658668868688886555555588886866\n", " 66688866468866855566655445555656888866\n", " 66688654888886868666555554556666666865\n", " 88688658688888888886666655556686688665\n", " 68888886666888888988888866666656686665\n", " 66888888845686888999888886666556866655\n", " 66688888862456668866666654431268686655\n", " 68688898886689696666655655313668688655\n", " 68888898888668998998998885356888986655\n", " 68688889888866899999999866666668986655\n", " 68888888888866666888866666666688866655\n", " 56888888888686889986868655566688886555\n", " 36668888888868888868688666686688866655\n", " 26686888888888888888888666688688865654\n", " 28688888888888888888668666668686666555\n", " 28666688888888888868668668688886665548\n", "```\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Our visual system is so optimized to spot faces that you might see the\n", "face even in this blob of numbers! But we hope our point is made. Also,\n", "you might want to look for the \"Faces In Things\" Tumblr, which demonstrates\n", "the face-finding optimization of our visual systems far more humorously.\n", "\n", "At any rate, the challenge is to make sense of those numbers, and where the\n", "boundaries lie that divide the different parts of the image. A popular\n", "approach is to find small regions (called superpixels) that\n", "you're *sure* belong in the same segment, and then merge those according\n", "to some more sophisticated rule.\n", "\n", "As a simple example, suppose you want to segment out the tiger in this\n", "picture, from the Berkeley Segmentation Dataset (BSDS):\n", "\n", "![BSDS-108073 tiger](http://www.eecs.berkeley.edu/Research/Projects/CS/vision/bsds/BSDS300/html/images/plain/normal/color/108073.jpg)\n", "\n", "A clustering algorithm, simple linear iterative clustering (SLIC) [^slic], can give\n", "us a decent starting point. It is available in the scikit-image library." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "url = ('http://www.eecs.berkeley.edu/Research/Projects/CS/vision/'\n", " 'bsds/BSDS300/html/images/plain/normal/color/108073.jpg')\n", "tiger = io.imread(url)\n", "from skimage import segmentation\n", "seg = segmentation.slic(tiger, n_segments=30, compactness=40.0,\n", " enforce_connectivity=True, sigma=3)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Scikit-image also has a function to *display* segmentations, which we use to\n", "visualize the result of SLIC:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from skimage import color\n", "io.imshow(color.label2rgb(seg, tiger));" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "This shows that the body of the tiger has been split in three parts, with the\n", "rest of the image in the remaining segments.\n", "\n", "A region adjacency graph (RAG) is a graph in which every node represents one\n", "of the above regions, and an edge connects two nodes when they touch. For a\n", "taste of what it looks like before we build one, we'll use the `show_rag` function\n", "from scikit-image — indeed, the library that contains this chapter's code snippet!" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from skimage.future import graph\n", "\n", "g = graph.rag_mean_color(tiger, seg)\n", "graph.show_rag(seg, g, tiger);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "Here, you can see the nodes corresponding to each segment, and the edges\n", "between adjacent segments. These are colored with the YlGnBu (yellow-green-blue)\n", "colormap from matplotlib, according to the difference in color between the\n", "two nodes.\n", "\n", "The figure also shows the magic of thinking of segmentations as graphs: you can\n", "see that edges between nodes within the tiger and those outside of it are brighter\n", "(higher-valued) than edges within the same object. Thus, if we can cut the\n", "graph along those edges, we will get our segmentation. We have chosen an easy\n", "example for color-based segmentation, but the same principles hold true for\n", "graphs with more complicated pairwise relationships.\n", "\n", "## Elegant ndimage: how to build graphs from image regions\n", "\n", "All the pieces are in place: you know about numpy arrays, image filtering,\n", "generic filters, graphs, and region adjacency graphs. Let's build one to pluck\n", "the tiger out of that picture!\n", "\n", "The obvious approach is to use two nested for-loops to iterate over every pixel\n", "of the image, look at the neighboring pixels, and checking for different labels:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import networkx as nx\n", "def build_rag(labels, image):\n", " g = nx.Graph()\n", " nrows, ncols = labels.shape\n", " for row in range(nrows):\n", " for col in range(ncols):\n", " current_label = labels[row, col]\n", " if not current_label in g:\n", " g.add_node(current_label)\n", " g.node[current_label]['total color'] = np.zeros(3, dtype=np.float)\n", " g.node[current_label]['pixel count'] = 0\n", " if row < nrows - 1 and labels[row + 1, col] != current_label:\n", " g.add_edge(current_label, labels[row + 1, col])\n", " if col < ncols - 1 and labels[row, col + 1] != current_label:\n", " g.add_edge(current_label, labels[row, col + 1])\n", " g.node[current_label]['total color'] += image[row, col]\n", " g.node[current_label]['pixel count'] += 1\n", " return g" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Whew! This works, but if you want to segment a 3D image, you'll have to write a\n", "different version:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import networkx as nx\n", "def build_rag_3d(labels, image):\n", " g = nx.Graph()\n", " nplns, nrows, ncols = labels.shape\n", " for pln in range(nplns):\n", " for row in range(nrows):\n", " for col in range(ncols):\n", " current_label = labels[pln, row, col]\n", " if not current_label in g:\n", " g.add_node(current_label)\n", " g.node[current_label]['total color'] = np.zeros(3, dtype=np.float)\n", " g.node[current_label]['pixel count'] = 0\n", " if pln < nplns - 1 and labels[pln + 1, row, col] != current_label:\n", " g.add_edge(current_label, labels[pln + 1, row, col])\n", " if row < nrows - 1 and labels[pln, row + 1, col] != current_label:\n", " g.add_edge(current_label, labels[pln, row + 1, col])\n", " if col < ncols - 1 and labels[pln, row, col + 1] != current_label:\n", " g.add_edge(current_label, labels[pln, row, col + 1])\n", " g.node[current_label]['total color'] += image[pln, row, col]\n", " g.node[current_label]['pixel count'] += 1\n", " return g" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Both of these are pretty ugly and unwieldy, too. And difficult to extend:\n", "if we want to count diagonally neighboring pixels as adjacent (that is,\n", "[row, col] is \"adjacent to\" [row + 1, col + 1]), the code becomes even\n", "messier. And if we want to analyze 3D video, we need yet another\n", "dimension, and another level of nesting. It's a mess!\n", "\n", "Enter Vighnesh's insight: SciPy's `generic_filter` function already does\n", "this iteration for us! We used it above to compute an arbitrarily\n", "complicated function on the neighborhood of every element of a numpy\n", "array. Only now we don't want a filtered image out of the function: we\n", "want a graph. It turns out that `generic_filter` lets you pass additional\n", "arguments to the filter function, and we can use that to build the graph:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import networkx as nx\n", "import numpy as np\n", "from scipy import ndimage as ndi\n", "\n", "def add_edge_filter(values, graph):\n", " center = values[len(values) // 2]\n", " for neighbor in values:\n", " if neighbor != center and not graph.has_edge(center, neighbor):\n", " graph.add_edge(center, neighbor)\n", " # float return value is unused but needed by `generic_filter`\n", " return 0.0\n", "\n", "def build_rag(labels, image):\n", " g = nx.Graph()\n", " footprint = ndi.generate_binary_structure(labels.ndim, connectivity=1)\n", " _ = ndi.generic_filter(labels, add_edge_filter, footprint=footprint,\n", " mode='nearest', extra_arguments=(g,))\n", " for n in g:\n", " g.nodes[n]['total color'] = np.zeros(3, np.double)\n", " g.nodes[n]['pixel count'] = 0\n", " for index in np.ndindex(labels.shape):\n", " n = labels[index]\n", " g.nodes[n]['total color'] += image[index]\n", " g.nodes[n]['pixel count'] += 1\n", " return g" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here are a few reasons this is a brilliant piece of code:\n", "\n", "- `ndi.generic_filter` iterates over array elements *with their neighbors*.\n", " (Use `numpy.ndindex` to simply iterate over array indices.)\n", "- We return \"0.0\" from the filter function because `generic_filter` requires\n", " the filter function to return a float. However, we ignore the filter\n", " output (which is zero everywhere), and use it only for its \"side effect\" of adding edges to the graph.\n", "- The loops are not nested several levels deep. This makes the code more\n", " compact, easier to take in in one go.\n", "- The code works identically for 1D, 2D, 3D, or even 8D images!\n", "- If we want to add support for diagonal connectivity, we just need to\n", " change the `connectivity` parameter to `ndi.generate_binary_structure`\n", "\n", "## Putting it all together: mean color segmentation\n", "\n", "Now, we can use everything we've learned to segment the tiger in the image above:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "g = build_rag(seg, tiger)\n", "for n in g:\n", " node = g.nodes[n]\n", " node['mean'] = node['total color'] / node['pixel count']\n", "for u, v in g.edges():\n", " d = g.nodes[u]['mean'] - g.nodes[v]['mean']\n", " g[u][v]['weight'] = np.linalg.norm(d)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Each edge holds the difference between the average color of each segment.\n", "We can now threshold the graph:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def threshold_graph(g, t):\n", " to_remove = [(u, v) for (u, v, d) in g.edges(data=True)\n", " if d['weight'] > t]\n", " g.remove_edges_from(to_remove)\n", "threshold_graph(g, 80)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally, we use the numpy index-with-an-array trick we learned in chapter 2:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "map_array = np.zeros(np.max(seg) + 1, int)\n", "for i, segment in enumerate(nx.connected_components(g)):\n", " for initial in segment:\n", " map_array[int(initial)] = i\n", "segmented = map_array[seg]\n", "plt.imshow(color.label2rgb(segmented, tiger));" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "Oops! Looks like the cat lost its tail!\n", "\n", "Still, we think that's a nice demonstration of the capabilities of RAGs...\n", "And the beauty with which SciPy and NetworkX make it feasible.\n", "Many of these functions are available in the scikit-image library. If you\n", "are interested in image analysis, look it up!\n", "\n", "[^alvyraysmith]: A Pixel Is Not A Little Square. Alvy Ray Smith, 1995, Technical\n", " Memo. http://alvyray.com/Memos/CG/Microsoft/6_pixel.pdf\n", "\n", "[^coins-source]: http://www.brooklynmuseum.org/opencollection/archives/image/15641/image\n", "\n", "[^openworm]: http://www.openworm.org\n", "\n", "[^file-url]: https://github.com/scikit-image/scikit-image/tree/master/skimage/io/util.py\n", "\n", "[^nxdoc]: http://networkx.github.io/documentation/latest/reference/index.html\n", "\n", "[^bwcdoc]: http://networkx.github.io/documentation/latest/reference/generated/networkx.algorithms.centrality.betweenness_centrality.html\n", "\n", "[^bsdstiger]: http://www.eecs.berkeley.edu/Research/Projects/CS/vision/bsds/BSDS300/html/dataset/images/color/108073.html\n", "\n", "[^slic]: http://ivrg.epfl.ch/research/superpixels\n", "\n", "" ] } ], "metadata": {}, "nbformat": 4, "nbformat_minor": 4 }