C#（OpenGL MathNet）处理Gauss光斑图像

C# 处理高斯光束的光斑图像

1 基础操作
2 图片截取
3 转灰度图
4 SharpGL画三维点云图
5 MathNet拟合，OxyPlot作图

1 基础操作

.Net平台必备VS，新建WinForm项目，项目名称Gauss，位置任选。然后就会进入窗口编辑页面，我们开始拖控件。

考虑到我们的需求无非是

读取图片
图片转灰度
展示图像灰度
高斯拟合

所以先排上四个按钮，分别用于这四种需求，如下图所示

然后我们首先想办法满足第一个功能，即打开图片。双击打开图片的按钮，进入函数编辑页面，对btnOpen_Click进行编辑。这个函数顾名思义，即当我们点击按钮btnOpen时将要执行的操作。

我们首先通过一个文件对话框OpenFileDialog来读取文件路径。

为了验证我们的代码，可以在Form1面板上添加一个textBox，其添加方式与button如出一辙。添加之后，将其(Name)改为textMessage，用于存放调试信息。

则得到代码如下

private void btnOpen_Click(object sender, EventArgs e)
{string strFileName = "";    //用于存储图片路径OpenFileDialog ofd = new OpenFileDialog();ofd.Filter = "图片文件(*.bmp)|*.bmp";if (ofd.ShowDialog() == DialogResult.OK)strFileName = ofd.FileName;textMessage.Text = strFileName;
}

按下F5运行，然后点击打开图像按钮，可以进入我们熟悉的文件对话框，选择图片后可以看到在textMessage中返回了图片路径。

接下来，开始进行读图操作，在C#中，提供了Bitmap的读取操作，我们可以直接在获取图像路径之后，通过Image.FromFile()函数进行图片的读取。

而图片读取之后，则需要进行展示，为此我们同样从工具箱中拖取PictureBox控件，来展示我们打开的图像，将其(Name)改为Facula，然后修改btnOpen_Click函数如下

private void btnOpen_Click(object sender, EventArgs e)
{string strFileName = "";OpenFileDialog ofd = new OpenFileDialog();ofd.Filter = "图片文件(*.bmp)|*.bmp";if (ofd.ShowDialog() == DialogResult.OK)strFileName = ofd.FileName;Bitmap facula = (Bitmap)Image.FromFile(strFileName);Facula.Image = facula.Clone() as Image;
}

运行之后，选择我们将要处理的图像，运行结果为

2 图片截取

由于光斑图片大部分是背景，故需截取感兴趣的区域。我们希望实现一种操作逻辑，即通过鼠标框选感兴趣区域，从而实现截图功能。在此可以把框选分为两个过程，即鼠标按下，然后拖动鼠标，最后鼠标松开。

所以需要封装两个函数，双击解决方案资源管理器中的Form.Designer.cs，进入代码编辑页面，Ctrl+F找到Facula的位置，在下面添加

this.Facula.MouseDown += new System.Windows.Forms.MouseEventHandler(this.Facule_Down);
this.Facula.MouseUp += new System.Windows.Forms.MouseEventHandler(this.Facule_Up);

其中，MouseDown表示按下鼠标左键触发的行为；MouseUp表示松开鼠标左键时的行为。只要我们分别可以记下按下和松开时的鼠标位置，就可以返回一个矩形框。

在C#中，Control.MousePosition表示当前鼠标相对于屏幕左上角的位置；this.Location表示窗体左上角相对于屏幕左上角的位置；Facula.Location表示我们创建的控件Facula相对于窗体左上角的位置。

考虑到我们需要在两个动作中调用同一参数，所以最好新建两个全局变量，存储按下鼠标时的坐标。然后新建一个函数，用于处理当前鼠标位置，并将其转化为图像坐标。

Bitmap facula;      //存放图像
int xAxis,yAxis;    //在后面调用
private float getAxis(char flag)
{float axis;if(flag == 'x')axis = (int)(Control.MousePosition.X - this.Location.X - Facula.Location.X-10)/ Facula.Width * facula.Width;elseaxis = (float)(Control.MousePosition.Y - this.Location.Y - Facula.Location.Y-38) / Facula.Height* facula.Height;return axis;
}//鼠标按下时的操作
private void Facule_Down(object sender, EventArgs e)
{this.xAxis = getAxis('x');this.yAxis = getAxis('y');
}
//鼠标抬起时的操作
private void Facule_Up(object sender, EventArgs e)
{//先放在这一会儿再写
}

