akima 插值拟合算法 Python/C++版本

news2024/11/25 10:34:46

目录

  • 前言
    • Akima简介
    • Akima优势
  • 算法的代码实现
    • python版
    • C++ 版
      • 代码解析1
      • 代码解析2
      • 代码解析3
    • 结果测试

前言

鉴于CSDN上Akima算法文章大部分要VIP观看或者下载,即使是付费也有质量不佳,浪费Money也浪费时间。
笔者更具查到的资料分享给大家。

Akima简介

Akima 拟合算法是 Hiroshi Akima 于 1970 年开发的一种插值和曲线拟合方法。Akima 插值算法对于构造通过给定数据点集的平滑曲线特别有用。它广泛应用于各个领域,包括计算机图形学、图像处理和数值分析。

Akima 拟合算法不同于传统的插值方法,例如线性或多项式插值,它提供了更稳健和视觉上令人愉悦的数据表示。它侧重于最大限度地减少在其他插值技术中经常观察到的振荡或摆动。此算法的工作原理是将给定数据分成小区间,然后将分段三次曲线拟合到每个区间。该方法确保生成的曲线在区间边界处是连续的并且准确地逼近数据。Akima 的方法同时考虑了数据点的斜率和曲率,从而产生更平滑和更具视觉吸引力的插值。

Akima优势

Akima 拟合的优势之一是它能够处理间隔不均匀的数据点,使其适用于不规则采样的数据。该算法还解决了数据点突然变化或不连续的情况。

算法的代码实现

鉴于很多是嵌入式上用Akima算法,这里在python版之外还提供了C++版。

python版

需要scipy包,里面直接有Akima拟合函数。
x,y是自己定义的需要Akima拟合的曲线的坐标,把这些坐标放到akima_interpolator里去。
然后将一个新的x输入akima_interpolator就能得到拟合的y了,是不是很简单。

import numpy as np
from scipy.interpolate import Akima1DInterpolator

# Generate sample data
x = np.array([1, 2, 3, 4, 5])
y = np.array([1, 3, 2, 5, 4])

# Perform Akima fitting
akima_interpolator = Akima1DInterpolator(x, y)

# Generate new x-values for interpolation
x_new = np.linspace(1, 5, num=100)

# Interpolate y-values using Akima fitting
y_new = akima_interpolator(x_new)

# Print the interpolated values
for i in range(len(x_new)):
    print(f"x: {x_new[i]}, y: {y_new[i]}")

C++ 版

前面的python版全程用akima包,细节看不到,C++没有这种包,但我们能清楚看到里面细节。

#include <iostream>
#include <vector>
#include <cmath>

// Akima插值函数
double akimaInterpolation(double x, const std::vector<double>& xData, const std::vector<double>& yData) {
    int n = xData.size();
    
    int index = 0;
    
    // Find the interval index
    for (int i = 0; i < n - 1; ++i) {
        if (x >= xData[i] && x <= xData[i + 1]) {
            index = i;
            break;
        }
    }
    
    // 计算斜率
    std::vector<double> slopes(n - 1); //初始化n-1个默认值为0的元素
    for (int i = 0; i < n - 1; ++i) {
        double dx = xData[i + 1] - xData[i];
        double dy = yData[i + 1] - yData[i];
        slopes[i] = dy / dx; //计算每段之间的斜率
    }
    
    // 计算权重
    std::vector<double> weights(n - 1); //初始化n-1个默认值为0的元素
    for (int i = 2; i < n - 2; ++i) {
        weights[i] = std::abs(slopes[i + 1] - slopes[i - 1]); //计算这些权重的目的是确定每个间隔附近的斜坡的“强度”。这些权重随后用于插值公式中,以确保插值曲线的平滑性和连续性。
    }
    
    // 计算插值
    double dx = xData[index + 1] - xData[index];
    double t = (x - xData[index]) / dx;  //参数 t 表示区间内的归一化位置,取值范围为 0 到 1
    
    //m0、m1、p0和p1是 A​​kima 插值公式中用于计算插值的系数
    double m0 = slopes[index] * dx; //详见代码解析1
    double m1 = slopes[index + 1] * dx;
    double p0 = (3 * weights[index] - 2 * m0 - m1) / dx; //这里的3,2系数是怎么来的详见代码解析2
    double p1 = (3 * weights[index + 1] - m0 - 2 * m1) / dx;
    //interpolatedValue 这个公式用于计算最终插值结果,详见代码解析3
    double interpolatedValue =
        yData[index] * (1 - t) * (1 - t) * (1 + 2 * t) +
        yData[index + 1] * t * t * (3 - 2 * t) +
        p0 * t * (1 - t) * (1 - t) +
        p1 * t * t * (t - 1);
    
    return interpolatedValue;
}

