Python计算机视觉——第二章局部图像描绘子

最新推荐文章于 2023-07-20 15:17:04 发布

Eaton_chen

最新推荐文章于 2023-07-20 15:17:04 发布

阅读量472

点赞数

分类专栏：计算机视觉文章标签： Harris角点检测 SIFT 特征匹配尺度不变性图像处理

本文链接：https://blog.csdn.net/weixin_44395365/article/details/118161999

版权

计算机视觉专栏收录该内容

10 篇文章 8 订阅

订阅专栏

文章目录

2.1 Harris角点检测器的原理
- 2.1.1 实验部分
2.2 SIFT（尺度不变特征变换）
- 2.2.1 实验部分
2.3 SIFT匹配描述子

2.1 Harris角点检测器的原理

在图像中，如果像素的周围显示存在多于一个方向的边，我们就认为该店为兴趣点，也称为角点。

角点具有如下的特点：

轮廓之间的交点；
对于同一场景，即使视角发生变化，通常具备稳定性质的特征；
该点附近区域的像素点无论在梯度方向上还是其梯度幅值上有着较大变化；

我们可以由一个函数来判断是否是角点：
$E(u,v)=\sum _{x,y} w(x,y)[I(x+u,y+v)-I(x,y)]^2$ w(x,y)表示矩形窗口，I(x,y)表示该点像素的强度，u与v表示窗口的偏移量，若在任意方向都存在着较大灰度的变化，则判断该点是角点。我们可以使用高斯权重矩阵来突出变化较大的点，赋予更大的权值。

我们可以使用一阶泰勒公式来简化上述的函数：
$E(u,v)\approx \sum|I(x,y)+uI_x+vI_y-I(x,y)|^2=\sum u^2I_x^2+2uvI_xI_y+v^2I_y^2$

用矩阵表示： $v](\sum[\left[ \begin{matrix} I_x^2 & I_xI_y\\ \\ I_xI_y\ & I_y^2 \end{matrix} \right]])[u,v]^T$

令矩阵M为中间的矩阵，求解特征值 $|\lambda E-M|=0$ ,可得特征值1= $I_x^2+I_y^2$ ,特征值2为0。

卷积： $\bar{M_t}=W*M$ W是高斯滤波器的权重函数,W的宽度决定了像素X周围的感兴趣区域。

我们可以通过两张图来总结一下特征值和角点的关系
在这里插入图片描述
角点区域的两个特征值都是很大的，而边缘部分是一大一小，平滑部分两个特征值均很小。

2.1.1 实验部分

Harris算法，每个函数的注释均已写出。

# 引入所有函数，绘图板块
from pylab import *
from numpy import *
#n维图像库，我们这里使用filters函数
from scipy.ndimage import filters

# 使用高斯滤波器进行卷积，标准差为3
# 这个函数得到Harris响应函数值的图像
def compute_harris_response(im,sigma=3):
    """ Compute the Harris corner detector response function 
        for each pixel in a graylevel image. """
    
    # derivatives
    imx = zeros(im.shape)
    # x方向上的高斯导数
    filters.gaussian_filter(im, (sigma,sigma), (0,1), imx)
    imy = zeros(im.shape)
    # y方向上的高斯导数
    filters.gaussian_filter(im, (sigma,sigma), (1,0), imy)
    
    # compute components of the Harris matrix
    # 导数运算的高斯模糊值，第一个参数表示input，通过这样的方式得到矩阵的分量
    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
# 从响应图像中获得harris角点检测结果，min_dist表示最少数目的分割角点，threshold是阈值 
def get_harris_points(harrisim,min_dist=10,threshold=0.1):
    """ Return corners from a Harris response image
        min_dist is the minimum number of pixels separating 
        corners and image boundary. """
    
    # find top corner candidates above a threshold
    # 寻找高于阈值的候选角点，.max是numpy的函数
    corner_threshold = harrisim.max() * threshold
    harrisim_t = (harrisim >corner_threshold) * 1
    
    # get coordinates of candidates
    #nonzeros(a)返回数组a中值不为零的元素的下标，它的返回值是一个长度为a.ndim(数组a的轴数)的元组,.T是矩阵转置操作
    coords = array(harrisim_t.nonzero()).T
    
    # ...and their values，harrisim的响应值
    candidate_values = [harrisim[c[0],c[1]] for c in coords]
    
    # 从小到大输出，注意输出的是下标
    index = argsort(candidate_values)
    
    # 标记分割角点范围的坐标
    allowed_locations = zeros(harrisim.shape)
    allowed_locations[min_dist:-min_dist,min_dist:-min_dist] = 1
    
    # select the best points taking min_distance into account,选择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
    
# 绘制检测角点，灰度图 ,上个函数的角点需要满足的条件，1：角点范围之内；2：高于harrisim响应值的阈值   
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],'ro')
    axis('off')
    show()

# 读入图像
im = array(Image.open('pic/empire.jpg').convert('L'))
 
