Homework 3 : Harris Corner Detect.

Kyunghoon Kim ( kyunghoon@unist.ac.kr )

Mathematical Sciences, UNIST

2014.11.02.

Class : Computer Vision (Prof. Jae-Young Sim)

Environment



In [194]:

    
from PIL import Image # Load PIL(Python Image Library)
import matplotlib.pyplot as plt # Load pyplot library
import numpy as np # Load Numerical Library
import copy
import timeit

Test for korean letter 'biup'



In [2]:

    
im = Image.open('biup.jpg') # Load the image file
imarray = np.array(im) # Construct the array
print im.size



In [164]:

    
imgray = im.convert('L') # Convert the image to gray scale
plt.imshow(imgray, cmap=plt.get_cmap('gray'), interpolation='nearest') # interpolation='nearest' simply display the image without try to interpolate betwen pixels if the display resolution is not the same as the image resolution (which is most often the case). It will results in an image in which is image pixel is displayed as a square of multiple display pixels.

imgray_array = np.array(imgray) # Let imgray_array be the array of imgray.
print imgray_array # As we know, each pixels have a color value.









    



[[255 255 255 ..., 255 255 255]
 [255 255 255 ..., 255 255 255]
 [255 255 255 ..., 255 255 255]
 ..., 
 [255 255 255 ..., 255 255 255]
 [255 255 255 ..., 255 255 255]
 [255 255 255 ..., 255 255 255]]



In [165]:

    
plt.hist(imgray_array.flatten(), 256, [0,256]);

$c(x, y)=\sum_W [I(x_i , y_i)-I(x_i +\Delta x , y_i + \Delta y )]^2$

We need moving function($\Delta$) for index range of image such as (0, 100)



In [6]:

    
I = imgray_array
# 0 ~ 49
maxx = I.shape[0]-1
maxy = I.shape[1]-1

def moving(value, maxvalue=100):
    """ index from 0 to max """
    temp = copy.copy(value) # assign different address for keeping value
    for i in value:
        if i < 0:
            temp.pop(temp.index(i))
            temp.append(max(temp)+1)
        if i > maxvalue:
            temp.pop(temp.index(i))
            temp.append(min(min(temp), maxvalue+1)-1)
    temp.sort()
    return temp



In [7]:

    
# test
xwidth = range(38,38+5)
xwidth = moving(xwidth, 49)
ywidth = range(40,40+5)
ywidth = moving(ywidth, 49)
print xwidth









    



[38, 39, 40, 41, 42]

Now calculate the matrix A

$A=\begin{bmatrix}

   \sum_W (I_x (x_i , y_j )^2) & \sum_W I_x (x_i , y_j )I_y (x_i , y_j ) \\[0.3em]
   \sum_W I_x (x_i , y_j )I_y (x_i , y_j ) & \sum_W (I_y (x_i , y_j )^2)
 \end{bmatrix}$



In [8]:

    
x=20
y=20
# x=0
# y=0
size = 4 # have to be even
xwidth = moving(range(x-int(size/2), x+int(size/2)), maxx)
ywidth = moving(range(y-int(size/2), y+int(size/2)), maxy)

print xwidth
print ywidth









    



[18, 19, 20, 21]
[18, 19, 20, 21]



In [9]:

    
matrix = []
for i in xwidth:
    temp = []
    for j in ywidth:
        temp.append(I[i, j])
    matrix.append(temp)
print matrix









    



[[0, 68, 255, 255], [1, 68, 255, 255], [0, 68, 255, 254], [0, 68, 255, 254]]



In [42]:

    
def minus(row):
    length = len(row)
    a = row[:length-1]
    b = row[1:length]

#     temp = np.array([], dtype='int64')
#     index = 0
#     for i in a:
#         np.append(temp, [i - b[index]])
#         index += 1
    temp = np.array(b)-np.array(a)
    return list(temp)



In [53]:

    
for row in matrix:
    print row
    print minus(row)
    print [i**2 for i in minus(row)]









    



