B样条插入节点

参考《计算几何算法与实现》–孔令德

通过插入新节点，可以提高B样条的局部修改性。插入一个节点后，需要重新计算节点分布及新控制点，以三次样条绘制为例，绘制上述形状：

#pragma once
#include"P2.h"

class Knot_Insertion_BSpline
{
public:
	Knot_Insertion_BSpline(void);
	~Knot_Insertion_BSpline(void);

public:
	//计算基函数
	 double BasicFunctionValue(double t, int i, int k, double knot[]);

	 //初始化控制点，节点矢量
	 void InitPoint();

	 //绘制样条线
	 void DrawCurve(CDC*pDC, double knot[], CP2 point[]);

	 //绘制控制点
	 void DrawControlPoint(CDC*pDC, CP2 point[]);

private:
	double knot1[12];	//定义节点矢量数组
	double knot2[13];   //插入后的节点矢量
	int n;              //控制点个数减1
	int k;              //次数
	int index;          //左下标索引
	CP2 P[8];           //控制点坐标
	CP2 Q[9];          //插入后的新控制点坐标
public:
	//计算插入后的节点
	int KnotInsertion(double theNewKnot);

	//计算插入节点后的新控制点
	void CalculateControlPnt(double theNewKnot, int r);

	//绘制
	void Draw(CDC* pDC, int theChoice);
};

#include "StdAfx.h"
#include "Knot_Insertion_BSpline.h"
#include "math.h"
#define ROUND(d) int(d+0.5)//四舍五入宏定义

Knot_Insertion_BSpline::Knot_Insertion_BSpline(void)
{
	InitPoint();
}


Knot_Insertion_BSpline::~Knot_Insertion_BSpline(void)
{
}

/*----------------------------------------------------
input   : t, i, k分别为参数值，支撑区间左端节点的下标，曲线次数，节点
function: 计算第i个k次B样条基函数值
-----------------------------------------------------*/

double Knot_Insertion_BSpline::BasicFunctionValue(double t, int i, int k, double knot[])
{
	//次数为0时的基函数值
	if(k==0)
	{
		if(t>=knot[i]&&t<knot[i+1])
			return 1.0;
		else 
			return 0.0;
	}

	//次数不为0时的基函数计算
	double value1, value2, value;
	if(k>0)
	{
		//节点矢量位于定义域之外
		if(t<knot[i]||t>knot[i+k+1])
			return 0.0;

		//节点矢量位于定义域之外
		else
		{
			double coff1, coff2;    //系数
			double denominator=0.0; //分母

			//递推公式第一项
			denominator=knot[i+k]-knot[i];
			if(denominator==0.0)
				coff1=0.0;                  //重节点分母为0
			else
				coff1=(t-knot[i])/denominator;

			//递推公式第二项
			denominator=knot[i+k+1]-knot[i+1];
			if(denominator==0.0)
				coff2=0.0;                  //重节点分母为0
			else
				coff2=(knot[i+k+1]-t)/denominator;

			//递推公式第一、二项的值
			value1=coff1*BasicFunctionValue(t,i,k-1,knot);
			value2=coff2*BasicFunctionValue(t,i+1,k-1,knot);

			value=value1+value2;
		}
	}
	return value;
}

/*----------------------------------------------------
input   : 
function: 初始化节点与控制点
-----------------------------------------------------*/
void Knot_Insertion_BSpline::InitPoint()
{
	n=7, k=3;

	//节点初始化
	knot1[0]=0;   knot1[1]=0;    knot1[2]=0;  
	knot1[3]=0;	  knot1[4]=0.2;  knot1[5]=0.4; 
	knot1[6]=0.6; knot1[7]=0.8;  knot1[8]=1.0;
	knot1[9]=1.0; knot1[10]=1.0; knot1[11]=1.0; 

	//控制点初始化
	P[0]=CP2(-300, -80);
	P[1]=CP2(-200, -20);
	P[2]=CP2(-100, -160);
	P[3]=CP2(170, -160); 
	P[4]=CP2(250, 0);
	P[5]=CP2(150, 160);
	P[6]=CP2(-80, 160);
	P[7]=CP2(-160, 40);
}

/*----------------------------------------------------
input   : 节点、控制点
function: 绘制相应节点与控制点的曲线
-----------------------------------------------------*/
void Knot_Insertion_BSpline::DrawCurve(CDC*pDC, double knot[], CP2 point[])
{
	//创建并选取画笔
	CPen NewPen, *pOldPen;
	NewPen.CreatePen(PS_SOLID,2,RGB(0,0,255));
	pOldPen=pDC->SelectObject(&NewPen);

	//等距取参数点计算值
	double tStep=0.01;
	for(double t=0.0; t<=1.0; t+=tStep)
	{
		CP2 p(0,0);//离散点
		for(int i=0; i<=n; i++)
		{
			double BValue=BasicFunctionValue(t, i, k, knot);
			p+=point[i]*BValue;
		}
		if(t==0)
			pDC->MoveTo(ROUND(p.x), ROUND(p.y));
		else
			pDC->LineTo(ROUND(p.x), ROUND(p.y));
	}

	//还原画笔，删除创建的画笔
	pDC->SelectObject(pOldPen);	
	NewPen.DeleteObject();      
}

