一、基本原理

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

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

MATLAB版本

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

高斯拟合法

% 输入：原始灰度图像，光斑二值蒙版
% 输出：中心坐标
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>0tmp_A(k) = pSrc*log(pSrc);endtmp_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)];
endend

重心法

% 输入：原始灰度图像，光斑二值蒙版
% 输出：中心坐标
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];
endend

中心法

% 输入：光斑二值蒙版
% 输出：中心坐标
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)];
endend

测试代码

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

第一，定位中心的前提是分割，调用算法之前，请先完成分割。不要直接抄我的“获取蒙版”这一步，除非我们图一模一样。
第二，分割完之后，先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）

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

C版本

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

核心代码

int gausscenter(float *x,                        // （out）xfloat *y,                        // （out）yfloat *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;
}

调用方法


#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）如果报错“内存炸了”，请检查自己的二值图像，大多数人的电脑是无法处理面积大于16000个像素的光斑。不过话也说回来，你的光斑真的那么大那么清晰，也用不着这个算法了，直接质心。
2）如果报错“S（1:5）巴拉巴拉”，检查你的二值图像，是不是面积太小了，小到没有了。
3）如何生成光斑图，最快捷的方式就是用PS软件自己画。
4）光斑中心定位的前提是光斑分割，如果分割这一步（也就是生成二值图）都做错了，后面一定会错的。遇到问题优先检查二值图，确定没问题还报错，再问我就好。

（MATLAB/C）高斯拟合法求光斑中心相关推荐

“剥皮”法求区域中心
"剥皮"法求区域中心本文作者徐庆荣(武汉大学) 在已出版的有关图形图像处理的书刊中,几乎都没有专门论述求区域中心的方法,对区域中心也没有明确的定义,然而它在图像处理和分析中有着 ...
“剥皮”法求区域中心 (转)
"剥皮"法求区域中心 (转)[@more@] "剥皮"法求区域中心XML:namespace prefix = o ns = "urn:schemas ...
二维高斯曲面拟合法求取光斑中心及算法的C++实现
(1)二维高斯去曲面拟合推导一个二维高斯方程可以写成如下形式: 其中,G为高斯分布的幅值,,为x,y方向上的标准差,对式(1)两边取对数,并展开平方项,整理后为: 假如参与拟合的数据点有N个,则将这 ...
二维高斯曲面拟合法求取光斑中心
(1)二维高斯去曲面拟合推导一个二维高斯方程可以写成如下形式: 其中,G为高斯分布的幅值,,为x,y方向上的标准差,对式(1)两边取对数,并展开平方项,整理后为: 假如参与拟合的数据点有N个,则将这 ...
matlab 半高斯拟合,高斯曲线拟合求半宽高
请问下我想用matlab求拉曼光谱的拟合曲线的半宽高,用origin的高斯拟合可以直接求出一系列的数值,如下图: QQ图片20180626113405.png (83.35 KB, 下载次数: 1) ...
matlab 计算变异系数,变异系数法求权重matlab 代码
利用matlab编程,很容易根据变异系数法,求得多指标综合评价模型的权重. 代码如果有不懂的地方,可以联系我. 变异系数法求权重matlab 代码 clear;clc; [data1,header1] ...
matlab 计算变异系数,变异系数法求权重matlab代码
<变异系数法求权重matlab代码>由会员分享,可在线阅读,更多相关<变异系数法求权重matlab代码(1页珍藏版)>请在读根文库上搜索. 1.变异系数法求权重 matlab ...
用matlab编程节点电压法求电路,MATLAB在电路中的应用
<MATLAB在电路中的应用>由会员分享,可在线阅读,更多相关<MATLAB在电路中的应用(59页珍藏版)>请在人人文库网上搜索. 1.MATLAB应用(三) Matlab在电 ...
matlab e 精确到,matlab中用0.618法求minf(x)=e^(-x)+x^2在区间(0,1)上的极小值,精确到0.03....
共回答了15个问题采纳率:86.7% clc clear all; elp=0.03; tao=0.618; N=fix(log(elp)/log(tao))+1; k=1; a(k)=0; b(k) ...
matlab怎么输入积分公式,在matlab下用梯形法求函数的积分
函数是网上看到的,加上了一点我自己的理解 %用来就数值积分 %fx是由syms定义的函数表达式 function Trapezia(a,b,fx,E,Nfprintf('\n***********st ...

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

高斯拟合法求光斑中心

一、基本原理

MATLAB版本

C版本

其他

（MATLAB/C）高斯拟合法求光斑中心相关推荐

最新文章

热门文章