自然三次样条曲线Java实现_雨林木风博客

http://blog.sina.com.cn/u/5308127060

首页博文目录关于我

个人资料

微博

加好友发纸条

写留言加关注

博客等级：
博客积分：

博客访问：
关注人气：
获赠金笔：0支
赠出金笔：0支
荣誉徽章：

正文字体大小：大中小

自然三次样条曲线Java实现

(2017-01-09 17:09:32)

标签：

自然三次样条曲线

java

分类：软件开发

自然三次样条曲线是三次样条曲线的一种，比较常用，自然三次样条是对三次样条曲线做了一些限制条件。简单来说，自然三次样条就是针对两点之间的曲线，本来两点之间是直接画线段的，为了平滑，现在要把这条线段变成曲线，这条曲线用的就是自然三次样条曲线。如果一个图中有n个点，那么就需要n-1个曲线方程。

S1(x)=y1+b1(x-x1)+c1(x-x1)^2+d1(x-x1)^3 [x1,x2]

S2(x)=y2+b2(x-x2)+c2(x-x2)^2+d2(x-x2)^3 [x2,x3]

.....

Sn-1=yn-1 + bn-1 (x-xn-1)+ cn-1 (x-xn-1)^2 + dn-1 (x-xn-1)^3 [xn-1, xn]

自然三次样条的三个性质：

（1）Si(xi)=yi, S(xi+1)=yi+1 , i=1,...,n-1

（2）S'i-1(xi) = S'i(xi), i=2,....,n-1

（3）S''i-1(xi) = Si''(xi), i=2,...,n-1

性质1保证S(x)插值数据点。

性质2使得相邻的样条段在他们相遇的地方斜率相等。

性质3保证在两条样条段相邻的地方曲率相同。该曲率由二阶导数表示。

比如有A(1, 2), B(2, 1), C(4, 4), D(5, 3)这4个数据点，则需要

S1(x)=2-13/8*(x-1)+0*(x-1)^2+5/8*(x-1)^3 A(1, 2)和B(2, 1)

S2(x)=1+1/4*(x-2)+15/8*(x-2)^2-5/8*(x-2)^3 B(2, 1)和C(4, 4)

S3(x)=4+1/4*(x-4)-15/8*(x-4)^2+5/8*(x-4)^3 C(4, 4)和D(5, 3)

读者可以试着把这四个点代入指定的方程内验证一下。

现在主要的任务就是求解这n-1个方程的bi, ci, di系数，求解出来的话，这n-1个方程就确定下来了，那么画曲线方程就不是事了。

求解3*（n-1）个系数，就需要3*（n-1）个方程。

性质1可以提供n-1个方程（S1(x2)=y2, ..., Sn-1(xn)=yn，共计n-1个）

性质2可以提供n-2个方程（S1'(x2)-S2'(x2)=0，...，Sn-2(xn-1)-Sn-1(xn-1)=0，共n-2个)）

性质3可以提供n-2个方程（S1''(x2)-S2''(x2)=0，...，Sn-2''(xn-1)-Sn-1''(xn-1)=0，共n-2个）

好像还差2个方程才能求解吧。那么就用S1''(x1)=0，S''n-1(xn)=0这两个方程来填补空缺，因此有这两个额外特殊方程的三次样条叫做自然三次样条。

具体化简结果请看数值分析的书籍，我也是跑到读书馆看了数值分析看了一下午，在本子上慢慢的一点点推导出来的，好多东西在博客里不好展示，我只是把流程给你过一遍，请谅解。

下面给出自然三次样条的伪代码：

给定x=[x1,x2,....,xn] ,其中x1<...

for i=1, ..., n-1

{

ai=yi

Ai=xi+1 - xi

Bi=yi+1 - yi

}

先求解ci，ci系数矩阵是一个严格对角占优矩阵，利用追赶法进行求解

求解矩阵，得到c1，c2...，cn（注意c1=0，cn=0）

for i=1,...,n-1

{

di=(ci+1 - ci)/3/Ai

bi=(Bi-Ai)-Ai/3*(2ci+ ci+1)

}

自然三次样条曲线方程如下：