/*----------------------------------------------------
input   : 
function: 绘制控制点
-----------------------------------------------------*/
void Knot_Insertion_BSpline::DrawControlPoint(CDC*pDC, CP2 point[])
{
	//初始化画笔
	CPen NewPen,*pOldPen;
	NewPen.CreatePen(PS_SOLID,3,RGB(0,0,0));
	pOldPen=pDC->SelectObject(&NewPen);

	//初始化画刷
	CBrush NewBrush,*pOldBrush;
	NewBrush.CreateSolidBrush(RGB(0,0,0));
	pOldBrush=pDC->SelectObject(&NewBrush);

	//加粗控制点
	for(int i=0;i<=n;i++)
	{
		if(0==i)
		{
			pDC->MoveTo(ROUND(point[i].x),ROUND(point[i].y));
			pDC->Ellipse(ROUND(point[i].x)-5,ROUND(point[i].y)-5,ROUND(point[i].x)+5,ROUND(point[i].y)+5);
		}
		else
		{
			pDC->LineTo(ROUND(point[i].x),ROUND(point[i].y));
			pDC->Ellipse(ROUND(point[i].x)-5,ROUND(point[i].y)-5,ROUND(point[i].x)+5,ROUND(point[i].y)+5);
		}
	}

	//还原画笔画刷
	pDC->SelectObject(pOldBrush);
	pDC->SelectObject(pOldPen);
	NewPen.DeleteObject(); 
	NewBrush.DeleteObject(); 
}

/*----------------------------------------------------
input   : 新节点
output  : 插入节点区间的左下标索引号
function: 插入新节点后，重新计算节点
-----------------------------------------------------*/
int Knot_Insertion_BSpline::KnotInsertion(double theNewKnot)
{
	//原节点数据赋值给新节点
	for(int i=0; i<=n+k+1; i++)
	{
		knot2[i]=knot1[i];
	}

	//在定义域内查找插入节点的位置
	for(int i=k; i<=n+1; i++)  
	{
		if(theNewKnot<knot1[i])
		{
			index=i;            //插入节点区间右端索引号
			break;
		}
	}

	//插入新节点
	int j=index;
	knot2[j++]=theNewKnot;

	//后续节点重新编号
	for(int i=index; i<=n+k+1; i++,j++)
	{
		knot2[j]=knot1[i];
	}

	//返回插入节点区间左端索引号
	index--;
	return index-1;
}

/*----------------------------------------------------
input   : 新节点,节点原来的重复度
output  : 
function: 插入新节点后，重新计算控制点
-----------------------------------------------------*/
void Knot_Insertion_BSpline::CalculateControlPnt(double theNewKnot, int r)
{
	//前不变控制点
	for(int i=0; i<=index-k; i++)
		Q[i]=P[i];

	//插入后的新控制点
	for(int i=index-k+1; i<=index-r; i++)
	{
		//系数计算
		double numerator=theNewKnot-knot1[i];
		double denominator=knot1[i+k]-knot1[i];

		double alpha=0.0;
		if(numerator==0||denominator==0)
			alpha=0.0;
		alpha=numerator/denominator;
		Q[i]=alpha*P[i]+(1-alpha)*P[i-1];
	}

	//后不变控制点
	for(int i=index-r+1; i<=n+1;i++)
		Q[i]=P[i-1];
}

/*----------------------------------------------------
input   : 绘制的类型
output  : 
function: 插入节点前后的样条线绘制
-----------------------------------------------------*/
void Knot_Insertion_BSpline::Draw(CDC* pDC, int theChoice)
{
	switch(theChoice)
	{
	case 0:
		{
			//插入节点前
			DrawControlPoint(pDC, P);
			DrawCurve(pDC,knot1, P);
			break;
		}
	case 1:
		{
			//插入非重复节点后
			KnotInsertion(0.5);
			CalculateControlPnt(0.5, 0);
			n=8;
			DrawControlPoint(pDC, Q);
			DrawCurve(pDC,knot2, Q);
			break;
		}
	case 2:
		{
			//插入重复节点
			KnotInsertion(0.4);
			CalculateControlPnt(0.4, 1);
			n=8;
			DrawControlPoint(pDC, Q);
			DrawCurve(pDC,knot2, Q);
			break;
		}
	}
}