【一种用opencv实现高斯曲线拟合的方法】

背景:

项目中需要实现数据的高斯拟合,进而提取数据中标准差,手头只有opencv库,经过资料查找验证,总结该方法。

基础知识:

1、opencv中solve可以实现对矩阵参数的求解;
2、线的拟合就是对多项式参数求解的过程,多项式可表示为矩阵形式;
3、高斯公式中的指数幂,可以通过取对数的方式转变成多项式的形式;
求解思路:
高斯公式->多项式公式->矩阵参数->调用solve求解;

实现过程及代码

1、确定所选的高斯公式形式

G(x)=a*exp(-((x-b)/c)^2);

2、对于给定的输入x1 ~ xn,有对输出y1 ~ yn。可以形成如下等式:

高斯公式及等式组

对等式左右两边取对数,并进行变换,可形成如下形式
对等式左右取自然对数

在这里插入图片描述
这里,就形成了AX^2+BX+C=Y的形式,其中
在这里插入图片描述
用A,B,C替换后后,原等式可写作
在这里插入图片描述
此时,我们只需要计算出A,B,C的值,再通过ABC与abc的关系即可得到abc的值。(请读者自行推导abc的公式,或见代码部分)

得到如上的多项式的形式后,直接构造参数矩阵,调用cv::solve(X,Y,A‘)接口,即可得到参数矩阵A’,其中即含有A,B,C的值。

上代码:

基础定义:

typedef struct StructMultinomialParamt
{
    double dB0;//多项式拟合的参数,数字表示幂次
    double dB1;
    double dB2;
}S_MULTNMNL_PARAMT;
typedef struct StructGaussParamT
{
    double dA;//指定的高斯参数
    double dB;//中心点
    double dC;//标准差
}S_GAUS_PARAMT;
void Gauss(S_GAUS_PARAMT sGsParamm, cv::Mat mX, cv::Mat& mY)
{
    cv::Mat mRslt = Mat::zeros(mX.size(), mX.type());
    double dx = 0;
    for (double i = 0.; i < mX.cols; i++)
    {
        for (double j = 0.; j < mX.rows; j++)
        {
            dx = mX.at<double>(j, i);
            mRslt.at<double>(j, i) = sGsParamm.dA * exp(-(pow((dx - sGsParamm.dB) / sGsParamm.dC, 2)));
        }
    }
    mY = mRslt;
    return;
}

高斯参数求解函数

void GaussFitT(cv::Mat mX, cv::Mat mY, S_GAUS_PARAMT* psGsParamm)
{
    //step1 构造参数矩阵mx与my
    cv::Mat X = Mat::zeros(mX.rows, 3, CV_64FC1);
    for (size_t i = 0; i < mX.rows; i++)
    {
        for (size_t J = 0; J < 3; J++)
        {
            X.at<double>(i, J) = pow(mX.at<double>(i, 0), 2 - J);
        }
    }
    cv::log(mY, mY);//对结果取对数
    //step2 多项式拟合
    cv::Mat A;//参数矩阵
    cv::solve(X, mY, A, cv::DECOMP_SVD);
    S_MULTNMNL_PARAMT sBparam;
    sBparam.dB2 = A.at<double>(0);
    sBparam.dB1 = A.at<double>(1);
    sBparam.dB0 = A.at<double>(2);

    //step3 高斯参数计算ABC-》abc
    psGsParamm->dA = exp(sBparam.dB0 - pow(sBparam.dB1, 2) / (4 * sBparam.dB2));
    psGsParamm->dB = -sBparam.dB1 / (2 * sBparam.dB2);
    psGsParamm->dC = sqrt(-1 / sBparam.dB2);

    return;
}

# 测试代码

double dX[50];//输入数据X
double dY[50];//输入数据Y
std::vector<cv::Point> pointsOri;

for (int i = 0; i < 50; i++)
{

    dX[i] = double(i);
    dY[i] = -0.5 * pow((dX[i] - 25), 2) + 320 + i;
    pointsOri.push_back(cv::Point(dX[i], dY[i]));
}
//转换成求解函数输入需要的数据格式
cv::Mat mGsInputX = Mat::zeros(50, 1, CV_64FC1);
cv::Mat mGsInputY = Mat::zeros(50, 1, CV_64FC1);
for (size_t i = 0; i < 50; i++)
{
    mGsInputX.at<double>(i) = dX[i];
    mGsInputY.at<double>(i) = dY[i];
}

S_GAUS_PARAMT sGsParamm;//求解结果
GaussFitT(mGsInputX, mGsInputY, &sGsParamm);

