OpenGL由已知控制点绘制模拟曲面地形

本人原创，欢迎转载，转载请注明出处http://www.cnblogs.com/zhouchanwen

主要要点：1.将离散的数据点网格化曲面

　　　　　2.对3d模型的鼠标控制，如虚拟球的实现

由已知控制点通过曲面拟合方法，将不规则的数据分布转换成规则的网格分布，然后绘制三维地形曲面图。即如图所示：



数据规则网格化（简称网格化）。网格化实际是一种曲面拟合方法。关于曲面拟合算法有很多，下面我们采用曲面样条方法实现网格化。我们定义曲面样条函数：



下面为c++/c代码：

3D_MAP.cpp

#include <iostream>
#include "tools.h"
#include "Grids.h"
#include "3dmap.h"

using namespace std;

// 3D_MAP.cpp : 定义控制台应用程序的入口点。

int xFar=0.0f,yFar=0.0f,zFar=0.0f;
int wWidth=1300,wHeight=700;
int oldX,oldY;
bool gIsStartTrackBall = false;
bool gIsMoveMap=false;
TrackBall trackball;
_3dMap map;

void displayEvent()
{
glClearColor(0, 0, 0.1f, 1);
glClear(GL_COLOR_BUFFER_BIT|GL_DEPTH_BUFFER_BIT);

glMatrixMode(GL_MODELVIEW);
glLoadIdentity();
glTranslatef(xFar, yFar,zFar);

trackball.makeRolate();
map.drawMap();

glFlush();
glutSwapBuffers();  // 交换缓存
}

void mouseMoveEvent(int x, int y)
{
if(gIsStartTrackBall)
{
trackball.MouseMove(x,y);
glutPostRedisplay();
}
if(gIsMoveMap)
{
xFar-=oldX-x;
yFar+=oldY-y;
oldX=x;
oldY=y;
glutPostRedisplay();
}
}
// 鼠标事件函数
void mouseEvent(int button, int state, int x, int y)
{
if(button == GLUT_LEFT_BUTTON)
{
if(state == GLUT_DOWN)
{
oldX=x;
oldY=y;
trackball.setXY(x,y);
gIsStartTrackBall = true;
}
else if(state == GLUT_UP)
{
oldX=x;
oldY=y;
gIsStartTrackBall = false;
}
glutPostRedisplay();
}else if(button == GLUT_RIGHT_BUTTON)
{
if(state == GLUT_DOWN)
{
oldX=x;
oldY=y;
gIsMoveMap=true;
}else if(state == GLUT_UP)
{
oldX=x;
oldY=y;
gIsMoveMap = false;
}
}
}
// 窗体尺寸变化事件
void resizeEvent(int w, int h)
{
wWidth=w;
wHeight=h;
zFar=0.0f;
xFar=0.0f;
yFar=0.0f;
glViewport(0,0, w, h);
glMatrixMode(GL_PROJECTION);
glLoadIdentity();
h = h > 0 ? h : 1;
float aspect = (float)w / (float)h;
gluPerspective(45,aspect,1.0,1500.0);
glTranslatef(0,0,-300.0f);

trackball.resize();

glutPostRedisplay();
}
void processSpecialKeys(int key, int x, int y) {
if(key==101)
{
zFar+=4;
glutPostRedisplay();
}
if(key==103)
{//
zFar-=4;
glutPostRedisplay();
}
printf("key:%d\n",key);
}

