大家注意，采点采的是Bitmap的像素坐标，不用转换到地理坐标！大家一定注意！舍友de了一天bug刚刚de出来快哭晕了~ ——2022/11/12修改
同学们，看完文章能不能点个赞，Part3、4、5加起来一千多的阅读量，没有一个人点赞😭😭。
还有，大家可以私信或者留言，不然我也不清楚下期要更什么。
好，今天是在赵同学强烈建议下的使用仿射变换进行图像配准，我们先来看demo。
一、Demo
这是图像配准前的两个图像，我们的目的呢，是把左图弄成右图的样式。接下来进行选点：
这里选了四个点，选好点之后，计算了他们的损失，就可以进行配准（这里使用最邻近像元法进行重采样）。
配准完成，因为有一些像素点重采样会采样到原图之外，所以采样到原图之外的地方我们设为黑色。
好，这就是图像配准的Demo。
二、思想
(能力强的同学可以根据思想自己手敲代码了，当然，解决方案不唯一)
从上面的Demo我们可以看到分为几个步骤：选点、计算损失、重采样、输出图像。
2.1 选点
选点这个部分呢在思想上没什么要讲的，无非就是读取鼠标的坐标嘛。但是有一点要记住，一定要把点按顺序对应起来，左边的第一个对的就是右边的第一个，不要搞混。
2.2 计算损失
其实这一“计算损失”的标题并不正确，这一部应该叫做“求解仿射变换参数并计算每个点的损失”。
首先，我们介绍一下仿射变换。
还是这个图，我们假设序号为1的点坐标如下（XL指的就是左图，R同理）：
| 左图点1 | 右图点1 | |
| X坐标 | XL1 | XR1 |
| Y坐标 | YL1 | YR1 |
我们知道左图是待矫正的点，想要它变成右图的样子。
那么本例中仿射变换的公式(仿射变换怎么来的我就不跟大家解释了哈)如下：
我们可以看到 变换后的(右图)XY坐标 分别由 变换前的XY坐标 和 一个常量 计算得出，在这个公式中，有abcdef六个参数，也就是说有6个未知数。那么一对点提供两个方程，六个未知数至少也要三对点。
接下来的问题就是，如果我们有三对以上的点，我们怎么通过C#编程，来实现最佳拟合求解参数呢？这就要有请我们的最小二乘了：
其中A就是我们 三对或以上XL与YL 组成的矩阵，b就是 XR和YR组成 的矩阵。
好，通过最小二乘，我们能得到了仿射变换六参数。接下来是计算损失，其实计算损失很简单，就以某一对点的XL和YL用求得的仿射变换六参数作变换，再以变换后的结果和这一对点的XR与YR作比较，即可得到差值，即是误差loss。
2.3 重采样
重采样呢，我这儿用了两种方法，一种是最邻近像元法，一种是双线性差值。
首先为了下面方便，我们把上面所说的待矫正的图称为“待矫正图”（废话），也就是上面说的左图；然后作为参考的图，也就是上面说的右图，称为“参考图”。那我们要生成一幅新图来接收图像配准的结果，我们就叫它“新图”。
新图应该是一片空白的，我们遍历新图上的每一个像素，通过仿射变换的逆变换，获得在待矫正图的坐标，然后在待矫正图上进行像素的采样，这就是重采样。
需要各种重采样的原因是：
由于我们逆变换获得的坐标存在小数，如图紫色点所示，而像素的位置都是整数（图中各个矩形的中心），所以我们获得的像素值应该取多少呢？
最邻近法就是取一个最近的像素的像素值。
双线性差值法就是以距离为权重，近的权重大，远的权重小，参考周围四个点来获得新的像素值。
至此，重采样完成。
2.4 输出图像
好，获得了新的图之后我们就可以输出，原理到此结束。
三、代码实现
3.1 最小二乘拟合类
上面有提到，我们需要最小二乘拟合，所以，在此制作了一个相关的类用来解决n元一次方程拟合的问题。另外，大家请先在C#的Nuget上安装MathNet，一个数学运算库。
class MultiX_1time_polynomial这是最小二乘拟合类的类名，多元_一次_多项式，当然大家也可以自己起一个自己喜欢的名字。接下来我们介绍它的类成员变量。（加上using System.Collections.Generic;）
private int NumFactors;
private List<double> Fts = new List<double>(); //存储系数的List
private List<double[]> Xs = new List<double[]>();
private List<double> Y = new List<double>(); //存储X和YNumFactors ：该方程的系数个数，比如在下面的方程中，系数就是abc，系数个数3个。
Fts ：Factors的意思，用以存储解出来的系数，比如存储上式的abc系数。
Xs ：存储输入的（X，Y）坐标，也就是变换前点的XY坐标，上式的XL、YL。
Ys ：存储输入的已知变换后的坐标，上式的XR。
这些变量就可以用来作最小二乘拟合：
A就是我们输入的XL和YL（也就是线性代数里说的系数阵），b就是XR、YR（常数阵），x就是解出来的待定系数abcdef。
好，这就是类成员变量，然后我们再来看类成员函数：
public MultiX_1time_polynomial(int NumFactors)//输入待定系数个数
{
this.NumFactors = NumFactors;
for (int i = 1; i <= NumFactors; i++)
{
this.Fts.Add(0.0);
}
}这是类的构造函数，其中在构造时需要输入一个参数，即待定系数个数，构造例子如：
MultiX_1time_polynomial MX1 = new MultiX_1time_polynomial(3) //3个待定系数
下一个函数：
public void Add_Xs_Y(double[] Xs,double Y)
{
if(Xs.Length != NumFactors)
{
System.Windows.Forms.MessageBox.Show("这是一个编程错误！多元一次多项式类中传入的X数量与系数数量不一致！", "错误！");
}
this.Xs.Add(Xs);
this.Y.Add(Y);
}这个函数就是添加X和Y的值，存储到成员变量中。假设我们有一组数据是：
52.3 = a * 14.5 + b * 20.1 + c
那这个函数就是这么用的(MX1定义引用上面说到的)：
MX1.Add_Xs_Y(new double[3]{14.5, 20.1, 1}, 52.3);//多一个1是因为c的系数已知为1
一组一组添加数据即可。来看下一个函数：
private void Matrix_solve()
{
var Y = new MathNet.Numerics.LinearAlgebra.Complex.DenseMatrix(this.Y.Count, 1);
var X = new MathNet.Numerics.LinearAlgebra.Complex.DenseMatrix(this.Xs.Count, this.NumFactors);
for (int i = 0; i < this.Y.Count; i++)
{
Y[i, 0] = this.Y[i];
for (int j = 0; j < this.NumFactors; j++)
{
X[i, j] = this.Xs[i][j];
}
}
var Paras = (X.Transpose().Multiply(X)).Inverse().Multiply(X.Transpose()).Multiply(Y);
for (int i = 0; i < this.NumFactors; i++)
{
this.Fts[i] = Paras.At(i, 0).Real;
}
}这个就是我们求解的核心了，你可以看到有不少的矩阵运算，这里就不给大家解释了，毕竟这也不是公共的接口，同学们用不到这个函数。来，最后一个Solve函数：
public double[] Solve()
{
if (this.Xs.Count < this.NumFactors ||this.Y.Count < this.NumFactors)
{
System.Windows.Forms.MessageBox.Show("这是一个编程错误！多元一次多项式数据个数不足无法拟合参数！", "错误！");
}
this.Matrix_solve(); //最小二乘拟合
double[] result = new double[this.Fts.Count];
for (int i = 0; i < this.Fts.Count; i++)
result[i] = this.Fts[i];
return result;
}这个函数呢，就是以double数组的方式求解并返回我们要求的待定参数啦！用法如下：
double[] paras = MX1.Solve()//在调用Solve前记得添加好适当个数的XY值哦
在我的程序中我是定义了两个MultiX_1time_polynomial来分别求解abc和def的，大家有更好的想法也可以修改一下代码。
好，最小二乘拟合类到此为止！
3.2 损失计算
在根据采点（这里界面的交互设计我就不讲了，每个人都有每个人的设计思路）求解出仿射变换六参数后，我们来看看使用仿射变换六参数进行损失的计算（似乎并没有什么好讲的，直接上代码）。
private double Compute_Loss(Point p1,Point p2,double[] TsF)//计算loss
{
double lossX = (TsF[0] * p1.X + TsF[1] * p1.Y + TsF[2]) - p2.X;
double lossY = (TsF[3] * p1.X + TsF[4] * p1.Y + TsF[5]) - p2.Y;
return System.Math.Sqrt(lossX * lossX + lossY * lossY);
}这里p1与p2就是我们所采的点，TsF就是仿射变换六参数的数组（按顺序abcdef），结果平方和开方即可。
3.3重采样
在合适的Loss的情况下我们进行矫正，重头戏就是重采样，这里介绍两种重采样方法的代码：
private Bitmap Execute_registration(Bitmap Pre_Bitmap,double[] TsF,int method=0)
{
//Pre_Bitmap是待矫正图像，TsF是仿射变换六参数，method为0是最邻近，1是双线性
Bitmap new_bitmap = new Bitmap(Pre_Bitmap);
for(int y = 0;y<new_bitmap.Height;y++)
for(int x = 0;x<new_bitmap.Width;x++)
{
//仿射变换逆变换
double[] preXY = this.Get_PreBMP_coor(x, y, TsF);
if(preXY[0] < 0 || preXY[0] > Pre_Bitmap.Width-1
|| preXY[1] < 0 || preXY[1] > Pre_Bitmap.Height - 1)
{//超出范围设为黑色
new_bitmap.SetPixel(x, y, Color.Black);
}
else
{//重采样
new_bitmap.SetPixel(x, y, this.Resample(preXY,Pre_Bitmap,method));
}
}
return new_bitmap;
}我们先以待矫正图像新建立一个Bitmap（其实图像配准应该输出的是和原Dataset一样的数据格式，我这边还没改过来，所以讲的就是输出新的BMP图像），遍历NewBitmap上的每一个像素，经过仿射变换逆变换（函数1）获得在原图中的坐标，若超出范围则采样为黑色，没有超出范围则进行重采样（函数2）。所以这里需要讲两个函数：1、仿射变换逆变换，2、重采样。先来第一个：
private double[] Get_PreBMP_coor(int preX,int preY,double[] TsF)//获得现在点在原图的坐标
{
var postXY = new MathNet.Numerics.LinearAlgebra.Complex.DenseMatrix(2, 1);
postXY[0, 0] = (double)preX;
postXY[1, 0] = (double)preY;
var cf = new MathNet.Numerics.LinearAlgebra.Complex.DenseMatrix(2, 1); //常数项
cf[0, 0] = TsF[2];
cf[1, 0] = TsF[5];
var abde = new MathNet.Numerics.LinearAlgebra.Complex.DenseMatrix(2, 2); //系数项
abde[0, 0] = TsF[0];abde[0, 1] = TsF[1];
abde[1, 0] = TsF[3];abde[1, 1] = TsF[4];
var preXY = abde.Inverse() * (postXY - cf); //计算得出原图坐标
return new double[2] { preXY.At(0,0).Real, preXY.At(1, 0).Real };
}这里的preX和preY就是新图的坐标，返回的就是在待矫正图上的坐标，其实在数学上挺简单的，纸笔写一下就知道要怎么算了，这里就不给大家细讲了，放进去用就OK。接下来就是Resample，我会讲一个函数分为多个部分来讲：
private Color Resample(double[] preXY,Bitmap bitmap,int method)
首先注意，这里返回的是Color类型，preXY就是逆变换求得的待矫正图坐标，bitmap就是待矫正图，method是选择方法。
if(method == 0)//最邻近像元法
{
int preX = ((int)System.Math.Round(preXY[0]));
int preY = ((int)System.Math.Round(preXY[1]));
return bitmap.GetPixel(preX, preY);
}最邻近像元法很简单，用Math.Round四舍五入就获得最邻近的值了，然后使用GetPixel获取Color就完成了。
else//双线性
{
int[] preX = new int[2] { (int)System.Math.Floor(preXY[0]), (int)System.Math.Ceiling(preXY[0]) };
int[] preY = new int[2] { (int)System.Math.Floor(preXY[1]), (int)System.Math.Ceiling(preXY[1]) };
double[] RGB = new double[3] { 0, 0, 0 };//存储三通道的结果
for(int i = 0;i<2;i++)
for(int j = 0;j<2;j++)
{
this.Bilinear_interpolation(bitmap.GetPixel(preX[i], preY[j]), new int[2] { preX[i], preY[j] }, preXY, ref RGB);
}
return Color.FromArgb(
(int)System.Math.Round(RGB[0]),
(int)System.Math.Round(RGB[1]),
(int)System.Math.Round(RGB[2])
);
}双线性插值法，首先计算与坐标邻近的四个像素的坐标，比如逆变换获得的X=12.3,Y=45.6，那么我们对他们进行取大和取小就有X=[12,13],Y=[45,46],组合起来就获得了2*2=4个坐标，也就是相邻的四个像素坐标。这一步就是在上面preX和preY实现。
然后我们可以看到两个循环，其实就是遍历X=[12,13]与Y=[45,46]，从而达到遍历周围的四个像素的效果，对每一个像素调用Bilinear_interpolation双线性插值函数，我们来看看这个函数：
private void Bilinear_interpolation(Color PixColor,int[] PixXY,double[] XY,ref double[] RGB)
{
//PixColor是该像素的Color值
//PixXY是像素的坐标，也就是说是两个整数；
//XY是逆变换获得的坐标，为小数。
//计算权重
double w_x = 1.0 - System.Math.Abs(XY[0] - (double)PixXY[0]);
double w_y = 1.0 - System.Math.Abs(XY[1] - (double)PixXY[1]);
//将结果放入RGB数组中
RGB[0] += PixColor.R * w_x * w_y;
RGB[1] += PixColor.G * w_x * w_y;
RGB[2] += PixColor.B * w_x * w_y;
}应该不难理解吧？类似于下图，XY就是黑点，PixXY就是粉色方块（像素）的中心坐标，PixColor就是该像素的Color值，我们通过计算黑点与像素的距离来求出采样的权重，然后将Color值乘以权重返回到RGB数组中，每一个坐标需要对周围四个像素进行采样。
到这儿，我们就完成了重采样的部分，Resample完成代码如下：
private Color Resample(double[] preXY,Bitmap bitmap,int method)
{
if(method == 0)//最邻近像元法
{
int preX = ((int)System.Math.Round(preXY[0]));
int preY = ((int)System.Math.Round(preXY[1]));
return bitmap.GetPixel(preX, preY);
}
else//双线性
{
int[] preX = new int[2] { (int)System.Math.Floor(preXY[0]), (int)System.Math.Ceiling(preXY[0]) };
int[] preY = new int[2] { (int)System.Math.Floor(preXY[1]), (int)System.Math.Ceiling(preXY[1]) };
double[] RGB = new double[3] { 0, 0, 0 };
for(int i = 0;i<2;i++)
for(int j = 0;j<2;j++)
{
this.Bilinear_interpolation(bitmap.GetPixel(preX[i], preY[j]), new int[2] { preX[i], preY[j] }, preXY, ref RGB);
}
return Color.FromArgb(
(int)System.Math.Round(RGB[0]),
(int)System.Math.Round(RGB[1]),
(int)System.Math.Round(RGB[2])
);
}
}好，Resample结束后，我们的图像配准也基本完成。
后话
当然这只是Bitmap的图像配准，输出的也只是Bitmap图像，我们需要做的应该是输入什么图像就输出相同格式的图像，这个我也还没来得及改，这两天在忙着做图像的融合，如果有时间改了再放上来啦。
好，今天的分享就到这里。
下一步想看更新什么或者说有什么建议都可以私信我或者在评论区留言哦。
725
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
