Goal:學習實作Harris corner detector
通過 OpenCV 練習Harris corner detector
了解參數對Harris corner detector的影響
作業環境
Window 11
Visual Studio 2019
OpenCV 4.5.5
CMake 3.24.3
CPU:11th Gen Intel(R) Core(TM) i5-1135G7 @ 2.40GHz 2.42 GHz
RAM:16.0 GB
Harris corner detector
Harris_corner_detector.cpp
#include "opencv2/highgui.hpp"
#include "opencv2/imgproc.hpp"
#include <iostream>
using namespace cv;//簡化了後面程式碼的撰寫,不用寫cv::、std::
using namespace std;
Mat src, src_gray;//src存原圖、src_gray存灰階影像
void cornerHarris_demo(int, void*);
int main(int argc, char** argv)
{
src = imread("C:/Users/User/Desktop/CV_testimage/P4/harris3.jpg");//彩色影像
imshow("Original", src);
if (src.empty())//讀不到影像就退出並顯示錯誤
{
cout << "Could not open or find the image!\n" << endl;
cout << "Usage: " << argv[0] << " <Input image>" << endl;
return -1;
}
cvtColor(src, src_gray, COLOR_BGR2GRAY);//將原始圖像轉成灰階影像
imshow("Graylevel", src_gray);
cornerHarris_demo(0, 0);//角點偵測
waitKey();
return 0;
}
//角點偵測核心功能Harris corner detector
void cornerHarris_demo(int, void*)
{
int thresh = 200;
int blockSize = 2;
int apertureSize = 3;
double k = 0.04;
Mat harris = Mat::zeros(src.size(), CV_32FC1);
cornerHarris(src_gray, harris, blockSize, apertureSize, k);//角點偵測
Mat harris_norm, harris_norm_scaled;
normalize(harris, harris_norm, 0, 255, NORM_MINMAX, CV_32FC1, Mat());//將檢測結果normalize到區間0~255
convertScaleAbs(harris_norm, harris_norm_scaled);//函數轉換,得到正確的值,將圖片轉換成CV_8U的格式
for (int i = 0; i < harris_norm.rows; i++)
{
for (int j = 0; j < harris_norm.cols; j++)
{
if ((int)harris_norm.at<float>(i, j) > thresh)
{
circle(harris_norm_scaled, Point(j, i), 5, Scalar(0), 2, 8, 0);
circle(src, Point(j, i), 5, Scalar(0,0,200), 2, 8, 0);
}
}
}
imshow("Harris", harris_norm_scaled);
imshow("Harris_o", src);
}
src:輸入圖像(8-bit unsigned, 16-bit unsigned...)。
dst:與 src 有相同大小和深度的輸出圖像。
code:顏色空間轉換(ColorConversionCodes)。
dstCn:目標圖像中的通道數(如果參數為 0,通道數自動從 src 和code中導出)
src :輸入單通道 8-bit 或浮點圖像。
dst: 存儲 Harris detector responses的圖像
(類型:CV_32FC1,大小與 src 相同)。
blockSize:Neighborhood size(cornerEigenValsAndVecs)。
ksize :Aperture parameter for the Sobel operator(只能取1、3、5、7)。
k:Harris detector 自由參數。 參見下面的公式。
borderType :Pixel extrapolation method,請參見BorderType(不支持 BORDER_WRAP)。
cornerEigenValsAndVecs():
λ1,λ2 是 M 的未排序特徵值
x1,y1 是對應於 λ1 的特徵向量
x2,y2 是對應於 λ2 的特徵向量
cornerHarris()在圖像上運行 Harris detector。 與 cornerMinEigenVal 和 cornerEigenValsAndVecs 類似,對於每個像素 (x,y),它在 blockSize×blockSize 鄰域上計算一個 2×2 梯度協方差矩陣 M(x,y)。
可以找到圖像中的角點作為此響應圖的局部最大值。
縮放、計算絕對值並將結果轉換為 8-bit。
src:輸入 array.
dst:輸出 array.
alpha:optional scale factor.
beta:optional delta added to the scaled values.
circle() 繪製一個具有給定圓心和半徑的簡單圓或實心圓。
img: 畫圓的圖像。
center: 圓的中心。
radius: 圓的半徑。
color :圓圈的顏色。
thickness: 圓形輪廓的厚度(正值),繪製實心圓(負值)。
lineType: 圓邊界的類型。
shift :中心坐標和半徑值中的小數位數。
Scalar(a,b,c)
a,b,c(藍,綠,紅):0~255
原圖
thresh = 200;blockSize = 2;apertureSize = 3;
circle(src, Point(j, i), 5, Scalar(0,0,200), 2, 8, 0);
thresh = 200;blockSize = 2;apertureSize = 3;
circle(harris_norm_scaled, Point(j, i), 5, Scalar(0), 2, 8, 0);
k = 0.06
k = 0.05
k = 0.04
k = 0.06
k = 0.05
k = 0.04
k 值增加,偵測到的角點增加。
kSize = 3
kSize = 5
kSize = 7
kSize = 3
kSize = 5
kSize = 7
kSize變高,被圈出來的角點變多。
blockSize = 2
blockSize = 3
blockSize = 7
blockSize = 2
blockSize = 3
blockSize = 7
blockSize變高,角點數量減少。
thresh = 100
thresh = 150
thresh = 200
thresh = 250
thresh = 100
thresh = 150
thresh = 200
thresh = 250
thresh變低:較多被認為可能是角點的地方被圈出來
thresh變高:被圈出的角點明顯變少
Bonus 1: cornerEigenValsAndVecs(), cornerMinEigenVal(), and cornerSubPix()
cornerEigenValsAndVecs()
src :輸入單通道 8-bit 或浮點圖像。
dst: 存儲 Harris detector responses的圖像
(類型:CV_32FC1,大小與 src 相同)。
blockSize:Neighborhood size。
ksize :Aperture parameter for the Sobel operator。
borderType :Pixel extrapolation method,請參見BorderType(不支持 BORDER_WRAP)。
計算圖像塊的特徵值和特徵向量以進行corner detection。對每個像素,cornerEigenValsAndVecs() 考慮 blockSize × blockSize 鄰域 S(p)。它將鄰域內導數的協方差矩陣計算為:
(用 Sobel 計算導數)
cornerEigenValsAndVecs.cpp
#include "opencv2/highgui/highgui.hpp"
#include "opencv2/imgproc/imgproc.hpp"
#include <iostream>
#include <stdio.h>
#include <stdlib.h>
using namespace cv;
using namespace std;
Mat src;
Mat srcGray, srcCopy;
Mat dstHarris; //dstHarris用於存儲角點檢測的結果
Mat resHarris; //resHarris用於存儲特徵點選擇後的結果
int HarrisQualityLevel = 50;
int maxQualityLevel = 100;
double HarrisMinVal = 0.0;
double HarrisMaxVal = 0.0;
// 角點檢測函數聲明
void HarrisFunction(int, void*);
int main(int argc, char** argv)
{
src = imread("C:/Users/User/Desktop/CV_testimage/P4/harris3.jpg", 1);
cvtColor(src, srcGray, COLOR_BGR2GRAY);
int blockSize = 3;
int apertureSize = 3;
// 使用cornerEigenValsAndVecs()函數檢測角點
dstHarris = Mat::zeros(srcGray.size(), CV_32FC(6));
resHarris = Mat::zeros(srcGray.size(), CV_32FC1);
cornerEigenValsAndVecs(srcGray, dstHarris, blockSize, apertureSize, BORDER_DEFAULT);
for (int j = 0; j < srcGray.rows; j++)
{
for (int i = 0; i < srcGray.cols; i++)
{
float* lambda = dstHarris.ptr<float>(j, i);
float lambda1 = lambda[0];
float lambda2 = lambda[1];
resHarris.at<float>(j, i) = lambda1 * lambda2 - 0.04 * pow((lambda1 + lambda2), 2);
}
}
minMaxLoc(resHarris, &HarrisMinVal, &HarrisMaxVal, 0, 0, Mat());
HarrisFunction(0, 0);
waitKey(0);
return(0);
}
// 角點檢測函數
void HarrisFunction(int, void*)
{
srcCopy = src.clone();
if (HarrisQualityLevel < 1)
HarrisQualityLevel = 1;
for (int j = 0; j < srcGray.rows; j++)
{
for (int i = 0; i < srcGray.cols; i++)
{
if (resHarris.at<float>(j, i) > HarrisMinVal + (HarrisMaxVal - HarrisMinVal) * HarrisQualityLevel / maxQualityLevel)
circle(srcCopy, Point(i, j), 2, Scalar(0, 0, 255), -1, 8, 0);
}
}
imshow("Harris_cornerEigen", srcCopy);
cout << "Harris Quality Level: " << HarrisQualityLevel << endl;
}
blockSize = 3;
kSize = 3;
blockSize = 3;
kSize = 5;
blockSize = 3;
kSize = 7;
跟Harris corner detector相同,kSize變大,被圈出來的角點變多。
cornerMinEigenVal()
src :輸入單通道 8-bit 或浮點圖像。
dst: 存儲 Harris detector responses的圖像
(類型:CV_32FC1,大小與 src 相同)。
blockSize:Neighborhood size(cornerEigenValsAndVecs)。
ksize :Aperture parameter for the Sobel operator。
borderType :Pixel extrapolation method,請參見BorderType(不支持 BORDER_WRAP)。
計算用於角點檢測的梯度矩陣的最小特徵值,該函數與cornerEigenValsAndVecs類似,但根據cornerEigenValsAndVecs描述中的公式計算並存儲導數協方差矩陣的最小特徵值min(λ1,λ2)。
cornerMinEigenVal.cpp
#include "opencv2/highgui/highgui.hpp"
#include "opencv2/imgproc/imgproc.hpp"
#include <iostream>
#include <stdio.h>
#include <stdlib.h>
using namespace cv;
using namespace std;
Mat src;
Mat srcGray, srcCopy; //srcCopy用於繪圖,srcGray用於檢測角點
Mat dstShiTomasi; //dstShiTomasi用於存儲角點檢測的結果
int ShiTomasiQualityLevel = 50;
int maxQualityLevel = 100;
double ShiTomasiMinVal;
double ShiTomasiMaxVal;
void ShiTomasiFunction(int, void*);
int main(int argc, char** argv)
{
src = imread("C:/Users/User/Desktop/CV_testimage/P4/harris3.jpg", 1);
cvtColor(src, srcGray, COLOR_BGR2GRAY);
int blockSize = 3;
int apertureSize = 7;
// 使用cornerMinEigenVal()函數檢測角點
dstShiTomasi = Mat::zeros(srcGray.size(), CV_32FC1);
cornerMinEigenVal(srcGray, dstShiTomasi, blockSize, apertureSize, BORDER_DEFAULT);
minMaxLoc(dstShiTomasi, &ShiTomasiMinVal, &ShiTomasiMaxVal, 0, 0, Mat());
ShiTomasiFunction(0, 0);
waitKey(0);
return(0);
}
void ShiTomasiFunction(int, void*)
{
srcCopy = src.clone();
if (ShiTomasiQualityLevel < 1)
ShiTomasiQualityLevel = 1;
for (int j = 0; j < srcGray.rows; j++)
{
for (int i = 0; i < srcGray.cols; i++)
{
if (dstShiTomasi.at<float>(j, i) > ShiTomasiMinVal + (ShiTomasiMaxVal - ShiTomasiMinVal) * ShiTomasiQualityLevel / maxQualityLevel)
circle(srcCopy, Point(i, j), 5, Scalar(0, 0, 255), 2, 8, 0);
}
}
imshow("corner_min", srcCopy);
}
kSize = 3;
kSize = 5;
kSize = 7;
跟Harris corner detector不同,kSize變大,被圈出來的角點變少。
cornerSubPix()
細化corner位置。
image:輸入單通道、8-bit或float圖像。
corners:輸入角的初始坐標和為輸出提供完善的坐標。
winSize:search windows邊長的一半。 例如,如果 winSize=Size(5,5) ,則使用 (5*2+1)×(5*2+1)=11×11 search windows。
zeroZone:搜索區中間死區大小的一半,在該區域上未進行以下公式中的求和。 它有時用於避免自相關矩陣可能出現的奇點。 (-1,-1) 的值表示沒有這樣的大小。
criteria:角點細化迭代過程的終止標準。 也就是說,角位置細化過程在 criteria.maxCount 迭代之後或在某個迭代中角位置移動小於 criteria.epsilon 時停止。
該函數迭代以找到corner或radial saddle points 的Sub-pixel精確位置,如右圖。
Sub-pixel corner detector基於以下觀察:
從中心 q 到位於 q 鄰域內的點 p 的每個向量都與受圖像和測量雜訊影響的 p 處的圖像梯度正交。 表示為:
其中 DIpi 是 q 鄰域中點 pi 之一的圖像梯度。要找到 q 的值,以使 εi 最小化。 可以建立一個方程組,其中 εi 設置為零:
其中梯度在 q 的鄰域(search window)內求和。 用第一個梯度項 G 和第二個梯度項 b 給出:
該算法將鄰域窗口的中心設置在這個新的中心 q 處,然後迭代直到該中心保持在設定的threshold內。
cornerSubPix.cpp
#include "opencv2/highgui.hpp"
#include "opencv2/imgproc.hpp"
#include <iostream>
using namespace cv;
using namespace std;
Mat src, src_gray;
int maxCorners = 10;
int maxTrackbar = 25;
RNG rng(12345);
const char* source_window = "image";
void goodFeaturesToTrack_Demo(int, void*);
int main(int argc, char** argv)
{
src = imread("C:/Users/User/Desktop/CV_testimage/P4/harris3.jpg");
if (src.empty())
{
cout << "Could not open or find the image!\n" << endl;
cout << "Usage: " << argv[0] << " <Input image>" << endl;
return -1;
}
cvtColor(src, src_gray, COLOR_BGR2GRAY);
imshow(source_window, src);
goodFeaturesToTrack_Demo(0, 0);
waitKey();
return 0;
}
void goodFeaturesToTrack_Demo(int, void*)
{
maxCorners = MAX(maxCorners, 1);
vector<Point2f> corners;//創建一個 Point2f 對象的向量 corners,保存檢測到的角點
double qualityLevel = 0.01;
double minDistance = 10;//角點之間的最小距離
int blockSize = 3, gradientSize = 3;
bool useHarrisDetector = false;//是否使用 Harris corner detector
double k = 0.04;
Mat copy = src.clone();//clone()創建原始圖像的副本
goodFeaturesToTrack(src_gray,
corners,
maxCorners,
qualityLevel,
minDistance,
Mat(),
blockSize,
gradientSize,
useHarrisDetector,
k);
cout << "** Number of corners detected: " << corners.size() << endl;
int radius = 4;
for (size_t i = 0; i < corners.size(); i++)
{
circle(copy, corners[i], radius, Scalar(rng.uniform(0, 255), rng.uniform(0, 256), rng.uniform(0, 256)), FILLED);
}
namedWindow(source_window);
imshow(source_window, copy);
Size winSize = Size(5, 5);
Size zeroZone = Size(-1, -1);
TermCriteria criteria = TermCriteria(TermCriteria::EPS + TermCriteria::COUNT, 40, 0.001);
cornerSubPix(src_gray, corners, winSize, zeroZone, criteria);//cornerSubPix()精確調整角點的位置
for (size_t i = 0; i < corners.size(); i++)
{
cout << " -- Refined Corner [" << i << "] (" << corners[i].x << "," << corners[i].y << ")" << endl;
circle(copy, corners[i], radius, Scalar(0,0,255), FILLED);
}
imshow("cornerSubPix", copy);
}
經過cornerSubPix()精確調整角點的位置 ,會更精確
Bonus 2 : Matlab Version of Harris Detector
harris_corner.m
function [ x, y, scores, Ix, Iy ] = harris_corner( image )
% with a high degree of 'cornerness' from
%RGB image matrix of type uint8
% Input - image = NxMx3 RGB image matrix
% Output - x = nx1 vector denoting the x location of each of n
% detected keypoints
% y = nx1 vector denoting the y location of each of n
% detected keypoints
% scores = an nx1 vector that contains the value (R) to which a
% a threshold was applied, for each keypoint
% Ix = A matrix with the same number of rows and columns as the
% input image, storing the gradients in the x-direction at each
% pixel
% Iy = A matrix with the same nuimber of rwos and columns as the
% input image, storing the gradients in the y-direction at each
% pixel
% compute the gradients, re-use code from HW2P, use window size of 5px
% convert image to grayscale first
image = imread('C:\Users\User\Desktop\CV_testimage\P4\harris3.jpg');
G = rgb2gray(image);
% convert to double
G2 = im2double(G);
% create X and Y Sobel filters
horizontal_filter = [1 0 -1; 2 0 -2; 1 0 -1];
vertical_filter = [1 2 1; 0 0 0 ; -1 -2 -1];
% using imfilter to get our gradient in each direction
filtered_x = imfilter(G2, horizontal_filter);
filtered_y = imfilter(G2, vertical_filter);
% store the values in our output variables, for clarity
Ix = filtered_x;
Iy = filtered_y;
% Compute the values we need for the matrix...
% Using a gaussian blur, because I get more positive values after applying
% it, my values all skew negative for some reason...
f = fspecial('gaussian');
Ix2 = imfilter(Ix.^2, f);
Iy2 = imfilter(Iy.^2, f);
Ixy = imfilter(Ix.*Iy, f);
% set empirical constant between 0.04-0.06
k = 0.04;
num_rows = size(image,1);
num_cols = size(image,2);
% create a matrix to hold the Harris values
H = zeros(num_rows, num_cols);
% % get our matrix M for each pixel
for y = 6:size(image,1)-6
for x = 6:size(image,2)-6
% build a matrix using what we computed earlier
% M = zeros(2,2);
% M(1,1) = Ix2(y,x);
% M(2,1) = Ixy(y,x);
% M(1,2) = Ixy(y,x);
% M(2,2) = Iy2(y,x);
%
% calculate means
% Ix2 mean
Ix2_matrix = Ix2(y-2:y+2,x-2:x+2);
Ix2_mean = mean(Ix2_matrix(:));
% Iy2 mean
Iy2_matrix = Iy2(y-2:y+2,x-2:x+2);
Iy2_mean = mean(Iy2_matrix(:));
% Ixy mean
Ixy_matrix = Ixy(y-2:y+2,x-2:x+2);
Ixy_mean = mean(Ixy_matrix(:));
Matrix = [Ix2_mean, Ixy_mean;
Ixy_mean, Iy2_mean];
R1 = det(Matrix) - (k * trace(Matrix)^2);
% compute R, using te matrix we just created
%R = det(M) - ( k * trace(M)^2 );
%alt
%R = (Ix2.*Iy2 - Ixy.^2) - k(Ix2+Iy2).^2;
% store the R values in our Harris Matrix
% H(y,x) = R;
H(y,x) = R1;
end
end
% set threshold of 'cornerness' to 5 times average R score
avg_r = mean(mean(H))
threshold = abs(5 * avg_r)
[row, col] = find(H > threshold);
scores = [];
%get all the values
for index = 1:size(row,1)
%see what the values are
r = row(index);
c = col(index);
%store the scores
%score(index) = H(r,c);
scores = cat(2, scores,H(r,c));
end
y = row;
x = col;
% This needs to be LAST
% non max suppression
% http://www.mathworks.com/help/images/ref/imregionalmax.html
% BW = imregionalmax(H);
figure; imshow(image)
hold on
for i = 1:size(scores,2)
plot(x(i), y(i), 'ro', 'MarkerSize', scores(i) * 2); % you may need to play with this multiplier or divisor based on your image
% I've used –> (/1000) to (* 10)
end
end
Input
im:待處理的圖像(灰階圖像)。
sigma:在計算結構張量時使用的平滑高斯標準偏差。
k:Harris測量中與結構張量的跡和行列式相關的參數(0.04-0.06)。
varargin:可選的半徑、thresh、N、subpixel 和display 。
Output
cim:角強度圖像。
r:角點的行坐標。
c:角點的列坐標。
轉換成灰階圖像,並轉換成double
建立X方向和Y方向的Sobel filters
計算derivatives 和結構張量的elements,用高斯去雜訊
計算Harris corner detection response
輸出結果
Reference
OpenCV: Feature Detection,OpenCV,Dec 28,2022,https://docs.opencv.org/4.7.0/dd/d1a/group__imgproc__feature.html#ga354e0d7c86d0d9da75de9b9701a9a87e
OpenCV: Detecting corners location in subpixels,OpenCV,Dec 25 2021,https://docs.opencv.org/4.5.5/dd/d92/tutorial_corner_subpixels.html
OpenCV - Detecting corner location in subpixel(cornerSubPix()) using C++,Another Techs,https://anothertechs.com/programming/cpp/opencv/corernsubpixel/
OpenCV中feature2D学习——自定义角点检测函数,holybin,Nov 10,2014,https://blog.csdn.net/holybin/article/details/40984955
Harris_Corner_Detector-Matlab,GitHub,adpoe, Nov 27, 2016,https://github.com/adpoe/Harris_Corner_Detector-Matlab
Harris Corner detection - File Exchange - MATLAB Central,MathWorks,27 Feb 2014,https://www.mathworks.com/matlabcentral/fileexchange/45704-harris-corner-detection