code.py

import cv2
import numpy
import numpy as np
import math
import copy

target = 'target.jpg'
test   = '3.jpg'


def atan2X(y,x):
    tmp = math.atan2(y,x)
    if tmp<0:
        return tmp + math.pi*2 
    else: return tmp

def getPart(y,x):
    return int(atan2X(y,x) / (math.pi * 2 /36))

def getPart(theta):
    return int(theta / (math.pi * 2 / 36))

def getNearByTheta(x,y):
    print (img[x][y+1], img[x][y-1], img[x+1][y], img[x-1][y],"delta=", int(img[x][y+1]) - int(img[x][y-1]), int(img[x+1][y]) - int(img[x-1][y]))
    return atan2X(int(img[x][y+1]) - int(img[x][y-1]), 
        int(img[x+1][y])- int(img[x-1][y]))

def getNearByLength(x,y):
    return math.hypot(int(img[x][y+1])-int(img[x][y-1]), int(img[x+1][y]) - int(img[x-1][y]))

def transf(x,y, theta):
    return int( math.cos(theta) * float(x) - math.sin(theta) * float(y)), int(math.sin(theta) * float(x) + math.cos(theta) * float(y))
    # tmptheta = math.atan2(y,x)
    # if tmptheta<0 : tmptheta += math.pi * 2
    # tmptheta += theta
    # if tmptheta > math.pi * 2 :
        # tmptheta -= math.pi * 2

    # tup = ((1,0), (1,1), (0,1), (-1,1), (-1,0), (-1,-1), (0,-1), (1,-1), (1,0))
    # for i in range(9):
        # if tmptheta < (math.pi / 4 * i + math.pi / 8:
            # return tup[i]

def getDescriptor():
    res = cv2.goodFeaturesToTrack(img, 500, 0.005, 10, blockSize=3, k=0.04, useHarrisDetector= True)
    ans = []
    tp = copy.deepcopy(img)
    for p in res:
        po = int(p[0][0]), int(p[0][1])
        if po[0] < 30 or po[0] > len(img)-30 or po[1] < 30 or po[1] > len(img[0])-30:
            continue

        a = {}
        for dx in range(-8, 8):
            for dy in range(-8, 8):
                theta = getNearByTheta(po[0]+dx, po[1]+dy)
                if theta == -999:
                    continue
                m = getNearByLength(po[0]+dx, po[1]+dy)
                a[getPart(theta)] = a.get(getPart(theta), 0.0) + m
                print (theta, getPart(theta))
        
        ma = 0.0
        mai = 0
        for i in a:
            if (a[i]>ma):
                ma = a[i]
                mai = i
        offTheta = mai * (math.pi / 18.0) 
        print ("result", po[0], po[1], offTheta, mai, ma)
        cv2.line(tp, (po[0], po[1]), (po[0] + int(math.cos(offTheta) * 50), po[1] + int(math.sin(offTheta) * 50)), (100,255,100))
        cv2.circle(tp, (po[0], po[1]), 10, (100,155,100))

        descriptor = []
        for i in range(16):
            descriptor.append({})

        cnt = 0
        cur = 0
        for i in range(-8, 8):
            for j in range(-8,8):
                cur = (i+8) / 4 * 4 + (j+8)/4
                vec1b = transf(i+1, j, offTheta)
                vec1a = transf(i-1, j, offTheta)
                vec2b = transf(i,j+1, offTheta)
                vec2a = transf(i, j-1, offTheta)
                
                t = math.atan2(int(img[po[0] + vec2b[0]][po[1] + vec2b[1]]) - int(img[po[0]+vec2a[0]][po[1] + vec2a[1]]),
                        int(img[po[0]+vec1b[0]][po[1]+vec1b[1]]) - int(img[po[0] + vec1a[0]][po[1] + vec1a[1]]))

                # print int(img[po[0] + vec2[0]][po[1] + vec2[1]]) - int(img[po[0]-vec2[0]][po[1] - vec2[1]]), int(img[po[0]+vec1[0]][po[1]+vec1[1]]) - int(img[po[0] - vec1[0]][po[1] - vec1[1]])
                
                # print math.hypot(int(img[po[0] + vec2[0]][po[1] + vec2[1]]) - int(img[po[0]-vec2[0]][po[1] - vec2[1]]), int(img[po[0]+vec1[0]][po[1]+vec1[1]]) - int(img[po[0] - vec1[0]][po[1] - vec1[1]]))

                while t<0: t += math.pi * 2
                while t>math.pi*2: t-= math.pi * 2

                result = int(t / (math.pi / 4))
                descriptor[cur][result] = descriptor[cur].get(result, 0) + 1
                # int(math.hypot(int(img[po[0] + vec2b[0]][po[1] + vec2b[1]]) - int(img[po[0]+vec2a[0]][po[1] + vec2a[1]]), int(img[po[0]+vec1b[0]][po[1]+vec1b[1]]) - int(img[po[0] + vec1a[0]][po[1] + vec1a[1]])))
        print ("des", descriptor)

        descriptorVec = []
        sum = 0
        for i in range(16):
            for j in range(8):
                sum += descriptor[i].get(j, 0)**2
        sum = math.sqrt(sum)

        if sum > 10:
            for i in range(16):
                for j in range(8):
                    descriptorVec.append(descriptor[i].get(j,0) * 1.0 / sum)

            ans.append((po[0], po[1], descriptorVec))
    return ans,tp

img =  cv2.imread(target, 0)
img1 = copy.deepcopy(img)
base,cc1 = getDescriptor()
cv2.imwrite('circle_1.jpg', cc1)

img = cv2.imread(test, 0)
img2 = copy.deepcopy(img)
vec,cc2 = getDescriptor()
cv2.imwrite('circle_2.jpg', cc2)


rows1 = img1.shape[0]
cols1 = img1.shape[1]
rows2 = img2.shape[0]
cols2 = img2.shape[1]

out = np.zeros((max([rows1,rows2]),cols1+cols2,3), dtype='uint8')

# Place the first image to the left
out[:rows1,:cols1] = np.dstack([img1, img1, img1])

# Place the next image to the right of it
out[:rows2,cols1:] = np.dstack([img2, img2, img2])

print base, vec

anscnt = 0
while True:
    mi = 0
    mi1 = 0
    mi2 = 0
    for c1 in range(len(base)):
        v1 = base[c1][2]
        for c2 in range(len(vec)):
            v2 = vec[c2][2]
            s=0
            for i in range(len(v2)):
                s+=v1[i] * v2[i]
            if s>mi:
                mi = s
                mi1 = c1
                mi2 = c2
    print mi
    if anscnt<5 :
        cv2.circle(out, (base[mi1][0], base[mi1][1]), 10, (100, 100, 255))
        cv2.circle(out, (vec[mi2][0]+cols1, vec[mi2][1]), 10, (100, 100, 100))
        cv2.line(out, (int(base[mi1][0]), int(base[mi1][1])), (int(vec[mi2][0])+cols1, int(vec[mi2][1])), (255,0,0), 1)
        del base[mi1]
        del vec[mi2]
        anscnt += 1
    else:
        break





cv2.imwrite('out.jpg', out)