# 检测harris角点
harrisim = compute_harris_response(im)
 
figure()
imshow(harrisim)
axis('off')
threshold = [0.01, 0.05, 0.1]
for i, thres in enumerate(threshold):
    filtered_coords = get_harris_points(harrisim, 6, thres)
    plot_harris_points(im,filtered_coords)
show()

在这里插入图片描述

实验分析：不同的阈值，检测出来角点的数量会不同，阈值区域的判定是用.max()再乘以threshold，第二张图片的阈值是0.01，相当的小，符合该阈值的区域应该相当得大，显示出来的被标记的红色圆圈是检测到的Harris角点，由于Harris检测器检测原理很简单，在稳健性上能力不足，有一些点可能并不准确，随着阈值的增大，检测出来的角点准确性将升高。

实验第二部分：通过比较各兴趣点，找到匹配的角点，方法是在每个点上加入描绘子信息，具体见下面的代码，注释均已写出。

# 每个点加入描述子信息，比较兴趣点来寻找匹配的角点，通过互相关矩阵，向量乘积越大，相似度越高，对于每个返回的点，返回点周围2*wid+1个像素的值
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)
    # matches_12>=0的部分不变，其他向量中的元素置0
    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]
    
    # 依据rows1/rows2的长度，选择填充的方式，纵轴方向
    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]
    # i是下标，m是元素
    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')


im1 = array(Image.open("pic/crans_1_small.jpg").convert("L"))
im2 = array(Image.open("pic/crans_2_small.jpg").convert("L"))

wid = 5
harrisim = compute_harris_response(im1, 5)
filtered_coords1 = get_harris_points(harrisim, wid+1)
d1 = get_descriptors(im1, filtered_coords1, wid)
 
harrisim =compute_harris_response(im2, 5)
filtered_coords2 =get_harris_points(harrisim, wid+1)
d2 = get_descriptors(im2, filtered_coords2, wid)
 
matches = match_twosided(d1, d2)
 
figure()
gray() 
harris.plot_matches(im1, im2, filtered_coords1, filtered_coords2, matches)
show()

在这里插入图片描述
实验分析：由于算法采用互相关矩阵的判定方式，具有较弱的稳定性，所以当两张不同视角的图片进行比对时，会出现一些不正确的匹配，在图片中可以显示为蓝色直线相对应物体的不同部分，通常采取一些更加稳健的方式来描述匹配。

2.2 SIFT（尺度不变特征变换）

SIFT是过去十年中最成功的图像局部描绘子之一，SIFT在兴趣点检测器和描述子中都具有很强的稳健性，SIFT特征对于尺度、旋转和亮度都具有不变性，因此可以用于三维视角和噪声的可靠匹配。

Lowe将SIFT算法分解为如下四步：

尺度空间极值检测：搜索所有尺度上的图像位置。通过高斯微分函数来识别潜在的对于尺度和旋转不变的兴趣点。
关键点定位：在每个候选的位置上，通过一个拟合精细的模型来确定位置和尺度。关键点的选择依据于它们的稳定程度。
方向确定：基于图像局部的梯度方向，分配给每个关键点位置一个或多个方向。所有后面的对图像数据的操作都相对于关键点的方向、尺度和位置进行变换，从而提供对于这些变换的不变性。
关键点描述：在每个关键点周围的邻域内，在选定的尺度上测量图像局部的梯度。这些梯度被变换成一种表示，这种表示允许比较大的局部形状的变形和光照变化。

尺度空间：在图像信息处理模型中引入一个被视为尺度的参数，通过连续变化尺度参数获得多尺度下的尺度空间表示序列，对这些序列进行尺度空间主轮廓的提取，并以该主轮廓作为一种特征向量，实现边缘、角点检测和不同分辨率上的特征提取等.

$L(x,y,\sigma)=G(x,y,\sigma)*I(x,y)$ 其中 $G(x,y,\sigma)=\frac{1}{2 \pi \sigma^2}e-\frac{(x-m/2)^2+(y-n/2)^2}{2\sigma ^2}$

尺度空间用高斯金字塔来描述，用来检测稳定的关键点。

对图像做不同尺度的高斯模糊；
对图像做降采样(隔点采样)。

SIFT利用高斯差分函数来定位兴趣点：

$D(x,\sigma)=[G_{k\sigma}(x)-G_\sigma(x)*I(x)]$ k是决定相差尺度的常数。利用差分来代替微分： $\frac{\alpha G}{\alpha \sigma}=\frac{G(x,y,k\sigma)-G(x,y,\sigma)}{k\sigma-\sigma}$

关键点是由DOG空间的局部极值点组成的，关键点的初步探查是通过同一组内各DoG相邻两层图像之间比较完成的,精确定位是通过拟合函数确定，这里略。

方向确定：我们通过计算梯度和方向，通过直方图的统计，确定主方向。
在这里插入图片描述
SIFT描述子是关键点邻域高斯图像梯度统计结果的一种表示。通过对关键点周围图像区域分块，计算块内梯度直方图，生成具有独特性的向量，这个向量是该区域图像信息的一种抽象，具有唯一性。