int main() {
    std::vector<double> xData = {1, 2, 3, 4, 5};
    std::vector<double> yData = {1, 3, 2, 5, 4};
    
    // 假设输入一个x=2.5,y输出多少?
    double interpolatedValue = akimaInterpolation(2.5, xData, yData);
    std::cout << "Interpolated value at x = 2.5: " << interpolatedValue << std::endl;
    
    return 0;
}

解释一下上面的斜率和权重,斜率是通过相邻点之间 k=dy/dx 来计算。而权重是区间附近斜率对这个区间影响的权重,将点i的左侧斜率slopes[i - 1]和右侧斜率slopes[i + 1]相减得到,存在weights[i]里。权重随后用于插值公式中,以确保插值曲线的平滑性和连续性。

这里展开讲一下:
在 Akima 插值中,插值曲线是通过将分段三次曲线拟合到连续数据点之间的每个区间来构建的。这些三次曲线的斜率在确定插值曲线的形状和行为方面起着至关重要的作用。目标是确保曲线连续并遵循数据的总体趋势,同时避免过度振荡。

通过计算权重算法考虑了相邻区间之间斜率的变化。权重通过捕获数据的局部行为并影响插值过程中每个斜率的“强度”。较大的权重表示区域斜率变化明显,而较小的权重表示区域较平滑。

在执行插值时,将权重合并到插值公式中以调整相邻斜率的贡献。权重充当控制不同区间斜率之间平衡的系数。此调整有助于平滑插值曲线并减少由异常值或噪声数据点引起的突然变化。

代码解析1

m0、m1、p0和p1是 A​​kima 插值公式中用于计算插值的系数:

m0:此变量表示左相邻区间的调整斜率。将当前区间 (slopes[index]) 的斜率乘以区间宽度 ( dx)得到。

m1:此变量表示右相邻区间的调整斜率。将下一个区间 (slopes[index + 1]) 的斜率乘以区间的宽度 (dx)得到。

p0:此变量表示左相邻区间的调整权重。它是使用当前区间 ( weights[index]) 的权重、左侧区间的调整斜率 ( m0) 和右侧区间的调整斜率 (m1) 计算得到。由公式(3 * weights[index] - 2 * m0 - m1) / dx确定左相邻区间对插值的贡献。

p1:此变量表示右相邻区间的调整权重。它是使用下一个区间的权重 ( weights[index + 1])、左侧区间的调整斜率 (m0) 和右侧区间的调整斜率 (m1) 计算得到。由公式(3 * weights[index + 1] - m0 - 2 * m1) / dx确定右相邻区间对插值的贡献。

代码解析2

公式(3 * weights[index] - 2 * m0 - m1) / dx 和 (3 * weights[index + 1] - m0 - 2 * m1) / dx 是基于Akima插值方案推导出来的。
为了理解推导,让我们考虑 Akima 插值方案的一般形式:

y(x) = p0(x) * y0 + p1(x) * y1 + q0(x) * m0 + q1(x) * m1

在此等式中,y(x)表示特定坐标处的插值x。y0和y1是x两侧的数据点, m0和m1是与数据点关联的斜率。项p0(x)和p1(x)是数据点的权重系数,q0(x)和q1(x)是数据点关联的斜率的权重系数。
为了确定p0(x)和p1(x),Akima 拟合使用三次多项式来确保平滑性和连续性。这些权重系数由斜率的局部行为决定。

通过考虑Akima插值方案,我们可以推导出代码中使用的具体权重公式:

对于p0(x):权重函数p0(x)决定了左邻域的贡献。在代码中,(3 * weights[index] - 2 * m0 - m1) / dx代表p0(x).
选择特定系数3、-2和1是为了平衡斜率的影响并确保间隔边界处的连续性。这些系数是通过数学分析和优化确定的。

对于p1(x):权重函数p1(x)决定了右邻区间的贡献。在代码中,(3 * weights[index + 1] - m0 - 2 * m1) / dx代表p1(x).同样,选择系数3、-1和-2以实现插值曲线的连续性和平滑性。

导出这些公式中的特定系数是为了最大限度地减少插值误差并保持曲线的连续性。它们是通过数学分析和优化技术确定的,以确保生成的曲线与基础数据点紧密匹配。

