（MATLAB/C）高斯拟合法求光斑中心 |
您所在的位置：网站首页 › 曲面拟合法 › （MATLAB/C）高斯拟合法求光斑中心 |
（MATLAB/C）高斯拟合法求光斑中心
|
高斯拟合法求光斑中心
一、基本原理MATLAB版本C版本其他
光斑图、阵列图、灰度图圆形等目标中心定位方法。分享高斯拟合法和更为简单的中心、重心法MATLAB代码，以及基于Eigen库的高斯拟合法C代码。互助互助 by HPC_ZY 一、基本原理大多数光斑其明暗分布情况都是中心最亮，往四周慢慢变暗，就类似二维高斯模型（如下图）。所以我们利用二维高斯模型去拟合光斑，从而得到光斑中心等参数。 由于本人不是数学大佬，就不推导数学公式了，直接上代码。 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>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 重心法 % 输入：原始灰度图像，光斑二值蒙版 % 输出：中心坐标 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 中心法 % 输入：光斑二值蒙版 % 输出：中心坐标 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 测试代码先别抄代码，听我说一句 先别抄代码，听我说一句 先别抄代码，听我说一句 第一，定位中心的前提是分割，调用算法之前，请先完成分割。 不要直接抄我的“获取蒙版”这一步，除非我们图一模一样。 第二，分割完之后，先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坐标。 实验结果 可以看出对于完整的光斑，三种方法结果几乎一致。 但对于边缘处光斑（不完整），高斯拟合法的优势就体现出来了。
主要利用Eigen的矩阵运算库 （Eigen下载） 核心代码 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 for (int j = 0; j xtmp[pn] = i; ytmp[pn] = j; pn++; } } } // 统计完毕 // 存入矩阵 MatrixXf X(pn, 1); MatrixXf Y(pn, 1); for (int i = 0; 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 qr; qr.compute(B); MatrixXf R = qr.matrixQR().triangularView(); MatrixXf Q = qr.householderQ(); // 特征解 MatrixXf S; S = Q.transpose()*A; MatrixXf S0(5, 1); MatrixXf R0(5, 5); for (int i = 0; i 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 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）光斑中心定位的前提是光斑分割，如果分割这一步（也就是生成二值图）都做错了，后面一定会错的。遇到问题优先检查二值图，确定没问题还报错，再问我就好。 |
【本文地址】
今日新闻 |
推荐新闻 |