PART 1

Load & crop image


In [3]:
%matplotlib inline

import matplotlib as mpl
mpl.rcParams['image.interpolation'] = 'none'
import matplotlib.pyplot as plt
from PIL import Image, ImageDraw
import numpy as np
import math
from scipy import ndimage
from scipy.misc import imresize

#LOAD IMAGES
img = Image.open("./img/picture.png")

#WARNING SIZE AND SHAPE INDICES ARE UPSIDE DOWN
size1_size=int(round(img.size[1]/10)) 
size2_size=int(round(img.size[0]/10))
img = imresize(img, ( size1_size,size2_size),interp='bilinear').astype('float')
img_np = np.array(img)

x1 = 20; y1 = 20;
x2 = 25; y2 = 333;
x3 = 24; y3 = 647;
w = 336; h = 302;

im1 = img_np[y1-2:y1+h-2,x1-1:x1+w-1]; 
im2 = img_np[y2-2:y2+h-2,x2-1:x2+w-1];
im3 = img_np[y3-2:y3+h-2,x3-1:x3+w-1];

I1 = 256*im1.astype('double')/im1.max(); 
I2 = 256*im2.astype('double')/im2.max(); 
I3 = 256*im3.astype('double')/im3.max();

RGB=np.zeros([h,w,3],dtype=type(img_np[0,0]))
RGB[:,:,0]=I3 
RGB[:,:,1]=I2 
RGB[:,:,2]=I1

plt.imshow(RGB.astype('uint8'))


Out[3]:
<matplotlib.image.AxesImage at 0x1a172fd0>

Extract image edges


In [4]:
##########################
## Functions
##########################


from scipy import fftpack 
# Filter in frequency domain
def filter_RGB_fft(img, filt, fftsize):
    # Transform filt to freq
    filtered = np.zeros(img.shape)
    fil_fft = fftpack.fft2(filt, (fftsize, fftsize))
    
    #Iterate over RGB
    for i in range(0,3):
        im_fft = fftpack.fft2(img[:,:,i], (fftsize, fftsize))
        im_fil_fft = im_fft * fil_fft
        im_fil = fftpack.ifft2(im_fil_fft)
        SZ = 50
        hs=np.floor(SZ/2.)

        filtered[:,:,i] = im_fil[hs:img.shape[0]+hs, hs:img.shape[1]+hs]
    
    return filtered

In [17]:
#Define sobel filter
sobelx = np.array([[-1,0,1],[-2,0,2],[-1,0,1]])
sobely = np.transpose(sobelx)

fftsize = 512
filt_x = filter_RGB_fft(RGB, sobelx, fftsize)
filt_y = filter_RGB_fft(RGB, sobely, fftsize)
# Generate filtered on xy
filt_xy = np.sqrt((filt_x*filt_x)+(filt_y*filt_y))

# Normalize betwen 0 and 1.
filt_xy *= (1./filt_xy.max())

# Show edges on channel R
plt.imshow(filt_xy[:,:,0], cmap = 'gray')
plt.colorbar()


Out[17]:
<matplotlib.colorbar.Colorbar instance at 0x000000002307AB88>

Select interest region


In [6]:
def crop(im, start, end):
    cr = np.array([end[0]-start[0],end[1]-start[1],im.shape[2]])
    cr = im[start[0]:end[0],start[1]:end[1],:]
    return cr

In [15]:
template = crop(filt_xy,[129,85], [160,140])

fig = plt.figure()
plt.subplot(1,3,1)
plt.imshow(template[:,:,0], cmap = 'gray')
plt.axis('off')
plt.subplot(1,3,2)
plt.imshow(template[:,:,1], cmap = 'gray')
plt.axis('off')
plt.subplot(1,3,3)
plt.imshow(template[:,:,0], cmap = 'gray')
plt.axis('off')
plt.gcf().set_size_inches((15,15))


Apply xCorrelation


In [8]:
import scipy as sci
import scipy.signal as signal

corrR_R = signal.correlate2d(template[:,:,0], template[:,:,0], mode='same', boundary='fill', fillvalue=0)
corrR_G = signal.correlate2d(template[:,:,1], template[:,:,0], mode='same', boundary='fill', fillvalue=0)
corrR_B = signal.correlate2d(template[:,:,2], template[:,:,0], mode='same', boundary='fill', fillvalue=0)

In [9]:
# Show figures
plt.subplot(1,3,1)
plt.imshow(corrR_R)
plt.axis('off')
plt.subplot(1,3,2)
plt.imshow(corrR_G)
plt.axis('off')
plt.subplot(1,3,3)
plt.imshow(corrR_B)
plt.axis('off')
plt.gcf().set_size_inches((15,15))


Get the shifting of the channels


In [10]:
#Get the position of the max
Rat = np.array(np.where(corrR_R == np.max(corrR_R)))
Gat = np.array(np.where(corrR_G == np.max(corrR_G)))
Bat = np.array(np.where(corrR_B == np.max(corrR_B)))

Rat = np.array([Rat[0][0],Rat[1][0]])
Gat = np.array([Gat[0][0],Gat[1][0]])
Bat = np.array([Bat[0][0],Bat[1][0]])

slideG = Rat-Gat
slideB = Rat-Bat

print '-> Assume channel R is our center'
print '-> Channel B is shifted',slideG
print '-> Channel B is shifted',slideB


-> Assume channel R is our center
-> Channel B is shifted [ 0 -1]
-> Channel B is shifted [ 3 -8]

Slide channels to fit the image


In [11]:
########################
## Function ShiftMatrix
########################

def shiftMatrix(mat, shx, shy):
    catx = np.zeros([abs(shx), mat.shape[1]])
    caty = np.zeros([mat.shape[0]+abs(shx), abs(shy)])
    shifted = mat
    if shx < 0:
        shifted = np.concatenate((shifted, catx), axis=0)
        shx = abs(shx)
    elif shx > 0:
        shifted = np.concatenate((catx, shifted), axis=0)
        shx = 0
    if shy < 0:
        shifted = np.concatenate((shifted, caty), axis=1)
        shy = abs(shy)
    elif shy > 0:
        shifted = np.concatenate((caty, shifted), axis=1)
        shy = 0
    return shifted[shx:mat.shape[0]+shx,shy:mat.shape[1]+shy]

In [12]:
#Generate final image
final = np.zeros(RGB.shape)
final[:,:,0] = RGB[:,:,0]
final[:,:,1] = shiftMatrix(RGB[:,:,1], slideG[0], slideG[1])
final[:,:,2] = shiftMatrix(RGB[:,:,2], slideB[0], slideB[1])

In [13]:
mpl.rcParams['image.interpolation'] = 'gaussian'
plt.imshow(final.astype('uint8'))
plt.gcf().set_size_inches((10,10))