%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.image as mpimg
import sys,os
ia898path = os.path.abspath('../../')
if ia898path not in sys.path:
    sys.path.append(ia898path)
import ia898.src as ia

"""
Explanation: Demo magnify
Illustrate the interpolation of magnified images
Description
In this demonstration, the interpolation of magnified image is explained using the frequency domain. The interpolation is a low pass filter that can be applied either in the spatial domain, which is the methods known as nearest neighbor or pixel replication and bi-linear interpolation. Or in the frequency domain, with an ideal filter or a butterworth filter.
End of explanation
"""


fin = mpimg.imread('../data/lenina.pgm')
nb = ia.nbshow(3)
#nb.nbshow(fin,'%s'% (fin.shape,))
froi = fin[137:137+64,157:157+64]
nb.nbshow(froi,'%s'% (froi.shape,))
nb.nbshow()
print(froi.shape)

"""
Explanation: Reading and ROI selection
The image is read and a 64x64 ROI is selected and displayed
End of explanation
"""


import numpy as np

fd = froi.astype(float)
F = np.fft.fft2(fd)
nb = ia.nbshow(2)
nb.nbshow(froi)
nb.nbshow(ia.dftview(F))
nb.nbshow()

"""
Explanation: DFT
The DFT of the small image is taken and its spectrum displayed
End of explanation
"""


fx4 = np.zeros(4*np.array(froi.shape), 'uint8')
fx4[::4,::4] = froi
nb.nbshow(froi)
nb.nbshow(fx4)
nb.nbshow()

"""
Explanation: Expansion by 4 without interpolation
The image is expanded by 4, but filling the new pixels with 0
End of explanation
"""


fdx4 = fx4.astype(np.float)
Fx4 = np.fft.fft2(fdx4)

nb.nbshow(fx4)
nb.nbshow(ia.dftview(Fx4))
nb.nbshow()

"""
Explanation: DFT of the expansion without interpolation
Using the expansion propertie of the DFT (only valid for the discrete case), the resulting DFT is a periodical replication of the original DFT.
End of explanation
"""


k = np.ones((4,4))
fx4nn = ia.pconv(fdx4, k)
nb.nbshow(fx4)
nb.nbshow(fx4nn.astype(np.uint8))
nb.nbshow()
print(fx4nn[:9,:9])

"""
Explanation: Filtering by mean filtering - nearest neighbor
Filtering the expanded image using an average filter of size 4x4 is equivalent of applying a nearest neighbor interpolator. The zero pixels are replaced by the nearest non-zero pixel. This is equivalent to interpolation by pixel replication.
End of explanation
"""


kzero = np.zeros(fx4.shape)
kzero[0:4,0:4] = k
K = np.fft.fft2(kzero)
nb.nbshow(ia.dftview(K))
Fx4nn = K * Fx4
nb.nbshow(ia.dftview(Fx4nn))
nb.nbshow()

nb.nbshow(ia.dftview(Fx4))
nb.nbshow(ia.dftview(Fx4nn))
nb.nbshow()

"""
Explanation: Interpretation of the mean filtering in the frequency domain
Filtering by the average filter in space domain is equivalent to filter in the frequency domain by the sync filter.
End of explanation
"""


klinear = np.array([[1,2,3,4,3,2,1]])/4.
k2dlinear = klinear.T.dot(klinear)
print('k2dlinear=\n',k2dlinear)
fx4li = ia.pconv(fdx4, k2dlinear)
nb.nbshow(fx4)
nb.nbshow(fx4li.astype(np.uint8))
nb.nbshow()

"""
Explanation: Filtering by pyramidal kernel, linear interpolation
Filtering by a pyramidal kernel in space domain is equivalent to make a bi-linear interpolation. The zero pixels are replaced by a weighted sum of the neighbor pixels, the weight is inversely proportional to the non-zero pixel distance.
End of explanation
"""


klizero = np.zeros(fx4.shape)
klizero[0:7,0:7] = k2dlinear
Klinear = np.fft.fft2(klizero)
nb.nbshow(ia.dftview(Klinear))
Fx4li = Klinear * Fx4
nb.nbshow(ia.dftview(Fx4li))
nb.nbshow()

nb.nbshow(ia.dftview(Fx4))
nb.nbshow(ia.dftview(Fx4li))
nb.nbshow()

"""
Explanation: Interpretation of the pyramid filtering in the frequency domain
Filtering by the pyramid filter in space domain is equivalent to filter in the frequency domain by the square of the sync filter.
End of explanation
"""


FI = np.zeros_like(Fx4)
FI[:FI.shape[0]//8,:FI.shape[1]//8] = 1
FI[-FI.shape[0]//8:,:FI.shape[1]//8] = 1
FI[:FI.shape[0]//8,-FI.shape[1]//8:] = 1
FI[-FI.shape[0]//8:,-FI.shape[1]//8:] = 1
Fx4ideal = Fx4 * FI
nb.nbshow(ia.dftview(Fx4))
nb.nbshow(ia.dftview(Fx4ideal))
nb.nbshow()

nb = ia.nbshow(2)
fx4ideal_c = np.fft.ifft2(Fx4ideal)
print(fx4ideal_c.imag.sum())
fx4ideal= 16*fx4ideal_c.real
print(fx4ideal.min(),fx4ideal.max())
print(froi.min(),froi.max())
nb.nbshow(ia.normalize(fx4ideal))
nb.nbshow()

nb = ia.nbshow(3)
#nb.nbshow(froi)
nb.nbshow(fx4li.astype(np.uint8))
nb.nbshow(ia.normalize(fx4ideal))
nb.nbshow()

"""
Explanation: Using an ideal filter
Filtering by cutoff period of 8
End of explanation
"""


