机器视觉径向畸变校正分步实现（标定，一，harris角点）

高考的那倒计时100天起始，父亲说，把高考的大目标，分成小目标，安排在这100天里，逐步完成。

目标，列提纲，再细化。（高考）

要识别的目标，粗匹配，细匹配。（机器视觉识别）

门，纲，目，科，属，种。（生物分类）

原来都是如此，但请不要忘记，把目标实现，限定在时光中。

我们前面的博客已经完成了标定目标确定和列提纲，现在要做的是再细化。

本节重点，harris角点（确定任意叉（×）整数角点A集合，在A集合中寻找正交（+）亚像素角点）。

需要说的是，harris角点，我们需要的是亚像素角点，第一步求的是任意×的整数角点，第二步，是在整数角点基础上求正交+亚像素角点。图中按钮“全部harris”代码实现：

1，高斯平滑

   for (int j = 1; j < (h - 1); j++)
            {
                for (int i = 1; i < (w - 1); i++)
                {
                    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);
                }
            }

2，任意×的整数角点

     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 特征 = (temp - m_Theta2[i, j] * m_Theta2[i, j]);//*-*
                    int 迹 = temp2 * temp2;

                  //m_k=0.04
                    m_Harris[i, j] = (int)(特征 - m_k * 迹);

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

3,  聚集的角点分组，亚像素角点必在其中

    List<List<Point>> fenzu = new List<List<Point>>();
            染色算法(ref  fenzu, ref tyimgshuzu, W, H);//此算法参考找斑//自己的东西，自己首先要发掘使用它

4，找出分组中，harris响应最大的，亚像素必在其附近

  List<Point> fenzuMaxPt = new List<Point>();
            List<int> fenzumaxVal = new List<int>();

            for (int i = 0; i < fenzu.Count; i++)
            {
                if (fenzu[i].Count < 3) continue;
                List<Point> temp2 = new List<Point>();

                temp2.AddRange(fenzu[i]);

                int max = Math.Abs(m_Harris[temp2[0].Y, temp2[0].X]);
                int maxI = 0, maxJ = 0;
                for (int ii = 0; ii < temp2.Count; ii++)
                {
                    int tempAbsVal = Math.Abs(m_Harris[temp2[ii].Y, temp2[ii].X]);
                    if (tempAbsVal >= max)
                    {
                        max = tempAbsVal;
                        maxI = temp2[ii].X;
                        maxJ = temp2[ii].Y;
                    }
                }
               
                fenzuMaxPt.Add(new Point(maxI, maxJ));//响应的最大值点
                fenzumaxVal.Add(max);//响应的最大值
            }

5，在分组最大响应角点处找亚像素角点

    for (int m = 0; m < fenzuMaxPt.Count; m++)
                {
                    int x = fenzuMaxPt[m].X;
                    int y = fenzuMaxPt[m].Y;
                    int[] xNum = new int[] { 1, 1, 0, -1, -1, -1, 0, 1 };
                    int[] yNum = new int[] { 0, -1, -1, -1, 0, 1, 1, 1 };
                    double  点积和  = 0;//初始化为0不妥，为什么？202005101007

//点积和入门，最大harris响应角点处其八个邻域内的点积和
                   for (int k = 0; k < 8; k++)
                    {
                        int xx = x + xNum[k];
                        int yy = y + yNum[k];
                        int x导数 = tyimg[(yy) * W + xx + 1] - tyimg[(yy) * W + xx];
                        int y导数 = tyimg[(yy) * W + xx + W] - tyimg[(yy) * W + xx];
                        int 梯度值 = Math.Abs(x导数) + Math.Abs(y导数);
                        float 梯度角 = (float)(jiaoduAndxiangxian(x导数, y导数) * 57.3);
                        float 向量X = (float)(梯度值 * Math.Cos(jiaoduAndxiangxian(x导数, y导数)));
                        float 向量Y = (float)(梯度值 * Math.Sin(jiaoduAndxiangxian(x导数, y导数)));
                        float 点积 = 向量X * (x - xx) + 向量Y * (y - yy);
                        点积和 += 点积;
                   }

//最大harris响应角点处，8*8矩形窗口中找亚像素角点（正交角点）

                    PointF tempcorner = tysubcorners1per8(fenzuMaxPt[m].X - 4, fenzuMaxPt[m].Y - 4, 8, 8, W, tyimg);              
                    tyCA_subcornersflitCali.Add(tempcorner);                                  
                }

结果如下：

注：jiaoduAndxiangxian函数在线图像工具设计中出现过；

tysubcorners1per8函数在亚像素角点最小二乘法实现中已经出现。