注意，我们所有的代码都写在public partial class Form1 : Form这个类里面。其中getAxis可以通过三元表达式写成更简洁的形式

private int getAxis(char flag)
{return flag=='x' ? (Control.MousePosition.X - this.Location.X - Facula.Location.X - 10)* facula.Width / Facula.Width: (Control.MousePosition.Y - this.Location.Y - Facula.Location.Y - 38)* facula.Height / Facula.Height;
}

在C#中，提供了Clone方法，可以实现Bitmap类的裁剪功能，其输入参数为一个Rectangle，即矩形，代表我们将要分割的感兴趣区域。

private void Facule_Up(object sender, EventArgs e)
{int xStart = Math.Min(getAxis('x'), xAxis);int yStart = Math.Min(getAxis('y'), yAxis);int width = Math.Abs(getAxis('x') - xAxis);int height = Math.Abs(getAxis('y') - yAxis);//此即我们感兴趣的区域Rectangle roi = new Rectangle(xStart,yStart,width,height);facula = facula.Clone(roi,facula.PixelFormat);//图像裁剪Facula.Image = facula.Clone() as Image;        //更改图片显示
}

由于裁剪前后图像尺寸相去甚远，所以需要将我们创建的PictureBox控件Facula的SizeMode设为ScratchImage。

其效果如图所示

3 转灰度图

我们看到的光斑图像虽然是黑白的，但并不是一个灰度图，因为这张图片仍然有四个通道，分别存放r, g, b以及alpha四个值，只不过透明度为0，且r, g, b相等，所以自然没有色彩。

为了将其转化为真·灰度图，我们需要稍微了解一下位图的编码方式。一张图片，主要包含两个部分，即文件头与数据，文件头一般占据54个字节，声明这是一张图片，并告诉我们这张图片的长宽信息；而其色彩信息则线性地存放在数据区里。图片软件通过文件头，获知这是一张图片，再根据其长宽，对数据进行扫描，呈现到屏幕上，我们就看到图了。

这里面有两个问题，

由于图片的像素值的大小是有限的，一般是8位无符号整型；但文件在操作系统中可能以64位为一个最小的存储单元，所以图片中的一行未必能够正好塞满64的倍数个bit，所以需要有一个数据对齐的问题。
红绿蓝三原色如何反演出灰度图像，直接取平均是否可行？对于我们的光斑图像来说是无所谓的，但比较流行的方案是gray=0.11r+0.59g+0.3b\text{gray}=0.11r+0.59g+0.3bgray=0.11r+0.59g+0.3b，我们的代码中也采取这个形式。

接下来，双击转灰度按钮，进入btnGray_Click函数的编辑位置，其代码为

byte[,] matFacula;        //创建一个全局变量来存储强度数据
private void btnGray_Click(object sender, EventArgs e)
{int width = facula.Width;     //图片宽度int height = facula.Height;   //图片高度//创建数据交换区，用来读取facula中的数据var bmdata = facula.LockBits(new Rectangle(Point.Empty, facula.Size),ImageLockMode.ReadOnly, PixelFormat.Format32bppArgb);byte[] buffer = new byte[bmdata.Stride* bmdata.Height]; Marshal.Copy(bmdata.Scan0, buffer, 0, buffer.Length); //复制图像数据facula.UnlockBits(bmdata);  //bmdata解锁matFacula = new byte[height, width]; //新建一个8位灰度图像，宽高分别是width,heightfacula = new Bitmap(width, height, PixelFormat.Format8bppIndexed);//这回用来写facula中的数据bmdata = facula.LockBits(new Rectangle(Point.Empty, facula.Size),ImageLockMode.WriteOnly, PixelFormat.Format8bppIndexed);int stride = bmdata.Stride;      //扫描行字节数byte[] grayValues = new byte[stride * height];   // 为图像数据分配内存int op = 0;    //用来遍历rgb与alpha通道for (int i = 0; i < height; i++)for (int j = 0; j < width; j++){matFacula[i, j] = (byte)(buffer[op++] * 0.11 + buffer[op++] * 0.59 + buffer[op++] * 0.3);op++;//跳过alpha通道grayValues[i*stride+j] = matFacula[i, j];}//参数分别是被复制数据,初始位置,目标指针，目标长度Marshal.Copy(grayValues, 0, bmdata.Scan0, stride * height);facula.UnlockBits(bmdata);  // 解锁内存区域  // 修改生成位图的索引表ColorPalette palette = facula.Palette;for (int i = 0; i < 256; i++)palette.Entries[i] = Color.FromArgb(i, i, i);facula.Palette = palette;Facula.Image = facula.Clone() as Image;
}

