计算机视觉--局部图像描述子：Harris角点检测算法、SIFT(尺度不变特征变换)

目录

一、Harris角点检测算法

1.1角点是什么

1.2角点算法的基本原理

1.3数学方法刻画角点特征

1.4角点响应函数

1.5Harris角点检测步骤

1.6角点检测实例

1.7图像中寻找对应点

二、SIFT(尺度不变特征变换)

2.1SIFT简介

2.2SIFT特点

2.3SIFT算法

2.4SIFT算法实现

一、Harris角点检测算法

1.1角点是什么

>轮廓之间的交点；

>对于同一场景，即使视角发生变化，通常具备稳定性质的特征；

>该点附近区域的像素点无论在梯度方向上还是其梯度幅值上有着较大变化；

1.2角点算法的基本原理 

算法基本思想是使用一个固定窗口在图像上进行任意方向上的滑动，比较滑动前与滑动后两种情况，窗口中的像素灰度变化程度，如果存在任意方向上的滑动，都有着较大灰度变化，那么我们可以认为该窗口中存在角点。

从图像局部的小窗口观察图像特征；角点定义 <-窗口向任意方向的移动都导致图像灰度的明

显变化

代码实现：

Mat RGB2GRAY(const Mat &src){
    if(!src.data || src.channels()!=3){
        exit(1);
    }
    int row = src.rows;
    int col = src.cols;
    Mat gray = Mat(row, col, CV_8UC1);
    for(int i = 0; i < row; i++){
        for(int j = 0; j < col; j++){
            gray.at<uchar>(i, j) = (uchar)(0.114 * src.at<Vec3b>(i, j)[0] + 0.587 * src.at<Vec3b>(i, j)[1] + 0.299 * src.at<Vec3b>(i, j)[2]);
        }
    }
    return gray;
}

void SobelGradDirction(Mat &src, Mat &sobelX, Mat &sobelY){
    int row = src.rows;
    int col = src.cols;
    sobelX = Mat::zeros(src.size(), CV_32SC1);
    sobelY = Mat::zeros(src.size(), CV_32SC1);

    for(int i = 0; i < row-1; i++){
        for(int j = 0; j < col-1; j++){
            double gradY = src.at<uchar>(i+1, j-1) + src.at<uchar>(i+1, j) * 2 + src.at<uchar>(i+1, j+1) - src.at<uchar>(i-1, j-1) - src.at<uchar>(i-1, j) * 2 - src.at<uchar>(i-1, j+1);
            sobelY.at<float>(i, j) = abs(gradY);
            double gradX = src.at<uchar>(i-1, j+1) + src.at<uchar>(i, j+1) * 2 + src.at<uchar>(i+1, j+1) - src.at<uchar>(i-1, j-1) - src.at<uchar>(i, j-1) * 2 - src.at<uchar>(i+1, j-1);
            sobelX.at<float>(i, j) = abs(gradX);
        }
    }
    //将梯度数组转换为8位无符号整形
    convertScaleAbs(sobelX, sobelX);
    convertScaleAbs(sobelY, sobelY);
}

Mat SobelXX(const Mat src){
    int row = src.rows;
    int col = src.cols;
    Mat_<float> dst(row, col, CV_32FC1);
    for(int i = 0; i < row; i++){
        for(int j = 0; j < col; j++){
            dst.at<float>(i, j) = src.at<uchar>(i, j) * src.at<uchar>(i, j);
        }
    }
    return dst;
}

Mat SobelYY(const Mat src){
    int row = src.rows;
    int col = src.cols;
    Mat_<float> dst(row, col, CV_32FC1);
    for(int i = 0; i < row; i++){
        for(int j = 0; j < col; j++){
            dst.at<float>(i, j) = src.at<uchar>(i, j) * src.at<uchar>(i, j);
        }
    }
    return dst;
}

