专栏首页爱笑的架构师二维高斯曲面拟合法求取光斑中心及算法的C++实现

二维高斯曲面拟合法求取光斑中心及算法的C++实现

(1)二维高斯去曲面拟合推导

一个二维高斯方程可以写成如下形式:

其中,G为高斯分布的幅值,,为x,y方向上的标准差,对式(1)两边取对数,并展开平方项,整理后为:

假如参与拟合的数据点有N个,则将这个N个数据点写成矩阵的形式为:A = B C

其中:

A为N*1的向量,其元素为:

B为N*5的矩阵:

C为一个由高斯参数组成的向量:

(2)求解二维高斯曲线拟合

N个数据点误差的列向量为:E=A-BC,用最小二乘法拟合,使其N个数据点的均方差最小,即:

在图像数据处理时,数据量比较大,为减小计算量,将矩阵B进行QR分解,即:B=QR,分解后Q为一个N*N的正交矩阵,R为一个N*5的上三角矩阵,对E=A-BC进行如下推导:

由于Q为正交矩阵,可以得到:

令:

其中: S为一个5维列向量;T为一个N-5维列向量;R1为一个5*5的上三角方阵,则

上式中,当S = R1C时取得最小值,因此只需解出:

即可求出:

中的

这些参数,这里先给出:

这里:

(3)C++代码实现,算法的实现过程中由于涉及大量的矩阵运算,所以采用了第三方的开源矩阵算法Eigen,这里真正用于高斯拟合的函数是

bool GetCentrePoint(float& x0,float& y0)

关于Eigen的用法请参考:http://blog.csdn.net/hjx_1000/article/details/8490941

#include "stdafx.h"
#include "GaussAlgorithm.h"

VectorXf m_Vector_A;
MatrixXf m_matrix_B;
int m_iN =-1;

bool InitData(int pSrc[100][100], int iWidth, int iHeight)
{
	if (NULL == pSrc || iWidth <=0 || iHeight <= 0)
		return false;
	m_iN = iWidth*iHeight;
	VectorXf tmp_A(m_iN);
	MatrixXf tmp_B(m_iN, 5);
	int i =0, j=0, iPos =0;
	while(i<iWidth)
	{
		 j=0;
		while(j<iHeight)
		{
			tmp_A(iPos) = pSrc[i][j] * log((float)pSrc[i][j]);

			tmp_B(iPos,0) = pSrc[i][j] ;
			tmp_B(iPos,1) = pSrc[i][j] * i;
			tmp_B(iPos,2) = pSrc[i][j] * j;
			tmp_B(iPos,3) = pSrc[i][j] * i * i;
			tmp_B(iPos,4) = pSrc[i][j] * j * j;
			++iPos;
			++j;
		}
		++i;
	}
	m_Vector_A = tmp_A;
	m_matrix_B = tmp_B;

}
bool GetCentrePoint(float& x0,float& y0)
{
	if (m_iN<=0)
		return false;
	//QR分解
	HouseholderQR<MatrixXf> qr;
	qr.compute(m_matrix_B);
	MatrixXf R = qr.matrixQR().triangularView<Upper>();
	MatrixXf Q =  qr.householderQ();

	//块操作,获取向量或矩阵的局部
	VectorXf S;
	S = (Q.transpose()* m_Vector_A).head(5);
	MatrixXf R1;
	R1 = R.block(0,0,5,5);

	VectorXf C;
	C = R1.inverse() * S;

	x0 = -0.5 * C[1] / C[3];
	y0 = -0.5 * C[2] / C[4];

	return true;
}

其中:

函数 bool InitData(int pSrc[100][100], int iWidth, int iHeight)主要用于将数组转换成相应的矩阵,因为我的所有数据都在一个整形的二维数组中存着,所以需要做这种转换。

函数bool GetCentrePoint(float& x0,float& y0)主要用于对数据点进行二维高斯曲面拟合,并返回拟合的光点中心。

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • 2020-10-22在线识图搜索引擎

    最近在逛淘宝时发现了淘宝的图片搜索功能,可能是我太Low了这个技术点已经实现很长时间了。想想自己能不能实现这个功能,起初我是这么想的,对两张图片从左上角的第一个...

    爱笑的架构师
  • 【java设计模式系列】1. 工厂方法模式(Factory Method)

    2. 多个工厂方法模式 ,是对普通工厂方法模式的改进,在普通工厂方法模式中,如果传递的字符串出错,则不能正确创建对象,而多个工厂方法模式是提供多个工厂方法,...

    爱笑的架构师
  • 2020-10-21CUDA从入门到精通

    在老板的要求下,本博主从2012年上高性能计算课程开始接触CUDA编程,随后将该技术应用到了实际项目中,使处理程序加速超过1K,可见基于图形显示器的并行计算对于...

    爱笑的架构师
  • 恕我直言,我怀疑你并不会生成随机数

    有一次,我在逛 Stack Overflow 的时候,发现有这样一个问题:“Java 中如何产生一个指定范围内的随机数”,我心想,“就这破问题,竟然有 398 ...

    沉默王二
  • Go语言排序与接口实例分析

    本文实例讲述了Go语言排序与接口用法。分享给大家供大家参考。具体如下: import "fmt" type Sorter interface { Len()...

    李海彬
  • Java高并发之无锁与Atomic源码分析

    用户1216491
  • 小数化为分数

    题目 简单的说就是将有限循环小数和无限循环小数转化为分数形式。比如: 0.9 = 9/10 0.333(3) = 1/3,其...

    小二哥
  • Spring Boot实战-使用WebFlux进行响应式编程(2)

    以上代码使用Map完成了关于用户的增删改查操作。这也是我们在项目中运用最多的操作。如果只是想要在完成时给出完成信号,就可以使用 Mono<Void>。 接下来...

    我的小熊不见了丶
  • Java基本数据类型的包装类

    Java语言是一个面向对象的语言,但是Java中的基本数据类型却是不面向对象的,这在实际使用时存在很多的不便,为了解决这个不足,在设计类时为每个基本数据类型设计...

    用户5224393
  • LWC 109: 933.Number of Recent Calls

    版权声明:本文为博主原创文章,未经博主允许不得转载。 https://blog.csdn.n...

    用户1147447

扫码关注云+社区

领取腾讯云代金券