Harris角点检测
一、原理
角点，一般是边缘的交点，在图像匹配中经常被拿来作特征描述点。
从 Moravec 检测（最早的角点检测算法之一）出发，首先角点是局部特征，因此检测过程发生在局部窗口内，通过局部窗口在多个方向的移动 shift，观察亮度的变化，考虑三种情况：
1.平坦区域，多个方向移动，平均亮度几乎没变化 （灰度与亮度正相关）
2.边缘，某个方向上亮度变化剧烈，而垂直的方向上变化很小
3.角点，任意方向移动，平均亮度变化都很剧烈
这种方法的缺点：
滑动窗口缺点：窗口滑动只有8个方向，当边缘角落的角度不在这8个方向上则检测不准。
Harris角点检测法使用特征值的方式，使得任何方向上的角点都可以被检测出来。
设图像函数为w，像素坐标为（x,y），滑动窗口的滑动变量为（u,v）。
泰勒回顾：
根据泰勒展开，有
（横切面Q相对与Q=Y=AΔx²+BΔy²+2CΔxΔy）
二、特殊情况：
三、通常情况
四、总结
五、步骤
六、demo
七、实验
7.1 参数影响
α 从 0.04 增加到 0.06，检测到的点从 92 个降低到 78 个
7.2 灰度变化实验
I = 0.5 * I
I = I - 50
可见，检测到的点基本不变，Harris 对灰度变化不敏感，具有一定的亮度不变性。
7.3 旋转不变性实验
之前分类讨论过， E 在某个值得横切面是一个椭圆，椭圆转动时，变化的是特征向量，而特征值大小不变，所以，Harris 具有旋转不变性，如下图
逆时针旋转 90 度，检测到的点和之前都是 92 个，且相对位置不变。
但如果旋转导致图像内容发生了改变，减少了一部分，则那部分的点是没法检测出来的，如下：
7.4 尺度不变性实验
看起来好像，检测到的角点减少了，但是这也是由尺度和图像内容决定的，在高分辨率时的弱角点，到低分辨率可能就是响应强的角点。
可见，Harris 不具有尺度不变性，检测到的点个数明显有变化，但是响应较强的点还是检测出来了，不会突然失效之类的。
7.5 噪声
添加方差 0.5 的高斯噪声的 harris 检测结果，检测到 105 个点
如果噪声级别增大，检测到的角点也会增多：
添加方差 1.5 的高斯噪声的 harris 检测结果，检测到 203 个点
7.6 边缘
前面提过，Harris 也可以用于检测边缘，响应值 < 0 的即可，以下是验证结果
响应值 &amp;lt; -3
响应值 &amp;lt; -100
响应值 &amp;lt; -1000
可见，Harris 在边缘这块效果不大好，也可能是我实现的错误。
八、扩展
8.1 Shi-Tomasi 角点检测
Shi-Tomasi 角点检测结果
看起来效果不错，倒是没看出来比 Harris 好多少。个人以为，速度上，可能略微慢点，因为涉及到了特征值的计算，相对 harris 的简单加减乘除要复杂。
个人实验时发现，Shi-Tomasi 相比 harris，少了一个控制变量α ，这个是经验变量。最主要的是，判断 harris 响应的阈值，由于存在平方项，所以阈值会很高，十万百万级别，很敏感；而 Shi-Tomasi 都是一次，所以阈值会低的多，几十几百以内即可，更方便选取最佳阈值。
8.2 多尺度 Harris-Laplace
waiting
8.3 仿射变换 Harris-Affine
waiting
九、应用
9.1. 图像拼接
waiting
十、代码
完整 C++ 代码
//C++
#include <cmath>
#include <chrono>
#include <vector>
#include <cstdlib>
#include <iostream>
// opencv
#include <opencv2/core.hpp>
#include <opencv2/highgui.hpp>
#include <opencv2/imgproc.hpp>
// Eigen3
#include <Eigen/Core>
#include <Eigen/Dense>
namespace {
void run(const std::function<void()>& work=[]{}, const std::string message="") {
auto start = std::chrono::steady_clock::now();
work();
auto finish = std::chrono::steady_clock::now();
auto duration = std::chrono::duration_cast<std::chrono::milliseconds>(finish - start);
std::cout << message << " " << duration.count() << " ms" << std::endl;
}
void cv_info(const cv::Mat& one_image) {
std::cout << "高 : " << one_image.rows << "\n宽 : " << one_image.cols << "\n通道 : " << one_image.channels() << std::endl;
std::cout << "步长 : " << one_image.step << std::endl;
}
void cv_show(const cv::Mat& one_image, const char* info="") {
cv::imshow(info, one_image);
cv::waitKey(0);
cv::destroyAllWindows();
}
cv::Mat make_pad(const cv::Mat& one_image, const int pad_H, const int pad_W) {
cv::Mat padded_image;
cv::copyMakeBorder(one_image, padded_image, pad_H, pad_H, pad_W, pad_W, cv::BORDER_REPLICATE);
return padded_image;
}
cv::Mat cv_concat(const std::vector<cv::Mat> images, const bool v=false) {
cv::Mat result;
if(not v) cv::hconcat(images, result);
else cv::vconcat(images, result);
return result;
}
bool cv_write(const cv::Mat& source, const std::string save_path) {
return cv::imwrite(save_path, source, std::vector<int>({cv::IMWRITE_PNG_COMPRESSION, 0}));
}
template<typename T>
cv::Mat toint8(const std::vector<T>& source, const int H, const int W) {
cv::Mat result(H, W, CV_8UC1);
const int length = H * W;
for(int i = 0;i < length; ++i) result.data[i] = std::abs(source[i]);
return result;
}
cv::Mat get_rotated(const cv::Mat& source, const int angle, const cv::Size& _size, const cv::Point2f& center) {
cv::Mat rotated_image;
cv::Mat rot_mat = cv::getRotationMatrix2D(center, angle, 1.0);
cv::warpAffine(source, rotated_image, rot_mat, _size, cv::INTER_LINEAR);
return rotated_image;
}
cv::Mat my_rotate(const cv::Mat& source) {
const int H = source.rows;
const int W = source.cols;
cv::Mat res(W, H, CV_8UC3);
for(int i = 0;i < H; ++i) {
for(int j = 0;j < W; ++j) {
res.at<cv::Vec3b>(W - 1 - j, i)[0] = source.at<cv::Vec3b>(i, j)[0];
res.at<cv::Vec3b>(W - 1 - j, i)[1] = source.at<cv::Vec3b>(i, j)[1];
res.at<cv::Vec3b>(W - 1 - j, i)[2] = source.at<cv::Vec3b>(i, j)[2];
}
}
return res;
}
// 代码取自 https://blog.csdn.net/qq_34784753/article/details/69379135
double generateGaussianNoise(double mu, double sigma) {
const double epsilon = std::numeric_limits<double>::min();
static double z0, z1;
static bool flag = false;
flag = !flag;
if (!flag) return z1 * sigma + mu;
double u1, u2;
do {
u1 = std::rand() * (1.0 / RAND_MAX);
u2 = std::rand() * (1.0 / RAND_MAX);
} while (u1 <= epsilon);
z0 = std::sqrt(-2.0 * log(u1)) * cos(2 * CV_PI * u2);
z1 = std::sqrt(-2.0 * log(u1)) * sin(2 * CV_PI * u2);
return z0 * sigma + mu;
}
// 为图像添加高斯噪声
// 代码取自 https://blog.csdn.net/qq_34784753/article/details/69379135
cv::Mat add_gaussian_noise(const cv::Mat &source) {
cv::Mat res = source.clone();
int channels = res.channels();
int rows_number = res.rows;
int cols_number = res.cols * channels;
if (res.isContinuous()) {
cols_number *= rows_number;
rows_number = 1;
}
for (int i = 0; i < rows_number; i++) {
for (int j = 0; j < cols_number; j++) {
int val = res.ptr<uchar>(i)[j] +
generateGaussianNoise(2, 1.5) * 16;
res.ptr<uchar>(i)[j] = cv::saturate_cast<uchar>(val);
}
}
return res;
}
// 画一下 harris 响应值的热力图, 可以参考之前何恺明去雾的那次代码
}
// 高斯模糊
void gausssi_filter(int* const source, double* const target, const int radius, const int H, const int W) {
// 首先创建一维的高斯模板
const double sigma = radius * 1. / 3;
const int kernel_size = (radius << 1) + 1;
std::vector<double> weights(kernel_size, 0);
for(int i = -radius;i <= radius; ++i)
weights[i + radius] = 1. / sigma * std::exp(-i * i / (2 * sigma * sigma));
double weights_sum = 0.0;
for(int i = 0;i < kernel_size; ++i) weights_sum += weights[i];
if(std::abs(weights_sum - 0) > 1e-10) for(int i = 0;i < kernel_size; ++i) weights[i] /= weights_sum;
// 先开始 x 方向的高斯卷积
std::vector<double> temp(H * W);
const int W2 = W - radius;
for(int i = 0;i < H; ++i) {
int* const row_ptr = source + i * W;
double* const res_ptr = temp.data() + i * W;
for(int j = radius; j < W2; ++j) {
double cur_sum = 0.0;
for(int k = -radius; k <= radius; ++k)
cur_sum += weights[radius + k] * row_ptr[j + k];
res_ptr[j] = cur_sum;
}
}
// 继续做 y 方向上的高斯卷积
const int H2 = H - radius;
for(int j = radius; j < W2; ++j) {
for(int i = radius; i < H2; ++i) {
double cur_sum = 0.0;
for(int k = -radius; k <= radius; ++k)
cur_sum += weights[radius + k] * temp[i * W + j + k];
target[i * W + j] = cur_sum;
}
}
}
void box_filter(int* const source, double* const target, const int radius, const int H, const int W) {
// 储存这一行的结果
const int kernel_size = (radius << 1) + 1;
std::vector<double> buffer(W);
// 遍历前 kernel_size 行, 计算每一列的和
for(int i = 0; i < kernel_size; ++i) {
int* const row_ptr = source + i * W;
for(int j = 0; j < W; ++j) buffer[j] += row_ptr[j];
}
// 计算剩下每一行的
const int H2 = H - radius;
const int W2 = W - 2 * radius;
for(int i = radius;i < H2; ++i) {
// 计算第一个位置的和
double cur_sum = 0;
for(int j = 0;j < kernel_size; ++j) cur_sum += buffer[j];
// 记录这第一个位置的结果
const int _beg = i * W + radius;
target[_beg] = cur_sum;
// 开始向右挪动
for(int j = 1; j < W2; ++j) {
cur_sum = cur_sum - buffer[j - 1] + buffer[j - 1 + kernel_size];
target[_beg + j] = cur_sum;
}
// 这一行移动完毕, 换到下一行, 更新 buffer
if(i != H2 - 1) {
int* up_ptr = source + (i - radius) * W;
int* down_ptr = source + (i + radius + 1) * W;
for(int j = 0;j < W; ++j) buffer[j] = buffer[j] - up_ptr[j] + down_ptr[j];
}
}
const int length = H * W;
const int area = kernel_size * kernel_size;
for(int i = 0;i < length; ++i) target[i] /= area;
}
// 给定图像 source, 计算 Ix^2, Iy^2, IxIy 的加权结果
std::tuple< std::vector<double>, std::vector<double>, std::vector<double> > compute_weighted_IxIy(
const cv::Mat& source, const int radius=2, const bool use_gaussi=false) {
// 获取图像信息
const int H = source.rows;
const int W = source.cols;
const int length = H * W;
// 先获取 x, y 两个方向上的梯度, 这里用的是 Sobel
std::vector<int> gradients_x(length, 0);
std::vector<int> gradients_y(length, 0);
const int H_1 = H - 1, W_1 = W - 1;
for(int i = 1; i < H_1; ++i) {
const uchar* const row_ptr = source.data + i * W;
int* const x_ptr = gradients_x.data() + i * W;
int* const y_ptr = gradients_y.data() + i * W;
for(int j = 1; j < W_1; ++j) {
// 计算 Sobel 梯度
x_ptr[j] = 2 * row_ptr[j + 1] + row_ptr[j + 1 + W] + row_ptr[j + 1 - W] - (2 * row_ptr[j - 1] + row_ptr[j - 1 + W] + row_ptr[j - 1 - W]);
y_ptr[j] = 2 * row_ptr[j + W] + row_ptr[j + W + 1] + row_ptr[j + W - 1] - (2 * row_ptr[j - W] + row_ptr[j - W + 1] + row_ptr[j - W - 1]);
}
}
// 计算 xx, yy, xy
std::vector<int> gradients_xx(length, 0), gradients_yy(length, 0), gradients_xy(length, 0);
for(int i = 0;i < length; ++i) gradients_xx[i] = gradients_x[i] * gradients_x[i];
for(int i = 0;i < length; ++i) gradients_yy[i] = gradients_y[i] * gradients_y[i];
for(int i = 0;i < length; ++i) gradients_xy[i] = gradients_x[i] * gradients_y[i];
// 计算每一个点的加权之和, 先储存起来, 窗口函数可以是均值, 也可以是高斯加权
std::vector<double> xx_sum(length, 0), yy_sum(length, 0), xy_sum(length, 0);
const auto window_function = use_gaussi ? gausssi_filter : box_filter;
window_function(gradients_xx.data(), xx_sum.data(), radius, H, W);
window_function(gradients_yy.data(), yy_sum.data(), radius, H, W);
window_function(gradients_xy.data(), xy_sum.data(), radius, H, W);
return std::make_tuple(xx_sum, yy_sum, xy_sum);
}
using key_points_type = std::vector< std::tuple<double, int, int> >;
key_points_type harris_corner_detect(
const cv::Mat& source,
const int radius=2,
const double alpha=0.04,
const double threshold=1e5,
const int point_num=-1,
const bool use_gaussi=false) {
// 获取图像信息
const int H = source.rows;
const int W = source.cols;
const int length = H * W;
// 计算 Ix, Iy, IxIy 的加权结果
const auto gradients_info = compute_weighted_IxIy(source, radius, use_gaussi);
// 开始计算每一个点的 harris 响应值
std::vector<double> R(length, 0);
const int H_radius = H - radius;
const int W_radius = W - radius;
for(int i = radius; i < H_radius; ++i) {
const double* const xx = std::get<0>(gradients_info).data() + i * W;
const double* const yy = std::get<1>(gradients_info).data() + i * W;
const double* const xy = std::get<2>(gradients_info).data() + i * W;
double* const res_ptr = R.data() + i * W;
for(int j = radius; j < W_radius; ++j) {
// 计算这个点所在窗口的加权和
const double A = xx[j] / 255;
const double B = yy[j] / 255;
const double C = xy[j] / 255;
// 计算 λ1 和 λ2
const double det = A * B - C * C;
const double trace = A + B;
res_ptr[j] = det - alpha * (trace * trace);
}
}
// 准备一个结果
key_points_type detection;
// 需要进行局部非极大化抑制
const int H_1 = H - 1, W_1 = W - 1;
for(int i = 1; i < H_1; ++i) {
double* row_ptr = R.data() + i * W;
for(int j = 1; j < W_1; ++j) {
const double center = row_ptr[j];
if(center > threshold) {
if(center > row_ptr[j - 1] and center > row_ptr[j + 1] and
center > row_ptr[j - 1 - W] and center > row_ptr[j - W] and center > row_ptr[j + 1 - W] and
center > row_ptr[j - 1 + W] and center > row_ptr[j + W] and center > row_ptr[j + 1 + W])
detection.emplace_back(center, i, j);
}
}
}
// 取前 point_sum 个
if(point_num > 0 and detection.size() > point_num) {
// 按照响应值大小排序
std::sort(detection.begin(), detection.end());
std::reverse(detection.begin(), detection.end());
detection.erase(detection.begin() + point_num, detection.end());
detection.shrink_to_fit();
}
std::cout << "收集到 " << detection.size() << " 个角点 " << std::endl;
return detection;
}
// Shi - Tomasi 的角点检测
key_points_type shi_tomasi_corner_detect(
const cv::Mat& source,
const int radius=2,
const double threshold=1e5,
const int point_num=-1,
const bool use_gaussi=false) {
// 获取图像信息
const int H = source.rows;
const int W = source.cols;
const int length = H * W;
// 计算 Ix, Iy, IxIy 的加权结果
const auto gradients_info = compute_weighted_IxIy(source, radius, use_gaussi);
// 开始计算每一个点的 harris 响应值
std::vector<double> R(length, 0);
const int H_radius = H - radius;
const int W_radius = W - radius;
for(int i = radius; i < H_radius; ++i) {
const double* const xx = std::get<0>(gradients_info).data() + i * W;
const double* const yy = std::get<1>(gradients_info).data() + i * W;
const double* const xy = std::get<2>(gradients_info).data() + i * W;
double* const res_ptr = R.data() + i * W;
for(int j = radius; j < W_radius; ++j) {
// 构建 M 矩阵
Eigen::MatrixXd M(2, 2);
M << xx[j] / 255, xy[j] / 255, xy[j] / 255, yy[j] / 255;
// 求 M 的特征值
Eigen::EigenSolver<Eigen::MatrixXd> solver(M);
const auto& result = solver.eigenvalues().real();
res_ptr[j] = std::min(result[0], result[1]);
}
}
// 准备一个结果
key_points_type detection;
// 需要进行局部非极大化抑制
const int H_1 = H - 1, W_1 = W - 1;
for(int i = 1; i < H_1; ++i) {
double* row_ptr = R.data() + i * W;
for(int j = 1; j < W_1; ++j) {
const double center = row_ptr[j];
if(center > threshold) {
if(center > row_ptr[j - 1] and center > row_ptr[j + 1] and
center > row_ptr[j - 1 - W] and center > row_ptr[j - W] and center > row_ptr[j + 1 - W] and
center > row_ptr[j - 1 + W] and center > row_ptr[j + W] and center > row_ptr[j + 1 + W])
detection.emplace_back(center, i, j);
}
}
}
// 取前 point_sum 个
if(point_num > 0 and detection.size() > point_num) {
// 按照响应值大小排序
std::sort(detection.begin(), detection.end());
std::reverse(detection.begin(), detection.end());
detection.erase(detection.begin() + point_num, detection.end());
detection.shrink_to_fit();
}
std::cout << "收集到 " << detection.size() << " 个角点 " << std::endl;
return detection;
}
void demo_1() {
const std::string save_dir("./images/output/1/");
std::string origin_path("../images/detail/harris_demo_1.png"); // harris_demo_1.png
const auto origin_image = cv::imread(origin_path);
if(origin_image.empty()) {
std::cout << "读取图像 " << origin_path << " 失败 !" << std::endl;
return;
}
// 转成灰度图
cv::Mat origin_gray;
cv::cvtColor(origin_image, origin_gray, cv::COLOR_BGR2GRAY);
// 写一个展示用的函数
auto show_harris = [](const cv::Mat& source, const key_points_type& harris_result, const std::string save_name, const int radius=2, const int thickness=4)
-> void {
// 画出来
cv::Mat display = source.clone();
for(const auto& item : harris_result)
cv::circle(display, cv::Point(std::get<2>(item), std::get<1>(item)), radius, cv::Scalar(0, 0, 255), thickness);
// 展示
cv_show(display);
cv_write(display, save_name);
};
// 检测 Harris 角点
auto harris_result = harris_corner_detect(origin_gray, 2, 0.04, 4e3, 0);
show_harris(origin_image, harris_result, save_dir + "original.png");
const int best_size = harris_result.size();
// 增大局部窗口的半径
harris_result = harris_corner_detect(origin_gray, 5, 0.04, 4e3, 0);
show_harris(origin_image, harris_result, save_dir + "original_radius.png");
// 增加参数 α
harris_result = harris_corner_detect(origin_gray, 2, 0.06, 4e3, 0);
show_harris(origin_image, harris_result, save_dir + "original_alpha.png");
// 证明对亮度或对比度改变的影响
cv::Mat darken_image = 0.5 * origin_image, darken_gray;
darken_image.convertTo(darken_image, CV_8UC3);
cv::cvtColor(darken_image, darken_gray, cv::COLOR_BGR2GRAY);
harris_result = harris_corner_detect(darken_gray, 2, 0.04, 2e2, best_size);
show_harris(darken_image, harris_result, save_dir + "darken.png");
cv::Mat darken_image_2 = origin_image - 50, darken_gray_2;
darken_image_2.convertTo(darken_image_2, CV_8UC3);
cv::cvtColor(darken_image, darken_gray_2, cv::COLOR_BGR2GRAY);
harris_result = harris_corner_detect(darken_gray_2, 2, 0.04, 2e2, best_size);
show_harris(darken_image_2, harris_result, save_dir + "darken_2.png");
// 证明旋转不变性
const auto my_rotated_image = my_rotate(origin_image);
cv::Mat my_rotated_gray;
cv::cvtColor(my_rotated_image, my_rotated_gray, cv::COLOR_BGR2GRAY);
harris_result = harris_corner_detect(my_rotated_gray, 2, 0.04, 4e3, 0);
show_harris(my_rotated_image, harris_result, save_dir + "rotated_1.png");
// 内容有损的旋转
const int H = origin_image.rows;
const int W = origin_image.cols;
const auto rotated_image = get_rotated(origin_image, 45, cv::Size(W, H), cv::Point2f(W / 2, H / 2));
cv::Mat rotated_gray;
cv::cvtColor(rotated_image, rotated_gray, cv::COLOR_BGR2GRAY);
harris_result = harris_corner_detect(rotated_gray, 2, 0.04, 4e3, 0);
show_harris(rotated_image, harris_result, save_dir + "rotated_2.png");
// 镜像
cv::Mat flipped_image, flipped_gray;
cv::flip(origin_image, flipped_image, 1);
cv::cvtColor(flipped_image, flipped_gray, cv::COLOR_BGR2GRAY);
harris_result = harris_corner_detect(flipped_gray, 2, 0.04, 4e3, 0);
show_harris(flipped_image, harris_result, save_dir + "flip_1.png");
// 证明不满足 尺度不变性
cv::Mat resized_image, resized_gray;
cv::resize(origin_image, resized_image, cv::Size(W / 8, H / 8));
cv::cvtColor(resized_image, resized_gray, cv::COLOR_BGR2GRAY);
harris_result = harris_corner_detect(resized_gray, 2, 0.04, 4e3, 0);
show_harris(resized_image, harris_result, save_dir + "scaled.png", 1, 1);
// 考虑抗噪性
const auto noise_image = add_gaussian_noise(origin_image);
cv::Mat noise_gray;
cv::cvtColor(noise_image, noise_gray, cv::COLOR_BGR2GRAY);
harris_result = harris_corner_detect(noise_gray, 2, 0.04, 4e3, 0);
show_harris(noise_image, harris_result, save_dir + "noisy.png");
// Shi-Tomasi 角点检测方法
harris_result = shi_tomasi_corner_detect(origin_gray, 2, 30, best_size);
show_harris(origin_image, harris_result, save_dir + "shi_tomasi.png");
}
int main() {
std::cout << "opencv : " << CV_VERSION << std::endl;
demo_1();
return 0;
}
Cmakelist:
cmake_minimum_required(VERSION 3.12)
project(harris)
# 添加编译选项
SET(CMAKE_CXX_FLAGS "-std=c++14 -O1")
set(PROJECT_SOURCE_DIR ./)
# 生成文件的输出目录
set(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${PROJECT_SOURCE_DIR}/bin)
# include 目录
include_directories(SYSTEM ${CMAKE_CURRENT_SOURCE_DIR}/include/)
# 添加可执行文件
add_executable(harris src/harris_demo.cpp )
# 设置 opencv
set(OpenCV_DIR "D:/environments/C++/OpenCV/opencv-4.5.2/build_no_qt/install")
# 寻找 opencv 的库
find_package(OpenCV REQUIRED)
# 添加 opencv 的头文件
include_directories(${OpenCV_INCLUDE_DIRS})
# 链接 opencv 动态库路径
link_directories(${OpenCV_LIBRARY_DIRS})
# 设置 Eigen3
include_directories(D:/environments/C++/3rdparty/Eigen3/eigen-3.3.9/installed/include/eigen3)
# 设置库文件
target_link_libraries(harris ${OpenCV_LIBS})
使用的图片
十一、参考
1.Harris C, Stephens M. A combined corner and edge detector[C]//Alvey vision conference. 1988, 15(50): 10-5244.
2.Th.Jonarch：[Athos’ CV] Chapter3.局部特征-1: 角点检测和Harris角点检测
3.Harris角点检测原理详解_lwzkiller的专栏-CSDN博客_harris角点检测
4.Shi J. Good features to track[C]//1994 Proceedings of IEEE conference on computer vision and pattern recognition. IEEE, 1994: 593-600.
8881
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
