Python图像处理（二）局部图像描述字
（仅学习摘抄）
1. Harris 角点检测器
Harris 角点检测算法是一个极为简单的角点检测算法。主要思想：如果像素周围显示存在多于一个方向的边，我们认为该点为兴趣点，这个点就称为角点。
角点，在通常意义来说，就是极值点，在某方面属性特别突出的点，是在某些属性上强度最大或者最小的孤立点、线段的终点。
① 一阶导数（即灰度的梯度）的局部最大所对应的像素点；
② 两条及两条以上边缘的交点；
③ 图像中梯度值和梯度方向的变化速率都很高的点；
④ 角点处的一阶导数最大、二阶导数为零，指示物体边缘变化不连续的地方。
先计算 Harris 矩阵，然后计算出角点：
def compute_harris_response(im,sigma=5):
'''在一幅灰度图像中，对每个像素计算Harris角点检测器响应函数'''
# 计算导数
imx = zeros(im.shape)
filters.gaussian_filter(im,(sigma,sigma),(0,1),imx)
imy = zeros(im.shape)
filters.gaussian_filter(im,(sigma,sigma),(1,0),imy)
# 计算Harris矩阵的分量
Wxx = filters.gaussian_filter(imx*imx,sigma)
Wxy = filters.gaussian_filter(imx*imy,sigma)
Wyy = filters.gaussian_filter(imy*imy,sigma)
# 计算特征值和迹
Wdet = Wxx*Wyy - Wxy**2
Wtr = Wxx + Wyy
return Wdet/Wtr
def get_harris_points(harrisim,min_dist=10,threshold=0.1):
'''从一幅Harris响应图像中返回角点。min_dist为分割角点和图像边界的最少像素数目'''
# 寻找高于阈值的候选角点
corner_threshold = harrisim.max() * threshold
harrisim_t = (harrisim > corner_threshold) * 1
# 得到候选点的坐标
coords = array(harrisim_t.nonzero()).T
# 以及它们的Harris响应值
candidate_values = [harrisim[c[0],c[1]] for c in coords]
# 对候选点按照Harris响应值进行排序
index = argsort(candidate_values)
# 将可行点的位置保存到数组中
allowed_locations = zeros(harrisim.shape)
allowed_locations[min_dist:-min_dist,min_dist:-min_dist] = 1
# 按照min_distance原则，选择最佳Harris点
filtered_coords = []
for i in index:
if allowed_locations[coords[i,0],coords[i,1]] == 1:
filtered_coords.append(coords[i])
allowed_locations[(coords[i,0]-min_dist):(coords[i,0]+min_dist),(coords[i,1]-min_dist):(coords[i,1]+min_dist)] = 0
return filtered_coords
def plot_harris_points(image,filtered_coords):
'''绘制图像中检测到的角点'''
figure()
gray()
imshow(image)
plot([p[1] for p in filtered_coords],[p[0] for p in filtered_coords],'*')
axis('off')
show()
设置不同的阈值，将会得到不同数目的角点。
对比两幅图像的兴趣点，找到匹配角点。需要在每个点上加入描述子信息。
兴趣点描述子是分配给兴趣点的一个向量，描述该点附近的图像的表观信息。描述子越好，寻找到的对应点越好。
def get_descriptors(image,filtered_coords,wid=5):
'''对于每个返回点，返回点周围 2*wid+1 个像素的值（假设选取点的min_distance>wid）'''
desc = []
for coords in filtered_coords:
path = image[coords[0]-wid:coords[0]+wid+1,coords[1]-wid:coords[1]+wid+1].flatten()
desc.append(path)
return desc
def match(desc1,desc2,threshold=0.5):
'''对于第一幅图像中的每个角点描述子，使用归一化互相关，选取它在第二幅图像中的匹配角点'''
n = len(desc1[0])
# 点对的距离
d = -ones((len(desc1),len(desc2)))
for i in range(len(desc1)):
for j in range(len(desc2)):
d1 = (desc1[i] - mean(desc1[i])) / std(desc1[i])
d2 = (desc2[j] - mean(desc2[j])) / std(desc2[j])
ncc_value = sum(d1 * d2) / (n - 1)
if ncc_value > threshold:
d[i,j] = ncc_value
ndx = argsort(-d)
matchscores = ndx[:,0]
return matchscores
def match_twosided(desc1,desc2,threshold=0.5):
'''两边对称版本的match()'''
matches_12 = match(desc1,desc2,threshold)
matches_21 = match(desc2,desc1,threshold)
ndx_12 = where(matches_12 >= 0)[0]
# 去除非对称的匹配
for n in ndx_12:
if matches_21[matches_12[n]] != n:
matches_12[n] = -1
return matches_12
def appendimage(im1,im2):
'''返回将两幅图像并排拼接成的一幅新图像'''
# 选取具有最少行数的图像，然后填充足够的空行
rows1 = im1.shape[0]
rows2 = im2.shape[0]
if rows1 < rows2:
im1 = concatenate((im1,zeros((rows2-rows1,im1.shape[1]))),axis=0)
elif rows1 > rows2:
im2 = concatenate((im2,zeros((rows2-rows1,im2.shape[1]))),axis=0)
# 如果这些情况都没有，那么它们的行数相同，不需要进行填充
return concatenate((im1,im2), axis=1)
def plot_matches(im1,im2,locs1,locs2,matchscores,show_below=True):
'''显示一幅带有连接匹配之间连线的图片
输入：im1,im2(数组图像)，locs1,locs2(特征位置)，matchscores(match()的输出)，
show_below(如果图像应该显示在匹配的下方)'''
im3 = appendimage(im1,im2)
if show_below:
im3 = vstack((im3,im3))
imshow(im3)
cols1 = im1.shape[1]
for i,m in enumerate(matchscores):
if m>0:
plot([locs1[i][1],locs2[m][1]+cols1],[locs1[i][0],locs2[m][0]],'c')
axis('off')
存在不正确的匹配，不具有尺度不变性和旋转不变性，而且算法中像素块的大小也会影响对应匹配的结果。
2. SIFT（尺度不变特性变换）
SITF 特征包括兴趣点检测和描述子。SITF 具有非常强的稳健性。SITF 特征对于尺度、旋转和亮度都具有不变性，它可以用于三维视角和噪声的可靠匹配。
2.1 兴趣点
SIFT 特征使用高斯差分函数来定位兴趣点：
其中，Gσ 时二维高斯核，Iσ 时使用 Gσ 模糊的灰度图像，k 是决定相差尺度的常数。兴趣点是在图像位置和尺度变化下 D(x,σ) 的最大值和最小值点。这些候选位置点通过滤波去除不稳定点。基于一些准则，比如认为低于对比度和位于边上的点不是兴趣点，可以去除一些候选兴趣点。
2.2 描述子
兴趣点（关键点）位置描述子给出了兴趣点的位置和尺度信息。为了实现旋转不变性，基于每个点周围图像梯度的方向和大小，SIFT 描述子又引入了参考方向。SIFT 描述子使用主方向描述参考方向。主方向使用方向直方图（以大小为权重）来度量。
为了对图像亮度具有稳健性，SIFT描述子使用图像梯度(之前Harris描述子使用图像亮度信息计算归-一化互相关矩阵)。SIFT 描述子在每个像素点附近选取子区域网格，在每个子区域内计算图像梯度方向直方图。每个子区域的直方图拼接起来组成描述子向量。SIFT 描述子的标准设置使用4x4的子区域，每个子区域使用8个小区间的方向直方图，会产生共128个小区间的直方图(4x4x8=128)。
2.3 检测兴趣点
from PIL import Image
from numpy import *
from pylab import *
import os
def process_image(imagename,resultname,params="--edge-thresh 10 --peak-thresh 5"):
'''处理一幅图像，然后将结果保存在文件中'''
if imagename[-3:] != 'pgm':
# 创建一个pgm文件
im = Image.open(imagename).convert('L')
im.save('tmp.pgm')
imagename = 'tmp.pgm'
cmmd = str("sift "+imagename+" --output="+resultname+" "+params)
os.system(cmmd)
print('processd',imagename,'to',resultname)
def read_features_from_file(filename):
'''读取特征属性值，然后将其以矩阵的形式返回'''
f = loadtxt(filename)
return f[:,:4],f[:,4:] # 特征位置，描述子
def write_features_to_file(filename,locs,desc):
'''将特征位置和描述子保存到文件中'''
savetxt(filename,hstack((locs,desc)))
def plot_features(im,locs,circle=False):
'''显示带有特征的图像
输入：im（数组图像），locs（每个特征的行、列、尺寸和朝向）'''
def draw_circle(c,r):
t = arange(0,1.01,.01)*2*pi
x = r*cos(t) + c[0]
y = r*sin(t) + c[1]
plot(x,y,'b',linewidth=2)
imshow()
if circle:
for p in locs:
draw_circle(p[:2],p[2])
else:
plot(locs[:,0],locs[:,1],'ob')
axis('off')
在原始图像中使用蓝色的点标记出 SIFT 特征点的位置。参数 circle 设为 True，将会绘制出圆圈，圆圈的半径是特征的尺度。
2.4 匹配描述子
将一幅图像中的特征匹配到另一副图像的特征，一种稳健的准则是使用这两个特征距离和两个最匹配特征距离的比率。和 Harris 匹配类似，匹配两张图片具有相同特征的点。
