OpenCV_Mark

最新推荐文章于 2024-03-29 14:07:45 发布

Domo1036

最新推荐文章于 2024-03-29 14:07:45 发布

阅读量237

点赞数

原文链接：https://blog.csdn.net/u014020344/article/details/50196635?utm_source=blogxgwz16

版权

Sobel 实现C++ 备忘

#include "cv.h" 
#include "highgui.h" 
#include "math.h" 
#include "stdio.h" 
#include "malloc.h"
IplImage *image; //声明IplImage指针 
int height, width; 
CvScalar s;
int sobel_y[9] = { 1, 2, 1, 0, 0, 0, -1, -2, -1 };  //y方向sobel算子 
int sobel_x[9] = { 1, 0, -1, 2, 0, -2, 1, 0, -1 };  //x方向sobel算子 
void sobel() 
{ 
    int i, j, k; 
    int grayx = 0, grayy = 0, gray; 
    int *data; 
    int a[9]; 
    data = (int *)malloc(height*width*sizeof(int)); 
    for (i = 0; i<height; i++) 
    { 
        for (j = 0; j<width; j++) 
        { 
            s = cvGet2D(image, i, j); 
            gray = (int)s.val[0]; 
            data[i*width + j] = gray; 
        } 
    }
    for (i = 1; i<height - 1; i++) 
    { 
        for (j = 1; j<width - 1; j++) 
        { 
            grayx = 0; 
            grayy = 0; 
            s = cvGet2D(image, i, j); 
            a[0] = data[width*(i - 1) + j - 1]; 
            a[1] = data[width*(i - 1) + j]; 
            a[2] = data[width*(i - 1) + j + 1];
            a[3] = data[width*i + j - 1]; 
            a[4] = data[width*i + j]; 
            a[5] = data[width*i + j + 1];
            a[6] = data[width*(i + 1) + j - 1]; 
            a[7] = data[width*(i + 1) + j]; 
            a[8] = data[width*(i + 1) + j + 1]; 
            for (k = 0; k<9; k++) 
                grayy += a[k] * sobel_y[k]; 
            for (k = 0; k<9; k++) 
                grayx += a[k] * sobel_x[k];
            s.val[0] = (abs(grayx) + abs(grayy));//此处需要取绝对值 
            cvSet2D(image, i, j, s); 
        } 
    } 
    free(data);
}
int main(int argc, char** argv) 
{
    image = cvLoadImage("image1.jpg", 0); 
    IplImage* image1 = cvLoadImage("image1.jpg", 1); 
    height = image->height; 
    width = image->width;
    sobel();
    cvNamedWindow("Image", 1);//创建窗口 
    cvNamedWindow("Sobel", 1);//创建窗口 
    cvShowImage("Image", image1);//显示图像 
    cvShowImage("Sobel", image);//显示图像 
    cvWaitKey(0); //等待按键
    cvDestroyWindow("Image");//销毁窗口 
    cvReleaseImage(&image); //释放图像 
    return 0;
}

