{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "from __future__ import division, print_function\n", "import numpy as np\n", "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Overview\n", "\n", "- http://scikit-image.org/docs/stable/api/skimage.transform.html\n", "- http://scikit-image.org/docs/stable/api/skimage.transform.html#skimage.transform.warp\n", "- http://scikit-image.org/docs/stable/api/skimage.transform.html#skimage.transform.AffineTransform (and other similar classes)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Image rotation from scratch\n", "\n", "The following code shows how to rotate an image using the skimage (scikit-image) library." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "from skimage import transform, data\n", "\n", "camera = data.camera()\n", "rotated = transform.rotate(camera, 30)\n", "\n", "f, (ax0, ax1) = plt.subplots(1, 2, figsize=(10, 5))\n", "ax0.imshow(camera, cmap='gray')\n", "ax1.imshow(rotated, cmap='gray');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Exercise:** Write an algorithm from scratch that will\n", "do the same (i.e., take an input image as an ndarray, and rotate it).\n", "\n", "If you feel creative, you can also write code to magnify (zoom) the image.\n", "

\n", "You may need: http://en.wikipedia.org/wiki/Polar_coordinate_system\n", "

\n", "A (bad) solution is given below--but try it yourself before looking!" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### A problematic approach" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from skimage import color" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def rotate(image, theta):\n", " theta = np.deg2rad(theta)\n", " \n", " height, width = image.shape[:2]\n", " out = np.zeros_like(image)\n", " \n", " centre_x, centre_y = width / 2., height / 2.\n", " \n", " for x in range(width):\n", " for y in range(height):\n", " \n", " x_c = x - centre_x\n", " y_c = y - centre_y\n", " \n", " # Determine polar coordinate of pixel\n", " radius = np.sqrt(x_c**2 + y_c**2)\n", " angle = np.arctan2(y_c, x_c)\n", " \n", " new_angle = angle + theta\n", " \n", " new_x = radius * np.cos(new_angle)\n", " new_y = radius * np.sin(new_angle)\n", " \n", " new_x = new_x + centre_x\n", " new_y = new_y + centre_y\n", " \n", " if (new_x >= width) or (new_x < 0) or\\\n", " (new_y >= height) or (new_y < 0):\n", " continue\n", " else:\n", " out[int(new_y), int(new_x)] = image[y, x]\n", " \n", " return out\n", "\n", "rotated = rotate(camera, 40)\n", " \n", "plt.imshow(rotated, cmap='gray', interpolation='nearest');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## And while we can attempt to fix the problem...\n", "\n", "...this is not an optimal approach" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Attempt at fixing the holes using a median filter\n", "# -- it works, sort of, but it's not the best approach.\n", "\n", "height, width = rotated.shape[:2]\n", "\n", "out = rotated.copy()\n", "\n", "for x in range(1, width - 1):\n", " for y in range(1, height - 1):\n", " if out[y, x] == 0:\n", " out[y, x] = np.median([out[y, x-1],\n", " out[y, x+1],\n", " out[y+1, x],\n", " out[y-1, x]])\n", " \n", "plt.imshow(out, cmap='gray', interpolation='nearest');" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "A = np.array([[4, 2], [1, 6]])\n", "print(A)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "plt.imshow(A, cmap='gray', interpolation='nearest');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# For later discussion: interpolation\n", "\n", "## Bi-linear interpolation\n", "\n", "\n", "
\n", "\n", "Also see [bilinear interpolation on Wikipedia](http://en.wikipedia.org/wiki/Bilinear_interpolation)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Some warping experiments!\n", "\n", "## Fish-eye" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from skimage import transform, data, io\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", "# Load face\n", "face = io.imread('../images/stefan.jpg')\n", "\n", "# Get the eye nicely in the middle\n", "face = face[:185, 15:]" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "plt.imshow(face)\n", "plt.plot([face.shape[1]/2.], [face.shape[0]/2.], 'or', markersize=14, alpha=0.4)\n", "plt.axis('image');" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Define a transformation on the x-y coordinates\n", "\n", "def fisheye(xy):\n", " center = np.mean(xy, axis=0)\n", " xc, yc = (xy - center).T\n", "\n", " # Polar coordinates\n", " r = np.sqrt(xc**2 + yc**2)\n", " theta = np.arctan2(yc, xc)\n", "\n", " r = 0.8 * np.exp(r**(1/2.1) / 1.8)\n", "\n", " return np.column_stack((r * np.cos(theta), r * np.sin(theta))) + center" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Warp and display\n", "\n", "out = transform.warp(face, fisheye)\n", "\n", "f, (ax0, ax1) = plt.subplots(1, 2, figsize=(10, 5))\n", "ax0.imshow(face)\n", "ax0.set_axis_off()\n", "\n", "ax1.imshow(out)\n", "ax1.set_axis_off()\n", "\n", "plt.title('Knock! Knock!')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Run the following scripts for fun:\n", "\n", "(Open up the terminal in the \"scripts\" directory first)\n", "\n", "- **deswirl.py** (run using: ``python deswirl.py``)\n", "\n", " In the UK, a criminal tried to hide his identity by posting\n", " swirled pictures of his face online. Here, we use the\n", " Mona Lisa to illustrate what he did. Can you restore\n", " her face back to normal? (Note that you can adjust the\n", " position of the red dot, as well as move the sliders.)\n", " \n", " \n", "- **clock_deblur.py**\n", "\n", " I took a picture of a wall clock while moving the camera. Or perhaps the clock moved.\n", " Either way, now I cannot read the time! I've implemented a deblurring\n", " algorithm--can you adjust its parameters to help me pin-point\n", " the time?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Here's code for a swirl transform:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from skimage import transform\n", "\n", "def swirl(xy, center=[0, 0], strength=1, radius=100, rotation=0):\n", " \"\"\"Compute the coordinate mapping for a swirl transformation.\n", "\n", " \"\"\"\n", " x, y = xy.T\n", " x0, y0 = center\n", " rho = np.sqrt((x - x0)**2 + (y - y0)**2)\n", "\n", " # Ensure that the transformation decays to approximately 1/1000-th\n", " # within the specified radius.\n", " radius = radius / 5 * np.log(2)\n", "\n", " theta = rotation + strength * \\\n", " np.exp(-rho / radius) + \\\n", " np.arctan2(y - y0, x - x0)\n", "\n", " xy[..., 0] = x0 + rho * np.cos(theta)\n", " xy[..., 1] = y0 + rho * np.sin(theta)\n", "\n", " return xy\n", "\n", "\n", "h, w = face.shape[:2]\n", "\n", "parameters = {'center': [w/2., h/2.],\n", " 'strength': 8,\n", " 'radius': 90,\n", " 'rotation': 0}\n", "\n", "out = transform.warp(face, swirl, parameters)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "f, (ax0, ax1) = plt.subplots(1, 2, figsize=(8, 4))\n", "\n", "ax0.imshow(face)\n", "ax1.imshow(out);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Can you come up with an even better distortion?\n", "\n", "## Start with this template:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def my_warp(xy):\n", " x = xy[:, 0]\n", " y = xy[:, 1]\n", " \n", " x = x + 1.5 * np.sin(y / 3)\n", " \n", " return np.hstack((x, y))\n", "\n", "image = plt.imread('../images/stefan.jpg')\n", "out = transform.warp(image, my_warp)\n", "\n", "f, (ax0, ax1) = plt.subplots(1, 2, figsize=(8, 4))\n", "ax0.imshow(image)\n", "ax1.imshow(out);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Composing Transformations\n", "\n", "scikit-image allows you to compose several transformations. For example:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from skimage import data\n", "\n", "cat = data.chelsea()\n", "horizontal_shift = transform.SimilarityTransform(translation=[20, 0])\n", "\n", "multiple_shifts = horizontal_shift + horizontal_shift + horizontal_shift\n", "\n", "f, (ax0, ax1, ax2) = plt.subplots(1, 3, figsize=(15, 5))\n", "ax0.imshow(cat)\n", "ax1.imshow(transform.warp(cat, horizontal_shift.inverse)) # Note the inverse!\n", "ax2.imshow(transform.warp(cat, multiple_shifts.inverse));" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The `transform` module allows us to rotate images. The inner workings is something like this:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def my_rotate(image, angle):\n", " rotation_tf = transform.SimilarityTransform(rotation=np.deg2rad(angle))\n", " return transform.warp(image, rotation_tf.inverse)\n", "\n", "plt.imshow(my_rotate(cat, 30))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that this rotates the cat around the origin (top-left).\n", "\n", "**Can you modify `my_rotate` to rotate the image around the center?**\n", "\n", "*Hint:*\n", "\n", "1. Shift the image (see above) so that the center of the image lies at (0, 0)\n", "2. Rotate the image\n", "3. Shift the image back---the opposite of what you did in step 1\n", "\n", "All of this can be achieved by composing transformations and calling `warp` once." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Advanced challenge: rectifying an image\n", "\n", "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We know the above tiles are laid out in a square--can you transform\n", "the image so that the tiles are displayed as if you were viewing them from above?\n", "\n", "The centre-points of the corner circles are, given as (row, column) coordinates:\n", "\n", "```\n", "(72, 129) -- top left\n", "(76, 302) -- top right\n", "(185, 90) -- bottom left\n", "(193, 326) -- bottom right\n", "```\n", "\n", "Hint: there is a linear transformation matrix, $H$, such that\n", "\n", "$H \\mathbf{x} = \\mathbf{x}'$\n", "\n", "where $\\mathbf{x}$ is the *homogeneous* coordinate in the original image and\n", "$\\mathbf{x}'$ is the *homogeneous* coordinate in the rectified image (with *homogeneous*\n", "we simply mean that we add an extra 1 at the end, e.g. (72, 129) becomes (72, 129, 1).\n", "The values for $\\mathbf{x}$ and their new values, $\\mathbf{x}'$,\n", "are therefore:\n", "\n", "```\n", "x = (72, 129, 1), x' = (0, 0, 1)\n", "x = (76, 302, 1), x' = (0, 400, 1)\n", "x = (185, 90, 1), x' = (400, 0, 1)\n", "x = (193, 326, 1) x' = (400, 400, 1)\n", "```\n", "\n", "(You can choose any output size you like--I chose $400 \\times 400$)\n", "\n", "Why do we need homogeneous coordinates? It allows us to have *translation* as part of H:\n", "\n", "$\n", "\\left[\\begin{array}{ccc}\n", "H_{00} & H_{01} & H_{02}\\\\\n", "H_{10} & H_{11} & H_{12}\\\\\n", "H_{20} & H_{21} & 1\n", "\\end{array}\\right]\\left[\\begin{array}{c}\n", "x\\\\\n", "y\\\\\n", "1\n", "\\end{array}\\right]=\\left[\\begin{array}{c}\n", "H_{00}x+H_{01}y+H_{02}\\\\\n", "H_{10}x+H_{11}y+H_{12}\\\\\n", "H_{20}x+H_{21}y+H_{22}\n", "\\end{array}\\right]\n", "$\n", "\n", "Note that each element of the output coordinate is of the form $ax + by + c$. Without the 1 in the last position of the coordinate, there would have been no $+ c$ and therefore no translation!\n", "\n", "The question on how to determine $H$ is left for another day. If you are curious, \n", "the [answer can be found here](homography.pdf).\n", "\n", "In the meantime, I provide some code to calculate $H$:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from skimage.transform import estimate_transform\n", "\n", "source = np.array([(129, 72),\n", " (302, 76),\n", " (90, 185),\n", " (326, 193)])\n", "\n", "target = np.array([[0, 0],\n", " [400, 0],\n", " [0, 400],\n", " [400, 400]])\n", "\n", "tf = estimate_transform('projective', source, target)\n", "H = tf.params\n", "print(H)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Using the code in the cell above, you can compute the target coordinate of any position in the original image." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Verify that the top left corner maps to (0, 0)\n", "\n", "x = np.array([[129, 72, 1]])\n", "\n", "z = np.dot(H, x.T)\n", "z /= z[2]\n", "\n", "print(z)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Here's a template solution:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def rectify(xy):\n", " x = xy[:, 0]\n", " y = xy[:, 1]\n", " \n", " # We need to provide the backward mapping, from the target\n", " # image to the source image.\n", " HH = np.linalg.inv(H)\n", " \n", " # You must fill in your code here to take\n", " # the matrix HH (given above) and to transform\n", " # each coordinate to its new position.\n", " # \n", " # Hint: handy functions are\n", " #\n", " # - np.dot (matrix multiplication)\n", " # - np.ones_like (make an array of ones the same shape as another array)\n", " # - np.column_stack\n", " # - A.T -- type .T after a matrix to transpose it\n", " # - x.reshape -- reshapes the array x\n", " \n", " # ... your code\n", " # ... your code\n", " # ... your code\n", " # ... your code\n", " \n", " return ...\n", "\n", " \n", "image = plt.imread('images/chapel_floor.png')\n", "out = transform.warp(image, rectify, output_shape=(400, 400))\n", "\n", "f, (ax0, ax1) = plt.subplots(1, 2, figsize=(8, 4))\n", "ax0.imshow(image)\n", "ax1.imshow(out)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
 
\n", "\n", "
\n", "The solution to the above problem is provided as [solutions/tile_rectify.py](solutions/tile_rectify.py). Only look at it after you've attempted the problem yourself!\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# For more fun examples see http://scikit-image.org/docs/dev/auto_examples" ] }, { "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.5.1+" }, "widgets": { "state": {}, "version": "1.1.2" } }, "nbformat": 4, "nbformat_minor": 0 }