In [25]:
%matplotlib inline
import matplotlib.pyplot as plt
import matplotlib as mpl
from PIL import Image, ImageDraw
import numpy as np
import math
from scipy import ndimage
from scipy.misc import imresize
from scipy import ndimage as nd
import scipy as sci
import scipy.signal as signal
In [26]:
#LOAD IMAGES
img = Image.open("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
In [27]:
def gaussKernel(sigma):
siz = 6*sigma+1
gauss = np.zeros([siz, siz])
for i in range(siz):
for j in range (siz):
gauss[i,j] = (1./(2*math.pi*sigma*sigma))*math.exp(-(math.pow((i-(float(siz)/2)),2)+math.pow((j-(float(siz)/2)),2))/(2*sigma*sigma))
return gauss
def filter_RGB_conv(img, filt):
filtered = np.zeros(img.shape)
try:
for i in range(0, 3):
filtered[:,:,i] = nd.convolve(img[:, :, i], filt, mode = 'constant', cval=0.0)
except:
filtered[:,:] = nd.convolve(img,filt,mode='constant',cval=0.0)
return filtered
def cropImage(mat, cent, area):
minx = max(cent[0]-(area[0]/2), 0)
miny = max(cent[1]-(area[1]/2), 0)
maxx = min(cent[0]+abs(minx-cent[0])+1, mat.shape[0]-1)
maxy = min(cent[1]+abs(cent[1]-miny)+1, mat.shape[1]-1)
return crop1d(mat, [minx,miny], [maxx,maxy])
def crop1d(im,start,end):
cr = im[start[0]:end[0]+1,start[1]:end[1]+1]
return np.array(cr)
########################
## 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]
def matToOdd(mat):
if (mat.shape[0]%2==0):
mat = crop1d(mat, [0,0], [mat.shape[0]-2,mat.shape[1]])
if (mat.shape[1]%2==0):
mat = crop1d(mat, [0,0], [mat.shape[0],mat.shape[1]-2])
return mat
In [28]:
#Deprecated functions, I leave it here just because I invested too much effort on them.
def getCorrelation2(temp, pat, shift):
##MAGIC COMES HERE!
st = np.array(shift)
end = np.array(shift) + (np.array(pat.shape))
corr = np.zeros([3,3])
for i in range(0,3):
for j in range(0,3):
st_t = st+np.array([i,j])-1
end_t = end+np.array([i,j])-1
tmp = crop1d(temp, st_t, end_t)
corr[i,j] = normalizedCrossCorrelation(pat,tmp)
return corr
def normalizedCrossCorrelation(m2,m1):
if (m1.shape != m2.shape):
print 'shape1', m1.shape, 'shape2', m2.shape
m1 = np.resize(m1, m2.shape)
mean1 = np.mean(m1)
mean2 = np.mean(m2)
std1 = np.std(m1)
std2 = np.std(m2)
xcorr = 0
n = (m1.shape[0]*m1.shape[1])
for x in xrange(m1.shape[0]):
for y in xrange(m1.shape[1]):
xcorr+= ((m1[x][y]*m2[x][y])-(mean1*mean2))
#xcorr*= (1./n)
xcorr/= (std1*std2)
return xcorr
This function returns the correlation given a template, a pattern and a shifting. It crops the template to fit the size of the pattern +- sampling/2
In [29]:
def getCorrelation(temp, pat, shift, sampling, show=False):
st = (getCenter(temp) + np.array(shift))-(getCenter(pat)) - sampling/2
end = (getCenter(temp) + np.array(shift)) + (getCenter(pat)) + sampling/2
t = crop1d(temp,st,end)
corr = signal.correlate2d(t, pat, mode='valid')
if show:
mod = mpl.rcParams['image.interpolation']
mpl.rcParams['image.interpolation'] = 'none'
plt.figure()
plt.subplot(131)
plt.imshow(pat)
plt.title('Pattern')
plt.axis('off')
plt.subplot(132)
plt.imshow(t)
plt.title('Template')
plt.axis('off')
plt.subplot(133)
plt.imshow(corr)
plt.title('Correlation')
plt.axis('off')
plt.gcf().set_size_inches((10,10))
mpl.rcParams['image.interpolation'] = mod
return corr
def getCenter(mat):
return np.array([mat.shape[0]/2,mat.shape[1]/2])
This function returns the shifting of a pattern considering a template. The shifting is calculated recursively.
There is an optional parameter called show that allows you to see the process.
In [30]:
#Parameter 'show' allows you to see all function steps
def pyramid_matching(template, pat, blurr, it, sampling, show=False):
# Change to odd numbers -- To have a valid center
template = matToOdd(template)
pat = matToOdd(pat)
if it == 0:
corr = signal.correlate2d(template, pat, mode='same')
maxAt = np.where(corr == np.max(corr))
maxAt = np.array([maxAt[0][0], maxAt[1][0]])
return maxAt-getCenter(corr)
template_samp = filter_RGB_conv(template, blurr)[::sampling, ::sampling]
pat_samp = filter_RGB_conv(pat, blurr)[::sampling, ::sampling]
#Get lower shifting
shift = sampling*pyramid_matching(template_samp, pat_samp, blurr,it-1, sampling, show)
#Calculate correlation
corr = getCorrelation(template,pat,shift, sampling, show)
#Get the max point of correlation
maxAt = np.where(corr == np.max(corr))
maxAt = np.array([maxAt[0][0], maxAt[1][0]])
#Add to shift
shift = np.array(shift+(maxAt-getCenter(corr)))
return np.array(shift)
In [31]:
#Define sobel filter
sobelx = np.array([[-1,0,1],[-2,0,2],[-1,0,1]])
sobely = np.transpose(sobelx)
filt_x = filter_RGB_conv(RGB, sobelx)
filt_y = filter_RGB_conv(RGB, sobely)
# 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())
In [32]:
filt_xy_R = crop1d(filt_xy[:,:,0], [140,100], [200,200])
patR = crop1d(filt_xy[:,:,0], [160,110], [185,160])
patG = crop1d(filt_xy[:,:,1], [160,110], [185,160])
patB = crop1d(filt_xy[:,:,2], [160,110], [185,160])
In [33]:
R = pyramid_matching(filt_xy_R, patR, gaussKernel(2),4,2)
G = pyramid_matching(filt_xy_R, patG, gaussKernel(2),4,2, show=True)
B = pyramid_matching(filt_xy_R, patB, gaussKernel(2),4,2)
slideR = R-R
slideG = G-R
slideB = B-R
In [34]:
#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 [36]:
# And show it!
mpl.rcParams['image.interpolation'] = 'gaussian'
plt.imshow(final.astype('uint8'))
plt.gcf().set_size_inches((10,10))
In [ ]: