Kstars

terrainrenderer.cpp
1 /*
2  SPDX-FileCopyrightText: 2021 Hy Murveit <[email protected]>
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.
21 TerrainRenderer * TerrainRenderer::_terrainRenderer = nullptr;
22 TerrainRenderer *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.
31 class 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.
62 class 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 
196 TerrainRenderer::TerrainRenderer()
197 {
198 }
199 
200 // Put degrees in the range of 0 -> 359.99999999
201 double 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.
214 double 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.
226 QRgb 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.
311 bool 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.zoomFactor == savedViewParams.zoomFactor &&
327  view.useRefraction == savedViewParams.useRefraction &&
328  view.useAltAz == savedViewParams.useAltAz &&
329  view.fillGround == savedViewParams.fillGround;
330  const double azDiff = fabs(savedAz - az);
331  const double altDiff = fabs(savedAlt - alt);
332  if (!forceRefresh && ok && azDiff < .0001 && altDiff < .0001)
333  return true;
334 
335  // Store the view
336  savedViewParams = view;
337  savedViewParams.focus = nullptr;
338  savedAz = az;
339  savedAlt = alt;
340  return false;
341 }
342 
343 bool TerrainRenderer::render(uint16_t w, uint16_t h, QImage *terrainImage, const Projector *proj)
344 {
345  // This is used to force a re-render, e.g. when the image is changed.
346  bool dirty = false;
347 
348  // Load the image if necessary.
349  QString filename = Options::terrainSource();
350  if (!initialized || filename != sourceFilename)
351  {
352  dirty = true;
353  QImage image;
354  if (image.load(filename))
355  {
357  qCDebug(KSTARS) << QString("Read terrain file %1 x %2").arg(sourceImage.width()).arg(sourceImage.height());
358  sourceFilename = filename;
359  initialized = true;
360  }
361  else
362  {
363  if (filename.isEmpty())
364  KStars::Instance()->statusBar()->showMessage(i18n("Failed to load terrain. Set terrain file in Settings."));
365  else
366  KStars::Instance()->statusBar()->showMessage(i18n("Failed to load terrain image (%1). Set terrain file in Settings.",
367  filename));
368  initialized = false;
369  Options::setShowTerrain(false);
371  }
372  }
373 
374  if (!initialized)
375  return false;
376 
377  // Check if parameters have changed.
378  if ((terrainDownsampling != Options::terrainDownsampling()) ||
379  (terrainSmoothPixels != Options::terrainSmoothPixels()) ||
380  (terrainSkipSpeedup != Options::terrainSkipSpeedup()) ||
381  (terrainTransparencySpeedup != Options::terrainTransparencySpeedup()) ||
382  (terrainSourceCorrectAz != Options::terrainSourceCorrectAz()) ||
383  (terrainSourceCorrectAlt != Options::terrainSourceCorrectAlt()))
384  dirty = true;
385 
386  terrainDownsampling = Options::terrainDownsampling();
387  terrainSmoothPixels = Options::terrainSmoothPixels();
388  terrainSkipSpeedup = Options::terrainSkipSpeedup();
389  terrainTransparencySpeedup = Options::terrainTransparencySpeedup();
390  terrainSourceCorrectAz = Options::terrainSourceCorrectAz();
391  terrainSourceCorrectAlt = Options::terrainSourceCorrectAlt();
392 
393  if (sameView(proj, dirty))
394  {
395  // Just return the previous image if the input view hasn't changed.
396  *terrainImage = savedImage.copy();
397  return true;
398  }
399 
400  QElapsedTimer timer;
401  timer.start();
402 
403  // Only compute the pixel's az and alt values for every Nth pixel.
404  // Get the other pixel az and alt values by interpolation.
405  // This saves a lot of time.
406  const int sampling = Options::terrainDownsampling();
407  InterpArray interp(w, h, sampling);
408  QElapsedTimer setupTimer;
409  setupTimer.start();
410  setupLookup(w, h, sampling, proj, interp.azimuthLookup(), interp.altitudeLookup());
411 
412  const double setupTime = setupTimer.elapsed() / 1000.0; ///////////////////
413 
414  // Another speedup. If true, our calculations are downsampled by 2 in each dimension.
415  const bool skip = Options::terrainSkipSpeedup() || SkyMap::IsSlewing();
416  int increment = skip ? 2 : 1;
417 
418  // Assign transparent pixels everywhere by default.
419  terrainImage->fill(0);
420 
421  // Go through the image, and for each pixel, using the previously computed az and alt values
422  // get the corresponding pixel from the terrain image.
423  bool lastTransparent = false;
424  for (int j = 0; j < h; j += increment)
425  {
426  bool notLastRow = j != h - 1;
427  for (int i = 0; i < w; i += increment)
428  {
429  if (lastTransparent && Options::terrainTransparencySpeedup())
430  {
431  // Speedup--if the last pixel was transparent, then this
432  // one is assumed transparent too (but next is calculated).
433  lastTransparent = false;
434  continue;
435  }
436 
437  const QPointF imgPoint(i, j);
438  bool equiRectangular = (proj->type() == Projector::Equirectangular);
439  bool usable = equiRectangular ? !dynamic_cast<const EquirectangularProjector*>(proj)->unusablePoint(imgPoint)
440  : !proj->unusablePoint(imgPoint);
441  if (usable)
442  {
443  float az, alt;
444  interp.get(i, j, &az, &alt);
445  const QRgb pixel = getPixel(az, alt);
446  terrainImage->setPixel(i, j, pixel);
447  lastTransparent = (pixel == 0);
448 
449  if (skip)
450  {
451  // If we've skipped, fill in the missing pixels.
452  bool notLastCol = i != w - 1;
453  if (notLastCol)
454  terrainImage->setPixel(i + 1, j, pixel);
455  if (notLastRow)
456  terrainImage->setPixel(i, j + 1, pixel);
457  if (notLastRow && notLastCol)
458  terrainImage->setPixel(i + 1, j + 1, pixel);
459  }
460  }
461  // Otherwise terrainImage was already filled with transparent pixels
462  // so i,j will be transparent.
463  }
464  }
465 
466  savedImage = terrainImage->copy();
467 
468  QFile f(sourceFilename);
469  QFileInfo fileInfo(f.fileName());
470  QString fName(fileInfo.fileName());
471  QString dbgMsg(QString("Terrain rendering: %1px, %2s (%3s) %4 ds %5 skip %6 trnsp %7 pan %8 smooth %9")
472  .arg(w * h)
473  .arg(timer.elapsed() / 1000.0, 5, 'f', 3)
474  .arg(setupTime, 5, 'f', 3)
475  .arg(fName)
476  .arg(Options::terrainDownsampling())
477  .arg(Options::terrainSkipSpeedup() ? "T" : "F")
478  .arg(Options::terrainTransparencySpeedup() ? "T" : "F")
479  .arg(Options::terrainPanning() ? "T" : "F")
480  .arg(Options::terrainSmoothPixels() ? "T" : "F"));
481  //qCDebug(KSTARS) << dbgMsg;
482  //fprintf(stderr, "%s\n", dbgMsg.toLatin1().data());
483 
484  dirty = false;
485  return true;
486 }
487 
488 // Goes through every Nth input pixel position, finding their azimuth and altitude
489 // and storing that for future use in the interpolations above.
490 // This is the most time-costly part of the computation.
491 void TerrainRenderer::setupLookup(uint16_t w, uint16_t h, int sampling, const Projector *proj, TerrainLookup *azLookup,
492  TerrainLookup *altLookup)
493 {
494  const auto &lst = KStarsData::Instance()->lst();
495  const auto &lat = KStarsData::Instance()->geo()->lat();
496  for (int j = 0, js = 0; j < h; j += sampling, js++)
497  {
498  for (int i = 0, is = 0; i < w; i += sampling, is++)
499  {
500  const QPointF imgPoint(i, j);
501  bool equiRectangular = (proj->type() == Projector::Equirectangular);
502  bool usable = equiRectangular ? !dynamic_cast<const EquirectangularProjector*>(proj)->unusablePoint(imgPoint)
503  : !proj->unusablePoint(imgPoint);
504  if (usable)
505  {
506  SkyPoint point = equiRectangular ?
507  dynamic_cast<const EquirectangularProjector*>(proj)->fromScreen(imgPoint, lst, lat, true)
508  : proj->fromScreen(imgPoint, lst, lat, true);
509  const double az = rationalizeAz(point.az().Degrees());
510  const double alt = rationalizeAlt(point.alt().Degrees());
511  azLookup->set(is, js, az);
512  altLookup->set(is, js, alt);
513  }
514  }
515  }
516 }
517 
const dms & alt() const
Definition: skypoint.h:281
Format_ARGB32_Premultiplied
void fill(uint pixelValue)
Stores dms coordinates for a point in the sky. for converting between coordinate systems.
Definition: skypoint.h:44
CachingDms * lst()
Definition: kstarsdata.h:223
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
bool load(QIODevice *device, const char *format)
virtual Q_INVOKABLE Projection type() const =0
Return the type of this projection.
static KStars * Instance()
Definition: kstars.h:125
void setPixel(int x, int y, uint index_or_rgb)
void showMessage(const QString &message, int timeout)
QString i18n(const char *text, const TYPE &arg...)
void syncOps()
Sync Options to GUI, if any.
const CachingDms * lat() const
Definition: geolocation.h:70
virtual bool unusablePoint(const QPointF &p) const
Check if the current point on screen is a valid point on the sky.
Definition: projector.cpp:390
bool isEmpty() const const
GeoLocation * geo()
Definition: kstarsdata.h:229
QImage convertToFormat(QImage::Format format, Qt::ImageConversionFlags flags) const &const
This is just a container that holds information needed to do projections.
Definition: projector.h:36
qint64 elapsed() const const
bool fillGround
If the ground is filled, then points below horizon are invisible.
Definition: projector.h:43
QStatusBar * statusBar() const const
QString arg(qlonglong a, int fieldWidth, int base, QChar fillChar) const const
const double & Degrees() const
Definition: dms.h:141
QImage copy(const QRect &rectangle) const const
virtual QVariant get(ScriptableExtension *callerPrincipal, quint64 objId, const QString &propName)
const dms & az() const
Definition: skypoint.h:275
This file is part of the KDE documentation.
Documentation copyright © 1996-2022 The KDE developers.
Generated on Fri Aug 12 2022 04:00:58 by doxygen 1.8.17 written by Dimitri van Heesch, © 1997-2006

KDE's Doxygen guidelines are available online.