Opencv实现的三次样条曲线(Cubic Spline)插值

2023-11-02 21:15
文章标签 实现 opencv 三次 插值 曲线 样条 cubic spline

本文主要是介绍Opencv实现的三次样条曲线(Cubic Spline)插值，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

1.样条曲线简介

样条曲线(Spline)本质是分段多项式实函数，在实数范围内有： S:[a,b]→R ，在区间 [a,b] 上包含 k 个子区间[ti−1,ti]，且有：

a=t0<t1<⋯<tk−1<tk=b(1)

对应每一段区间 i 的存在多项式： Pi:[ti−1,ti]→R，且满足于：

S(t)=P1(t) , t0≤t<t1,S(t)=P2(t) , t1≤t<t2,⋮S(t)=Pk(t) , tk−1≤t≤tk.(2)

其中， Pi(t) 多项式中最高次项的幂，视为样条的阶数或次数（Order of spline），根据子区间 [ti−1,ti] 的区间长度是否一致分为均匀（Uniform）样条和非均匀（Non-uniform）样条。

满足了公式 (2) 的多项式有很多，为了保证曲线在 S 区间内具有据够的平滑度，一条n次样条，同时应具备处处连续且可微的性质：

P(j)i(ti)=P(j)i+1(ti);(3)

其中 i=1,…,k−1;j=0,…,n−1 。

2.三次样条曲线

2.1曲线条件

按照上述的定义，给定节点：

t:z:a=t0z0<t1z1<⋯⋯<tk−1zk−1<tkzk=b(4)

三次样条曲线满足三个条件：

  1. 在每段分段区间 [ti,ti+1],i=0,1,…,k−1 上， S(t)=Si(t) 都是一个三次多项式；
  2. 满足 S(ti)=zi,i=1,…,k−1 ;
  3. S(t) 的一阶导函数 S′(t) 和二阶导函数 S′′(t) 在区间 [a,b] 上都是连续的，从而曲线具有光滑性。

则三次样条的方程可以写为：

Si(t)=ai+bi(t−ti)+ci(t−ti)2+di(t−ti)3,(5)

其中， ai,bi,ci,di 分别代表 n 个未知系数。

  • 曲线的连续性表示为：

Si(ti)=zi,(6)

Si(ti+1)=zi+1,(7)

其中 i=0,1,…,k−1 。

  • 曲线微分连续性：

S′i(ti+1)=S′i+1(ti+1),(8)

S′′i(ti+1)=S′′i+1(ti+1),(9)

其中 i=0,1,…,k−2 。

  • 曲线的导函数表达式：

S′i=bi+2ci(t−ti)+3di(t−ti)2,(10)

S′′i(x)=2ci+6di(t−ti),(11)

令区间长度 hi=ti+1−ti ，则有：

  1. 由公式 (6) ，可得： ai=zi ；

  2. 由公式 (7) ，可得： ai+bihi+cih2i+dih3i=zi+1 ；

  3. 由公式 (8) ，可得：
    S′i(ti+1)=bi+2cihi+3dih2i ;
    S′i+1(ti+1)=bi+1 ；
    ⇒bi+2cihi+3dih2i−bi+1=0 ；

  4. 由公式 (9) ，可得：
    S′′i(ti+1)=2ci+6dihi ；
    S′′i+1(ti+1)=2ci+1 ；
    ⇒2ci+6dihi=2ci+1 ；

    设 mi=S′′i(xi)=2ci ，则：

    A. mi+6dihi−mi+1=0⇒
    di=mi+1−mi6hi ；

    B.将 ci,di 代入 zi+bihi+cih2i+dih3i=zi+1⇒
    bi=zi+1−zihi−hi2mi−hi6(mi+1−mi) ；

    C.将 bi,ci,di 代入 bi+2cihi+3dih2i=bi+1⇒

    himi+2(hi+hi+1)mi+1+hi+1mi+2=6[zi+2−zi+1hi+1−zi+1−zihi].(12)

2.2端点条件

在上述分析中，曲线段的两个端点 t0 和 tk 是不适用的，有一些常用的端点限制条件，这里只讲解自然边界。
在自然边界下，首尾两端的二阶导函数满足 S′′=0 ，即 m0=0 和 mk=0 。

3.三次样条插值类的实现