2.2.1 实验部分

注意要装好VLfeat的sift.exe和VL.dll，解压在python文件中。

from PIL import Image
import os
from numpy import *
from pylab import *


def process_image(imagename,resultname,params="--edge-thresh 10 --peak-thresh 5"):
    """ Process an image and save the results in a file. """

    if imagename[-3:] != 'pgm':
        # create a pgm file
        im = Image.open(imagename).convert('L')
        im.save('tmp.pgm')
        imagename = 'tmp.pgm'

    cmmd = str("sift.exe "+imagename+" --output="+resultname+
            " "+params)
    # 将字符串转化成命令在服务器上运行
    os.system(cmmd)
    print 'processed', imagename, 'to', resultname


def read_features_from_file(filename):
    """ Read feature properties and return in matrix form. """
    #读取文件
    f = loadtxt(filename)
    # 返回特征位置，描述子
    return f[:,:4],f[:,4:] # feature locations, descriptors

def write_features_to_file(filename,locs,desc):
    """ Save feature location and descriptor to file. """
    # 水平叠堆两个向量
    savetxt(filename,hstack((locs,desc)))
    

def plot_features(im,locs,circle=False):
    """ Show image with features. input: im (image as array), 
        locs (row, col, scale, orientation of each feature). """

    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)
    # 将参数circle的选项设置为True
    if circle:
        for p in locs:
            draw_circle(p[:2],p[2]) 
    else:
        plot(locs[:,0],locs[:,1],'ob')
    axis('off')

imname='pic/empire.jpg'
im1=array(Image.open(imname).convert('L'))
process_image(imname,'empire.sift')
l1,d1=read_features_from_file('empire.sift')

figure()
gray()
plot_features(im1,l1,circle=True)
show()

在这里插入图片描述

实验分析：上图是SIFT特征检测，下图是对照的Harris检测器，我们可以注意到第一张图的圆圈大小不同，为什么会出现这样的情况呢，通过程序的代码的算法思想，read_features_from_file这个函数所返回的参数是特征位置和描述子，二进制文件中前4个数值依次表示兴趣点的坐标、尺度和方向角度，由于尺度的不同，导致了圆圈大小不一样，所以SIFT特征检测比起Harris更具有描述子上的优势，特征点的位置检测地不同，算法更精确、稳健。

2.3 SIFT匹配描述子

SIFT的特征匹配，采取一种稳健的准则，使用两个特征距离和两个最匹配特征距离的比率，相比于其他特征，该准则能保证找到足够相似的唯一特征，使错误的匹配数降低。

def match(desc1,desc2):
    """ For each descriptor in the first image, 
        select its match in the second image.
        input: desc1 (descriptors for the first image), 
        desc2 (same for second image). """
    #linalg.norm(d)是范数
    desc1 = array([d/linalg.norm(d) for d in desc1])
    desc2 = array([d/linalg.norm(d) for d in desc2])
    
    dist_ratio = 0.6
    desc1_size = desc1.shape
    
    matchscores = zeros((desc1_size[0]),'int')
    desc2t = desc2.T # precompute matrix transpose
    for i in range(desc1_size[0]):
        dotprods = dot(desc1[i,:],desc2t) # vector of dot products
        dotprods = 0.9999*dotprods
        # inverse cosine and sort, return index for features in second image
        indx = argsort(arccos(dotprods))
        
        # check if nearest neighbor has angle less than dist_ratio times 2nd
        if arccos(dotprods)[indx[0]] < dist_ratio * arccos(dotprods)[indx[1]]:
            matchscores[i] = int(indx[0])
    
    return matchscores

#进一步增强稳健性

def match_twosided(desc1,desc2):
    """ Two-sided symmetric version of match(). """
    
    matches_12 = match(desc1,desc2)
    matches_21 = match(desc2,desc1)
    
    ndx_12 = matches_12.nonzero()[0]
    
    # remove matches that are not symmetric
    for n in ndx_12:
        if matches_21[int(matches_12[n])] != n:
            matches_12[n] = 0
    
    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]
    
    # 依据rows1/rows2的长度，选择填充的方式，纵轴方向
    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]
    # i是下标，m是元素
    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')

from PIL import Image
import os
from numpy import *
from pylab import *

imname1='pic/climbing_1_small.jpg'
im1=array(Image.open(imname))
process_image(imname1,'climbing1.sift')
l1,d1=read_features_from_file('climbing1.sift')



imname2='pic/climbing_2_small.jpg'
im2=array(Image.open(imname2))
process_image(imname2,'climbing2.sift')
l2,d2=read_features_from_file('climbing2.sift')


# d是描述子
matches = match_twosided(d1, d2)

figure(figsize=(10,10))
gray()
axis('equal')
#l是特征位置
plot_matches(im1, im2, l1, l2, matches, show_below=True)
show()