#!/usr/bin/env python
# coding: utf-8

# In[ ]:


# Introduction to Image Processing


# This numerical tour explores some basic image processing tasks.

# *Important:* Please read the [installation page](http://gpeyre.github.io/numerical-tours/installation_python/) for details about how to install the toolboxes.
# $\newcommand{\dotp}[2]{\langle #1, #2 \rangle}$
# $\newcommand{\enscond}[2]{\lbrace #1, #2 \rbrace}$
# $\newcommand{\pd}[2]{ \frac{ \partial #1}{\partial #2} }$
# $\newcommand{\umin}[1]{\underset{#1}{\min}\;}$
# $\newcommand{\umax}[1]{\underset{#1}{\max}\;}$
# $\newcommand{\umin}[1]{\underset{#1}{\min}\;}$
# $\newcommand{\uargmin}[1]{\underset{#1}{argmin}\;}$
# $\newcommand{\norm}[1]{\|#1\|}$
# $\newcommand{\abs}[1]{\left|#1\right|}$
# $\newcommand{\choice}[1]{ \left\{  \begin{array}{l} #1 \end{array} \right. }$
# $\newcommand{\pa}[1]{\left(#1\right)}$
# $\newcommand{\diag}[1]{{diag}\left( #1 \right)}$
# $\newcommand{\qandq}{\quad\text{and}\quad}$
# $\newcommand{\qwhereq}{\quad\text{where}\quad}$
# $\newcommand{\qifq}{ \quad \text{if} \quad }$
# $\newcommand{\qarrq}{ \quad \Longrightarrow \quad }$
# $\newcommand{\ZZ}{\mathbb{Z}}$
# $\newcommand{\CC}{\mathbb{C}}$
# $\newcommand{\RR}{\mathbb{R}}$
# $\newcommand{\EE}{\mathbb{E}}$
# $\newcommand{\Zz}{\mathcal{Z}}$
# $\newcommand{\Ww}{\mathcal{W}}$
# $\newcommand{\Vv}{\mathcal{V}}$
# $\newcommand{\Nn}{\mathcal{N}}$
# $\newcommand{\NN}{\mathcal{N}}$
# $\newcommand{\Hh}{\mathcal{H}}$
# $\newcommand{\Bb}{\mathcal{B}}$
# $\newcommand{\Ee}{\mathcal{E}}$
# $\newcommand{\Cc}{\mathcal{C}}$
# $\newcommand{\Gg}{\mathcal{G}}$
# $\newcommand{\Ss}{\mathcal{S}}$
# $\newcommand{\Pp}{\mathcal{P}}$
# $\newcommand{\Ff}{\mathcal{F}}$
# $\newcommand{\Xx}{\mathcal{X}}$
# $\newcommand{\Mm}{\mathcal{M}}$
# $\newcommand{\Ii}{\mathcal{I}}$
# $\newcommand{\Dd}{\mathcal{D}}$
# $\newcommand{\Ll}{\mathcal{L}}$
# $\newcommand{\Tt}{\mathcal{T}}$
# $\newcommand{\si}{\sigma}$
# $\newcommand{\al}{\alpha}$
# $\newcommand{\la}{\lambda}$
# $\newcommand{\ga}{\gamma}$
# $\newcommand{\Ga}{\Gamma}$
# $\newcommand{\La}{\Lambda}$
# $\newcommand{\si}{\sigma}$
# $\newcommand{\Si}{\Sigma}$
# $\newcommand{\be}{\beta}$
# $\newcommand{\de}{\delta}$
# $\newcommand{\De}{\Delta}$
# $\newcommand{\phi}{\varphi}$
# $\newcommand{\th}{\theta}$
# $\newcommand{\om}{\omega}$
# $\newcommand{\Om}{\Omega}$
# 

# In[5]:


from nt_toolbox.general import *
from nt_toolbox.signal import *
get_ipython().run_line_magic('pylab', 'inline')
get_ipython().run_line_magic('matplotlib', 'inline')


# ## Image Loading and Displaying

# Several functions are implemented to load and display images.
# 
# 
# First we load an image.
# 
# path to the images

# In[18]:


name = 'nt_toolbox/data/hibiscus.bmp'
n = 256
M = load_image(name, n)


# We can display it. It is possible to zoom on it, extract pixels, etc.

# In[19]:


m = int(n/2)
imageplot(M[m-25:m+25,m-25:m+25], 'Zoom', [1, 2, 2])


# ## Image Modification

# An image is a 2D array, that can be modified as a matrix.

# In[20]:


imageplot(-M, '-M', [1,2,1])
imageplot(M[::-1,:], 'Flipped', [1,2,2])


# Blurring is achieved by computing a convolution with a kernel. 

# Compute the low pass Gaussian kernel. Warning, the indexes needs to be modulo $n$ in order to use FFTs.

# In[21]:


sigma = 5
t = concatenate( (arange(0,n/2+1), arange(-n/2,-1)) )
[Y,X] = np.meshgrid(t,t)
h = exp( -(X**2+Y**2)/(2.0*float(sigma)**2) )
h = h/sum(h)
imageplot(fftshift(h))


# Compute the periodic convolution ussing FFTs

# In[22]:


Mh = real( ifft2(fft2(M) * fft2(h)) )


# Display

# In[23]:


imageplot(M, 'Image', [1, 2, 1])
imageplot(Mh, 'Blurred', [1, 2, 2])


# Several differential and convolution operators are implemented.

# In[24]:


G = grad(M)
imageplot(G[:,:,0], 'd/ dx', [1, 2, 1])
imageplot(G[:,:,1], 'd/ dy', [1, 2, 2])


# ## Fourier Transform

# The 2D Fourier transform can be used to perform low pass approximation
# and interpolation (by zero padding).
# 
# Compute and display the Fourier transform (display over a log scale).
# The function fftshift is useful to put the 0 low frequency in the
# middle. After fftshift, the zero frequency is located at position
# $(n/2+1,n/2+1)$.

# In[25]:


Mf = fft2(M)
Lf = fftshift(log(abs(Mf) + 1e-1))
imageplot(M, 'Image', [1, 2, 1])
imageplot(Lf, 'Fourier transform', [1, 2, 2])


# **Exercise 1:** To avoid boundary artifacts and estimate really the frequency content of
# the image (and not of the artifacts!), one needs to multiply M by a
# smooth windowing function h and compute fft2(M*h). Use a sine
# windowing function. Can you interpret the resulting filter ?

# In[26]:


run -i nt_solutions/introduction_3_image/exo1


# **Exercise 2:** Perform low pass filtering by removing the high frequencies of the
# spectrum. What do you oberve ?

# In[28]:


run -i nt_solutions/introduction_3_image/exo2

<script>
  $(document).ready(function(){
      $('div.prompt').hide();
  });
</script>