bsplinedoc.cpp

B3次样条拟合

// BSplineDoc.cpp : implementation of the CBSplineDoc class
//

#include "stdafx.h"
#include "BSpline.h"
#include "math.h"

#include "BSplineDoc.h"

#ifdef _DEBUG
#define new DEBUG_NEW
#undef THIS_FILE
static char THIS_FILE[] = __FILE__;
#endif

/////////////////////////////////////////////////////////////////////////////
// CBSplineDoc

IMPLEMENT_DYNCREATE(CBSplineDoc, CDocument)

BEGIN_MESSAGE_MAP(CBSplineDoc, CDocument)
	//{{AFX_MSG_MAP(CBSplineDoc)
		// NOTE - the ClassWizard will add and remove mapping macros here.
		//    DO NOT EDIT what you see in these blocks of generated code!
	//}}AFX_MSG_MAP
END_MESSAGE_MAP()

/////////////////////////////////////////////////////////////////////////////
// CBSplineDoc construction/destruction

CBSplineDoc::CBSplineDoc()
{
	// TODO: add one-time construction code here
	m_type = 2;
	m_subtype = 1;
	m_degree = 3;
	m_length = 0.01;
	m_ptArray.SetSize(0,10);
	m_lthArray.SetSize(1,10);

	r.SetSize(0,5);
	t.SetSize(0,5);
}

CBSplineDoc::~CBSplineDoc()
{
}

BOOL CBSplineDoc::OnNewDocument()
{
	if (!CDocument::OnNewDocument())
		return FALSE;

	// TODO: add reinitialization code here
	
	// (SDI documents will reuse this document)

	return TRUE;
}



/////////////////////////////////////////////////////////////////////////////
// CBSplineDoc serialization

void CBSplineDoc::Serialize(CArchive& ar)
{
	if (ar.IsStoring())
	{
		// TODO: add storing code here
	}
	else
	{
		// TODO: add loading code here
	}
}

/////////////////////////////////////////////////////////////////////////////
// CBSplineDoc diagnostics

#ifdef _DEBUG
void CBSplineDoc::AssertValid() const
{
	CDocument::AssertValid();
}

void CBSplineDoc::Dump(CDumpContext& dc) const
{
	CDocument::Dump(dc);
}
#endif //_DEBUG

/////////////////////////////////////////////////////////////////////////////
// CBSplineDoc commands

CPoint CBSplineDoc::Dbordv(double u)
/*
    输入参数:
	k表示次数，l表示相异节点个数，u表示参数值，t[l]表示相异节点数组，r[l]表示重复度数组，
	d[dim][kl]表示所涉及的控制顶点。
	输出参数:
	p表示B样条曲线上参数为u值的一点
*/
{
	int k = Knot_Index(u) - m_degree;  //相关节点的起始点索引号
	if(k < 0) return CPoint(0,0);

	CArray<double,double> dt0,dt1;
	dt0.SetSize(m_degree + 1);
	dt1.SetSize(m_degree + 1);

	for(int i = 0;i <= m_degree;i++)
	{
	    dt0[i] = m_ptArray[k + i].x;
		dt1[i] = m_ptArray[k + i].y;
	}
	

    for(int II = 1;II <= m_degree;II++)
	{
		for(int j = 0;j <= m_degree - II;j++)
		{
			double denom,alpha;
			denom = Knot_Value(j+m_degree+1+k) - Knot_Value(j+II+k);
			
			if(fabs(denom) <= 1e-5)
				alpha = 0.0;
			else 
				alpha = (u - Knot_Value(j+II+k))/denom;

			dt0[j] = (1.0 - alpha) * dt0[j] + alpha * dt0[j+1];
			dt1[j] = (1.0 - alpha) * dt1[j] + alpha * dt1[j+1];
		}
	}

	CPoint p;
	p.x = (int)(dt0.GetAt(0)+0.5);
	p.y = (int)(dt1.GetAt(0)+0.5);

	return p;

}

double CBSplineDoc::Knot_Value(int j)
/*
功能: 给定节点矢量的一个下标值,求得该节点的节点值
输入参数: l表示相异节点个数,x表示索引号为该值的节点值,
          t[l]表示相异节点数组，r[l]表示重复度数组
*/
{
	double u;
	int index = 0;
	int l = t.GetSize();
	for(int i = 0;i < l;i++)
	{
		index = index + r[i];
		if((index - 1) >= j)
		{
			u = t[i];
			break;
		}
	}
	return u;
}

int CBSplineDoc::Knot_Index(double u)
/*
功能: 给定u值,求得该u值在哪两个节点之间
*/
{

	int index = 0; 
	int l = t.GetUpperBound();

	if(l == -1 || u < 0 || u > t.GetAt(l)) return -1;  //t内无值或者u〈0 时，返回-1

	if(u == t.GetAt(l))           //当u超过最大值时
	{
		for(int i = 0;i < r.GetSize()-1;i++) index += r[i];
		return index - 1;
	}

	for(int i = 0;t[i] <= u && i < t.GetSize()-1;i++)
		index += r[i];
	return index - 1;
}