[0, 68, 255, 255]
[68, 187, 0]
[4624, 34969, 0]
[1, 68, 255, 255]
[67, 187, 0]
[4489, 34969, 0]
[0, 68, 255, 254]
[68, 187, 255]
[4624, 34969, 65025]
[0, 68, 255, 254]
[68, 187, 255]
[4624, 34969, 65025]



In [54]:

    
for row in matrix:
    gradient = minus(row)
    temp = [i**2 for i in gradient]
    print np.sum(temp)

FUNCTION TEST



In [105]:

    
def minus(row):
    length = len(row)
    a = row[:length-1]
    b = row[1:length]

#     temp = np.array([], dtype='int64')
#     index = 0
#     for i in a:
#         np.append(temp, [i - b[index]])
#         index += 1
    temp = np.array(b)-np.array(a)
    return list(temp)

def submatrix(x, y, size):
    xwidth = moving(range(x-int(size/2), x+int(size/2)), maxx)
    ywidth = moving(range(y-int(size/2), y+int(size/2)), maxy)

    matrix = []
    for i in xwidth:
        temp = []
        for j in ywidth:
            temp.append(I[i, j])
        matrix.append(temp)
    return matrix

def gradientx(matrix):
    result = 0
    for row in matrix:
        gradient = minus(row)
        temp = [i**2 for i in gradient]
        result += np.sum(temp)
    return result

def gradienty(matrix):
    result = 0
    matrix = np.array(matrix).transpose().tolist()
    for row in matrix:
        gradient = minus(row)
        temp = [i**2 for i in gradient]
        result += np.sum(temp)
    return result

def gradientxy(matrix):
    result = 0
    
    gradient_x = []
    for row in matrix:
        gradient_x.append(minus(row))
    matrix = np.array(matrix).transpose().tolist()
    gradient_y = []
    for row in matrix:
        gradient_y.append(minus(row))

    multiply = np.array(gradient_x)*np.array(gradient_y)
    for i in multiply:
        result += np.sum(i)
    return result



In [110]:

    
x=20
y=20
# x=0
# y=0
size = 4 # have to be even

matrix = submatrix(x, y, size)
print gradientx(matrix)
print gradienty(matrix)
print gradientxy(matrix)



In [111]:

    
A = np.zeros([2,2])
A[0,0] = gradientx(matrix)
A[0,1] = gradientxy(matrix)
A[1,0] = gradientxy(matrix)
A[1,1] = gradienty(matrix)
print A









    



[[  2.88287000e+05  -3.06000000e+02]
 [ -3.06000000e+02   3.00000000e+00]]



In [115]:

    
k=0.04
np.linalg.det(A)-k*np.trace(A)









    Out[115]:





759693.39999999979



In [116]:

    
np.linalg.eig(A)









    Out[116]:





(array([  2.88287325e+05,   2.67519566e+00]),
 array([[ 0.99999944,  0.00106145],
        [-0.00106145,  0.99999944]]))

FUNCTION

Window size is 4



In [169]:

    
size = 4 # have to be even
k=0.04

corner = np.zeros([maxx, maxy])
for x in range(maxx+1):
    for y in range(maxy+1):
        matrix = submatrix(x, y, size)
        A = np.zeros([2,2])
        A[0,0] = gradientx(matrix)
        A[0,1] = gradientxy(matrix)
        A[1,0] = gradientxy(matrix)
        A[1,1] = gradienty(matrix)
        R = np.linalg.det(A)-k*np.trace(A)
        if R>0:
            corner[x,y] = R



In [170]:

    
plt.imshow(corner, interpolation='nearest')









    Out[170]:





<matplotlib.image.AxesImage at 0x114207910>



In [171]:

    
plt.imshow(corner>1000000, interpolation='nearest')









    Out[171]:





<matplotlib.image.AxesImage at 0x114320990>

Window size is 2



In [172]:

    
size = 2 # have to be even
k=0.04

