一、B样条曲面定义
二、双三次均匀B样条曲面片
2.1 理论
已知曲面控制顶点Pij, i=0,1,2,3; j= 0,1,2,3,即m=n=3,参数(u, v) [0, 1)x[0, 1], p=q=3。依次用线段连接点Pij阵列中的相邻两点,所形成的空间网格称为控制网格。由m=n=3构成了4×4=16个顶点的控制网格,其相应的曲面称为双三次B样条曲面。
具体曲面定义表达式,如下:
2.2 代码实现
2.2.1 定义曲面类
#pragma once
#include"P3.h"
class CBicubicUniformBSplinePatch
{
public:
CBicubicUniformBSplinePatch(void);
~CBicubicUniformBSplinePatch(void);
void ReadControlPoint(CP3 P[4][4]);//读入16个控制点
void DrawCurvedPatch(CDC* pDC);//绘制双三次均匀B样条曲面片
void DrawControlGrid(CDC* pDC);//绘制控制网格
private:
void LeftMultiplyMatrix(double M[][4],CP3 P[][4]);//左乘顶点矩阵
void RightMultiplyMatrix(CP3 P[][4],double M[][4]);//右乘顶点矩阵
void TransposeMatrix(double M[][4]);//转置矩阵
CP2 ObliqueProjection(CP3 Point3);//斜二测投影
public:
CP3 P[4][4];//三维控制点
};
#include "StdAfx.h"
#include "BicubicUniformBSplinePatch.h"
#include<math.h>
#define ROUND(d) int(d+0.5)
CBicubicUniformBSplinePatch::CBicubicUniformBSplinePatch(void)
{
}
CBicubicUniformBSplinePatch::~CBicubicUniformBSplinePatch(void)
{
}
void CBicubicUniformBSplinePatch::ReadControlPoint(CP3 P[4][4])
{
for (int i=0;i<4;i++)
for(int j=0;j<4;j++)
this->P[i][j]=P[i][j];
}
void CBicubicUniformBSplinePatch::DrawCurvedPatch(CDC* pDC)
{
CPen NewPen, *pOldPen;
NewPen.CreatePen(PS_SOLID, 1, RGB(255, 0, 0));
pOldPen = pDC->SelectObject(&NewPen);
double M[4][4];//系数矩阵Mb
M[0][0]=-1;M[0][1]=3; M[0][2]=-3;M[0][3]=1;
M[1][0]=3; M[1][1]=-6;M[1][2]=3; M[1][3]=0;
M[2][0]=-3;M[2][1]=0; M[2][2]=3; M[2][3]=0;
M[3][0]=1; M[3][1]=4; M[3][2]=1; M[3][3]=0;
CP3 P3[4][4];//曲线计算用控制点数组
for(int i=0;i<4;i++)
for(int j=0;j<4;j++)
P3[i][j]=P[i][j];
LeftMultiplyMatrix(M,P3);//数字矩阵左乘三维点矩阵
TransposeMatrix(M);//计算转置矩阵
RightMultiplyMatrix(P3,M);//数字矩阵右乘三维点矩阵
double Step=0.1;//步长
double u0,u1,u2,u3,v0,v1,v2,v3;//u,v参数的幂
for(double u=0;u<=1;u+=Step)
for(double v=0;v<=1;v+=Step)
{
u3=u*u*u;u2=u*u;u1=u;u0=1;v3=v*v*v;v2=v*v;v1=v;v0=1;
CP3 pt=(u3*P3[0][0]+u2*P3[1][0]+u1*P3[2][0]+u0*P3[3][0])*v3
+(u3*P3[0][1]+u2*P3[1][1]+u1*P3[2][1]+u0*P3[3][1])*v2
+(u3*P3[0][2]+u2*P3[1][2]+u1*P3[2][2]+u0*P3[3][2])*v1
+(u3*P3[0][3]+u2*P3[1][3]+u1*P3[2][3]+u0*P3[3][3])*v0;
pt=pt/36;
CP2 Point2=ObliqueProjection(pt);//斜投影
if(v==0)
pDC->MoveTo(ROUND(Point2.x),ROUND(Point2.y));
else
pDC->LineTo(ROUND(Point2.x),ROUND(Point2.y));
}
for(double v=0;v<=1;v+=Step)
for(double u=0;u<=1;u+=Step)
{
u3=u*u*u;u2=u*u;u1=u;u0=1;v3=v*v*v;v2=v*v;v1=v;v0=1;
CP3 pt=(u3*P3[0][0]+u2*P3[1][0]+u1*P3[2][0]+u0*P3[3][0])*v3
+(u3*P3[0][1]+u2*P3[1][1]+u1*P3[2][1]+u0*P3[3][1])*v2
+(u3*P3[0][2]+u2*P3[1][2]+u1*P3[2][2]+u0*P3[3][2])*v1
+(u3*P3[0][3]+u2*P3[1][3]+u1*P3[2][3]+u0*P3[3][3])*v0;
pt=pt/36;
CP2 Point2=ObliqueProjection(pt);//斜投影
if(0==u)
pDC->MoveTo(ROUND(Point2.x),ROUND(Point2.y));
else
pDC->LineTo(ROUND(Point2.x),ROUND(Point2.y));
}
pDC->SelectObject(pOldPen);
NewPen.DeleteObject();
}
void CBicubicUniformBSplinePatch::LeftMultiplyMatrix(double M[][4],CP3 P[][4])//左乘矩阵M*P
{
CP3 T[4][4];//临时矩阵
for(int i=0;i<4;i++)
for(int j=0;j<4;j++)
T[i][j]=M[i][0]*P[0][j]+M[i][1]*P[1][j]+M[i][2]*P[2][j]+M[i][3]*P[3][j];
for(int i=0;i<4;i++)
for(int j=0;j<4;j++)
P[i][j]=T[i][j];
}
void CBicubicUniformBSplinePatch::RightMultiplyMatrix(CP3 P[][4],double M[][4])//右乘矩阵P*M
{
CP3 T[4][4];//临时矩阵
for(int i=0;i<4;i++)
for(int j=0;j<4;j++)
T[i][j]=P[i][0]*M[0][j]+P[i][1]*M[1][j]+P[i][2]*M[2][j]+P[i][3]*M[3][j];
for(int i=0;i<4;i++)
for(int j=0;j<4;j++)
P[i][j]=T[i][j];
}
void CBicubicUniformBSplinePatch::TransposeMatrix(double M[][4])//转置矩阵
{
double T[4][4];//临时矩阵
for(int i=0;i<4;i++)
for(int j=0;j<4;j++)
T[j][i]=M[i][j];
for(int i=0;i<4;i++)
for(int j=0;j<4;j++)
M[i][j]=T[i][j];
}
CP2 CBicubicUniformBSplinePatch::ObliqueProjection(CP3 Point3)//斜二测投影
{
CP2 Point2;
Point2.x=Point3.x-Point3.z*sqrt(2.0)/2.0;
Point2.y=Point3.y-Point3.z*sqrt(2.0)/2.0;
return Point2;
}
void CBicubicUniformBSplinePatch::DrawControlGrid(CDC *pDC)//绘制控制网格
{
CP2 P2[4][4];//二维控制点
for(int i=0;i<4;i++)
for(int j=0;j<4;j++)
P2[i][j]=ObliqueProjection(P[i][j]);
CPen NewPen,*pOldPen;
NewPen.CreatePen(PS_SOLID,2,RGB(0,0,255));
pOldPen=pDC->SelectObject(&NewPen);
for(int i=0;i<4;i++)
{
pDC->MoveTo(ROUND(P2[i][0].x),ROUND(P2[i][0].y));
for(int j=1;j<4;j++)
pDC->LineTo(ROUND(P2[i][j].x),ROUND(P2[i][j].y));
}
for(int j=0;j<4;j++)
{
pDC->MoveTo(ROUND(P2[0][j].x),ROUND(P2[0][j].y));
for(int i=1;i<4;i++)
pDC->LineTo(ROUND(P2[i][j].x),ROUND(P2[i][j].y));
}
pDC->SelectObject(pOldPen);
NewPen.DeleteObject();
}
2.2.2 读取16个控制点
CP3 P[4][4];//控制点
void CGeometricfiguretestView::ReadPoint()//读入控制多边形顶点
{
P[0][0].x=20; P[0][0].y=0; P[0][0].z=200;
P[0][1].x=0; P[0][1].y=100;P[0][1].z=150;
P[0][2].x=-130;P[0][2].y=100;P[0][2].z=50;
P[0][3].x=-250;P[0][3].y=50; P[0][3].z=0;
P[1][0].x=100; P[1][0].y=100;P[1][0].z=150;
P[1][1].x= 30; P[1][1].y=100;P[1][1].z=100;
P[1][2].x=-40; P[1][2].y=100;P[1][2].z=50;
P[1][3].x=-110;P[1][3].y=100;P[1][3].z=0;
P[2][0].x=280; P[2][0].y=90; P[2][0].z=140;
P[2][1].x=110; P[2][1].y=120;P[2][1].z=80;
P[2][2].x=0; P[2][2].y=130;P[2][2].z=30;
P[2][3].x=-100;P[2][3].y=150;P[2][3].z=-50;
P[3][0].x=350; P[3][0].y=30; P[3][0].z=150;
P[3][1].x=200; P[3][1].y=150;P[3][1].z=50;
P[3][2].x=50; P[3][2].y=200;P[3][2].z=0;
P[3][3].x=0; P[3][3].y=100;P[3][3].z =-70;
}
2.2.3 调用绘制
void CGeometricfiguretestView::OnDraw(CDC* pDC)
{
CTestDoc* pDoc = GetDocument();
ASSERT_VALID(pDoc);
if (!pDoc)
return;
// TODO: add draw code for native data here
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);
DrawGraph(pDC);//绘制图形
}
void CGeometricfiguretestView::DrawGraph(CDC* pDC)//绘制图形
{
CBicubicUniformBSplinePatch patch;
patch.ReadControlPoint(P);
patch.DrawCurvedPatch(pDC);
patch.DrawControlGrid(pDC);
}
编译,运行。
三、非均匀B样条曲面
3.1 理论
3.2 代码实现
3.2.1 定义曲面类
#pragma once
#include"P3.h"
class CBSplineSurface
{
public:
CBSplineSurface(void);
~CBSplineSurface(void);
void Initialize(CP3** ptr,int m,int p,int n,int q);//初始化
void GetKnotVector(double* T,int nCount,int num,int order, BOOL bU);//获取节点矢量
void DrawCurvedSurface(CDC* pDC);//绘制曲面
double BasisFunctionValue(double u,int i,int k,double* T);//根据节点矢量T,计算基函数
void DrawControlGrid(CDC* pDC);//绘制控制网格
CP2 ObliqueProjection(CP3 Point3);//斜二测投影
public:
int m,p;//u向的顶点数减1与次数
int n,q;//v向的顶点数减1与次数
double **V;//u向节点矢量数组
double **U;//v向节点矢量数组
CP3 **P;//三维控制点
};
#include "StdAfx.h"
#include "BSPlineSurface.h"
#include<math.h>
#define ROUND(d) int(d+0.5)
#define PI 3.1415926
CBSplineSurface::CBSplineSurface(void)
{
P = NULL;
U = NULL;
V = NULL;
}
CBSplineSurface::~CBSplineSurface(void)
{
if(NULL != P)
{
for(int i=0;i<n+1;i++)
{
delete []P[i];
P[i]=NULL;
}
delete []P;
P=NULL;
}
if(NULL!=U)
{
for(int i=0;i<n+1;i++)
{
delete []U[i];
U[i]=NULL;
}
delete []U;
U=NULL;
}
if(NULL!=V)
{
for(int i=0;i<m+1;i++)
{
delete []V[i];
V[i]=NULL;
}
delete []V;
V=NULL;
}
}
void CBSplineSurface::Initialize(CP3** ptr,int n,int q,int m,int p)//初始化
{
P=new CP3*[n+1];//建立动态二维数组
for(int i=0;i<n+1;i++)
P[i]=new CP3[m+1];
for(int i=0;i<n+1;i++)//二维数组初始化
for(int j=0;j<m+1;j++)
P[i][j]=ptr[i][j];
this->m=m,this->p=p;
this->n=n,this->q=q;
U=new double*[n+1];//建立u向节点矢量动态数组
for(int i=0;i<n+1;i++)
U[i]=new double[m+p+2];
V=new double*[m+1];//建立v向节点矢量动态数组
for(int i=0;i<m+1;i++)
V[i]=new double[n+q+2];
}
void CBSplineSurface::DrawCurvedSurface(CDC* pDC)//绘制曲面
{
for(int i=0;i<m+1;i++)
GetKnotVector(V[i],i,n,q, FALSE);//获取V
for(int i=0;i<n+1;i++)
GetKnotVector(U[i],i,m,p, TRUE);//获取U
CPen NewPen(PS_SOLID,1,RGB(255,0,0));//曲线颜色
CPen* pOldPen=pDC->SelectObject(&NewPen);
double Step=0.1;//步长
for(double u=0.0;u<=1.0;u+=Step)//绘制v向曲线
{
for(double v=0.0;v<=1.0;v+=Step)
{
CP3 pt(0,0,0);
for(int i=0;i<n+1;i++)
{
for(int j=0;j<m+1;j++)
{
double BValueU=BasisFunctionValue(u,j,p,U[i]);
double BValueV=BasisFunctionValue(v,i,q,V[j]);
pt=pt+P[i][j]*BValueU*BValueV;
}
}
CP2 Point2=ObliqueProjection(pt);//斜投影
if (v==0.0)
pDC->MoveTo(ROUND(Point2.x),ROUND(Point2.y));
else
pDC->LineTo(ROUND(Point2.x),ROUND(Point2.y));
}
}
for(double v=0.0;v<=1.0;v+=Step)//绘制u向曲线
{
for(double u=0.0;u<=1.0;u+=Step)
{
CP3 pt(0,0,0);
for(int i=0;i<n+1;i++)
{
for(int j=0;j<m+1;j++)
{
double BValueU=BasisFunctionValue(u,j,p,U[i]);
double BValueV=BasisFunctionValue(v,i,q,V[j]);
pt=pt+P[i][j]*BValueU*BValueV;
}
}
CP2 Point2=ObliqueProjection(pt);//斜投影
if (u==0.0)
pDC->MoveTo(ROUND(Point2.x),ROUND(Point2.y));
else
pDC->LineTo(ROUND(Point2.x),ROUND(Point2.y));
}
}
pDC->SelectObject(pOldPen);
NewPen.DeleteObject();
}
void CBSplineSurface::DrawControlGrid(CDC* pDC)//绘制控制网格
{
CP2** P2=new CP2*[n+1];
for(int i=0;i<n+1;i++)
P2[i]=new CP2[m+1];
for(int i=0;i<n+1;i++)
for(int j=0;j<m+1;j++)
P2[i][j]=ObliqueProjection(P[i][j]);
CPen NewPen,*pOldPen;
NewPen.CreatePen(PS_SOLID,2,RGB(0,0,255));
pOldPen=pDC->SelectObject(&NewPen);
for(int i=0;i<n+1;i++)
{
pDC->MoveTo(ROUND(P2[i][0].x),ROUND(P2[i][0].y));
for(int j=1;j<m+1;j++)
pDC->LineTo(ROUND(P2[i][j].x),ROUND(P2[i][j].y));
}
for(int j=0;j<m+1;j++)
{
pDC->MoveTo(ROUND(P2[0][j].x),ROUND(P2[0][j].y));
for(int i=1;i<n+1;i++)
pDC->LineTo(ROUND(P2[i][j].x),ROUND(P2[i][j].y));
}
pDC->SelectObject(pOldPen);
NewPen.DeleteObject();
if(NULL!=P2)
{
for(int i=0;i<n+1;i++)
{
delete []P2[i];
P2[i] = NULL;
}
delete []P2;
P2 = NULL;
}
}
void CBSplineSurface::GetKnotVector(double* T,int nCount,int num,int order, BOOL bU)//Hartley-Judd算法获取节点矢量数组
{
for(int i=0;i<=order;i++) //小于曲线次数order的节点值为0
T[i]=0.0;
for(int i=num+1;i<=num+order+1;i++)//大于控制点个数num的节点值为1
T[i]=1.0;
//计算num-order个内节点
for(int i=order+1;i<=num;i++)
{
double sum=0.0;
for(int j=order+1;j<=i;j++)
{
double numerator=0.0;
for(int loop=j-order;loop<=j-1;loop++)
{
if(bU)
numerator+=(P[nCount][loop].x-P[nCount][loop-1].x)*(P[nCount][loop].x-P[nCount][loop-1].x)+(P[nCount][loop].y-P[nCount][loop-1].y)*(P[nCount][loop].y-P[nCount][loop-1].y);
else
numerator+=(P[loop][nCount].x-P[loop-1][nCount].x)*(P[loop][nCount].x-P[loop-1][nCount].x)+(P[loop][nCount].y-P[loop-1][nCount].y)*(P[loop][nCount].y-P[loop-1][nCount].y);
}
double denominator=0.0;//计算分母
for(int loop1=order+1;loop1<=num+1;loop1++)
{
for(int loop2=loop1-order;loop2<=loop1-1;loop2++)
{
if(bU)
denominator+=(P[nCount][loop2].x-P[nCount][loop2-1].x)*(P[nCount][loop2].x-P[nCount][loop2-1].x)+(P[nCount][loop2].y-P[nCount][loop2-1].y)*(P[nCount][loop2].y-P[nCount][loop2-1].y);
else
denominator+=(P[loop2][nCount].x-P[loop2-1][nCount].x)*(P[loop2][nCount].x-P[loop2-1][nCount].x)+(P[loop2][nCount].y-P[loop2-1][nCount].y)*(P[loop2][nCount].y-P[loop2-1][nCount].y);
}
}
sum+=numerator/denominator;
}
T[i]=sum;
}
}
double CBSplineSurface::BasisFunctionValue(double t,int i,int order,double* T)//计算B样条基函数
{
double value1,value2,value;
if(order==0)
{
if(t>=T[i]&&t<T[i+1])
return 1.0;
else
return 0.0;
}
if(order>0)
{
if(t<T[i]||t>T[i+order+1])
return 0.0;
else
{
double coffcient1,coffcient2;//凸组合系数1,凸组合系数2
double denominator=0.0;//分母
denominator=T[i+order]-T[i];//递推公式第一项分母
if(denominator==0.0)//约定0/0
coffcient1=0.0;
else
coffcient1=(t-T[i])/denominator;
denominator=T[i+order+1]-T[i+1];//递推公式第二项分母
if(0.0==denominator)//约定0/0
coffcient2=0.0;
else
coffcient2=(T[i+order+1]-t)/denominator;
value1=coffcient1*BasisFunctionValue(t,i,order-1,T);//递推公式第一项的值
value2=coffcient2*BasisFunctionValue(t,i+1,order-1,T);//递推公式第二项的值
value=value1+value2;//基函数的值
}
}
return value;
}
CP2 CBSplineSurface::ObliqueProjection(CP3 Point3)//斜二测投影
{
CP2 Point2;
Point2.x=Point3.x-Point3.z*sqrt(2.0)/2.0;
Point2.y=Point3.y-Point3.z*sqrt(2.0)/2.0;
return Point2;
}
3.2.2 初始化数据
CP3** P3;//二维控制网格
int m,n;//u向和v向的控制点数-1
int p,q;//u向和v向的次数
CBSplineSurface BSurface;//B样条曲面对象
CGeometricfiguretestView::CGeometricfiguretestView()
{
// TODO: add construction code here
m=4, p=3;//u向顶点数5个,曲线次数为3
n=3, q=2;//v向顶点数4个,曲线次数为2
P3=new CP3*[n+1];//建立动态二维数组
for(int i=0;i<n+1;i++)
P3[i]=new CP3[m+1];
P3[0][0]=CP3(20,0,200), P3[0][1]=CP3(100,100,150),P3[0][2]=CP3(230,130,140), P3[0][3]=CP3(350,100,150),P3[0][4]=CP3(450,30,150);
P3[1][0]=CP3(0,100,150), P3[1][1]=CP3(30,100,120), P3[1][2]=CP3(110,120,80), P3[1][3]=CP3(200,150,50), P3[1][4]=CP3(300,150,50);
P3[2][0]=CP3(-130,100,50),P3[2][1]=CP3(-40,120,50), P3[2][2]=CP3(30,150,30), P3[2][3]=CP3(50,150,0), P3[2][4]=CP3(150,200,0);
P3[3][0]=CP3(-250,50,0), P3[3][1]=CP3(-150,100,0), P3[3][2]=CP3(-100,150,-50), P3[3][3]=CP3(0,150,-70), P3[3][4]=CP3(80,220,-70);
BSurface.Initialize(P3,n,q,m,p);
}
3.2.3 调用绘制
void CGeometricfiguretestView::OnDraw(CDC* pDC)
{
CTestDoc* pDoc = GetDocument();
ASSERT_VALID(pDoc);
if (!pDoc)
return;
// TODO: add draw code for native data here
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);
DrawGraph(pDC);//绘制图形
}
void CGeometricfiguretestView::DrawGraph(CDC* pDC)//绘制B样条曲面
{
BSurface.DrawCurvedSurface(pDC);
BSurface.DrawControlGrid(pDC);
}
编译,运行。
© 版权声明
文章版权归作者所有,未经允许请勿转载。
THE END
暂无评论内容