其中PixelFormat.Format32bppArgb代表32位Argb图像，A即alpha，代表透明度；rgb分别代表红绿蓝，四个通道每个通道都是8位，合在一起就是32位；PixelFormat.Format8bppIndexed代表8位索引图像，即只有一个通道，但这个通道可以根据颜色表进行对照，我们想做出灰度图像，所以颜色表中的rgb值应该是1:1:11:1:11:1:1。

Marshal是System.Runtime.InteropServices中的一个类，可进行一些内存操作，其中Copy函数即将一个内存区域复制给另一个变量。以Marshal.Copy(bmdata.Scan0, buffer, 0, buffer.Length);为例，其复制的起始位置为bmdata.Scan0，将复制给buffer，偏移量为0，复制buffer.Length个字节。

ColorPalette是位于System.Drawing.Imaging中的一个类，可以定义调色板的颜色的数组，可以理解为位图的一个属性。

最终效果如下图所示

在以上代码中，有一个貌似没用上的变量matFacula，我们将用这个二维数组进行一些数学计算，并绘制强度图。

4 SharpGL画三维点云图

这里所谓的三维图，目的是把光斑的强度分布展示出来，很遗憾C#在三维作图这一点上远远不如专用的科学计算语言，但方法总比问题多，我们选择Sharp GL，即OpenGL的C#版本，来进行三维图的绘制。

点击VS菜单栏的工具→NuGet包管理器→管理解决方案的Nuget包，选择浏览，然后搜索SharpGL，选择搜索结果中的SharpGL以及SharpGL.WinForms进行安装。

安装成功后，工具箱中会出现SharpGL.WinForms组，里面有一些SharpGL控件。

细心的朋友可能会发现我刚刚展示的图片中多了一个选项卡，即光斑图像和画图，这个看上去很高大上的东西也能在工具箱中找到，是Continer中的TabControl，一时之间找不到的话可以直接在工具箱中搜索。将选项卡控件拖放到窗口中后，可以看到这个控件右上角有一个侧三角图标，可以进行添加或删除Page的操作，很符合人的直觉；改名字之类的操作与其他控件相同，不再赘述。

我们把之前的PictureBox控件拖入光斑图像选项卡中，然后在画图选项卡中拖入OpenGLControl(SharpGL)控件，取个名字叫glMesh。然后点击三维图按钮，进入btnMesh_Click函数的编辑位置。

  private void btnMesh_Click(object sender, EventArgs e){OpenGL gl = this.glMesh.OpenGL;   //新建openGL对象//清除缓冲区gl.Clear(OpenGL.GL_COLOR_BUFFER_BIT | OpenGL.GL_DEPTH_BUFFER_BIT);gl.MatrixMode(OpenGL.GL_PROJECTION);//设置工作模式gl.LoadIdentity();                    //生成单位阵gl.Perspective(65, 1.0, 1.0, 1000);  // 投影矩阵（ y 45度，纵横比1：1，near截平面1，far截平面10000.）gl.LookAt(-10, -10, 150, 0.5*(double)facula.Width, 0.5 * (double)facula.Height, 0.0f, 1.0f, 1.0f, -1.0f); //defines a viewing transformationfloat z;//用来存放像素点的强度信息//查找matFacula的最大值for (int i = 0; i < facula.Height; i++)for (int j = 0; j < facula.Width; j++)zMax = zMax > matFacula[i,j] ? zMax : matFacula[i,j];gl.Begin(OpenGL.GL_POINTS);   //开始画图for (int x = 0; x < facula.Height; x += 1)for (int y = 0; y < facula.Width; y += 1){z = (float)matFacula[x,y];gl.Color(z/zMax,0.3,0);        //设置点的颜色gl.Vertex(x,y, z);                //绘制点的位置}gl.End();         //画图完毕gl.Flush();       //更新图像
}

