# Basis for grayscale images

## Introduction

Consider the set of real-valued matrices of size $M\times N$; we can turn this into a vector space by defining addition and scalar multiplication in the usual way:

\begin{align}
\mathbf{A} + \mathbf{B} &=  
    \left[ 
        \begin{array}{ccc} 
            a_{0,0} & \dots & a_{0,N-1} \\ 
            \vdots & & \vdots \\ 
            a_{M-1,0} & \dots & b_{M-1,N-1} 
        \end{array}
    \right]
    + 
    \left[ 
        \begin{array}{ccc} 
            b_{0,0} & \dots & b_{0,N-1} \\ 
            \vdots & & \vdots \\ 
            b_{M-1,0} & \dots & b_{M-1,N-1} 
        \end{array}
    \right]
    \\
    &=
    \left[ 
        \begin{array}{ccc} 
            a_{0,0}+b_{0,0} & \dots & a_{0,N-1}+b_{0,N-1} \\ 
            \vdots & & \vdots \\ 
            a_{M-1,0}+b_{M-1,0} & \dots & a_{M-1,N-1}+b_{M-1,N-1} 
        \end{array}
    \right]     
    \\ \\ \\
\beta\mathbf{A} &=  
    \left[ 
        \begin{array}{ccc} 
            \beta a_{0,0} & \dots & \beta a_{0,N-1} \\ 
            \vdots & & \vdots \\ 
            \beta a_{M-1,0} & \dots & \beta a_{M-1,N-1}
        \end{array}
    \right]
\end{align}


As a matter of fact, the space of real-valued $M\times N$ matrices is completely equivalent to $\mathbb{R}^{MN}$ and we can always "unroll" a matrix into a vector. Assume we proceed column by column; then the matrix becomes

$$
    \mathbf{a} = \mathbf{A}[:] = [
        \begin{array}{ccccccc}
            a_{0,0} & \dots & a_{M-1,0} & a_{0,1} & \dots & a_{M-1,1} & \ldots & a_{0, N-1} & \dots & a_{M-1,N-1}
        \end{array}]^T
$$

Although the matrix and vector forms represent exactly the same data, the matrix form allows us to display the data in the form of an image. Assume each value in the matrix is a grayscale intensity, where zero is black and 255 is white; for example we can create a checkerboard pattern of any size with the following function:


```python
# usual python bookkeeping...
%pylab inline
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
import IPython
from IPython.display import Image
import math
from __future__ import print_function

# ensure all images will be grayscale
gray();
```

    Populating the interactive namespace from numpy and matplotlib



    <Figure size 432x288 with 0 Axes>


#### (?1) `gray()  # ensure all images will be grayscale`? What?

#### (R1)


```python
gray.__module__
```




    'matplotlib.pyplot'




```python
help(gray)
```

    Help on function gray in module matplotlib.pyplot:
    
    gray()
        Set the colormap to "gray".
        
        This changes the default colormap as well as the colormap of the current
        image if there is one. See ``help(colormaps)`` for more information.
    


#### (!1) Worth digging deep: `colormaps`


