主要要點：1.將離散的資料點網格化曲面

　　　　　2.對3d模型的滑鼠控制，如虛擬球的實現

由已知控制點通過曲面擬合方法，將不規則的資料分佈轉換成規則的網格分佈，然後繪製三維地形曲面圖。即如圖所示：

資料規則網格化（簡稱網格化）。網格化實際是一種曲面擬合方法。關於曲面擬合算法有很多，下面我們採用曲面樣條方法實現網格化。我們定義曲面樣條函式：

下面為c++/c程式碼：

3D_MAP.cpp

1 #include <iostream>
  2 #include "tools.h"
  3 #include "Grids.h"
  4 #include "3dmap.h"
  5 
  6 using namespace std;
  7 
  8 // 3D_MAP.cpp : 定義控制檯應用程式的入口點。
  9 
 10 int xFar=0.0f,yFar=0.0f,zFar=0.0f;
 11 int wWidth=1300,wHeight=700;
 12 int oldX,oldY;
 13 bool gIsStartTrackBall = false;
 14 bool gIsMoveMap=false;
 15 TrackBall trackball;
 16 _3dMap map;
 17 
 18 void displayEvent()
 19 {
 20     glClearColor(0, 0, 0.1f, 1);
 21     glClear(GL_COLOR_BUFFER_BIT|GL_DEPTH_BUFFER_BIT);
 22     
 23     glMatrixMode(GL_MODELVIEW);
 24     glLoadIdentity();
 25     glTranslatef(xFar, yFar,zFar);
 26     
 27     trackball.makeRolate();
 28     map.drawMap();
 29 
 30     glFlush();
 31     glutSwapBuffers();  // 交換快取 
 32 }
 33 
 34 void mouseMoveEvent(int x, int y)
 35 {
 36     if(gIsStartTrackBall)
 37     {
 38         trackball.MouseMove(x,y);
 39         glutPostRedisplay();
 40     }
 41     if(gIsMoveMap)
 42     {
 43         xFar-=oldX-x;
 44         yFar+=oldY-y;
 45         oldX=x;
 46         oldY=y;
 47         glutPostRedisplay();
 48     }
 49 }
 50 // 滑鼠事件函式 
 51 void mouseEvent(int button, int state, int x, int y)
 52 {
 53     if(button == GLUT_LEFT_BUTTON)
 54     {
 55         if(state == GLUT_DOWN)
 56         {
 57             oldX=x;
 58             oldY=y;
 59             trackball.setXY(x,y);
 60             gIsStartTrackBall = true;
 61         }
 62         else if(state == GLUT_UP)
 63         {  
 64             oldX=x;
 65             oldY=y;
 66             gIsStartTrackBall = false;
 67         }            
 68         glutPostRedisplay();
 69     }else if(button == GLUT_RIGHT_BUTTON)
 70     {
 71         if(state == GLUT_DOWN)
 72         {
 73             oldX=x;
 74             oldY=y;
 75             gIsMoveMap=true;
 76         }else if(state == GLUT_UP)
 77         {
 78             oldX=x;
 79             oldY=y;
 80             gIsMoveMap = false;
 81         }
 82     }
 83 }       
 84 // 窗體尺寸變化事件 
 85 void resizeEvent(int w, int h)
 86 {
 87     wWidth=w;
 88     wHeight=h;
 89     zFar=0.0f;
 90     xFar=0.0f;
 91     yFar=0.0f;
 92     glViewport(0,0, w, h);
 93     glMatrixMode(GL_PROJECTION);
 94     glLoadIdentity();
 95     h = h > 0 ? h : 1;
 96     float aspect = (float)w / (float)h;     
 97     gluPerspective(45,aspect,1.0,1500.0);
 98     glTranslatef(0,0,-300.0f);
 99 
100     trackball.resize();
101 
102     
103     glutPostRedisplay();
104 }   
105 void processSpecialKeys(int key, int x, int y) { 
106     if(key==101)
107     {
108         zFar+=4;
109         glutPostRedisplay();
110     }
111     if(key==103)
112     {//
113         zFar-=4;
114         glutPostRedisplay();
115     }
116     printf("key:%d\n",key);
117 } 
118 
119 void MenuFunc(int data)
120 {
121     switch(data)
122     {
123     case 1:
124         map.setLineOrFill();break;
125     default:break;
126     }
127     glutPostRedisplay();
128 }
129 void glInit()
130 {
131     glShadeModel(GL_FLAT);//SMOOTH//GL_FLAT
132     glClearColor(1.0f, 1.0f, 1.0f, 1.0f);
133     glClearDepth(1.0f);
134 
135     glEnable(GL_NORMALIZE);
136 
137     glEnable ( GL_DEPTH_TEST );
138     glAlphaFunc(GL_GREATER,0);
139     glDepthFunc(GL_LEQUAL);
140     glHint(GL_PERSPECTIVE_CORRECTION_HINT, GL_NICEST);
141 
142     glEnable( GL_BLEND); 
143     glBlendFunc( GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
144 
145     glEnable(GL_POINT_SMOOTH);
146     glEnable(GL_LINE_SMOOTH);
147     glEnable (GL_POLYGON_SMOOTH);
148 
149     glHint(GL_POINT_SMOOTH_HINT, GL_NICEST); // Make round points, not square points
150     glHint(GL_LINE_SMOOTH_HINT, GL_NICEST);  // Antialias the lines
151     glHint (GL_POLYGON_SMOOTH_HINT, GL_NICEST);
152 
153     
154     //glClearColor(1.0,1.0,1.0,1.0);  //視窗背景設定為白色
155     glMatrixMode(GL_MODELVIEW); //設定投影引數
156 
157     glEnable( GL_COLOR_MATERIAL) ;
158     glColorMaterial(GL_FRONT_AND_BACK,GL_AMBIENT_AND_DIFFUSE);
159 
160 }
161 
162 int main(int argc,char* argv[])
163 {
164     map.initMap();
165 
166     glutInit(&argc,argv);                                             //初始化GLUT
167     glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB|GLUT_DEPTH | GLUT_MULTISAMPLE);  //設定顯示模式
168     glutInitWindowPosition(0,0);  //設定顯示視窗的左上角位置
169     glutInitWindowSize(wWidth,wHeight);         //設定視窗的長和高
170     glutCreateWindow("3DMap" );     //創造顯示視窗
171 
172     glInit();                                       //開始初始化過程
173     glutReshapeFunc(resizeEvent);
174     glutDisplayFunc(displayEvent);
175     glutMouseFunc(mouseEvent);
176     glutSpecialFunc(processSpecialKeys);
177     glutMotionFunc(mouseMoveEvent);
178     glutCreateMenu(MenuFunc);
179     glutAddMenuEntry("填充/網格",1);
180     glutAttachMenu(GLUT_MIDDLE_BUTTON);
181 
182     glutMainLoop();    //顯示所有並等候
183 
184     getchar();
185     return 0;
186 }
 