由方向校正Sobel系数 Mark

	//生成核//根据向量方向，修改sobel算子
	std::array<cv::Mat, 2> colorWish::genKernel( cv::Vec4f  &axis, int kerRadius, float order,bool isRefine)
	{
		std::array<cv::Mat, 2>  mKs;
		cv::Mat   mK1, mK2, A0, Ar, ArVis;
		const double pi = 3.141592653;
 
		//计算偏离角度-使用第一个方向//注意角度，一定限制在180度之内
		double alpha = arccos(0, 0, axis[0], axis[1]);
		alpha = alpha > pi ? 2 * pi - alpha : alpha;
		//alpha = 0;
		double alphaVis = arccos(0, 0, axis[2], axis[3]);
		alphaVis = alphaVis > pi ? 2 * pi - alphaVis : alphaVis;
 
		//测试-标准核-0角度
		A0 = cv::Mat::zeros(kerRadius, kerRadius, CV_32FC1);
		const float r = kerRadius / 2;//直接取整
		const int c = r;
		A0.at<float>(c, c) = 0;//只取第一象限，四象限根据对称规则来取
		for (int i = 0; i < r; ++i) {
			for (int j = 0; j < r; ++j) {
				A0.at<float>(c + 1 + j, c + 1 + i) = (j + 1)*(c - i) / r;//使用线性核
			}
		}
 
		//旋转//根据角度进行判断，不能使用一个公式//注意下标
		Ar = cv::Mat::zeros(kerRadius, kerRadius, CV_32FC1);
		ArVis = cv::Mat::zeros(kerRadius, kerRadius, CV_32FC1);
 
		// 注意坐标旋转方向--象限的对应关系
		double x, y;
		for (int i = 0; i < kerRadius; ++i) {
			y = i - c;
			double d = sqrt(y *y + y*y);// *sqrt(2) / 2;//平方核
			//d = abs(y) * 1;//线性核//修改距离判断，使符合线性叠加！
			//d = r;//不修改线性，
			for (int j = 0; j < kerRadius; ++j) {
				x = j - c;
				d = sqrt(x *x + y*y);// 
				d = pow(d, order);
				//d*=sqrt(2) / 2;
				double Beta = arccos(c, c, i, j);//已约束y值为正方向//需要修改顺序
				double AReal = Beta + alpha;
				double ARealVis = Beta + alphaVis;;// Beta - alphaVis;
				 //使用线性核//y= sin ,x = cos;
				if (AReal > 0 && AReal < pi / 2){
					Ar.at< float >(j, i) = (d*cos(AReal) + 0)*(r + 1 - d*sin(AReal)) / r;//同步修改公式
				}
				else {
					if (AReal > pi / 2 && AReal < pi){
						AReal = pi - AReal;
						//std::cout << d*cos(AReal) << "   "; std::cout << d*sin(AReal) << "   ";
						Ar.at< float >(j, i) = 0 - (d*cos(AReal) + 0)*(r + 1 - d*sin(AReal)) / r;// 注意坐标旋转方向--象限的对应关系
					}else {
						if (AReal > pi  && AReal < 3 * pi / 2){
							AReal = AReal - pi;
							//std::cout << d*cos(AReal) << "   "; std::cout << d*sin(AReal) << "   ";
							Ar.at< float >(j, i) = 0 - (d*cos(AReal) + 0)*(r + 1 - d*sin(AReal)) / r;
						}else {
							if (AReal > 2 * pi){
								AReal = AReal - pi * 2;
								//std::cout << d*cos(AReal) << "   "; std::cout << d*sin(AReal) << "   ";
								Ar.at< float >(j, i) = (d*cos(AReal) + 0)*(r + 1 - d*sin(AReal)) / r;
							}else {
								AReal = pi * 2 - AReal;
								//std::cout << d*cos(AReal) << "   "; std::cout << d*sin(AReal) << "   ";
								Ar.at< float >(j, i) = (d*cos(AReal) + 0)*(r + 1 - d*sin(AReal)) / r;
							}
						}//if ( AReal > pi  && AReal < 3*pi/2 )
					}//if (AReal > pi / 2 && AReal < pi)
				}//if (AReal > pi / 2 && AReal < pi )
				 //std::cout << Ar.at<float>(j, i) << "   ";
				if (ARealVis > 0 && ARealVis < pi / 2){
					ArVis.at< float >(j, i) = d*cos(ARealVis)*(r + 1 - d*sin(ARealVis)) / r;
				}else {
					if (ARealVis > pi / 2 && ARealVis < pi){
						ARealVis = pi - ARealVis;
						ArVis.at< float >(j, i) = 0 - d*cos(ARealVis)*(r + 1 - d*sin(ARealVis)) / r;
					}else {
						if (ARealVis > pi  && ARealVis < 3 * pi / 2){
							ARealVis = ARealVis - pi;
							ArVis.at< float >(j, i) = 0 - d*cos(ARealVis)*(r + 1 - d*sin(ARealVis)) / r;
						}else {
							if (ARealVis > 2 * pi)
							{
								ARealVis = ARealVis - pi * 2;
								ArVis.at< float >(j, i) = d*cos(ARealVis)*(r + 1 - d*sin(ARealVis)) / r;//(3*pi/2,pi)
							}
							else {
								ARealVis = pi * 2 - ARealVis;
								ArVis.at< float >(j, i) = d*cos(ARealVis)*(r + 1 - d*sin(ARealVis)) / r;//(3*pi/2,pi)
							}
						}//if ( AReal > pi  && AReal < 3*pi/2 )
					}//if (AReal > pi / 2 && AReal < pi)
				}//if (AReal > pi / 2 && AReal < pi )
				 //std::cout << Ar.at<float>(j, i) << "   ";
			}//for
		}//for
 
		Ar.at<float>(c, c) = 0;
		ArVis.at<float>(c, c) = 0;
 
		//转换坐标关系
		cv::transpose(Ar, Ar);
		cv::transpose(ArVis, ArVis);
 
		for (int i = 0; i < kerRadius; ++i) {
			std::cout << std::endl;
			for (int j = 0; j < kerRadius; ++j) {
				//std::swap(Ar.at<float>(j, i), Ar.at<float>(i, j));
				std::cout << Ar.at<float>(j, i) << "  ";
			}
		}
		std::cout << std::endl;
		for (int i = 0; i < kerRadius; ++i) {
			std::cout << std::endl;
			for (int j = 0; j < kerRadius; ++j) {
				std::cout << ArVis.at<float>(j, i) << "  ";
			}
		}
		std::cout << std::endl;
 
		mKs[0] = Ar.clone();
		mKs[1] = ArVis.clone();
		return mKs;
	}//genKernel(cv::Vec4f  &axis, int kerRadius)

拟合圆

