1、发文初心

样条曲线插值应用挺多的，其实这方面的代码不少。

不过都有一些缺憾，大家使用起来不便利，因而特意改编一个供初学者使用。

本文的程序改编自。。。

上面的算法与程序有很大的局限性：

（1）仅仅适用于相对简单的、只能从左到右、不能折弯的曲线插值，比如 ECharts 内的曲线插值。

（2）仅限于二维；

（3）计算效率比较低；

（4）数据点分布比较复杂的时候无法计算；

工业化的插值应该选用更好的、更高效的算法与代码，比如双三次B样条，T样条及 NURBS 等等。

2、三次样条曲线插值的算法与源代码

using System;
using System.Text;
using System.Collections;
using System.Collections.Generic;

namespace Legalsoft.Truffer
{
/// <summary>
/// 三次样条曲线插值（2D）
/// </summary>
public class Spline_Interpolation
{
private double[] Xi { get; set; }
private double[] Yi { get; set; }
private double[] A { get; set; }
private double[] B { get; set; }
private double[] C { get; set; }
private double[] H { get; set; }
private double[] Lamda { get; set; }
private double[] Mu { get; set; }
private double[] G { get; set; }
private int N { get; set; } = 0;
private int n { get; set; } = 0;
private double[] M;

public Spline_Interpolation()
{
}

/// <summary>
/// 应用演示代码
/// </summary>
public static void Drive()
{
#region 构建实验数据
Random rnd = new Random();
int num = 10;
double[] xx = new double[num];
double[] yy = new double[num];
for (int i = 0; i < num; i++)
{
xx[i] = i + rnd.NextDouble() * 0.1;
yy[i] = rnd.NextDouble();
}
#endregion

#region 曲线插值与绘制
Spline_Interpolation spline = new Spline_Interpolation();
spline.Initialize(xx, yy);
double xv = xx[0];
double last_xv = xv;
double last_yv = yy[0];
while (xv < xx[num - 1])
{
double yv = spline.Interpolate(xv);
if (last_xv > xv)
{
// 绘制线段
//g.DrawLine(pen, last_xv, last_yv, xv, yv);
}
last_xv = xv;
last_yv = yv;
xv += 0.01;
}
#endregion
}

public bool Initialize(double[] Xi, double[] Yi)
{
if (Xi.Length != Yi.Length)
{
return false;
}
if (Xi.Length == 0)
{
return false;
}

// 根据给定的Xi,Yi元素的数目来确定N的大小。
N = Xi.Length;
n = N - 1;

// 根据实际大小给各个成员变量分配内存
A = new double[N - 1];
B = new double[N];
C = new double[N - 1];

this.Xi = new double[N];
this.Yi = new double[N];

H = new double[N - 1];
Lamda = new double[N - 1];
Mu = new double[N - 1];

G = new double[N];

// 把输入数据点的数值赋给成员变量。
for (int i = 0; i <= n; i++)
{
this.Xi[i] = Xi[i];
this.Yi[i] = Yi[i];
}

// 求出hi,Lamda(i),Mu(i),gi
GetH();
GetLamda_Mu_G();
GetABC();

// 调用追赶法求出系数矩阵M
Chasing_Solver chase = new Chasing_Solver();
if (chase.Initialize(A, B, C, G))
{
if (chase.Solve(out M))
{
return true;
}
}
return false;
}

private void GetH()
{
for (int i = 0; i < n; i++)
{
H[i] = Xi[i + 1] - Xi[i];
}
}

private void GetLamda_Mu_G()
{
for (int i = 1; i < n; i++)
{
Lamda[i] = H[i] / (H[i] + H[i - 1]);
Mu[i] = 1.0 - Lamda[i];

double t1 = (Yi[i] - Yi[i - 1]) / H[i - 1];
double t2 = (Yi[i + 1] - Yi[i]) / H[i];
G[i] = 3.0 * (Lamda[i] * t1 + Mu[i] * t2);
}
G[0] = 3.0 * (Yi[1] - Yi[0]) / H[0];
G[n] = 3.0 * (Yi[n] - Yi[n - 1]) / H[n - 1];
Mu[0] = 1.0;
Lamda[0] = 0.0;
}

private void GetABC()
{
for (int i = 1; i < n; i++)
{
A[i - 1] = Lamda[i];
C[i] = Mu[i];
}
A[n - 1] = 1.0;
C[0] = 1.0;
for (int i = 0; i <= n; i++)
{
B[i] = 2.0;
}
}

private double fai0(double x)
{
double t1 = 2.0 * x + 1.0;
double t2 = (x - 1.0) * (x - 1.0);
return t1 * t2;
}

private double fai1(double x)
{
return x * (x - 1.0) * (x - 1.0);
}

private int GetSection(double x)
{
for (int i = 0; i < n; i++)
{
if (x <= Xi[i])
{
return i - 1;
}
}
return -999999;// iNum;
}

/// <summary>
/// 按x值计算y值（插值计算）
/// </summary>
/// <param name="x"></param>
/// <returns></returns>
public double Interpolate(double x)
{
double t = x;
int iNum = GetSection(x);
if (iNum == -1)
{
return Yi[0];
}
else if (iNum == -999999)
{
return Yi[n];
}
double P1 = (t - Xi[iNum]) / H[iNum];
double P2 = (Xi[iNum + 1] - t) / H[iNum];
return Yi[iNum] * fai0(P1) +
Yi[iNum + 1] * fai0(P2) +
M[iNum] * H[iNum] * fai1(P1) -
M[iNum + 1] * H[iNum] * fai1(P2);
}

}
}

