{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Discrete Cosine Transform" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exercise 2: Simple JPEG Encoder" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The Discrete Cosine Transform is the key to JPEG compression. Blocks of $8 \\times 8$ pixels are transformed to the frequency domain to be in a format which is more suited for compression. The code template provides you with the pixel values of the given $8 \\times 8$ block of the well-known test image \\emph{Lena}.\n", "\n", "a) In a greyscale image, the pixel values are usually encoded with 8 bit with values in $[0, 255]$. The DCT, however, works on $[-128, 127]$. Write a function that normalises the pixel values.\n", "\n", "b) The DCT transforms an $8 \\times 8$ block from the spatial domain to the frequency domain. Use\n", "\n", "$$\$$F_{uv} = \\frac{1}{\\sqrt{2N}} c_u c_v \\sum_{x=0}^{N-1} \\sum_{y=0}^{N-1} f_{xy} \\cos \\left(\\frac{(2x+1)u\\pi}{2N}\\right) \\cos \\left(\\frac{(2y+1)v\\pi}{2N}\\right)\n", "\$$$$\n", "\n", "where \n", "\n", "$$\$$\n", "c_{u;v} = \\begin{cases} \\frac{1}{\\sqrt{2}} & u;v = 0 \\\\ 1 & \\text{otherwise.}\\end{cases}\n", "\$$$$\n", "What is the complexity of your DCT routine? What is the meaning of $F_{00}$? \n", "\n", "c) The coefficients $F_{uv}$ are divided by quantisation values $Q_{uv}$ and rounded to the nearest integer. Quantisation is a lossy process; high quantisation coefficients result in a high compression factor, though at the expense of image quality. Implement the quantisation step. A common choice for the quantisation matrix is given in the code template. As mentioned in the lecture, a full JPEG encoder would now apply further compression techniques. \n", "\n", "d) Derive the IDCT and write a decoder for our simple JPEG block encoder. Compare the original image with the subsequently encoded and then decoded image." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPgAAAD8CAYAAABaQGkdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAC7dJREFUeJzt3d2LXfUVxvHnybyYFxOTtJMSMyHJRQjRQo2EgKQIjVhi\nFe1FLxJQqBRypSgtiPau/4DYiyJI1AqmShsVRKxi8a1Ca03GpDUvljQoSdAmWmJ0YjLJzOrFnMjo\nRGZPzv7tc2bx/cCQOWc2Z60zyZPfnj177+WIEICcZnW6AQDlEHAgMQIOJEbAgcQIOJAYAQcSI+BA\nYgQcSIyAA4n1FnnR3t7o7+8v8dKTjI2NNVJHkmbNavb/w6a+h5LU21vkn8JFjY6OpqwlSU2dGXrm\nzBmdO3fOU21X5G+1v79fa9asKfHSk5w+fbqROpI0Z86cxmpJ0ooVKxqrtXDhwsZqDQ8PN1br5MmT\njdWSpHPnzjVSZ2hoqNJ27KIDiRFwIDECDiRGwIHECDiQGAEHEiPgQGIEHEisUsBtb7b9vu1Dtu8v\n3RSAekwZcNs9kn4n6SZJV0naavuq0o0BaF+VFXyDpEMRcTgiRiQ9Lem2sm0BqEOVgC+TdGTC46Ot\n5wB0udouNrG9TdI2Serr66vrZQG0ocoKfkzS8gmPB1vPfU1EPBIR6yNifZOXHgL4dlUC/o6k1bZX\n2e6XtEXS82XbAlCHKZfaiDhv+y5JL0vqkfRYROwr3hmAtlXal46IFyW9WLgXADXjTDYgMQIOJEbA\ngcQIOJAYAQcSI+BAYgQcSIyAA4kVO2m8qREuIyMjjdTphLNnz3a6hSKavFahp6ensVqSdP78+Ubr\nTYUVHEiMgAOJEXAgMQIOJEbAgcQIOJAYAQcSI+BAYgQcSKzKZJPHbB+3/V4TDQGoT5UV/PeSNhfu\nA0ABUwY8It6U9L8GegFQM34GBxJjdBGQWG0rOKOLgO7DLjqQWJVfkz0l6W+S1tg+avsX5dsCUIcq\ns8m2NtEIgPqxiw4kRsCBxAg4kBgBBxIj4EBiBBxIjIADiRFwILEiJ43PmjVLs2fPLvHSkzQ53ifz\nOfZNfh/PnDnTWK2xsbHGaknj//abYLvSdqzgQGIEHEiMgAOJEXAgMQIOJEbAgcQIOJAYAQcSI+BA\nYgQcSKzKTReX237N9n7b+2zf00RjANpX5eTq85J+FRFDtudL2m37lYjYX7g3AG2qMpvso4gYan3+\nuaQDkpaVbgxA+6Z1eZTtlZLWSXr7Il/7anRRf39/Da0BaFflg2y2L5f0jKR7I+LUN78+cXQRs8mA\n7lAp4Lb7NB7uHRHxbNmWANSlylF0S3pU0oGIeLB8SwDqUmUF3yjpDkmbbO9pffykcF8AalBlNtlb\nkqrdHwZAV+FMNiAxAg4kRsCBxAg4kBgBBxIj4EBiBBxIjIADiRUZttXT06PFixeXeOlJFi1a1Egd\nSY3NW7tgzpw5jdVqcl7YqVOTrlUqpukrGwcGBhqpU/WCLlZwIDECDiRGwIHECDiQGAEHEiPgQGIE\nHEiMgAOJEXAgsSo3XZxt+x+297ZGF/2micYAtK/KqapnJW2KiC9at09+y/afI+LvhXsD0KYqN10M\nSV+0Hva1PqJkUwDqUXXwQY/tPZKOS3olIi46usj2Ltu7RkZG6u4TwCWoFPCIGI2IayQNStpg+/sX\n2ear0UXMJgO6w7SOokfESUmvSdpcph0AdapyFH3A9sLW53Mk3SjpYOnGALSvylH0pZKesN2j8f8Q\n/hgRL5RtC0AdqhxF/6fGZ4IDmGE4kw1IjIADiRFwIDECDiRGwIHECDiQGAEHEiPgQGJFRhfZrjxa\npV3z5s1rpI4kLViwoLFakjR+pW4zmhwn1OTVhkuWLGmsliStXbu2kTpvvPFGpe1YwYHECDiQGAEH\nEiPgQGIEHEiMgAOJEXAgMQIOJEbAgcQqB7x1b/R3bXM/NmCGmM4Kfo+kA6UaAVC/qpNNBiXdLGl7\n2XYA1KnqCv6QpPskjRXsBUDNqgw+uEXS8YjYPcV2zCYDukyVFXyjpFttfyDpaUmbbD/5zY2YTQZ0\nnykDHhEPRMRgRKyUtEXSqxFxe/HOALSN34MDiU3rji4R8bqk14t0AqB2rOBAYgQcSIyAA4kRcCAx\nAg4kRsCBxAg4kBgBBxIrMrpodHRUn332WYmXnuTLL79spI4knT59urFa0vgIqKY0+d6avFZhYGCg\nsVqStHz58kbqVP0esoIDiRFwIDECDiRGwIHECDiQGAEHEiPgQGIEHEiMgAOJVTqTrXVH1c8ljUo6\nHxHrSzYFoB7TOVX1RxHxSbFOANSOXXQgsaoBD0l/sb3b9raSDQGoT9Vd9B9GxDHbSyS9YvtgRLw5\ncYNW8LdJ0mWXXVZzmwAuRaUVPCKOtf48Luk5SRsuss1Xo4v6+vrq7RLAJakyfHCe7fkXPpf0Y0nv\nlW4MQPuq7KJ/T9JzrZsP9Er6Q0S8VLQrALWYMuARcVjSDxroBUDN+DUZkBgBBxIj4EBiBBxIjIAD\niRFwIDECDiRGwIHEiowu6u3t1eLFi0u89CQnTpxopI4kffrpp43VkqQFCxY0VqvJET9Lly5trNbq\n1asbqyVJV155ZSN1ql7vwQoOJEbAgcQIOJAYAQcSI+BAYgQcSIyAA4kRcCAxAg4kVingthfa3mn7\noO0Dtq8r3RiA9lU9VfW3kl6KiJ/Z7pc0t2BPAGoyZcBtXyHpekk/l6SIGJE0UrYtAHWosou+StIJ\nSY/bftf29tb90QF0uSoB75V0raSHI2KdpGFJ939zI9vbbO+yvevs2bM1twngUlQJ+FFJRyPi7dbj\nnRoP/NdMHF3EbDKgO0wZ8Ij4WNIR22taT90gaX/RrgDUoupR9Lsl7WgdQT8s6c5yLQGoS6WAR8Qe\nSesL9wKgZpzJBiRGwIHECDiQGAEHEiPgQGIEHEiMgAOJEXAgMQIOJFZkNtncuXO1bt26Ei89yd69\nexupI0kffvhhY7Ukaf78+Y3VuvrqqxurtXbt2sZqDQ4ONlZLkhYtWtRInf7+/krbsYIDiRFwIDEC\nDiRGwIHECDiQGAEHEiPgQGIEHEiMgAOJTRlw22ts75nwccr2vU00B6A9U56qGhHvS7pGkmz3SDom\n6bnCfQGowXR30W+Q9J+IaPakbACXZLoB3yLpqYt9YeLoouHh4fY7A9C2ygFvDT24VdKfLvb1iaOL\n5s1jNiHQDaazgt8kaSgi/luqGQD1mk7At+pbds8BdKdKAW/NA79R0rNl2wFQp6qzyYYlfadwLwBq\nxplsQGIEHEiMgAOJEXAgMQIOJEbAgcQIOJAYAQcSc0TU/6L2CUnTvaT0u5I+qb2Z7pD1vfG+OmdF\nRAxMtVGRgF8K27siYn2n+ygh63vjfXU/dtGBxAg4kFg3BfyRTjdQUNb3xvvqcl3zMziA+nXTCg6g\nZl0RcNubbb9v+5Dt+zvdTx1sL7f9mu39tvfZvqfTPdXJdo/td22/0Ole6mR7oe2dtg/aPmD7uk73\n1I6O76K37rX+b43fMeaopHckbY2I/R1trE22l0paGhFDtudL2i3ppzP9fV1g+5eS1ktaEBG3dLqf\nuth+QtJfI2J760ajcyPiZKf7ulTdsIJvkHQoIg5HxIikpyXd1uGe2hYRH0XEUOvzzyUdkLSss13V\nw/agpJslbe90L3WyfYWk6yU9KkkRMTKTwy11R8CXSToy4fFRJQnCBbZXSlon6e3OdlKbhyTdJ2ms\n043UbJWkE5Ieb/34sb11P8IZqxsCnprtyyU9I+neiDjV6X7aZfsWSccjYneneymgV9K1kh6OiHWS\nhiXN6GNC3RDwY5KWT3g82HpuxrPdp/Fw74iILHek3SjpVtsfaPzHqU22n+xsS7U5KuloRFzY09qp\n8cDPWN0Q8Hckrba9qnVQY4uk5zvcU9tsW+M/yx2IiAc73U9dIuKBiBiMiJUa/7t6NSJu73BbtYiI\njyUdsb2m9dQNkmb0QdFKt00uKSLO275L0suSeiQ9FhH7OtxWHTZKukPSv2zvaT3364h4sYM9YWp3\nS9rRWmwOS7qzw/20peO/JgNQTjfsogMohIADiRFwIDECDiRGwIHECDiQGAEHEiPgQGL/B6xowHRU\nzGExAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import numpy as np\n", "from matplotlib import pyplot as plt\n", "from numpy import dtype\n", "\n", "\n", "Cosines = np.zeros((8,8))\n", "Coeffs = np.ones((8,8))\n", "qtable = np.matrix( [[ 16, 11, 10, 16, 24, 40, 51, 61],\n", " [ 12, 12, 14, 19, 26, 58, 60, 55],\n", " [ 14, 13, 16, 24, 40, 57, 69, 56],\n", " [ 14, 17, 22, 29, 51, 87, 80, 62],\n", " [ 18, 22, 37, 56, 68, 109, 103, 77],\n", " [ 24, 35, 55, 64, 81, 104, 113, 92],\n", " [ 49, 64, 78, 87, 103, 121, 120, 101],\n", " [ 72, 92, 95, 98, 112, 100, 103, 99]])\n", "\n", "\n", "def normalise(block):\n", " '''\n", " Shift pixel values of an 8-by-8 matrix to [-128, 127]\n", " '''\n", " return block - 128*np.ones((8,8), dtype=np.int)\n", "\n", "\n", "def denormalise(block):\n", " '''\n", " Shift pixel values to [0, 255]\n", " '''\n", " return block + 128*np.ones((8,8), dtype=np.int)\n", "\n", "\n", "def dct(block):\n", " '''\n", " DCT matrix for an 8-by-8 block\n", " '''\n", " dctBlock = np.zeros((8,8)) \n", " \n", " for i in range(0,8):\n", " for j in range(0,8):\n", " blockSum = 0.\n", " for x in range(0,8):\n", " for y in range(0,8):\n", " pixel = block[x,y]\n", " blockSum = blockSum + pixel * Cosines[x,i] * Cosines[y,j]\n", " dctBlock[i,j] = int(0.25 * Coeffs[i,j] * blockSum)\n", " \n", " return dctBlock\n", "\n", "\n", "def idct(block):\n", " '''\n", " inverse DCT for an 8-by-8 block\n", " '''\n", " image = np.zeros((8,8)) \n", " \n", " for x in range(0,8):\n", " for y in range(0,8):\n", " blockSum = 0.\n", " for u in range(0,8):\n", " for v in range(0,8):\n", " blockSum = blockSum + Coeffs[u,v] * block[u,v] * Cosines[x,u] * Cosines[y,v]\n", " image[x,y] = int(0.25 * blockSum)\n", " \n", " return image\n", "\n", " \n", "def quantise(block):\n", " '''\n", " Lossy compression through quantisation of DCT-coefficients\n", " ''' \n", " return np.round(block/qtable, 0)\n", "\n", "\n", "def dequantise(block):\n", " '''\n", " Dequantise the DCT coefficients\n", " '''\n", " return np.multiply(qtable, block)\n", "\n", "\n", "def encodeBlock(block):\n", " '''\n", " JPEG-encode an 8-by-8 block (no lossless compression)\n", " '''\n", " normalisedBlock = normalise(block)\n", " dctBlock = dct(normalisedBlock)\n", " quantisedBlock = quantise(dctBlock)\n", " return quantisedBlock\n", "\n", "\n", "def decodeBlock(block):\n", " '''\n", " Decode an 8-by-8 block for your simple JPEG encoder\n", " '''\n", " return denormalise(idct(dequantise(block)))\n", "\n", "\n", "# values in [0, 255]\n", "image = np.matrix([[ 37, 41, 53, 68, 89, 79, 54, 68],\n", " [ 36, 28, 38, 65, 100, 94, 62, 65],\n", " [ 46, 36, 46, 66, 84, 76, 59, 78],\n", " [ 84, 67, 73, 84, 86, 73, 64, 87],\n", " [ 77, 86, 116, 113, 77, 57, 79, 129],\n", " [ 73, 62, 76, 74, 58, 70, 109, 155],\n", " [ 82, 54, 55, 57, 63, 105, 154, 189],\n", " [121, 73, 65, 86, 126, 181, 203, 200]], dtype=int)\n", "\n", "plt.imshow(image, cmap='gray', interpolation='nearest', vmin=0, vmax=255)\n", "plt.savefig('origblock8x8.jpg')\n", "#plt.show()\n", "\n", "# compute lookup tables\n", "for x in range(0,8):\n", " for y in range(0,8):\n", " Cosines[x,y] = np.cos((2.*x+1)*y*np.pi/16)\n", " \n", "vec = np.array([1/np.sqrt(2), 1, 1, 1, 1, 1, 1, 1])\n", "Coeffs = np.outer(vec, vec)\n", "\n", "# encode and decode the given block \n", "encodedBlock = encodeBlock(image)\n", "decodedImage = decodeBlock(encodedBlock)\n", "plt.imshow(decodedImage, cmap='gray', interpolation='nearest', vmin=0, vmax=255)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exercise 3: Discrete Cosine Transform" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We start with a dataset $f_{-N+1}, \\ldots, f_N$, which fulfills the following\n", "symmetry constraint:\n", "\$$\n", " f_{-n} = f_n \\qquad \\mbox{for n=1,\\ldots, N-1}\n", "\$$\n", "\n", "a) Show that the corresponding Fourier coefficients\n", " $$\$$\n", " F_k = \\frac{1}{2N}\n", " \\sum\\limits_{n=-N+1}^{N} f_n \\omega_{2N}^{-kn}\n", " \$$$$\n", " are real values only and can be written as:\n", "\n", " $$\$$\n", " F_k = \\frac{1}{N} \\left(\\frac{1}{2} f_0 + \\sum\\limits_{n=1}^{N-1} f_n \\cos\\left( \\frac{\\pi nk}{N} \\right)+ \\frac{1}{2} f_N \\cos(\\pi k) \\right).\n", " \$$$$\n", "\n", "b) Show that the $F_k$ is symmetric too.\n", "\n", "c)\n", "Let $FFT(f,N)$ be a procedure that computes the coefficients $F_k$\n", "efficiently (according to dft equation from a field $f$\n", "which consists of $2N$ values $f_n$. (The result is written back to field $f$)\n", "\n", "Write a short procedure $DCT(g,N)$ which uses procedure $FFT$\n", "to compute the coefficients $F_k$ for $k=0,\\ldots, N$ from dft equation for the \n", "(non-symmetrical) data $f_0,\\ldots,f_N$, stored in the parameter field\n", "$g$." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "f = [0.0, -1.0, 2.0, 4.0, 5.0, 6.0, 5.0, 4.0, 2.0, -1.0]\n", "naive_dct: [ 26. -16.94427191 -1.23606798 0.94427191 3.23606798 2. ]\n", "dct : [ 26. -16.94427191 -1.23606798 0.94427191 3.23606798 2. ]\n" ] } ], "source": [ "from scipy.fftpack import *\n", "\n", "def naive_dct(g):\n", " f = g + g[-2:0:-1]\n", " print(\"f = \" + str(f))\n", " F = fft(f)\n", " return F[:len(g)]\n", " #return F\n", "\n", "g = [0.0, -1.0, 2.0, 4.0, 5.0, 6.0]\n", "\n", "print(\"naive_dct: \" + str(naive_dct(g).real))\n", "print(\"dct : \" + str(dct(g, type=1)))\n" ] } ], "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.6.4" } }, "nbformat": 4, "nbformat_minor": 1 }