Si(x)=ai+bi(x-xi)+ci(x-xi)^2+di(x-xi)^3, [xi, xi+1], i=1, ...., n-1

java代码如下：

package DrawLine;

public class SplineInterpolation {

double x[], y[];//采样点的(x,y)

int n = 0;//采样点的数目

//s1(x)=y1+b1(x-x1)+c1(x-x1)^2+d1(x-x1)^3

//只需要求出a[], b[], c[], d[]就可以得出n-1个自然三次样条曲线方程

double a[], b[], c[], d[];

//y1=a1

//xx=xi+1 - xi

//yy=yi+1 - yi

double xx[], yy[];

public SplineInterpolation(double x[], double y[], int n)

{

this.n = n;

this.x = new double[this.n];

this.y = new double[this.n];

this.x = x;

this.y = y;

this.a = new double[this.n-1];

this.b = new double[this.n-1];

this.c = new double[this.n];

this.d = new double[this.n-1];

this.xx = new double[this.n-1];

this.yy = new double[this.n-1];

}

public void Interpolation()

{

//解三角线方程组 Ax=f 可化为求解两个三角形方程组 => A=PQ

//Py=f Qx=y

//p1=b1 qi=ci/pi(i=1,2...,n-1) pi=bi-ai*qi-1(i=2,3...,n)

// bi相当于(xx[i-1]+xx[i])*2 ci相当于xx[i] ai相当于xx[i-1]

//此处的a, b, c和上边的方程系数不是一回事

double a2[], b2[], c2[];

a2 = new double[n];

b2 = new double[n];

c2 = new double[n-1];

double f[];

f = new double[n];

for(i=0; i

{

if(i==0)

{

b2[i] = 1;

c2[i] = 0;

f[i] = 0;

}else if(i==n-1)

{

a2[i] = 0;

b2[i] = 1;

f[i] = 0;

}else{

a2[i] = xx[i-1];

b2[i] = 2 * (xx[i-1] + xx[i]);

c2[i] = xx[i];

f[i] = 3 * (yy[i]/xx[i] - yy[i-1]/xx[i-1]);

}

double p[], q[];

p = new double[n];

q = new double[n-1];

p[0] = b2[0];

for(i=0; i < n; i++)

{

if(i > 0)

{

p[i] = b2[i] - a2[i]*q[i-1];

}

if(i < n-1)

{

q[i] = c2[i]/p[i];

}

//求解Py=f

//y1=f1/p1

//yi=(fi-ai*yi-1)/pi (i=2,3...,n)

double y2[];

y2 = new double[n];

y2[0] = f[0]/p[0];

for(i=1; i

{

y2[i] = (f[i] - a2[i]*y2[i-1])/p[i];

}

//求解Qx=y x就是要求的c[]

//xn=yn

//xi=yi-qi*xi+1

c[n-1] = y2[n-1];

for(i=n-2; i>=0; i--)

{

c[i] = y2[i] - q[i]*c[i+1];

}

public double computeY(double inputX)

{

int i = 0;

for(i=0; i

{

if(inputX >= x[i] && inputX < x[i+1])

{

break;

}

//计算相应曲线中的Y坐标

double returnY = 0.0;

double inputX_x_1 = inputX-x[i];

double inputX_x_2 = inputX_x_1 * inputX_x_1;

double inputX_x_3 = inputX_x_1 * inputX_x_2;

returnY = a[i] + b[i]*inputX_x_1 + c[i]*inputX_x_2 + d[i]*inputX_x_3;

return returnY;

}

package DrawLine;

import java.awt.BasicStroke;

import java.awt.Color;

import java.awt.Graphics;

import java.awt.Graphics2D;

import java.awt.event.MouseEvent;

import java.awt.event.MouseListener;

import java.awt.event.MouseMotionListener;

import javax.swing.JFrame;

import javax.swing.JPanel;

@SuppressWarnings("serial")

public class DrawSplineLine extends JPanel implements MouseListener, MouseMotionListener{

public SplineInterpolation sp;//声明对象

int n = 0;//原采样点集合的大小

int maxNum = 10000; //扩大后的采样点集合的大小，目前设为10000

double maxY[]; //扩大后的采样点y集合

double maxX[]; //扩大后的采样点x集合

double x[]; //源采样点x集合

double y[]; //源采样点y集合

public DrawSplineLine(double x[], double y[], int n, SplineInterpolation sp)

{

this.sp = sp;

this.n = n;

this.x = new double[n];

this.y = new double[n];

this.x = x;

this.y = y;

//鼠标事件监听，一定不要忘加这两个函数

addMouseListener(this);

addMouseMotionListener(this);

}

public void Paint(Graphics2D g)

{

g.setColor(Color.RED);

g.setStroke(new BasicStroke(4.0f));

for(int i=0; i

{

//将每个坐标点扩大为自身的100倍，之后在屏幕上显示

g.drawLine((int)(maxX[i]*100), (int)(maxY[i]*100), (int)(maxX[i+1]*100), (int)(maxY[i+1]*100));

}

@Override

protected void paintComponent(Graphics arg0) {

Graphics2D g = (Graphics2D) arg0;

g.setColor(Color.WHITE);

g.fill3DRect(0, 0, getWidth(), getHeight(), true);

Paint(g);

}

public void ChangeXY()

{

maxNum = 10000;

maxX = new double[maxNum];

maxY = new double[maxNum];

int maxFlag = 0;

for(int i=0; i

{

maxX[maxFlag] = x[i];

maxY[maxFlag] = y[i];

if(i

{

while(maxX[maxFlag]+0.05 < x[i+1])

{

maxFlag ++;

maxX[maxFlag] = maxX[maxFlag-1] + 0.05;

maxY[maxFlag] = sp.computeY(maxX[maxFlag]);

}

maxFlag ++;

}

maxNum = maxFlag;

}

@Override

public void mouseDragged(MouseEvent e) {

int i = 0;

for(i=0; i

{

if(Math.abs(e.getX()-x[i]*100)<20 && Math.abs(e.getY()-y[i]*100)<20)

{

//生成新的采样点

x[i] = (double)e.getX()/100;

y[i] = (double)e.getY()/100;

break;

}

//输出鼠标焦点的坐标

System.out.println("x="+e.getX()+":y="+e.getY());

//实时的更新大采样点集合

ChangeXY();

//重画调用函数

repaint();

}

public static void main(String[] args) {

//测试数据 x[] y[]

double x[] = {1, 2, 4, 5};

double y[] = {2, 1, 4, 3};

int n = 4; //测试数据集的大小

SplineInterpolation sp = new SplineInterpolation(x, y, n);

sp.Interpolation();

//输出n-1个曲线方程的a,b,c,d参数

//曲线方程 Sn(x)=a+b(x-xn)+c(x-xn)^2+d(x-xn)^3

for(int i=0; i

{

System.out.println(sp.a[i]+" "+sp.b[i]+" "+sp.c[i]+" "+sp.d[i]);

}

//*************绘制曲线部分***********************

DrawSplineLine ds = new DrawSplineLine(x, y, n, sp);

//根据源采样点集合生成大的数据集

ds.ChangeXY();

//============画图老套路===========

JFrame jf = new JFrame("自然三次样条插值曲线(鼠标可拖拽版)");

jf.setSize(800, 800);

jf.setDefaultCloseOperation(JFrame.EXIT_ON_CLOSE);

jf.setVisible(true);

jf.add(ds);

//===============================

//*************************************

}

//***********************无用函数***********************

@Override

public void mouseMoved(MouseEvent e) {

}

@Override

public void mouseClicked(MouseEvent e) {

}

@Override

public void mousePressed(MouseEvent e) {

}

@Override

public void mouseReleased(MouseEvent e) {

}

@Override

public void mouseEntered(MouseEvent e) {

}

@Override

public void mouseExited(MouseEvent e) {

}

运行结果：（曲线支持鼠标拖动）

http://s10/mw690/005NemkAzy77QUnYgm5a9&690

http://s14/mw690/005NemkAzy77QUonUlv5d&690

阅读┊ 收藏 ┊ 喜欢 ▼ ┊打印┊举报/Report

前一篇：机器学习实战-手写数字识别

后一篇：分类算法之朴素贝叶斯

新浪BLOG意见反馈留言板　欢迎批评指正