Mat SobelXY(const Mat src1, const Mat &src2){
    int row = src1.rows;
    int col = src1.cols;
    Mat_<float> dst(row, col, CV_32FC1);
    for(int i = 0; i < row; i++){
        for(int j = 0; j < col; j++){
            dst.at<float>(i, j) = src1.at<uchar>(i, j) * src2.at<uchar>(i, j);
        }
    }
    return dst;
}

void separateGaussianFilter(Mat_<float> &src, Mat_<float> &dst, int ksize, double sigma){
    CV_Assert(src.channels()==1 || src.channels() == 3); //只处理单通道或者三通道图像
    //生成一维的
    double *matrix = new double[ksize];
    double sum = 0;
    int origin = ksize / 2;
    for(int i = 0; i < ksize; i++){
        double g = exp(-(i-origin) * (i-origin) / (2 * sigma * sigma));
        sum += g;
        matrix[i] = g;
    }
    for(int i = 0; i < ksize; i++) matrix[i] /= sum;
    int border = ksize / 2;
    copyMakeBorder(src, dst, border, border, border, border, BORDER_CONSTANT);
    int channels = dst.channels();
    int rows = dst.rows - border;
    int cols = dst.cols - border;
    //水平方向
    for(int i = border; i < rows - border; i++){
        for(int j = border; j < cols - border; j++){
            float sum[3] = {0};
            for(int k = -border; k<=border; k++){
                if(channels == 1){
                    sum[0] += matrix[border + k] * dst.at<float>(i, j+k);
                }else if(channels == 3){
                    Vec3f rgb = dst.at<Vec3f>(i, j+k);
                    sum[0] += matrix[border+k] * rgb[0];
                    sum[1] += matrix[border+k] * rgb[1];
                    sum[2] += matrix[border+k] * rgb[2];
                }
            }
            for(int k = 0; k < channels; k++){
                if(sum[k] < 0) sum[k] = 0;
                else if(sum[k] > 255) sum[k] = 255;
            }
            if(channels == 1)
                dst.at<float>(i, j) = static_cast<float>(sum[0]);
            else if(channels == 3){
                Vec3f rgb = {static_cast<float>(sum[0]), static_cast<float>(sum[1]), static_cast<float>(sum[2])};
                dst.at<Vec3f>(i, j) = rgb;
            }
        }
    }
    //竖直方向
    for(int i = border; i < rows - border; i++){
        for(int j = border; j < cols - border; j++){
            float sum[3] = {0};
            for(int k = -border; k<=border; k++){
                if(channels == 1){
                    sum[0] += matrix[border + k] * dst.at<float>(i+k, j);
                }else if(channels == 3){
                    Vec3f rgb = dst.at<Vec3f>(i+k, j);
                    sum[0] += matrix[border+k] * rgb[0];
                    sum[1] += matrix[border+k] * rgb[1];
                    sum[2] += matrix[border+k] * rgb[2];
                }
            }
            for(int k = 0; k < channels; k++){
                if(sum[k] < 0) sum[k] = 0;
                else if(sum[k] > 255) sum[k] = 255;
            }
            if(channels == 1)
                dst.at<float>(i, j) = static_cast<float>(sum[0]);
            else if(channels == 3){
                Vec3f rgb = {static_cast<float>(sum[0]), static_cast<float>(sum[1]), static_cast<float>(sum[2])};
                dst.at<Vec3f>(i, j) = rgb;
            }
        }
    }
    delete [] matrix;
}

Mat harrisResponse(Mat_<float> GaussXX, Mat_<float> GaussYY, Mat_<float> GaussXY, float k){
    int row = GaussXX.rows;
    int col = GaussXX.cols;
    Mat_<float> dst(row, col, CV_32FC1);
    for(int i = 0; i < row; i++){
        for(int j = 0; j < col; j++){
            float a = GaussXX.at<float>(i, j);
            float b = GaussYY.at<float>(i, j);
            float c = GaussXY.at<float>(i, j);
            dst.at<float>(i, j) = a * b - c * c - k * (a + b) * (a + b);
        }
    }
    return dst;
}