代码解析3

  • yData[index] * (1 - t) * (1 - t) * (1 + 2 * t):这一项代表左边数据点(yData [index]) 对插值的贡献。它乘以三次多项式“(1 - t) * (1 - t) * (1 + 2 * t)”,该多项式取决于参数“t”,范围从 0 到 1。多项式旨在确保左侧数据点的平滑过渡和适当加权。
  • yData[index + 1] * t * t * (3 - 2 * t):此项表示右侧数据点 (yData[index + 1]) 对插值的贡献。它乘以三次多项式“t * t * (3 - 2 * t)”。与上面类似,这个多项式确保了右侧数据点的平滑过渡和适当加权。
  • p0 * t * (1 - t) * (1 - t):此项表示左侧相邻区间的调整权重 (p0) 对插值的贡献。它乘以三次多项式“t * (1 - t) * (1 - t)”。该多项式表示左侧相邻区间对插值的影响。
  • p1 * t * t * (t - 1):此项表示右相邻区间的调整权重 (p1) 对插值的贡献。它乘以三次多项式“t * t * (t - 1)”。该多项式表示右侧邻区间对插值的影响。
    该方程结合了相邻数据点的贡献及其相应的权重来计算最终的插值。参数 t 表示区间内的归一化位置,取值范围为 0 到 1。它决定了相邻数据点及其对应区间的相对权重。应用于数据点和权重的三次多项式确保插值曲线的平滑性和连续性。

把上面这些影响因素加一起就是插值点的函数值interpolatedValue了。

结果测试

笔者使用Dev C++ 运行代码,插值在合理范围内,大家可以多找几个点试试。
在这里插入图片描述

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

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈,一经查实,立即删除!

相关文章

C++技能系列 ( 5 ) - 详解函数入参/返回参使用(值传递/引用传递/指针传递/智能指针传递)

系列文章目录 C技能系列 Linux通信架构系列 C高性能优化编程系列 深入理解软件架构设计系列 高级C并发线程编程 期待动动小手&#xff0c;点击关注哦&#xff01;&#xff01;&#xff01; 当你休息的时候&#xff0c;一定要想到别人还在奔跑。 When you rest, we must thin…

数据库相关

1、主要考点思维导图 2、如何设计一个关系型数据库 存储管理&#xff1a;数据逻辑关系转为物理存储关系。 缓存机制&#xff1a;优化执行效率。 SQL解析&#xff1a;将Sql语句进行解析。 日志管理&#xff1a;记录操作。 权限划分&#xff1a;多用户管理。 容灾机制&…

2、瑞丽-伯纳德对流的拉格朗日拟序结构(FTLE场结果对比)

在上篇博客中&#xff0c;我简单比较了瑞丽伯纳德对流的FTLE场&#xff0c;但是因为粒子追踪采用的是欧拉方法&#xff0c;所以精度不是很高&#xff0c; 因此与文献中的结果还是有些差别。 下面放一张文献中的FTLE场&#xff0c;参数与上篇文章是一致的&#xff0c;Ra 1e8;Pr…

《面试1v1》SpringBean生命周期

&#x1f345; 作者简介&#xff1a;王哥&#xff0c;CSDN2022博客总榜Top100&#x1f3c6;、博客专家&#x1f4aa; &#x1f345; 技术交流&#xff1a;定期更新Java硬核干货&#xff0c;不定期送书活动 &#x1f345; 王哥多年工作总结&#xff1a;Java学习路线总结&#xf…

python:使用Scikit-image对遥感影像进行边缘检测(edges)

作者:CSDN @ _养乐多_ 本文将介绍使用Scikit-image库中用于边缘检测特征提取的一些方法及其代码。方法包括 Canny边缘检测(canny),Sobel边缘检测(sobel),Scharr边缘检测(scharr), Roberts边缘检测(roberts),Prewitt边缘检测(prewitt),Farid边缘检测(farid),…

HJ48 从单向链表中删除指定值的节点

描述 输入一个单向链表和一个节点的值&#xff0c;从单向链表中删除等于该值的节点&#xff0c;删除后如果链表中无节点则返回空指针。 链表的值不能重复。 构造过程&#xff0c;例如输入一行数据为: 6 2 1 2 3 2 5 1 4 5 7 2 2 则第一个参数6表示输入总共6个节点&am…

Day24 实战篇 ——Jmeter通过JDBC测试实战

Day24 实战篇 ——Jmeter通过JDBC测试实战 文章目录 Day24 实战篇 ——Jmeter通过JDBC测试实战1、**业务级脚本开发**2、**接口级脚本开发**3、**JDBC脚本开发**4、**JMS Point-to-Poibt脚本开发**5、**Jmeter轻量级接口自动化测试框架**(了解就行)1、业务级脚本开发 登录脚…

《2023 年 React 生态》

大家好&#xff0c;我是 Chocolate。 前不久看到一篇不错的内容&#xff0c;来自于 The React Ecosystem in 2023&#xff0c;也结合自己今年使用的 React 生态总结一下。 本文并非视频演讲稿&#xff0c;和视频内容还是有一点点区别&#xff0c;视频内容相对来说会更加详细一…

ffmpeg(一) ffmpeg+QT开发环境搭建