头文件
/**Cubic spline interpolation class.*
*/#ifndef CUBICSPLINEINTERPOLATION_H
#pragma once
#define CUBICSPLINEINTERPOLATION_H#include <iostream>
#include <vector>
#include <math.h>#include <opencv2/opencv.hpp>using namespace std;
using namespace cv;/* Cubic spline interpolation coefficients */
class CubicSplineCoeffs
{
public:CubicSplineCoeffs( const int &count ){a = std::vector<double>(count);b = std::vector<double>(count);c = std::vector<double>(count);d = std::vector<double>(count);}~CubicSplineCoeffs(){std::vector<double>().swap(a);std::vector<double>().swap(b);std::vector<double>().swap(c);std::vector<double>().swap(d);}public:std::vector<double> a, b, c, d;
};enum CubicSplineMode
{CUBIC_NATURAL,    // NaturalCUBIC_CLAMPED,    // TODO: ClampedCUBIC_NOT_A_KNOT  // TODO: Not a knot
};enum SplineFilterMode
{CUBIC_WITHOUT_FILTER, // without filterCUBIC_MEDIAN_FILTER  // median filter
};/* Cubic spline interpolation */
class CubicSplineInterpolation
{
public:CubicSplineInterpolation() {}~CubicSplineInterpolation() {}public:/*Calculate cubic spline coefficients.- node list x (input_x);- node list y (input_y);- output coefficients (cubicCoeffs);- ends mode (splineMode).*/void calCubicSplineCoeffs( std::vector<double> &input_x,std::vector<double> &input_y, CubicSplineCoeffs *&cubicCoeffs,CubicSplineMode splineMode = CUBIC_NATURAL,SplineFilterMode filterMode = CUBIC_MEDIAN_FILTER );/*Cubic spline interpolation for a list.- input coefficients (cubicCoeffs);- input node list x (input_x);- output node list x (output_x);- output node list y (output_y);- interpolation step (interStep).*/void cubicSplineInterpolation( CubicSplineCoeffs *&cubicCoeffs,std::vector<double> &input_x, std::vector<double> &output_x,std::vector<double> &output_y, const double interStep = 0.5 );/*Cubic spline interpolation for a value.- input coefficients (cubicCoeffs);- input a value(x);- output interpolation value(y);*/void cubicSplineInterpolation2( CubicSplineCoeffs *&cubicCoeffs,std::vector<double> input_x, double x, double &y );/*calculate  tridiagonal matrices with Thomas Algorithm(TDMA) :example:| b1 c1 0  0  0  0  |  |x1 |   |d1 || a2 b2 c2 0  0  0  |  |x2 |   |d2 || 0  a3 b3 c3 0  0  |  |x3 | = |d3 || ...         ...   |  |...|   |...|| 0  0  0  0  an bn |  |xn |   |dn |Ci = ci/bi , i=1; ci / (bi - Ci-1 * ai) , i = 2, 3, ... n-1;Di = di/bi , i=1; ( di - Di-1 * ai )/(bi - Ci-1 * ai) , i = 2, 3, ..., n-1xi = Di - Ci*xi+1 , i = n-1, n-2, 1;*/bool caltridiagonalMatrices( cv::Mat_<double> &input_a,cv::Mat_<double> &input_b, cv::Mat_<double> &input_c,cv::Mat_<double> &input_d, cv::Mat_<double> &output_x );/* Calculate the curve index interpolation belongs to */int calInterpolationIndex( double &pt, std::vector<double> &input_x );/* median filtering */void cubicMedianFilter( std::vector<double> &input, const int filterSize = 5 );double cubicSort( std::vector<double> &input );// double cubicNearestValue( std::vector );
};#endif // CUBICSPLINEINTERPOLATION_H
实现文件（cpp）
/** CubicSplineInterpolation.cpp*/#include "cubicsplineinterpolation.h"void CubicSplineInterpolation::calCubicSplineCoeffs(std::vector<double> &input_x,std::vector<double> &input_y,CubicSplineCoeffs *&cubicCoeffs,CubicSplineMode splineMode /* = CUBIC_NATURAL */,SplineFilterMode filterMode /*= CUBIC_MEDIAN_FILTER*/ )
{int sizeOfx = input_x.size();int sizeOfy = input_y.size();if ( sizeOfx != sizeOfy ){std::cout << "Data input error!" << std::endl <<"Location: CubicSplineInterpolation.cpp" <<" -> calCubicSplineCoeffs()" << std::endl;return;}/*hi*mi + 2*(hi + hi+1)*mi+1 + hi+1*mi+2=  6{ (yi+2 - yi+1)/hi+1 - (yi+1 - yi)/hi }so, ignore the both ends:| -     -     -        0           ...             0     |  |m0 || h0 2(h0+h1) h1       0           ...             0     |  |m1 || 0     h1    2(h1+h2) h2 0        ...                   |  |m2 ||         ...                      ...             0     |  |...|| 0       ...           0 h(n-2) 2(h(n-2)+h(n-1)) h(n-1) |  |   || 0       ...                      ...             -     |  |mn |*/std::vector<double> copy_y = input_y;if ( filterMode == CUBIC_MEDIAN_FILTER ){cubicMedianFilter(copy_y, 5);}const int count  = sizeOfx;const int count1 = sizeOfx - 1;const int count2 = sizeOfx - 2;const int count3 = sizeOfx - 3;cubicCoeffs = new CubicSplineCoeffs( count1 );std::vector<double> step_h( count1, 0.0 );// for m matrixcv::Mat_<double> m_a(1, count2, 0.0);cv::Mat_<double> m_b(1, count2, 0.0);cv::Mat_<double> m_c(1, count2, 0.0);cv::Mat_<double> m_d(1, count2, 0.0);cv::Mat_<double> m_part(1, count2, 0.0);cv::Mat_<double> m_all(1, count, 0.0);// initial step hifor ( int idx=0; idx < count1; idx ++ ){step_h[idx] = input_x[idx+1] - input_x[idx];}// initial coefficientsfor ( int idx=0; idx < count3; idx ++ ){m_a(idx) = step_h[idx];m_b(idx) = 2 * (step_h[idx] + step_h[idx+1]);m_c(idx) = step_h[idx+1];}// initial dfor ( int idx =0; idx < count3; idx ++ ){m_d(idx) = 6 * ((copy_y[idx+2] - copy_y[idx+1]) / step_h[idx+1] -(copy_y[idx+1] - copy_y[idx]) / step_h[idx] );}//cv::Mat_<double> matOfm( count2,  )bool isSucceed = caltridiagonalMatrices(m_a, m_b, m_c, m_d, m_part);if ( !isSucceed ){std::cout<<"Calculate tridiagonal matrices failed!"<<std::endl<<"Location: CubicSplineInterpolation.cpp -> " <<"caltridiagonalMatrices()"<<std::endl;return;}if ( splineMode == CUBIC_NATURAL ){m_all(0)      = 0.0;m_all(count1) = 0.0;for ( int i=1; i<count1; i++ ){m_all(i) = m_part(i-1);}for ( int i=0; i<count1; i++ ){cubicCoeffs->a[i] = copy_y[i];cubicCoeffs->b[i] = ( copy_y[i+1] - copy_y[i] ) / step_h[i] -step_h[i]*( 2*m_all(i) + m_all(i+1) ) / 6;cubicCoeffs->c[i] = m_all(i) / 2.0;cubicCoeffs->d[i] = ( m_all(i+1) - m_all(i) ) / ( 6.0 * step_h[i] );}}else{std::cout<<"Not define the interpolation mode!"<<std::endl;}
}void CubicSplineInterpolation::cubicSplineInterpolation(CubicSplineCoeffs *&cubicCoeffs,std::vector<double> &input_x,std::vector<double> &output_x,std::vector<double> &output_y,const double interStep )
{const int count = input_x.size();double low  = input_x[0];double high = input_x[count-1];double interBegin = low;for ( ; interBegin < high; interBegin += interStep ){int index = calInterpolationIndex(interBegin, input_x);if ( index >= 0 ){double dertx = interBegin - input_x[index];double y = cubicCoeffs->a[index] + cubicCoeffs->b[index] * dertx +cubicCoeffs->c[index] * dertx * dertx +cubicCoeffs->d[index] * dertx * dertx * dertx;output_x.push_back(interBegin);output_y.push_back(y);}}
}void CubicSplineInterpolation::cubicSplineInterpolation2(CubicSplineCoeffs *&cubicCoeffs,std::vector<double> input_x, double x, double &y)
{const int count = input_x.size();double low  = input_x[0];double high = input_x[count-1];if ( x<low || x>high ){std::cout<<"The interpolation value is out of range!"<<std::endl;}else{int index = calInterpolationIndex(x, input_x);if ( index >= 0 ){double dertx = x - input_x[index];y = cubicCoeffs->a[index] + cubicCoeffs->b[index] * dertx +cubicCoeffs->c[index] * dertx * dertx +cubicCoeffs->d[index] * dertx * dertx * dertx;}else{std::cout<<"Can't find the interpolation range!"<<std::endl;}}
}bool CubicSplineInterpolation::caltridiagonalMatrices(cv::Mat_<double> &input_a,cv::Mat_<double> &input_b,cv::Mat_<double> &input_c,cv::Mat_<double> &input_d,cv::Mat_<double> &output_x )
{int rows = input_a.rows;int cols = input_a.cols;if ( ( rows == 1 && cols > rows ) ||(cols == 1 && rows > cols ) ){const int count = ( rows > cols ? rows : cols ) - 1;output_x = cv::Mat_<double>::zeros(rows, cols);cv::Mat_<double> cCopy, dCopy;input_c.copyTo(cCopy);input_d.copyTo(dCopy);if ( input_b(0) != 0 ){cCopy(0) /= input_b(0);dCopy(0) /= input_b(0);}else{return false;}for ( int i=1; i < count; i++ ){double temp = input_b(i) - input_a(i) * cCopy(i-1);if ( temp == 0.0 ){return false;}cCopy(i) /= temp;dCopy(i) = ( dCopy(i) - dCopy(i-1)*input_a(i) ) / temp;}output_x(count) = dCopy(count);for ( int i=count-2; i > 0; i-- ){output_x(i) = dCopy(i) - cCopy(i)*output_x(i+1);}return true;}else{return false;}
}int CubicSplineInterpolation::calInterpolationIndex(double &pt, std::vector<double> &input_x )
{const int count = input_x.size()-1;int index = -1;for ( int i=0; i<count; i++ ){if ( pt > input_x[i] && pt <= input_x[i+1] ){index = i;return index;}}return index;
}void CubicSplineInterpolation::cubicMedianFilter(std::vector<double> &input, const int filterSize /* = 5 */ )
{const int count = input.size();for ( int i=filterSize/2; i<count-filterSize/2; i++ ){std::vector<double> temp(filterSize, 0.0);for ( int j=0; j<filterSize; j++ ){temp[j] = input[i+j - filterSize/2];}input[i] = cubicSort(temp);std::vector<double>().swap(temp);}for ( int i=0; i<filterSize/2; i++ ){std::vector<double> temp(filterSize, 0.0);for ( int j=0; j<filterSize; j++ ){temp[j] = input[j];}input[i] = cubicSort(temp);std::vector<double>().swap(temp);}for ( int i=count-filterSize/2; i<count; i++ ){std::vector<double> temp(filterSize, 0.0);for ( int j=0; j<filterSize; j++ ){temp[j] = input[j];}input[i] = cubicSort(temp);std::vector<double>().swap(temp);}
}double CubicSplineInterpolation::cubicSort( std::vector<double> &input )
{int iCount = input.size();for ( int j=0; j<iCount-1; j++ ){for ( int k=iCount-1; k>j; k-- ){if ( input[k-1] > input[k] ){double tp  = input[k];input[k]   = input[k-1];input[k-1] = tp;}}}return input[iCount/2];
}