tools.h

tools.h中定義了TrackBall類，用於對3D模型進行滑鼠的旋轉，世界中心點的移動等控制

  1 #include <cmath>
  2 #include <gl/glut.h>
  3 
  4 using namespace std;
  5 
  6 struct GPoint3d{
  7     double mX, mY, mZ;
  8     double x() { return mX; }
  9     double y() { return mY; }
 10     double z() {return mZ; } 
 11     void setX(double x)  { mX = x; } 
 12     void setY(double y)  { mY = y; }
 13     void setZ(double z)     { mZ = z; } 
 14     void set(double x, double y, double z) { mX = x; mY = y; mZ = z;}
 15 };
 16 
 17 class TrackBall
 18 {
 19     int OldX;
 20     int OldY;
 21     double mMatrix[16];
 22 public:    
 23     TrackBall(){}
 24         // 向量的點積
 25     double dotMult(GPoint3d v1, GPoint3d v2);
 26     // 向量的叉積
 27     GPoint3d crossMult(GPoint3d v1, GPoint3d v2);
 28         // 將滑鼠二維點對映為球面向量(用於滑鼠追蹤球) 
 29     GPoint3d gMousePtToSphereVec(int x, int y, int w, int h);
 30     void makeRolate();
 31 
 32     void MouseMove(int x, int y);
 33     void resize()
 34     {
 35         glGetDoublev(GL_MODELVIEW_MATRIX, mMatrix);  // 返回當前模型矩陣
 36     }
 37     void setXY(int x, int y){OldX = x;OldY = y;}
 38     void setP(double *v)//{x1,y1,z1,  x2,y2,z2,  x3,y3,z3}
 39     {
 40         GPoint3d v1,v2,v3;
 41         v1.setX(v[3]-v[0]);
 42         v1.setY(v[4]-v[1]);
 43         v1.setZ(v[5]-v[2]);
 44 
 45         v2.setX(v[6]-v[0]);
 46         v2.setY(v[7]-v[1]);
 47         v2.setZ(v[8]-v[2]);
 48 
 49         v3=crossMult(v1,v2);
 50 
 51         glNormal3f(v3.x(),v3.y(),v3.z());
 52     }
 53 };
 54 
 55 double TrackBall::dotMult(GPoint3d v1, GPoint3d v2)
 56 {
 57     double angle;
 58     angle = v1.x()*v2.x()+v1.y()*v2.y()+v1.z()*v2.z();
 59     return angle;
 60 }
 61 GPoint3d TrackBall::crossMult(GPoint3d v1, GPoint3d v2)
 62 {
 63     GPoint3d v;
 64     v.setX(v1.y()*v2.z()-v1.z()*v2.y());
 65     v.setY(v1.z()*v2.x()-v1.x()*v2.z());
 66     v.setZ(v1.x()*v2.y()-v1.y()*v2.x());
 67     return v;
 68 }
 69 
 70         // 將滑鼠二維點對映為球面向量(用於滑鼠追蹤球) 
 71 GPoint3d TrackBall::gMousePtToSphereVec(int x, int y, int w, int h)
 72 {
 73     double x1,y1,z1,r,len;
 74     GPoint3d vec;
 75     x1 = (2.0*x - w) / w;
 76     y1 = (h - 2.0*y) / h;
 77     r=x1*x1+y1*y1;
 78     if(r > 1) r = 1;
 79     z1 = sqrt(1 - r);
 80     len = sqrt(x1*x1+y1*y1+z1*z1);
 81     vec.setX(x1/len);
 82     vec.setY(y1/len);
 83     vec.setZ(z1/len);
 84     return vec;
 85 }
 86 void TrackBall::makeRolate()
 87 {
 88     glMultMatrixd(mMatrix);
 89 }
 90 void TrackBall::MouseMove(int x, int y)
 91 {
 92     if(x != OldX || y != OldY)
 93     {
 94         int wWidth,wHeight;
 95         wWidth = glutGet(GLUT_WINDOW_WIDTH);
 96         wHeight = glutGet(GLUT_WINDOW_HEIGHT);
 97         GPoint3d lastVec = gMousePtToSphereVec(OldX, OldY, wWidth, wHeight);
 98         GPoint3d currentVec = gMousePtToSphereVec(x, y, wWidth, wHeight);
 99         OldX = x;        OldY = y;
100                 // 求旋轉角度
101         double rotAngle = acos(dotMult(lastVec,currentVec))*57.29577958;
102         // 求旋轉向量軸
103         GPoint3d axis = crossMult(lastVec,currentVec);
104         glMatrixMode(GL_MODELVIEW);
105         glPushMatrix();
106         glLoadIdentity();
107         glRotated(rotAngle, axis.x(), axis.y(), axis.z()); // 旋轉
108             
109         glMultMatrixd(mMatrix);
110         glGetDoublev(GL_MODELVIEW_MATRIX, mMatrix);  // 返回當前模型矩陣
111             
112         glPopMatrix();
113     }
114 }
 