void MenuFunc(int data)
{
switch(data)
{
case 1:
map.setLineOrFill();break;
default:break;
}
glutPostRedisplay();
}
void glInit()
{
glShadeModel(GL_FLAT);//SMOOTH//GL_FLAT
glClearColor(1.0f, 1.0f, 1.0f, 1.0f);
glClearDepth(1.0f);

glEnable(GL_NORMALIZE);

glEnable ( GL_DEPTH_TEST );
glAlphaFunc(GL_GREATER,0);
glDepthFunc(GL_LEQUAL);
glHint(GL_PERSPECTIVE_CORRECTION_HINT, GL_NICEST);

glEnable( GL_BLEND);
glBlendFunc( GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);

glEnable(GL_POINT_SMOOTH);
glEnable(GL_LINE_SMOOTH);
glEnable (GL_POLYGON_SMOOTH);

glHint(GL_POINT_SMOOTH_HINT, GL_NICEST); // Make round points, not square points
glHint(GL_LINE_SMOOTH_HINT, GL_NICEST);  // Antialias the lines
glHint (GL_POLYGON_SMOOTH_HINT, GL_NICEST);

//glClearColor(1.0,1.0,1.0,1.0);  //窗口背景设置为白色
glMatrixMode(GL_MODELVIEW); //设置投影参数

glEnable( GL_COLOR_MATERIAL) ;
glColorMaterial(GL_FRONT_AND_BACK,GL_AMBIENT_AND_DIFFUSE);

}

int main(int argc,char* argv[])
{
map.initMap();

glutInit(&argc,argv);                                             //初始化GLUT
glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB|GLUT_DEPTH | GLUT_MULTISAMPLE);  //设置显示模式
glutInitWindowPosition(0,0);  //设置显示窗口的左上角位置
glutInitWindowSize(wWidth,wHeight);         //设置窗口的长和高
glutCreateWindow("3DMap" );     //创造显示窗口

glInit();                                       //开始初始化过程
glutReshapeFunc(resizeEvent);
glutDisplayFunc(displayEvent);
glutMouseFunc(mouseEvent);
glutSpecialFunc(processSpecialKeys);
glutMotionFunc(mouseMoveEvent);
glutCreateMenu(MenuFunc);
glutAddMenuEntry("填充/网格",1);
glutAttachMenu(GLUT_MIDDLE_BUTTON);

glutMainLoop();    //显示所有并等候

getchar();
return 0;
}

tools.h

tools.h中定义了TrackBall类，用于对3D模型进行鼠标的旋转，世界中心点的移动等控制

#include <cmath>
#include <gl/glut.h>

using namespace std;

struct GPoint3d{
double mX, mY, mZ;
double x() { return mX; }
double y() { return mY; }
double z() {return mZ; }
void setX(double x)  { mX = x; }
void setY(double y)  { mY = y; }
void setZ(double z)     { mZ = z; }
void set(double x, double y, double z) { mX = x; mY = y; mZ = z;}
};

class TrackBall
{
int OldX;
int OldY;
double mMatrix[16];
public:
TrackBall(){}
// 向量的点积
double dotMult(GPoint3d v1, GPoint3d v2);
// 向量的叉积
GPoint3d crossMult(GPoint3d v1, GPoint3d v2);
// 将鼠标二维点映射为球面向量(用于鼠标追踪球)
GPoint3d gMousePtToSphereVec(int x, int y, int w, int h);
void makeRolate();

void MouseMove(int x, int y);
void resize()
{
glGetDoublev(GL_MODELVIEW_MATRIX, mMatrix);  // 返回当前模型矩阵
}
void setXY(int x, int y){OldX = x;OldY = y;}
void setP(double *v)//{x1,y1,z1,  x2,y2,z2,  x3,y3,z3}
{
GPoint3d v1,v2,v3;
v1.setX(v[3]-v[0]);
v1.setY(v[4]-v[1]);
v1.setZ(v[5]-v[2]);

v2.setX(v[6]-v[0]);
v2.setY(v[7]-v[1]);
v2.setZ(v[8]-v[2]);

v3=crossMult(v1,v2);

glNormal3f(v3.x(),v3.y(),v3.z());
}
};