//结果对比
Mat mGsOutputY;
Gauss(sGsParamm, mGsInputX, mGsOutputY);
std::vector<cv::Point> pointsNew;//拟合结果
for (int i = 0; i < 50; i++)
{
    pointsNew.push_back(cv::Point(dX[i], mGsOutputY.at<double>(i)));
}
cv::Mat img(450, 60, CV_8UC3, cv::Scalar(0, 0, 0));
cv::polylines(img, std::vector<std::vector<cv::Point>>{pointsOri}, false, cv::Scalar(0, 0, 255), 2);
cv::polylines(img, std::vector<std::vector<cv::Point>>{pointsNew}, false, cv::Scalar(255, 255, 255), 0.5);

// 显示图像
cv::imshow("Line Chart", img);
cv::waitKey(0);

运行输出

在这里插入图片描述

红色的为原始数据分布,白色的为拟合计算结果。
而我需要的标准差,则为sGsParamm.dC。

参考:https://blog.csdn.net/guangjie2333/article/details/115629152
https://blog.csdn.net/KYJL888/article/details/103073956
https://blog.csdn.net/qq_35097289/article/details/103910984

后记:

调用solve的接口求解时,OPENCV提供了以下六种方式以对应不同的情况。对于多项式的求解,也可以采用最小二乘法的逼近,不再调用solve方法,这块后面再填坑吧。

cv::DECOMP_LU 高斯消元法(LU分解)
cv::DECOMP_SVD 奇异值分解(SVD)
cv::DECOMP_CHOLESKY 对于对称正定矩阵
cv::DECOMP_EIG 特征值分解,只用于对称矩阵
cv::DECOMP_QR QR因式分解
cv::DECOMP_NORMAL 可选附加标志,表示要求解标准方程

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:/a/252829.html

如若内容造成侵权/违法违规/事实不符,请联系我们进行投诉反馈qq邮箱809451989@qq.com,一经查实,立即删除!

相关文章

【深度强化学习】确定性策略梯度算法 DDPG

前面讲到如 REINFORCE&#xff0c;Actor-Critic&#xff0c;TRPO&#xff0c;PPO 等算法&#xff0c;它们都是随机性策略梯度算法&#xff08;Stochastic policy&#xff09;&#xff0c;在广泛的任务上表现良好&#xff0c;因为这类方法鼓励了算法探索&#xff0c;给出的策略是…

探索 Vim:一个强大的文本编辑器

引言&#xff1a; Vim&#xff08;Vi IMproved&#xff09;是一款备受推崇的文本编辑器&#xff0c;拥有强大的功能和高度可定制性&#xff0c;提供丰富的编辑和编程体验。本文将探讨 Vim 的基本概念、使用技巧以及为用户带来的独特优势。 简介和发展 1. Vim 的简介和历史 V…

【二分查找】自写二分函数的总结

作者推荐 【动态规划】【广度优先搜索】LeetCode:2617 网格图中最少访问的格子数 本文涉及的基础知识点 二分查找算法合集 自写二分函数 的封装 我暂时只发现两种&#xff1a; 一&#xff0c;在左闭右开的区间寻找最后一个符合条件的元素&#xff0c;我封装成FindEnd函数。…

力扣刷题-二叉树-平衡二叉树

110 平衡二叉树 给定一个二叉树&#xff0c;判断它是否是高度平衡的二叉树。 本题中&#xff0c;一棵高度平衡二叉树定义为&#xff1a;一个二叉树每个节点 的左右两个子树的高度差的绝对值不超过1。 示例 1: 给定二叉树 [3,9,20,null,null,15,7] 返回 true 。 给定二叉树 [1…

AUTOSAR ComM模块配置以及代码

ComM模块配置以及代码执行流程 1、基本的一个通道的配置列表 ComMNmVariant 概念的个人理解&#xff1a; FULL&#xff1a; 完全按照AUTOSAR NM方式进行调用 LIGHT &#xff1a;设置一个超时时间&#xff0c;在请求停止通信的时候开始计时&#xff0c;超时之后才会进入FULLCOM…

运维实践|采集MySQL数据出现many connection errors

文章目录 问题出现问题分析当前环境问题分析 解决方案1 检查调度事件任务是否开启2 开启调度事件任务3 创建一张日志表4 创建函数存储过程5 创建事件定时器6 开启事件调度任务7 检查核实是否创建 总结 问题出现 最近在做OGG结构化数据采集工作&#xff0c;在数据采集过程中&am…

将博客搬至微信公众号了

一、博客搬家通知 各位码友们好&#xff0c;大家是不是基本很少看 CSDN 了呢&#xff1f;平时开发是不都依靠着 chatGPT 来解决工作中的技术问题了&#xff0c;不过我觉得在工作中的使用场景各式各样的&#xff0c;具体问题还得自己具体去梳理逻辑&#xff0c;再考虑用什么样的…

C语言:求和1+1/2-1/3+1/4-1/5+……-1/99+1/100

#include<stdio.h> int main() {int i 0;double sum 0.0;int flag 1;for (i 1;i < 100;i){sum 1.0 / i * flag;flag -flag;}printf("sum%lf\n", sum);return 0; }