void CBSplineDoc::UpdateData()   //根据控制点及其他参数重新设置节点相异数组
{
	t.RemoveAll();
	r.RemoveAll();

	int i,l = m_ptArray.GetSize();
	double lsum = 0.0,lt = 0.0,ut = 0.0;  //非均匀B样条计算边长用到的参数
	if(l < m_degree+1) return;   //控制点数 少于曲线次数加一时,不进行更新

	switch(m_type) {
	case 1: //均匀B样条的节点划分
		for(i = 0;i <= l + m_degree;i++)
		{
			t.Add(i);
			r.Add(1);
		}
		break;
	case 2: //准均匀B样条的节点划分
		t.Add(0);
		r.Add(m_degree+1);
		for(i = 1;i <= l - 1 - m_degree;i++)
		{
			t.Add(i);
			r.Add(1);
		}
		t.Add(l - m_degree);
		r.Add(m_degree+1);
		break;
	case 3:  //分段贝奇尔曲线的节点划分	
		if((l - 1) % m_degree != 0) return;
		t.Add(0);
		r.Add(m_degree+1);
		for(i = 1;i < (l - 1) / m_degree;i++)
		{
			t.Add(i);
			r.Add(m_degree);
		}
		t.Add((l - 1) / m_degree);
		r.Add(m_degree+1);
		break;
	case 4:   //  非均匀B样条的节点划分
		m_lthArray.RemoveAll();
		m_lthArray.Add(0); 
		for(i = 1;i <= l - 1;i++)
		{
			m_lthArray.Add(Length(i));
		}
		if(m_subtype == 1)  //里森费尔德方法
		{
			for(i = 1;i <= l - 1;i++)	
			{
				lsum += m_lthArray[i];
			}
			t.Add(0);
			r.Add(m_degree+1);
			if(m_degree % 2 == 0)
			{
				for(int k = 1;k < m_degree / 2;k++)
					lt += m_lthArray[k];
				for(int j = 0;j <= l - m_degree - 2;j++)
				{
					lt += m_lthArray[m_degree / 2 + j];
					t.Add((lt + m_lthArray[m_degree / 2 + j + 1]) / lsum);
					r.Add(1);
				}
			}
			else
			{
				for(int k = 1;k < (m_degree + 1) / 2;k++)
					lt += m_lthArray[k];
				for(int j = 0;j <= l - m_degree - 2;j++)
				{
					lt += m_lthArray[(m_degree + 1) / 2 + j];
					t.Add(lt / lsum);
					r.Add(1);
				}
			}
			t.Add(1);
			r.Add(m_degree+1);
		}
		else if(m_subtype == 2) //哈特利-贾德方法
		{
			lsum = 0.0;	
			for(i = m_degree + 1;i <= l;i++)
			{
				for(int j = i - m_degree;j <= i - 1;j++)
				{
					lsum += m_lthArray[j];
				}
			}
			t.Add(0);
			r.Add(m_degree+1);
			for(i = m_degree + 1;i <= l - 1;i++)
			{
				lt = 0.0;
				for(int j = i - m_degree;j <= i - 1;j++)
				{
					lt += m_lthArray[j];
				}
				ut += (lt / lsum);
				t.Add(ut);
				r.Add(1);
			}
			t.Add(1);
			r.Add(m_degree+1);
		}
		break;
	default:
		break;
	}
}

double CBSplineDoc::Basic_Spline(int j,int k,double u)  
/*
功能：给定u值，输出次数为k，序号为j的B样条函数对应的值（根据书本p239的方法推导）
参数说明：j表示下标值，k表示次数，j应该取值在index-k和index之间 共k+1个值
*/
{
	double Value1,Value2,alphl1,alphl2,BSt;
	int index = Knot_Index(u);
	if(j < index - k || j > index) return 0.0;

	if(j == index && k == 0) return 1.0;
	else
	{
		Value1 = Knot_Value(j+k) - Knot_Value(j);
		Value2 = Knot_Value(j+k+1) - Knot_Value(j+1);
		if(fabs(Value1) < 1e-5) alphl1 = 0;
		else alphl1 = (u - Knot_Value(j))/Value1;
		if(fabs(Value2) < 1e-5) alphl2 = 0;
		else alphl2 = (Knot_Value(j+k+1) - u)/Value2;	
		BSt = Basic_Spline(j,k-1,u) * alphl1 + Basic_Spline(j+1,k-1,u) * alphl2;
		return BSt;
	}
}

double CBSplineDoc::Length(int j)  //计算两个控制点之间的距离
{
	double lx,ly;
	lx = m_ptArray[j - 1].x - m_ptArray[j].x;
	ly = m_ptArray[j - 1].y - m_ptArray[j].y;
	return sqrt(pow(lx,2) + pow(ly,2));
}