计算机视觉二

一.角点检测器

什么是角点？在任意方向的一个微小变动都会引起灰度很大变化的点

1. 角点是一阶导数（即灰度的梯度）的局部最大所对应的像素点；
2. 角点是两条及两条以上边缘的交点;
3. 角点指示了物体边缘变化不连续的方向;
4. 角点处的一阶导数最大，二阶导数为0；
5. 角点是指图像中梯度值和梯度方向的变化速率都很高的点。

目前的角点检测算法可以归结为以下三类

基于灰度图的角点检测
基于二值图的角点检测
基于轮廓曲线的角点检测
基于灰度图的角点检测又可以分为基于梯度、基于模板和基于模板梯度组合的三类方法。

其中基于模板的方法主要考虑像素领域点的灰度变化，即图像亮度的变化，将与邻点亮度对比足够大的点定义为角点，Harris角点检测算法就是基于模板的角点检测方法。

 代码

# -*- coding:utf-8 -*-
from PIL import Image
from numpy import *
from pylab import *
from scipy.ndimage import filters

def compute_harris_response(im, sigma=3):
    """ 在一幅灰度图像中，对每一个像素计算Harris角点检测器响应函数
     im:(数组图像）  sigma=3:标准差为3
    """

    #计算x方向的导数
    imx = zeros(im.shape) #zeros()生成全零数组
    filters.gaussian_filter(im, (sigma,sigma),(0,1),imx)

    #计算x方向的导数
    imy = zeros(im.shape) #zeros()生成全零数组
    filters.gaussian_filter(im, (sigma,sigma),(0,1),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(harrism, min_dist=10, threshold=0.1):
    """从一幅Harris响应图像中返回角点。min_dist为分割角点和图像边界的最少像素数目"""

    #寻找高于阈值的候选角点
    corner_threshold = harrism.max()*threshold
    harrisim_t = (harrisim > corner_threshold) * 1

    #得到候选点的坐标
    coords = array(harrisim_t.nonzero()).T

    #以及它们的Harris响应值
    candidate_values = [harrism[c[0],c[1]] for c in coords]

    #对侯选点按照Harris响应值进行排序
    index = argsort(candidate_values)

    #将可行点的位置保存到数组中
    allowed_locations = zeros(harrism.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()

if __name__ == "__main__":

    im = array(Image.open("/Users/apple/Desktop/集大.jpeg").convert("L"))
    harrisim = compute_harris_response(im)
    filtered_coords = get_harris_points(harrisim,6)
    plot_harris_points(im, filtered_coords)

结果

原图

 二.sift

SIFT全称尺度不变特征变换，SIFT算子是把图像中检测到的特征点用一个128维的特征向量进行描述，因此一幅图像经过SIFT算法后表示为一个128维的特征向量集，该特征向量集具有对图像缩放，平移，旋转不变的特征，对于光照、仿射和投影变换也有一定的不变性，是一种非常优秀的局部特征描述算法。
SIFT算法步骤分为：
（1）尺度空间极点检测
（2）关键点精确定位
（3）确定关键点方向
（4）生成特征向量

尺度空间
特征点的检测就是寻找特征点的位置和尺度，需要位置的原因显而易见，而需要尺度的原因则是因为真实世界中的物体只有在一定尺度下才有意义。我们寻找的特征点就是要找到在连续的尺度空间下位置不发生改变的点。

构建尺度空间的目的就是找到在尺度变化中具有不变性的位置，可以使用连续的尺度变化，即在尺度空间中所有可能的尺度变化中找到稳定的特征点，通过这种方式找到的极点可以保证在图像缩放和旋转变化中具有不变性。

经过前人证明，尺度空间内核是高斯函数。因此假设I(x,y)是原始图像，G(x,y,σ)是尺度空间可变的高斯函数，则一个图像的尺度空间可以定义为：
其中，∗表示的是卷积运算，σ表示尺度空间的大小，σ越大则表示越模糊，表示图像的概貌，σ越小则表示越清晰，表示图像的细节。高斯函数G(x,y,σ)定义为：

经过一系列的尺度空间变换和二倍的下采样，最终得到高斯金字塔。

高斯差分（DOG）
高斯差分可以在尺度空间中找到稳定不变的极值点，高斯差分(DOG)函数D(x,y,σ)定义为：

其中kσ和σ是连续的两个图像的平滑尺度，所得到的差分图像再高斯差分金字塔中。

代码

import numpy as np
import cv2
from matplotlib import pyplot as plt

imgname1 = '/Users/apple/Desktop/集大.jpeg'
imgname2 = '/Users/apple/Desktop/集大2.jpeg'

sift = cv2.xfeatures2d.SIFT_create()

# FLANN 参数设计
FLANN_INDEX_KDTREE = 0
index_params = dict(algorithm = FLANN_INDEX_KDTREE, trees = 5)
search_params = dict(checks=50)
flann = cv2.FlannBasedMatcher(index_params,search_params)

img1 = cv2.imread(imgname1)
gray1 = cv2.cvtColor(img1, cv2.COLOR_BGR2GRAY) #灰度处理图像
kp1, des1 = sift.detectAndCompute(img1,None)#des是描述子
rows, cols = img1.shape[:2] #获取img1的高度、宽度

img2 = cv2.imread(imgname2)
img2 = cv2.resize(img2,(cols,rows),interpolation=cv2.INTER_CUBIC) #放大图像
gray2 = cv2.cvtColor(img2, cv2.COLOR_BGR2GRAY)
kp2, des2 = sift.detectAndCompute(img2,None)

hmerge = np.hstack((gray1, gray2)) #水平拼接
cv2.imshow("gray", hmerge) #拼接显示为gray
cv2.waitKey(0)

img3 = cv2.drawKeypoints(img1,kp1,img1,color=(255,0,255))
img4 = cv2.drawKeypoints(img2,kp2,img2,color=(255,0,255))

hmerge = np.hstack((img3, img4)) #水平拼接
cv2.imshow("point", hmerge) #拼接显示为gray
cv2.waitKey(0)
matches = flann.knnMatch(des1,des2,k=2)
matchesMask = [[0,0] for i in range(len(matches))]

good = []
for m,n in matches:
    if m.distance < 0.7*n.distance:
        good.append([m])

img5 = cv2.drawMatchesKnn(img1,kp1,img2,kp2,good,None,flags=2)

cv2.imshow("FLANN", img5)
cv2.waitKey(0)
cv2.destroyAllWindows()

结果

 三.匹配地理标记图像

from PIL import Image
import os
from numpy import *
from pylab import *
from scipy.ndimage import filters
def imresize(im,sz):
    pil_im = Image.fromarray(uint8(im))
    
    return array(pil_im.resize(sz))

def compute_harris_response(im,sigma=3):
    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)
    
    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):
    corner_threshold = harrisim.max() * threshold
    harrisim_t = (harrisim > corner_threshold) * 1
    
    coords = array(harrisim_t.nonzero()).T
    
    candidate_values = [harrisim[c[0],c[1]] for c in coords]
    
    index = argsort(candidate_values)[::-1]
    
    allowed_locations = zeros(harrisim.shape)
    allowed_locations[min_dist:-min_dist,min_dist:-min_dist] = 1
    
    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):
    """ Plots corners found in image. """
    
    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):
    """ For each point return pixel values around the point
        using a neighbourhood of width 2*wid+1. (Assume points are 
        extracted with min_distance > wid). """
    
    desc = []
    for coords in filtered_coords:
        patch = image[coords[0]-wid:coords[0]+wid+1,
                            coords[1]-wid:coords[1]+wid+1].flatten()
        desc.append(patch)
    
    return desc