由于电脑屏幕无论如何也只能进行二维的显示，所以模型做好之后，需要用一个“相机”把模型拍下来，使之成为一张图片显示在屏幕上。

其中，gl.Perspective的用途是定义相机的性能，输入的四个参数分别代表：视场角、纵横比，后面两个代表焦深，即相机能看到的距离范围。

gl.LookAt代表了相机的位置，在OpenGL中的定义为

void gluLookAt(GLdouble eyex, GLdouble eyey, GLdouble eyez,GLdouble centerx, GLdouble centery, GLdouble centerz,GLdouble upx, GLdouble upy, GLdouble upz)

其中，eyex, eyey, eyez代表相机所在位置；centerx, centery, centerz则代表相机正对着的位置；upx, upy, upz表示相机平面的角度。

我们可以根据光斑图像的坐标范围来确定这些参数，由于像素点的值最大不超过255，而图片坐标的起始点为(0,0)，所以我们选定(-10,-10,150)作为相机位置；相机正对着的点为图像中心处，即其z向坐标为0，x,y向坐标为图像长宽的一半；最后的视角，我们选择斜向下45度，即(1,1,-1)。

gl.Color为设置点的颜色，需要注意的是，这里的输入值虽然可以是rgb，但其取值范围是0到1，所以我们先搜索出图像矩阵的最大值，然后根据当前值与最大值的比值来设置图像颜色。

根据我们的设置gl.Color(z/zMax,0.3,0)，当像素的灰度值越大时，图像越红，否则越绿。

运行之后，拖动裁剪，然后转灰度，再绘制三维图，得到的结果为

5 MathNet拟合，OxyPlot作图

拟合曲线之前，我们需要找到一群可以被拟合的点。由于已经得到了光斑的灰度矩阵，所以我们可以从每一行选出一个最大值来代表这行光斑，从而得到一个一维的序列。

方法很简单

//private void btnFit_Click(object sender, EventArgs e)
double[] xArray = new double[facula.Height];
double[] pixArray = new double[facula.Height];for (int i = 0; i < facula.Height; i++)
{xArray[i] = i + 1;pixArray[i] = 0;for (int j = 0; j < facula.Width; j++)pixArray[i] = Math.Max(pixArray[i], matFacula[i, j]);
}

接下来就是问题之关键——如何将xArray和pixArray逆合成一条高斯曲线。但在这个关键问题之前，我们还很好奇那些被选出来的点到底是什么样的，希望把这两组数先画出来看一看。

尽管SharpGL完全可以胜任这一工作，但相比之下，我们可以选择更适合的画图工具来完成，比如OxyPlot。那么接下来就是同样的套路

在NuGet包管理器中搜索OxyPlot，然后安装OxyPlot.Core以及OxyPlot.WindowsForms
在tabPages中添加一个选项卡，起个名字叫拟合
从工具箱中选择PlotView添加到拟合选项卡上，将其更名为curveView
using OxyPlot，using OxyPlot.Series

然后在btnFit_Click函数中继续写

//接在这一行后面
//pixArray[i] = Math.Max(pixArray[i], matFacula[i, j]);
//}
curveView.Model = new PlotModel { };
var scatters = new ScatterSeries()
{Title = "原始数据",MarkerType = MarkerType.Circle,MarkerSize = 1
};
for (int i = 0; i < facula.Height; i++)scatters.Points.Add(new ScatterPoint(xArray[i], pixArray[i]));curveView.Model.Series.Add(scatters);
curveView.Model.InvalidatePlot(true);

运行结果如图所示

看上去的确有一点高斯的样子，所谓高斯函数，其表达形式为

y=a⋅exp⁡(−(x−bc)2)y=a⋅\exp{(−(\frac{x−b}{c})^2)}y=a⋅exp(−(cx−b)2)

