from PIL import Image
import numpy as np
import matplotlib.pyplot as plt
from scipy.ndimage import filters

## 1.找出角點(diǎn)
### 1.1 步驟：img求導(dǎo) --> imx ,imy 濾波-->Wxx ,Wxy,Wyy --> （Wxx*Wyy - Wxy**2）/（Wxx + Wyy）
def compute_harris_response(img,sigma=3):
? ? # 求梯度
? ? imgx,imgy ?= np.zeros(img.shape),np.zeros(img.shape)
? ? filters.gaussian_filter(img,(sigma,sigma),(0,1),imgx)
? ? filters.gaussian_filter(img,(sigma,sigma),(1,0),imgy) # [260,263]
? ? # 對梯度進(jìn)行高斯濾波
? ? wxx = filters.gaussian_filter(imgx**2,sigma)
? ? wyy = filters.gaussian_filter(imgy**2,sigma) # [260,263]
? ? wxy = filters.gaussian_filter(imgx*imgy,sigma)# [260,263]
? ? ## 求行列式和跡
? ? wdet = wxx*wyy -wxy**2
? ? wtr = wxx + wyy
? ? return ?wdet/wtr

## 2 篩選角點(diǎn)
### 2.1 步驟：根據(jù)閾值篩選角點(diǎn)--> 取得角點(diǎn)的坐標(biāo) -->根據(jù)角點(diǎn)坐標(biāo)得到角點(diǎn)所在的行 ?-->
# --> 在角點(diǎn)周圍，刪除掉其他角點(diǎn)
def get_harris_points(harri,min_dist=4,threshold=0.1):
? ? corner_thre = harri.max()*threshold ?# 角點(diǎn)閾值
? ? mask = (harri > corner_thre)*1 ?# 取出大于閾值的點(diǎn)為候選角點(diǎn)
? ? cords = np.array(mask.nonzero()).T ?# 取候選角點(diǎn)的坐標(biāo)

? ? values = [harri[i[0],i[1]] for i in cords] ?# 候選角點(diǎn)的值
? ? cls = np.argsort(values) ? ? # 對角點(diǎn)排序得到排序后的序列號，序列號也是候選角點(diǎn)所在的行

? ? loc = np.zeros(harri.shape) ? # 劃出可行性區(qū)域
? ? loc[min_dist:-min_dist,min_dist:-min_dist] = 1

? ? re_cords = []
? ? for i in cls: ?# 篩選角點(diǎn)。先取出角點(diǎn)，角點(diǎn)周圍的點(diǎn)不再取出
? ? ? ? if loc[cords[i,0],cords[i,1]] == 1 :
? ? ? ? ? ? re_cords.append(cords[i])
? ? ? ? ? ? loc[cords[i,0]-min_dist:cords[i,0]+min_dist,cords[i,1]-min_dist:cords[i,1]+min_dist]=0
? ? return re_cords

def plot_harri(img,cords):
? ? plt.figure()
? ? plt.gray()
? ? plt.imshow(img)
? ? plt.plot([i[1] for i in cords],[i[0] for i in cords],'.')
? ? plt.axis('off')
? ? plt.show()

## 3 測試
if __name__ == '__main__':
? ? img = np.array(Image.open('luna.png').convert('L'))
? ? harri = compute_harris_response(img)
? ? re_cords = get_harris_points(harri)
? ? plot_harri(img,re_cords)