3、求解严格对角占优的对角方程组的追赶法

using System;
using System.Text;
using System.Collections;
using System.Collections.Generic;

namespace Legalsoft.Truffer
{
/// <summary>
/// 求解严格对角占优的对角方程组的追赶法
/// </summary>
public class Chasing_Solver
{
/// <summary>
/// Dimension of Martrix Ax=d
/// </summary>
private int N { get; set; }
/// <summary>
/// Ax=d
/// </summary>
private double[] d { get; set; }
/// <summary>
/// a in A
/// </summary>
private double[] Aa { get; set; }
/// <summary>
/// b in A
/// </summary>
private double[] Ab { get; set; }
/// <summary>
/// c in A
/// </summary>
private double[] Ac { get; set; }
/// <summary>
/// LU-->L
/// </summary>
private double[] L { get; set; }
/// <summary>
/// LU-->U
/// </summary>
private double[] U { get; set; }
/// <summary>
/// store the result
/// </summary>
private double[] S { get; set; }

public Chasing_Solver()
{
}

/// <summary>
/// 数据初始化
/// </summary>
/// <param name="a"></param>
/// <param name="b"></param>
/// <param name="c"></param>
/// <param name="d"></param>
/// <returns></returns>
public bool Initialize(double[] a, double[] b, double[] c, double[] d)
{
int na = a.Length;
int nb = b.Length;
int nc = c.Length;
int nd = d.Length;

if (nb < 3)
{
return false;
}
N = nb;

if (na != N - 1 || nc != N - 1 || nd != N)
{
return false;
}
S = new double[N];
L = new double[N - 1];
U = new double[N];

Aa = new double[N - 1];
Ab = new double[N];
Ac = new double[N - 1];
this.d = new double[N];

for (int i = 0; i <= N - 2; i++)
{
Ab[i] = b[i];
this.d[i] = d[i];
Aa[i] = a[i];
Ac[i] = c[i];
}

Ab[N - 1] = b[N - 1];
this.d[N - 1] = d[N - 1];
return true;
}

/// <summary>
/// 追赶法求解器
/// </summary>
/// <param name="R"></param>
/// <returns></returns>
public bool Solve(out double[] R)
{
#region A=LU
R = new double[Ab.Length];
U[0] = Ab[0];
for (int i = 2; i <= N; i++)
{
// L[i] = Aa[i] / U[i - 1];
L[i - 2] = Aa[i - 2] / U[i - 2];
//U[i]=Ab[i]-Ac[i-1]*L[i];
U[i - 1] = Ab[i - 1] - Ac[i - 2] * L[i - 2];
}
#endregion

#region Ly=d
double[] Y = new double[d.Length];
Y[0] = d[0];

for (int i = 2; i <= N; i++)
{
//Y[k]=d[k]-L[k]*Y[k-1];
Y[i - 1] = d[i - 1] - (L[i - 2]) * (Y[i - 2]);

}
#endregion

#region Ux=Y
//X[n]=Y[n]/U[n];
R[N - 1] = Y[N - 1] / U[N - 1];
//X[k]=(Y[k]-C[k]*X[k+1])/U[k];(n-1,,.....1)
for (int i = N - 1; i >= 1; i--)
{
R[i - 1] = (Y[i - 1] - Ac[i - 1] * R[i]) / U[i - 1];
}
#endregion

for (int i = 0; i < R.Length; i++)
{
if (double.IsInfinity(R[i]) || double.IsNaN(R[i]))
{
return false;
}
}
return true;
}
}
}