int dir[8][2] = {0, 1, 0, -1, 1, 0, -1, 0, 1, 1, 1, -1, -1, 1, -1, -1};

void MaxLocalValue(Mat_<float>&resultData, Mat srcGray, Mat &resultImage, int kSize){
    int r = kSize / 2;
    resultImage = srcGray.clone();
    int row = resultImage.rows;
    int col = resultImage.cols;
    for(int i = r; i < row - r; i++){
        for(int j = r; j < col - r; j++){
            bool flag = true;
            for(int k = 0; k < 8; k++){
                int tx = i + dir[k][0];
                int ty = j + dir[k][1];
                if(resultData.at<float>(i, j) < resultData.at<float>(tx, ty)){
                    flag = false;
                }
            }
            if(flag && (int)resultData.at<float>(i, j) > 18000){
                circle(resultImage, Point(i, j), 5, Scalar(0, 0, 255), 2, 8, 0);
            }
        }
    }
}

效果图：

 1.3数学方法刻画角点特征

将图像窗口平移 [ u,v ] 产生灰度变化 E ( u,v )

 

 

共可分为三种情况：

   （1） 图像中的直线：一个特征值大，另一个特征值小，λ1>λ2或λ2>λ1。自相关函数值在某一方向上大，在其他方向上小。

   （2）图像中的平面：两个特征值都小，且近似相等；自相关函数数值在各个方向上都小。

   （3）图像中的角点：两个特征值都大，且近似相等，自相关函数在所有方向都增大。

 1.4角点响应函数

由于我们是通过M的两个特征值的大小对图像进行分类，所以，定义角点相应函数R： 

 

 其中k为经验常数，一般取k=0.04~0.06。增大k的值将减小响应值R，降低检测的灵敏性，减少被检测角点的数量；减小k值，将增大角点响应值R，增加被检测角点的数量。为了去除加权常数κ，我们通常使用商数detM/(traceM)2作为指示器。所以，上图可以转化为：

•  R 只与M的特征值有关

               • 角点：R 为大数值正数（λ1和λ2都是很大的正数）

               • 边缘：R 为大数值负数

               • 平坦区：R 为小数值（如果λ1≈λ2≈0，该区域为空）

         在判断角点的时候，对角点响应函数R进行阈值处理：R > threshold，提取R的局部极大值。
 

1.5Harris角点检测步骤

Harris角点检测可以分为5个步骤：

（1）计算图像I(x,y)I(x,y)在xx和yy两个方向的梯度IxIx,IyIy;

（2）计算图像两个方向梯度的乘积;

（3）使用高斯函数对Ix2、Iy2、IxIy进行高斯加权(取σ=2,ksize=3)，计算中心点为(x,y)(x,y)的窗口W对应的矩阵M；

（4）计算每个像素点(x,y)处的(x,y)处的Harris响应值R； 

 （5）过滤大于某一阈值t的R值；

 1.6角点检测实例

（1）代码实现:

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

def compute_harris_response(im, sigma=3):
    # 在一幅灰度图像中，对每个像素计算Harris角点检测器响应函数
    
    # 计算导数
    k = 0.04
    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
    
    # 返回像素值为 Harris 响应函数值的一幅图像
    print(Wdet-(k*Wtr)) # 输出角点响应函数R
    return Wdet / Wtr  # 此处可消除参数k的影响
    
# 读入图像(读取图像到数组中)
im = array(Image.open('p3.jpg').convert('L'))  # 括号内：读取一幅图像，并将其转换成灰度图像

# 检测harris角点
harrisim = compute_harris_response(im)
                # compute_harris_response(im, sigma):
                #     在一幅灰度图像中，对每一个像素计算Harris角点检测器响应函数
                #     im:(数组图像）  sigma:标准差默认为3
                
# Harris响应函数
# print (harrisim)
harrisim1 = 255 - harrisim
figure()
gray()