其中，a的值表示该函数的最大值；b表示其中心值，c表示当y值降到1e2\frac{1}{e^2}e21分之一处时x距离中心的位置。

为了拟合这个函数，我们还需要再安装一个包，MathNet.Numerics，这个包里有一个可以对"任意"函数进行非线性拟合的工具，但悲催的是，目前只支持两个参数的拟合，而我们想要拟合的函数里有三个参数。

所以我们有两种方案，其一是做一些预处理，先消掉一个参数，例如对数据做归一化，从而消去强度项a，或者对数据取质心，从而将其移到原点，消去参数b。另一个方案则是通过某种变形，使得高斯函数变成某种多项式的形式。一个最直观的想法就是取对数。

ln⁡y=ln⁡a−(x−b)2c2ln⁡y=−x2c2+2bxc2+ln⁡a−b2c2\begin{aligned} \ln y&=\ln a-\frac{(x-b)^2}{c^2}\\ \ln y&=-\frac{x^2}{c^2}+\frac{2bx}{c^2}+\ln a-\frac{b^2}{c^2} \end{aligned} lnylny=lna−c2(x−b)2=−c2x2+c22bx+lna−c2b2

令Y=ln⁡y,A=−1c2,B=2bc2,C=ln⁡a−b2c2Y=\ln y, A=-\frac{1}{c^2}, B = \frac{2b}{c^2}, C=\ln a-\frac{b^2}{c^2}Y=lny,A=−c21,B=c22b,C=lna−c2b2，则Gauss函数变为

Y=Ax2+Bx+CY=Ax^2+Bx+C Y=Ax2+Bx+C

待拟合成功后，有c=1−A,b=−B2A,a=exp⁡B2+4AC4Ac=\frac{1}{\sqrt{-A}},b=-\frac{B}{2A},a=\exp\frac{B^2+4AC}{4A}c=−A1,b=−2AB,a=exp4AB2+4AC。

又考虑到我们截图的过程必然不会非常完美，于是会产生一些噪声，所以在正式做数据处理之前，需要先去下噪声，这里只用最简单的阈值方式，即只有大于阈值的点才可以参与到拟合中来。

最后将btnFit_Click函数分拆整理一下

private void btnFit_Click(object sender, EventArgs e)
{var data = showOrigin();curveFit(data);
}private double[] showOrigin()
{int pNum = facula.Height;            int[] xArray = new int[pNum];double[] pixArray = new double[pNum];for (int i = 0; i < pNum; i++){xArray[i] = i;pixArray[i] = 0;for (int j = 0; j < facula.Width; j++)pixArray[i] = Math.Max(pixArray[i], matFacula[i, j]);}curveView.Model = new PlotModel { };var scatters = new ScatterSeries(){Title = "原始数据",MarkerType = MarkerType.Circle,MarkerSize = 1};foreach (var item in xArray)scatters.Points.Add(new ScatterPoint(item, pixArray[item]));curveView.Model.Series.Add(scatters);curveView.Model.InvalidatePlot(true);return pixArray;
}private void curveFit(double[] yArray)
{int pNum = 0;int thres = 10;foreach (var item in yArray)if (item> thres)pNum += 1;double[] xFit = new double[pNum];double[] yFit = new double[pNum];int j = 0;for (int i = 0; i < yArray.Length; i++)if (yArray[i]>thres){xFit[j] = i;yFit[j++] = Math.Log(yArray[i]);}double[] fit = MathNet.Numerics.Fit.Polynomial(xFit, yFit, 2);var plotView = new LineSeries() { };plotView.Title = "拟合结果";plotView.Color = OxyColors.Red;for (int i = 0; i < yArray.Length; i++)plotView.Points.Add(new DataPoint(i, Math.Exp(fit[0] + fit[1] * i + fit[2] * Math.Pow(i, 2))));double c = -1 / Math.Sqrt(-fit[2]);double b = -fit[1] / fit[2] / 2;double a = Math.Exp((fit[1] * fit[1] + 4 * fit[0] * fit[2]) / fit[2] / 4);textMessage.AppendText(string.Format("y={0:F2}*exp(-((x-{1:F2})/{2:F2})^2)", a, b, c));curveView.Model.Series.Add(plotView);curveView.Model.InvalidatePlot(true);
}

结果为