第三神奇----harris角点之亚像素角点最小二乘法实现（感兴趣区域（roi）变化为角点）

最新推荐文章于 2024-08-21 20:30:00 发布

工业机器视觉设计和实现

最新推荐文章于 2024-08-21 20:30:00 发布

阅读量809

点赞数

分类专栏：机器视觉自我实现（三）

本文链接：https://blog.csdn.net/ganggangwawa/article/details/88232825

版权

机器视觉自我实现（三）专栏收录该内容

59 篇文章 26 订阅

订阅专栏

已经明白两个向量点积为零的道理（x1*x2+y1*y2=0;即k1*k2=-1），亚像素角点即就是整数角点间的插值。

时间回到2016.07.25

今天使用聚集区harris整数角点做为向量P(i),求向量q（亚像素点）。

在角点聚集区域再一次harris角点迭代仍然是整数角点，意义应该不大。

两个向量点积为零是一个非常好听的名字，即让角点聚集区非负点积接近零的插值，不就是最小二乘法吗？

这个想法实在是太伟大了。

两个向量：

(p(i).x-q.x,p(i).y-q.y)和p(i)的梯度向量（x(i),y(i)）,

两个向量的点积为零：x(i)*(p(i).x-q.x)+y(i)*(p(i).y-q.y)=0

点积为零集最小即所有|x(i)*(p(i).x-q.x)+y(i)*(p(i).y-q.y)|和最小------>0

Q(q.x,q.y)= $\sum$ $\binom{n}{i=1}$ [x(i)*(p(i).x-q.x)+y(i)*(p(i).y-q.y)]^2/最小二乘法

令Q'q.x=0=q.x* $\sum$ $\binom{n}{i=1}$ xi*yi+q.y* $\sum$ $\binom{n}{i=1}$ yi^2- $\sum$ $\binom{n}{i=1}$ (xi*yi*pi.x+yi^2*pi.y)

令Q'q.y=0=q.x* $\sum$ $\binom{n}{i=1}$ xi^2+q.y* $\sum$ $\binom{n}{i=1}$ yi*xi- $\sum$ $\binom{n}{i=1}$ (xi^2*pi.y+xi*yi*pi.y)

得解(q.x,q.y)。

另外一种矩阵求解也验证此结果，显然可靠，即

p(i)的梯度向量（x(i),y(i)）；使其组成n*2的列向量A，

$A^{T}$ =| x1 x2 x3 x4 x5...|

| y1 y2 y3 y4 y5... |

$A^{T}$ *A=| $\sum$ $\binom{n}{i=1}$ xi^2 $\sum$ $\binom{n}{i=1}$ xi*yi|

| $\sum$ $\binom{n}{i=1}$ yi*xi $\sum$ $\binom{n}{i=1}$ yi^2|

令误差ε(i)=梯度向量(x(i),y(i))* (q-p(i))

化简得：梯度向量(x(i),y(i))*q= 误差ε(i)+ 梯度向量(x(i),y(i))* p(i)

A=梯度向量(x(i),y(i))，X= $\binom{q.x}{q.y}$ ，b= 梯度向量(x(i),y(i))* p(i)

它实质是矩阵 AX-b（即误差ε(i)=0）最小化问题，矩阵求解得X= $\binom{q.x}{q.y}$ = (A的转置*A)的逆*A的转置*b

(即 $\binom{q.x}{q.y}$ = $(A^{T}A)^{-1}$ * $A^{T}$ *b)。

经推导与上面最小二乘法结果一致。

人有时并不知道自己走的路是否正确，但是两条腿走路是可靠的。

整数harris角点分组后，在每组中求亚像素角点，以下是亚像素角点算法C#实现：

private PointF tysubcorners1per8(int left, int top, int w, int h, int picW, byte[] img)//仅仅梯度幅度不行,M和theta可以计算出梯度向量
{
double D, A, B, E, F, G, H; A = B = D = E = F = G = H = 0;
double _D, _A, _B; _A = _B = _D = 0;
for (int i = 0; i < h; i++)
for (int j = 0; j < w; j++)
{
int xdaoshu = img[(top + i) * picW + left + j + 1] - img[(top + i) * picW + left + j];//无需考虑边界
int ydaoshu = img[(i + top) * picW + j + left + picW] - img[(top + i) * picW + left + j];
double TiDuvalue = Math.Sqrt(xdaoshu * xdaoshu + ydaoshu * ydaoshu);
float MVX = (float)(TiDuvalue * Math.Cos(jiaoduAndxiangxian(xdaoshu, ydaoshu)));//梯度向量X
float MVY = (float)(TiDuvalue * Math.Sin(jiaoduAndxiangxian(xdaoshu, ydaoshu)));//梯度向量Y
_A += (MVX) * (MVX);
_B += (MVX) * (MVY);
_D += (MVY) * (MVY);
E += (MVX) * (MVX) * (left + j);//x
F += (MVX) * (MVY) * (top + i);//y
G += (MVX) * (MVY) * (left + j);//x
H += (MVY) * (MVY) * (top + i);//y
}
double result = _A * _D - _B * _B;
A = _D / (result * 1.0); B = -_B / (result * 1.0); D = _A / (result * 1.0);
double x = A * (E + F) + B * (G + H);
double y = B * (E + F) + D * (G + H);
return new PointF((float)(x), (float)(y));
}