最小二乘拟合圆曲线方程为：R²=（X-A）²+（Y-B）²

       /// <summary>
       /// 拟合圆
       /// </summary>
        private struct Circle
        {
            public double X;//圆心X
            public double Y;//圆心Y
            public double R;//半径R
        }
 
         /// <summary>
        /// 定义的拟合点
        /// </summary>
        private struct PixelPoint
        {
            public double x;
            public double y;
        }
        /// <summary>
        /// 拟合圆程序
        /// </summary>
        /// <param name="pPointList">要拟合点集</param>
        /// <returns>返回圆对象</returns>
        private Circle FittingCircle(List<PixelPoint> pPointList)
        {
            Circle pCircle=new Circle();
            if (pPointList.Count < 3)
            {
                XtraMessageBox.Show("最少需要三个点进行拟合");
                return pCircle;
            }
            double X1 = 0;
            double Y1 = 0;
            double X2 = 0;
            double Y2 = 0;
            double X3 = 0;
            double Y3 = 0;
            double X1Y1 = 0;
            double X1Y2 = 0;
            double X2Y1 = 0;
            for (int i = 0; i < pPointList.Count; i++)
            {
                X1 = X1 + pPointList[i].x;
                Y1 = Y1 + pPointList[i].y;
                X2 = X2 + pPointList[i].x * pPointList[i].x;
                Y2 = Y2 + pPointList[i].y * pPointList[i].y;
                X3 = X3 + pPointList[i].x * pPointList[i].x * pPointList[i].x;
                Y3 = Y3 + pPointList[i].y * pPointList[i].y * pPointList[i].y;
                X1Y1 = X1Y1 + pPointList[i].x * pPointList[i].y;
                X1Y2 = X1Y2 + pPointList[i].x * pPointList[i].y * pPointList[i].y;
                X2Y1 = X2Y1 + pPointList[i].x * pPointList[i].x * pPointList[i].y;
            }
            double C, D, E, G, H, N;
            double a, b, c;
            N = pPointList.Count;
            C = N * X2 - X1 * X1;
            D = N * X1Y1 - X1 * Y1;
            E = N * X3 + N * X1Y2 - (X2 + Y2) * X1;
            G = N * Y2 - Y1 * Y1;
            H = N * X2Y1 + N * Y3 - (X2 + Y2) * Y1;
            a = (H * D - E * G) / (C * G - D * D);
            b = (H * C - E * D) / (D * D - G * C);
            c = -(a * X1 + b * Y1 + X2 + Y2) / N;
            pCircle.X = a / (-2);
            pCircle.Y = b / (-2);
            pCircle.R = Math.Sqrt(a * a + b * b - 4 * c) / 2;
            return pCircle;
        }

多个点拟合圆并给出圆心坐标

public static PointF FitCenter(List<PointF> pts, double epsilon = 0.1)
        {
            double totalX = 0, totalY = 0;
            int setCount = 0;

            for (int i = 0; i < pts.Count; i++)
            {
                for (int j = 1; j < pts.Count; j++)
                {
                    for (int k = 2; k < pts.Count; k++)
                    {
                        double delta = (pts[k].X - pts[j].X) * (pts[j].Y - pts[i].Y) - (pts[j].X - pts[i].X) * (pts[k].Y - pts[j].Y);

                        if (Math.Abs(delta) > epsilon)
                        {
                            double ii = Math.Pow(pts[i].X, 2) + Math.Pow(pts[i].Y, 2);
                            double jj = Math.Pow(pts[j].X, 2) + Math.Pow(pts[j].Y, 2);
                            double kk = Math.Pow(pts[k].X, 2) + Math.Pow(pts[k].Y, 2);

                            double cx = ((pts[k].Y - pts[j].Y) * ii + (pts[i].Y - pts[k].Y) * jj + (pts[j].Y - pts[i].Y) * kk) / (2 * delta);
                            double cy = -((pts[k].X - pts[j].X) * ii + (pts[i].X - pts[k].X) * jj + (pts[j].X - pts[i].X) * kk) / (2 * delta);

                            totalX += cx;
                            totalY += cy;

                            setCount++;
                        }
                    }
                }
            }

            if (setCount == 0)
            {
                //failed
                return PointF.Empty;
            }

            return new PointF((float)totalX / setCount, (float)totalY / setCount);
        }