corner = np.zeros([maxx, maxy])
for x in range(maxx+1):
    for y in range(maxy+1):
        matrix = submatrix(x, y, size)
        A = np.zeros([2,2])
        A[0,0] = gradientx(matrix)
        A[0,1] = gradientxy(matrix)
        A[1,0] = gradientxy(matrix)
        A[1,1] = gradienty(matrix)
        R = np.linalg.det(A)-k*np.trace(A)
        if R>0:
            corner[x,y] = R
plt.imshow(corner, interpolation='nearest')









    Out[172]:





<matplotlib.image.AxesImage at 0x114560510>



In [173]:

    
plt.imshow(corner>1000000, cmap='Greys',  interpolation='nearest')









    Out[173]:





<matplotlib.image.AxesImage at 0x114665b50>

Demo 1 : UNIST Logo



In [312]:

    
im = Image.open('unist_logo.jpg') # Load the image file
imarray = np.array(im) # Construct the array
print 'Image Size : ', im.size
imgray = im.convert('L') # Convert the image to gray scale
plt.imshow(imgray, cmap=plt.get_cmap('gray'), interpolation='nearest')
"""interpolation='nearest' simply display the image without try to interpolate betwen pixels if the display resolution is not the same as the image resolution (which is most often the case).
It will results in an image in which is image pixel is displayed as a square of multiple display pixels."""

I = np.array(imgray) # Let imgray_array be the array of imgray.

maxx = I.shape[0]-1
maxy = I.shape[1]-1









    



Image Size :  (600, 600)



In [191]:

    
plt.hist(imgray_array.flatten(), 256, [0,256]);



In [332]:

    
size = 2 # have to be even
k=0.04

tic = timeit.default_timer() # for time check
start = tic
timeindex = 0
corner = np.zeros([maxx, maxy])
for x in range(maxx):
    for y in range(maxy):
        matrix = submatrix(x, y, size)
        A = np.zeros([2,2])
        A[0,0] = gradientx(matrix)
        A[0,1] = gradientxy(matrix)
        A[1,0] = gradientxy(matrix)
        A[1,1] = gradienty(matrix)
        R = np.linalg.det(A)-k*np.trace(A)
        if R>0:
            corner[x,y] = R
    if timeindex % 100 == 0:
        toc = timeit.default_timer()
        print 'row num : ', timeindex, 'sec:', toc-tic
        tic = toc
    timeindex += 1
print "FINISHED : ", timeit.default_timer()-start









    



row num :  0 sec: 0.164710998535
row num :  100 sec: 15.0754880905
row num :  200 sec: 15.4575419426
row num :  300 sec: 15.7959108353
row num :  400 sec: 15.7063159943
row num :  500 sec: 15.6441760063
FINISHED :  93.1917200089



In [260]:

    
plt.imshow(corner, interpolation='nearest')









    Out[260]:





<matplotlib.image.AxesImage at 0x13145f490>



In [333]:

    
fig, (ax1, ax2, ax3, ax4) = plt.subplots(1,4)

fig.suptitle('Harris Corner Detect(size=2), Demo1 : UNIST Logo', fontsize=20)

# Line plots
ax1.set_title('Original Gray Scale Image')
ax1.imshow(imgray, cmap=plt.get_cmap('gray'), interpolation='nearest')

ax2.set_title('Harris Corner Detect with R values')
ax2.imshow(corner, interpolation='nearest')

ax3.set_title('R > 10')
ax3.imshow(corner>10, cmap='Greys',  interpolation='nearest')

ax4.set_title('R > 200000')
ax4.imshow(corner>200000, cmap='Greys',  interpolation='nearest')

fig.set_size_inches(18,5)



In [313]:

    
size = 10 # have to be even
k=0.04

tic = timeit.default_timer() # for time check
start = tic
timeindex = 0
corner = np.zeros([maxx, maxy])
for x in range(maxx):
    for y in range(maxy):
        matrix = submatrix(x, y, size)
        A = np.zeros([2,2])
        A[0,0] = gradientx(matrix)
        A[0,1] = gradientxy(matrix)
        A[1,0] = gradientxy(matrix)
        A[1,1] = gradienty(matrix)
        R = np.linalg.det(A)-k*np.trace(A)
        if R>0:
            corner[x,y] = R
    if timeindex % 100 == 0:
        toc = timeit.default_timer()
        print 'row num : ', timeindex, 'sec:', toc-tic
        tic = toc
    timeindex += 1