1、开发库的选择 &#xff08;1&#xff09;音视频开发库 每个主流平台基本都有自己的音视频开发库&#xff08;API&#xff09;&#xff0c;用以处理音视频数据&#xff0c;比如&#xff1a; iOS&#xff1a;AVFoundation、AudioUnit 等 Android&#xff1a;MediaPlayer、Me…

Flink 学习九 Flink 程序分布式运行部署

Flink 学习九 Flink 程序分布式运行部署 1.Job 执行计划 层级说明备注StreamGraph用户代码生成的最初的图程序的运行流程图JobGraph将多个符合条件的节点多个符合条件的节点合并,减少序列化和反序列化ExecutionGraphJobGraph 的并行化调度层的核心数据结构PhysicalGraphJobMa…

【计算机组成原理】信息编码与数据表示

目录 一、进位计数制 二、信息编码 三、定点数的表示 四、校验码 五、浮点数的表示 一、进位计数制 整数部分&#xff1a; 二进制、八进制、十六进制 ---> 十进制&#xff1a;加权求和二进制 ---> 八进制&#xff1a;每三位分为一组&#xff0c;转为八进制…

CloFormer实战:使用CloFormer实现图像分类任务(一)

文章目录 摘要安装包安装efficientnet_pytorch安装timm安装 grad-cam 数据增强Cutout和MixupEMA项目结构计算mean和std生成数据集 摘要 论文翻译&#xff1a;https://blog.csdn.net/m0_47867638/article/details/131161083 官方源码&#xff1a;https://github.com/qhfan/CloF…

faceswap安装教程图文详解

Faceswap是一种人脸识别技术&#xff0c;可以将一个人的面部特征与另一个人的面部特征进行交换&#xff0c;从而创建出一个看起来像是两个人融合在一起的图像或视频。这项技术可以用于各种目的&#xff0c;包括艺术创作、电影制作、虚拟现实、安全监控等领域。Faceswap的实现方…

UE特效案例 —— 寒冰武器

一&#xff0c;环境配置 创建默认地形Landscape&#xff0c;如给地形上材质需确定比例&#xff1b;添加环境主光源DirectionalLight&#xff0c;设置相应的强度和颜色&#xff1b;PostProcessVolume设置曝光&#xff0c;设置Min/Max Brightness为1&#xff1b; 与关闭Game Sett…

怎样开始用selenium进行自动化测试?

如果您刚开始使用 Selenium 进行自动化测试&#xff0c;以下是建议的步骤。 1、安装 Selenium 首先&#xff0c;您需要安装 Selenium。Selenium 支持多种编程语言&#xff0c;如 Python、Java、C# 等。可以通过 pip 命令在 Python 中安装 Selenium&#xff1a; pip install …

CloFormer实战:使用CloFormer实现图像分类任务(二)

文章目录 训练部分导入项目使用的库设置随机因子设置全局参数图像预处理与增强读取数据设置Loss设置模型设置优化器和学习率调整算法设置混合精度&#xff0c;DP多卡&#xff0c;EMA定义训练和验证函数训练函数验证函数调用训练和验证方法 运行以及结果查看测试热力图可视化展示…

秀米编辑器(xiumi)+百度编辑器(Ueditor) 集成 :解决集成问题,秀米编辑器导出到百度编辑器格式问题,图片保存到自己的服务器(阿里云OSS)

1.集成前提条件&#xff1a; 1. 需要集成百度编辑器到环境中 2.https环境下才可以导出数据到百度编辑器&#xff0c;如果不是https环境&#xff0c;会出现错误 然后我们开始讲解如何集成&#xff1a; 2.引入资源&#xff1a; //百度编辑器需要修改的文件&#xff08;配置与原始…

测试入门第一步------编写接口测试用例

自动化始终只是辅助测试工作的一个手段&#xff0c;对于测试人员而言&#xff0c;测试基础和测试用例的设计才是核心。如果测试用例的覆盖率或者质量不高&#xff0c;那将这部分用例实现为自动化用例的意义也就不大了。 那么&#xff0c;接口测试用例应该怎么编写呢&#xff1…

Spring boot集成RabbitMq

Spring boot集成RabbitMq 一、搭建RabbitMq1.1 参考1.2 配置erlong的环境变量1.3 RabbitMQ对应的在注册表中的位置 二、使用教程2.1 打开服务端2.2 注意的问题2.3 Queue的包 三、git命令3.1 git remote3.2 git remote add origin "xxxx"3.3 git push -u origin maste…

使用esp32+micropython+microdot搭建web(http+websocket)服务器(超详细)第一部分

使用esp32micropythonmicrodot搭建web(httpwebsocket)服务器&#xff08;超详细&#xff09;第一部分 microdot文档速查 什么是Microdot?Microdot是一个可以在micropython中搭建物联网web服务器的框架micropyton文档api速查 Quick reference for the ESP32 先来个小demo先体…