cad.cpp

用参数样条曲线拟合型值点的算法。实现将11个型值点的拟和。在CAD中有广泛应用。用OpenGL绘制出图形。

////////////////////////////////////////////////////////////////////////////////////////
//                                                                                    //
//                                                                                    //
//                                                                                    //
//                                                                                    //
//    Microsoft Visual C++ 6.0 , Windows XP 编译通过                                  //
//                                                                                    //
//    注：编译该程序时须将 opengl32.lib glu32.lib glaux.lib 加入工程中                //
//                                                                                    //
//    菜单选取操作project->setting->link->general->object/library modules，           //
//    在object/library modules编辑框中加入opengl32.lib glu32.lib glaux.lib 。         //
//                                                                                    //
////////////////////////////////////////////////////////////////////////////////////////

#include <windows.h>
#include <GL/gl.h>
#include <GL/glu.h>
#include <GL/glaux.h>
#include <math.h>

#define STEP 0.01                 //步长
#define PI 3.14159265359          //π
#define K 15                      //图形放大倍数


///////////////////////////////////////////////////////////////////////////////////////
//                                                                                   //
//                              用参数样条曲线拟合型值点                             //
//                                                                                   //
///////////////////////////////////////////////////////////////////////////////////////

struct PolarPoint
{//极坐标点类型
	float rho;
	float theta;
};

PolarPoint P[11];              //存放型值点极坐标
float m[11];                   //m
float h[11];                   //h
float lambda[10];              //λ
float mu[10];                  //μ
float c[10];                   //c

///////////////////////////////////////////////////////////////////////////////
//用追赶法解三对角型线性方程组
//a[i],b[i],c[i]分别为三条对角线元素，n为方程个数，x为解向量

void Pursue(float *a, float *b, float *c, float *d, const int n ,float *x)
{
	int k;
	float* r = new float[n];
	float* y = new float[n];

	r[1] = c[1] / b[1];
	y[1] = d[1] / b[1];

	for(k=2; k<=n-1; ++k)
	{
		r[k] = c[k] / (b[k] - r[k-1]*a[k]);
		y[k] = ( (d[k] - y[k-1]*a[k]) / (b[k] - r[k-1]*a[k]) );
	}
	y[n] = ( (d[n] - y[n-1]*a[n]) / (b[n] - r[n-1]*a[n]) );

	x[n] = y[n];
	for(k=n-1; k>=1; --k)
	{
		x[k] = y[k] - r[k]*x[k+1];
	}
		

//	delete []r;       
//	delete []y;
}

///////////////////////////////////////////////////////////////////////////////
//三次样条插值
void Spline()
{
	int i;
	
	for(i=0; i<=10; ++i)                            //输入11个型值点的极坐标
	{
		P[i].rho = 1.0;
		P[i].theta = i * 10.0;
	}  
	
	m[0] = 0;                                       //已知条件
	m[9] = 0;                                       //已知条件

	for(i=1; i<=10; ++i)                            //计算h[i]
	{
		h[i] = P[i].theta - P[i-1].theta;
	}

	for(i=1; i<=9; ++i)                            //计算lambda[i]及mu[i]
	{
		lambda[i] = h[i+1] / ( h[i] + h[i+1] );

		mu[i] = 1 - lambda[i];
	}

	for(i=1; i<=9; ++i)                            //计算c[i]
	{
		c[i] = 3 * ( lambda[i]*(P[i].rho - P[i-1].rho)/h[i] 
			   + mu[i]*(P[i+1].rho - P[i].rho)/h[i+1] );
	}


	float TWO[8];
	for(i=0; i<8; ++i)
	{
		TWO[i] = 2.0;
	}

	Pursue(lambda, TWO, mu, c, 8, m);            //追赶法解m[1]～m[8]

	m[10] = ( c[9] - lambda[9]*m[8] - 2*m[9] ) / mu[9];    //计算m[10]

}

////////////////////////////////////////////////////////////////////////////////////
//Hermite基函数

float F0(float u)
{
	return ( 2*u*u*u - 3*u*u + 1 );
}