接下來是3dmap.h，其中定義了_3dMap類，主要用來完成3d模型的主要繪製任務

1 class _3dMap
  2 {
  3     double **data;
  4     int M,N;
  5     int dip;
  6     double scale;
  7     bool showBaseLine;
  8     double max,min;
  9     bool LineMode;
 10 public :
 11     _3dMap()
 12     {
 13         M=320;
 14         N=220;
 15         showBaseLine=true;
 16         LineMode=false;
 17         scale=20.0;
 18         dip=1;
 19         max=-1000;
 20         min=10000;
 21     }
 22     ~_3dMap(){}
 23     void initMap();
 24     void drawBaseLine();
 25     void drawMap();
 26     void _3dMap::setLineOrFill()
 27     {    
 28         LineMode=!LineMode;
 29     }
 30     double getAver(double * arr)
 31     {
 32         double num=0;
 33         int n=3;//(int)(sizeof(arr)/sizeof(double))/3;
 34         for(int i=0;i<n;i++)
 35             num+=arr[3*i+2];
 36         return num/n;
 37     }
 38     
 39     void setColor(double z1,double z2,double z3)
 40     {
 41         float r,g,b;
 42         double temp=(max+min)/2;
 43         double aver=(z1+z2+z3)/3;
 44         /*printf("%lf  ",aver);*/
 45         if(aver>temp)
 46         {
 47             r=(aver-temp)/(temp-min);
 48             b=0;
 49         }else {
 50             r=0;
 51             b=(temp-aver)/(temp-min);
 52         }
 53         g=1-((abs(temp-aver))/(temp-min));
 54         glColor3f(r,g,b);
 55     }
 56 };
 57 
 58 void _3dMap::initMap()
 59 {
 60 
 61     int temp[]={
 62         229,219,199,216,235,255,266,285,272,241,246,281,284,275,261,273,
 63         221,214,195,216,234,258,273,289,281,249,259,278,287,272,275,277,
 64         213,203,196,206,221,232,259,293,294,277,258,285,287,283,288,286,
 65         204,195,200,201,209,218,231,259,288,306,286,291,301,311,319,298,
 66         196,207,201,211,239,234,241,259,294,315,317,321,325,322,325,341,
 67         208,218,204,214,235,260,239,268,298,291,331,313,281,280,280,280,
 68         216,231,218,196,220,255,271,253,264,303,322,312,276,243,238,239,
 69         236,242,218,198,200,215,224,238,261,294,324,312,280,255,220,200,
 70         255,241,219,211,206,225,252,275,284,285,305,316,271,237,208,191,
 71         245,218,207,198,214,241,261,256,273,276,291,298,281,238,197,175,
 72         225,215,205,195,208,221,235,252,262,271,301,275,245,212,181,171
 73     };
 74     int len=11*16;
 75     int i,j;
 76     double** src=new double*[3];
 77     for(i=0;i<3;i++)
 78         src[i]=new double[len];    
 79     for(i=0;i<11;i++)
 80     {//y
 81         for(j=0;j<16;j++)
 82         {//x
 83             src[0][i*16+j]=j*10;
 84             src[1][i*16+j]=i*10;
 85             src[2][i*16+j]=temp[i*16+j];
 86         }
 87     }
 88     Grids g1(3,len,src);
 89     double** des=new double*[2];
 90     for(i=0;i<2;i++)
 91         des[i]=new double[N*M];    
 92     for(i=0;i<N;i++)
 93     {//y
 94         for(j=0;j<M;j++)
 95         {//x
 96             des[0][i*M+j]=j/2;
 97             des[1][i*M+j]=i/2;
 98         }
 99     }
100     Grids g2(2,N*M,des);
101     data=g1.Surface_Fitting(&g1,&g2);
102 
103     for(i=0;i<N;i++)
104     {//height
105         for(int j=0;j<M;j++)
106         {//width
107             data[0][i*M+j]*=2;
108             data[1][i*M+j]*=2;
109             data[2][i*M+j]/=5;
110 
111             if(max<data[2][i*M+j])
112                 max=data[2][i*M+j];
113             if(min>data[2][i*M+j])
114                 min=data[2][i*M+j];
115         }
116     }
117     printf("max=%lf,min=%lf\n",max,min);
118 }
119 void _3dMap::drawBaseLine()
120 {
121         glBegin(GL_LINES); 
122             glColor3f(1, 0, 0);
123             glVertex2i(0,0);//x line
124             glVertex2i(100,0);
125 
126             glColor3f(0, 1, 0);
127             glVertex2i(0,0);//y line
128             glVertex2i(0,100);
129 
130             glColor3f(0, 0, 1);
131             glVertex3f(0,0,0);//
132             glVertex3f(0,0,500);
133         glEnd();  
134 }
135 void _3dMap::drawMap()
136 {
137 
138     if(showBaseLine)
139         drawBaseLine();
140 
141     if(LineMode)
142         glPolygonMode(GL_FRONT_AND_BACK,GL_LINE);
143     else 
144         glPolygonMode(GL_FRONT_AND_BACK,GL_FILL);
145     
146     glLineWidth(0.4f);
147     for(int i=0;i<N-1;i++)
148     {//height22
149         
150         for(int j=0;j<M-1;j++)
151         {//width32
152             glBegin(GL_TRIANGLE_FAN);
153                 int t=i*M+j;
154                 setColor(data[2][t],data[2][t+1],data[2][t+1+M]);
155                 
156                 glVertex3f(data[0][t],data[1][t],data[2][t]);
157                 glVertex3f(data[0][t+1],data[1][t+1],data[2][t+1]);
158                 glVertex3f(data[0][t+1+M],data[1][t+1+M],data[2][t+1+M]);
159                 setColor(data[2][t+M],data[2][t+1],data[2][t+1+M]);
160                 glVertex3f(data[0][t+M],data[1][t+M],data[2][t+M]);
161             glEnd();
162         }
163     }
164 }
 