SpringIOC之@Primary

博主介绍&#xff1a;✌全网粉丝5W&#xff0c;全栈开发工程师&#xff0c;从事多年软件开发&#xff0c;在大厂呆过。持有软件中级、六级等证书。可提供微服务项目搭建与毕业项目实战&#xff0c;博主也曾写过优秀论文&#xff0c;查重率极低&#xff0c;在这方面有丰富的经验…

力扣刷题-二叉树-找树左下角的值

513 找树左下角的值 给定一个二叉树的 根节点 root&#xff0c;请找出该二叉树的 最底层 最左边 节点的值。 假设二叉树中至少有一个节点。 示例 1&#xff1a; 示例 2&#xff1a; 思路 层序遍历 直接层序遍历&#xff0c;因为题目说了是最底层&#xff0c;最左边的值&a…

【数据结构—队列的实现】

提示&#xff1a;文章写完后&#xff0c;目录可以自动生成&#xff0c;如何生成可参考右边的帮助文档 文章目录 前言 一、队列 1.1队列的概念及结构 二、队列的实现 2.1头文件的实现—Queue.h 2.2源文件的实现—Queue.c 2.3源文件的测试—test.c 三、测试队列实际数据的展示 3.…

mysql使用st_distance_sphere函数报错Incorrect arguments to st_distance_sphere

前言 最近使用空间点位查询数据时函数报错Incorrect arguments to st_distance_sphere报错。 发现问题 因为之前是没有问题的&#xff0c;所以把问题指向了数据&#xff0c;因为是外部数据&#xff0c;不是通过系统打点获取&#xff0c;发现是因为经纬度反了&#xff0c;loc…

VRRP(虚拟路由冗余协议)

一.VRRP简介 1.VRRP是什么 Virtual route Redundancy Protocol&#xff0c;也叫虚拟路由器冗余协议。 利用VRRP&#xff0c;一组路由器协同工作&#xff0c;单只有一个处于Master状态&#xff0c;处于该状态的路由器&#xff08;的接口&#xff09;承担实际的数据流量转发任…

MapReduce序列化实例代码

1 &#xff09;需求&#xff1a;统计每个学号该月的超市消费、食堂消费、总消费 2 &#xff09;输入数据格式 序号 学号 超市消费 食堂消费 18 202200153105 8.78 12 3 &#xff09;期望输出格式 key &#xff08;学号&#xff09; value &#xff08; bean 对象&#xf…

二分查找算法的概念、原理、效率以及使用C语言循环和数组的简单实现

二分查找的概念 二分查找也称折半查找&#xff08;Binary Search&#xff09;&#xff0c;它是一种效率较高的查找方法。但是&#xff0c;折半查找要求线性表必须采用顺序存储结构&#xff0c;而且表中元素按关键字有序排列。 实现原理 首先&#xff0c;假设表中元素是按升序…

深度学习项目实战:垃圾分类系统

简介&#xff1a; 今天开启深度学习另一板块。就是计算机视觉方向&#xff0c;这里主要讨论图像分类任务–垃圾分类系统。其实这个项目早在19年的时候&#xff0c;我就写好了一个版本了。之前使用的是python搭建深度学习网络&#xff0c;然后前后端交互的采用的是java spring …

VLAN间的通讯---三层交换

一.三层交换 1.概念 使用三层交换技术实现VLAN间通信 三层交换二层交换 三层转发 2.基于CEF的MLS CEF是一种基于拓补转发的模型 转发信息库&#xff08;FIB&#xff09;临接关系表 转发信息库&#xff08;FIB&#xff09;可以理解为路由表 邻接关系表可以理解为MAC地址表…

js获取日期

目录 Date 对象 1. 获取当前时间 2. 获取特定日期时间 Date 对象的方法 1. 获取各种日期时间组件 2. 获取星期几 3. 获取时间戳 格式化日期时间 1. 使用 toLocaleString() 方法 2. 使用第三方库 UNIX 时间戳 内部表示 时区 Date 对象 JavaScript中内置的 Date 对象…

【sqli靶场】第六关和第七关通关思路

目录 前言 一、sqli靶场第六关 1.1 判断注入类型 1.2 观察报错 1.3 使用extractvalue函数报错 1.4 爆出数据库中的表名 二、sqli靶场第七关 1.1 判断注入类型 1.2 判断数据表中的字段数 1.3 提示 1.4 构造poc爆库名 1.5 构造poc爆表名 1.6 构造poc爆字段名 1.7 构造poc获取账…

MySQL 报错 You can‘t specify target table for update in FROM clause解决办法

You can’t specify target table for update in FROM clause 其含义是&#xff1a;不能在同一表中查询的数据作为同一表的更新数 单独执行复合查询是正常的&#xff0c;如下&#xff1a; 但是当执行子查询删除命令时&#xff0c;报如下错误 DELETE FROM abpusers WHERE Id I…