（MATLAB/C）高斯拟合法求光斑中心

光斑图、阵列图、灰度图圆形等目标中心定位方法。分享高斯拟合法和更为简单的中心、重心法MATLAB代码，以及基于Eigen库的高斯拟合法C代码。互助互助

by HPC_ZY

一、基本原理

大多数光斑其明暗分布情况都是中心最亮，往四周慢慢变暗，就类似二维高斯模型（如下图）。所以我们利用二维高斯模型去拟合光斑，从而得到光斑中心等参数。

由于本人不是数学大佬，就不推导数学公式了，直接上代码。

MATLAB版本

三种方法都需要提供一个叫mask的矩阵，它是一个二值图像，描述的是被判定而光斑的像素。务必提供争取或你认为对的mask，中心只会在mask的范围内搜索。

1. 高斯拟合法

% 输入：原始灰度图像，光斑二值蒙版
% 输出：中心坐标
function coor = gausscenter(im,mask)

% 连通域
[label,num] = bwlabel(mask);

% 计算
coor = zeros(num,2);
for n = 1:num
    [x,y] = find(label==n);
    % 生成计算矩阵
    m_iN = length(x);
    tmp_A = zeros(m_iN,1);
    tmp_B = zeros(m_iN,5);    
    for k = 1:m_iN        
        pSrc = im(x(k),y(k));
        if pSrc>0
            tmp_A(k) = pSrc*log(pSrc);
        end
        tmp_B(k,1) = pSrc ;
        tmp_B(k,2) = pSrc*x(k);
        tmp_B(k,3) = pSrc*y(k);
        tmp_B(k,4) = pSrc*x(k)*x(k);
        tmp_B(k,5) = pSrc*y(k)*y(k);
    end
    
    % QR分解
    Vector_A = tmp_A;
    matrix_B = tmp_B;
    [Q,R] = qr(matrix_B);
    
    % 求解中心
    S = Q'*Vector_A;
    S = S(1:5);
    R1 = R(1:5,1:5);
    C = R1\S;   
    coor(n,:) = [-0.5*C(2)/C(4),-0.5*C(3)/C(5)];
end

end

2. 重心法

% 输入：原始灰度图像，光斑二值蒙版
% 输出：中心坐标
function coor = gravitycenter(im,mask)

% 连通域
[M,~] = size(im);
[label,num] = bwlabel(mask);
coor = zeros(num,2);

for n = 1:num
    [x,y] = find(label==n);
    % 取出对应点灰度
    idx = (y-1)*M+x; % matlab是列优先
    imtmp = im(idx);
    % 计算权值
    w = (imtmp/sum(imtmp))';   
    coor(n,:) = [w*x,w*y];
end

end

3. 中心法

% 输入：光斑二值蒙版
% 输出：中心坐标
function coor = geometriccenter(mask)

% 连通域
[label,num] = bwlabel(mask);

% 计算
coor = zeros(num,2);
for n = 1:num
    [x,y] = find(label==n);
    coor(n,:) = [mean(x),mean(y)];
end

end

4. 测试代码

先别抄代码，听我说一句
先别抄代码，听我说一句
先别抄代码，听我说一句

第一，定位中心的前提是分割，调用算法之前，请先完成分割。 不要直接抄我的“获取蒙版”这一步，除非我们图一模一样。
第二，分割完之后，先imshow一个你的蒙版mask，看看是不是跟你想的一样（是不是所有光斑都是白色，其他都是黑色）
第三，主要太多人报错来问我，都是二值图有错误造成的，统一回复了。这篇文章目前没有包含图像分割的内容，更多报错请看文末“其他”

%% 读取图像

im = imread('test2.png'); % 这是我的测试图片，换成你的，如果是彩色你还得转为灰度
% 注意注意，我的图是灰度图，所以这样就ok了
im = im2double(im);
% 如果你是彩色图  请  im = im2double(rgb2gray(im));

%% 获取蒙版
% 要选择适合你图像的方法，对于复杂图像这一步并不简单
% 请务必使用合适的图像分割算法，有必要的请imshow(mask)，看看你的mask对不对
mask = im>0.5; % 由于我的图比较单纯，直接阈值。 请勿照搬，否则基本失败

%% 光斑中心定位
% 高斯拟合法
coor1 = gausscenter(im,mask);

% 重心法
coor2 = gravitycenter(im,mask);

% 中心法
coor3 = geometriccenter(mask);

%% 显示结果
imshow(im),hold on
plot(coor1(:,2),coor1(:,1),'r+','MarkerSize',15) % 高斯法
plot(coor2(:,2),coor2(:,1),'g.','MarkerSize',15) % 重心法
plot(coor3(:,2),coor3(:,1),'bo','MarkerSize',15) % 中心法
legend({'高斯法','重心法','中心法'})