double TrackBall::dotMult(GPoint3d v1, GPoint3d v2)
{
double angle;
angle = v1.x()*v2.x()+v1.y()*v2.y()+v1.z()*v2.z();
return angle;
}
GPoint3d TrackBall::crossMult(GPoint3d v1, GPoint3d v2)
{
GPoint3d v;
v.setX(v1.y()*v2.z()-v1.z()*v2.y());
v.setY(v1.z()*v2.x()-v1.x()*v2.z());
v.setZ(v1.x()*v2.y()-v1.y()*v2.x());
return v;
}

// 将鼠标二维点映射为球面向量(用于鼠标追踪球)
GPoint3d TrackBall::gMousePtToSphereVec(int x, int y, int w, int h)
{
double x1,y1,z1,r,len;
GPoint3d vec;
x1 = (2.0*x - w) / w;
y1 = (h - 2.0*y) / h;
r=x1*x1+y1*y1;
if(r > 1) r = 1;
z1 = sqrt(1 - r);
len = sqrt(x1*x1+y1*y1+z1*z1);
vec.setX(x1/len);
vec.setY(y1/len);
vec.setZ(z1/len);
return vec;
}
void TrackBall::makeRolate()
{
glMultMatrixd(mMatrix);
}
void TrackBall::MouseMove(int x, int y)
{
if(x != OldX || y != OldY)
{
int wWidth,wHeight;
wWidth = glutGet(GLUT_WINDOW_WIDTH);
wHeight = glutGet(GLUT_WINDOW_HEIGHT);
GPoint3d lastVec = gMousePtToSphereVec(OldX, OldY, wWidth, wHeight);
GPoint3d currentVec = gMousePtToSphereVec(x, y, wWidth, wHeight);
OldX = x;        OldY = y;
// 求旋转角度
double rotAngle = acos(dotMult(lastVec,currentVec))*57.29577958;
// 求旋转向量轴
GPoint3d axis = crossMult(lastVec,currentVec);
glMatrixMode(GL_MODELVIEW);
glPushMatrix();
glLoadIdentity();
glRotated(rotAngle, axis.x(), axis.y(), axis.z()); // 旋转

glMultMatrixd(mMatrix);
glGetDoublev(GL_MODELVIEW_MATRIX, mMatrix);  // 返回当前模型矩阵

glPopMatrix();
}
}

接下来是3dmap.h，其中定义了_3dMap类，主要用来完成3d模型的主要绘制任务

class _3dMap
{
double **data;
int M,N;
int dip;
double scale;
bool showBaseLine;
double max,min;
bool LineMode;
public :
_3dMap()
{
M=320;
N=220;
showBaseLine=true;
LineMode=false;
scale=20.0;
dip=1;
max=-1000;
min=10000;
}
~_3dMap(){}
void initMap();
void drawBaseLine();
void drawMap();
void _3dMap::setLineOrFill()
{
LineMode=!LineMode;
}
double getAver(double * arr)
{
double num=0;
int n=3;//(int)(sizeof(arr)/sizeof(double))/3;
for(int i=0;i<n;i++)
num+=arr[3*i+2];
return num/n;
}

void setColor(double z1,double z2,double z3)
{
float r,g,b;
double temp=(max+min)/2;
double aver=(z1+z2+z3)/3;
/*printf("%lf  ",aver);*/
if(aver>temp)
{
r=(aver-temp)/(temp-min);
b=0;
}else {
r=0;
b=(temp-aver)/(temp-min);
}
g=1-((abs(temp-aver))/(temp-min));
glColor3f(r,g,b);
}
};