#画出Harris响应图
subplot(141)
imshow(harrisim1)
print harrisim1.shape
axis('off')
axis('equal')

threshold = [0.01, 0.05, 0.1]
for i, thres in enumerate(threshold):
    # enumerate()函数用于将一个可遍历的数据对象(如列表、元组或字符串)组合为一个索引序列，同时列出数据和数据下标，一般用在for 循环当中
        filtered_coords = harris.get_harris_points(harrisim, 6, thres)
                                # get_harris_points(harrisim, min_dist=10, threshold=0.1):
                                #       从一幅Harris 响应图像harrisim中返回角点。min_dist 为分割角点和图像边界的最少像素数目（默认为10）。
                                #       threshold 为阀值，大于阀值的harris响应函数值被认为是可能的角点，并在harrism_t矩阵中相应的位置1,其余地方置0
        subplot(1, 4, i+2)
        imshow(im)
        print im.shape
        plot([p[1] for p in filtered_coords], [p[0] for p in filtered_coords], '*') # 标出角点
        axis('off')
        
show()

上面代码中重写compute_harris_response()函数，以便将角点响应函数R输出。

截图：

 （2）输入图片特征为纹理平坦丰富：

 1.7图像中寻找对应点

Harris 角点检测器仅仅能够检测出图像中的兴趣点，但是没有给出通过比较图像间的兴趣点来寻找匹配角点的方法。我们需要在每个点上加入描述子信息，并给出一 个比较这些描述子的方法。兴趣点描述子是分配给兴趣点的一个向量，描述该点附近的图像的表观信息。描述子越好，寻找到的对应点越好。我们用对应点或者点的对应来描述相同物体和场景点在不同图像上形成的像素点。使用归一化的互相关矩阵进行两幅图像的比较，代码如下：
 

