{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# Image Processing SS 16 - Assignment - 02\n",
    "\n",
    "### Deadline is 4.5.2016 at 16:00 o'clock\n",
    "\n",
    "Please solve the assignments together with a partner.\n",
    "I will run every notebook. Make sure the code runs through. Select `Kernel` -> `Restart & Run All` to test it.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercise 1 - 10 Points\n",
    "\n",
    "Implement affine transformation with [bicubic interpolation](https://en.wikipedia.org/wiki/Bicubic_interpolation)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# display the plots inside the notebook\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pylab\n",
    "pylab.rcParams['figure.figsize'] = (12, 12)   # This makes the plot bigger"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The [skimage](http://scikit-image.org/) library comes with multiple useful test images.  Let's start with an image of an astronaut. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from skimage.data import astronaut\n",
    "from skimage.color import rgb2gray"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# We use a gray image. All the algorithms should work with color images too.\n",
    "img = rgb2gray(astronaut() / 255.)\n",
    "plt.imshow(img, cmap='gray')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def derive_y(image):\n",
    "    \"\"\"Computes the derivative of the image w.r.t the y coordinate\"\"\"\n",
    "    derived_image = np.zeros_like(image)\n",
    "    for x in range(image.shape[0]):\n",
    "        for y in range(image.shape[1]):\n",
    "            if y + 1 < image.shape[1] and y - 1 > 0:\n",
    "                derived_image[x,y] = image[x, y - 1] - image[x, y + 1]\n",
    "    return derived_image\n",
    "\n",
    "def derive_x(image):\n",
    "    \"\"\"Computes the derivative of the image w.r.t the x coordinate\"\"\"\n",
    "    derived_image = np.zeros_like(image)\n",
    "    for x in range(image.shape[0]):\n",
    "        for y in range(image.shape[1]):\n",
    "            if x + 1 < image.shape[1] and x - 1 > 0:\n",
    "                derived_image[x,y] = image[x - 1, y] - image[x + 1, y]\n",
    "    return derived_image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "dx_img = derive_x(img)\n",
    "dy_img = derive_y(img)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(18, 12))\n",
    "plt.subplot(131)\n",
    "plt.imshow(img, cmap='gray')\n",
    "plt.subplot(132)\n",
    "plt.imshow(dx_img, cmap='gray')\n",
    "plt.subplot(133)\n",
    "plt.imshow(dy_img, cmap='gray')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# The derivatives are no longer in the range [0,1].\n",
    "print(\"min: {}, max: {}\".format(dx_img.min(), dx_img.max()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "indicies = np.indices(img.shape).reshape(2, -1)\n",
    "indicies.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# pick some random index"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "indicies[:, 123456]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "indicies_hg = np.concatenate([\n",
    "        indicies, np.ones((1, indicies.shape[1]))], axis=0)\n",
    "indicies_hg.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "indicies_hg[:, 123456]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "T_scale = np.array([\n",
    "    [0.75, 0, 0],\n",
    "    [0, 0.75, 0],\n",
    "    [0, 0, 1],\n",
    "])\n",
    "# np.dot(T_affine, indicies_hg).shape, for python < 3.5\n",
    "(T_scale @ indicies_hg).shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "T_affine = np.array([\n",
    "    [0.75, 0.2, 100],\n",
    "    [-0.2, 0.75, 100],\n",
    "    [0, 0.001, 1],\n",
    "])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# you can use this function to invert the matricies\n",
    "np.linalg.inv(T_scale)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def affine_transformation(img, matrix):\n",
    "    indicies = np.indices(img.shape).reshape(2, -1)\n",
    "    indicies_hg = matrix @ np.concatenate([indicies, np.ones((1, indicies.shape[1]))], axis=0)\n",
    "    \n",
    "    res_img = np.zeros((int(np.ceil(np.max(indicies_hg[0,:]))), int(np.ceil(np.max(indicies_hg[1,:])))))\n",
    "    indicies_res = np.indices(res_img.shape).reshape(2, -1)\n",
    "    indicies_res_hg = np.concatenate([indicies_res, np.ones((1, indicies_res.shape[1]))], axis=0)\n",
    "    indicies_interpolation = np.linalg.inv(matrix) @ indicies_res_hg\n",
    "    \n",
    "    return bicubic_interpolation(img, indicies_interpolation, matrix, res_img)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def bicubic_interpolation(img, indicies, matrix, result):\n",
    "    dx_img = derive_x(img)\n",
    "    dy_img = derive_y(img)\n",
    "    dxy_img = derive_x(dy_img)\n",
    "    inv_matrix = np.array([\n",
    "        [ 1,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0],\n",
    "        [ 0,  0,  0,  0,  1,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0],\n",
    "        [-3,  3,  0,  0, -2, -1,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0],\n",
    "        [ 2, -2,  0,  0,  1,  1,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0],\n",
    "        [ 0,  0,  0,  0,  0,  0,  0,  0,  1,  0,  0,  0,  0,  0,  0,  0],\n",
    "        [ 0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  1,  0,  0,  0],\n",
    "        [ 0,  0,  0,  0,  0,  0,  0,  0, -3,  3,  0,  0, -2, -1,  0,  0],\n",
    "        [ 0,  0,  0,  0,  0,  0,  0,  0,  2, -2,  0,  0,  1,  1,  0,  0],\n",
    "        [-3,  0,  3,  0,  0,  0,  0,  0, -2,  0, -1,  0,  0,  0,  0,  0],\n",
    "        [ 0,  0,  0,  0, -3,  0,  3,  0,  0,  0,  0,  0, -2,  0, -1,  0],\n",
    "        [ 9, -9, -9,  9,  6,  3, -6, -3,  6, -6,  3, -3,  4,  2,  2,  1],\n",
    "        [-6,  6,  6, -6, -3, -3,  3,  3, -4,  4, -2,  2, -2, -2, -1, -1],\n",
    "        [ 2,  0, -2,  0,  0,  0,  0,  0,  1,  0,  1,  0,  0,  0,  0,  0],\n",
    "        [ 0,  0,  0,  0,  2,  0, -2,  0,  0,  0,  0,  0,  1,  0,  1,  0],\n",
    "        [-6,  6,  6, -6, -4, -2,  4,  2, -3,  3, -3,  3, -2, -1, -2, -1],\n",
    "        [ 4, -4, -4,  4,  2,  2, -2, -2,  2, -2,  2, -2,  1,  1,  1,  1]\n",
    "    ])\n",
    "\n",
    "    x_size = img.shape[0]\n",
    "    y_size = img.shape[1]\n",
    "    \n",
    "    for i in range(indicies.shape[-1]):\n",
    "        point = indicies[:, i]\n",
    "        x_val_floor = int(np.floor(point[0]))\n",
    "        x_val_ceil = int(np.ceil(point[0]))\n",
    "        y_val_floor = int(np.floor(point[1]))\n",
    "        y_val_ceil = int(np.ceil(point[1]))\n",
    "\n",
    "        if 0 < x_val_floor < x_size and 0 < x_val_ceil < x_size and 0 < y_val_floor < y_size and 0 < y_val_ceil < y_size:\n",
    "            func_values = np.array([\n",
    "                img[x_val_floor][y_val_floor], img[x_val_ceil][y_val_floor], img[x_val_floor][y_val_ceil], img[x_val_ceil][y_val_ceil],\n",
    "                dx_img[x_val_floor][y_val_floor], dx_img[x_val_ceil][y_val_floor], dx_img[x_val_floor][y_val_ceil], dx_img[x_val_ceil][y_val_ceil],\n",
    "                dy_img[x_val_floor][y_val_floor], dy_img[x_val_ceil][y_val_floor], dy_img[x_val_floor][y_val_ceil], dy_img[x_val_ceil][y_val_ceil],\n",
    "                dxy_img[x_val_floor][y_val_floor], dxy_img[x_val_ceil][y_val_floor], dxy_img[x_val_floor][y_val_ceil], dxy_img[x_val_ceil][y_val_ceil]\n",
    "            ])\n",
    "            alpha = inv_matrix @ func_values\n",
    "            res_point = matrix @ point\n",
    "            result[int(np.rint(res_point[0]))][int(np.rint(res_point[1]))] = calc_value(point[0] - x_val_floor, point[1] - y_val_floor, alpha)\n",
    "\n",
    "    return result\n",
    "\n",
    "def calc_value(x, y, alpha):\n",
    "    x_pow2 = x ** 2\n",
    "    x_pow3 = x ** 3\n",
    "    y_pow2 = y ** 2\n",
    "    y_pow3 = y ** 3\n",
    "    \n",
    "    return alpha[0] + alpha[4] * y +  alpha[8] * y_pow2 + alpha[12] * y_pow3\\\n",
    "        + (alpha[1] + alpha[5] * y +  alpha[9] * y_pow2 + alpha[13] * y_pow3) * x\\\n",
    "        + (alpha[2] + alpha[6] * y + alpha[10] * y_pow2 + alpha[14] * y_pow3) * x_pow2\\\n",
    "        + (alpha[3] + alpha[7] * y + alpha[11] * y_pow2 + alpha[15] * y_pow3) * x_pow3\\"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "img_scale = affine_transformation(img, T_scale)\n",
    "img_affine = affine_transformation(img, T_affine)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.imshow(img_scale, cmap='gray')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "plt.imshow(img_affine, cmap='gray')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "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+"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}