print "FINISHED : ", timeit.default_timer()-start









    



row num :  0 sec: 0.837815999985
row num :  100 sec: 84.2350687981
row num :  200 sec: 83.5878801346
row num :  300 sec: 80.0482399464
row num :  400 sec: 83.1345701218
row num :  500 sec: 83.2046859264
FINISHED :  496.835048914



In [331]:

    
fig, (ax1, ax2, ax3, ax4) = plt.subplots(1,4)

fig.suptitle('Harris Corner Detect(size=10), Demo1 : UNIST Logo', fontsize=20)

# Line plots
ax1.set_title('Original Gray Scale Image')
ax1.imshow(imgray, cmap=plt.get_cmap('gray'), interpolation='nearest')

ax2.set_title('Harris Corner Detect with R values')
ax2.imshow(corner, interpolation='nearest')

ax3.set_title('R > 1000000000')
ax3.imshow(corner>1000000000, cmap='Greys',  interpolation='nearest')

ax4.set_title('R > 100000000000')
ax4.imshow(corner>100000000000, cmap='Greys',  interpolation='nearest')

fig.set_size_inches(18,5)



In [328]:

    
plt.imshow(corner>220000000000, cmap='Greys',  interpolation='nearest')









    Out[328]:





<matplotlib.image.AxesImage at 0x1247cc450>

Demo 2



In [246]:

    
im = Image.open('IU.jpg') # Load the image file
imarray = np.array(im) # Construct the array
print 'Image Size : ', im.size
imgray = im.convert('L') # Convert the image to gray scale
plt.imshow(imgray, cmap=plt.get_cmap('gray'), interpolation='nearest')

I = np.array(imgray) # Let imgray_array be the array of imgray.

maxx = I.shape[0]-1
maxy = I.shape[1]-1









    



Image Size :  (480, 300)



In [247]:

    
size = 2 # have to be even
k=0.04

tic = timeit.default_timer() # for time check
start = tic
timeindex = 0
corner = np.zeros([maxx, maxy])
for x in range(maxx):
    for y in range(maxy):
        matrix = submatrix(x, y, size)
        A = np.zeros([2,2])
        A[0,0] = gradientx(matrix)
        A[0,1] = gradientxy(matrix)
        A[1,0] = gradientxy(matrix)
        A[1,1] = gradienty(matrix)
        R = np.linalg.det(A)-k*np.trace(A)
        if R>0:
            corner[x,y] = R
    if timeindex % 100 == 0:
        toc = timeit.default_timer()
        print 'row num : ', timeindex, 'sec:', toc-tic
        tic = toc
    timeindex += 1
print "FINISHED : ", timeit.default_timer()-start









    



row num :  0 sec: 0.13964009285
row num :  100 sec: 12.2578580379
row num :  200 sec: 12.3582618237
FINISHED :  36.8223719597



In [257]:

    
fig, (ax1, ax2, ax3, ax4) = plt.subplots(1,4)

fig.suptitle('Harris Corner Detect, Demo2 : IU Face (Person Image)', fontsize=20)

# Line plots
ax1.set_title('Original Gray Scale Image')
ax1.imshow(imgray, cmap=plt.get_cmap('gray'), interpolation='nearest')

ax2.set_title('Harris Corner Detect with R values')
ax2.imshow(corner, interpolation='nearest')

ax3.set_title('R > 100000')
ax3.imshow(corner>100000, cmap='Greys',  interpolation='nearest')

ax4.set_title('R > 1000000')
ax4.imshow(corner>1000000, cmap='Greys',  interpolation='nearest')

fig.set_size_inches(18,4)

Demo 3



In [235]:

    
im = Image.open('IU_text.jpg') # Load the image file
imarray = np.array(im) # Construct the array
print 'Image Size : ', im.size
imgray = im.convert('L') # Convert the image to gray scale
plt.imshow(imgray, cmap=plt.get_cmap('gray'), interpolation='nearest')