void _3dMap::initMap()
{

int temp[]={
229,219,199,216,235,255,266,285,272,241,246,281,284,275,261,273,
221,214,195,216,234,258,273,289,281,249,259,278,287,272,275,277,
213,203,196,206,221,232,259,293,294,277,258,285,287,283,288,286,
204,195,200,201,209,218,231,259,288,306,286,291,301,311,319,298,
196,207,201,211,239,234,241,259,294,315,317,321,325,322,325,341,
208,218,204,214,235,260,239,268,298,291,331,313,281,280,280,280,
216,231,218,196,220,255,271,253,264,303,322,312,276,243,238,239,
236,242,218,198,200,215,224,238,261,294,324,312,280,255,220,200,
255,241,219,211,206,225,252,275,284,285,305,316,271,237,208,191,
245,218,207,198,214,241,261,256,273,276,291,298,281,238,197,175,
225,215,205,195,208,221,235,252,262,271,301,275,245,212,181,171
};
int len=11*16;
int i,j;
double** src=new double*[3];
for(i=0;i<3;i++)
src[i]=new double[len];
for(i=0;i<11;i++)
{//y
for(j=0;j<16;j++)
{//x
src[0][i*16+j]=j*10;
src[1][i*16+j]=i*10;
src[2][i*16+j]=temp[i*16+j];
}
}
Grids g1(3,len,src);
double** des=new double*[2];
for(i=0;i<2;i++)
des[i]=new double[N*M];
for(i=0;i<N;i++)
{//y
for(j=0;j<M;j++)
{//x
des[0][i*M+j]=j/2;
des[1][i*M+j]=i/2;
}
}
Grids g2(2,N*M,des);
data=g1.Surface_Fitting(&g1,&g2);

for(i=0;i<N;i++)
{//height
for(int j=0;j<M;j++)
{//width
data[0][i*M+j]*=2;
data[1][i*M+j]*=2;
data[2][i*M+j]/=5;

if(max<data[2][i*M+j])
max=data[2][i*M+j];
if(min>data[2][i*M+j])
min=data[2][i*M+j];
}
}
printf("max=%lf,min=%lf\n",max,min);
}
void _3dMap::drawBaseLine()
{
glBegin(GL_LINES);
glColor3f(1, 0, 0);
glVertex2i(0,0);//x line
glVertex2i(100,0);

glColor3f(0, 1, 0);
glVertex2i(0,0);//y line
glVertex2i(0,100);

glColor3f(0, 0, 1);
glVertex3f(0,0,0);//
glVertex3f(0,0,500);
glEnd();
}
void _3dMap::drawMap()
{

if(showBaseLine)
drawBaseLine();

if(LineMode)
glPolygonMode(GL_FRONT_AND_BACK,GL_LINE);
else
glPolygonMode(GL_FRONT_AND_BACK,GL_FILL);

glLineWidth(0.4f);
for(int i=0;i<N-1;i++)
{//height22

for(int j=0;j<M-1;j++)
{//width32
glBegin(GL_TRIANGLE_FAN);
int t=i*M+j;
setColor(data[2][t],data[2][t+1],data[2][t+1+M]);

glVertex3f(data[0][t],data[1][t],data[2][t]);
glVertex3f(data[0][t+1],data[1][t+1],data[2][t+1]);
glVertex3f(data[0][t+1+M],data[1][t+1+M],data[2][t+1+M]);
setColor(data[2][t+M],data[2][t+1],data[2][t+1+M]);
glVertex3f(data[0][t+M],data[1][t+M],data[2][t+M]);
glEnd();
}
}
}

下面是Grids.h文件，定义了Grids类，用于对控制点的数据进行网格化

#include <iostream>

#include <cmath>

using namespace std;

class Grids
{
int row;
int col;
double **mat;
public:
Grids();
Grids(int r, int c);
Grids(int r, int c, double **m);
~Grids();
void Mul_Mat(Grids* m1, Grids* m2, Grids* m3);
void Inv_Mat(Grids* m1, Grids* m2);
double** Surface_Fitting(Grids* m1, Grids* m2);
};

Grids::Grids()
{
row = 0;
col = 0;
mat = 0;
}

Grids::Grids(int r, int c)
{
row = r;
col = c;
mat = 0;
}

Grids::Grids(int r, int c, double **m)
{
row = r;
col = c;
mat = m;
}

Grids::~Grids()
{
mat=0;
row=0;
col=0;
}

void Grids::Mul_Mat(Grids* m1, Grids* m2, Grids* m3)
{
int i=0,j=0,p=0;
double sum=0;
if (m1->col != m2->row) {
cout<<"\n行、列数不匹配!";
exit(0);
}
m3->row=m1->row;
m3->col=m2->col;
m3->mat=new double*[m1->row];
if (NULL==m3->mat) {
cout<<"ERROR!\n";
exit(0);
}
for (i=0;i<m3->row;i++) {
m3->mat[i]=new double[m3->col];
for (j=0;j<m3->col;j++) {
for (m3->mat[i][j]=0,p=0;p<m1->col;p++) {
m3->mat[i][j]+=m1->mat[i][p]*m2->mat[p][j];
}
}
}
}

void Grids::Inv_Mat(Grids* m1, Grids* m2)
{
int i,j,n,*is,*js,k;
double d,p;
if(m1->row!=m1->col) {
cout<<"ERROR! 必须是方阵才能求逆!\n";
exit(1);
}
m2->mat=new double*[m1->row];//申请行指针数组
if(NULL==m2->mat){
cout<<"ERROR! 申请内存出错!\n";
exit(1);
}
for (i=0;i<m1->row;i++) {
m2->mat[i]=new double[m1->col];//申请行
for (j=0;j<m1->col;j++)
m2->mat[i][j]=m1->mat[i][j];
}
n=m1->row;
m2->row=m1->row;
m2->col=m1->col;
is=new int
;
js=new int
;
if (NULL==is || NULL==js){
cout<<"ERROR! 申请内存出错!\n";
exit(1);
}
for (k=0;k<=n-1;k++){ //全选主元
d=0.000;
for (i=k;i<=n-1;i++){
for (j=k;j<=n-1;j++){
p=fabs(m2->mat[i][j]);
if (p>d){
d=p;
is[k]=i;
js[k]=j;
}
}
}
if(1.0==d+1.0){
delete []is;
delete []js;
cout<<"ERROR ! 矩阵求逆出错!\n";
exit(1);
}
if(is[k]!=k){ /*行交换*/
for (j=0;j<=n-1;j++){
p=m2->mat[k][j];
m2->mat[k][j]=m2->mat[is[k]][j];
m2->mat[is[k]][j]=p;
}
}
if(js[k] != k) { /*列交换*/
for (i=0;i<=n-1;i++) {
p=m2->mat[i][k];
m2->mat[i][k]=m2->mat[i][js[k]];
m2->mat[i][js[k]]=p;
}
}
m2->mat[k][k]=1/m2->mat[k][k];
for (j=0;j<=n-1;j++){
if (j != k){
m2->mat[k][j]=m2->mat[k][j]*m2->mat[k][k];
}
}
for (i=0;i<=n-1;i++){
if(i!=k){
for (j=0;j<=n-1;j++){
if(j!=k){
m2->mat[i][j]=m2->mat[i][j]-m2->mat[i][k]*m2->mat[k][j];
}
}
}
}
for (i=0;i<=n-1;i++){
if (i!=k){
m2->mat[i][k]=-m2->mat[i][k]*m2->mat[k][k];
}
}
}
for (k=n-1;k>=0;k--){
if (js[k]!=k){
for (j=0;j<=n-1;j++){
p=m2->mat[k][j];
m2->mat[k][j]=m2->mat[js[k]][j];
m2->mat[js[k]][j]=p;
}
}
if (is[k] != k){
for (i=0;i<=n-1;i++){
p=m2->mat[i][k];
m2->mat[i][k]=m2->mat[i][is[k]];
m2->mat[i][is[k]]=p;
}
}
}
delete []is;
delete []js;
}