def match(desc1,desc2,threshold=0.5):
    """ For each corner point descriptor in the first image, 
        select its match to second image using
        normalized cross correlation. """
    
    n = len(desc1[0])
    
    # pair-wise distances
    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):
    """ Two-sided symmetric version of match(). """
    
    matches_12 = match(desc1,desc2,threshold)
    matches_21 = match(desc2,desc1,threshold)
    
    ndx_12 = where(matches_12 >= 0)[0]
    
    # remove matches that are not symmetric
    for n in ndx_12:
        if matches_21[matches_12[n]] != n:
            matches_12[n] = -1
    
    return matches_12
def appendimages(im1,im2):
    """ Return a new image that appends the two images side-by-side. """
    
    # select the image with the fewest rows and fill in enough empty rows
    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((rows1-rows2,im2.shape[1]))),axis=0)
    # if none of these cases they are equal, no filling needed.
    
    return concatenate((im1,im2), axis=1)
def plot_matches(im1,im2,locs1,locs2,matchscores,show_below=True):
    """ Show a figure with lines joining the accepted matches 
        input: im1,im2 (images as arrays), locs1,locs2 (feature locations), 
        matchscores (as output from 'match()'), 
        show_below (if images should be shown below matches). """
    
    im3 = appendimages(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')
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("D:\Python36\Lib\win64\sift.exe "+imagename+" --output="+resultname+" "+params)
    os.system(cmmd)

def read_features_from_file(filename):
    """读取特征属性值，然后将其以矩阵的形式返回"""
    f = loadtxt(filename)
    return f[:,:4], f[:,4:] #特征位置，描述子

def write_featrues_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(im)
    if circle:
        for p in locs:
            draw_circle(p[:2], p[2])
    else: 
        plot(locs[:,0], locs[:,1], 'ob')
    axis('off')

def get_imlist(path):
        
    return [os.path.join(path,f) for f in os.listdir(path) if f.endswith('.jpeg')]

from pylab import *
from PIL import Image
import pydot
import os


download_path = "/Users/apple/Desktop/集大"  
path = "/Users/apple/Desktop/集大"  

imlist = get_imlist(download_path)
nbr_images = len(imlist)

featlist = [imname[:-3] + 'sift' for imname in imlist]
for i, imname in enumerate(imlist):
    process_image(imname, featlist[i])

matchscores = zeros((nbr_images, nbr_images))

for i in range(nbr_images):
    for j in range(i, nbr_images):  # only compute upper triangle
        print ('comparing ', imlist[i], imlist[j])
        l1, d1 = read_features_from_file(featlist[i])
        l2, d2 = read_features_from_file(featlist[j])
        matches = match_twosided(d1, d2)
        nbr_matches = sum(matches > 0)
        print ('number of matches = ', nbr_matches)
        matchscores[i, j] = nbr_matches
print ("The match scores is: \n", matchscores)

for i in range(nbr_images):
    for j in range(i + 1, nbr_images):  # no need to copy diagonal
        matchscores[j, i] = matchscores[i, j]
threshold = 2  # min number of matches needed to create link

g = pydot.Dot(graph_type='graph')  # don't want the default directed graph

for i in range(nbr_images):
    for j in range(i + 1, nbr_images):
        if matchscores[i, j] > threshold:
            # first image in pair
            im = Image.open(imlist[i])
            im.thumbnail((100, 100))
            filename = path + str(i) + '.png'
            im.save(filename)  # need temporary files of the right size
            g.add_node(pydot.Node(str(i), fontcolor='transparent', shape='rectangle', image=filename))

            im = Image.open(imlist[j])
            im.thumbnail((100, 100))
            filename = path + str(j) + '.png'
            im.save(filename)  # need temporary files of the right size
            g.add_node(pydot.Node(str(j), fontcolor='transparent', shape='rectangle', image=filename))

            g.add_edge(pydot.Edge(str(i), str(j)))
g.write_png('whitehouse.png')

因为使用的系统为mac，无法执行exe

 导致后续操作中文件夹内没有需要的sift文件，导致报错