这篇关于Opencv实现的三次样条曲线(Cubic Spline)插值的文章就介绍到这儿，希望我们推荐的文章对编程师们有所帮助！


http://www.chinasem.cn/article/333357。

相关文章

Java中使用Java Mail实现邮件服务功能示例

Java中使用Java Mail实现邮件服务功能示例

《Java中使用JavaMail实现邮件服务功能示例》：本文主要介绍Java中使用JavaMail实现邮件服务功能的相关资料,文章还提供了一个发送邮件的示例代码,包括创建参数类、邮件类和执行结... 目录前言一、历史背景二编程、pom依赖三、API说明（一）Session （会话）（二）Message编程客
阅读更多...
Java中List转Map的几种具体实现方式和特点

Java中List转Map的几种具体实现方式和特点

《Java中List转Map的几种具体实现方式和特点》：本文主要介绍几种常用的List转Map的方式,包括使用for循环遍历、Java8StreamAPI、ApacheCommonsCollect... 目录前言1、使用for循环遍历：2、Java8 Stream API：3、Apache Commons
阅读更多...
C#提取PDF表单数据的实现流程

C#提取PDF表单数据的实现流程

《C#提取PDF表单数据的实现流程》PDF表单是一种常见的数据收集工具,广泛应用于调查问卷、业务合同等场景,凭借出色的跨平台兼容性和标准化特点,PDF表单在各行各业中得到了广泛应用,本文将探讨如何使用... 目录引言使用工具C# 提取多个PDF表单域的数据C# 提取特定PDF表单域的数据引言PDF表单是一
阅读更多...
使用Python实现高效的端口扫描器

