{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "# Exercises from Linear algebra, signal processing, and wavelets. A unified approach.\\\\ Python version\n", "\n", " **Øyvind Ryan**\n", "\n", "Date: **Oct 24, 2018**\n", "\n", "\n", "\n", "\n", "# Digital images\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "## Example 8.2: Extracting the different colors\n", "
\n", "\n", "\n", "If we have a color image $P=(r_{i,j},g_{i,j},b_{i,j})_{i,j=1}^{m,n}$, it is often useful to manipulate the three color components separately as the three images" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "P_r=(r_{i,j})_{i,j=1}^{m,n},\\quad P_g=(g_{i,j})_{i,j=1}^{m,n},\\quad P_b=(b_{i,j})_{i,j=1}^{m,n}.\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As an example, let us first see how we can produce three separate images, showing the red, green, and blue color components of the Lena image. First we need to read the file:" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import os, sys\n", "sys.path.append(os.path.join(os.getcwd(), 'python'))" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline\n", "\n", "from images import *\n", "from fft import *\n", "from scipy.fftpack import dct, idct\n", "from numpy import *\n", "import matplotlib.pyplot as plt\n", "from sound import filter_impl\n", "from forward_compress_reverse import *" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [], "source": [ "img = double(imread('images/lena.png', 'png'))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The returned object has three dimensions, the first two representing the spatial directions (row/column), the third the color component. \n", "By using different values for the third dimension we can thus obtain the different color components as follows:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [], "source": [ "X1 = zeros_like(img)\n", "X1[:, :, 0] = img[:, :, 0]\n", "\n", "X2 = zeros_like(img)\n", "X2[:, :, 1] = img[:, :, 1]\n", "\n", "X3 = zeros_like(img)\n", "X3[:, :, 2] = img[:, :, 2]" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [], "source": [ "imshow(uint8(X1))\n", "imshow(uint8(X2))\n", "imshow(uint8(X3))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "\n", "\n", "\n", "\n", "\n", "## Example 8.3: Converting from color to grey-level\n", "
\n", "\n", "\n", "A color image can easily be converted to a grey-level image: At each point in the image we have to replace the three color values $(r,g,b)$ by a single value $q$ that will represent the grey level. \n", "If we want the grey-level image to be a reasonable representation of the color image, the value $q$ should somehow reflect the intensity of the image at the point. There are several ways to do this.\n", "If the image is $P=(r_{i,j},g_{i,j},b_{i,j})_{i,j=1}^{m,n}$, we can (for all $i$ and $j$) either\n", "* set $q_{i,j}=\\max(r_{i,j},g_{i,j},b_{i,j})$, i.e. use the largest of the three color components,\n", "\n", "* set $q_{i,j}=(r_{i,j}+g_{i,j}+b_{i,j})/3$, i.e. use the average of the three values,\n", "\n", "* set $q_{i,j}=\\sqrt{r_{i,j}^2+g_{i,j}^2+b_{i,j}^2}$, i.e. use the \"length\" of the color vector $(r,g,b)$.\n", "\n", "The three methods can be implemented as follows." ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false }, "outputs": [], "source": [ "mx = maximum(img[:, :, 0], img[:, :, 1])\n", "X1 = maximum(img[:, :, 2], mx)\n", "\n", "X2 = (img[:, :, 0] + img[ :, :, 1] + img[ :, :, 2])/3.\n", "\n", "X3 = sqrt(img[:,:,0]**2 + img[:, :, 1]**2 + img[:, :, 2]**2)\n", "map_to_01(X3)\n", "X3 *= 255" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For the third method we had to be careful, since $\\sqrt{r^2+g^2+b^2}$ may be outside $[0,1]$, even if $r$, $g$ and $b$ all lie in $[0,1]$. We thus need to normalize the values as explained above.\n", "In practice one of the last two methods are preferred, perhaps with a preference for the last method, but the actual choice depends on the application." ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false }, "outputs": [], "source": [ "imshow(uint8(X1))\n", "imshow(uint8(X2))\n", "imshow(uint8(X3))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "\n", "\n", "\n", "\n", "\n", "## Example 8.4: Computing the negative image\n", "
\n", "\n", "\n", "In film-based photography a negative image was obtained when the film was developed, and then a positive image was created from the negative. \n", "We can easily simulate this and compute a negative digital image.\n", "\n", "Suppose we have a grey-level image $P=(p_{i,j})_{i,j=1}^{m,n}$ with intensity values in the interval $[0,1]$. Here intensity value 0 corresponds to black and 1 \n", "corresponds to white. To obtain the negative image we just have to replace an intensity $p$ by its 'mirror value' $1-p$.\n", "This is also easily translated to code as above." ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false }, "outputs": [], "source": [ "X1 = 255 - X1\n", "X2 = 255 - X2\n", "X3 = 255 - X3\n", "\n", "imshow(uint8(X1))\n", "imshow(uint8(X2))\n", "imshow(uint8(X3))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "\n", "\n", "\n", "\n", "\n", "## Example 8.5: Adjusting the contrast\n", "
\n", "\n", "\n", "A common problem with images is that the contrast often is not good enough. \n", "This typically means that a large proportion of the grey values are concentrated in a rather small subinterval of $[0,1]$. \n", "The obvious solution is to somehow spread out the values. \n", "This can be accomplished by applying a monotone function $f$ which maps $[0,1]$ onto $[0,1]$. \n", "If we choose $f$ so that its derivative is large in the area where many intensity values are concentrated, we obtain the desired effect. \n", "We will consider two such families of functions:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\begin{align*}\n", "f_n(x) &= \\frac{\\arctan(n(x-1/2))}{2\\arctan(n/2)}+ \\frac{1}{2} \\\\ \n", "g_\\epsilon(x) &= \\frac{\\ln(x+\\epsilon)-\\ln \\epsilon}{\\ln(1+\\epsilon)-\\ln \\epsilon}.\n", "\\end{align*}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The first type of functions have quite large derivatives near $x=0.5$ and will therefore increase the contrast in images with a concentration of intensities near 0.5. \n", "The second type of functions have a large derivative near $x=0$ and will therefore increase the contrast in images with a large proportion of small intensity values, i.e., very dark images. \n", "The following code plots the functions $f_4$, $f_{10}$, and $f_{100}$." ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def contrast_adjust0(X,n):\n", " # Assumes that the values are in [0,255]\n", " X /= 255.\n", " X -= 1/2.\n", " X *= n\n", " arctan(X, X)\n", " X /= (2*arctan(n/2.)) \n", " X += 1/2.0\n", " X *= 255 # Maps the values back to [0,255]" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": false }, "outputs": [], "source": [ "x=linspace(0, 1, 100)\n", "\n", "n=4\n", "f4 = arctan(n*(x-1/2.))/(2*arctan(n/2.)) + 1/2.\n", "n=10\n", "f10=arctan(n*(x-1/2.))/(2*arctan(n/2.)) + 1/2.\n", "n=100\n", "f100=arctan(n*(x-1/2.))/(2*arctan(n/2.)) + 1/2.\n", "plt.plot(x,f4,'r',x,f10,'g',x,f100,'b');\n", "plt.axis([0, 1, 0, 1])\n", "plt.legend(['n=4','n=10','n=100'], loc='lower right', frameon=False)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The following code plots the functions $g_{0.1}$, $g_{0.01}$, and $g_{0.001}$." ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": false }, "outputs": [], "source": [ "eps=0.1\n", "g01=(log(x+eps)-log(eps))/(log(1+eps)-log(eps))\n", "eps=0.01\n", "g001=(log(x+eps)-log(eps))/(log(1+eps)-log(eps))\n", "eps=0.001\n", "g0001=(log(x+eps)-log(eps))/(log(1+eps)-log(eps))\n", "plt.plot(x,g01,'r',x,g001,'g',x,g0001,'b')\n", "plt.axis([0, 1, 0, 1])\n", "plt.legend(['$\\epsilon$=0.1','$\\epsilon$=0.01','$\\epsilon$=0.001'], loc='lower right', frameon=False)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The following code applies $f_{10}$ to the image in the right part of Figure 8.4." ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": false }, "outputs": [], "source": [ "X1 = img.copy()\n", "contrast_adjust0(X1, 10)\n", "imshow(uint8(X1))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The following code applies $g_{0.01}$ to the image in the right part of Figure 8.4." ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": false }, "outputs": [], "source": [ "X2 = img.copy()\n", "contrast_adjust(X2, 0.01)\n", "imshow(uint8(X2))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Since the image was quite well balanced, $f_{10}$ made the dark areas too dark and the bright areas too bright.\n", "$g_{0.01}$ on the other hand has made the image as a whole too bright.\n", "\n", "Increasing the contrast is easy to implement. The following function uses the contrast adjusting function $g_\\epsilon$, with $\\epsilon$ as parameter." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " def contrast_adjust(X,epsilon):\n", " # Assumes that the values are in [0,255]\n", " X /= 255.\n", " X += epsilon\n", " log(X, X) \n", " X -= log(epsilon)\n", " X /= (log(1+epsilon)-log(epsilon))\n", " X *= 255\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "\n", "\n", "\n", "\n", "\n", "## Exercise 8.6: Generating black and white images\n", "
\n", "\n", "\n", "Black and white images can be generated from grayscale images (with values between $0$ and $255$) by replacing each pixel value with the one of $0$ and $255$ which\n", "is closest. Use this strategy to generate the black and white image shown in \n", "the right part of Figure 8.2.\n", "\n", "\n", "\n", "**Solution.**\n", "The following code can be used" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": false }, "outputs": [], "source": [ "X2 = 255*(X1 >= 128)\n", "imshow(uint8(X2))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "## Exercise 8.7: Adjusting contrast\n", "
\n", "\n", "\n", "\n", "**a)**\n", "Write a function which instead uses the function $f_n$ from Example 8.5 to adjust the contrast (rather than $g_{\\epsilon}$). $n$ should be a parameter to your function.\n", "\n", "**b)**\n", "Generate the left and right images in Figure 8.7 on your own by using code which calls the function you have written, as well as `contrast_adjust`.\n", "\n", "\n", "\n", "**Solution.**\n", "The following code can be used." ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "collapsed": false }, "outputs": [], "source": [ "X1 = img.copy()\n", "contrast_adjust0(X1, 10)\n", "imshow(uint8(X1))" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "collapsed": false }, "outputs": [], "source": [ "X2 = img.copy()\n", "contrast_adjust(X2, 0.01)\n", "imshow(uint8(X2))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "## Exercise 8.8: Adjusting contrast with another function\n", "
\n", "\n", "\n", "\n", "**a)**\n", "Show that the function $h_n: \\mathbb{R} \\rightarrow \\mathbb{R}$ given by" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "h_n(x) = x^n,\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "maps the interval $[0,1] \\rightarrow [0,1]$ for any $n$, and that $h_n^\\prime (1) \\rightarrow \\infty$ as $n \\rightarrow \\infty$.\n", "\n", "\n", "\n", "**Solution.**\n", "We have that $h_n^\\prime(x)=nx^{n-1}$, so that $h_n^\\prime(1)=n\\to\\infty$.\n", "\n", "\n", "\n", "**b)**\n", "The image `secret.jpg`, shown in Figure 8.8, \n", "contains some information that is nearly invisible to the naked eye on most monitors. Use\n", "the function $h_n$, to reveal the secret message.\n", "\n", "\n", "\n", "\n", "
\n", "\n", "

