开通博客赚积分
发布资源赚积分
分类

源码开发语言/平台
当前位置：首页 > 源码/资料 > 通讯/手机编程 > 传真(Fax)编程 > 查看源码
Earth3D.cpp
资源名称：GTIFFSample.zip　[点击查看]
上传用户：lbr_007
上传日期：2019-05-31
资源大小：282k
文件大小：28k
源码类别：

传真(Fax)编程

开发平台：

Visual C++

Earth3D.cpp：源码内容
  1.  // ogl3d.cpp
  2. //
  4. #include "stdafx.h"
  5. #include "earth3d.h"
  7. #include <math.h>
  8. #include <vector>
  9. #include <arraycontainer.h>
  10. #include <geo.h>
  12. namespace E3D {
  13. const INT32 BaseMapList = 111;
  14. const INT32 OverlayList = 112;
  15. const INT32 HeaderList = 113;
  16. const INT32 PickList = 114;
  17. const INT32 AircraftList = 115;
  18. const UINT32 CellSize = 8; // 8x8 rectangles
  19. const UINT32 PickCellSize = 2;
  20. const double hdrSize = 0.2;
  21. const double AircraftLength = 0.000050;
  22. }
  24. class V3D {
  25. public:
  26. double x;
  27. double y;
  28. double z;
  29. V3D(void){ x = y = z = 0.0;}
  30. V3D(const V3D& v){ x = y = z = 0.0; Copy(v);}
  31. V3D(double * p){ x = p[0]; y = p[1]; z = p[2];}
  32. ~V3D(void){}
  33. void Create(double x1, double x2, double x3){ x = x1; y = x2; z = x3;}
  34. void Create(double * p){ x = p[0]; y = p[1]; z = p[2];}
  35. void Copy(const V3D& v)
  36. {
  37. x = v.x;
  38. y = v.y;
  39. z = v.z;
  40. }
  41. V3D& operator = (const V3D& v)
  42. {
  43. Copy(v);
  44. return *this;
  45. }
  46. void Init(void){ x = y = z = 0.0;}
  47. };
  49. //
  50. // arbitrary rotation of the point p theta radians around the axis r to
  51. // produce the point q
  52. //
  53. void ArbitraryRotate(V3D& p, double theta, V3D& r, V3D& q)
  54. {
  55.    double costheta,sintheta;
  57.    costheta = cos(theta);
  58.    sintheta = sin(theta);
  60.    q.x = (costheta + (1 - costheta) * r.x * r.x) * p.x;
  61.    q.x += ((1 - costheta) * r.x * r.y - r.z * sintheta) * p.y;
  62.    q.x += ((1 - costheta) * r.x * r.z + r.y * sintheta) * p.z;
  64.    q.y = ((1 - costheta) * r.x * r.y + r.z * sintheta) * p.x;
  65.    q.y += (costheta + (1 - costheta) * r.y * r.y) * p.y;
  66.    q.y += ((1 - costheta) * r.y * r.z - r.x * sintheta) * p.z;
  68.    q.z = ((1 - costheta) * r.x * r.z - r.y * sintheta) * p.x;
  69.    q.z += ((1 - costheta) * r.y * r.z + r.x * sintheta) * p.y;
  70.    q.z += (costheta + (1 - costheta) * r.z * r.z) * p.z;
  71. }
  73. void GEOTex::CreateTexture(DIBSection& dib, UTMProjection * pj)
  74. {
  75. UINT32 width = dib.Width();
  76. UINT32 height = dib.Height();
  78. UINT32 w, h;
  80. if (width >= 1024) w = 1024;
  81. else if (width >= 512) w = 512;
  82. else if (width >= 256) w = 256;
  83. else if (width >= 128) w = 128;
  84. else if (width >= 64) w = 64;
  85. else if (width >= 32) w = 32;
  86. else w = 16;
  88. if (height >= 1024) h = 1024;
  89. else if (height >= 512) h = 512;
  90. else if (height >= 256) h = 256;
  91. else if (height >= 128) h = 128;
  92. else if (height >= 64) h = 64;
  93. else if (height >= 32) h = 32;
  94. else h = 16;
  96. DIBSection tmpDIB;
  97. dib.ResizeImage(tmpDIB, w, h);
  98. UINT32 totalSize = w * h * 3; // 24-bit
  99. m_imgData.CreateArray(totalSize);
  100. memcpy(m_imgData.GetData(), tmpDIB.GetBits(), totalSize);
  102. m_hScale = (double)width / (double)w;
  103. m_vScale = (double)height / (double)h;
  105. m_gr = dib.GetConstGeoReference();
  106. m_pj = pj;
  107. m_texWidth = w;
  108. m_texHeight = h;
  109. }
  111. void GEOTex::DrawTexture(void)
  112. {
  113. }
  115. Earth3D::Earth3D(void)
  116. {
  117. InitEarth3D();
  118. }
  120. Earth3D::~Earth3D(void)
  121. {
  122. }
  124. void Earth3D::GeoToXYZ(double lat, double lon, double& x, double& y, double& z, double R)
  125. {
  126. //
  127. // its necessary to use the negative of the longitude because the OpenGL
  128. // (x,y,z) coordinate system is essentially backwards from the lat/long coordinate system
  129. // where we assume negative longitudes mean west and positive means eastern hemisphere
  130. //
  131. x = R * cos((-lon*GEO::DE2RA)) * cos((lat*GEO::DE2RA));
  132. y = R * sin((lat*GEO::DE2RA));
  133. z = R * sin((-lon*GEO::DE2RA)) * cos((lat*GEO::DE2RA));
  134. }
  136. void Earth3D::SetGeoImage(DIBSection * dib, UTMProjection * pj)
  137. {
  138. m_gt.CreateTexture(*dib, pj);
  140. double x1, y1;
  141. double EL = m_model_matrix.XRot();
  142. double elev = m_eyeDist * sin(EL * GEO::DE2RA);
  144. double southLat, southLon, eastLat, eastLon, M;
  145. double southPt[3], eastPt[3];
  147. x1 = m_gt.m_gr.m_tieX + (m_gt.m_gr.m_resX * (m_gt.m_gr.m_cols/2.0));
  148. y1 = m_gt.m_gr.m_tieY - (m_gt.m_gr.m_resY * (m_gt.m_gr.m_rows/2.0));
  150. m_gt.m_pj->inverse(x1,y1,m_centerLat,m_centerLon);
  152. m_eyeDist = 1.5 * sqrt((y1 - m_gt.m_gr.m_tieY)*(y1 - m_gt.m_gr.m_tieY)) / 1000.0;
  154. GeoToXYZ(m_centerLat, m_centerLon, m_centerPt[0], m_centerPt[1], m_centerPt[2]);
  156. // setting up the rectangular block at the northern extent of the
  157. // image that says "North"
  158. m_northHdr.Create(64, 32, 24);
  160. CRect r;
  161. m_northHdr.GetDC()->DrawText("North", -1, &r, DT_CALCRECT);
  162. m_northHdr.GetDC()->PatBlt(0,0,64,32,WHITENESS);
  163. m_northHdr.GetDC()->TextOut(32-r.Width()/2, 16, "North", 5);
  165. // calculate south and eastern coordinates
  166. GeoCalc::GCDistanceAzimuth(m_centerLat, m_centerLon, m_eyeDist, 180, southLat, southLon);
  167. GeoCalc::GCDistanceAzimuth(m_centerLat, m_centerLon, m_eyeDist, 90, eastLat, eastLon);
  169. // calculate southern and eastern points in the model coordinate system
  170. GeoToXYZ(southLat,southLon, southPt[0],southPt[1],southPt[2]);
  171. GeoToXYZ(eastLat,eastLon, eastPt[0],eastPt[1],eastPt[2]);
  173. // make unitSouthPt a unit vector in the southern direction (model coord sys)
  174. m_unitSouthPt[0] = southPt[0] - m_centerPt[0];
  175. m_unitSouthPt[1] = southPt[1] - m_centerPt[1];
  176. m_unitSouthPt[2] = southPt[2] - m_centerPt[2];
  178. M = sqrt(m_unitSouthPt[0]*m_unitSouthPt[0] + m_unitSouthPt[1]*m_unitSouthPt[1] +
  179. m_unitSouthPt[2]*m_unitSouthPt[2]);
  180. m_unitSouthPt[0] /= M;
  181. m_unitSouthPt[1] /= M;
  182. m_unitSouthPt[2] /= M;
  184. // make unitEastPt a unit vector in the eastern direction (model coord sys)
  185. m_unitEastPt[0] = eastPt[0] - m_centerPt[0];
  186. m_unitEastPt[1] = eastPt[1] - m_centerPt[1];
  187. m_unitEastPt[2] = eastPt[2] - m_centerPt[2];
  189. M = sqrt(m_unitEastPt[0]*m_unitEastPt[0] + m_unitEastPt[1]*m_unitEastPt[1] +
  190. m_unitEastPt[2]*m_unitEastPt[2]);
  191. m_unitEastPt[0] /= M;
  192. m_unitEastPt[1] /= M;
  193. m_unitEastPt[2] /= M;
  195. // make unitCenterPt a unit vector pointing from the origin to the center point
  196. // (model coord sys)
  197. m_unitCenterPt[0] = m_centerPt[0];
  198. m_unitCenterPt[1] = m_centerPt[1];
  199. m_unitCenterPt[2] = m_centerPt[2];
  201. M = sqrt(m_centerPt[0]*m_centerPt[0] + m_centerPt[1]*m_centerPt[1] +
  202. m_centerPt[2]*m_centerPt[2]);
  203. m_unitCenterPt[0] /= M;
  204. m_unitCenterPt[1] /= M;
  205. m_unitCenterPt[2] /= M;
  206. }
  208. void Earth3D::SetOrientation(void)
  209. {
  210. if (m_flyMode && m_flyFirstPerson)
  211. {
  212. //
  213. // if we're flying then the orientation is simple: our eye point is the
  214. // current fly point and we're looking at the next fly point, our up
  215. // vector is the next fly point
  216. //
  217. gluLookAt(m_flyPt[0], m_flyPt[1], m_flyPt[2],
  218. m_flyTo[0], m_flyTo[1], m_flyTo[2],
  219. m_flyTo[0], m_flyTo[1], m_flyTo[2]);
  220. }
  221. else
  222. {
  223. double upVec[3];
  225. // dividing by the Earth radius essentially scales from a world
  226. // coordinate system to this OpenGL coordinate system where the
  227. // radius of the Earth equals 1.0
  228. double D = (m_eyeDist / m_model_matrix.XScale())/GEO::ERAD;
  230. V3D p(m_unitSouthPt);
  231. V3D u(m_unitCenterPt);
  232. V3D q, v, w, x;
  234. // rotate the south vector around the center vector
  235. ArbitraryRotate(p,(-m_model_matrix.YRot())*GEO::DE2RA,u,q);
  237. v.Create(m_unitEastPt);
  239. // rotate the east vector around the center vector
  240. ArbitraryRotate(v,(-m_model_matrix.YRot())*GEO::DE2RA,u,w);
  242. p = q;
  244. // rotate the (formerly) south vector around the (formerly) east vector
  245. ArbitraryRotate(p,-m_model_matrix.XRot()*GEO::DE2RA,w,q);
  247. // calculate the eye point
  248. m_eyePt[0] = m_centerPt[0] + D * q.x;
  249. m_eyePt[1] = m_centerPt[1] + D * q.y;
  250. m_eyePt[2] = m_centerPt[2] + D * q.z;
  252. p.Create(m_unitCenterPt);
  254. // rotate the center vector around the (formerly) east vector
  255. ArbitraryRotate(p,-m_model_matrix.XRot()*GEO::DE2RA,w,q);
  257. upVec[0] = q.x;
  258. upVec[1] = q.y;
  259. upVec[2] = q.z;
  261. gluLookAt( m_eyePt[0], m_eyePt[1], m_eyePt[2],
  262. m_centerPt[0], m_centerPt[1], m_centerPt[2],
  263. upVec[0], upVec[1], upVec[2]);
  264. }
  265. }
  267. void RenderCylinder(double * startPt, double * endPt, double r = 0.0000075, UINT32 steps = 100)
  268. {
  269. if (steps > 1)
  270. {
  271. int i, j;
  272. double d[3], n[3], c1[3], c2[3], u[3], v[3], w[3];
  273. double dIncr[3], angle1, D, N;
  275. d[0] = endPt[0] - startPt[0];
  276. d[1] = endPt[1] - startPt[1];
  277. d[2] = endPt[2] - startPt[2];
  279. D = sqrt(d[0]*d[0] + d[1]*d[1] + d[2]*d[2]);
  281. dIncr[0] = (endPt[0] - startPt[0]) / (double)(steps - 1);
  282. dIncr[1] = (endPt[1] - startPt[1]) / (double)(steps - 1);
  283. dIncr[2] = (endPt[2] - startPt[2]) / (double)(steps - 1);
  285. n[0] = (endPt[1]*startPt[2] - endPt[2]*startPt[1]);
  286. n[1] = (endPt[2]*startPt[0] - endPt[0]*startPt[2]);
  287. n[2] = (endPt[0]*startPt[1] - endPt[1]*startPt[0]);
  289. // making the n-vector a unit vector
  290. N = sqrt(n[0]*n[0] + n[1]*n[1] + n[2]*n[2]);
  292. n[0] /= N;
  293. n[1] /= N;
  294. n[2] /= N;
  296. n[0] *= r;
  297. n[1] *= r;
  298. n[2] *= r;
  300. const int numAngles = 36;
  301. double perp[numAngles][3];
  303. V3D p, t, q;
  304. p.x = n[0]; p.y = n[1]; p.z = n[2];
  305. t.x = d[0]/D; t.y = d[1]/D; t.z = d[2]/D;
  307. double dTheta = 360.0 / (numAngles - 1);
  309. for (i = 0; i < numAngles; i++){
  310. angle1 = (double)i*dTheta;
  311. ArbitraryRotate(p,angle1*GEO::DE2RA,t,q);
  312. perp[i][0] = q.x;
  313. perp[i][1] = q.y;
  314. perp[i][2] = q.z;
  315. }
  318. glBegin( GL_TRIANGLES );
  320. for (j = 0; j < (steps - 1); j++){
  322. c1[0] = startPt[0] + (double)j * dIncr[0];
  323. c1[1] = startPt[1] + (double)j * dIncr[1];
  324. c1[2] = startPt[2] + (double)j * dIncr[2];
  326. c2[0] = startPt[0] + (double)(j+1) * dIncr[0];
  327. c2[1] = startPt[1] + (double)(j+1) * dIncr[1];
  328. c2[2] = startPt[2] + (double)(j+1) * dIncr[2];
  330. for (i = 0; i < (numAngles-1); i++){
  332. u[0] = c1[0] + perp[i][0];
  333. u[1] = c1[1] + perp[i][1];
  334. u[2] = c1[2] + perp[i][2];
  336. v[0] = c2[0] + perp[i][0];
  337. v[1] = c2[1] + perp[i][1];
  338. v[2] = c2[2] + perp[i][2];
  340. w[0] = c1[0] + perp[i+1][0];
  341. w[1] = c1[1] + perp[i+1][1];
  342. w[2] = c1[2] + perp[i+1][2];
  344. glNormal3dv(perp[i]);
  345. glVertex3dv(u);
  346. glVertex3dv(v);
  347. glVertex3dv(w);
  349. u[0] = c2[0] + perp[i+1][0];
  350. u[1] = c2[1] + perp[i+1][1];
  351. u[2] = c2[2] + perp[i+1][2];
  353. glNormal3dv(perp[i+1]);
  354. glVertex3dv(w);
  355. glVertex3dv(v);
  356. glVertex3dv(u);
  357. }
  358. }
  360. glEnd();
  361. }
  362. }
  364. void Earth3D::SetTexture(DIBSection& dib)
  365. {
  366. if (dib.IsCreated())
  367. {
  368. GenTexture(dib);
  369. }
  370. }
  372. void Earth3D::GenTexture(DIBSection& dib)
  373. {
  374. UINT32 width = dib.Width();
  375. UINT32 height = dib.Height();
  377. int w, h;
  379. UINT32 max_dimension = width > height ? width : height;
  381. if (width >= 512) w = 512;
  382. else if (width >= 256) w = 256;
  383. else if (width >= 128) w = 128;
  384. else if (width >= 64) w = 64;
  385. else if (width >= 32) w = 32;
  386. else w = 16;
  388. if (height >= 512) h = 512;
  389. else if (height >= 256) h = 256;
  390. else if (height >= 128) h = 128;
  391. else if (height >= 64) h = 64;
  392. else if (height >= 32) h = 32;
  393. else h = 16;
  395. dib.ResizeImage(m_texture_dib,w,h);
  396. }
  398. void Earth3D::RenderPrimaryImage(void)
  399. {
  400. glPolygonMode(GL_FRONT_AND_BACK,GL_FILL);
  401. if (m_pickMode)
  402. {
  403. glCallList( E3D::PickList );
  405. m_pickPoint = false;
  407. if (m_dib.IsCreated())
  408. {
  409. COLORREF clr;
  410. m_dib.GetPixel(m_pickX, m_pickY, clr);
  412. INT32 horzCells = (m_gt.m_texWidth - E3D::PickCellSize)/E3D::PickCellSize;
  413. INT32 vertCells = (m_gt.m_texHeight - E3D::PickCellSize)/E3D::PickCellSize;
  415. INT32 row = clr / horzCells;
  416. INT32 col = clr - (row * horzCells);
  418. if ((row >= 0) && (row < vertCells) && (col >= 0) && (col < horzCells))
  419. {
  420. double x = m_gt.m_gr.m_tieX + (double)col * m_gt.m_gr.m_resX * m_gt.m_hScale*E3D::PickCellSize;
  421. x += m_gt.m_gr.m_tieX + (double)(col+E3D::PickCellSize) * m_gt.m_gr.m_resX * m_gt.m_hScale*E3D::PickCellSize;
  422. x /= 2.0;
  424. double y = m_gt.m_gr.m_tieY - (double)row * m_gt.m_gr.m_resY * m_gt.m_vScale*E3D::PickCellSize;
  425. y += m_gt.m_gr.m_tieY - (double)(row+E3D::PickCellSize) * m_gt.m_gr.m_resY * m_gt.m_vScale*E3D::PickCellSize;
  426. y /= 2.0;
  428. m_gt.m_pj->inverse(x,y,m_pickLat,m_pickLon);
  429. m_pickPoint = true;
  430. }
  431. }
  432. }
  433. else
  434. {
  435. if (m_flyMode)
  436. {
  437. RenderDirectionHeaders(m_northHdr, 0);
  438. RenderEarthImage(); // shows the earth underneath image
  439. RenderGEOImage();
  440. if (!m_flyFirstPerson)
  441. {
  442. RenderAircraftImage();
  443. }
  444. }
  445. else
  446. {
  447. RenderPickImage();
  448. RenderDirectionHeaders(m_northHdr, 0);
  449. RenderEarthImage(); // shows the earth underneath image
  450. RenderGEOImage();
  451. }
  452. }
  453. }
  455. void Earth3D::SetFlyPosition(double u[3], double v[3], bool firstPerson)
  456. {
  457. m_flyPt[0] = u[0];
  458. m_flyPt[1] = u[1];
  459. m_flyPt[2] = u[2];
  460. m_flyTo[0] = v[0];
  461. m_flyTo[1] = v[1];
  462. m_flyTo[2] = v[2];
  463. m_flyFirstPerson = firstPerson;
  464. }
  466. void Earth3D::RenderGEOImage(void)
  467. {
  468. if (m_gt.IsCreated() && m_gt.m_pj)
  469. {
  470. if (m_rebuildTexture)
  471. {
  472. glNewList( E3D::OverlayList, GL_COMPILE);
  474. UINT32 row, col;
  475. double x1, x2, y1, y2;
  476. double lat[4], lon[4];
  477. GLdouble u[3], v[3], w[3], n[3], tx1, tx2, ty1, ty2;
  478. GLdouble R = 1.0000;
  480. glColor3ub(32,32,32);
  482. glTexImage2D(GL_TEXTURE_2D,0,3, m_gt.m_texWidth, m_gt.m_texHeight,
  483. 0,GL_BGR_EXT,GL_UNSIGNED_BYTE,(GLvoid *)m_gt.m_imgData.GetData());
  485. glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_REPEAT);
  486. glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_REPEAT);
  487. glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST);
  488. glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST);
  489. glTexEnvf(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_DECAL);
  491. glEnable( GL_TEXTURE_2D );
  492. {
  493. glBegin( GL_TRIANGLES );
  494. {
  495. //
  496. // Notes:
  497. // 1. each texture mapping point is like a spot weld so
  498. // most of the time we don't need to texture map every point,
  499. // instead we will break up the image into 2D cells.
  500. // 2. all four points of the cell must be un-projected for
  501. // accuracy...inverting a projection is slow...another reason
  502. // we don't want to texture map every point.
  503. // 3. each cell will be broken into 2 triangles and we have to be
  504. // careful that the normal vector for the 2 triangles are both
  505. // pointing the same direction and that direction must be up.
  506. // 4. one optimization i use below is to re-use the v & w points
  507. // from the first triangle in the second triangle because those
  508. // points are used in both. you may be able to come up with much
  509. // more re-use of points and improve performance of this code, but
  510. // keep in mind we are only building an OpenGL call list here so
  511. // this code does not get executed on each and every re-draw, the
  512. // non-opengl code only gets executed when this model is re-built.
  513. // this fact pushes us more towards using high resolution up to the
  514. // point that we are creating too much opengl code.
  515. // 5. set the texture coordinate before setting the vertex. the
  516. // texture coordinate is handled like a state when the vertex coordinates
  517. // get set.
  518. // 6. upper left texture coordinate is (1.0,1.0) as you go right and down
  519. // the texture coordinate moves towards (0.0, 0.0)
  520. //
  521. for (row = 0; row < (m_gt.m_texHeight - E3D::CellSize); row += E3D::CellSize){
  523. //
  524. // x1, y1, x2, y2 are based off the coordinate system of the 
  525. // original geotiff image
  526. //
  527. y1 = m_gt.m_gr.m_tieY - m_gt.m_gr.m_resY * ((double)row * m_gt.m_vScale);
  528. y2 = m_gt.m_gr.m_tieY - m_gt.m_gr.m_resY * ((double)(row+E3D::CellSize) * m_gt.m_vScale);
  530. //
  531. // tx1, ty1, tx2, ty2 are based off the image data that was re-sized to
  532. // meet OpenGL requirements for image dimensions that are powers of two.
  533. //
  534. ty1 = 1.0 - (double)row/(double)m_gt.m_texHeight;
  535. ty2 = 1.0 - (double)(row + E3D::CellSize)/(double)m_gt.m_texHeight;
  537. for (col = 0; col < (m_gt.m_texWidth - E3D::CellSize); col += E3D::CellSize){
  539. x1 = m_gt.m_gr.m_tieX + m_gt.m_gr.m_resX * ((double)col * m_gt.m_hScale);
  540. x2 = m_gt.m_gr.m_tieX + m_gt.m_gr.m_resX * ((double)(col + E3D::CellSize) * m_gt.m_hScale);
  542. tx1 = (double)col/(double)(m_gt.m_texWidth);
  543. tx2 = (double)(col + E3D::CellSize)/(double)(m_gt.m_texWidth);
  545. //
  546. // going clockwise around this rectangle for inverse projections
  547. //
  548. m_gt.m_pj->inverse( x1, y1, lat[0], lon[0]);
  549. m_gt.m_pj->inverse( x2, y1, lat[1], lon[1]);
  550. m_gt.m_pj->inverse( x1, y2, lat[2], lon[2]);
  551. m_gt.m_pj->inverse( x2, y2, lat[3], lon[3]);
  553. GeoToXYZ(lat[0],lon[0],u[0],u[1],u[2]);
  554. GeoToXYZ(lat[1],lon[1],v[0],v[1],v[2]);
  555. GeoToXYZ(lat[2],lon[2],w[0],w[1],w[2]);
  557. //
  558. // i'm switching the order of the coordinates u,w,v instead of
  559. // u,v,w so that we get counter-clockwise direction. using the right
  560. // hand rule this means that the normal vector will point up. here
  561. // we do upper-left, lower-left, upper-right
  562. //
  563. NormalVector(u,w,v,n);
  565. glNormal3dv(n);
  566. glTexCoord2d(tx1, ty1);
  567. glVertex3dv(u);
  568. glTexCoord2d(tx1, ty2);
  569. glVertex3dv(w);
  570. glTexCoord2d(tx2, ty1);
  571. glVertex3dv(v);
  573. GeoToXYZ(lat[3],lon[3],u[0],u[1],u[2]);
  575. //
  576. // again, we want to go counter clockwise for the second triangle
  577. // (the lower right triangle) to keep our normal vector going up.
  578. // so we want the latitude & longitude of points lower-left then
  579. // lower-right, then upper-right...and this completes our cell.
  580. //
  581. NormalVector(w,u,v,n);
  582. glNormal3dv(n);
  583. glTexCoord2d(tx1, ty2);
  584. glVertex3dv(w);
  585. glTexCoord2d(tx2, ty2);
  586. glVertex3dv(u);
  587. glTexCoord2d(tx2, ty1);
  588. glVertex3dv(v);
  589. }
  590. }
  591. }
  592. glEnd();
  593. }
  594. glDisable( GL_TEXTURE_2D );
  596. glEndList();
  598. m_rebuildTexture = false;
  599. }
  601. glCallList( E3D::OverlayList );
  602. }
  603. }
  605. void Earth3D::RenderAircraftImage(void)
  606. {
  607. glPushMatrix();
  608. {
  609. //
  610. // we're going to draw a simple aircraft made up of
  611. // 3 triangles: a fuselage, wings, and rudder. First
  612. // we have to get unit vectors for the aircraft's orientation
  613. // and direction.
  614. //
  615. GLdouble u[3], v[3], w[3], n1[3], n2[3], n3[3], norm[3], m;
  616. u[0] = m_flyTo[0] - m_flyPt[0];
  617. u[1] = m_flyTo[1] - m_flyPt[1];
  618. u[2] = m_flyTo[2] - m_flyPt[2];
  620. m = sqrt(u[0]*u[0] + u[1]*u[1] + u[2]*u[2]);
  622. // u is the unit vector for the direction of the aircraft
  623. u[0] /= m;
  624. u[1] /= m;
  625. u[2] /= m;
  627. v[0] = m_flyPt[0];
  628. v[1] = m_flyPt[1];
  629. v[2] = m_flyPt[2];
  631. m = sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]);
  633. // v is the unit vector pointing up
  634. v[0] /= m;
  635. v[1] /= m;
  636. v[2] /= m;
  638. // cross product of the up vector and the forward vector should
  639. // the vector to the left wing.
  640. w[0] = v[1]*u[2] - v[2]*u[1];
  641. w[1] = v[2]*u[0] - v[0]*u[2];
  642. w[2] = v[0]*u[1] - v[1]*u[0];
  644. m = sqrt(w[0]*w[0] + w[1]*w[1] + w[2]*w[2]);
  646. // w is the unit vector from the flight point to the left wing
  647. w[0] /= m;
  648. w[1] /= m;
  649. w[2] /= m;
  651. glTranslated( m_flyTo[0], m_flyTo[1], m_flyTo[2]);
  653. //if (m_buildAircraft)
  654. //{
  655. //glNewList( E3D::AircraftList, GL_COMPILE);
  656. //m_buildAircraft = false;
  658. glBegin( GL_TRIANGLES );
  659. {
  660. glColor3ub(64,64,224); // light blue...i think???
  661. // first six vertices will be the aircraft fuselage
  662. n1[0] = E3D::AircraftLength*u[0];
  663. n1[1] = E3D::AircraftLength*u[1];
  664. n1[2] = E3D::AircraftLength*u[2];
  666. n2[0] = E3D::AircraftLength*w[0]*0.20;
  667. n2[1] = E3D::AircraftLength*w[1]*0.20;
  668. n2[2] = E3D::AircraftLength*w[2]*0.20;
  670. n3[0] = E3D::AircraftLength*u[0]*0.20;
  671. n3[1] = E3D::AircraftLength*u[1]*0.20;
  672. n3[2] = E3D::AircraftLength*u[2]*0.20;
  674. NormalVector(n1,n2,n3,norm);
  675. glNormal3dv(norm);
  676. glVertex3dv(n1);
  677. glVertex3dv(n2);
  678. glVertex3dv(n3);
  680. n1[0] = E3D::AircraftLength*u[0];
  681. n1[1] = E3D::AircraftLength*u[1];
  682. n1[2] = E3D::AircraftLength*u[2];
  684. n2[0] = E3D::AircraftLength*(-w[0])*0.20;
  685. n2[1] = E3D::AircraftLength*(-w[1])*0.20;
  686. n2[2] = E3D::AircraftLength*(-w[2])*0.20;
  688. NormalVector(n1,n3,n2,norm);
  689. glNormal3dv(norm);
  690. glVertex3dv(n1);
  691. glVertex3dv(n3);
  692. glVertex3dv(n2);
  694. // next lets do the wings
  695. glColor3ub(64,64,128);
  696. n1[0] = E3D::AircraftLength*u[0]*0.75;
  697. n1[1] = E3D::AircraftLength*u[1]*0.75;
  698. n1[2] = E3D::AircraftLength*u[2]*0.75;
  700. n2[0] = -E3D::AircraftLength*w[0]*0.50 + E3D::AircraftLength*u[0]*0.20;
  701. n2[1] = -E3D::AircraftLength*w[1]*0.50 + E3D::AircraftLength*u[1]*0.20;
  702. n2[2] = -E3D::AircraftLength*w[2]*0.50 + E3D::AircraftLength*u[2]*0.20;
  704. n3[0] = E3D::AircraftLength*(w[0])*0.50 + E3D::AircraftLength*u[0]*0.20;
  705. n3[1] = E3D::AircraftLength*(w[1])*0.50 + E3D::AircraftLength*u[1]*0.20;
  706. n3[2] = E3D::AircraftLength*(w[2])*0.50 + E3D::AircraftLength*u[2]*0.20;
  708. NormalVector(n1,n2,n3,norm);
  709. glNormal3dv(norm);
  710. glVertex3dv(n1);
  711. glVertex3dv(n2);
  712. glVertex3dv(n3);
  714. // next lets do the rudders
  715. glColor3ub(192,192,192);
  716. n1[0] = E3D::AircraftLength*u[0]*0.50;
  717. n1[1] = E3D::AircraftLength*u[1]*0.50;
  718. n1[2] = E3D::AircraftLength*u[2]*0.50;
  720. n2[0] = 0.0;
  721. n2[1] = 0.0;
  722. n2[2] = 0.0;
  724. n3[0] = E3D::AircraftLength*v[0]*0.15;
  725. n3[1] = E3D::AircraftLength*v[1]*0.15;
  726. n3[2] = E3D::AircraftLength*v[2]*0.15;
  728. NormalVector(n1,n2,n3,norm);
  729. glNormal3dv(norm);
  730. glVertex3dv(n1);
  731. glVertex3dv(n2);
  732. glVertex3dv(n3);
  733. }
  734. glEnd();
  736. //glEndList();
  737. //}
  738. //else
  739. //{
  740. // glCallList( E3D::AircraftList );
  741. //}
  742. }
  743. glPopMatrix();
  744. }
  746. void Earth3D::RenderPickImage(void)
  747. {
  748. if (m_gt.IsCreated() && m_gt.m_pj)
  749. {
  750. if (m_rebuildTexture)
  751. {
  752. glNewList( E3D::PickList, GL_COMPILE);
  754. glShadeModel(GL_FLAT);
  755. glDisable(GL_LIGHT0);
  756. glDisable(GL_LIGHTING);
  757. glDisable( GL_TEXTURE_2D );
  759. UINT32 row, col;
  760. double x1, x2, y1, y2;
  761. double lat[4], lon[4];
  762. GLdouble u[3], v[3], w[3];
  764. COLORREF clr = 0;
  766. glBegin( GL_TRIANGLES );
  767. {
  768. for (row = 0; row < (m_gt.m_texHeight - E3D::PickCellSize); row+=E3D::PickCellSize){
  770. //
  771. // x1, y1, x2, y2 are based off the coordinate system of the 
  772. // original geotiff image
  773. //
  774. y1 = m_gt.m_gr.m_tieY - m_gt.m_gr.m_resY * ((double)row * m_gt.m_vScale);
  775. y2 = m_gt.m_gr.m_tieY - m_gt.m_gr.m_resY * ((double)(row+E3D::PickCellSize) * m_gt.m_vScale);
  777. for (col = 0; col < (m_gt.m_texWidth - E3D::PickCellSize); col+=E3D::PickCellSize){
  779. x1 = m_gt.m_gr.m_tieX + m_gt.m_gr.m_resX * ((double)col * m_gt.m_hScale);
  780. x2 = m_gt.m_gr.m_tieX + m_gt.m_gr.m_resX * ((double)(col + E3D::PickCellSize) * m_gt.m_hScale);
  782. //
  783. // going clockwise around this rectangle for inverse projections
  784. //
  785. m_gt.m_pj->inverse( x1, y1, lat[0], lon[0]);
  786. m_gt.m_pj->inverse( x2, y1, lat[1], lon[1]);
  787. m_gt.m_pj->inverse( x1, y2, lat[2], lon[2]);
  788. m_gt.m_pj->inverse( x2, y2, lat[3], lon[3]);
  790. GeoToXYZ(lat[0],lon[0],u[0],u[1],u[2]);
  791. GeoToXYZ(lat[1],lon[1],v[0],v[1],v[2]);
  792. GeoToXYZ(lat[2],lon[2],w[0],w[1],w[2]);
  794. glColor3ub(GetRValue(clr),GetGValue(clr),GetBValue(clr));
  795. clr++;
  797. glVertex3dv(u);
  798. glVertex3dv(w);
  799. glVertex3dv(v);
  801. GeoToXYZ(lat[3],lon[3],u[0],u[1],u[2]);
  803. glVertex3dv(w);
  804. glVertex3dv(u);
  805. glVertex3dv(v);
  806. }
  807. }
  808. }
  809. glEnd();
  811. glEndList();
  812. }
  813. }
  814. }
  816. void Earth3D::RenderDirectionHeaders(DIBSection& dib, int direction)
  817. {
  818. if (m_gt.IsCreated() && m_gt.m_pj)
  819. {
  820. if (m_rebuildTexture)
  821. {
  822. glNewList( E3D::HeaderList, GL_COMPILE);
  824. double x1, x2, y1, y2;
  825. double lat[4], lon[4];
  826. GLdouble u[3], v[3], w[3], n[3], e1[3], e2[3];
  827. GLdouble R = 1.0;
  829. glTexImage2D(GL_TEXTURE_2D,0,3, dib.Width(), dib.Height(),
  830. 0,GL_BGR_EXT,GL_UNSIGNED_BYTE,(GLvoid *)dib.GetBits());
  832. glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_REPEAT);
  833. glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_REPEAT);
  834. glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST);
  835. glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST);
  836. glTexEnvf(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_DECAL);
  838. glEnable( GL_TEXTURE_2D );
  839. {
  840. glBegin( GL_TRIANGLES );
  841. {
  842. x1 = m_gt.m_gr.m_tieX;
  843. x2 = x1 + m_gt.m_gr.m_resX * m_gt.m_gr.m_cols;
  844. y1 = m_gt.m_gr.m_tieY;
  845. y2 = y1 - m_gt.m_gr.m_resY * m_gt.m_gr.m_rows;
  847. m_gt.m_pj->inverse( x1, y1, lat[0], lon[0]);
  848. m_gt.m_pj->inverse( x2, y1, lat[1], lon[1]);
  849. m_gt.m_pj->inverse( x1, y2, lat[2], lon[2]);
  850. m_gt.m_pj->inverse( x2, y2, lat[3], lon[3]);
  852. NormalVector(u,w,v,n);
  853. glNormal3dv(n);
  855. switch (direction){
  856. case 0: // north
  857. {
  858. GeoToXYZ(lat[0], lon[0],u[0],u[1],u[2]);
  859. GeoToXYZ(lat[1], lon[1],v[0],v[1],v[2]);
  860. GeoToXYZ(lat[2], lon[2],w[0],w[1],w[2]);
  862. glTexCoord2d(0.0, 0.0);
  863. glVertex3dv(u);
  864. glTexCoord2d(1.0, 0.0);
  865. glVertex3dv(v);
  867. e1[0] = u[0] + E3D::hdrSize * (u[0] - w[0]);
  868. e1[1] = u[1] + E3D::hdrSize * (u[1] - w[1]);
  869. e1[2] = u[2] + E3D::hdrSize * (u[2] - w[2]);
  871. GeoToXYZ(lat[3], lon[3],w[0],w[1],w[2]);
  873. e2[0] = v[0] + E3D::hdrSize * (v[0] - w[0]);
  874. e2[1] = v[1] + E3D::hdrSize * (v[1] - w[1]);
  875. e2[2] = v[2] + E3D::hdrSize * (v[2] - w[2]);
  877. glTexCoord2d(1.0, 1.0);
  878. glVertex3dv(e2);
  880. glTexCoord2d(0.0, 0.0);
  881. glVertex3dv(u);
  882. glTexCoord2d(1.0, 1.0);
  883. glVertex3dv(e2);
  884. glTexCoord2d(0.0, 1.0);
  885. glVertex3dv(e1);
  886. }
  887. break;
  888. case 1: // south
  889. {
  890. }
  891. break;
  892. case 2: // east
  893. {
  894. }
  895. break;
  896. case 3: // west
  897. {
  898. }
  899. break;
  900. }
  901. }
  902. glEnd();
  903. }
  904. glDisable( GL_TEXTURE_2D );
  906. glEndList();
  907. }
  909. glCallList( E3D::HeaderList );
  910. }
  911. }
  913. void Earth3D::RenderEarthImage(void)
  914. {
  915. if (m_rebuildEarth)
  916. {
  917. glNewList( E3D::BaseMapList, GL_COMPILE);
  919. const int NumLatitudes = 90;
  920. const int NumLongitudes = 180;
  922. double start_lat, start_lon;
  923. double theta1, phi1, theta2, phi2, lat_incr, lon_incr;
  925. GLdouble u[3], v[3], w[3], n[3];
  926. GLdouble R;
  928. glColor3ub(224,224,224);
  930. /*
  931. if (m_texture_dib.IsCreated())
  932. {
  933. glTexImage2D(GL_TEXTURE_2D,0,3,m_texture_dib.Width(),m_texture_dib.Height(),
  934. 0,GL_BGR_EXT,GL_UNSIGNED_BYTE,(GLvoid *)m_texture_dib.GetBits());
  936. glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_REPEAT);
  937. glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_REPEAT);
  938. glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST);
  939. glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST);
  940. glTexEnvf(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_DECAL);
  942. glEnable( GL_TEXTURE_2D );
  943. }
  944. else
  945. {
  946. glDisable( GL_TEXTURE_2D );
  947. }
  948. */
  950. start_lat = -90.0;
  951. start_lon = -180.0;
  952. R = 1.0;
  954. lat_incr = 180.0 / (double)NumLatitudes;
  955. lon_incr = 360.0 / (double)NumLongitudes;
  957. int row, col;
  958. //GLdouble tx1, tx2, ty1, ty2;
  960. glBegin( GL_TRIANGLES );
  961. {
  962. for (col = 0; col < NumLongitudes; col++){
  963. phi1 = (start_lon + col * lon_incr);
  964. phi2 = (start_lon + (col + 1) * lon_incr);
  966. //tx1 = (180.0 + phi1)/360.0;
  967. //tx2 = (180.0 + phi2)/360.0;
  969. for (row = 0; row < NumLatitudes; row++){
  970. theta1 = (start_lat + row * lat_incr);
  971. theta2 = (start_lat + (row + 1) * lat_incr);
  973. //ty1 = (90.0 + theta1)/180.0;
  974. //ty2 = (90.0 + theta2)/180.0;
  976. GeoToXYZ(theta1,phi1,u[0],u[1],u[2]);
  977. GeoToXYZ(theta2,phi1,v[0],v[1],v[2]);
  978. GeoToXYZ(theta2,phi2,w[0],w[1],w[2]);
  980. NormalVector(u,v,w,n);
  981. glNormal3dv(n);
  982. //glTexCoord2d(tx1,ty1);
  983. glVertex3dv(u);
  984. //glTexCoord2d(tx1,ty2);
  985. glVertex3dv(v);
  986. //glTexCoord2d(tx2,ty2);
  987. glVertex3dv(w);
  989. GeoToXYZ(theta1,phi2,v[0],v[1],v[2]);
  991. NormalVector(u,w,v,n);
  992. glNormal3dv(n);
  993. //glTexCoord2d(tx1,ty1);
  994. glVertex3dv(u);
  995. //glTexCoord2d(tx2,ty2);
  996. glVertex3dv(w);
  997. //glTexCoord2d(tx2,ty1);
  998. glVertex3dv(v);
  999. }
  1000. }
  1001. }
  1002. glEnd();
  1004. //glDisable( GL_TEXTURE_2D );
  1006. glEndList();
  1008. m_rebuildEarth = false;
  1010. }
  1012. glCallList( E3D::BaseMapList );
  1013. }