使用Python实现高效的端口扫描器

《使用Python实现高效的端口扫描器》在网络安全领域,端口扫描是一项基本而重要的技能,通过端口扫描,可以发现目标主机上开放的服务和端口,这对于安全评估、渗透测试等有着不可忽视的作用,本文将介绍如何使... 目录1. 端口扫描的基本原理2. 使用python实现端口扫描2.1 安装必要的库2.2 编写端口扫
阅读更多...
PyCharm接入DeepSeek实现AI编程的操作流程

PyCharm接入DeepSeek实现AI编程的操作流程

《PyCharm接入DeepSeek实现AI编程的操作流程》DeepSeek是一家专注于人工智能技术研发的公司,致力于开发高性能、低成本的AI模型,接下来,我们把DeepSeek接入到PyCharm中... 目录引言效果演示创建API key在PyCharm中下载Continue插件配置Continue引言
阅读更多...
MySQL分表自动化创建的实现方案

MySQL分表自动化创建的实现方案

《MySQL分表自动化创建的实现方案》在数据库应用场景中,随着数据量的不断增长,单表存储数据可能会面临性能瓶颈,例如查询、插入、更新等操作的效率会逐渐降低,分表是一种有效的优化策略,它将数据分散存储在... 目录一、项目目的二、实现过程（一）mysql 事件调度器结合存储过程方式1. 开启事件调度器2. 创
阅读更多...
使用Python实现操作mongodb详解

