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

# ##  Moebius transformation applied to a discrete image ##

#  The aim of this notebook is to illustrate how a Moebius transformation, defined by its action on three points, distorts a  discrete image.
# A Mobius transformation acts as a geometric transformation of an image, defined via the function
# `geometric_transform` from `scipy.ndimage.interpolation`.
# 
# In image processing, a geometric transformation of an image consists in modifying
# the positions of pixels in that image, but keeping
# their colors unchanged.
# 
# In order to define a geometric transformation, as it is implemented in `scipy.ndimage`, we fix some background:
# 
# A 3 or 4 channels image of resolution, $m\times n$, is interpreted as a mapping
# $$img:\{0,1, \ldots, m-1\}\times \{0,1, \ldots, n-1\}\to\mathbb{R}^3 (\mathbb{R}^4),$$
# defined by $img(k,l)=(r,g,b)$ or (r, g, b, a), i.e. to the pixel in the row k, column l, one 
# associates its color code (r, g, b) or (r, g, b, a).  
#     
# The  position of a pixel, $(row=k, col=l)$, is also given by its  coordinates $(x=l, y=k)$,
#  with respect     to the image system of axes, $Oxy$, where $O$ is the upper-left corner, $Ox$ points horizontally to the right,  and $Oy$ downwards.
# 
# The geometric transformation, as a mapping from image to a planar region, is  given in the world coordinate system, XO'Y,  with O' usually chosen as being  the lower left corner of the image.
# 
# We have the following relationship between a pixel world coordinates, $(X,Y)$, and image coordinates, $(x,y)$:
#     
# $$ \begin{array}{lll} 
# X&=&x\\
# Y&=&m-1-y
# \end{array}
# $$
# 
# Denoting by $D$ the rectangular region that covers  our image, 
# a geometric transformation  is an invertible map $T:D\to D'\subset \mathbb{R}^2$. The transformation $T$ is chosen such that  the transformed image has common regions with the original one, i.e. $T(D)\cap D\neq \emptyset$.
#     

#  The geometric transformation implemented in `scipy.ndimage`  works as [follows](https://docs.scipy.org/doc/scipy-0.18.1/reference/tutorial/ndimage.html): theoretically one defines the output image (i.e. the deformed image) as an image of the same resolution (or not) with the input image. For each point P in the output image  one evaluates $T^{-1}(P)$. 
#  - If $T^{-1}(P)$ is a pixel   in the input image,  one assigns to P the color of $T^{-1}(P)$.
#  - If $T^{-1}(P)$ is not just a  pixel ( a point of integer coordinates),  one assigns a color to P, through a spline interpolation of the colors of neighbouring pixels of $T^{-1}(P)$.  
#  - If $T^{-1}(P)$  is outside the rectangle that covers the input image, then $P$  is filled by a prescribed method. 
#  

# The Python function that implements $T^{-1}$, let us call it `inverse_map`,  has as a first mandatory parameter, a tuple of length equal to the output array (image) rank, and returns  a tuple of length equal to the input array (image) rank (recall that the rank of a `ndarray` is the length of its shape).

# The `scipy.ndimage.interpolation.geometric_transform` that applies the transformation $T$  to an image, `img`, is defined as follows:
# 
# `geometric_transform(img, inverse_mapping, output_shape=None,  output=None, order=3, mode='constant',
#                      cval=0.0, prefilter=True, extra_arguments=(),extra_keywords={})`    
# 
# - `order` sets the order of the spline interpolation;
# - `mode` is a key that sets the method of filling the points P, for which $T^{-1}(P)$ is not in img.
#    The options are: 'constant' (all such points are colored with the same color), 'nearest', 'reflect' or 'wrap';
#    default is 'constant'.
# - `cval` has effect when mode='constant'. It gives the  grey color code, between 0 and 255 (for jpg images), to fill the regions consisting in points that are not mapped by $T^{-1}$ in img.
# - `extra_arguments=()` is a tuple defining the arguments of inverse_mapping, other than the mandatory one, defined above;
# - for other keywords see [scipy docs](https://docs.scipy.org/doc/scipy-0.15.1/reference/generated/scipy.ndimage.interpolation.geometric_transform.html).