I = np.array(imgray) # Let imgray_array be the array of imgray.

maxx = I.shape[0]-1
maxy = I.shape[1]-1









    



Image Size :  (700, 467)



In [236]:

    
size = 2 # have to be even
k=0.04

tic = timeit.default_timer() # for time check
start = tic
timeindex = 0
corner = np.zeros([maxx, maxy])
for x in range(maxx):
    for y in range(maxy):
        matrix = submatrix(x, y, size)
        A = np.zeros([2,2])
        A[0,0] = gradientx(matrix)
        A[0,1] = gradientxy(matrix)
        A[1,0] = gradientxy(matrix)
        A[1,1] = gradienty(matrix)
        R = np.linalg.det(A)-k*np.trace(A)
        if R>0:
            corner[x,y] = R
    if timeindex % 100 == 0:
        toc = timeit.default_timer()
        print 'row num : ', timeindex, 'sec:', toc-tic
        tic = toc
    timeindex += 1
print "FINISHED : ", timeit.default_timer()-start









    



row num :  0 sec: 0.197516918182
row num :  100 sec: 18.1818010807
row num :  200 sec: 18.1433038712
row num :  300 sec: 18.2068171501
row num :  400 sec: 18.1124298573
FINISHED :  84.6004538536



In [244]:

    
fig, (ax1, ax2, ax3, ax4) = plt.subplots(1,4)

fig.suptitle('Harris Corner Detect, Demo3 : IU with text', fontsize=20)

# Line plots
ax1.set_title('Original Gray Scale Image')
ax1.imshow(imgray, cmap=plt.get_cmap('gray'), interpolation='nearest')

ax2.set_title('Harris Corner Detect with R values')
ax2.imshow(corner, interpolation='nearest')

ax3.set_title('R > 10')
ax3.imshow(corner>10, cmap='Greys',  interpolation='nearest')

ax4.set_title('R > 5000000')
ax4.imshow(corner>5000000, cmap='Greys',  interpolation='nearest')

fig.set_size_inches(18,4)

Demo 4



In [270]:

    
im = Image.open('map.jpg') # Load the image file
imarray = np.array(im) # Construct the array
print 'Image Size : ', im.size
imgray = im.convert('L') # Convert the image to gray scale
plt.imshow(imgray, cmap=plt.get_cmap('gray'), interpolation='nearest')

I = np.array(imgray) # Let imgray_array be the array of imgray.

maxx = I.shape[0]-1
maxy = I.shape[1]-1









    



Image Size :  (1309, 585)



In [271]:

    
size = 2 # have to be even
k=0.04

tic = timeit.default_timer() # for time check
start = tic
timeindex = 0
corner = np.zeros([maxx, maxy])
for x in range(maxx):
    for y in range(maxy):
        matrix = submatrix(x, y, size)
        A = np.zeros([2,2])
        A[0,0] = gradientx(matrix)
        A[0,1] = gradientxy(matrix)
        A[1,0] = gradientxy(matrix)
        A[1,1] = gradienty(matrix)
        R = np.linalg.det(A)-k*np.trace(A)
        if R>0:
            corner[x,y] = R
    if timeindex % 100 == 0:
        toc = timeit.default_timer()
        print 'row num : ', timeindex, 'sec:', toc-tic
        tic = toc
    timeindex += 1
print "FINISHED : ", timeit.default_timer()-start









    



row num :  0 sec: 0.344913959503
row num :  100 sec: 34.0332839489
row num :  200 sec: 34.2081160545
row num :  300 sec: 33.7014529705
row num :  400 sec: 32.9766590595
row num :  500 sec: 32.9560739994
FINISHED :  196.576879025



In [291]:

    
fig, (ax1, ax2, ax3, ax4) = plt.subplots(1,4)

fig.suptitle('Harris Corner Detect, Demo4 : Google Maps', fontsize=20)

# Line plots
ax1.set_title('Original Gray Scale Image')
ax1.imshow(imgray, cmap=plt.get_cmap('gray'), interpolation='nearest')

