Kstars

terrainrenderer.cpp
1/*
2 SPDX-FileCopyrightText: 2021 Hy Murveit <hy@murveit.com>
3
4 SPDX-License-Identifier: GPL-2.0-or-later
5*/
6
7#include "terrainrenderer.h"
8
9#include "kstars_debug.h"
10#include "Options.h"
11#include "skymap.h"
12#include "skyqpainter.h"
13#include "projections/projector.h"
14#include "projections/equirectangularprojector.h"
15#include "skypoint.h"
16#include "kstars.h"
17
18#include <QStatusBar>
19
20// This is the factory that builds the one-and-only TerrainRenderer.
21TerrainRenderer * TerrainRenderer::_terrainRenderer = nullptr;
22TerrainRenderer *TerrainRenderer::Instance()
23{
24 if (_terrainRenderer == nullptr)
25 _terrainRenderer = new TerrainRenderer();
26
27 return _terrainRenderer;
28}
29
30// This class implements a quick 2D float array.
31class TerrainLookup
32{
33 public:
34 TerrainLookup(int width, int height) :
35 valPtr(new float[width * height]), valWidth(width)
36 {
37 memset(valPtr, 0, width * height * sizeof(float));
38 }
39 ~TerrainLookup()
40 {
41 delete[] valPtr;
42 }
43 inline float get(int w, int h)
44 {
45 return valPtr[h * valWidth + w];
46 }
47 inline void set(int w, int h, float val)
48 {
49 valPtr[h * valWidth + w] = val;
50 }
51 private:
52 float *valPtr;
53 int valWidth = 0;
54};
55
56// Samples 2-D array and returns interpolated values for the unsampled elements.
57// Used to speed up the calculation of azimuth and altitude values for all pixels
58// of the array to be rendered. Instead we compute every nth pixel's coordinates (e.g. n=4)
59// and interpolate the azimuth and alitutde values (so 1/16 the number of calculations).
60// This should be more than accurate enough for this rendering, and this part
61// of the code definitely needs speed.
62class InterpArray
63{
64 public:
65 // Constructor calculates the downsampled size and allocates the 2D arrays
66 // for azimuth and altitude.
67 InterpArray(int width, int height, int samplingFactor) : sampling(samplingFactor)
68 {
69 int downsampledWidth = width / sampling;
70 if (width % sampling != 0)
71 downsampledWidth += 1;
72 lastDownsampledCol = downsampledWidth - 1;
73
74 int downsampledHeight = height / sampling;
75 if (height % sampling != 0)
76 downsampledHeight += 1;
77 lastDownsampledRow = downsampledHeight - 1;
78
79 azLookup = new TerrainLookup(downsampledWidth, downsampledHeight);
80 altLookup = new TerrainLookup(downsampledWidth, downsampledHeight);
81 }
82
83 ~InterpArray()
84 {
85 delete azLookup;
86 delete altLookup;
87 }
88
89 // Get the azimuth and altitude values from the 2D arrays.
90 // Inputs are a full-image position
91 inline void get(int x, int y, float *az, float *alt)
92 {
93 const bool rowSampled = y % sampling == 0;
94 const bool colSampled = x % sampling == 0;
95 const int xSampled = x / sampling;
96 const int ySampled = y / sampling;
97
98 // If the requested position has been calculated, the use the stored value.
99 if (rowSampled && colSampled)
100 {
101 *az = azLookup->get(xSampled, ySampled);
102 *alt = altLookup->get(xSampled, ySampled);
103 return;
104 }
105 // The desired position has not been calculated, but values on its row have been calculated.
106 // Interpolate between the nearest values on the row.
107 else if (rowSampled)
108 {
109 if (xSampled >= lastDownsampledCol)
110 {
111 // If the requested value is on or past the last calculated value,
112 // just use that last value.
113 *az = azLookup->get(xSampled, ySampled);
114 *alt = altLookup->get(xSampled, ySampled);
115 return;
116 }
117 // Weight the contributions of the 2 calculated values on the row by their distance to the desired position.
118 const float weight2 = static_cast<float>(x) / sampling - xSampled;
119 const float weight1 = 1 - weight2;
120 *az = weight1 * azLookup->get(xSampled, ySampled) + weight2 * azLookup->get(xSampled + 1, ySampled);
121 *alt = weight1 * altLookup->get(xSampled, ySampled) + weight2 * altLookup->get(xSampled + 1, ySampled);
122 return;
123 }
124 // The desired position has not been calculated, but values on its column have been calculated.
125 // Interpolate between the nearest values on the column.
126 else if (colSampled)
127 {
128 if (ySampled >= lastDownsampledRow)
129 {
130 // If the requested value is on or past the last calculated value,
131 // just use that last value.
132 *az = azLookup->get(xSampled, ySampled);
133 *alt = altLookup->get(xSampled, ySampled);
134 return;
135 }
136 // Weight the contributions of the 2 calculated values on the column by their distance to the desired position.
137 const float weight2 = static_cast<float>(y) / sampling - ySampled;
138 const float weight1 = 1 - weight2;
139 *az = weight1 * azLookup->get(xSampled, ySampled) + weight2 * azLookup->get(xSampled, ySampled + 1);
140 *alt = weight1 * altLookup->get(xSampled, ySampled) + weight2 * altLookup->get(xSampled, ySampled + 1);
141 return;
142 }
143 else
144 {
145 // The desired position has not been calculated, and no values on its column nor row have been calculated.
146 // Interpolate between the nearest 4 values.
147 if (xSampled >= lastDownsampledCol || ySampled >= lastDownsampledRow)
148 {
149 // If we're near the edge, just use a close value. Harder to interpolate and not too important.
150 *az = azLookup->get(xSampled, ySampled);
151 *alt = altLookup->get(xSampled, ySampled);
152 return;
153 }
154 // The section uses distance along the two nearest calculated rows to come up with an estimate.
155 const float weight2 = static_cast<float>(x) / sampling - xSampled;
156 const float weight1 = 1 - weight2;
157 const float azWval = (weight1 * ( azLookup->get(xSampled, ySampled) + azLookup->get(xSampled, ySampled + 1)) +
158 weight2 * ( azLookup->get(xSampled + 1, ySampled) + azLookup->get(xSampled + 1, ySampled + 1)));
159 const float altWval = (weight1 * (altLookup->get(xSampled, ySampled) + altLookup->get(xSampled, ySampled + 1)) +
160 weight2 * (altLookup->get(xSampled + 1, ySampled) + altLookup->get(xSampled + 1, ySampled + 1)));
161
162 // The section uses distance along the two nearest calculated columns to come up with an estimate.
163 const float hweight2 = static_cast<float>(y) / sampling - ySampled;
164 const float hweight1 = 1 - hweight2;
165 const float azHval = (hweight2 * ( azLookup->get(xSampled, ySampled) + azLookup->get(xSampled + 1, ySampled)) +
166 hweight1 * ( azLookup->get(xSampled, ySampled + 1) + azLookup->get(xSampled + 1, ySampled + 1)));
167 const float altHval = (hweight2 * (altLookup->get(xSampled, ySampled) + altLookup->get(xSampled + 1, ySampled)) +
168 hweight1 * (altLookup->get(xSampled, ySampled + 1) + altLookup->get(xSampled + 1, ySampled + 1)));
169 // We average the 2 estimates (note above actually estimated twice the value).
170 *az = (azWval + azHval) / 4.0;
171 *alt = (altWval + altHval) / 4.0;
172 return;
173 }
174 }
175 TerrainLookup *azimuthLookup()
176 {
177 return azLookup;
178 }
179 TerrainLookup *altitudeLookup()
180 {
181 return altLookup;
182 }
183
184 private:
185 // These are the indeces of the last rows and columns (in the downsampled sized space) that were filled.
186 // This is needed because the downsample factor might not be an even multiple of the image size.
187 int lastDownsampledCol = 0;
188 int lastDownsampledRow = 0;
189 // The downsample factor.
190 int sampling = 0;
191 // The azimuth and altitude values are stored in these 2D arrays.
192 TerrainLookup *azLookup = nullptr;
193 TerrainLookup *altLookup = nullptr;
194};
195
196TerrainRenderer::TerrainRenderer()
197{
198}
199
200// Put degrees in the range of 0 -> 359.99999999
201double rationalizeAz(double degrees)
202{
203 if (degrees < -1000 || degrees > 1000)
204 return 0;
205
206 while (degrees < 0)
207 degrees += 360.0;
208 while (degrees >= 360.0)
209 degrees -= 360.0;
210 return degrees;
211}
212
213// Checks that degrees in the range of -90 -> 90.
214double rationalizeAlt(double degrees)
215{
216 if (degrees > 90.0)
217 return 90.0;
218 if (degrees < -90)
219 return -90;
220 return degrees;
221}
222
223// Assumes the source photosphere has rows which, left-to-right go from AZ=0 to AZ=360
224// and columns go from -90 altitude on the bottom to +90 on top.
225// Returns the pixel for the desired azimuth and altitude.
226QRgb TerrainRenderer::getPixel(double az, double alt) const
227{
228 az = rationalizeAz(az + Options::terrainSourceCorrectAz());
229 // This may make alt > 90 (due to a negative sourceCorrectAlt).
230 // If so, it returns 0, which is a transparent pixel.
231 alt = alt - Options::terrainSourceCorrectAlt();
232 if (az < 0 || az >= 360 || alt < -90 || alt > 90)
233 return(0);
234
235 // shift az to be -180 to 180
236 if (az > 180)
237 az = az - 360.0;
238 const int width = sourceImage.width();
239 const int height = sourceImage.height();
240
241 if (!Options::terrainSmoothPixels())
242 {
243 // az=0 should be the middle of the image.
244 int pixX = width / 2 + (az / 360.0) * width;
245 if (pixX > width - 1)
246 pixX = width - 1;
247 else if (pixX < 0)
248 pixX = 0;
249 int pixY = ((alt + 90.0) / 180.0) * height;
250 if (pixY > height - 1)
251 pixY = height - 1;
252 pixY = (height - 1) - pixY;
253 return sourceImage.pixel(pixX, pixY);
254 }
255
256 // Get floating point pixel positions so we can interpolate.
257 float pixX = width / 2 + (az / 360.0) * width;
258 if (pixX > width - 1)
259 pixX = width - 1;
260 else if (pixX < 0)
261 pixX = 0;
262 float pixY = ((alt + 90.0) / 180.0) * height;
263 if (pixY > height - 1)
264 pixY = height - 1;
265 pixY = (height - 1) - pixY;
266
267 // Don't bother interpolating for transparent pixels.
268 constexpr int lowAlpha = 0.1 * 255;
269 if (qAlpha(sourceImage.pixel(pixX, pixY)) < lowAlpha)
270 return sourceImage.pixel(pixX, pixY);
271
272 // Instead of just returning the pixel at the truncated position as above,
273 // below we interpolate the pixel RGBA values based on the floating-point pixel position.
274 int x1 = static_cast<int>(pixX);
275 int y1 = static_cast<int>(pixY);
276
277 if ((x1 >= width - 1) || (y1 >= height - 1))
278 return sourceImage.pixel(x1, y1);
279
280 // weights for the x & x+1, and y & y+1 positions.
281 float wx2 = pixX - x1;
282 float wx1 = 1.0 - wx2;
283 float wy2 = pixY - y1;
284 float wy1 = 1.0 - wy2;
285
286 // The pixels we'll interpolate.
287 QRgb c11(qUnpremultiply(sourceImage.pixel(pixX, pixY)));
288 QRgb c12(qUnpremultiply(sourceImage.pixel(pixX, pixY + 1)));
289 QRgb c21(qUnpremultiply(sourceImage.pixel(pixX + 1, pixY)));
290 QRgb c22(qUnpremultiply(sourceImage.pixel(pixX + 1, pixY + 1)));
291
292 // Weights for the above pixels.
293 float w11 = wx1 * wy1;
294 float w12 = wx1 * wy2;
295 float w21 = wx2 * wy1;
296 float w22 = wx2 * wy2;
297
298 // Finally, interpolate each component.
299 int red = w11 * qRed(c11) + w12 * qRed(c12) + w21 * qRed(c21) + w22 * qRed(c22);
300 int green = w11 * qGreen(c11) + w12 * qGreen(c12) + w21 * qGreen(c21) + w22 * qGreen(c22);
301 int blue = w11 * qBlue(c11) + w12 * qBlue(c12) + w21 * qBlue(c21) + w22 * qBlue(c22);
302 int alpha = w11 * qAlpha(c11) + w12 * qAlpha(c12) + w21 * qAlpha(c21) + w22 * qAlpha(c22);
303 return qPremultiply(qRgba(red, green, blue, alpha));
304}
305
306// Checks to see if the view is the same as the last call to render.
307// If true, render (below) will skip its computations and return the same image
308// as was previously calculated.
309// If the view is not the same, this method stores away the details of the current view
310// so that it may make this comparison again in the future.
311bool TerrainRenderer::sameView(const Projector *proj, bool forceRefresh)
312{
313 ViewParams view = proj->viewParams();
314 SkyPoint point = *(view.focus);
315 const auto &lst = KStarsData::Instance()->lst();
316 const auto &lat = KStarsData::Instance()->geo()->lat();
317
318 // We compare az and alt, not RA and DEC, as the RA changes,
319 // but the rendering is tied to azimuth and altitude which may not change.
320 // Thus we convert point to horizontal coordinates.
321 point.EquatorialToHorizontal(lst, lat);
322 const double az = rationalizeAz(point.az().Degrees());
323 const double alt = rationalizeAlt(point.alt().Degrees());
324
325 bool ok = view.width == savedViewParams.width &&
326 view.height == savedViewParams.height &&
327 view.zoomFactor == savedViewParams.zoomFactor &&
328 view.rotationAngle == savedViewParams.rotationAngle &&
329 view.useRefraction == savedViewParams.useRefraction &&
330 view.useAltAz == savedViewParams.useAltAz &&
331 view.fillGround == savedViewParams.fillGround;
332 const double azDiff = fabs(savedAz - az);
333 const double altDiff = fabs(savedAlt - alt);
334 if (!forceRefresh && ok && azDiff < .0001 && altDiff < .0001)
335 return true;
336
337 // Store the view
338 savedViewParams = view;
339 savedViewParams.focus = nullptr;
340 savedAz = az;
341 savedAlt = alt;
342 return false;
343}
344
345bool TerrainRenderer::render(uint16_t w, uint16_t h, QImage *terrainImage, const Projector *proj)
346{
347 // This is used to force a re-render, e.g. when the image is changed.
348 bool dirty = false;
349
350 // Load the image if necessary.
351 QString filename = Options::terrainSource();
352 if (!initialized || filename != sourceFilename)
353 {
354 dirty = true;
355 QImage image;
356 if (image.load(filename))
357 {
359 qCDebug(KSTARS) << QString("Read terrain file %1 x %2").arg(sourceImage.width()).arg(sourceImage.height());
360 sourceFilename = filename;
361 initialized = true;
362 }
363 else
364 {
365 if (filename.isEmpty())
366 KStars::Instance()->statusBar()->showMessage(i18n("Failed to load terrain. Set terrain file in Settings."));
367 else
368 KStars::Instance()->statusBar()->showMessage(i18n("Failed to load terrain image (%1). Set terrain file in Settings.",
369 filename));
370 initialized = false;
371 Options::setShowTerrain(false);
373 }
374 }
375
376 if (!initialized)
377 return false;
378
379 // Check if parameters have changed.
380 if ((terrainDownsampling != Options::terrainDownsampling()) ||
381 (terrainSmoothPixels != Options::terrainSmoothPixels()) ||
382 (terrainSkipSpeedup != Options::terrainSkipSpeedup()) ||
383 (terrainTransparencySpeedup != Options::terrainTransparencySpeedup()) ||
384 (terrainSourceCorrectAz != Options::terrainSourceCorrectAz()) ||
385 (terrainSourceCorrectAlt != Options::terrainSourceCorrectAlt()))
386 dirty = true;
387
388 terrainDownsampling = Options::terrainDownsampling();
389 terrainSmoothPixels = Options::terrainSmoothPixels();
390 terrainSkipSpeedup = Options::terrainSkipSpeedup();
391 terrainTransparencySpeedup = Options::terrainTransparencySpeedup();
392 terrainSourceCorrectAz = Options::terrainSourceCorrectAz();
393 terrainSourceCorrectAlt = Options::terrainSourceCorrectAlt();
394
395 if (sameView(proj, dirty))
396 {
397 // Just return the previous image if the input view hasn't changed.
398 *terrainImage = savedImage.copy();
399 return true;
400 }
401
402 QElapsedTimer timer;
403 timer.start();
404
405 // Only compute the pixel's az and alt values for every Nth pixel.
406 // Get the other pixel az and alt values by interpolation.
407 // This saves a lot of time.
408 const int sampling = Options::terrainDownsampling();
409 InterpArray interp(w, h, sampling);
410 QElapsedTimer setupTimer;
411 setupTimer.start();
412 setupLookup(w, h, sampling, proj, interp.azimuthLookup(), interp.altitudeLookup());
413
414 const double setupTime = setupTimer.elapsed() / 1000.0; ///////////////////
415
416 // Another speedup. If true, our calculations are downsampled by 2 in each dimension.
417 const bool skip = Options::terrainSkipSpeedup() || SkyMap::IsSlewing();
418 int increment = skip ? 2 : 1;
419
420 // Assign transparent pixels everywhere by default.
421 terrainImage->fill(0);
422
423 // Go through the image, and for each pixel, using the previously computed az and alt values
424 // get the corresponding pixel from the terrain image.
425 bool lastTransparent = false;
426 for (int j = 0; j < h; j += increment)
427 {
428 bool notLastRow = j != h - 1;
429 for (int i = 0; i < w; i += increment)
430 {
431 if (lastTransparent && Options::terrainTransparencySpeedup())
432 {
433 // Speedup--if the last pixel was transparent, then this
434 // one is assumed transparent too (but next is calculated).
435 lastTransparent = false;
436 continue;
437 }
438
439 const QPointF imgPoint(i, j);
440 bool equiRectangular = (proj->type() == Projector::Equirectangular);
441 bool usable = equiRectangular ? !dynamic_cast<const EquirectangularProjector*>(proj)->unusablePoint(imgPoint)
442 : !proj->unusablePoint(imgPoint);
443 if (usable)
444 {
445 float az, alt;
446 interp.get(i, j, &az, &alt);
447 const QRgb pixel = getPixel(az, alt);
448 terrainImage->setPixel(i, j, pixel);
449 lastTransparent = (pixel == 0);
450
451 if (skip)
452 {
453 // If we've skipped, fill in the missing pixels.
454 bool notLastCol = i != w - 1;
455 if (notLastCol)
456 terrainImage->setPixel(i + 1, j, pixel);
457 if (notLastRow)
458 terrainImage->setPixel(i, j + 1, pixel);
459 if (notLastRow && notLastCol)
460 terrainImage->setPixel(i + 1, j + 1, pixel);
461 }
462 }
463 // Otherwise terrainImage was already filled with transparent pixels
464 // so i,j will be transparent.
465 }
466 }
467
468 savedImage = terrainImage->copy();
469
470 QFile f(sourceFilename);
471 QFileInfo fileInfo(f.fileName());
472 QString fName(fileInfo.fileName());
473 QString dbgMsg(QString("Terrain rendering: %1px, %2s (%3s) %4 ds %5 skip %6 trnsp %7 pan %8 smooth %9")
474 .arg(w * h)
475 .arg(timer.elapsed() / 1000.0, 5, 'f', 3)
476 .arg(setupTime, 5, 'f', 3)
477 .arg(fName)
478 .arg(Options::terrainDownsampling())
479 .arg(Options::terrainSkipSpeedup() ? "T" : "F")
480 .arg(Options::terrainTransparencySpeedup() ? "T" : "F")
481 .arg(Options::terrainPanning() ? "T" : "F")
482 .arg(Options::terrainSmoothPixels() ? "T" : "F"));
483 //qCDebug(KSTARS) << dbgMsg;
484 //fprintf(stderr, "%s\n", dbgMsg.toLatin1().data());
485
486 dirty = false;
487 return true;
488}
489
490// Goes through every Nth input pixel position, finding their azimuth and altitude
491// and storing that for future use in the interpolations above.
492// This is the most time-costly part of the computation.
493void TerrainRenderer::setupLookup(uint16_t w, uint16_t h, int sampling, const Projector *proj, TerrainLookup *azLookup,
494 TerrainLookup *altLookup)
495{
496 const auto &lst = KStarsData::Instance()->lst();
497 const auto &lat = KStarsData::Instance()->geo()->lat();
498 for (int j = 0, js = 0; j < h; j += sampling, js++)
499 {
500 for (int i = 0, is = 0; i < w; i += sampling, is++)
501 {
502 const QPointF imgPoint(i, j);
503 bool equiRectangular = (proj->type() == Projector::Equirectangular);
504 bool usable = equiRectangular ? !dynamic_cast<const EquirectangularProjector*>(proj)->unusablePoint(imgPoint)
505 : !proj->unusablePoint(imgPoint);
506 if (usable)
507 {
508 SkyPoint point = equiRectangular ?
509 dynamic_cast<const EquirectangularProjector*>(proj)->fromScreen(imgPoint, lst, lat, true)
510 : proj->fromScreen(imgPoint, lst, lat, true);
511 const double az = rationalizeAz(point.az().Degrees());
512 const double alt = rationalizeAlt(point.alt().Degrees());
513 azLookup->set(is, js, az);
514 altLookup->set(is, js, alt);
515 }
516 }
517 }
518}
519
Implememntation of Equirectangular projection
const CachingDms * lat() const
Definition geolocation.h:70
CachingDms * lst()
Definition kstarsdata.h:224
GeoLocation * geo()
Definition kstarsdata.h:230
static KStars * Instance()
Definition kstars.h:121
void syncOps()
Sync Options to GUI, if any.
The Projector class is the primary class that serves as an interface to handle projections.
Definition projector.h:58
virtual bool unusablePoint(const QPointF &p) const
Check if the current point on screen is a valid point on the sky.
virtual Q_INVOKABLE Projection type() const =0
Return the type of this projection.
The sky coordinates of a point in the sky.
Definition skypoint.h:45
void EquatorialToHorizontal(const CachingDms *LST, const CachingDms *lat)
Determine the (Altitude, Azimuth) coordinates of the SkyPoint from its (RA, Dec) coordinates,...
Definition skypoint.cpp:77
const dms & az() const
Definition skypoint.h:275
const dms & alt() const
Definition skypoint.h:281
This is just a container that holds information needed to do projections.
Definition projector.h:37
bool fillGround
If the ground is filled, then points below horizon are invisible.
Definition projector.h:44
const double & Degrees() const
Definition dms.h:141
QString i18n(const char *text, const TYPE &arg...)
qint64 elapsed() const const
Format_ARGB32_Premultiplied
QImage convertToFormat(Format format, Qt::ImageConversionFlags flags) &&
QImage copy(const QRect &rectangle) const const
void fill(Qt::GlobalColor color)
int height() const const
bool load(QIODevice *device, const char *format)
QRgb pixel(const QPoint &position) const const
void setPixel(const QPoint &position, uint index_or_rgb)
int width() const const
QStatusBar * statusBar() const const
void showMessage(const QString &message, int timeout)
QString arg(Args &&... args) const const
bool isEmpty() const const
This file is part of the KDE documentation.
Documentation copyright © 1996-2024 The KDE developers.
Generated on Mon Nov 4 2024 16:38:44 by doxygen 1.12.0 written by Dimitri van Heesch, © 1997-2006

KDE's Doxygen guidelines are available online.