注意：1）函数输入图像是灰度图，如果你的是彩色图请转为灰度图。2）Mask是指被你认定为是光斑的像素，请使用合适的图像分割方法获取mask，我代码中写的是阈值分割法（且使用了一个很简单的值0.5）

2023年3月21日补充：上面的coor1,coor2,coor3就是计算出的中心坐标，它们是N*2的向量（X,Y）。也就是行数是找到的光斑数量，每一行的第一列是x坐标，第二列是y坐标。

5. 实验结果
   可以看出对于完整的光斑，三种方法结果几乎一致。
   但对于边缘处光斑（不完整），高斯拟合法的优势就体现出来了。

C版本

主要利用Eigen的矩阵运算库 （Eigen下载）

1. 核心代码

int gausscenter(
	float *x,						// （out）x
	float *y,						// （out）y
	float *pimg,					// （in）原始采集图像
	int *imglabel,					// （in）原图标记
	int labeln,						// （in）光斑数量
	int imgWidth, 					// （in）图像宽
	int imgHeight					// （in）图像高
)
{
	// 预分配内存（缓存每个光斑像素位置，根据需要修改）
	int *xtmp = new int[2048]; // 因为我的光斑尺寸小，所以2048足以
	int *ytmp = new int[2048];

	// 遍历所有光斑
	for (int k = 1; k <= labeln; k++) 
	{
		int pn = 0;
		// 统计单个光斑坐标
		for (int i = 0; i < imgHeight; i++)
		{
			for (int j = 0; j < imgWidth; j++)
			{
				if (imglabel[i*imgWidth + j] == k)
				{
					xtmp[pn] = i;
					ytmp[pn] = j;
					pn++;
				}
			}
		} // 统计完毕

		// 存入矩阵
		MatrixXf X(pn, 1);
		MatrixXf Y(pn, 1);
		for (int i = 0; i < pn; i++)
		{
			X(i, 0) = xtmp[i];
			Y(i, 0) = ytmp[i];
		}
		// 生成计算矩阵
		MatrixXf A(pn, 1);
		MatrixXf B(pn, 5);
		for (int i = 0; i < pn; i++)
		{
			float img = pimg[(int)X(i, 0)*imgWidth + (int)Y(i, 0)];
			A(i, 0) = img*log(img);
			B(i, 0) = img;
			B(i, 1) = img*X(i, 0);
			B(i, 2) = img*Y(i, 0);
			B(i, 3) = img*X(i, 0)*X(i, 0);
			B(i, 4) = img*Y(i, 0)*Y(i, 0);
		}
		// QR分解
		HouseholderQR<MatrixXf> qr;
		qr.compute(B);
		MatrixXf R = qr.matrixQR().triangularView<Eigen::Upper>();
		MatrixXf Q = qr.householderQ();
		// 特征解
		MatrixXf S;
		S = Q.transpose()*A;
		MatrixXf S0(5, 1);
		MatrixXf R0(5, 5);
		for (int i = 0; i < 5; i++)
		{
			S0(i, 0) = S(i, 0);
			for (int j = 0; j < 5; j++)
			{
				R0(i, j) = R(i, j);
			}
		}
		MatrixXf C = R0.inverse()*S0;

		x[k-1] = -0.5 * C(1) / C(3);
		y[k-1] = -0.5 * C(2) / C(4);
	}
	
	return 0;
}

2. 调用方法

#include <Eigen/Dense>
using namespace Eigen;

int main()
{

float *img = new float[widt*height];	
int *label = new int[width*height]; 	// 连通域结果
int pn = 0; 							// 用来保存连通域个数

// 省略读图+计算连通域
// 加油
// 我就不写了

float *x = new float[pn];
float *y = new float[pn];
getlightcenter(x, y, img, label, pn, width, height);

return 0;
}

其他

1. 代码都公开了，就不上传文件了
2. 如果我的代码有帮助到小伙伴们，不妨点赞留言
3. 答疑（在提问前可先看以下内容）
   1）如果报错“内存炸了”，请检查自己的二值图像，大多数人的电脑是无法处理面积大于16000个像素的光斑。不过话也说回来，你的光斑真的那么大那么清晰，也用不着这个算法了，直接质心。
   2）如果报错“S（1:5）巴拉巴拉”，检查你的二值图像，是不是面积太小了，小到没有了。
   3）如何生成光斑图，最快捷的方式就是用PS软件自己画。
   4）光斑中心定位的前提是光斑分割，如果分割这一步（也就是生成二值图）都做错了，后面一定会错的。遇到问题优先检查二值图，确定没问题还报错，再问我就好。