第三神奇----harris角点的理解和实现（感兴趣区域（roi）变化为角点）

首先，要有这样的意识，梯度的方向会垂直边缘。那么x的梯度和y的梯度相交，也即两个边缘相交就会形成角点。

那么，什么是梯度？在图像处理上，可以这样表达：Gx=I(x+1,y)-I(x,y);Gy=I(x,y+1)-I(x,y).I是图像阵列

那么x，y处的梯度这样表达，▽I（x，y）=|Gx|+|Gy|;梯度的方向θ=arctg(Gx/Gy).

Gx，Gy很像一阶导数，虽然从数学上来讲此种表达有待探讨，但是他并不影响解决问题。

再延伸一下，Gxx=Gx-G(x-1)=I(x+1,y)-I(x,y)-[I(x-1+1,y)-I(x-1,y)]=I(x+1,y)-2*I(x,y)+I(x-1,y)

Gxy=Gyx=I(x+1,y+1)+I(x,y)-I(x,y+1)-I(x+1,y)

Gyy=I(x,y+1)-2*I(x,y)+I(x,y-1)

其次，图像I做高斯平滑，

            for (int j = 1; j < (h - 1); j++)
            {             |  1      2     1 |
                for (int i = 1; i < (w - 1); i++)/     1/16*|   2      4     2 |-------------->高斯模板
                {///              |  1       2     1 |
                    int n0 = (j * w + i);
            tyimg[n0] = (byte)((tongyongimg33gaos[n0 - w - 1] + 2 * tongyongimg33gaos[n0 - w] + tongyongimg33gaos[n0 - w + 1] +
                                            2 * tongyongimg33gaos[n0 - 1] + 4 * tongyongimg33gaos[n0] + 2 * tongyongimg33gaos[n0 + 1] +
                               tongyongimg33gaos[n0 + w - 1] + 2 * tongyongimg33gaos[n0 + w] + tongyongimg33gaos[n0 + w + 1]) / 16);
                }
            }

最后，矩阵A=a11（B=Gxx）    a12（C=Gxy）    ,A的迹=a11+a22=T

                       a21（C=Gxy）   a22（D=Gyy）

echo=（B*D-C*C）-T*T*系数k,所以

            int[,] m_Theta2 = new int[h, w];
            int[,] m_P2 = new int[h, w]; int[,] m_Q2 = new int[h, w];
            int[,] m_Harris = new int[h, w];
            for (int i = 1; i < (h - 1); i++)
            {
                for (int j = 1; j < (w - 1); j++)
                {
                    int n0 = (i * w + j);
                    m_P2[i, j] = tyimg[n0 + 1] + tyimg[n0 - 1] - 2 * tyimg[n0];
                    m_Q2[i, j] = tyimg[n0 + w] + tyimg[n0 - w] - 2 * tyimg[n0];
                    m_Theta2[i, j] = tyimg[n0] - tyimg[n0 + 1] - tyimg[n0 + w] + tyimg[n0 + w + 1];
                    int temp = m_P2[i, j] * m_Q2[i, j];
                    int temp2 = m_P2[i, j] + m_Q2[i, j];//拉普拉斯
                    int tezheng = (temp - m_Theta2[i, j] * m_Theta2[i, j]);
                    int trace = temp2 * temp2;
                    m_Harris[i, j] = (int)(tezheng - m_k * trace);//echo=（B*D-C*C）-T*T*系数k

                    if (Math.Abs(m_Harris[i, j]) > thres && i > 3 && j > 3 && i < h - 3 && j < w - 3)//已经更改201702121652
                    {
                        tyimgshuzu[i * w + j] = 0;
                    }
                    else
                    {
                        tyimgshuzu[i * w + j] = 255;
                    }
                }
            }

那么，角点就被255标记出来了。/睁眼说瞎话，明明是0嘛！

现在得到的角点是整数角点。

我最初的版本很简单，也很直白。当你往前走了很久后，回头再看，你发现harris自相关矩阵中的w相当于高斯平滑，再说的深一点，就是波长函数，当你豁然时，你会发现其中的道理。

另外，人家的自相关矩阵，是一个i*j的窗口，而我的窗口是1*1的。i*j的好处是找亚像素角点时，仍然使用i*j窗口，不用变来变去得了。

还有一个猜测没有验证，找到的角点两个交叉边缘不一定是垂直的。而亚像素角点两个交叉边缘是垂直的，原因是，亚像素角点原理是点积为零（k1*k2=-1,即相互垂直）。

识别特征固然重要，但能度量更重要。度量需要准确，所以会用到亚像素角点，这便是机器视觉中标定的基础。

                                                  刘刚（有光）  ganggangwa   ganggangwawa

 13摄氏度      晴天       201903061012