float F1(float u)
{
	return ( -2*u*u*u + 3*u*u );
}

float G0(float u)
{
	return ( u*(u-1)*(u-1) );
}

float G1(float u)
{
	return ( u*u*(u-1) );
}

////////////////////////////////////////////////////////////////////////////////////////
//参数样条曲线方程
float SplineFunction(float theta)
{
	float rho, u; 

	int i = int(theta / 10) + 1;        //根据theta值判断此点应用哪一段方程计算

	u = (theta - P[i-1].theta) / h[i];
	
	rho = P[i-1].rho * F0(u) + P[i].rho * F1(u)
		  +h[i] * ( m[i-1] * G0(u) + m[i] * G1(u) );

	return rho;
}



///////////////////////////////////////////////////////////////////////////////////
//                                                                               //
//                               用OpenGL绘制曲线                                //
//                                                                               //
///////////////////////////////////////////////////////////////////////////////////

void myinit(void);
void DrawMyObjects(void);
void CALLBACK myReshape(GLsizei w,GLsizei h);
void CALLBACK display(void);

////////////////////////////////////////////////////////////////////////////////////

void myinit(void)
{
	glClearColor(0.0,0.0,0.0,0.0);
	glClear(GL_COLOR_BUFFER_BIT);
	glShadeModel(GL_FLAT);
}

void CALLBACK myReshape(GLsizei w,GLsizei h)
{
	glViewport(0,0,w,h);
    glMatrixMode(GL_PROJECTION);
    glLoadIdentity();
	
	if(w<=h)
		glOrtho(-20.0,20.0,-20.0*(GLfloat)h/(GLfloat)w,
		20.0*(GLfloat)h/(GLfloat)w,-50.0,50.0);
    else
		glOrtho(-20.0*(GLfloat)w/(GLfloat)h,
		20.0*(GLfloat)w/(GLfloat)h,-20.0,20.0,-50.0,50.0);
	
    glMatrixMode(GL_MODELVIEW);
    glLoadIdentity();
}

void CALLBACK display(void)
{//绘制和显示图形
	glColor3f(1.0, 0.0, 0.0);
	DrawMyObjects();
	glFlush();
}

void DrawMyObjects(void)
{//绘制图形
	float theta;
	int i;

	//绘制z轴(横轴)
	glBegin(GL_LINE_STRIP);
    glColor3f(0.0, 0.0, 1.0);                     //设置为蓝色
    glVertex2f(-200.0, 0.0);
    glVertex2f(200.0, 0.0);
    glEnd();

	//绘制y轴(纵轴)
	glBegin(GL_LINE_STRIP);
    glColor3f(0.0, 0.0, 1.0);                     //设置为蓝色
    glVertex2f(0.0, -200.0);
    glVertex2f(0.0, 200.0);
    glEnd();

	//绘制参数样条曲线
	glBegin(GL_LINE_STRIP);                       //用折线连接点
	glColor3f(1.0, 0.0, 0.0);                     //设置为红色
	for(theta=0.0; theta<=100.0; theta+=STEP)
	{
		glVertex2f( K*SplineFunction(theta)*sin(theta*PI/180), 
		            K*SplineFunction(theta)*cos(theta*PI/180) );
	}
	glEnd();

	//标记型值点
	glBegin(GL_POINTS);
    glColor3f(1.0, 1.0, 0.0);                     //设置为黄色
    for(i=0; i<11; ++i)
	{
		glVertex2f( K*P[i].rho * sin(P[i].theta*PI/180), 
			K*P[i].rho * cos(P[i].theta*PI/180) );
	}
	
    glEnd();
}

//////////////////////////////////////////////////////////////////////////////////////

void main(void)
{
	Spline();               //对型值点进行三次参数样条插值
	
	auxInitDisplayMode(AUX_SINGLE|AUX_RGBA);	
	auxInitPosition(0,0,640,480);
	auxInitWindow("三次参数样条插值");
	myinit();
	auxReshapeFunc(myReshape);
	auxMainLoop(display);
}