# Next we show how to apply a Moebius transform, $T(z)=(az+b)/(cz+d)$, $ad-bc\neq0$, to a color image.
# The inverse map is defined by  $T^{-1}(z)=(dz-b)/(-cz+a)$.

# In[1]:


import skimage.io as sio
import numpy as np
from scipy.ndimage import geometric_transform
import plotly.express as px
from plotly.offline import download_plotlyjs, init_notebook_mode,  iplot
init_notebook_mode(connected=True)


# Our color image has the shape $(m, n, 3)$. Hence the `inverse_map` has as a mandatory first parameter, a tuple of len(3):

# In[2]:


def inv_Moebius_transform(index, a, b, c, d, img): 
    #index[0] gives the row of a pixel in the  output image, 
    #index[1] the column, and index[2], the color channel
    #a, b, c, d are the Moebius transform coeffcients and img is the output image (ndarray)                                        
    
    z = index[1] + 1j*(img.shape[0]-1-index[0])#the complex number associated to a pixel (to the corresponding 
                                               # point expressed in coordinates X,Y)
    w = (d*z-b)/(-c*z+a) #T^{-1}(z)
    return     img.shape[0]-1-np.imag(w), np.real(w), index[2] #returns the "approx"  row, and column 
                                                               # for  T^{-1}(z)


# Read the image, img:

# In[3]:


img = sio.imread('https://raw.githubusercontent.com/empet/Math/master/Imags/lena_color.jpg')
fig = px.imshow(img)
fig.update_layout(width=500, height=500)
iplot(fig)


# Let us define the Moebius transformation that maps three points,  zp[i], from Lena's image, referenced to XO'Y, to three points, wp[i], i=0,1,2, referenced to the same system of coordinates:

# In[4]:


z1 = 300+1j*230
w1 = np.exp(1j*np.pi/11)*z1  #z1 is rotated counter-clockwise about O' with  pi/11 to get w1
zp = [z1, 20+250*1j, 400+1j*180 ]#zp[1] is a fixed point, zp[2] is translated vertically
wp = [w1, 20+250*1j, 400+1j*210]


# The coefficients, a, b, c, d, of the corresponding Moebius transformation are computed as
# [follows](https://en.wikipedia.org/wiki/M%C3%B6bius_transformation#Specifying_a_transformation_by_three_points):

# In[5]:


a = np.linalg.det([[zp[0]*wp[0], wp[0], 1], 
                   [zp[1]*wp[1], wp[1], 1], 
                   [zp[2]*wp[2], wp[2], 1]])

b = np.linalg.det([[zp[0]*wp[0], zp[0], wp[0]], 
                   [zp[1]*wp[1], zp[1], wp[1]], 
                   [zp[2]*wp[2], zp[2], wp[2]]])         

c = np.linalg.det([[zp[0], wp[0], 1], 
                   [zp[1], wp[1], 1], 
                   [zp[2], wp[2], 1]])

d = np.linalg.det([[zp[0]*wp[0], zp[0], 1], 
                   [zp[1]*wp[1], zp[1], 1], 
                   [zp[2]*wp[2], zp[2], 1]])


# Apply this transformation to the Lena's image:

# In[6]:


transformed_img = geometric_transform(img, inv_Moebius_transform, extra_arguments=(a, b, c, d, img), cval=240) 


# In[7]:


tr_fig = px.imshow(transformed_img)
tr_fig.update_layout(width=500, height=500)
iplot(tr_fig)


# Now let us apply the inverse, $T^{-1}$, to the initial image. Its inverse is just $T$:

# In[8]:


def inverse_map(index, a, b, c , d, img): 
    z = index[1] + 1j*(img.shape[0]-1-index[0]) 
    w = (a*z+b)/(c*z+d) 
    return  img.shape[0]-1-np.imag(w), np.real(w), index[2]   


# In[9]:


inv_transformed_img = geometric_transform(img, inverse_map, extra_arguments=(a, b,c, d, img), cval=240) 
tr2_fig =px.imshow(inv_transformed_img)
tr2_fig.update_layout(width=500, height=500)
iplot(tr2_fig)


# In[ ]: