1 /***************************************************************************
2 * This file is part of Tecorrec. *
3 * Copyright 2008 James Hogan <james@albanarts.com> *
5 * Tecorrec is free software: you can redistribute it and/or modify *
6 * it under the terms of the GNU General Public License as published by *
7 * the Free Software Foundation, either version 2 of the License, or *
8 * (at your option) any later version. *
10 * Tecorrec is distributed in the hope that it will be useful, *
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
13 * GNU General Public License for more details. *
15 * You should have received a copy of the GNU General Public License *
16 * along with Tecorrec. If not, write to the Free Software Foundation, *
17 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. *
18 ***************************************************************************/
22 * @brief Manages data for a globe.
26 #include "tcGeoImageData.h"
27 #include "tcLandsatData.h"
28 #include "tcProcessingData.h"
29 #include "tcCustomSunDirection.h"
30 #include "tcSrtmModel.h"
39 * Constructors + destructor
42 /// Primary constructor.
43 tcGlobe::tcGlobe(double meanRadius
)
44 : m_meanRadius(meanRadius
)
45 , m_elevation(new tcSrtmModel
[2])
48 , m_elevationInterpolation(1.0f
)
49 , m_colourCoding(NoColourCoding
)
52 for (int i
= 0; i
< 2; ++i
)
54 m_elevationMode
[i
] = CorrectedElevation
;
56 for (int i
= 0; i
< 3; ++i
)
58 m_colourMapping
[i
][0] = 0;
59 m_colourMapping
[i
][1] = 2-i
;
61 for (int i
= 3; i
< 6; ++i
)
63 m_colourMapping
[i
][0] = -1;
64 m_colourMapping
[i
][1] = -1;
66 QList
<tcShadowClassifyingData
*> datasets
;
67 // this one is a bit rubbish, the sun is too high
68 //datasets << new tcLandsatData("/home/james/cs/pro/data/LE71950281999206EDC01/", m_elevation);
70 // high azimuth, low elevation
71 datasets
<< new tcLandsatData("/home/james/cs/pro/data/LE71950282001307EDC00/", m_elevation
);
72 datasets
<< new tcLandsatData("/home/james/cs/pro/data/LE71950282002006EDC00/", m_elevation
);
73 datasets
<< new tcLandsatData("/home/james/cs/pro/data/LE71950282002310EDC00/", m_elevation
);
75 datasets
<< new tcLandsatData("/home/james/cs/pro/data/LE71950282000081EDC00/", m_elevation
);
76 datasets
<< new tcLandsatData("/home/james/cs/pro/data/etp195r28_5t19900910/", m_elevation
);
78 #define NEW_DATASET(AZ, EL) \
79 datasets << new tcCustomSunDirection(m_elevation+1, tcGeo((180.0-AZ)*M_PI/180, (EL)*M_PI/180), \
80 tcGeo(7.0*M_PI/180, 46.0*M_PI/180))
82 for (int i = 0; i < 360; i += 5)
84 for (int j = 30; j <= 30; j += 2)
86 NEW_DATASET((double)i, (double)j);
91 foreach (tcShadowClassifyingData
* dataset
, datasets
)
95 addImagery(new tcProcessingData(datasets
, m_elevation
));
101 foreach (tcGeoImageData
* data
, m_imagery
)
105 delete [] m_elevation
;
112 /// Add some imagery data.
113 void tcGlobe::addImagery(tcGeoImageData
* data
)
115 m_imagery
.push_back(data
);
118 /// Get the imagery data.
119 const QList
<tcGeoImageData
*>& tcGlobe::imagery() const
124 /// Get the elevation model.
125 tcSrtmModel
* tcGlobe::dem() const
134 /// Draw a line of latitude.
135 void tcGlobe::drawLineOfLatitude(double latitude
) const
137 double z
= sin(latitude
) * m_meanRadius
;
138 double xy
= cos(latitude
) * m_meanRadius
;
139 glBegin(GL_LINE_LOOP
);
141 for (int lon
= 0; lon
< 360; ++lon
)
143 glVertex3d(xy
*sin(M_PI
/180*lon
), xy
*cos(M_PI
/180*lon
), z
);
150 void tcGlobe::renderCell(tcObserver
* const observer
, const tcGeo
& swCorner
, const tcGeo
& neCorner
, int samples
, bool normals
,
151 bool northEdge
, bool eastEdge
, bool southEdge
, bool westEdge
) const
153 // Sample at a sensible level
154 tcSrtmModel::RenderState state
;
155 m_elevation
->sampleAlign(swCorner
, neCorner
, &state
, samples
);
157 if (state
.moreAvailableLon
|| state
.moreAvailableLat
)
159 // Find the square distances to each corner
161 glGetv(GL_MODELVIEW_MATRIX
, &modelview
);
162 tcGeo geoCorners
[4] = {
163 tcGeo(swCorner
.lon(), swCorner
.lat()),
164 tcGeo(swCorner
.lon(), neCorner
.lat()),
165 tcGeo(neCorner
.lon(), neCorner
.lat()),
166 tcGeo(neCorner
.lon(), swCorner
.lat())
168 GLvec3d cartCorners
[4];
170 double altitudeMean
= m_meanRadius
+ altitudeAt((geoCorners
[0] + geoCorners
[1] + geoCorners
[2] + geoCorners
[3])/4);
171 for (int i
= 0; i
< 4; ++i
)
173 cartCorners
[i
] = (GLvec3d
)geoCorners
[i
] * altitudeMean
;
174 toCorners
[i
] = (modelview
*(cartCorners
[i
], 1.0)).slice
<0,3>().sqr();
175 // Cull faces which are roughly backfacing
178 if ((modelview
*(cartCorners
[i
], 0.0))[2] <= 0.0)
185 // Decide whether to subdivide
186 float diagonal
= ( (cartCorners
[0]-cartCorners
[2]).sqr()
187 + (cartCorners
[3]-cartCorners
[1]).sqr())/2*4;
188 // If it is disproportionately tall, only subdivide horizontally
189 bool tall
= (cartCorners
[1] - cartCorners
[0]).sqr() > (cartCorners
[3] - cartCorners
[0]).sqr()*4.0
190 || (cartCorners
[2] - cartCorners
[3]).sqr() > (cartCorners
[2] - cartCorners
[1]).sqr()*4.0;
191 // If it is disproportionately wide, only subdivide vertically
192 bool wide
= (cartCorners
[3] - cartCorners
[0]).sqr() > (cartCorners
[1] - cartCorners
[0]).sqr()*4.0
193 || (cartCorners
[2] - cartCorners
[1]).sqr() > (cartCorners
[2] - cartCorners
[3]).sqr()*4.0;
194 bool subdivide
= true;
195 for (int i
= 0; i
< 4; ++i
)
197 if (toCorners
[i
] > diagonal
)
209 renderCell(observer
, geoCorners
[0], (geoCorners
[3] + geoCorners
[2]) * 0.5, samples
, normals
,
210 false, eastEdge
, southEdge
, westEdge
);
212 renderCell(observer
, (geoCorners
[0] + geoCorners
[1]) * 0.5, geoCorners
[2], samples
, normals
,
213 northEdge
, eastEdge
, false, westEdge
);
215 else if (wide
&& !tall
)
218 renderCell(observer
, geoCorners
[0], (geoCorners
[1] + geoCorners
[2]) * 0.5, samples
, normals
,
219 northEdge
, false, southEdge
, westEdge
);
221 renderCell(observer
, (geoCorners
[0] + geoCorners
[3]) * 0.5, geoCorners
[2], samples
, normals
,
222 northEdge
, eastEdge
, southEdge
, false);
227 renderCell(observer
, geoCorners
[0], (geoCorners
[0] + geoCorners
[2]) * 0.5, samples
, normals
,
228 false, false, southEdge
, westEdge
);
230 renderCell(observer
, (geoCorners
[0] + geoCorners
[3]) * 0.5, (geoCorners
[3] + geoCorners
[2]) * 0.5, samples
, normals
,
231 false, eastEdge
, southEdge
, false);
233 renderCell(observer
, (geoCorners
[0] + geoCorners
[1]) * 0.5, (geoCorners
[1] + geoCorners
[2]) * 0.5, samples
, normals
,
234 northEdge
, false, false, westEdge
);
236 renderCell(observer
, (geoCorners
[0] + geoCorners
[2]) * 0.5, geoCorners
[2], samples
, normals
,
237 northEdge
, eastEdge
, false, false);
246 glColor3f(0.0f
, 0.0f
, 1.0f
);
250 glColor3f(1.0f
, 0.5f
, 0.0f
);
253 glBegin(GL_LINE_LOOP
);
254 for (int i
= 0; i
< 4; ++i
)
256 glVertex3(cartCorners
[i
]);
261 #define EDGE_KERNEL_1 \
262 glColor4f(0.5f, 0.3f, 0.2f, 1.0f); \
263 glVertex3(dir * (m_meanRadius+alt)); \
264 glColor4f(0.5f, 0.5f, 0.5f, 1.0f); \
265 glVertex3(dir * m_meanRadius);
266 #define EDGE_KERNEL_2 \
267 glColor4f(0.5f, 0.5f, 1.0f, 1.0f); \
268 glVertex3(dir * m_meanRadius); \
269 glColor4f(0.5f, 0.5f, 1.0f, 1.0f); \
270 glVertex3(dir * (m_meanRadius-(accurate ? 100.0 : 1000.0)));
271 #define EDGE_KERNEL_3 \
272 glColor4f(0.5f, 0.5f, 0.5f, 1.0f); \
273 glVertex3(dir * (m_meanRadius-(accurate ? 100.0 : 1000.0))); \
274 glColor4f(0.5f, 0.5f, 0.5f, 1.0f); \
275 glVertex3(dir * (m_meanRadius-5000.0));
276 #define EDGE_KERNEL_4 \
277 glColor4f(0.5f, 0.5f, 0.5f, 1.0f); \
278 glVertex3(dir * (m_meanRadius-5000.0)); \
279 glColor4f(1.0f, 0.0f, 0.0f, 0.5f); \
280 glVertex3(dir * (m_meanRadius-8000.0));
281 #define EDGE_KERNEL(STAGE, EDGE, SAMPLES, LON, LAT) \
284 glBegin(GL_TRIANGLE_STRIP); \
285 for (int i = 0; i < SAMPLES; ++i) \
288 bool accurate = true; \
289 double alt = altitudeAt(state, (LON), (LAT), &coord, &accurate); \
290 GLvec3d dir = coord; \
291 EDGE_KERNEL_##STAGE \
295 #define EDGE(STAGE) \
296 EDGE_KERNEL(STAGE, north, state.samplesLon, i, state.samplesLat-1); \
297 EDGE_KERNEL(STAGE, east, state.samplesLat, state.samplesLon-1, i); \
298 EDGE_KERNEL(STAGE, south, state.samplesLon, i, 0); \
299 EDGE_KERNEL(STAGE, west, state.samplesLat, 0, i);
304 float normColours
[3][2][3] = {
305 { { 1.0f
, 1.0f
, 1.0f
}, { 1.0f
, 0.0f
, 0.0f
} },
306 { { 1.0f
, 1.0f
, 1.0f
}, { 0.0f
, 1.0f
, 0.0f
} },
307 { { 1.0f
, 1.0f
, 1.0f
}, { 0.0f
, 0.0f
, 1.0f
} }
310 for (int i
= 0; i
< state
.samplesLon
; ++i
)
312 for (int j
= 0; j
< state
.samplesLat
; ++j
)
315 bool accurate
= true;
316 double alt
= altitudeAt(state
, i
, j
, &coord
, &accurate
);
317 for (int n
= 0; n
< 3; ++n
)
319 GLvec3f norm
= normalAt(n
, coord
);
324 double rad
= m_meanRadius
+ alt
;
325 glColor3fv(normColours
[n
][0]);
326 glVertex3(dir
* rad
);
327 glColor3fv(normColours
[n
][1]);
328 glVertex3(dir
* rad
+ norm
*50.0f
);
338 // Render the solid rock walls
339 glDisable(GL_CULL_FACE
);
344 glEnable(GL_CULL_FACE
);
348 /// @todo cache one edge of strip to save time on other
349 for (int i
= 0; i
< state
.samplesLon
-1; ++i
)
351 glBegin(GL_TRIANGLE_STRIP
);
353 for (int j
= 0; j
< state
.samplesLat
; ++j
)
355 for (int k
= 0; k
< 2; ++k
)
358 bool accurate
= true;
359 double alt
= altitudeAt(state
, i
+k
, j
, &coord
, &accurate
);
362 foreach (tcGeoImageData
* imagery
, m_imagery
)
364 imagery
->texCoord(coord
);
369 glColor4f(0.5f
, 0.5f
, 1.0f
, 1.0f
);
373 glColor4f(1.0f
, 1.0f
, 1.0f
, 1.0f
);
374 //glColor4f(1.0f, 1.0f-(float)alt/3278.0f, 0.0f, 1.0f);
376 // Colour code if applicable
377 if (m_colourCoding
== ElevationSampleAlignment
)
379 if (state
.moreAvailableLon
&& state
.moreAvailableLat
)
381 glColor3f(1.0f
, 0.0f
, 1.0f
);
383 else if (state
.moreAvailableLon
)
385 glColor3f(1.0f
, 0.0f
, 0.0f
);
387 else if (state
.moreAvailableLat
)
389 glColor3f(0.0f
, 0.0f
, 1.0f
);
393 glColor3f(0.0f
, 1.0f
, 0.0f
);
397 double rad
= m_meanRadius
+ alt
;
398 glVertex3(dir
* rad
);
407 /// Render from the POV of an observer.
408 void tcGlobe::render(tcObserver
* const observer
, bool adaptive
, const tcGeo
* extent
)
410 /// @todo use a really simple fragment shader to cull backfacing lines
413 glColor3f(0.0f
, 1.0f
, 0.0f
);
414 for (int lat
= -75; lat
<= 75; lat
+= 15)
418 drawLineOfLatitude(M_PI
/180*lat
);
422 double tropic
= (23.0 + 26.0/60 + 22.0/3600) * M_PI
/180;
424 glColor3f(1.0f
, 0.0f
, 0.0f
);
425 drawLineOfLatitude(0.0);
426 // Tropics (Capricorn and Cancer)
427 glColor3f(1.0f
, 0.0f
, 1.0f
);
428 drawLineOfLatitude(-tropic
);
429 glColor3f(1.0f
, 1.0f
, 0.0f
);
430 drawLineOfLatitude(+tropic
);
431 // Arctic and Antarctic Circles
432 glColor3f(1.0f
, 1.0f
, 1.0f
);
433 drawLineOfLatitude(+M_PI
/2 - tropic
);
434 drawLineOfLatitude(-M_PI
/2 + tropic
);
436 // Lines of longitude
437 for (int lon
= 0; lon
< 360; lon
+= 15)
439 double x
= sin(M_PI
/180*lon
) * m_meanRadius
;
440 double y
= -cos(M_PI
/180*lon
) * m_meanRadius
;
446 glColor3f(1.0f
, lon
/180, 0.0f
);
453 glBegin(GL_LINE_STRIP
);
455 for (int lat
= minLat
; lat
<= maxLat
; ++lat
)
457 double z
= cos(M_PI
/180*lat
) * m_meanRadius
;
458 double xy
= sin(M_PI
/180*lat
);
459 glVertex3d(xy
*x
, xy
*y
, z
);
466 glColor3f(0.0f
, 1.0f
, 0.0f
);
471 // Draw data diagramatically
472 foreach (tcGeoImageData
* imagery
, m_imagery
)
474 imagery
->renderSchematic(m_meanRadius
, observer
);
479 int colourMapping
[6];
480 for (int band
= 0; band
< m_imagery
.size(); ++band
)
482 for (int i
= 0; i
< 6; ++i
)
484 colourMapping
[i
] = (m_colourMapping
[i
][0] == band
? m_colourMapping
[i
][1] : -1);
486 tcGeoImageData
* imagery
= m_imagery
[band
];
487 imagery
->setupThumbnailRendering(6, colourMapping
);
492 for (int lon
= -180; lon
< 180; lon
+= dlon
)
494 for (int lat
= -90; lat
< 90; lat
+= dlat
)
496 tcGeo
sw(M_PI
/180 * (lon
), M_PI
/180 * (lat
));
497 tcGeo
ne(M_PI
/180 * (lon
+dlon
), M_PI
/180 * (lat
+dlat
));
498 renderCell(observer
, sw
, ne
, adaptive
? 5 : 0, false);
501 foreach (tcGeoImageData
* imagery
, m_imagery
)
503 imagery
->setupNormalRendering();
505 for (int lon
= -180; lon
< 180; lon
+= dlon
)
507 for (int lat
= -90; lat
< 90; lat
+= dlat
)
509 tcGeo
sw(M_PI
/180 * (lon
), M_PI
/180 * (lat
));
510 tcGeo
ne(M_PI
/180 * (lon
+dlon
), M_PI
/180 * (lat
+dlat
));
511 renderCell(observer
, sw
, ne
, adaptive
? 5 : 0, true);
517 tcGeo sw
= extent
[0];
518 tcGeo ne
= extent
[1];
519 if (sw
.lon() > ne
.lon())
522 ne
.setLon(extent
[0].lon());
524 if (sw
.lat() > ne
.lat())
527 ne
.setLat(extent
[0].lat());
529 int colourMapping
[6];
530 for (int band
= 0; band
< m_imagery
.size(); ++band
)
532 for (int i
= 0; i
< 6; ++i
)
534 colourMapping
[i
] = (m_colourMapping
[i
][0] == band
? m_colourMapping
[i
][1] : -1);
536 tcGeoImageData
* imagery
= m_imagery
[band
];
537 imagery
->setupDetailedRendering(6, colourMapping
, sw
, ne
);
539 /// @todo If it is really big, split it
540 renderCell(observer
, sw
, ne
, adaptive
? 16 : 0, false, true, true, true, true);
541 foreach (tcGeoImageData
* imagery
, m_imagery
)
543 imagery
->setupNormalRendering();
545 renderCell(observer
, sw
, ne
, adaptive
? 16 : 0, true, true, true, true, true);
547 foreach (tcGeoImageData
* imagery
, m_imagery
)
549 imagery
->finishRendering();
553 /// Set the elevation mode to render in.
554 void tcGlobe::setElevationMode(int dem
, ElevationMode mode
)
556 Q_ASSERT(dem
>= 0 && dem
< 2);
557 m_elevationMode
[dem
] = mode
;
560 /// Set the level of interpolation.
561 void tcGlobe::setElevationInterpolation(float interpolation
)
563 m_elevationInterpolation
= interpolation
;
566 /// Set the elevation data set name.
567 void tcGlobe::setElevationDataSet(int dem
, const QString
& name
)
569 Q_ASSERT(dem
>= 0 && dem
< 2);
570 m_elevation
[dem
].setDataSet(name
);
573 /// Set colour coding method.
574 void tcGlobe::setColourCoding(ColourCoding colourCoding
)
576 m_colourCoding
= colourCoding
;
579 /// Adjust the mapping between bands and colour channels.
580 void tcGlobe::setColourMapping(int outputChannel
, int inputBand
, int inputGroup
)
582 if (outputChannel
< 6)
584 m_colourMapping
[outputChannel
][0] = inputGroup
;
585 m_colourMapping
[outputChannel
][1] = inputBand
;
593 /// Get the mean radius.
594 double tcGlobe::meanRadius() const
599 /// Get the altitude above sea level at a sample in a render state.
600 double tcGlobe::altitudeAt(const tcSrtmModel::RenderState
& state
, int x
, int y
, tcGeo
* outCoord
, bool* isAccurate
) const
611 m_elevationInterpolation
< 1.0f
,
612 m_elevationInterpolation
> 0.0f
614 for (int i
= 0; i
< 2; ++i
)
616 switch (m_elevationMode
[i
])
622 alt
[i
] = m_elevation
[i
].altitudeAt(state
, x
, y
, outCoord
, false, &accurate
[i
]);
629 case CorrectedElevation
:
632 alt
[i
] = m_elevation
[i
].altitudeAt(state
, x
, y
, outCoord
, true, &accurate
[i
]);
640 double result
= alt
[0]*(1.0-m_elevationInterpolation
) + alt
[1]*m_elevationInterpolation
;
643 *isAccurate
= accurate
[0] && accurate
[1];
648 /// Get the altitude above sea level at a coordinate.
649 double tcGlobe::altitudeAt(const tcGeo
& coord
, bool* isAccurate
) const
660 m_elevationInterpolation
< 1.0f
,
661 m_elevationInterpolation
> 0.0f
663 for (int i
= 0; i
< 2; ++i
)
665 switch (m_elevationMode
[i
])
671 alt
[i
] = m_elevation
[i
].altitudeAt(coord
, false, &accurate
[i
]);
678 case CorrectedElevation
:
681 alt
[i
] = m_elevation
[i
].altitudeAt(coord
, true, &accurate
[i
]);
689 double result
= alt
[0]*(1.0-m_elevationInterpolation
) + alt
[1]*m_elevationInterpolation
;
692 *isAccurate
= accurate
[0] && accurate
[1];
697 /// Get the radius at a coordinate.
698 double tcGlobe::radiusAt(const tcGeo
& coord
) const
700 return m_meanRadius
+ altitudeAt(coord
);
703 /// Get the texture coordinate of the effective texture at a geographical coordinate.
704 maths::Vector
<2,double> tcGlobe::textureCoordOfGeo(const tcGeo
& coord
) const
706 /// @todo Reimplement tcGlobe::textureCoordOfGeo with multiple imagery
707 Q_ASSERT(0 && "Reimplement tcGlobe::TextureCoordOfGeo with multiple imagery");
708 return m_imagery
[0]->geoToEffectiveTex() * coord
;
711 /// Get the current normal at a coordinate.
712 maths::Vector
<3,float> tcGlobe::normalAt(int norm
, const tcGeo
& coord
) const
714 Q_ASSERT(norm
>= 0 && norm
< 3);
715 if (m_colourMapping
[3+norm
][0] >= 0)
717 return m_imagery
[m_colourMapping
[3+norm
][0]]->normalAt(norm
, coord
);
721 return maths::Vector
<3,float>(0.0f
);