下面是Grids.h檔案，定義了Grids類，用於對控制點的資料進行網格化

 1 #include <iostream>
  2 
  3 #include <cmath>
  4 
  5 using namespace std;
  6 
  7 class Grids
  8 {
  9     int row;
 10     int col;
 11     double **mat;
 12 public:
 13     Grids();
 14     Grids(int r, int c);
 15     Grids(int r, int c, double **m);
 16     ~Grids();
 17     void Mul_Mat(Grids* m1, Grids* m2, Grids* m3);
 18     void Inv_Mat(Grids* m1, Grids* m2);
 19     double** Surface_Fitting(Grids* m1, Grids* m2);
 20 };
 21 
 22 Grids::Grids()
 23 {
 24     row = 0;
 25     col = 0;
 26     mat = 0;
 27 }
 28 
 29 Grids::Grids(int r, int c)
 30 {
 31     row = r;
 32     col = c;
 33     mat = 0;
 34 }
 35 
 36 Grids::Grids(int r, int c, double **m)
 37 {
 38     row = r;
 39     col = c;
 40     mat = m;
 41 }
 42 
 43 Grids::~Grids()
 44 {
 45     mat=0;
 46     row=0;
 47     col=0;
 48 }
 49 
 50 void Grids::Mul_Mat(Grids* m1, Grids* m2, Grids* m3)
 51 {
 52     int i=0,j=0,p=0;
 53     double sum=0;
 54     if (m1->col != m2->row) {
 55         cout<<"\n行、列數不匹配!";
 56         exit(0);
 57     }
 58     m3->row=m1->row;
 59     m3->col=m2->col;
 60     m3->mat=new double*[m1->row];
 61     if (NULL==m3->mat) {
 62         cout<<"ERROR!\n";
 63         exit(0);
 64     }
 65     for (i=0;i<m3->row;i++) {
 66         m3->mat[i]=new double[m3->col];
 67         for (j=0;j<m3->col;j++) {
 68             for (m3->mat[i][j]=0,p=0;p<m1->col;p++) {
 69                 m3->mat[i][j]+=m1->mat[i][p]*m2->mat[p][j];
 70             }
 71         }
 72     }
 73 }
 74 
 75 void Grids::Inv_Mat(Grids* m1, Grids* m2)
 76 {
 77     int i,j,n,*is,*js,k;
 78     double d,p;
 79     if(m1->row!=m1->col) {
 80         cout<<"ERROR! 必須是方陣才能求逆!\n"; 
 81         exit(1); 
 82     }
 83     m2->mat=new double*[m1->row];//申請行指標陣列
 84     if(NULL==m2->mat){
 85         cout<<"ERROR! 申請記憶體出錯!\n";
 86         exit(1);
 87     }
 88     for (i=0;i<m1->row;i++) {
 89         m2->mat[i]=new double[m1->col];//申請行
 90         for (j=0;j<m1->col;j++)
 91             m2->mat[i][j]=m1->mat[i][j];
 92     }
 93     n=m1->row;
 94     m2->row=m1->row;
 95     m2->col=m1->col;
 96     is=new int[n];
 97     js=new int[n];
 98     if (NULL==is || NULL==js){
 99         cout<<"ERROR! 申請記憶體出錯!\n";
100         exit(1);
101     }
102     for (k=0;k<=n-1;k++){ //全選主元
103         d=0.000;
104         for (i=k;i<=n-1;i++){
105             for (j=k;j<=n-1;j++){
106                 p=fabs(m2->mat[i][j]);
107                 if (p>d){
108                     d=p;
109                     is[k]=i;
110                     js[k]=j;
111                 }
112             }
113         }
114         if(1.0==d+1.0){
115             delete []is;
116             delete []js;
117             cout<<"ERROR ! 矩陣求逆出錯!\n";
118             exit(1);
119         }
120         if(is[k]!=k){ /*行交換*/
121             for (j=0;j<=n-1;j++){
122                 p=m2->mat[k][j];
123                 m2->mat[k][j]=m2->mat[is[k]][j];
124                 m2->mat[is[k]][j]=p;
125             }
126         }
127         if(js[k] != k) { /*列交換*/
128             for (i=0;i<=n-1;i++) {
129                 p=m2->mat[i][k];
130                 m2->mat[i][k]=m2->mat[i][js[k]];
131                 m2->mat[i][js[k]]=p;
132             }
133         }
134         m2->mat[k][k]=1/m2->mat[k][k];
135         for (j=0;j<=n-1;j++){
136             if (j != k){
137                 m2->mat[k][j]=m2->mat[k][j]*m2->mat[k][k];
138             }
139         }
140         for (i=0;i<=n-1;i++){
141             if(i!=k){
142                 for (j=0;j<=n-1;j++){
143                     if(j!=k){
144                         m2->mat[i][j]=m2->mat[i][j]-m2->mat[i][k]*m2->mat[k][j];
145                     }
146                 }
147             }
148         }
149         for (i=0;i<=n-1;i++){
150             if (i!=k){
151                 m2->mat[i][k]=-m2->mat[i][k]*m2->mat[k][k];
152             }
153         }
154     }
155     for (k=n-1;k>=0;k--){
156         if (js[k]!=k){
157             for (j=0;j<=n-1;j++){
158                 p=m2->mat[k][j];
159                 m2->mat[k][j]=m2->mat[js[k]][j];
160                 m2->mat[js[k]][j]=p;
161             }
162         }
163         if (is[k] != k){
164             for (i=0;i<=n-1;i++){
165                 p=m2->mat[i][k];
166                 m2->mat[i][k]=m2->mat[i][is[k]];
167                 m2->mat[i][is[k]]=p;
168             }
169         }
170     }
171     delete []is;
172     delete []js;
173 }
174 
175 double** Grids::Surface_Fitting(Grids* m1, Grids* m2)
176 {  // m1為已知陣列R,m2為自變數陣列
177     int i,j;
178     double ipso=0.05;
179     Grids *R=new Grids(m1->col+3,m1->col+3);
180     R->mat=new double*[R->row];
181     for(i=0;i<R->row;i++)
182         R->mat[i]=new double[R->col];
183     for(i=0;i<R->row;i++){
184         for(j=0;j<R->col;j++){
185             if(i<m1->col&&j<m1->col){
186                 if(i==j)
187                     R->mat[i][j]=ipso;
188                 else
189                     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));
190             }
191             if(i>=m1->col&&j<m1->col){
192                 if(i==m1->col)
193                     R->mat[i][j]=1;
194                 else
195                     R->mat[i][j]=m1->mat[i-m1->col-1][j];
196             }
197             if(i<m1->col&&j>=m1->col){
198                 if(j==m1->col)
199                     R->mat[i][j]=1;
200                 else
201                     R->mat[i][j]=m1->mat[j-(m1->col)-1][i];
202             }
203             if(i>=m1->col&&j>=m1->col){
204                 R->mat[i][j]=0;
205             }
206         }
207     }
208 
209     Grids *z1=new Grids(m1->col+3,1);
210     Grids *invm1=new Grids();
211     Grids *F=new Grids();
212     z1->mat=new double*[z1->row];
213     for(i=0;i<z1->row;i++){
214         z1->mat[i]=new double[z1->col];
215         for(j=0;j<z1->col;j++){
216             if(i<m1->col)
217                 z1->mat[i][j]=m1->mat[2][i];
218             else
219                 z1->mat[i][j]=0;
220         }
221     }
222     Inv_Mat(R, invm1);
223     
224     Mul_Mat(invm1, z1, F);
225 
226     Grids *z=new Grids();
227     Grids *r=new Grids(m2->col,F->row);
228     r->mat=new double*[r->row];
229     for(i=0;i<r->row;i++)
230         r->mat[i]=new double[r->col];
231     for(i=0;i<r->row;i++){
232         for(j=0;j<m1->col+3;j++){
233             if(j<m1->col){
234                 double temp=sqrt(pow(m2->mat[0][i] - m1->mat[0][j],2)+pow(m2->mat[1][i] - m1->mat[1][j],2));
235                 r->mat[i][j]=temp*log(temp+ipso);
236             }
237             else{
238                 if(j==m1->col)
239                     r->mat[i][j]=1;
240                 else
241                     r->mat[i][j]=m2->mat[j-m1->col-1][i];
242             }
243         }
244     }
245 /*    FF->mat=new double*[FF->row];
246     for(i=0;i<FF->row;i++)
247         FF->mat[i]=new double[FF->col];
248     for(i=0;i<FF->row;i++){
249         for(j=0;j<FF->col;j++){
250             if(i<m2->col)
251                 FF->mat[i][j]=F->mat[i/4][j/4];
252             else
253                 FF->mat[i][j]=F->mat[i-m2->col+m1->col][j];
254         }
255     }*/
256     Mul_Mat(r, F, z);
257 
258     Grids *xyz=new Grids(m2->row+1,m2->col);
259     xyz->mat=new double*[xyz->row];
260     for(i=0;i<xyz->row;i++)
261         xyz->mat[i]=new double[xyz->col];
262     for(i=0;i<xyz->row;i++){
263         for(j=0;j<xyz->col;j++){
264             if(i<m2->row)
265                 xyz->mat[i][j]=m2->mat[i][j];
266             else
267                 xyz->mat[i][j]=z->mat[j][0];
268         }
269     }
270     for(i=0;i<R->row;i++)
271         delete [] R->mat[i];
272     delete [] R->mat;
273     for(i=0;i<z1->row;i++)
274         delete [] z1->mat[i];
275     delete [] z1->mat;
276     for(i=0;i<invm1->row;i++)
277         delete [] invm1->mat[i];
278     delete [] invm1->mat;
279     for(i=0;i<F->row;i++)
280         delete [] F->mat[i];
281     delete [] F->mat;
282     for(i=0;i<z->row;i++)
283         delete [] z->mat[i];
284     delete [] z->mat;
285     for(i=0;i<r->row;i++)
286         delete [] r->mat[i];
287     delete [] r->mat;
288     return xyz->mat;
289 }

下面是最終效果：

操作：使用滑鼠左鍵進行虛擬球控制模型，鍵盤上下鍵進行拉近拉遠檢視，滑鼠右鍵移動檢視中心點，點選滑鼠滾輪彈出選單，選單中可以設定地圖繪製方式，網格狀和填充。

轉：https://www.cnblogs.com/zhouchanwen/archive/2012/03/13/2394078.html

 目標：