def get_descriptors(image,filtered_coords,wid=5):
    """对于每个返回的点，返回点周围2*wid+1个像素的值（假设选取点的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):
    """对于第一幅图像中的每个角点描述子，使用归一化互相关，选取它在第二幅图像中的匹配角点"""
    n = len(desc1[0])
    # 点对的距离
    d = -np.ones((len(desc1), len(desc2)))
    for i in range(len(desc1)):
        for j in range(len(desc2)):
            # mean:求平均值，std:求标准差
            d1 = (desc1[i] - np.mean(desc1[i])) / np.std(desc1[i])
            d2 = (desc2[j] - np.mean(desc2[j])) / np.std(desc2[j])
            ncc_value = sum(d1 * d2) / (n-1)
            if ncc_value > threshold:  # 数值较高的距离代表两个点
                d[i,j] = ncc_value
    # np.argsort():将数组中的元素从小到大排列，提取其对应的index(索引)，然后输出
    # 数值较高的距离代表两个点能够更好的匹配，因此对距离取相反数
    ndx = np.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)
    # np.where:找到n维数组中特定数值的索引
    ndx_12 = np.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 appendimages(im1,im2):
    """返回将两幅图像并排拼接成一幅图像"""
    # 选取具有最少行数的图像，然后填充足够的空行
    rows1 = im1.shape[0]
    rows2 = im2.shape[0]
 
    if rows1 < rows2:
        im1 = np.concatenate((im1, np.zeros((rows2 - rows1, im1.shape[1]))), axis=0)
    elif rows1 > rows2:
        im2 = np.concatenate((im2, np.zeros((rows1 - rows2, im2.shape[1]))), axis=0)
    # 如果无此两种情况，则他们行数相同，无需填充
    return np.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 = appendimages(im1,im2)
    if show_below:
        im3 = np.vstack((im3, im3))
    plt.imshow(im3)
    cols1 = im1.shape[1]
    for i,m in enumerate(matchscores):
        if m>0:
            plt.plot([locs1[i][1],locs2[m][1]+cols1],[locs1[i][0],locs2[m][0]],'c')
    plt.axis('off')

测试代码：

# -*- coding: utf-8 -*-
"""
@author: RRJ
@software: PyCharm
@file: harrisdemo3.py
@time: 2022/3/27 15:19
"""
import numpy as np
from PIL import Image
import matplotlib.pyplot as plt
from Localdescriptors import harris
from tools.imtools import imresize
 
# 读入图像
im1 = np.array(Image.open("../JMU_img/JMU_ZS1.jpg").convert("L"))
im2 = np.array(Image.open("../JMU_img/JMU_ZS2.jpg").convert("L"))
 
# resize加快匹配速度
im1 = imresize(im1, (im1.shape[1]//2, im1.shape[0]//2))
im2 = imresize(im2, (im2.shape[1]//2, im2.shape[0]//2))
 
wid = 5
harrisim = harris.compute_harris_response(im1, 5)
filtered_coords1 = harris.get_harris_points(harrisim, wid+1)
d1 = harris.get_descriptors(im1, filtered_coords1, wid)
 
harrisim = harris.compute_harris_response(im2, 5)
filtered_coords2 = harris.get_harris_points(harrisim, wid+1)
d2 = harris.get_descriptors(im2, filtered_coords2, wid)
 
print('starting matching')
matches = harris.match_twosided(d1, d2)
 
plt.figure()
plt.gray()
harris.plot_matches(im1, im2, filtered_coords1, filtered_coords2, matches)
plt.show()

运行截图：

二、SIFT(尺度不变特征变换)

2.1SIFT简介

SIFT特征是基于物体上的一些局部外观的兴趣点而与影像的大小和旋转无关。对于光线、噪声、微视角改变的容忍度也相当高。基于这些特性，它们是高度显著而且相对容易撷取，在母数庞大的特征数据库中，很容易辨识物体而且鲜有误认。使用SIFT特征描述对于部分物体遮蔽的侦测率也相当高，甚至只需要3个以上的SIFT物体特征就足以计算出位置与方位。在现今的电脑硬件速度下和小型的特征数据库条件下，辨识速度可接近即时运算。SIFT特征的信息量大，适合在海量数据库中快速准确匹配。

2.2 SIFT特点

1.SIFT特征是图像的局部特征，其对旋转、尺度缩放、亮度变化保持不变性，对视角变化、仿射变换、噪声也保持一定程度的稳定性；

2. 区分性（Distinctiveness）好，信息量丰富，适用于在海量特征数据库中进行快速、准确的匹配；

3. 多量性，即使少数的几个物体也可以产生大量的SIFT特征向量；

4.高速性，经优化的SIFT匹配算法甚至可以达到实时的要求；

5.可扩展性，可以很方便的与其他形式的特征向量进行联合。

2.3SIFT算法

SIFT算法可以解决的问题

（1）目标的旋转、缩放、平移；

（2）图像仿射或投影变换；

（3）弱光照影响以及部分目标遮挡；

（4）杂物场景以及噪声点影响。

SIFT算法实现特征匹配的三个流程：

1.提取关键点；

2.对关键点附加详细的信息（局部特征），即描述符；

3.通过特征点（附带上特征向量的关键点）的两两比较找出相互匹配的若干对特征点，建立景物间的对应关系。

2.4SIFT算法实现

代码实现;

#pragma once                    //防止头文件重复包含和下面作用一样
 
#include<iostream>
#include<vector>
#include<opencv2\core\core.hpp>
#include<opencv2\features2d\features2d.hpp>
 
using namespace std;
using namespace cv;
 
/*************************定义常量*****************************/
 
//高斯核大小和标准差关系，size=2*(GAUSS_KERNEL_RATIO*sigma)+1,经常设置GAUSS_KERNEL_RATIO=2-3之间
const double GAUSS_KERNEL_RATIO = 3;
 
const int MAX_OCTAVES = 8;					//金字塔最大组数
 
const float CONTR_THR = 0.04f;				//默认是的对比度阈值(D