C#，计算几何，计算机图形学，三次B样条曲线插值算法（B-Spline Curve Interpolation）的基本原理及源代码相关推荐

计算机图形学三(补充)：重心坐标(barycentric coordinates)详解及其作用
重心坐标(Barycentric Coordinates) 1 重心坐标的定义及求解 1.1 基础定义 1.2 几何面积角度求解 1.3 坐标系角度求解 2 重心坐标的运用 Reference (本篇 ...
计算机图形学--------充分理解B样条曲线
样条(spline)二字,从英文翻译过来的,让人费解.B样条的数学定义更是让人匪夷所思.看了好几本参考教材,还是把总结一下B样条这个概念. 一.解释什么是样条. 实际应用中,样条是一根富有弹性的细木条 ...
计算机图形学(三)-图形学中的基本变换（缩放、平移、旋转、剪切、镜像）
图形学中的基本变换 1. 二维变换 1.1 缩放变换 1.2.镜像变换 1.3 剪切变换 1.4 旋转变换 1.5 平移变换 1.5.1 什么是线性变换 1.5.2 平移变换(仿射变换) 1.5.3 ...
计算机图形学三：光栅化-Rasterization
文章目录什么是光栅化? 像素和屏幕直线光栅化(Linear Rasterization) DDA数值微分算法中点Bresenham算法三角形光栅化(Triangle Rasterization ...
计算机图形学(三) -- 3D 变换
文章目录 3D 变换缩放(Scale) 平移(Translation) 旋转(Rotation) 3D 旋转(3D Rotation) 什么是欧拉角罗德里格斯旋转公式(Rodrigues' Rot ...
计算机图形学学习笔记（九）：曲线曲面（一）：参数曲线、参数几何代数形式
接上文计算机图形学学习笔记(八):三维图形变换:三维几何变换,投影变换(平行/ 透视投影) 计算机图形学三大块内容:光栅图形显示(前面已经介绍完了 1-8).几何造型技术.真实感图形显示.光栅图 ...
计算机图形什么叫参数连续性,计算机图形学--参数三次插值样条曲线.ppt
计算机图形学--参数三次插值样条曲线参数三次插值样条曲线三次多项式方程是能表示曲线段的端点通过特定点且在连接处保持位置和斜率的连续性的最低阶次的方程.与更高次的多项式方程相比,三次样条只需要较少的 ...
计算机图形学——bazier曲线和B样条曲线
目录前言一.bazier曲线 1.bazier曲线的由来二.B-spline 1.为何要引入B-spline曲线,Bazier曲线有何不足 2.B-spline曲线方法前言写这篇文章是因为最 ...
三个基本原理和概念 - 计算机图形学、数据加密、数据挖掘
一. 计算机图形学最基本原理计算机屏幕由像素组成.一个像素点包括X和Y坐标. 高级语言有画基本图形的函数或语句,可以直接调用画图形.比如画线,画圆,画四方形. 但是最底层的编程接口, ...

C#，计算几何，计算机图形学，三次B样条曲线插值算法（B-Spline Curve Interpolation）的基本原理及源代码

1、发文初心

2、三次样条曲线插值的算法与源代码

3、求解严格对角占优的对角方程组的追赶法

C#，计算几何，计算机图形学，三次B样条曲线插值算法（B-Spline Curve Interpolation）的基本原理及源代码相关推荐

最新文章

热门文章