Secret message.

\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "**Hint.**\n", "Convert to a grayscale image before you adjust the contrast with $h_n$.\n", "\n", "\n", "\n", "\n", "\n", "**Solution.**\n", "The following code can be used" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "collapsed": false }, "outputs": [], "source": [ "n = 100\n", "X = double(imread('images/secret.jpg', 'jpg'))\n", "X = (X[:,:,0] + X[:,:,1]+ X[:,:,2])/3.\n", "map_to_01(X)\n", "X **= n\n", "X *= 255\n", "X = uint8(X)\n", "imshow(X)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "## Example 8.9: Smoothing an image\n", "
\n", "\n", "\n", "When we considered filtering of digital sound, low-pass filters dampened high frequencies. \n", "We will here similarly see that an image can be smoothed by applying a low-pass filters to the rows and the columns.\n", "Let us write down the corresponding computational molecules, and let us use low-pass filters with filter coefficients from Pascal's triangle. \n", "If we use the filter $S=\\frac{1}{4}\\{1,\\underline{2},1\\}$ (row 2 from Pascal's triangle), Theorem 2 gives the computational molecule" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "A=\n", "\\frac{1}{16}\n", "\\begin{pmatrix}\n", "1 & 2 & 1\\\\ \n", "2 & \\underline{4} & 2\\\\ \n", "1 & 2 & 1\n", "\\end{pmatrix}.\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If the pixels in the image are $p_{i,j}$, this means that we compute the new pixels by" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\begin{multline*}\n", "\\hat p_{i,j}=\\frac{1}{16}\\bigl(4 p_{i,j}+2(p_{i,j-1}+p_{i-1,j}+p_{i+1,j}+p_{i,j+1})\\\\ \n", "+p_{i-1,j-1}+p_{i+1,j-1}+p_{i-1,j+1}+p_{i+1,j+1}\\bigr).\n", "\\end{multline*}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If we instead use the filter $S=\\frac{1}{64}\\{1,6,15,\\underline{20},15,6,1\\}$ (row 6 from Pascal's triangle), we obtain the molecule" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\frac{1}{4096}\n", "\\begin{pmatrix}\n", "1 & 6 & 15 & 20 & 15 & 6 & 1 \\\\ \n", " 6 & 36 & 90 & 120 & 90 & 36 & 6 \\\\ \n", " 15 & 90 & 225 & 300 & 225 & 90 & 15 \\\\ \n", " 20 & 120 & 300 & \\underline{400} & 300 & 120 & 20 \\\\ \n", " 15 & 90 & 225 & 300 & 225 & 90 & 15 \\\\ \n", " 6 & 36 & 90 & 120 & 90 & 36 & 6 \\\\ \n", " 1 & 6 & 15 & 20 & 15 & 6 & 1\n", "\\end{pmatrix}.\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We anticipate that both molecules give a smoothing effect, but that the second molecule provides more smoothing. \n", "To make the smoothing effect visible, we have zoomed in on the face in the image. The following function was used to obtain this portion of the image, and will be used repeatedly in the following:" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def create_excerpt():\n", " img = double(imread('images/lena.png','png'))\n", " return img[128:384,128:384,:]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The result of applying the two molecules above to our grayscale image" ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "collapsed": false }, "outputs": [], "source": [ "excerpt=create_excerpt()\n", "X1 = excerpt.copy()\n", "imshow(uint8(X1))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "is shown below" ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def shortmolecule(x, bd_mode):\n", " filter_impl(array([1, 2, 1])/4., x, bd_mode)\n", " \n", "X2 = excerpt.copy()\n", "tensor2_impl(X2, shortmolecule, shortmolecule, 'symm')\n", "imshow(uint8(X2))" ] }, { "cell_type": "code", "execution_count": 21, "metadata": { "collapsed": false }, "outputs": [], "source": [ "X3 = excerpt.copy()\n", "\n", "S2 = convolve((1/8.)*array([1, 3, 3, 1]),(1/8.)*array([1, 3, 3, 1]))\n", "\n", "def longmolecule(x, bd_mode):\n", " filter_impl(S2, x, bd_mode)\n", " \n", "tensor2_impl(X3, longmolecule, longmolecule, 'symm')\n", "imshow(uint8(X3))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The smoothing effect is best visible in the second image. Smoothing effects are perhaps more visible if we use the simple chess pattern image\n", "generated as follows" ] }, { "cell_type": "code", "execution_count": 22, "metadata": { "collapsed": false }, "outputs": [], "source": [ "N=16\n", "img=255*ones((N, N))\n", "img[0:(int(N/2)),0:(int(N/2))] = 0\n", "img[(int(N/2)):N,(int(N/2)):N] = 0\n", " \n", "X1 = img.copy()\n", "imshow(uint8(X1))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The image can be smoothed out horizontally as follows," ] }, { "cell_type": "code", "execution_count": 23, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def donothing(x, bd_mode):\n", " a=1\n", " \n", "X2 = img.copy()\n", "tensor2_impl(X2, donothing, shortmolecule, 'per')\n", "imshow(uint8(X2))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "vertically as follows," ] }, { "cell_type": "code", "execution_count": 24, "metadata": { "collapsed": false }, "outputs": [], "source": [ "X3 = img.copy()\n", "tensor2_impl(X3, shortmolecule, donothing, 'per')\n", "imshow(uint8(X3))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "and both horizontally and vertically as follows." ] }, { "cell_type": "code", "execution_count": 25, "metadata": { "collapsed": false }, "outputs": [], "source": [ "X4 = img.copy()\n", "tensor2_impl(X4, shortmolecule, shortmolecule, 'per')\n", "imshow(uint8(X4))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Again we have used the filter $S=\\frac{1}{4}\\{1,\\underline{2},1\\}$. \n", "Here we also have shown what happens if we only smooth the image in one of the directions.\n", "In the right image we have smoothed in both directions. \n", "We then see the union of the two one-dimensional smoothing operations. \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "## Example 8.10: Edge detection\n", "
\n", "\n", "\n", "Another operation on images which can be expressed in terms of computational molecules is *edge detection*. \n", "An edge in an image is characterized by a large change in intensity values over a small distance in the image. \n", "For a continuous function this corresponds to a large derivative. \n", "We cannot compute exact derivatives for an image since it is only defined at isolated points, however,\n", "but we have a perfect situation for numerical differentiation techniques.\n", "\n", "\n", "### **Partial derivative in $x$-direction**\n", "\n", "Let us first compute the partial derivative $\\partial P/\\partial x$ at all points in the image.\n", "Note first that it is the second coordinate in an image which refers to the $x$-direction you are used to from plotting functions. This means that the much used symmetric Newton quotient\n", "approximation to the partial derivative [[knut]](#knut) takes the form" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\frac{\\partial P}{\\partial x}(i,j)\\approx \\frac{p_{i,j+1}-p_{i,j-1}}{2},\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "where we have used the convention $h=1$ (i.e. the derivative measures 'intensity per pixel'). \n", "This corresponds to applying the high-pass filter $S=\\frac{1}{2}\\{1,\\underline{0},-1\\}$ to all the rows, \n", "alternatively applying the tensor product $I\\otimes S$ to the image, i.e. the molecule" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\frac{1}{2}\n", "\\begin{pmatrix}\n", "0 & 0 & 0\\\\ \n", "1 & \\underline{0} & -1\\\\ \n", "0 & 0 & 0\n", "\\end{pmatrix}.\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We have included two rows of zeroes to make it clear how the computational molecule is to be applied when we place it over the pixels. \n", "If we apply this molecule to the usual excerpt of the Lena image" ] }, { "cell_type": "code", "execution_count": 26, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def diffxmolecule(x, bd_mode):\n", " filter_impl(array([1, 0, -1])/2., x, bd_mode)\n", " \n", "excerpt=create_excerpt()\n", "X1 = excerpt.copy()\n", "tensor2_impl(X1, donothing, diffxmolecule, 'per')\n", "imshow(uint8(X1))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This image shows many artifacts since the pixel values lie outside the legal range: many of the intensities are in fact negative. Let us therefore normalize and map all intensities to $[0,1]$." ] }, { "cell_type": "code", "execution_count": 27, "metadata": { "collapsed": false }, "outputs": [], "source": [ "X2 = excerpt.copy()\n", "tensor2_impl(X2, donothing, diffxmolecule, 'per')\n", "map_to_01(X2)\n", "X2 *= 255\n", "imshow(uint8(X2))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The predominant color of this image is an average grey, i.e. an intensity of about 0.5. \n", "To get more detail in the image we therefore try to increase the contrast by applying the function $f_{50}$ from Example 8.5 to each intensity value. \n", "This does indeed show more detail." ] }, { "cell_type": "code", "execution_count": 28, "metadata": { "collapsed": false }, "outputs": [], "source": [ "X3 = excerpt.copy()\n", "tensor2_impl(X3, donothing, diffxmolecule, 'per')\n", "map_to_01(X3)\n", "X3 *= 255\n", "contrast_adjust0(X3,50)\n", "imshow(uint8(X3))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It is important to understand the colors in these images. We have computed the derivative in the $x$-direction, and for this image it turns out that the smallest and largest value of the derivative have about the same absolute value, but with opposite signs. \n", "A value of 0 in the left image in Figure 8.11 corresponds to no change in intensity between the two pixels, but these will be mapped to about $0.5$. \n", "In the second image the edges (where the largest values of the derivative appear) have been mapped to black and white, while points with small or no changes in intensity have been mapped to a middle grey-tone. \n", "The middle image thus tells us that large parts of the image have little variation in intensity. \n", "\n", "\n", "### **Partial derivative in $y$-direction**\n", "\n", "The partial derivative $\\partial P/\\partial y$ can be computed analogously to $\\partial P/\\partial x$,\n", "i.e. we apply the filter $-S=\\frac{1}{2}\\{-1,\\underline{0},1\\}$ to all columns of the image \n", "(alternatively, apply the tensor product $-S\\otimes I$ to the image), where $S$ is the filter we\n", "used in the $x$-direction. The positive direction of this axis in an image is opposite to the direction of the $y$-axis we use when plotting functions, \n", "and this explains the additional minus sign.\n", "We now obtain the molecule" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\frac{1}{2}\n", "\\begin{pmatrix}\n", "0 & 1 & 0\\\\ \n", "0 & \\underline{0} & 0\\\\ \n", "0 & -1 & 0\n", "\\end{pmatrix}.\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let us compare the partial derivatives in both directions.\n", "First we reproduce the partial derivative in $x$-direction." ] }, { "cell_type": "code", "execution_count": 29, "metadata": { "collapsed": false }, "outputs": [], "source": [ "excerpt=create_excerpt() \n", "X1 = excerpt.copy()\n", "tensor2_impl(X1, donothing, diffxmolecule, 'per')\n", "map_to_01(X1)\n", "X1 *= 255\n", "contrast_adjust0(X1, 50)\n", "imshow(uint8(X1))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Then we compute the partial derivative in $y$-direction." ] }, { "cell_type": "code", "execution_count": 30, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def diffymolecule(x, bd_mode):\n", " filter_impl(array([-1, 0, 1])/2., x, bd_mode)\n", " \n", "X2 = excerpt.copy()\n", "tensor2_impl(X2, diffymolecule, donothing, 'per')\n", "map_to_01(X2)\n", "X2 *= 255\n", "contrast_adjust0(X2, 50)\n", "imshow(uint8(X2))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The intensities have been normalized and the contrast enhanced by the function $f_{50}$ from Example 8.5.\n", "\n", "### **The gradient**\n", "\n", "The gradient of a scalar function in two variables is defined by the vector" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\nabla P=\\Biggl(\\frac{\\partial P}{\\partial x}, \\frac{\\partial P}{\\partial y}\\Biggr).\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The length of the gradient is" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\vert \\nabla P\\vert = \\sqrt{\\Biggl(\\frac{\\partial P}{\\partial x}\\Biggr)^2 + \\Biggl(\\frac{\\partial P}{\\partial y}\\Biggr)^2}.\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "When the two first derivatives have been computed it is thus simple to compute the length of the gradient. \n", "As for the first order derivatives, it is possible for the length of the gradient to be outside the legal range of values. \n", "Let us first compute and show the values of the gradient." ] }, { "cell_type": "code", "execution_count": 31, "metadata": { "collapsed": false }, "outputs": [], "source": [ "excerpt=create_excerpt()\n", "X1 = excerpt.copy()\n", "tensor2_impl(X1, donothing, diffxmolecule, 'per')\n", "X2 = excerpt.copy()\n", "tensor2_impl(X2, diffymolecule, donothing, 'per')\n", "X1 = X1**2 + X2**2\n", "imshow(uint8(X1))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Then we map the intensities to the legal range." ] }, { "cell_type": "code", "execution_count": 32, "metadata": { "collapsed": false }, "outputs": [], "source": [ "X2[:] = X1\n", "map_to_01(X2)\n", "X2 *= 255\n", "imshow(uint8(X2))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Then we use the function $g_{0.01}$ to adjust the contrast." ] }, { "cell_type": "code", "execution_count": 33, "metadata": { "collapsed": false }, "outputs": [], "source": [ "X3[:] = X2\n", "contrast_adjust0(X3, 50)\n", "imshow(uint8(X3))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The image of the gradient looks quite different from the images of the two partial derivatives. The reason is that the numbers that represent the length of \n", "the gradient are (square roots of) sums of squares of numbers. This means that the parts of the image that have virtually constant intensity \n", "(partial derivatives close to 0) are colored black. In the images of the partial derivatives these values ended up in the middle of the range of intensity values, with a final color of grey, since there were both positive and negative values.\n", "To enhance the contrast for this image we should thus do something different, like applying a one of the functions in the right plot of Figure 8.6. \n", "\n", "The gradient contains information about both derivatives and therefore emphasizes edges in all directions. It also gives a simpler image since the sign of \n", "the derivatives has been removed.\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "## Example 8.11: Second-order derivatives\n", "
\n", "\n", "\n", "To compute the three second order derivatives we can combine the two computational molecules which we already have described. For the mixed second order derivative we get\n", "$(I\\otimes S)((-S)\\otimes I)=-S\\otimes S$. For the last two second order derivative $\\frac{\\partial^2P}{\\partial x^2}$, $\\frac{\\partial^2P}{\\partial y^2}$, we can also use the\n", "three point approximation to the second derivative [[knut]](#knut)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\frac{\\partial P}{\\partial x^2}(i,j)\\approx p_{i,j+1}-2p_{i,j}+p_{i,j-1}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "(again we have set $h=1$). This gives a smaller molecule than if we combine the two molecules for order one differentiation \n", "(i.e. $(I\\otimes S)(I\\otimes S)=(I\\otimes S^2)$ and $((-S)\\otimes I)((-S)\\otimes I)=(S^2\\otimes I)$), since \n", "$S^2=\\frac{1}{2}\\{1,\\underline{0},-1\\}\\frac{1}{2}\\{1,\\underline{0},-1\\}=\\frac{1}{4}\\{1,0,\\underline{-2},0,1\\}$. \n", "The second order derivatives of an image $P$ can thus be computed by applying the computational molecules" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\begin{align*}\n", "\\frac{\\partial^2P}{\\partial x^2}:&\\qquad\n", "\\begin{pmatrix}\n", "0 & 0 & 0\\\\ \n", "1 & \\underline{-2} & 1\\\\ \n", "0 & 0 & 0\n", "\\end{pmatrix}, \n", "\\\\ \n", "\\frac{\\partial^2P}{\\partial y \\partial x}:&\\qquad\n", "\\frac{1}{4}\n", "\\begin{pmatrix}\n", "-1 & 0 & 1\\\\ \n", "0 & \\underline{0} & 0\\\\ \n", "1 & 0 & -1\n", "\\end{pmatrix}, \n", "\\\\ \n", "\\frac{\\partial^2P}{\\partial y^2}:&\\qquad\n", "\\begin{pmatrix}\n", "0 & 1 & 0\\\\ \n", "0 & \\underline{-2} & 0\\\\ \n", "0 & 1 & 0\n", "\\end{pmatrix}. \n", "\\end{align*}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "With these molecules it is quite easy to compute the second-order derivatives.\n", "First the $xx$-derivative." ] }, { "cell_type": "code", "execution_count": 34, "metadata": { "collapsed": false }, "outputs": [], "source": [ "excerpt=create_excerpt() \n", "X1 = excerpt.copy()\n", "tensor2_impl(X1, donothing, diffxmolecule, 'per')\n", "tensor2_impl(X1, donothing, diffxmolecule, 'per')\n", "map_to_01(X1)\n", "X1 *= 255\n", "contrast_adjust0(X1, 100)\n", "imshow(uint8(X1))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Then the $xy$-derivative." ] }, { "cell_type": "code", "execution_count": 35, "metadata": { "collapsed": false }, "outputs": [], "source": [ "X2 = excerpt.copy()\n", "tensor2_impl(X2, diffymolecule, diffxmolecule, 'per')\n", "map_to_01(X2)\n", "X2 *= 255\n", "contrast_adjust0(X2, 100)\n", "imshow(uint8(X2))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Then the $yy$-derivative." ] }, { "cell_type": "code", "execution_count": 36, "metadata": { "collapsed": false }, "outputs": [], "source": [ "X3 = excerpt.copy()\n", "tensor2_impl(X3, diffymolecule, donothing, 'per')\n", "tensor2_impl(X3, diffymolecule, donothing, 'per')\n", "map_to_01(X3)\n", "X3 *= 255\n", "contrast_adjust0(X3, 100)\n", "imshow(uint8(X3))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The computed derivatives were first normalized and then the contrast enhanced with the function $f_{100}$ from Example 8.5.\n", "\n", "As for the first derivatives, the $xx$-derivative seems to emphasize vertical edges and the $yy$-derivative horizontal edges. However, we also see that the second derivatives are\n", "more sensitive to noise in the image (the areas of grey are less uniform). The mixed derivative behaves a bit differently from the other two, and not surprisingly it seems to pick\n", "up both horizontal and vertical edges.\n", "\n", "This procedure can be generalized to higher order derivatives also. To apply $\\frac{\\partial^{k+l}P}{\\partial x^k\\partial y^l}$ to an image we can compute $S_l\\otimes S_k$ where\n", "$S_r$ corresponds to any point method for computing the $r$'th order derivative. We can also compute $(S^l)\\otimes(S^k)$, where we iterate the\n", "filter $S=\\frac{1}{2}\\{1,\\underline{0},-1\\}$ for the first derivative. As pointed out, this gives longer filters.\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "## Example 8.12: A chess pattern image\n", "
\n", "\n", "\n", "Let us also apply the molecules for differentiation to the chess pattern test image produced by the code" ] }, { "cell_type": "code", "execution_count": 37, "metadata": { "collapsed": false }, "outputs": [], "source": [ "N = 128;\n", "img=zeros((N, N))\n", "for x in range(N):\n", " for y in range(N):\n", " img[x,y] = 255*( (mod(x-1,64)>=32) == (mod(y-1,64)>=32) )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "and looks as follows" ] }, { "cell_type": "code", "execution_count": 38, "metadata": { "collapsed": false }, "outputs": [], "source": [ "X11 = img.copy()\n", "imshow(uint8(X11))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We first apply $S\\otimes I$." ] }, { "cell_type": "code", "execution_count": 39, "metadata": { "collapsed": false }, "outputs": [], "source": [ "X12 = img.copy()\n", "tensor2_impl(X12, diffymolecule, donothing, 'per')\n", "map_to_01(X12)\n", "X12 *= 255\n", "imshow(uint8(X12))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Then we apply $I\\otimes S$." ] }, { "cell_type": "code", "execution_count": 40, "metadata": { "collapsed": false }, "outputs": [], "source": [ "X13 = img.copy()\n", "tensor2_impl(X13, donothing, diffxmolecule, 'per')\n", "map_to_01(X13)\n", "X13 *= 255\n", "imshow(uint8(X13))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Then we apply $S\\otimes S$." ] }, { "cell_type": "code", "execution_count": 41, "metadata": { "collapsed": false }, "outputs": [], "source": [ "X21 = img.copy()\n", "tensor2_impl(X21, diffymolecule, diffxmolecule, 'per')\n", "map_to_01(X21)\n", "X21 *= 255\n", "imshow(uint8(X21))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Then we apply $I\\otimes S^2$." ] }, { "cell_type": "code", "execution_count": 42, "metadata": { "collapsed": false }, "outputs": [], "source": [ "X22 = img.copy()\n", "tensor2_impl(X22, donothing, diffxmolecule, 'per')\n", "tensor2_impl(X22, donothing, diffxmolecule, 'per')\n", "map_to_01(X22)\n", "X22 *= 255\n", "imshow(uint8(X22))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The we apply $S^2\\otimes I$." ] }, { "cell_type": "code", "execution_count": 43, "metadata": { "collapsed": false }, "outputs": [], "source": [ "X23 = img.copy()\n", "tensor2_impl(X23, diffymolecule, donothing, 'per')\n", "tensor2_impl(X23, diffymolecule, donothing, 'per')\n", "map_to_01(X23)\n", "X23 *= 255\n", "imshow(uint8(X23))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "These images make it even clearer that \n", "* $S\\otimes I$ detects horizontal edges, \n", "\n", "* $I\\otimes S$ detects vertical edges, \n", "\n", "* $S\\otimes S$ detects all points where abrupt changes appear in both directions. \n", "\n", "We also see that the second order partial derivative detects exactly the same edges as the first order\n", "partial derivative, and that the edges detected with $I\\otimes S^2$ are wider than the ones detected with $I\\otimes S$. The reason is that the filter $S^2$\n", "has more filter coefficients than $S$. \n", "\n", "Edges are detected with different colors, reflecting whether the difference between the neighboring pixels is positive or negative. \n", "Note also the additional edges at the first and last rows/edges in the images. The reason is that the filter $S$ is defined by assuming that the pixels repeat periodically (i.e. it is a\n", "circulant Toeplitz matrix). \n", "\n", "Defining a two-dimensional filter by filtering columns and then rows is not the only way we can define a two-dimensional filter. Another possible way is to let the $MN\\times\n", "MN$-matrix itself be a filter. Unfortunately, this is a bad way to define filtering, since there are some undesirable effects near the boundaries between rows:\n", "in the vector we form, the last element of one row is followed by the first element of the next row. These boundary effects are unfortunate when a filter is applied.\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "## Exercise 8.13: Generating images\n", "
\n", "\n", "\n", "Write code which calls the function `tensor2_impl` with appropriate filters to produce the following images:\n", "\n", "\n", "**a)**\n", "The right image in Figure 8.9.\n", "\n", "\n", "\n", "**Solution.**\n", "The following code can be used" ] }, { "cell_type": "code", "execution_count": 44, "metadata": { "collapsed": false }, "outputs": [], "source": [ "X3 = excerpt.copy()\n", "\n", "S2 = convolve((1/8.)*array([1, 3, 3, 1]),(1/8.)*array([1, 3, 3, 1]))\n", "\n", "def longmolecule(x, bd_mode):\n", " filter_impl(S2, x, bd_mode)\n", " \n", "tensor2_impl(X3, longmolecule, longmolecule, 'symm')\n", "imshow(uint8(X3))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "**b)**\n", "The right image in Figure 8.11.\n", "\n", "\n", "\n", "**Solution.**\n", "The following code can be used" ] }, { "cell_type": "code", "execution_count": 45, "metadata": { "collapsed": false }, "outputs": [], "source": [ "X3 = excerpt.copy()\n", "tensor2_impl(X3, donothing, diffxmolecule, 'per')\n", "map_to_01(X3)\n", "X3 *= 255\n", "contrast_adjust0(X3,50)\n", "imshow(uint8(X3))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "**c)**\n", "The images in Figure 8.13.\n", "\n", "\n", "\n", "**Solution.**\n", "The following code can be used" ] }, { "cell_type": "code", "execution_count": 46, "metadata": { "collapsed": false }, "outputs": [], "source": [ "excerpt=create_excerpt()\n", "X1 = excerpt.copy()\n", "tensor2_impl(X1, donothing, diffxmolecule, 'per')\n", "X2 = excerpt.copy()\n", "tensor2_impl(X2, diffymolecule, donothing, 'per')\n", "X1 = X1**2 + X2**2\n", "imshow(uint8(X1))" ] }, { "cell_type": "code", "execution_count": 47, "metadata": { "collapsed": false }, "outputs": [], "source": [ "X2[:] = X1\n", "map_to_01(X2)\n", "X2 *= 255\n", "imshow(uint8(X2))" ] }, { "cell_type": "code", "execution_count": 48, "metadata": { "collapsed": false }, "outputs": [], "source": [ "X3[:] = X2\n", "contrast_adjust0(X3, 50)\n", "imshow(uint8(X3))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "**d)**\n", "The images in Figure 8.12.\n", "\n", "\n", "\n", "**Solution.**\n", "The following code can be used" ] }, { "cell_type": "code", "execution_count": 49, "metadata": { "collapsed": false }, "outputs": [], "source": [ "excerpt=create_excerpt() \n", "X1 = excerpt.copy()\n", "tensor2_impl(X1, donothing, diffxmolecule, 'per')\n", "map_to_01(X1)\n", "X1 *= 255\n", "contrast_adjust0(X1, 50)\n", "imshow(uint8(X1))" ] }, { "cell_type": "code", "execution_count": 50, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def diffymolecule(x, bd_mode):\n", " filter_impl(array([-1, 0, 1])/2., x, bd_mode)\n", " \n", "X2 = excerpt.copy()\n", "tensor2_impl(X2, diffymolecule, donothing, 'per')\n", "map_to_01(X2)\n", "X2 *= 255\n", "contrast_adjust0(X2, 50)\n", "imshow(uint8(X2))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "**e)**\n", "The images in Figure 8.14.\n", "\n", "\n", "\n", "**Solution.**\n", "The following code can be used" ] }, { "cell_type": "code", "execution_count": 51, "metadata": { "collapsed": false }, "outputs": [], "source": [ "excerpt=create_excerpt() \n", "X1 = excerpt.copy()\n", "tensor2_impl(X1, donothing, diffxmolecule, 'per')\n", "tensor2_impl(X1, donothing, diffxmolecule, 'per')\n", "map_to_01(X1)\n", "X1 *= 255\n", "contrast_adjust0(X1, 100)\n", "imshow(uint8(X1))" ] }, { "cell_type": "code", "execution_count": 52, "metadata": { "collapsed": false }, "outputs": [], "source": [ "X2 = excerpt.copy()\n", "tensor2_impl(X2, diffymolecule, diffxmolecule, 'per')\n", "map_to_01(X2)\n", "X2 *= 255\n", "contrast_adjust0(X2, 100)\n", "imshow(uint8(X2))" ] }, { "cell_type": "code", "execution_count": 53, "metadata": { "collapsed": false }, "outputs": [], "source": [ "X3 = excerpt.copy()\n", "tensor2_impl(X3, diffymolecule, donothing, 'per')\n", "tensor2_impl(X3, diffymolecule, donothing, 'per')\n", "map_to_01(X3)\n", "X3 *= 255\n", "contrast_adjust0(X3, 100)\n", "imshow(uint8(X3))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "## Example 8.22: Change of coordinates with the DFT\n", "
\n", "\n", "\n", "The DFT was defined as a change of coordinates from the standard basis to the Fourier basis. Let us substitute the DFT and the IDFT for $S_1$, $S_2$.\n", "\n", "Modern image standards do typically not apply a change of coordinates to the entire image. \n", "Rather the image is split into smaller squares of appropriate size, called blocks, and a\n", "change of coordinates is performed independently on each block. In this example we will split the image into blocks of size $8\\times 8$. \n", "\n", "Recall that the DFT values express frequency components. The same applies for the 2D DFT and thus for images, but \n", "frequencies are now represented in two different directions. Let us make a parallel to Example 2.5 for the Lena image, i.e. \n", "we will view the image after a 2D DFT, followed by discarding DFT coefficients below a given threshold, followed by a 2D IDFT.\n", "As for sound this should have little effect\n", "on the human perception of the image, if the threshold is chosen with care. \n", "DFT-coefficients in a matrix `X` below a threshold can be discarded with the following code:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " X *= (abs(X) >= threshold)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\\verb+abs(X)>=threshold+ returns a *threshold matrix* with $1$ and $0$ of the same size as `X`. \n", "\n", "The following code generates images using threshold values of 1200, 200, and 400, respectively." ] }, { "cell_type": "code", "execution_count": 54, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def dft_impl8(x, bd_mode):\n", " N = shape(x)[0]\n", " for n in range(0, N, 8):\n", " x[n:(n+8), :] = fft.fft(x[n:(n+8), :], axis=0)" ] }, { "cell_type": "code", "execution_count": 55, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def idft_impl8(x, bd_mode):\n", " N = shape(x)[0]\n", " for n in range(0, N, 8):\n", " x[n:(n+8), :] = fft.ifft(x[n:(n+8), :], axis=0)" ] }, { "cell_type": "code", "execution_count": 56, "metadata": { "collapsed": false }, "outputs": [], "source": [ "X = create_excerpt()\n", "X = X.astype(complex)\n", "X1 = X.copy()\n", "forw_comp_rev_2d(X1, dft_impl8, idft_impl8, 100)\n", "imshow(uint8(abs(X1)))" ] }, { "cell_type": "code", "execution_count": 57, "metadata": { "collapsed": false }, "outputs": [], "source": [ "X2 = X.copy()\n", "forw_comp_rev_2d(X2, dft_impl8, idft_impl8, 200)\n", "imshow(uint8(abs(X2)))" ] }, { "cell_type": "code", "execution_count": 58, "metadata": { "collapsed": false }, "outputs": [], "source": [ "X3 = X.copy()\n", "forw_comp_rev_2d(X3, dft_impl8, idft_impl8, 400)\n", "imshow(uint8(X3))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "When increasing the threshold, the image becomes more and more unclear, but the image is quite clear in \n", "the first case, where as much as more than $76.6\\%$ of the samples have been zeroed out. A blocking effect at the block boundaries is clearly visible.\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "## Example 8.23: Change of coordinates with the DCT\n", "
\n", "\n", "\n", "Similarly to the DFT, the DCT was the change of coordinates from the standard basis to what we called the DCT basis. Let us substitute the DCT and the IDCT for $S_1$, $S_2$.\n", "\n", "The DCT is used more than the DFT in image processing. In particular, the JPEG standard applies a two-dimensional DCT, rather than a two-dimensional DFT. \n", "With the JPEG standard, the blocks are always $8\\times 8$, as in the previous example. It is of course not a coincidence that a power of $2$ is chosen here: the DCT, as the\n", "DFT, has an efficient implementation for powers of $2$.\n", "\n", "If we follow the same strategy for the DCT as above, i.e. zero out DCT-coefficients below a given threshold[^jpegcomment], \n", "the following code generates images using threshold values of 30, 50, and 100, respectively." ] }, { "cell_type": "code", "execution_count": 59, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def dct_impl8(x, bd_mode):\n", " N = shape(x)[0]\n", " for n in range(0, N, 8):\n", " x[n:(n+8), :] = dct(x[n:(n+8), :], norm='ortho', axis=0)" ] }, { "cell_type": "code", "execution_count": 60, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def idct_impl8(x, bd_mode):\n", " N = shape(x)[0]\n", " for n in range(0, N, 8):\n", " x[n:(n+8), :] = idct(x[n:(n+8), :], norm='ortho', axis=0)" ] }, { "cell_type": "code", "execution_count": 61, "metadata": { "collapsed": false }, "outputs": [], "source": [ "X = create_excerpt()\n", "X1 = X.copy()\n", "forw_comp_rev_2d(X1, dct_impl8, idct_impl8, 30)\n", "imshow(uint8(X1))" ] }, { "cell_type": "code", "execution_count": 62, "metadata": { "collapsed": false }, "outputs": [], "source": [ "X2 = X.copy()\n", "forw_comp_rev_2d(X2, dct_impl8, idct_impl8, 50)\n", "imshow(uint8(X2))" ] }, { "cell_type": "code", "execution_count": 63, "metadata": { "collapsed": false }, "outputs": [], "source": [ "X3 = X.copy()\n", "forw_comp_rev_2d(X3, dct_impl8, idct_impl8, 100)\n", "imshow(uint8(X3))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It is also interesting to see what happens if we don't split the image into blocks. Of course, when we discard many of the DCT-coefficients, we\n", "should see some artifacts, but there is no reason to believe that these occur at the old block boundaries. \n", "The new artifacts can be seen \n", "with the following code, and take a completely different shape." ] }, { "cell_type": "code", "execution_count": 64, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def dct_impl_full(x, bd_mode):\n", " x[:] = dct(x, norm='ortho', axis=0)" ] }, { "cell_type": "code", "execution_count": 65, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def idct_impl_full(x, bd_mode):\n", " x[:] = idct(x, norm='ortho', axis=0)" ] }, { "cell_type": "code", "execution_count": 66, "metadata": { "collapsed": false }, "outputs": [], "source": [ "X1 = X.copy()\n", "forw_comp_rev_2d(X1, dct_impl_full, idct_impl_full, 30)\n", "imshow(uint8(X1))" ] }, { "cell_type": "code", "execution_count": 67, "metadata": { "collapsed": false }, "outputs": [], "source": [ "X2 = X.copy()\n", "forw_comp_rev_2d(X2, dct_impl_full, idct_impl_full, 50)\n", "imshow(uint8(X2))" ] }, { "cell_type": "code", "execution_count": 68, "metadata": { "collapsed": false }, "outputs": [], "source": [ "X3 = X.copy()\n", "forw_comp_rev_2d(X3, dct_impl_full, idct_impl_full, 100)\n", "imshow(uint8(X3))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the exercises you will be asked to write code which generates the images from these examples.\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "## Exercise 8.24: Implementing DFT and DCT on blocks\n", "
\n", "\n", "\n", "In this section we have used functions which apply the DCT and the DFT either to subblocks of size $8\\times 8$, or to the full image.\n", "\n", "\n", "**a)**\n", "Implement functions `dft_impl8`, `idft_impl8`, `dct_impl8`, and `idct_impl8` \n", "which apply the DFT, IDFT, DCT, and IDCT, to consecutive blocks of length $8$.\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "**b)**\n", "Implement the two-dimensional FFT, IFFT, DCT, and IDCT on images, with the help of their one-dimensional counterparts, as well as the function `tensor2_impl`.\n", "\n", "\n", "\n", "**Solution.**\n", "The following code can be used." ] }, { "cell_type": "code", "execution_count": 69, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def dft_impl_full(x, bd_mode):\n", " x[:] = fft.fft(x, axis=0)" ] }, { "cell_type": "code", "execution_count": 70, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def idft_impl_full(x, bd_mode):\n", " x[:] = fft.ifft(x, axis=0)" ] }, { "cell_type": "code", "execution_count": 71, "metadata": { "collapsed": false }, "outputs": [], "source": [ "tensor2_impl(X, dft_impl_full, dft_impl_full, 'symm')\n", "tensor2_impl(X, idft_impl_full, idft_impl_full, 'symm')" ] }, { "cell_type": "code", "execution_count": 72, "metadata": { "collapsed": false }, "outputs": [], "source": [ "tensor2_impl(X, dct_impl_full, dct_impl_full, 'symm')\n", "tensor2_impl(X, idct_impl_full, idct_impl_full, 'symm')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "**c)**\n", "The function `forw_comp_rev_2d` in the library applies given transforms to the rows and columns of an input image, \n", "sets the coefficients below a certain threshold to zero, and transforms back.\n", "Run `forw_comp_rev_2d` for different threshold parameters on the sample image, \n", "and with the functions `dct_impl8`, `idct_impl8`, as well as DCT and IDCT applied to the entire input, as parameters. \n", "Check that this reproduces the DCT test images of this section, \n", "and that the correct numbers of values which have been discarded (i.e. which are below the threshold) are printed on screen.\n", "\n", "\n", "\n", "**Solution.**\n", "The following code can be used." ] }, { "cell_type": "code", "execution_count": 73, "metadata": { "collapsed": false }, "outputs": [], "source": [ "X = create_excerpt()\n", "X1 = X.copy()\n", "forw_comp_rev_2d(X1, dct_impl8, idct_impl8, 30)\n", "imshow(uint8(X1))" ] }, { "cell_type": "code", "execution_count": 74, "metadata": { "collapsed": false }, "outputs": [], "source": [ "X2 = X.copy()\n", "forw_comp_rev_2d(X2, dct_impl8, idct_impl8, 50)\n", "imshow(uint8(X2))" ] }, { "cell_type": "code", "execution_count": 75, "metadata": { "collapsed": false }, "outputs": [], "source": [ "X3 = X.copy()\n", "forw_comp_rev_2d(X3, dct_impl8, idct_impl8, 100)\n", "imshow(uint8(X3))" ] }, { "cell_type": "code", "execution_count": 76, "metadata": { "collapsed": false }, "outputs": [], "source": [ "X1 = X.copy()\n", "forw_comp_rev_2d(X1, dct_impl_full, idct_impl_full, 30)\n", "imshow(uint8(X1))" ] }, { "cell_type": "code", "execution_count": 77, "metadata": { "collapsed": false }, "outputs": [], "source": [ "X2 = X.copy()\n", "forw_comp_rev_2d(X2, dct_impl_full, idct_impl_full, 50)\n", "imshow(uint8(X2))" ] }, { "cell_type": "code", "execution_count": 78, "metadata": { "collapsed": false }, "outputs": [], "source": [ "X3 = X.copy()\n", "forw_comp_rev_2d(X3, dct_impl_full, idct_impl_full, 100)\n", "imshow(uint8(X3))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "\n", "" ] } ], "metadata": {}, "nbformat": 4, "nbformat_minor": 2 }