.net core(c#)拟合圆测试
说明
很多时候，我们需要运动物体的转弯半径去描述其机器性能。但在大多数的现实条件下，我们只能够获取到运动物体的 GPS 位置点集，并不能直接得到转弯半径或者圆心位置。为此，我们可以利用拟合圆的方式得到圆坐标方程，由此得到转弯半径和圆心位置。

解决过程
关于拟合圆方程的方法有很多，曾经在这篇译文中获益良多代数逼近法、最小二乘法、正交距离回归法来拟合圆及其结果对比（Python）。此系列文中也给出了提及的三种方法的性能及效果对比，最终得出最优的解决方案就是最小二乘法。由于最近的学习中又进一步了解到，可以利用线性代数的方法去求解。本着大学课程中曾学过的《线性代数》知识，所以想着用此方法再加以解决该问题，以作最对比。

接下来，本文就最小二乘法和线性代数的方法求取圆方程作一论述。

准备
引用矩阵计算库MathNet.Numerics。该库是一个强大的科学计算库，遵循 .Net Standard，所以可以跨平台使用。

创建描述圆的类

public class Circle
{
    /// <summary>
    /// 圆心横坐标
    /// </summary>
    /// <value></value>
    public   double X { get; set; }
    /// <summary>
    /// 圆心纵坐标
    /// </summary>
    /// <value></value>
    public   double Y { get; set; }
    /// <summary>
    /// 圆半径
    /// </summary>
    /// <value></value>
    public  double R { get; set; }
}

画图，引用System.Drawing.Common库，以实现跨平台的图像生成。接下来，我们简单的实现一个图像帮助类来进行图像绘制

public class ImageHelp
{
    private Image _image;
    public ImageHelp(int width, int height)
    {
        _image = new Bitmap(width, height);
        var graph = Graphics.FromImage(_image);
        graph.Clear(Color.White);
    }
    public void DrawCicle(Circle circle, Brush brush)
    {
        var graph = Graphics.FromImage(_image);
        var count=200;
        var fitPoints = new Point[count+1];
        var step = 2 * Math.PI / count;
        for (int i = 0; i < count; i++)
        {
            //circle
            var p = new Point();
            p.X = (int)(circle.X + Math.Cos(i * step) * circle.R);
            p.Y = (int)(circle.Y + Math.Sin(i * step) * circle.R);
            fitPoints[i] = p;
        }
        fitPoints[count] = fitPoints[0];//闭合圆
        graph.DrawLines(new Pen(brush, 2), fitPoints);
        graph.Dispose();
    }
    public void DrawPoints(double[] X, double[] Y, Brush brush)
    {
        var graph = Graphics.FromImage(_image);
        for (int i = 0; i < X.Length; i++)
        {
            graph.DrawEllipse(new Pen(brush, 2), (int)X[i], (int)Y[i], 6, 6);
        }
        graph.Dispose();
    }
    public void SaveImage(string file)
    {
        _image.Save(file, System.Drawing.Imaging.ImageFormat.Png);
    }
}

模拟点集，由于现实中的数据采集存在着精度、数据记录等众多不确定因素的影像。模拟点集中也将加入一定程度的噪音。以下代码中 x 与 y 中存储着我们的点集数据：

var count = 50;
var step = 2 * Math.PI / 100;
var rd = new Random();
//参照圆
var x0 = 204.1;
var y0 = 213.1;
var r0 = 98.4;
//噪音绝对差
var diff = (int)(r0 * 0.1);
var x = new double[count];
var y = new double[count];
//输出点集
for (int i = 0; i < count; i++)
{
    //circle
    x[i] = x0 + Math.Cos(i * step) * r0;
    y[i] = y0 + Math.Sin(i * step) * r0;
    //noise
    x[i] += Math.Cos(rd.Next() % 2 * Math.PI) * rd.Next(diff);
    y[i] += Math.Cos(rd.Next() % 2 * Math.PI) * rd.Next(diff);
}

最小二乘法
网上有很多的原理解析，上文中提到的译文中也有提及，这里不在过多赘述。直接贴出 c#代码实现：

public Circle LeastSquaresFit(double[] X, double[] Y)
{
    if (X.Length < 3)
    {
        return null;
    }
    double cent_x = 0.0,
        cent_y = 0.0,
        radius = 0.0;
    double sum_x = 0.0f, sum_y = 0.0f;
    double sum_x2 = 0.0f, sum_y2 = 0.0f;
    double sum_x3 = 0.0f, sum_y3 = 0.0f;
    double sum_xy = 0.0f, sum_x1y2 = 0.0f, sum_x2y1 = 0.0f;
    int N = X.Length;
    double x, y, x2, y2;
    for (int i = 0; i < N; i++)
    {
        x = X[i];
        y = Y[i];
        x2 = x * x;
        y2 = y * y;
        sum_x += x;
        sum_y += y;
        sum_x2 += x2;
        sum_y2 += y2;
        sum_x3 += x2 * x;
        sum_y3 += y2 * y;
        sum_xy += x * y;
        sum_x1y2 += x * y2;
        sum_x2y1 += x2 * y;
    }
    double C, D, E, G, H;
    double a, b, c;
    C = N * sum_x2 - sum_x * sum_x;
    D = N * sum_xy - sum_x * sum_y;
    E = N * sum_x3 + N * sum_x1y2 - (sum_x2 + sum_y2) * sum_x;
    G = N * sum_y2 - sum_y * sum_y;
    H = N * sum_x2y1 + N * sum_y3 - (sum_x2 + sum_y2) * sum_y;
    a = (H * D - E * G) / (C * G - D * D);
    b = (H * C - E * D) / (D * D - G * C);
    c = -(a * sum_x + b * sum_y + sum_x2 + sum_y2) / N;
    cent_x = a / (-2);
    cent_y = b / (-2);
    radius = Math.Sqrt(a * a + b * b - 4 * c) / 2;
    var result = new Circle();
    result.X = cent_x;
    result.Y = cent_y;
    result.R = radius;
    return result;
}

线性代数
从标准圆方程(x-c1)^2+(y-c2)2=r^{2中进行方程变换得到2xc1+2yc2+(r}2−c1^2−c22)=x^2+y2，其中，我们c3替换常量值r^2−c12−c2^2，即：r2−c1^2−c22=c3。由此，我们得到2xc1+2yc2+c3=x^2+y2，将点集带入，方程就只剩三个未知数`c1,c2 和 c3。

简单起见，假设我们有四个点{[0,5],[0,-5],[5,0],[-5,0]},代入方程可得到四个方程：
0c1 + 10c2 + c3 = 25
0c1 - 10c2 + c3 = 25
10c1 + 0c2 + c3 = 25
-10c1 + 0c2 + c3 = 25
该方程组比较简单，一眼便能看出解。但用线性代数我们可以得到矩阵：

/*****************AB****C/
| 0c1 10c2 1c3| | 0 10 1| |c1| |25|
| 0c1 -10c2 1c3| = | 0 -10 1| * |c2| = |25|
| 10c1 0c2 1c3| | 10 0 1| |c3| |25|
|-10c1 0c2 1c3| |-10 0 1| |25|
在矩阵方程中AB=C，只需求出矩阵B即可得到方程组的解。c#中MathNet.Numerics可以轻松胜任这一工作：

public Circle LinearAlgebraFit(double[] X, double[] Y)
{
    if (X.Length < 3)
    {
        return null;
    }
    var count = X.Length;
    var a = new double[count, 3];
    var c = new double[count, 1];
    for (int i = 0; i < count; i++)
    {
        //matrix
        a[i, 0] = 2 * X[i];
        a[i, 1] = 2 * Y[i];
        a[i, 2] = 1;
        c[i, 0] = X[i] * X[i] + Y[i] * Y[i];
    }
    var A = DenseMatrix.OfArray(a);
    var C = DenseMatrix.OfArray(c);
    //A*B=C
    var B = A.Solve(C);
    double c1 = B.At(0, 0),
        c2 = B.At(1, 0),
        r = Math.Sqrt(B.At(2, 0) + c1 * c1 + c2 * c2);
    var result = new Circle();
    result.X = c1;
    result.Y = c2;
    result.R = r;
    return result;
}

总结

Console.WriteLine($"raw   c1:{x0}, c2:{y0}, r:{r0}");
var fit = new FitCircle();
var sth = new Stopwatch();
sth.Start();
var lsf = fit.LeastSquaresFit(x, y);![](https://img2018.cnblogs.com/blog/1214143/201908/1214143-20190804173821455-1022769486.jpg)
 
 
Console.WriteLine($"LeastSquaresFit   c1:{lsf.X}, c2:{lsf.Y}, r:{lsf.R}, time:{sth.Elapsed}");
sth.Restart();
var laf = fit.LinearAlgebraFit(x, y);
Console.WriteLine($"LinearAlgebraFit  c1:{laf.X}, c2:{laf.Y}, r:{laf.R}, time:{sth.Elapsed}");
var img = new ImageHelp(512, 512);
img.DrawPoints(x, y, Brushes.Red);
img.DrawCicle(lsf, Brushes.Green);
img.DrawCicle(laf, Brushes.Orange);
img.SaveImage("graph.jpeg");

控制台输出：

raw   c1:204.1, c2:213.1, r:98.4
LeastSquaresFit   c1:204.791071061878, c2:210.86075318831, r:100.436594821545, time:00:00:00.0011029
LinearAlgebraFit  c1:204.791071061878, c2:210.860753188315, r:100.436594821541, time:00:00:00.1691119

从结果中可以看出，两种方法的结果基本一样，在小数点后好几位才出现差别。但是其计算效率却差异巨大，最小二乘法比线性代数快上 100 多倍。
在最小二乘法中，只有一个及其简单的 for 循环，很少涉及内存写。但在线性代数中，需要进行矩阵的生成DenseMatrix.OfArray，以及矩阵运算，这二者都需要内存写。再者，矩阵计算有着繁重的计算量，这些都在影响着线性代数拟合圆的效率。最终的胜利还是属于最小二乘法。

修改特定位置像素值

#include<opencv2\opencv.hpp>
int main() 
{
	cv::Mat img = cv::imread("Lena.jpg");
	cv::Mat img_backup = img.clone();
	//如果是单通道，转为三通道
	if (1 == img.channels())
		cv::cvtColor(img, img, cv::COLOR_GRAY2BGR);
 
	/**
	核心：
	*	*(img.data + img.step[0] * i + img.step[1] * j + img.elemSize1() * c)=new_val
	*	这行代码是解析出来image的第i行第j列(即坐标为[j,i])第c通道，然后就可以对它进行赋值了。
	*/
	for (int i = 0; i < img.rows; ++i)
		for (int j = 0; j < img.cols; ++j)
		{
			//对第0个通道进行赋值操作，也就是blue通道
			*(img.data + img.step[0] * i + img.step[1] * j + img.elemSize1() * 0) = 255;
 
			//对第2个通道进行赋值操作，也就是green通道
			*(img.data + img.step[0] * i + img.step[1] * j + img.elemSize1() * 1) = 0;
 
			//对第3个通道进行赋值操作，也就是red通道
			*(img.data + img.step[0] * i + img.step[1] * j + img.elemSize1() * 2) = 0;
		}
 
	cv::imshow("img_src", img_backup);
	cv::imshow("img_rst", img);
	cv::waitKey();
}

OpencvSharp实现Mat定位设置值

//用C#中的get和set方法，速度较慢

for(int i = 1; i < rows -1; i++)
            {
                for(int j = 1; j < cols -1; j++)
                {
                    v = (int)Math.Abs(k[0, 0] * src.Get<byte>(i - 1, j - 1) + k[0, 2] * src.Get<byte>(i - 1, j + 1) + k[0, 1] * src.Get<byte>(i - 1, j) + k[2, 1] * src.Get<byte>(i + 1, j));
                    v =(int)(255 - 2 * v + 2 * k[1, 1] * src.Get<byte>(i, j));
                    v = v > 0 ? v : 0;
                    v = v < 255 ? v : 255;
                    dst.Set(i, j, v);
                }
             }

OpencvSharp 像素与通道操作
一：API：
1: AT();获取像素值，4个重载，用法大同小异

public T At<T>(int i0, int i1) where T : struct;  //返回指定数组元素的值。
public T At<T>(int i0) where T : struct;
public T At<T>(params int[] idx) where T : struct;
public T At<T>(int i0, int i1, int i2) where T : struct;

i0: Index along the dimension 0 （沿着空间维度0索引）
i1: Index along the dimension 1 （沿着空间维度1索引）

2: Get();与 AT方法一样

public T Get<T>(int i0, int i1) where T : struct;
public T Get<T>(int i0, int i1, int i2) where T : struct;
public T Get<T>(int i0) where T : struct;
public T Get<T>(params int[] idx) where T : struct;

3: Set();对指定位置的像素赋值 4个重载方法用法大同小异

//将值设置为指定的数组元素。
 public void Set<T>(int i0, int i1, T value) where T : struct;
 public void Set<T>(int i0, int i1, int i2, T value) where T : struct;
 public void Set<T>(int[] idx, T value) where T : struct;
 public void Set<T>(int i0, T value) where T : struct;

i0: Index along the dimension 0 （沿着空间维度0索引）
i1: Index along the dimension 1 （沿着空间维度1索引）
value: 像素值，泛型

二：单通道图像像素处理
代码：

static void Main(string[] args)
        {
            string imagePath = @"d:\ZG190016\1.jpg";
            ReadImageSingleChannelPixel(imagePath);
        }
        /// <summary>
        /// 单通道像素操作
        /// </summary>
        public static void ReadImageSingleChannelPixel(string path)
        {
            using (Mat src = new Mat(path, ImreadModes.AnyColor | ImreadModes.AnyDepth))
            using (Mat dst = new Mat())
            {
                //ColorConversionCodes 色彩空间转换枚举类型 BGRA2GRAY 转为灰度单通道
                Cv2.CvtColor(src, dst, ColorConversionCodes.BGRA2GRAY);
                Mat gray = new Mat();
                dst.CopyTo(gray);

                int height = dst.Rows; //行 获取图片的行叠加起来就是高度
                int width = dst.Cols;  //列  ... ...     宽度

                for (int row = 0; row < height; row++)
                {
                    for (int col = 0; col < width; col++)
                    {
                       byte p=  dst.At<byte>(row, col); //获对应矩阵坐标的取像素
                        byte value =byte.Parse( (255-p).ToString()); //反转像素值
                        dst.Set(row, col, value); //赋值
                    }
                }

                using(new Window("gray",WindowMode.FreeRatio,gray)) //单通道图像
                using(new Window("dst",WindowMode.FreeRatio,dst)) //反转后的单通道图像
                using (new Window("src", WindowMode.FreeRatio, src)) //源图
                {
                    Cv2.WaitKey(0);
                }
            }
        }

二：三通道图像像素处理
代码：

static void Main(string[] args)
        {
            string imagePath = @"d:\ZG190016\1.jpg";
            ReadImageThreeChannelsPixel(imagePath);
        }
        /// <summary>
        /// 三通道像素操作
        /// </summary>
        public static void ReadImageThreeChannelsPixel(string path)
        {
            using (Mat src = new Mat(path, ImreadModes.AnyColor | ImreadModes.AnyDepth))
            using (Mat dst = new Mat(src.Size(), src.Type()))
            {

                int height = src.Rows;
                int width = src.Cols;
                int cn = src.Channels(); //获取通道数

                #region 反转像素算法
                
                for (int row = 0; row < height; row++)
                {
                    for (int col = 0; col < width; col++)
                    {
                        if (cn == 1) //如果是单通道
                        {
                            byte p = dst.At<byte>(row, col); //获取像素
                            byte value = byte.Parse((255 - p).ToString()); //反转像素值
                            dst.Set(row, col, value); //赋值
                        }
                        else if(cn==3) //如果是三通道
                        {
                            //读取源图的像素
                            int b = src.At<Vec3b>(row, col)[0];
                            int g = src.At<Vec3b>(row, col)[1];
                            int r = src.At<Vec3b>(row, col)[2];

                            Vec3b color = new Vec3b
                            {
                                Item0 = (byte)(255 - b), //反转像素   (byte)( Math.Max(r, Math.Max(b, g)));
                                Item1 = (byte)(255 - g), //                 (byte)(Math.Max(r, Math.Max(b, g)));
                                Item2 = (byte)(255 - r) //                   (byte)(Math.Max(r, Math.Max(b, g)));
                            };
                            /*  
                             Vec3b color = new Vec3b //反转像素
                            {
                               Item0 = (byte)Math.Abs(src.Get<Vec3b>(row, col).Item0 - 255),
                               Item1 = (byte)Math.Abs(src.Get<Vec3b>(row, col).Item1 - 255),
                               Item2 = (byte)Math.Abs(src.Get<Vec3b>(row, col).Item2 - 255)
                            };
                            */
                            //Math.Max(r, Math.Max(b, g));取灰度，min也可以，但是亮度不同

                            //赋值
                            dst.Set<Vec3b>(row,col,color );
                        }
                       
                    }
                } 
                #endregion
                using (new Window("dst", WindowMode.FreeRatio, dst)) //反转
                using (new Window("src", WindowMode.FreeRatio, src)) //源图
                {
                    Cv2.WaitKey(0);
                }
            }
        }

Vec3b：通道顺序 blue grenn red 在C++中对应uchar类型，C#对应byte类型:
像素数据结构 BGR 包含三个像素值 R代表红，red; G代表绿，green; B代表蓝，blue。RGB模式就是，色彩数据模式，R在高位，G在中间，B在低位。BGR正好相反。例如，如果色彩数据是24位，对于RGB模式，就是高8位是R，中间8位是G，低8位是B。一个色彩数据共24位，3个字节。

Vec3f:对应三通道的float类型。
把CV_8UC1 转换到 CV32F1实现:

src.convertTo(dst,CV_32F1);

反转像素API：

 Cv2.BitwiseNot(src, dst); //反转像素函数，不需要操作像素，达到的效果一样

~ 取反符号：输出同样效果：

                using (new Window("dst", WindowMode.FreeRatio, ~src)) //反转
                using (new Window("src", WindowMode.FreeRatio, src)) //源图
                {
                    Cv2.WaitKey(0);
                }

OpenCV 基于距离变换的高精度轮廓匹配

轮廓匹配在定位测量应用中对其匹配的精度有更高的要求，通常的像素级的匹配结果难以满足其要求。本文给出了一种具有亚像素精度的快速轮廓匹配定位方法，其进行数学计算的基础为二值图像的距离变换。

二值图像距离变换的概念由Rosenfeld和Pfaltz于1966年其论文中提出，目前广泛应用于计算机图形学，计算机视觉及GIS空间分析等领域，其基本含义是计算一个图像中非零像素点到最近的零像素点的距离，也就是对每一各非零像素点计算其到零像素点的最短距离，并将该距离值赋值给该非零像素位置，从而将一幅二值图像变换为一幅距离图像。OpenCV通过cv：：distanceTransform()函数，给出了该功能的快速实现方法。并通过参数的方式给出了距离的几种定义方式，如欧式距离、马氏距离、倒角距离等。

高精度轮廓匹配的算法实现步骤：

假设待匹配的两个轮廓分别为A和B

（1）将轮廓A栅格化并生成为一幅二值图像Ia，其中轮廓点为255，背景点为0；

（2）将二值图Ia进行距离变换得到距离图Id；

（3）通过一个变换描述轮廓B和轮廓A之间的关系，如平移变换、刚体变换、仿射变换等；

（4）通过一定的搜索策略或者最优化方法搜索轮廓B在距离图Id中的最佳变换参数，并输出。

在上述过程中（4）是最关键的一步，具体的搜索策略可考虑使用Powell搜索、分支定界搜索或者Gauss-Newton方法。下面给出了使用Gauss-Newton方法实现求解平移变换的代码，仅供参考。其它复杂的变换可参考此过程完成。

1 数据栅格化，并距离变换

int GetRasterDistanceImage(cv::Point2f* refPoints, int pntsNum, double resolution, 
						   int& offX, int& offY, cv::Mat& rasterDistImg)
{
	double maxX, minX, maxY, minY;
 
	int i;
	for (i = 0; i < pntsNum; ++ i)
	{
		maxX = refPoints[i].x;
		minX = refPoints[i].x;
		maxY = refPoints[i].y;
		minY = refPoints[i].y;
 
		break;
	}
 
	for ( ; i < pntsNum; ++ i)
	{
		if (refPoints[i].x > maxX)
			maxX = refPoints[i].x;
		if (refPoints[i].x < minX)
			minX = refPoints[i].x;
 
		if (refPoints[i].y > maxY)
			maxY = refPoints[i].y;
		if (refPoints[i].y < minY)
			minY = refPoints[i].y;
	}
 
	int nX = int((maxX - minX) / resolution) + 200;
	int nY = int((maxY - minY) / resolution) + 200;
 
	cv::Mat rasterImg;
	rasterImg.create(nY, nX, CV_8UC1);
	memset(rasterImg.data, 0xFF, nY * nX);
 
	double realCenterX = (maxX + minX) / 2.0f;
	double realCenterY = (maxY + minY) / 2.0f;
 
	int offsetX = nX / 2 - int(realCenterX / resolution + 0.5f);
	int offsetY = nY / 2 - int(realCenterY / resolution + 0.5f);
 
	uchar* dat = (uchar*)(rasterImg.data);
	for (i = 0; i < pntsNum; ++ i)
	{
		if (validateMask[i] == 0)
			continue;
 
		int x = (int)floor(refPoints[i].x / resolution + 0.5f) + offsetX;
		int y = nY - ((int)floor(refPoints[i].y / resolution + 0.5f) + offsetY);
 
		*(dat + nX * y + x) = 0;
	}
 
	cv::distanceTransform(rasterImg, rasterDistImg, CV_DIST_C, 3);
 
	offX = offsetX;
	offY = offsetY;
 
	return 1;
}

2 高精度轮廓匹配，并输出平移量 offX，offY

void MatchReferencePoints(cv::Mat& rasterDistImg, cv::Point2f* rasterPoints, int pointsNum , double& offX, double& offY)
{
	offX = 0.0f;
	offY = 0.0f;
 
	cv::Mat imgB = rasterDistImg;
 
	double* pL = new double[pointsNum];
	double* pA = new double[pointsNum * 2];
 
	double a_b[2] = {0};
	double deltaX[2];
	double AL[2];
	double MM[4];
 
	double va0, va1, va2, va3, va4;
	double sx, sy;
	double pError = 0.0f;
 
	for (int i = 0; i < 100; ++ i)
	{
		double error = 0.0f;
		for (int j = 0; j < pointsNum; ++ j)
		{
			sx = myPoints[j].x + a_b[0];
			sy = myPoints[j].y + a_b[1];
 
			va0 = GetDataValue(imgB, sx, sy);
 
			va1 = GetDataValue(imgB, sx + 0.5, sy);
			va2 = GetDataValue(imgB, sx - 0.5, sy);
			va3 = GetDataValue(imgB, sx, sy + 0.5);
			va4 = GetDataValue(imgB, sx, sy - 0.5);
 
			pA[j] = va1 - va2;
			pA[j + pointsNum] = va3 - va4;
 
			pL[j] = 0.0f - va0;
 
			error += pL[j] * pL[j];
		}
 
		if (i == 0)
			pError = error;
		else if (error > pError)
			break;
 
		pError = error;
 
		int nOff_M = 0;
		for(int j = 0; j < 2; ++ j){
			for(int jj = 0; jj < 2; ++ jj){
				MM[nOff_M + jj] = 0.0;
				for(int kk = 0; kk < pointsNum; ++ kk)
					MM[nOff_M + jj] += pA[j * pointsNum + kk] * pA[jj * pointsNum + kk]; 
			}							// M = A * A(T)
			nOff_M += 2;
		}
 
		cv::Mat cvMatM(2, 2, CV_64FC1, MM);
		cv::invert(cvMatM, cvMatM);
 
		for (int j = 0; j < 2; ++ j)
		{
			AL[j] = 0.0f;
			for (int kk = 0; kk < pointsNum; ++ kk)
				AL[j] += pA[j * pointsNum + kk] * pL[kk];
		}
 
		for (int j = 0; j < 2; ++ j)
		{
			deltaX[j] = 0.0f;
			for (int jj = 0; jj < 2; ++ jj)
				deltaX[j] += MM[j * 2+ jj] * AL[jj];
		}
 
		int s = 0;
		for (; s < 2; ++ s)
		{
			if (fabs(deltaX[s]) > 0.0001)
				break;
		}
 
		if (s == 2)
			break;
 
		a_b[0] += deltaX[0];
		a_b[1] += deltaX[1];
	}
 
	offX = a_b[0];
	offY = a_b[1];
 
	delete[] pL;
	delete[] pA;
 
	delete[] myPoints;
 
	double temp = sqrt(pError / pointsNum);
}