ax2.set_title('Harris Corner Detect with R values')
ax2.imshow(corner, interpolation='nearest')

ax3.set_title('R > 10')
ax3.imshow(corner>10, cmap='Greys',  interpolation='nearest')

ax4.set_title('R > 5000000')
ax4.imshow(corner>5000000, cmap='Greys',  interpolation='nearest')

fig.set_size_inches(18,3)



In [292]:

    
fig, (ax1, ax2, ax3, ax4) = plt.subplots(1,4)

fig.suptitle('Harris Corner Detect, Demo4 : Google Maps', fontsize=20)

# Line plots
ax1.set_title('R > 10000000')
ax1.imshow(corner>10000000, cmap='Greys',  interpolation='nearest')

ax2.set_title('R > 40000000')
ax2.imshow(corner>40000000, cmap='Greys',  interpolation='nearest')

ax3.set_title('R > 70000000')
ax3.imshow(corner>70000000, cmap='Greys',  interpolation='nearest')

ax4.set_title('R > 100000000')
ax4.imshow(corner>100000000, cmap='Greys',  interpolation='nearest')

fig.set_size_inches(18,3)

Modified with size 10 before (window) size=2



In [293]:

    
size = 10 # have to be even
k=0.04

tic = timeit.default_timer() # for time check
start = tic
timeindex = 0
corner = np.zeros([maxx, maxy])
for x in range(maxx):
    for y in range(maxy):
        matrix = submatrix(x, y, size)
        A = np.zeros([2,2])
        A[0,0] = gradientx(matrix)
        A[0,1] = gradientxy(matrix)
        A[1,0] = gradientxy(matrix)
        A[1,1] = gradienty(matrix)
        R = np.linalg.det(A)-k*np.trace(A)
        if R>0:
            corner[x,y] = R
    if timeindex % 100 == 0:
        toc = timeit.default_timer()
        print 'row num : ', timeindex, 'sec:', toc-tic
        tic = toc
    timeindex += 1
print "FINISHED : ", timeit.default_timer()-start









    



row num :  0 sec: 1.90029597282
row num :  100 sec: 185.654971123
row num :  200 sec: 180.005522013
row num :  300 sec: 180.201139927
row num :  400 sec: 185.08220005
row num :  500 sec: 183.339393854
FINISHED :  1069.03331399



In [307]:

    
fig, (ax1, ax2, ax3, ax4) = plt.subplots(1,4)

fig.suptitle('Harris Corner Detect, Demo4 : Google Maps', fontsize=20)

# Line plots
ax1.set_title('Original Gray Scale Image')
ax1.imshow(imgray, cmap=plt.get_cmap('gray'), interpolation='nearest')

ax2.set_title('Harris Corner Detect with R values')
ax2.imshow(corner, interpolation='nearest')

ax3.set_title('R > 1000000')
ax3.imshow(corner>1000000, cmap='Greys',  interpolation='nearest')

ax4.set_title('R > 1000000000')
ax4.imshow(corner>1000000000, cmap='Greys',  interpolation='nearest')

fig.set_size_inches(18,3)



In [308]:

    
plt.imshow(corner>10000000000, cmap='Greys',  interpolation='nearest')









    Out[308]:





<matplotlib.image.AxesImage at 0x119749450>



In [309]:

    
plt.imshow(corner>100000000000, cmap='Greys',  interpolation='nearest')









    Out[309]:





<matplotlib.image.AxesImage at 0x12ded1c10>

Conclusion

For the image size (1309, 585), when we take 2 as a size of window, the computation time is 196 sec. And when we take 10 a size of window, the time is 1069 sec. ($2:196=10:1069$, $196*5=980$)

Harris Corner Detect[1] is a simple idea. It means interestingly it has many of holes. So we can add various ideas [2]. []

In this algorithm, the parameter k, size(window) can't be calculated. It's only depended experiences. It's big defect.

[1]: Harris, Chris, and Mike Stephens. "A combined corner and edge detector." Alvey vision conference. Vol. 15. 1988.



In [ ]: