{
"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",
"# Using tensor products to apply wavelets to images\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"## Example 9.1: Piecewise constant functions\n",
"\n",
"\n",
"\n",
"If $V_m$ is a resolution space of piecewise constant functions, \n",
"the resolution space $V_m\\otimes V_m$ consists of functions in two variables which are constant on any square of the form $[k_12^{-m},(k_1+1)2^{-m})\\times[k_22^{-m},(k_2+1)2^{-m})$: \n",
"Clearly $\\phi_{m,k_1}\\otimes\\phi_{m,k_2}$ is constant on such a square and $0$ elsewhere, and these functions form a basis for $V_m\\otimes V_m$. \n",
"\n",
"Let us compute the orthogonal projection of $\\phi_{1,k_1}\\otimes\\phi_{1,k_2}$ onto $V_0\\otimes V_0$. Since the Haar wavelet is orthonormal, the basis functions \n",
"in (9.4) are orthonormal, so that we can use the orthogonal decomposition formula to find this projection. \n",
"Clearly $\\phi_{1,k_1}\\otimes\\phi_{1,k_2}$ has different support from all except\n",
"one of $\\phi_{0,n_1}\\otimes\\phi_{0,n_2}$. Since"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\langle \\phi_{1,k_1}\\otimes\\phi_{1,k_2},\\phi_{0,n_1}\\otimes\\phi_{0,n_2}\\rangle\n",
"=\n",
"\\langle \\phi_{1,k_1},\\phi_{0,n_1}\\rangle \\langle \\phi_{1,k_2},\\phi_{0,n_2}\\rangle\n",
"=\n",
"\\frac{\\sqrt{2}}{2}\\frac{\\sqrt{2}}{2} = \\frac{1}{2}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"when the supports intersect, we obtain"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\text{proj}_{V_0\\otimes V_0} (\\phi_{1,k_1}\\otimes\\phi_{1,k_2}) =\n",
" \\begin{cases}\n",
" \\frac{1}{2}(\\phi_{0,k_1/2}\\otimes\\phi_{0,k_2/2}) & \\text{$k_1,k_2$ even} \\\\ \n",
" \\frac{1}{2}(\\phi_{0,k_1/2}\\otimes\\phi_{0,(k_2-1)/2}) & \\text{$k_1$ even, $k_2$ odd} \\\\ \n",
" \\frac{1}{2}(\\phi_{0,(k_1-1)/2}\\otimes\\phi_{0,k_2/2}) & \\text{$k_1$ odd, $k_2$ even} \\\\ \n",
" \\frac{1}{2}(\\phi_{0,(k_1-1)/2}\\otimes\\phi_{0,(k_2-1)/2}) & \\text{$k_1,k_2$ odd}\n",
" \\end{cases}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"So, in this case there are $4$ different formulas. Let us also compute the projection onto the orthogonal complement of $V_0\\otimes V_0$ in $V_1\\otimes V_1$. \n",
"Also here there are $4$ different formulas. When $k_1,k_2$ are both even we obtain"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\begin{align*} \n",
" \\lefteqn{\\phi_{1,k_1}\\otimes\\phi_{1,k_2} - \\text{proj}_{V_0\\otimes V_0}(\\phi_{1,k_1}\\otimes\\phi_{1,k_2})} \\\\ \n",
" &= \\phi_{1,k_1}\\otimes\\phi_{1,k_2} - \\frac{1}{2}(\\phi_{0,k_1/2}\\otimes\\phi_{0,k_2/2}) \\\\ \n",
" &= \\left(\\frac{1}{\\sqrt{2}}(\\phi_{0,k_1/2}+\\psi_{0,k_1/2})\\right)\\otimes\\left(\\frac{1}{\\sqrt{2}}(\\phi_{0,k_2/2}+\\psi_{0,k_2/2})\\right) \n",
" - \\frac{1}{2}(\\phi_{0,k_1/2}\\otimes\\phi_{0,k_2/2}) \\\\ \n",
" &= \\frac{1}{2}(\\phi_{0,k_1/2}\\otimes\\phi_{0,k_2/2}) + \\frac{1}{2}(\\phi_{0,k_1/2}\\otimes\\psi_{0,k_2/2}) \\\\ \n",
" & + \\frac{1}{2}(\\psi_{0,k_1/2}\\otimes\\phi_{0,k_2/2}) \n",
" + \\frac{1}{2}(\\psi_{0,k_1/2}\\otimes\\psi_{0,k_2/2}) - \\frac{1}{2}(\\phi_{0,k_1/2}\\otimes\\phi_{0,k_2/2}) \\\\ \n",
" &= \\frac{1}{2}(\\phi_{0,k_1/2}\\otimes\\psi_{0,k_2/2}) + \\frac{1}{2}(\\psi_{0,k_1/2}\\otimes\\phi_{0,k_2/2}) + \\frac{1}{2}(\\psi_{0,k_1/2}\\otimes\\psi_{0,k_2/2}).\n",
"\\end{align*}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here we have used the relation $\\phi_{1,k_i}=\\frac{1}{\\sqrt{2}}(\\phi_{0,k_i/2}+\\psi_{0,k_i/2})$, which we have previously obtained. \n",
"When either $k_1$ or $k_2$ is odd, similar formulas for the projection onto the orthogonal complement can be found. \n",
"In all cases, the formulas use the basis functions for $W_0^{(0,1)}$, $W_0^{(1,0)}$, $W_0^{(1,1)}$."
]
},
{
"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": "markdown",
"metadata": {},
"source": [
"Let us plot these."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"%matplotlib inline\n",
"\n",
"from numpy import *\n",
"import matplotlib.pyplot as plt\n",
"from mpl_toolkits.mplot3d import Axes3D\n",
"from images import *\n",
"from dwt import *\n",
"from forward_compress_reverse import *"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"First we define $\\phi$ and $\\psi$."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"t1 = linspace(-0.4, 1.4, 300)\n",
"t2 = linspace(-0.4, 1.4, 300)\n",
"x, y = meshgrid(t1,t2)\n",
"\n",
"phidef = lambda t: (t >= 0).astype(float)*(t <= 1).astype(float);\n",
"psidef = lambda t: (t >= 0).astype(float)*(t <= 0.5).astype(float) - (t >= 0.5).astype(float)*(t <= 1).astype(float);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Then $\\phi\\otimes\\phi$."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"fig = plt.figure()\n",
"ax = fig.gca(projection='3d')\n",
"ax.plot_surface(x, y, phidef(x)*phidef(y))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Then $\\phi\\otimes\\psi$."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"fig = plt.figure()\n",
"ax = fig.gca(projection='3d')\n",
"ax.plot_surface(x, y, phidef(x)*psidef(y))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Then $\\psi\\otimes\\phi$."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"fig = plt.figure()\n",
"ax = fig.gca(projection='3d')\n",
"ax.plot_surface(x, y, psidef(x)*phidef(y));"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Then $\\psi\\otimes\\psi$."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"fig = plt.figure()\n",
"ax = fig.gca(projection='3d')\n",
"ax.plot_surface(x, y, psidef(x)*psidef(y))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"## Example 9.2: Piecewise linear functions\n",
"\n",
"\n",
"\n",
"If we instead use any of the wavelets for piecewise linear functions, the wavelet basis functions are not orthogonal anymore, just as in the one-dimensional case.\n",
"Let us also plot these. First we define the new $\\phi$ and $\\psi$."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"t1 = linspace(-1, 2, 300)\n",
"t2 = linspace(-1, 2, 300)\n",
"x, y = meshgrid(t1, t2)\n",
"phidef = lambda t: (1-abs(t)).astype(float)*(abs(t)<=1).astype(float);\n",
"psidef = lambda t: (1-2*abs(t-0.5)).astype(float)*(t>=0).astype(float)*(t<=1).astype(float) - 0.25*phidef(t) - 0.25*phidef(t-1);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Then $\\phi\\otimes\\phi$."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"fig = plt.figure()\n",
"ax = fig.gca(projection='3d')\n",
"ax.plot_surface(x, y, phidef(x)*phidef(y))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Then $\\phi\\otimes\\psi$."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"fig = plt.figure()\n",
"ax = fig.gca(projection='3d')\n",
"ax.plot_surface(x, y, phidef(x)*psidef(y))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Then $\\psi\\otimes\\phi$."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"fig = plt.figure()\n",
"ax = fig.gca(projection='3d')\n",
"ax.plot_surface(x, y, psidef(x)*phidef(y))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Then $\\psi\\otimes\\psi$."
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"fig = plt.figure()\n",
"ax = fig.gca(projection='3d')\n",
"ax.plot_surface(x, y, psidef(x)*psidef(y))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"## Example 9.5: Applying the Haar wavelet to a chess pattern image\n",
"\n",
"\n",
"\n",
"Let us apply the Haar wavelet to the sample chess pattern example image from Figure 8.15, which was produced by the code"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"img=zeros((128,128))\n",
"for x in range(128):\n",
" for y in range(128):\n",
" img[x,y]=255*( (mod(x-1,64)>=32) == (mod(y-1,64)>=32) )"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The low-pass filter of the Haar wavelet was\n",
"essentially a low-pass filter with two coefficients. Also, as we have seen, the high-pass filter essentially computes an approximation to the partial derivative.\n",
"Clearly, abrupt changes in the vertical and horizontal directions appear here only at the edges in the chess pattern, and abrupt changes in both directions appear\n",
"only at the grid points in the chess pattern. After a 2D DWT Observation 8 thus states that we should see vertical edges from the chess pattern in one corner, \n",
"horizontal edges in another, and the grid pattern in yet another. \n",
"This can be verified with the following code,"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"dwt_impl(img, 'Haar', dims=2)\n",
"map_to_01(img)\n",
"img *= 255\n",
"imshow(uint8(img))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"## Example 9.6: Creating thumbnail images\n",
"\n",
"\n",
"\n",
"Let us apply the Haar wavelet to the Lena image. The following code computes the low-resolution approximation for $m=1$:"
]
},
{
"cell_type": "code",
"execution_count": 15,
"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": "code",
"execution_count": 16,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"X = create_excerpt()\n",
"dwt_impl(X, 'Haar')\n",
"X = X[0:(int(shape(X)[0]/2)), 0:(int(shape(X)[1]/2))]\n",
"map_to_01(X); X *= 255"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Also here it was necessary to map the result back to $[0,255]$. We also have omitted the third parameter here, \n",
"as $m=1$ is the default number of resolutions in `dwt_impl`. \n",
"Also, we have not entered the number of dimensions (2) in the parameter list in this case. \n",
"The reason is that this image, contrary to the previous one, is an RGB image: Since RGB images have data of dimension 3, the library \n",
"correctly assumes that the DWT should be parallelized on the third dimension (the color component), which is the desirable behavior. \n",
"\n",
"The following code can be used to get the result up to four resolutions."
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"img = create_excerpt()\n",
"l1, l2, l3 = shape(img)\n",
"X = 255*ones((l1/2,2*l2,l3))\n",
"newx = zeros_like(img)\n",
" \n",
"for m in range(1,5):\n",
" newx[:] = img\n",
" dwt_impl(newx, 'Haar', m=m)\n",
" X[(int(l1/4)-int(l1/2**(m+1))):(int(l1/4)+int(l1/2**(m+1))) , ((m-1)*int(l2/2) + int(l2/4)-int(l2/2**(m+1))):((m-1)*int(l2/2) + int(l2/4) + int(l2/2**(m+1))), :] = newx[:(int(l1/2**m)),:(int(l2/2**m)),:]/(2**m)\n",
"imshow(uint8(X))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The entire result after a one-level 2D DWT can be obtained as follows."
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"img = create_excerpt()\n",
"l1, l2, l3 = shape(img)\n",
" \n",
"X1 = img.copy()\n",
"dwt_impl(X1, 'Haar')\n",
"map_to_01(X1)\n",
"X1 *= 255\n",
"imshow(uint8(X1))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The entire result after a two-level 2D DWT can be obtained as follows."
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"X2 = img.copy()\n",
"dwt_impl(X2, 'Haar', m=2)\n",
"map_to_01(X2)\n",
"X2 *= 255\n",
"imshow(uint8(X2))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The first two thumbnail images can be seen as the upper left corners of the first two images. The other corners represent detail.\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"## Example 9.7: Detail and low-resolution approximations for different wavelets\n",
"\n",
"\n",
"\n",
"Let us take a closer look at the images generated when we use different wavelets, setting the detail coefficients to zero,\n",
"and viewing the result as an image of the same size. \n",
"In particular, let us discard the coefficients as pictured in Figure 9.4 and Figure 9.5.\n",
"We should expect that the lower order resolution approximations from $V_0$ are worse as $m$ increases.\n",
"\n",
"We will use the function `forw_comp_rev_dwt2` to generate the first four low-resolution approximations or detail for the usual excerpt of the sample image.\n",
"The low-level resolution approximations look as follows when we use the Haar wavelet."
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"X1 = img.copy()\n",
"forw_comp_rev_dwt2(X1, 1, 'Haar')\n",
"imshow(uint8(X1))"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"X2 = img.copy()\n",
"forw_comp_rev_dwt2(X2, 2, 'Haar')\n",
"imshow(uint8(X2))"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"X3 = img.copy()\n",
"forw_comp_rev_dwt2(X3, 3, 'Haar')\n",
"imshow(uint8(X3))"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"X4 = img.copy()\n",
"forw_comp_rev_dwt2(X4, 4, 'Haar')\n",
"imshow(uint8(X4))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The details look as follows when we use the Haar wavelet."
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"X1 = img.copy()\n",
"forw_comp_rev_dwt2(X1, 1, 'Haar', False)\n",
"imshow(uint8(X1))"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"X2 = img.copy()\n",
"forw_comp_rev_dwt2(X2, 2, 'Haar', False)\n",
"imshow(uint8(X2))"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"X3 = img.copy()\n",
"forw_comp_rev_dwt2(X3, 3, 'Haar', False)\n",
"imshow(uint8(X3))"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"X4 = img.copy()\n",
"forw_comp_rev_dwt2(X4, 4, 'Haar', False)\n",
"imshow(uint8(X4))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The low-level resolution approximations look as follows when we use the CDF 9/7 wavelet."
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"X1 = img.copy()\n",
"forw_comp_rev_dwt2(X1, 1, 'cdf97')\n",
"imshow(uint8(X1))"
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"X2 = img.copy()\n",
"forw_comp_rev_dwt2(X2, 2, 'cdf97')\n",
"imshow(uint8(X2))"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"X3 = img.copy()\n",
"forw_comp_rev_dwt2(X3, 3, 'cdf97')\n",
"imshow(uint8(X3))"
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"X4 = img.copy()\n",
"forw_comp_rev_dwt2(X4, 4, 'cdf97')\n",
"imshow(uint8(X4))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The details look as follows when we use the CDF 9/7 wavelet."
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"X1 = img.copy()\n",
"forw_comp_rev_dwt2(X1, 1, 'cdf97', False)\n",
"imshow(uint8(X1))"
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"X2 = img.copy()\n",
"forw_comp_rev_dwt2(X2, 2, 'cdf97', False)\n",
"imshow(uint8(X2))"
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"X3 = img.copy()\n",
"forw_comp_rev_dwt2(X3, 3, 'cdf97', False)\n",
"imshow(uint8(X3))"
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"X4 = img.copy()\n",
"forw_comp_rev_dwt2(X4, 4, 'cdf97', False)\n",
"imshow(uint8(X4))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Since black indicates values which are $0$, most of the coefficients must be small.\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"## Example 9.8: The Spline 5/3 wavelet and removing bands in the detail spaces\n",
"\n",
"\n",
"\n",
"The detail components in images are split into three bands, so let us see what happens when we discard some of them \n",
"(i.e. some of $W_m^{(1,1)}$, $W_m^{(1,0)}$, $W_m^{(0,1)}$). Let us use the Spline 5/3 wavelet. \n",
"\n",
"We will use the function `mmsubbands` to generate the low-resolution approximations obtained by zeroing out various numbers of the detail bands at level $m$."
]
},
{
"cell_type": "code",
"execution_count": 36,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"def mmsubbands(m):\n",
" img = create_excerpt()\n",
" l1, l2, l3 = shape(img)\n",
" X = zeros((3,l1,l2,l3))\n",
" \n",
" X[0] = img.copy()\n",
" dwt_impl(X[0], 'spline53', m=m)\n",
" X[0, (int(l1/2**(m-1))):, :, :] = 0\n",
" X[0, :(int(l1/2**(m-1))), (int(l2/2**(m-1))):, :] = 0\n",
" X[0, (int(l1/2**m)):(int(l1/2**(m-1))), (int(l2/2**m)):(int(l2/2**(m-1))), :] = 0\n",
" X[1] = X[0]\n",
" idwt_impl(X[0], 'spline53', m=m)\n",
" \n",
" X[1, (int(l1/2**m)):(int(l1/2**(m-1))), :(int(l2/2**m)), :] = 0\n",
" X[2] = X[1]\n",
" idwt_impl(X[1], 'spline53', m=m)\n",
" \n",
" X[2, :(int(l1/2**m)), (int(l2/2**m)):(int(l2/2**(m-1))), :] = 0\n",
" idwt_impl(X[2], 'spline53', m=m)\n",
" \n",
" for k in range(3):\n",
" map_to_01(X[k])\n",
" X[k] *= 255\n",
" return X[0], X[1], X[2]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The following code shows the resulting images when the bands on the first level indicated in Figure 9.4 are removed."
]
},
{
"cell_type": "code",
"execution_count": 37,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"X0, X1, X2 = mmsubbands(1)\n",
"imshow(uint8(X0))\n",
"imshow(uint8(X1))\n",
"imshow(uint8(X2))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let us also try this out at level 2."
]
},
{
"cell_type": "code",
"execution_count": 38,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"X0, X1, X2 = mmsubbands(2)\n",
"imshow(uint8(X0))\n",
"imshow(uint8(X1))\n",
"imshow(uint8(X2))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The image is seen still to resemble the original one, even after two levels of wavelet coefficients have been discarded. \n",
"This in itself is good for compression purposes, since we may achieve compression simply by dropping the given coefficients. \n",
"However, if we continue to discard more levels, the result will look poorer. \n",
"\n",
"With the following code we can discard also the third and fourth levels of detail."
]
},
{
"cell_type": "code",
"execution_count": 39,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"X3 = img.copy()\n",
"forw_comp_rev_dwt2(X3, 3, 'spline53')\n",
"imshow(uint8(X3))"
]
},
{
"cell_type": "code",
"execution_count": 40,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"X4 = img.copy()\n",
"forw_comp_rev_dwt2(X4, 4, 'spline53')\n",
"imshow(uint8(X4))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Although we still can see details in the image, the quality is definitely poorer.\n",
"\n",
"The following code displays the detail we obtain for DWT's with $1$, $2$, $3$, and $4$ levels."
]
},
{
"cell_type": "code",
"execution_count": 41,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"X1 = img.copy()\n",
"forw_comp_rev_dwt2(X1, 1, 'spline53', False)\n",
"imshow(uint8(X1))"
]
},
{
"cell_type": "code",
"execution_count": 42,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"X2 = img.copy()\n",
"forw_comp_rev_dwt2(X2, 2, 'spline53', False)\n",
"imshow(uint8(X2))"
]
},
{
"cell_type": "code",
"execution_count": 43,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"X3 = img.copy()\n",
"forw_comp_rev_dwt2(X3, 3, 'spline53', False)\n",
"imshow(uint8(X3))"
]
},
{
"cell_type": "code",
"execution_count": 44,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"X4 = img.copy()\n",
"forw_comp_rev_dwt2(X4, 4, 'spline53', False)\n",
"imshow(uint8(X4))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Clearly, more detail can be seen in the image when more of the detail is included.\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"## Exercise 9.10: Applying the Haar wavelet to another chess pattern image\n",
"\n",
"\n",
"\n",
"The following code produces a chess pattern type image almost identical to that from Example 9.5."
]
},
{
"cell_type": "code",
"execution_count": 45,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"img=zeros((128,128))\n",
"for x in range(128):\n",
" for y in range(128):\n",
" img[x,y]=255*( (mod(x,64)>=32) == (mod(y,64)>=32) )\n",
"imshow(uint8(img))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let us now apply a 2D DWT to this image as well with the Haar wavelet:"
]
},
{
"cell_type": "code",
"execution_count": 46,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"img2 = img.copy()\n",
"dwt_impl(img2, 'Haar', bd_mode='per', dims=2)\n",
"map_to_01(img2)\n",
"img2 *= 255\n",
"imshow(uint8(img2))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"There seem to be no detail components here, which is very different from \n",
"what you saw in Example 9.5, even though the images are very similar. \n",
"Attempt to explain what causes this to happen.\n",
"\n",
"\n",
"\n",
"**Hint.**\n",
"Compare with Exercise 4.14.\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"In Example 9.5, the borders in the chess pattern was chosen so that they occur at odd numbers. \n",
"This means that the image can not be represented\n",
"exactly in $V_{m-1}\\otimes V_{m-1}$, so that there is detail present in the image at all the borders in the chess pattern. \n",
"In the new image, the borders in the chess pattern was chosen so that they occur at even numbers. This means that the image can be represented\n",
"exactly in $V_{m-1}\\otimes V_{m-1}$, so that there is no detail components present in the image.\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"## Exercise 9.11: Implementing the fingerprint compression scheme\n",
"\n",
"\n",
"\n",
"Write code which generates the images shown in figures 9.20, \n",
"9.22, \n",
"and 9.23.\n",
"Use the functions `dwt_impl` and `idwt_impl` with the CDF 9/7 wavelet.\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"The following code generates the different parts of Figure 9.20."
]
},
{
"cell_type": "code",
"execution_count": 47,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"img=double(imread('images/fingerprint_500x500.png'))\n",
"imgnew=255*ones((512,512,3))\n",
"imgnew[0:500,0:500,:]=img\n",
"img=imgnew\n",
"l1, l2, l3 = 512, 512, 3\n",
" \n",
"X0 = img.copy()\n",
"dwt_impl(X0, 'cdf97')\n",
"X1 = X0.copy()\n",
"map_to_01(X0)\n",
"X0 *= 255\n",
"imshow(uint8(X0))"
]
},
{
"cell_type": "code",
"execution_count": 48,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"dwt_impl(X1[:(int(l1/2)), :(int(l2/2)), :], 'cdf97')\n",
"dwt_impl(X1[:(int(l1/2)), (int(l2/2)):, :], 'cdf97')\n",
"dwt_impl(X1[(int(l1/2)):, :(int(l2/2)), :], 'cdf97')\n",
"dwt_impl(X1[(int(l1/2)):, (int(l2/2)):512, :], 'cdf97')\n",
"X2 = X1.copy()\n",
"map_to_01(X1)\n",
"X1 *= 255\n",
"imshow(uint8(X1))"
]
},
{
"cell_type": "code",
"execution_count": 49,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"dwt_impl(X2[:(int(l1/4)), :(int(l2/4)), :], 'cdf97')\n",
"dwt_impl(X2[:(int(l1/4)), (int(l2/4)):(int(l2/2)), :], 'cdf97')\n",
"dwt_impl(X2[(int(l1/4)):(int(l1/2)), :(int(l2/4)), :], 'cdf97')\n",
"X3 = X2.copy()\n",
"map_to_01(X2)\n",
"X2 *= 255\n",
"imshow(uint8(X2))"
]
},
{
"cell_type": "code",
"execution_count": 50,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"dwt_impl(X3[:64, :64, :], 'cdf97')\n",
"dwt_impl(X3[:64, 64:128, :], 'cdf97')\n",
"dwt_impl(X3[:64, 128:192, :], 'cdf97')\n",
"dwt_impl(X3[:64, 192:256, :], 'cdf97')\n",
" \n",
"dwt_impl(X3[64:128, :64, :], 'cdf97')\n",
"dwt_impl(X3[64:128, 64:128, :], 'cdf97')\n",
"dwt_impl(X3[64:128, 128:192, :], 'cdf97')\n",
"dwt_impl(X3[64:128, 192:256, :], 'cdf97')\n",
" \n",
"dwt_impl(X3[128:192, :64, :], 'cdf97')\n",
"dwt_impl(X3[128:192, 64:128, :], 'cdf97')\n",
" \n",
"dwt_impl(X3[192:256, :64, :], 'cdf97')\n",
"dwt_impl(X3[192:256, 64:128, :], 'cdf97')\n",
"Z = X3.copy()\n",
"map_to_01(X3)\n",
"X3 *= 255\n",
"imshow(uint8(X3))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The following code generates Figure 9.22."
]
},
{
"cell_type": "code",
"execution_count": 51,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"dwt_impl(Z[:32,:32,:],'cdf97')\n",
"Y0 = Z.copy()\n",
"map_to_01(Z)\n",
"Z *= 255\n",
"Z = uint8(Z)\n",
"imshow(Z)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The following code generates the parts of Figure 9.23."
]
},
{
"cell_type": "code",
"execution_count": 52,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"Y0[16:, :, :] = 0\n",
"Y0[:16, 16:, :] = 0\n",
"idwt_impl(Y0, 'cdf97', m=5)\n",
"Y1 = Y0.copy()\n",
"map_to_01(Y0)\n",
"Y0 *= 255\n",
"imshow(uint8(Y0))"
]
},
{
"cell_type": "code",
"execution_count": 53,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"Y1 = img.copy()\n",
"dwt_impl(Y1, 'cdf97', m=5)\n",
"Y1[:16, :16, :] = 0\n",
"idwt_impl(Y1, 'cdf97', m=5)\n",
"map_to_01(Y1)\n",
"Y1 *= 255\n",
"imshow(uint8(Y1))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
""
]
}
],
"metadata": {},
"nbformat": 4,
"nbformat_minor": 2
}