MFC绘制三次样条曲线

参考孔令德老师的书籍《计算几何算法与实现》

算法部分在书中有详细的讲解，下面放上代码实现

1、新建MFC单文档项目
2、添加二维点类CP2：

#pragma once
class CP2
{
public:
	CP2(void);
	CP2(double x,double y);
	~CP2(void);
public:
	double x;
	double y;
	double w;
};

#include "StdAfx.h"
#include "P2.h"


CP2::CP2(void)
{
	x=0;
	y=0;
	w=1;
}
CP2::CP2(double x,double y)
{
	this->x=x;
	this->y=y;
}

CP2::~CP2(void)
{

}

3、在View类中包含头文件，并添加如下的成员属性

#include"P2.h"
#define ROUND(d) int(d+0.5)//四舍五入取整

public:
	CP2 P[7];
	void ReadPoint(void);
	void DrawDataPoint(CDC*pDC);
	void DrawCubicSpline(CDC*pDC);

void CcubicSplineView::ReadPoint(void)
{
	P[1].x=-340,P[1].y=-200;//读入型值点，按照横坐标从小到大排序
	P[2].x=-150,P[2].y=0;
	P[3].x=0,P[3].y=-50;
	P[4].x=100,P[4].y=-100;
	P[5].x=250,P[5].y=-100;
	P[6].x=350,P[6].y=-50;
}

void CcubicSplineView::DrawDataPoint(CDC*pDC)
{
	CBrush NewBrush,*pOldBrush;//加粗型值点
	NewBrush.CreateSolidBrush(RGB(0,0,0));
	pOldBrush=pDC->SelectObject(&NewBrush);
	for(int i=1;i<7;i++)
	{
		pDC->Ellipse(ROUND(P[i].x-5),ROUND(P[i].y-5),
			ROUND(P[i].x+5),ROUND(P[i].y+5));
	}
	pDC->SelectObject(pOldBrush);
}

void CcubicSplineView::DrawCubicSpline(CDC*pDC)//绘制三次样条线
{
	int n=6;
	const int dim=7;//二维数组维数
	double b1=10,bn=-10;//给出起点和终点的一阶导数
	double h[dim],lambda[dim],mu[dim],D[dim];//四个数组
	double l[dim],m[dim],u[dim];//追赶法参数
	double M[dim],K[dim];//追赶法过渡矩阵
	double a[dim],b[dim],c[dim],d[dim];//函数的系数
	for(int i=1;i<n;i++)//计算hi=xi+1-xi
	{
		h[i]=P[i+1].x-P[i].x;
	}
	for(int i=2;i<n;i++)
	{
		lambda[i]=h[i-1]/(h[i-1]+h[i]);
		mu[i]=h[i]/(h[i-1]+h[i]);
		D[i]=6/(h[i-1]+h[i])*((P[i+1].y-P[i].y)/h[i]-(P[i].y-P[i-1].y)/h[i-1]);//计算Di
	}
	D[1]=6*((P[2].y-P[1].y)/h[1]-b1)/h[1];//夹持端的D[1]
	D[n]=6*(bn-(P[n].y-P[n-1].y)/h[n-1])/h[n-1];//夹持端的D[n]
	mu[1]=1;
	lambda[n]=1;
	//追赶法求解三弯矩方程
	l[1]=2;
	u[1]=mu[1]/l[1];
	for(int i=2;i<=n;i++)
	{
		m[i]=lambda[i];
		l[i]=2-m[i]*u[i-1];
		u[i]=mu[i]/l[i];
	}
	K[1]=D[1]/l[1];//解LK=D方程
	for(int i=2;i<=n;i++)
	{
		K[i]=(D[i]-m[i]*K[i-1])/l[i];
	}
	M[n]=K[n];//解UM=K
   for(int i=n-1;i>=1;i--)
   {
	   M[i]=K[i]-u[i]*M[i+1];
   }
   //计算三次样条函数的系数
   for(int i=1;i<n;i++)
   {
	   a[i]=P[i].y;
	   b[i]=(P[i+1].y-P[i].y)/h[i]-h[i]*(M[i]/3+M[i+1]/6);
	   c[i]=M[i]/2;
	   d[i]=(M[i+1]-M[i])/(6*h[i]);
   }
   pDC->MoveTo(ROUND(P[1].x),ROUND(P[1].y));
   double xStep=0.5;//当前步长
   double x,y;//当前点
   for(int i=1;i<n;i++)//循环访问每个节点
   {
	   for(x=P[i].x;x<P[i+1].x;x+=xStep)//按照步长xStep进行计算
	   { 
		   y=a[i]+b[i]*(x-P[i].x)+c[i]*(x-P[i].x)*(x-P[i].x)+d[i]*(x-P[i].x)*(x-P[i].x)*(x-P[i].x);
		   pDC->LineTo(ROUND(x),ROUND(y));//绘制样条曲线
	   }
   }
}

4、构造函数初始化：

CcubicSplineView::CcubicSplineView()
{
	// TODO: 在此处添加构造代码
	ReadPoint();

}

5、OnDraw函数进行调用绘图

void CcubicSplineView::OnDraw(CDC* pDC)
{
	CcubicSplineDoc* pDoc = GetDocument();
	ASSERT_VALID(pDoc);
	if (!pDoc)
		return;

	// TODO: 在此处为本机数据添加绘制代码
	CRect rect;//定义客户区矩形
	GetClientRect(&rect);//获得客户区的大小
	pDC->SetMapMode(MM_ANISOTROPIC);//pDC自定义坐标系
	pDC->SetWindowExt(rect.Width(),rect.Height());//设置窗口范围
	pDC->SetViewportExt(rect.Width(),-rect.Height());//设置视区范围,x轴水平向右，y轴垂直向上
	pDC->SetViewportOrg(rect.Width()/2,rect.Height()/2);//客户区中心为原点
	rect.OffsetRect(-rect.Width()/2,-rect.Height()/2);

	DrawDataPoint(pDC);
	DrawCubicSpline(pDC);
}