光斑图、阵列图、灰度图圆形等目标中心定位方法。分享高斯拟合法和更为简单的中心、重心法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坐标。
- 实验结果
可以看出对于完整的光斑,三种方法结果几乎一致。
但对于边缘处光斑(不完整),高斯拟合法的优势就体现出来了。
C版本
主要利用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 <= 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)光斑中心定位的前提是光斑分割,如果分割这一步(也就是生成二值图)都做错了,后面一定会错的。遇到问题优先检查二值图,确定没问题还报错,再问我就好。