double** Grids::Surface_Fitting(Grids* m1, Grids* m2)
{  // m1为已知数组R,m2为自变量数组
int i,j;
double ipso=0.05;
Grids *R=new Grids(m1->col+3,m1->col+3);
R->mat=new double*[R->row];
for(i=0;i<R->row;i++)
R->mat[i]=new double[R->col];
for(i=0;i<R->row;i++){
for(j=0;j<R->col;j++){
if(i<m1->col&&j<m1->col){
if(i==j)
R->mat[i][j]=ipso;
else
R->mat[i][j]=sqrt(pow(m1->mat[0][i] - m1->mat[0][j],2)+pow(m1->mat[1][i] - m1->mat[1][j],2));
}
if(i>=m1->col&&j<m1->col){
if(i==m1->col)
R->mat[i][j]=1;
else
R->mat[i][j]=m1->mat[i-m1->col-1][j];
}
if(i<m1->col&&j>=m1->col){
if(j==m1->col)
R->mat[i][j]=1;
else
R->mat[i][j]=m1->mat[j-(m1->col)-1][i];
}
if(i>=m1->col&&j>=m1->col){
R->mat[i][j]=0;
}
}
}

Grids *z1=new Grids(m1->col+3,1);
Grids *invm1=new Grids();
Grids *F=new Grids();
z1->mat=new double*[z1->row];
for(i=0;i<z1->row;i++){
z1->mat[i]=new double[z1->col];
for(j=0;j<z1->col;j++){
if(i<m1->col)
z1->mat[i][j]=m1->mat[2][i];
else
z1->mat[i][j]=0;
}
}
Inv_Mat(R, invm1);

Mul_Mat(invm1, z1, F);

Grids *z=new Grids();
Grids *r=new Grids(m2->col,F->row);
r->mat=new double*[r->row];
for(i=0;i<r->row;i++)
r->mat[i]=new double[r->col];
for(i=0;i<r->row;i++){
for(j=0;j<m1->col+3;j++){
if(j<m1->col){
double temp=sqrt(pow(m2->mat[0][i] - m1->mat[0][j],2)+pow(m2->mat[1][i] - m1->mat[1][j],2));
r->mat[i][j]=temp*log(temp+ipso);
}
else{
if(j==m1->col)
r->mat[i][j]=1;
else
r->mat[i][j]=m2->mat[j-m1->col-1][i];
}
}
}
/*    FF->mat=new double*[FF->row];
for(i=0;i<FF->row;i++)
FF->mat[i]=new double[FF->col];
for(i=0;i<FF->row;i++){
for(j=0;j<FF->col;j++){
if(i<m2->col)
FF->mat[i][j]=F->mat[i/4][j/4];
else
FF->mat[i][j]=F->mat[i-m2->col+m1->col][j];
}
}*/
Mul_Mat(r, F, z);

Grids *xyz=new Grids(m2->row+1,m2->col);
xyz->mat=new double*[xyz->row];
for(i=0;i<xyz->row;i++)
xyz->mat[i]=new double[xyz->col];
for(i=0;i<xyz->row;i++){
for(j=0;j<xyz->col;j++){
if(i<m2->row)
xyz->mat[i][j]=m2->mat[i][j];
else
xyz->mat[i][j]=z->mat[j][0];
}
}
for(i=0;i<R->row;i++)
delete [] R->mat[i];
delete [] R->mat;
for(i=0;i<z1->row;i++)
delete [] z1->mat[i];
delete [] z1->mat;
for(i=0;i<invm1->row;i++)
delete [] invm1->mat[i];
delete [] invm1->mat;
for(i=0;i<F->row;i++)
delete [] F->mat[i];
delete [] F->mat;
for(i=0;i<z->row;i++)
delete [] z->mat[i];
delete [] z->mat;
for(i=0;i<r->row;i++)
delete [] r->mat[i];
delete [] r->mat;
return xyz->mat;
}

下面是最终效果：





操作：使用鼠标左键进行虚拟球控制模型，键盘上下键进行拉近拉远视图，鼠标右键移动视图中心点，点击鼠标滚轮弹出菜单，菜单中可以设置地图绘制方式，网格状和填充。

本人原创，欢迎转载，转载请注明出处http://www.cnblogs.com/zhouchanwen