```python
help(colormaps)
```

    Help on function colormaps in module matplotlib.pyplot:
    
    colormaps()
        Matplotlib provides a number of colormaps, and others can be added using
        :func:`~matplotlib.cm.register_cmap`.  This function documents the built-in
        colormaps, and will also return a list of all registered colormaps if
        called.
        
        You can set the colormap for an image, pcolor, scatter, etc,
        using a keyword argument::
        
          imshow(X, cmap=cm.hot)
        
        or using the :func:`set_cmap` function::
        
          imshow(X)
          pyplot.set_cmap('hot')
          pyplot.set_cmap('jet')
        
        In interactive mode, :func:`set_cmap` will update the colormap post-hoc,
        allowing you to see which one works best for your data.
        
        All built-in colormaps can be reversed by appending ``_r``: For instance,
        ``gray_r`` is the reverse of ``gray``.
        
        There are several common color schemes used in visualization:
        
        Sequential schemes
          for unipolar data that progresses from low to high
        Diverging schemes
          for bipolar data that emphasizes positive or negative deviations from a
          central value
        Cyclic schemes
          for plotting values that wrap around at the endpoints, such as phase
          angle, wind direction, or time of day
        Qualitative schemes
          for nominal data that has no inherent ordering, where color is used
          only to distinguish categories
        
        Matplotlib ships with 4 perceptually uniform color maps which are
        the recommended color maps for sequential data:
        
          =========   ===================================================
          Colormap    Description
          =========   ===================================================
          inferno     perceptually uniform shades of black-red-yellow
          magma       perceptually uniform shades of black-red-white
          plasma      perceptually uniform shades of blue-red-yellow
          viridis     perceptually uniform shades of blue-green-yellow
          =========   ===================================================
        
        The following colormaps are based on the `ColorBrewer
        <https://colorbrewer2.org>`_ color specifications and designs developed by
        Cynthia Brewer:
        
        ColorBrewer Diverging (luminance is highest at the midpoint, and
        decreases towards differently-colored endpoints):
        
          ========  ===================================
          Colormap  Description
          ========  ===================================
          BrBG      brown, white, blue-green
          PiYG      pink, white, yellow-green
          PRGn      purple, white, green
          PuOr      orange, white, purple
          RdBu      red, white, blue
          RdGy      red, white, gray
          RdYlBu    red, yellow, blue
          RdYlGn    red, yellow, green
          Spectral  red, orange, yellow, green, blue
          ========  ===================================
        
        ColorBrewer Sequential (luminance decreases monotonically):
        
          ========  ====================================
          Colormap  Description
          ========  ====================================
          Blues     white to dark blue
          BuGn      white, light blue, dark green
          BuPu      white, light blue, dark purple
          GnBu      white, light green, dark blue
          Greens    white to dark green
          Greys     white to black (not linear)
          Oranges   white, orange, dark brown
          OrRd      white, orange, dark red
          PuBu      white, light purple, dark blue
          PuBuGn    white, light purple, dark green
          PuRd      white, light purple, dark red
          Purples   white to dark purple
          RdPu      white, pink, dark purple
          Reds      white to dark red
          YlGn      light yellow, dark green
          YlGnBu    light yellow, light green, dark blue
          YlOrBr    light yellow, orange, dark brown
          YlOrRd    light yellow, orange, dark red
          ========  ====================================
        
        ColorBrewer Qualitative:
        
        (For plotting nominal data, `.ListedColormap` is used,
        not `.LinearSegmentedColormap`.  Different sets of colors are
        recommended for different numbers of categories.)
        
        * Accent
        * Dark2
        * Paired
        * Pastel1
        * Pastel2
        * Set1
        * Set2
        * Set3
        
        A set of colormaps derived from those of the same name provided
        with Matlab are also included:
        
          =========   =======================================================
          Colormap    Description
          =========   =======================================================
          autumn      sequential linearly-increasing shades of red-orange-yellow
          bone        sequential increasing black-white color map with
                      a tinge of blue, to emulate X-ray film
          cool        linearly-decreasing shades of cyan-magenta
          copper      sequential increasing shades of black-copper
          flag        repetitive red-white-blue-black pattern (not cyclic at
                      endpoints)
          gray        sequential linearly-increasing black-to-white
                      grayscale
          hot         sequential black-red-yellow-white, to emulate blackbody
                      radiation from an object at increasing temperatures
          jet         a spectral map with dark endpoints, blue-cyan-yellow-red;
                      based on a fluid-jet simulation by NCSA [#]_
          pink        sequential increasing pastel black-pink-white, meant
                      for sepia tone colorization of photographs
          prism       repetitive red-yellow-green-blue-purple-...-green pattern
                      (not cyclic at endpoints)
          spring      linearly-increasing shades of magenta-yellow
          summer      sequential linearly-increasing shades of green-yellow
          winter      linearly-increasing shades of blue-green
          =========   =======================================================
        
        A set of palettes from the `Yorick scientific visualisation
        package <https://dhmunro.github.io/yorick-doc/>`_, an evolution of
        the GIST package, both by David H. Munro are included:
        
          ============  =======================================================
          Colormap      Description
          ============  =======================================================
          gist_earth    mapmaker's colors from dark blue deep ocean to green
                        lowlands to brown highlands to white mountains
          gist_heat     sequential increasing black-red-orange-white, to emulate
                        blackbody radiation from an iron bar as it grows hotter
          gist_ncar     pseudo-spectral black-blue-green-yellow-red-purple-white
                        colormap from National Center for Atmospheric
                        Research [#]_
          gist_rainbow  runs through the colors in spectral order from red to
                        violet at full saturation (like *hsv* but not cyclic)
          gist_stern    "Stern special" color table from Interactive Data
                        Language software
          ============  =======================================================
        
        A set of cyclic color maps:
        
          ================  =================================================
          Colormap          Description
          ================  =================================================
          hsv               red-yellow-green-cyan-blue-magenta-red, formed by
                            changing the hue component in the HSV color space
          twilight          perceptually uniform shades of
                            white-blue-black-red-white
          twilight_shifted  perceptually uniform shades of
                            black-blue-white-red-black
          ================  =================================================
        
        Other miscellaneous schemes:
        
          ============= =======================================================
          Colormap      Description
          ============= =======================================================
          afmhot        sequential black-orange-yellow-white blackbody
                        spectrum, commonly used in atomic force microscopy
          brg           blue-red-green
          bwr           diverging blue-white-red
          coolwarm      diverging blue-gray-red, meant to avoid issues with 3D
                        shading, color blindness, and ordering of colors [#]_
          CMRmap        "Default colormaps on color images often reproduce to
                        confusing grayscale images. The proposed colormap
                        maintains an aesthetically pleasing color image that
                        automatically reproduces to a monotonic grayscale with
                        discrete, quantifiable saturation levels." [#]_
          cubehelix     Unlike most other color schemes cubehelix was designed
                        by D.A. Green to be monotonically increasing in terms
                        of perceived brightness. Also, when printed on a black
                        and white postscript printer, the scheme results in a
                        greyscale with monotonically increasing brightness.
                        This color scheme is named cubehelix because the (r, g, b)
                        values produced can be visualised as a squashed helix
                        around the diagonal in the (r, g, b) color cube.
          gnuplot       gnuplot's traditional pm3d scheme
                        (black-blue-red-yellow)
          gnuplot2      sequential color printable as gray
                        (black-blue-violet-yellow-white)
          ocean         green-blue-white
          rainbow       spectral purple-blue-green-yellow-orange-red colormap
                        with diverging luminance
          seismic       diverging blue-white-red
          nipy_spectral black-purple-blue-green-yellow-red-white spectrum,
                        originally from the Neuroimaging in Python project
          terrain       mapmaker's colors, blue-green-yellow-brown-white,
                        originally from IGOR Pro
          turbo         Spectral map (purple-blue-green-yellow-orange-red) with
                        a bright center and darker endpoints. A smoother
                        alternative to jet.
          ============= =======================================================
        
        The following colormaps are redundant and may be removed in future
        versions.  It's recommended to use the names in the descriptions
        instead, which produce identical output:
        
          =========  =======================================================
          Colormap   Description
          =========  =======================================================
          gist_gray  identical to *gray*
          gist_yarg  identical to *gray_r*
          binary     identical to *gray_r*
          =========  =======================================================
        
        .. rubric:: Footnotes
        
        .. [#] Rainbow colormaps, ``jet`` in particular, are considered a poor
          choice for scientific visualization by many researchers: `Rainbow Color
          Map (Still) Considered Harmful
          <https://ieeexplore.ieee.org/document/4118486/?arnumber=4118486>`_
        
        .. [#] Resembles "BkBlAqGrYeOrReViWh200" from NCAR Command
          Language. See `Color Table Gallery
          <https://www.ncl.ucar.edu/Document/Graphics/color_table_gallery.shtml>`_
        
        .. [#] See `Diverging Color Maps for Scientific Visualization
          <http://www.kennethmoreland.com/color-maps/>`_ by Kenneth Moreland.
        
        .. [#] See `A Color Map for Effective Black-and-White Rendering of
          Color-Scale Images
          <https://www.mathworks.com/matlabcentral/fileexchange/2662-cmrmap-m>`_
          by Carey Rappaport
    



```python
# let's create a checkerboard pattern
SIZE = 4
img = np.zeros((SIZE, SIZE))
for n in range(0, SIZE):
    for m in range(0, SIZE):
        if (n & 0x1) ^ (m & 0x1):
            ## recall that ^ is the exclusive or
            #img[n, m] = 255
            img[m, n] = 255

# now display the matrix as an image
plt.matshow(img); 
```


```python
SIZE = 4
img = np.zeros((SIZE, SIZE))
img.dtype  # float64, better uint8 for images
```




    dtype('float64')




```python
np.uint8
```




    numpy.uint8



`m & 0x1` and `n & 0x1` tests the parity: if `m` is odd, then `m & 0x1` outputs `1` (or `True`); if `m` even,
output `0`. Using XOR after that guarantees that only when exactly one of `m` and `n` is odd will the case
be colored white. This implies that every time one crosses from one case to any of its neighboring cases,
the color changes


Note how the axes of `plt.matshow()` are arranged diff from that of `plt.imshow()`
help(plt.matshow)

```python
plt.imshow(img);
```

Given the equivalence between the space of $M\times N$ matrices and $\mathbb{R}^{MN}$ we can easily define the inner product between two matrices in the usual way:

$$
\langle \mathbf{A}, \mathbf{B} \rangle = \sum_{m=0}^{M-1} \sum_{n=0}^{N-1} a_{m,n} b_{m, n}
$$

(where we have neglected the conjugation since we'll only deal with real-valued matrices); in other words, we can take the inner product between two matrices as the standard inner product of their unrolled versions. The inner product allows us to define orthogonality between images and this is rather useful since we're going to explore a couple of bases for this space.

## Actual images

Conveniently, using IPython, we can read images from disk in any given format and convert them to numpy arrays; let's load and display for instance a JPEG image:
img = np.array(plt.imread('cameraman.jpg'), dtype=int)
plt.matshow(img);---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-9-37d88f4789f6> in <module>()
----> 1 img = np.array(plt.imread('cameraman.jpg'), dtype=int)
      2 plt.matshow(img);

/home/phunc20/.virtualenvs/dsp-py2.7/lib/python2.7/site-packages/matplotlib/pyplot.pyc in imread(*args, **kwargs)
   2371 @docstring.copy_dedent(_imread)
   2372 def imread(*args, **kwargs):
-> 2373     return _imread(*args, **kwargs)
   2374 
   2375 

/home/phunc20/.virtualenvs/dsp-py2.7/lib/python2.7/site-packages/matplotlib/image.pyc in imread(fname, format)
   1351             raise ValueError('Only know how to handle extensions: %s; '
   1352                              'with Pillow installed matplotlib can handle '
-> 1353                              'more images' % list(handlers))
   1354         return im
   1355 

ValueError: Only know how to handle extensions: [u'png']; with Pillow installed matplotlib can handle more images!pip install pillow

```python
img = np.array(plt.imread('cameraman.jpg'), dtype=int)
plt.matshow(img);
```

The image is a $64\times 64$ low-resolution version of the famous "cameraman" test picture. Out of curiosity, we can look at the first column of this image, which is is a $64×1$ vector:


```python
img[:,0]
```




    array([156, 157, 157, 152, 154, 155, 151, 157, 152, 155, 158, 159, 159,
           160, 160, 161, 155, 160, 161, 161, 164, 162, 160, 162, 158, 160,
           158, 157, 160, 160, 159, 158, 163, 162, 162, 157, 160, 114, 114,
           103,  88,  62, 109,  82, 108, 128, 138, 140, 136, 128, 122, 137,
           147, 114, 114, 144, 112, 115, 117, 131, 112, 141,  99,  97])



The values are integers between zero and 255, meaning that each pixel is encoded over 8 bits (or 256 gray levels).

## The canonical basis

The canonical basis for any matrix space $\mathbb{R}^{M\times N}$ is the set of "delta" matrices where only one element equals to one while all the others are 0. Let's call them $\mathbf{E}_n$ with $0 \leq n < MN$. Here is a function to create the canonical basis vector given its index:


```python
def canonical(n, M=5, N=10):
    e = np.zeros((M, N))
    #e[(n % M), int(n / M)] = 1
    e[n % M, n // M] = 1
    return e
```


```python
!python --version
```

    Python 3.8.2



```python
10 / 3, 10 // 3, int(10/3)
```




    (3.3333333333333335, 3, 3)



Here are some basis vectors: look for the position of white pixel, which differentiates them and note that we enumerate pixels column-wise:


```python
plt.matshow(canonical(0));
plt.matshow(canonical(1));
plt.matshow(canonical(49));
```

Note how diff `matshow()` is from `imshow()`:
> In a jupyter cell, successive `matshow()`'s can draw any number of images, while `imshow()` will only do it for the last one


```python
plt.imshow(canonical(0));
plt.imshow(canonical(1));
plt.imshow(canonical(49));
```

##### Stopped here (2020/11/19 12h25)

## Transmitting images

Suppose we want to transmit the "cameraman" image over a communication channel. The intuitive way to do so is to send the  pixel values one by one, which corresponds to sending the coefficients of the decomposition of the image over the canonical basis. So far, nothing complicated: to send the cameraman image, for instance, we will send $64\times 64 = 4096$ coefficients in a row. 

Now suppose that a communication failure takes place after the first half of the pixels have been sent. The received data will allow us to display an approximation of the original image only. If we replace the missing data with zeros, here is what we would see, which is not very pretty:


```python
# unrolling of the image for transmission (we go column by column, hence "F")
tx_img = np.ravel(img, "F")

# oops, we lose half the data
tx_img[int(len(tx_img)/2):] = 0

# rebuild matrix
rx_img = np.reshape(tx_img, (64, 64), "F")
plt.matshow(rx_img);
```
help(np.ravel)
Can we come up with a trasmission scheme that is more robust in the face of channel loss? Interestingly, the answer is yes, and it involves a different, more versatile basis for the space of images. What we will do is the following: 

* describe the Haar basis, a new basis for the image space
* project the image in the new basis
* transmit the projection coefficients
* rebuild the image using the basis vectors

We know a few things: if we choose an orthonormal basis, the analysis and synthesis formulas will be super easy (a simple inner product and a scalar multiplication respectively). The trick is to find a basis that will be robust to the loss of some coefficients. 

One such basis is the **Haar basis**. We cannot go into too many details in this notebook but, for the curious, a good starting point is [here](https://chengtsolin.wordpress.com/2015/04/15/real-time-2d-discrete-wavelet-transform-using-opengl-compute-shader/). (An even better starting point is [these two papers](http://grail.cs.washington.edu/projects/wavelets/article/).) Mathematical formulas aside, the Haar basis works by encoding the information in a *hierarchical* way: the first basis vectors encode the broad information and the higher coefficients encode the detail. Let's have a look. 

First of all, to keep things simple, we will remain in the space of square matrices whose size is a power of two. The code to generate the Haar basis matrices is the following: first we generate a 1D Haar vector and then we obtain the basis matrices by taking the outer product of all possible 1D vectors (don't worry if it's not clear, the results are what's important):


```python
def haar1D(n, SIZE):
    # check power of two
    if math.floor(math.log(SIZE) / math.log(2)) != math.log(SIZE) / math.log(2):
        print("Haar defined only for lengths that are a power of two")
        return None
    if n >= SIZE or n < 0:
        print("invalid Haar index")
        return None
    
    # zero basis vector
    if n == 0:
        return np.ones(SIZE)
    
    # express n >= 1 as 2^p + q with p as large as possible;
    # then k = SIZE/2^p is the length of the support
    # and s = qk is the shift
    p = math.floor(math.log(n) / math.log(2))
    pp = int(pow(2, p))
    k = SIZE / pp
    s = (n - pp) * k
    
    h = np.zeros(SIZE)
    h[int(s):int(s+k/2)] = 1
    h[int(s+k/2):int(s+k)] = -1
    # these are not normalized
    return h


def haar2D(n, SIZE=8):
    # get horizontal and vertical indices
    hr = haar1D(n % SIZE, SIZE)
    hv = haar1D(int(n / SIZE), SIZE)
    # 2D Haar basis matrix is separable, so we can
    # just take the column-row product
    H = np.outer(hr, hv)
    # np.outer() is just column vector times row vector,
    # the 1st arg being the col vec, the 2nd the row vec.
    H = H / math.sqrt(np.sum(H * H))
    # the previous line just divides H by its Frobenius norm
    # so that the returned value of haar2D() has norm 1.
    return H
```
help(np.outer)
First of all, let's look at a few basis matrices; note that the matrices have **both positive and negative values**, so that the value of **zero** will be represented as **gray**:


```python
plt.matshow?
```


```python
plt.matshow(haar2D(0));
plt.matshow(haar2D(1));
plt.matshow(haar2D(10));
plt.matshow(haar2D(63));
```


```python
np.unique(haar2D(0))
```




    array([0.125])




```python
np.unique(haar2D(1))
```




    array([-0.125,  0.125])




```python
np.unique(haar2D(63))
```




    array([-0.5,  0. ,  0.5])




```python
white = np.unique(haar2D(0))[0]
black = -white
plt.matshow(haar2D(0), vmax=white, vmin=black);
```

We can notice two key properties

* each basis matrix has positive and negative values in some symmetric pattern: this means that the basis matrix will implicitly compute the difference between image areas
* low-index basis matrices take differences between large areas, while high-index ones take differences in smaller **localized** areas of the image

We can immediately verify that the Haar matrices are orthogonal:


```python
# let's use an 8x8 space; there will be 64 basis vectors
# compute all possible inner product and only print the nonzero results
for m in range(0,64):
    for n in range(0,64):
        r = np.sum(haar2D(m, 8) * haar2D(n, 8))
        if r != 0:
            print("[%dx%d -> %f] " % (m, n, r), end="")
```

    [0x0 -> 1.000000] [1x1 -> 1.000000] [2x2 -> 1.000000] [3x3 -> 1.000000] [4x4 -> 1.000000] [5x5 -> 1.000000] [6x6 -> 1.000000] [7x7 -> 1.000000] [8x8 -> 1.000000] [9x9 -> 1.000000] [10x10 -> 1.000000] [11x11 -> 1.000000] [12x12 -> 1.000000] [13x13 -> 1.000000] [14x14 -> 1.000000] [15x15 -> 1.000000] [16x16 -> 1.000000] [16x17 -> -0.000000] [17x16 -> -0.000000] [17x17 -> 1.000000] [18x18 -> 1.000000] [19x19 -> 1.000000] [20x20 -> 1.000000] [21x21 -> 1.000000] [22x22 -> 1.000000] [23x23 -> 1.000000] [24x24 -> 1.000000] [24x25 -> -0.000000] [25x24 -> -0.000000] [25x25 -> 1.000000] [26x26 -> 1.000000] [27x27 -> 1.000000] [28x28 -> 1.000000] [29x29 -> 1.000000] [30x30 -> 1.000000] [31x31 -> 1.000000] [32x32 -> 1.000000] [33x33 -> 1.000000] [34x34 -> 1.000000] [35x35 -> 1.000000] [36x36 -> 1.000000] [37x37 -> 1.000000] [38x38 -> 1.000000] [39x39 -> 1.000000] [40x40 -> 1.000000] [41x41 -> 1.000000] [42x42 -> 1.000000] [43x43 -> 1.000000] [44x44 -> 1.000000] [45x45 -> 1.000000] [46x46 -> 1.000000] [47x47 -> 1.000000] [48x48 -> 1.000000] [49x49 -> 1.000000] [50x50 -> 1.000000] [51x51 -> 1.000000] [52x52 -> 1.000000] [53x53 -> 1.000000] [54x54 -> 1.000000] [55x55 -> 1.000000] [56x56 -> 1.000000] [57x57 -> 1.000000] [58x58 -> 1.000000] [59x59 -> 1.000000] [60x60 -> 1.000000] [61x61 -> 1.000000] [62x62 -> 1.000000] [63x63 -> 1.000000] 

OK! Everything's fine. Now let's transmit the "cameraman" image: first, let's verify that it works


```python
# project the image onto the Haar basis, obtaining a vector of 4096 coefficients
# this is simply the analysis formula for the vector space with an orthogonal basis
tx_img = np.zeros(64*64)
for k in range(0, (64*64)):
    tx_img[k] = np.sum(img * haar2D(k, 64))

# now rebuild the image with the synthesis formula; since the basis is orthonormal
# we just need to scale the basis matrices by the projection coefficients
rx_img = np.zeros((64, 64))
for k in range(0, (64*64)):
    rx_img += tx_img[k] * haar2D(k, 64)

plt.matshow(rx_img);
```


```python
np.linalg.norm(img - rx_img, inf)
```




    2.112088282046898e-12


help(np.linalg.norm)
Cool, it works! Now let's see what happens if we lose the second half of the coefficients:


```python
# oops, we lose half the data
lossy_img = np.copy(tx_img);
lossy_img[int(len(tx_img)/2):] = 0

# rebuild matrix
rx_img = np.zeros((64, 64))
for k in range(0, (64*64)):
    rx_img += lossy_img[k] * haar2D(k, 64)

plt.matshow(rx_img);
```

That's quite remarkable, no? We've lost the same amount of information as before but the image is still acceptable. This is because we lost the coefficients associated to the fine details of the image but we retained the "broad strokes" encoded by the first half. 

Note that if we lose the first half of the coefficients, the result would look remarkably different:


```python
lossy_img = np.copy(tx_img);
lossy_img[0:int(len(tx_img)/2)] = 0

rx_img = np.zeros((64, 64))
for k in range(0, (64*64)):
    rx_img += lossy_img[k] * haar2D(k, 64)

plt.matshow(rx_img);
```

In fact, schemes like this one are used in *progressive encoding*: send the most important information first and add details if the channel permits it. You may have experienced this while browsing the internet over a slow connection. 

All in all, a great application of a change of basis!

## A few of my own questions

**(?)** About `# check power of two` of `haar1D()`


```python
help(math.log)
```

    Help on built-in function log in module math:
    
    log(...)
        log(x[, base])
        
        Return the logarithm of x to the given base.
        If the base not specified, returns the natural logarithm (base e) of x.
    



```python
math.log(2)
```




    0.6931471805599453




```python
math.log(2, 2)
```




    1.0




```python
math.log(2, 10)
```




    0.30102999566398114



The author is just checking whether
$
\log_{2} \texttt{SIZE} = \frac{\ln \texttt{SIZE}}{\ln 2}
$
is an integer.

**N.B.** Unlike `math.log`, numpy has
- `np.log`: Natural logarithm
- `np.log2`: base 2
- `np.log10`: base 10
- `np.log1p`: log(1+p)

**(?)** What is `n` in `haar2D()`?
Must `n` be bounded?

**(R)**
From reading the code, it seems that if `SIZE = k`, then `n = 0, 1, 2, ..., k**2 -1`, like the above examples
- when `SIZE=8`, `n = 0, 1, ..., 63`
- when `SIZE=64`, `n = 0, 1, ..., 64**2 - 1`

The `n` in `haar2D()` has to do with the `n` in `haar1D()`. To better understand how the function is written that way, readers would better have read [one of the papers](http://grail.cs.washington.edu/projects/wavelets/article/wavelet1.pdf) mentioned above. Briefly speaking, `Haar1D(n , SIZE)` will return a basis for the vector space $V^j$, where with `SIZE`$ = 2^j$. For example, with the box basis and the Haar wavelets described in the paper,

$$
\forall\, f \in V^3,\; \text{we have}\\
f = 
c_{0}^{0} \phi_{0}^{0}
+ d_{0}^{0} \psi_{0}^{0}
+ \left(d_{0}^{1} \psi_{0}^{1} + d_{1}^{1} \psi_{1}^{1}\right)
+ \left(d_{0}^{2} \psi_{0}^{2} + d_{1}^{2} \psi_{1}^{2} + d_{2}^{2} \psi_{2}^{2} + d_{3}^{2} \psi_{3}^{2}\right).
$$
<br>
What `haar1D(n, SIZE)` does is that it returns those $\phi_{s}^{t}$ and $\psi_{s}^{t}$.<br> As you can see from the $V^3$ example, there should be $2^3$ such basis vectors/functions; in general, there will be `SIZE`$= 2^j$ basis vectors and that's why the subscript (here in the `haar1D` function `n`) runs from `0` to `SIZE-1`.

**(?)** Is it correct for `haar2D` to normalize by Frobenius norm?

**(R)**


## Ref.
- [http://grail.cs.washington.edu/projects/wavelets/](http://grail.cs.washington.edu/projects/wavelets/)
    - [http://grail.cs.washington.edu/Research/](http://grail.cs.washington.edu/Research/)
- **exercise3.4 p.57** of the textbook written by the same authors of the course also talks about Haar basis.    


```python
type(pow(2,3)), type(pow(2,3.0))
```




    (int, float)




```python
pow(2,3), pow(2,3.0)
```




    (8, 8.0)




```python
2**3, 2**3.0
```




    (8, 8.0)




```python

```

069c785a4dedfd845f1cebf03b3abe4464474103
ipynb
Jupyter Notebook
epfl/2020/hw-ipynb/HaarBasis/hb.ipynb
phunc20/dsp
e7c496eb5fd4b8694eab0fc049cf98a5e3dfd886
MIT
2021-03-12T18:32:06.000Z
2021-03-12T18:32:06.000Z
epfl/2020/hw-ipynb/HaarBasis/hb.ipynb
phunc20/dsp
e7c496eb5fd4b8694eab0fc049cf98a5e3dfd886
MIT
epfl/2020/hw-ipynb/HaarBasis/hb.ipynb
phunc20/dsp
e7c496eb5fd4b8694eab0fc049cf98a5e3dfd886
MIT
Qwen/Qwen-72B
1. YES
2. YES