使用Python实现操作mongodb详解

《使用Python实现操作mongodb详解》这篇文章主要为大家详细介绍了使用Python实现操作mongodb的相关知识,文中的示例代码讲解详细,感兴趣的小伙伴可以跟随小编一起学习一下... 目录一、示例二、常用指令三、遇到的问题一、示例from pymongo import MongoClientf
阅读更多...
SQL Server使用SELECT INTO实现表备份的代码示例

SQL Server使用SELECT INTO实现表备份的代码示例

《SQLServer使用SELECTINTO实现表备份的代码示例》在数据库管理过程中,有时我们需要对表进行备份,以防数据丢失或修改错误,在SQLServer中,可以使用SELECTINT... 在数据库管理过程中，有时我们需要对表进行备份，以防数据丢失或修改错误。在 SQL Server 中，可以使用 SE
阅读更多...
基于Go语言实现一个压测工具

基于Go语言实现一个压测工具

《基于Go语言实现一个压测工具》这篇文章主要为大家详细介绍了基于Go语言实现一个简单的压测工具,文中的示例代码讲解详细,感兴趣的小伙伴可以跟随小编一起学习一下... 目录整体架构通用数据处理模块Http请求响应数据处理Curl参数解析处理客户端模块Http客户端处理Grpc客户端处理Websocket客户端
阅读更多...
Java CompletableFuture如何实现超时功能

Java CompletableFuture如何实现超时功能

《JavaCompletableFuture如何实现超时功能》：本文主要介绍实现超时功能的基本思路以及CompletableFuture（之后简称CF）是如何通过代码实现超时功能的,需要的... 目录基本思路CompletableFuture 的实现1. 基本实现流程2. 静态条件分析3. 内存泄露 bug
阅读更多...

Copyright China编程 23002807@qq.com

沪ICP备2023033072号-2