1 /*
2 SPDX-FileCopyrightText: 2010 Henry de Valence <hdevalence@gmail.com>
3
4 SPDX-License-Identifier: GPL-2.0-or-later
5 */
6
7 #include "projector.h"
8
9 #include "ksutils.h"
10 #ifdef KSTARS_LITE
11 #include "skymaplite.h"
12 #endif
13 #include "skycomponents/skylabeler.h"
14
15 namespace
16 {
toXYZ(const SkyPoint * p,double * x,double * y,double * z)17 void toXYZ(const SkyPoint *p, double *x, double *y, double *z)
18 {
19 double sinRa, sinDec, cosRa, cosDec;
20
21 p->ra().SinCos(sinRa, cosRa);
22 p->dec().SinCos(sinDec, cosDec);
23 *x = cosDec * cosRa;
24 *y = cosDec * sinRa;
25 *z = sinDec;
26 }
27 }
28
pointAt(double az)29 SkyPoint Projector::pointAt(double az)
30 {
31 SkyPoint p;
32 p.setAz(az);
33 p.setAlt(0.0);
34 p.HorizontalToEquatorial(KStarsData::Instance()->lst(), KStarsData::Instance()->geo()->lat());
35 return p;
36 }
37
Projector(const ViewParams & p)38 Projector::Projector(const ViewParams &p)
39 {
40 m_data = KStarsData::Instance();
41 setViewParams(p);
42 // Force clip polygon update
43 updateClipPoly();
44 }
45
setViewParams(const ViewParams & p)46 void Projector::setViewParams(const ViewParams &p)
47 {
48 m_vp = p;
49
50 /** Precompute cached values */
51 //Find Sin/Cos for focus point
52 m_sinY0 = 0;
53 m_cosY0 = 0;
54 if (m_vp.useAltAz)
55 {
56 // N.B. We explicitly check useRefraction and not use
57 // SkyPoint::altRefracted() here because it could be different
58 // from Options::useRefraction() in some future use-case
59 // --asimha
60 SkyPoint::refract(m_vp.focus->alt(), m_vp.useRefraction).SinCos(m_sinY0, m_cosY0);
61 }
62 else
63 {
64 m_vp.focus->dec().SinCos(m_sinY0, m_cosY0);
65 }
66
67 double currentFOV = m_fov;
68 //Find FOV in radians
69 m_fov = sqrt(m_vp.width * m_vp.width + m_vp.height * m_vp.height) / (2 * m_vp.zoomFactor * dms::DegToRad);
70 //Set checkVisibility variables
71 double Ymax;
72 m_xrange = 1.2 * m_fov / m_cosY0;
73 if (m_vp.useAltAz)
74 {
75 Ymax = fabs(SkyPoint::refract(m_vp.focus->alt().Degrees(), m_vp.useRefraction)) + m_fov;
76 }
77 else
78 {
79 Ymax = fabs(m_vp.focus->dec().Degrees()) + m_fov;
80 }
81 m_isPoleVisible = (Ymax >= 90.0);
82
83 // Only update clipping polygon if there is an FOV change
84 if (currentFOV != m_fov)
85 updateClipPoly();
86 }
87
fov() const88 double Projector::fov() const
89 {
90 return m_fov;
91 }
92
toScreen(const SkyPoint * o,bool oRefract,bool * onVisibleHemisphere) const93 QPointF Projector::toScreen(const SkyPoint *o, bool oRefract, bool *onVisibleHemisphere) const
94 {
95 return KSUtils::vecToPoint(toScreenVec(o, oRefract, onVisibleHemisphere));
96 }
97
onScreen(const QPointF & p) const98 bool Projector::onScreen(const QPointF &p) const
99 {
100 return (0 <= p.x() && p.x() <= m_vp.width && 0 <= p.y() && p.y() <= m_vp.height);
101 }
102
onScreen(const Eigen::Vector2f & p) const103 bool Projector::onScreen(const Eigen::Vector2f &p) const
104 {
105 return onScreen(QPointF(p.x(), p.y()));
106 }
107
clipLine(SkyPoint * p1,SkyPoint * p2) const108 QPointF Projector::clipLine(SkyPoint *p1, SkyPoint *p2) const
109 {
110 return KSUtils::vecToPoint(clipLineVec(p1, p2));
111 }
112
clipLineVec(SkyPoint * p1,SkyPoint * p2) const113 Eigen::Vector2f Projector::clipLineVec(SkyPoint *p1, SkyPoint *p2) const
114 {
115 /* ASSUMES p1 was not clipped but p2 was.
116 * Return the QPoint that barely clips in the line twixt p1 and p2.
117 */
118 //TODO: iteration = ceil( 0.5*log2( w^2 + h^2) )??
119 // also possibly rewrite this
120 // --hdevalence
121 int iteration = 15; // For "perfect" clipping:
122 // 2^iterations should be >= max pixels/line
123 bool isVisible = true; // so we start at midpoint
124 SkyPoint mid;
125 Eigen::Vector2f oMid;
126 double x, y, z, dx, dy, dz, ra, dec;
127 int newx, newy, oldx, oldy;
128 oldx = oldy = -10000; // any old value that is not the first omid
129
130 toXYZ(p1, &x, &y, &z);
131 // -jbb printf("\np1: %6.4f %6.4f %6.4f\n", x, y, z);
132
133 toXYZ(p2, &dx, &dy, &dz);
134
135 // -jbb printf("p2: %6.4f %6.4f %6.4f\n", dx, dy, dz);
136 dx -= x;
137 dy -= y;
138 dz -= z;
139 // Successive approximation to point on line that just clips.
140 while (iteration-- > 0)
141 {
142 dx *= .5;
143 dy *= .5;
144 dz *= .5;
145 if (!isVisible) // move back toward visible p1
146 {
147 x -= dx;
148 y -= dy;
149 z -= dz;
150 }
151 else // move out toward clipped p2
152 {
153 x += dx;
154 y += dy;
155 z += dz;
156 }
157
158 // -jbb printf(" : %6.4f %6.4f %6.4f\n", x, y, z);
159 // [x, y, z] => [ra, dec]
160 ra = atan2(y, x);
161 dec = asin(z / sqrt(x * x + y * y + z * z));
162
163 mid = SkyPoint(ra * 12. / dms::PI, dec * 180. / dms::PI);
164 mid.EquatorialToHorizontal(m_data->lst(), m_data->geo()->lat());
165
166 oMid = toScreenVec(&mid, false, &isVisible);
167 //AND the result with checkVisibility to clip things going below horizon
168 isVisible &= checkVisibility(&mid);
169 newx = (int)oMid.x();
170 newy = (int)oMid.y();
171
172 // -jbb printf("new x/y: %4d %4d", newx, newy);
173 if ((oldx == newx) && (oldy == newy))
174 {
175 break;
176 }
177 oldx = newx;
178 oldy = newy;
179 }
180 return oMid;
181 }
182
checkVisibility(const SkyPoint * p) const183 bool Projector::checkVisibility(const SkyPoint *p) const
184 {
185 //TODO deal with alternate projections
186 //not clear how this depends on projection
187 //FIXME do these heuristics actually work?
188
189 double dX, dY;
190
191 //Skip objects below the horizon if the ground is drawn
192 /* Is the cost of this conversion actually less than drawing it anyways?
193 * EquatorialToHorizontal takes 3 SinCos calls -- so 6 trig calls if not using GNU exts.
194 */
195 /*
196 if( m_vp.fillGround ) {
197 if( !m_vp.useAltAz )
198 p->EquatorialToHorizontal( m_data->lst(), m_data->geo()->lat() );
199 if( p->alt().Degrees() < -1.0 ) return false;
200 }
201 */ //Here we hope that the point has already been 'synchronized'
202 if (m_vp.fillGround /*&& m_vp.useAltAz*/ && p->alt().Degrees() <= SkyPoint::altCrit)
203 return false;
204
205 if (m_vp.useAltAz)
206 {
207 /** To avoid calculating refraction, we just use the unrefracted
208 altitude and add a 2-degree 'safety factor' */
209 dY = fabs(p->alt().Degrees() - m_vp.focus->alt().Degrees()) - 2.;
210 }
211 else
212 {
213 dY = fabs(p->dec().Degrees() - m_vp.focus->dec().Degrees());
214 }
215 if (m_isPoleVisible)
216 dY *= 0.75; //increase effective FOV when pole visible.
217 if (dY > m_fov)
218 return false;
219 if (m_isPoleVisible)
220 return true;
221
222 if (m_vp.useAltAz)
223 {
224 dX = fabs(p->az().Degrees() - m_vp.focus->az().Degrees());
225 }
226 else
227 {
228 dX = fabs(p->ra().Degrees() - m_vp.focus->ra().Degrees());
229 }
230 if (dX > 180.0)
231 dX = 360.0 - dX; // take shorter distance around sky
232
233 return dX < m_xrange;
234 }
235
236 // FIXME: There should be a MUCH more efficient way to do this (see EyepieceField for example)
findNorthPA(const SkyPoint * o,float x,float y) const237 double Projector::findNorthPA(const SkyPoint *o, float x, float y) const
238 {
239 //Find position angle of North using a test point displaced to the north
240 //displace by 100/zoomFactor radians (so distance is always 100 pixels)
241 //this is 5730/zoomFactor degrees
242 KStarsData *data = KStarsData::Instance();
243 double newDec = o->dec().Degrees() + 5730.0 / m_vp.zoomFactor;
244 if (newDec > 90.0)
245 newDec = 90.0;
246 SkyPoint test(o->ra().Hours(), newDec);
247 if (m_vp.useAltAz)
248 test.EquatorialToHorizontal(data->lst(), data->geo()->lat());
249 Eigen::Vector2f t = toScreenVec(&test);
250 float dx = t.x() - x;
251 float dy = y - t.y(); //backwards because QWidget Y-axis increases to the bottom
252 float north;
253 if (dy)
254 {
255 north = atan2f(dx, dy) * 180.0 / dms::PI;
256 }
257 else
258 {
259 north = (dx > 0.0 ? -90.0 : 90.0);
260 }
261
262 return north;
263 }
264
findPA(const SkyObject * o,float x,float y) const265 double Projector::findPA(const SkyObject *o, float x, float y) const
266 {
267 return (findNorthPA(o, x, y) + o->pa());
268 }
269
groundPoly(SkyPoint * labelpoint,bool * drawLabel) const270 QVector<Eigen::Vector2f> Projector::groundPoly(SkyPoint *labelpoint, bool *drawLabel) const
271 {
272 QVector<Eigen::Vector2f> ground;
273
274 static const QString horizonLabel = i18n("Horizon");
275 float marginLeft, marginRight, marginTop, marginBot;
276 SkyLabeler::Instance()->getMargins(horizonLabel, &marginLeft, &marginRight, &marginTop, &marginBot);
277
278 //daz is 1/2 the width of the sky in degrees
279 double daz = 90.;
280 if (m_vp.useAltAz)
281 {
282 daz = 0.5 * m_vp.width * 57.3 / m_vp.zoomFactor; //center to edge, in degrees
283 if (type() == Projector::Orthographic)
284 {
285 daz = daz * 1.4;
286 }
287 daz = qMin(qreal(90.0), daz);
288 }
289
290 double faz = m_vp.focus->az().Degrees();
291 double az1 = faz - daz;
292 double az2 = faz + daz;
293
294 bool allGround = true;
295 bool allSky = true;
296
297 double inc = 1.0;
298 //Add points along horizon
299 for (double az = az1; az <= az2 + inc; az += inc)
300 {
301 SkyPoint p = pointAt(az);
302 bool visible = false;
303 Eigen::Vector2f o = toScreenVec(&p, false, &visible);
304 if (visible)
305 {
306 ground.append(o);
307 //Set the label point if this point is onscreen
308 if (labelpoint && o.x() < marginRight && o.y() > marginTop && o.y() < marginBot)
309 *labelpoint = p;
310
311 if (o.y() > 0.)
312 allGround = false;
313 if (o.y() < m_vp.height)
314 allSky = false;
315 }
316 }
317
318 if (allSky)
319 {
320 if (drawLabel)
321 *drawLabel = false;
322 return QVector<Eigen::Vector2f>();
323 }
324
325 if (allGround)
326 {
327 ground.clear();
328 ground.append(Eigen::Vector2f(-10., -10.));
329 ground.append(Eigen::Vector2f(m_vp.width + 10., -10.));
330 ground.append(Eigen::Vector2f(m_vp.width + 10., m_vp.height + 10.));
331 ground.append(Eigen::Vector2f(-10., m_vp.height + 10.));
332 if (drawLabel)
333 *drawLabel = false;
334 return ground;
335 }
336
337 //In Gnomonic projection, or if sufficiently zoomed in, we can complete
338 //the ground polygon by simply adding offscreen points
339 //FIXME: not just gnomonic
340 if (daz < 25.0 || type() == Projector::Gnomonic)
341 {
342 ground.append(Eigen::Vector2f(m_vp.width + 10.f, ground.last().y()));
343 ground.append(Eigen::Vector2f(m_vp.width + 10.f, m_vp.height + 10.f));
344 ground.append(Eigen::Vector2f(-10.f, m_vp.height + 10.f));
345 ground.append(Eigen::Vector2f(-10.f, ground.first().y()));
346 }
347 else
348 {
349 double r = m_vp.zoomFactor * radius();
350 double t1 =
351 atan2(-1. * (ground.last().y() - 0.5 * m_vp.height), ground.last().x() - 0.5 * m_vp.width) / dms::DegToRad;
352 double t2 = t1 - 180.;
353 for (double t = t1; t >= t2; t -= inc) //step along circumference
354 {
355 dms a(t);
356 double sa(0.), ca(0.);
357 a.SinCos(sa, ca);
358 ground.append(Eigen::Vector2f(0.5 * m_vp.width + r * ca, 0.5 * m_vp.height - r * sa));
359 }
360 }
361
362 if (drawLabel)
363 *drawLabel = true;
364 return ground;
365 }
366
updateClipPoly()367 void Projector::updateClipPoly()
368 {
369 m_clipPolygon.clear();
370
371 double r = m_vp.zoomFactor * radius();
372 double t1 = 0;
373 double t2 = 360;
374 double inc = 1.0;
375 for (double t = t1; t <= t2; t += inc)
376 {
377 //step along circumference
378 dms a(t);
379 double sa(0.), ca(0.);
380 a.SinCos(sa, ca);
381 m_clipPolygon << QPointF(0.5 * m_vp.width + r * ca, 0.5 * m_vp.height - r * sa);
382 }
383 }
384
clipPoly() const385 QPolygonF Projector::clipPoly() const
386 {
387 return m_clipPolygon;
388 }
389
unusablePoint(const QPointF & p) const390 bool Projector::unusablePoint(const QPointF &p) const
391 {
392 //r0 is the angular size of the sky horizon, in radians
393 double r0 = radius();
394 //If the zoom is high enough, all points are usable
395 //The center-to-corner distance, in radians
396 double r = 0.5 * 1.41421356 * m_vp.width / m_vp.zoomFactor;
397 if (r < r0)
398 return false;
399 //At low zoom, we have to determine whether the point is beyond the sky horizon
400 //Convert pixel position to x and y offsets in radians
401 double dx = (0.5 * m_vp.width - p.x()) / m_vp.zoomFactor;
402 double dy = (0.5 * m_vp.height - p.y()) / m_vp.zoomFactor;
403 return (dx * dx + dy * dy) > r0 * r0;
404 }
405
fromScreen(const QPointF & p,dms * LST,const dms * lat,bool onlyAltAz) const406 SkyPoint Projector::fromScreen(const QPointF &p, dms *LST, const dms *lat, bool onlyAltAz) const
407 {
408 dms c;
409 double sinc, cosc;
410 /** N.B. We don't cache these sin/cos values in the inverse
411 * projection because it causes 'shaking' when moving the sky.
412 */
413 double sinY0, cosY0;
414 //Convert pixel position to x and y offsets in radians
415 double dx = (0.5 * m_vp.width - p.x()) / m_vp.zoomFactor;
416 double dy = (0.5 * m_vp.height - p.y()) / m_vp.zoomFactor;
417
418 double r = sqrt(dx * dx + dy * dy);
419 c.setRadians(projectionL(r));
420 c.SinCos(sinc, cosc);
421
422 if (m_vp.useAltAz)
423 {
424 dx = -1.0 * dx; //Azimuth goes in opposite direction compared to RA
425 SkyPoint::refract(m_vp.focus->alt(), m_vp.useRefraction).SinCos(sinY0, cosY0);
426 }
427 else
428 {
429 m_vp.focus->dec().SinCos(sinY0, cosY0);
430 }
431
432 double Y = asin(cosc * sinY0 + (r == 0 ? 0 : (dy * sinc * cosY0) / r));
433 double atop = dx * sinc;
434 double abot = r * cosY0 * cosc - dy * sinY0 * sinc;
435 double A = atan2(atop, abot);
436
437 SkyPoint result;
438 if (m_vp.useAltAz)
439 {
440 dms alt, az;
441 alt.setRadians(Y);
442 az.setRadians(A + m_vp.focus->az().radians());
443 if (m_vp.useRefraction)
444 alt = SkyPoint::unrefract(alt);
445 result.setAlt(alt);
446 result.setAz(az);
447 if (onlyAltAz)
448 return result;
449 result.HorizontalToEquatorial(LST, lat);
450 }
451 else
452 {
453 dms ra, dec;
454 dec.setRadians(Y);
455 ra.setRadians(A + m_vp.focus->ra().radians());
456 result.set(ra.reduce(), dec);
457 result.EquatorialToHorizontal(LST, lat);
458 }
459
460 return result;
461 }
462
toScreenVec(const SkyPoint * o,bool oRefract,bool * onVisibleHemisphere) const463 Eigen::Vector2f Projector::toScreenVec(const SkyPoint *o, bool oRefract, bool *onVisibleHemisphere) const
464 {
465 double Y, dX;
466 double sindX, cosdX, sinY, cosY;
467
468 oRefract &= m_vp.useRefraction;
469 if (m_vp.useAltAz)
470 {
471 if (oRefract)
472 Y = SkyPoint::refract(o->alt()).radians(); //account for atmospheric refraction
473 else
474 Y = o->alt().radians();
475 dX = m_vp.focus->az().radians() - o->az().radians();
476 }
477 else
478 {
479 dX = o->ra().radians() - m_vp.focus->ra().radians();
480 Y = o->dec().radians();
481 }
482
483 if (!(std::isfinite(Y) && std::isfinite(dX)))
484 {
485 return Eigen::Vector2f(0, 0);
486
487 // JM: Enable this again later when trying to find a solution for it
488 // As it is now creating too much noise in the log file.
489 /*
490 qDebug() << "Assert in Projector::toScreenVec failed!";
491 qDebug() << "using AltAz?" << m_vp.useAltAz << " Refract? " << oRefract;
492 const SkyObject *obj;
493 qDebug() << "Point supplied has RA0 = " << o->ra0().toHMSString() << " Dec0 = " << o->dec0().toDMSString() << "; alt = " << o->alt().toDMSString() << "; az = " << o->az().toDMSString();
494 if ( (obj = dynamic_cast<const SkyObject *>(o) ) ) {
495 qDebug() << "Point is object with name = " << obj->name() << " longname = " << obj->longname();
496 }
497 qDebug() << "dX = " << dX << " and isfinite(dX) is" << std::isfinite(dX);
498 qDebug() << "Y = " << Y << " and isfinite(Y) is" << std::isfinite(Y);
499
500 //Q_ASSERT( false );
501 */
502 }
503
504 dX = KSUtils::reduceAngle(dX, -dms::PI, dms::PI);
505
506 //Convert dX, Y coords to screen pixel coords, using GNU extension if available
507 #if (__GLIBC__ >= 2 && __GLIBC_MINOR__ >= 1)
508 sincos(dX, &sindX, &cosdX);
509 sincos(Y, &sinY, &cosY);
510 #else
511 sindX = sin(dX);
512 cosdX = cos(dX);
513 sinY = sin(Y);
514 cosY = cos(Y);
515 #endif
516
517 //c is the cosine of the angular distance from the center
518 double c = m_sinY0 * sinY + m_cosY0 * cosY * cosdX;
519
520 //If c is less than 0.0, then the "field angle" (angular distance from the focus)
521 //is more than 90 degrees. This is on the "back side" of the celestial sphere
522 //and should not be drawn.
523 if (onVisibleHemisphere)
524 *onVisibleHemisphere =
525 (c >
526 cosMaxFieldAngle()); // TODO: Isn't it more efficient to bypass the final calculation below if the object is not visible?
527
528 double k = projectionK(c);
529
530 double origX = m_vp.width / 2;
531 double origY = m_vp.height / 2;
532
533 double x = origX - m_vp.zoomFactor * k * cosY * sindX;
534 double y = origY - m_vp.zoomFactor * k * (m_cosY0 * sinY - m_sinY0 * cosY * cosdX);
535 #ifdef KSTARS_LITE
536 double skyRotation = SkyMapLite::Instance()->getSkyRotation();
537 if (skyRotation != 0)
538 {
539 dms rotation(skyRotation);
540 double cosT, sinT;
541
542 rotation.SinCos(sinT, cosT);
543
544 double newX = origX + (x - origX) * cosT - (y - origY) * sinT;
545 double newY = origY + (x - origX) * sinT + (y - origY) * cosT;
546
547 x = newX;
548 y = newY;
549 }
550 #endif
551 return Eigen::Vector2f(x, y);
552 }
553