projector.cpp

/*
    SPDX-FileCopyrightText: 2010 Henry de Valence <hdevalence@gmail.com>
 
    SPDX-License-Identifier: GPL-2.0-or-later
*/
 
#include "projector.h"
 
#include "ksutils.h"
#ifdef KSTARS_LITE
#include "skymaplite.h"
#endif
#include "skycomponents/skylabeler.h"
 
namespace
{
void toXYZ(const SkyPoint *p, double *x, double *y, double *z)
{
    double sinRa, sinDec, cosRa, cosDec;
 
    p->ra().SinCos(sinRa, cosRa);
    p->dec().SinCos(sinDec, cosDec);
    *x = cosDec * cosRa;
    *y = cosDec * sinRa;
    *z = sinDec;
}
}
 
SkyPoint Projector::pointAt(double az)
{
    SkyPoint p;
    p.setAz(az);
    p.setAlt(0.0);
    p.HorizontalToEquatorial(KStarsData::Instance()->lst(), KStarsData::Instance()->geo()->lat());
    return p;
}
 
Projector::Projector(const ViewParams &p)
{
    m_data = KStarsData::Instance();
    setViewParams(p);
    // Force clip polygon update
    updateClipPoly();
}
 
void Projector::setViewParams(const ViewParams &p)
{
    m_vp = p;
 
    /** Precompute cached values */
    //Find Sin/Cos for focus point
    m_sinY0 = 0;
    m_cosY0 = 0;
    if (m_vp.useAltAz)
    {
        // N.B. We explicitly check useRefraction and not use
        // SkyPoint::altRefracted() here because it could be different
        // from Options::useRefraction() in some future use-case
        // --asimha
        SkyPoint::refract(m_vp.focus->alt(), m_vp.useRefraction).SinCos(m_sinY0, m_cosY0);
    }
    else
    {
        m_vp.focus->dec().SinCos(m_sinY0, m_cosY0);
    }
 
    double currentFOV = m_fov;
    //Find FOV in radians
    m_fov = sqrt(m_vp.width * m_vp.width + m_vp.height * m_vp.height) / (2 * m_vp.zoomFactor * dms::DegToRad);
    //Set checkVisibility variables
    double Ymax;
    m_xrange = 1.2 * m_fov / m_cosY0;
    if (m_vp.useAltAz)
    {
        Ymax     = fabs(SkyPoint::refract(m_vp.focus->alt().Degrees(), m_vp.useRefraction)) + m_fov;
    }
    else
    {
        Ymax     = fabs(m_vp.focus->dec().Degrees()) + m_fov;
    }
    m_isPoleVisible = (Ymax >= 90.0);
 
    // Only update clipping polygon if there is an FOV change
    if (currentFOV != m_fov)
        updateClipPoly();
}
 
double Projector::fov() const
{
    return m_fov;
}
 
QPointF Projector::toScreen(const SkyPoint *o, bool oRefract, bool *onVisibleHemisphere) const
{
    return KSUtils::vecToPoint(toScreenVec(o, oRefract, onVisibleHemisphere));
}
 
bool Projector::onScreen(const QPointF &p) const
{
    return (0 <= p.x() && p.x() <= m_vp.width && 0 <= p.y() && p.y() <= m_vp.height);
}
 
bool Projector::onScreen(const Eigen::Vector2f &p) const
{
    return onScreen(QPointF(p.x(), p.y()));
}
 
QPointF Projector::clipLine(SkyPoint *p1, SkyPoint *p2) const
{
    return KSUtils::vecToPoint(clipLineVec(p1, p2));
}
 
Eigen::Vector2f Projector::clipLineVec(SkyPoint *p1, SkyPoint *p2) const
{
    /* ASSUMES p1 was not clipped but p2 was.
     * Return the QPoint that barely clips in the line twixt p1 and p2.
     */
    //TODO: iteration = ceil( 0.5*log2( w^2 + h^2) )??
    //      also possibly rewrite this
    //     --hdevalence
    int iteration = 15; // For "perfect" clipping:
    // 2^iterations should be >= max pixels/line
    bool isVisible = true; // so we start at midpoint
    SkyPoint mid;
    Eigen::Vector2f oMid;
    double x, y, z, dx, dy, dz, ra, dec;
    int newx, newy, oldx, oldy;
    oldx = oldy = -10000; // any old value that is not the first omid
 
    toXYZ(p1, &x, &y, &z);
    // -jbb printf("\np1: %6.4f %6.4f %6.4f\n", x, y, z);
 
    toXYZ(p2, &dx, &dy, &dz);
 
    // -jbb printf("p2: %6.4f %6.4f %6.4f\n", dx, dy, dz);
    dx -= x;
    dy -= y;
    dz -= z;
    // Successive approximation to point on line that just clips.
    while (iteration-- > 0)
    {
        dx *= .5;
        dy *= .5;
        dz *= .5;
        if (!isVisible) // move back toward visible p1
        {
            x -= dx;
            y -= dy;
            z -= dz;
        }
        else // move out toward clipped p2
        {
            x += dx;
            y += dy;
            z += dz;
        }
 
        // -jbb printf("  : %6.4f %6.4f %6.4f\n", x, y, z);
        // [x, y, z] => [ra, dec]
        ra  = atan2(y, x);
        dec = asin(z / sqrt(x * x + y * y + z * z));
 
        mid = SkyPoint(ra * 12. / dms::PI, dec * 180. / dms::PI);
        mid.EquatorialToHorizontal(m_data->lst(), m_data->geo()->lat());
 
        oMid = toScreenVec(&mid, false, &isVisible);
        //AND the result with checkVisibility to clip things going below horizon
        isVisible &= checkVisibility(&mid);
        newx = (int)oMid.x();
        newy = (int)oMid.y();
 
        // -jbb printf("new x/y: %4d %4d", newx, newy);
        if ((oldx == newx) && (oldy == newy))
        {
            break;
        }
        oldx = newx;
        oldy = newy;
    }
    return oMid;
}
 
bool Projector::checkVisibility(const SkyPoint *p) const
{
    //TODO deal with alternate projections
    //not clear how this depends on projection
    //FIXME do these heuristics actually work?
 
    double dX, dY;
 
    //Skip objects below the horizon if the ground is drawn
    /* Is the cost of this conversion actually less than drawing it anyways?
     * EquatorialToHorizontal takes 3 SinCos calls -- so 6 trig calls if not using GNU exts.
     */
    /*
    if( m_vp.fillGround  ) {
        if( !m_vp.useAltAz )
            p->EquatorialToHorizontal( m_data->lst(), m_data->geo()->lat() );
        if( p->alt().Degrees() < -1.0 ) return false;
    }
    */ //Here we hope that the point has already been 'synchronized'
    if (m_vp.fillGround /*&& m_vp.useAltAz*/ && p->alt().Degrees() <= SkyPoint::altCrit)
        return false;
 
    if (m_vp.useAltAz)
    {
        /** To avoid calculating refraction, we just use the unrefracted
            altitude and add a 2-degree 'safety factor' */
        dY = fabs(p->alt().Degrees() - m_vp.focus->alt().Degrees()) - 2.;
    }
    else
    {
        dY = fabs(p->dec().Degrees() - m_vp.focus->dec().Degrees());
    }
    if (m_isPoleVisible)
        dY *= 0.75; //increase effective FOV when pole visible.
    if (dY > m_fov)
        return false;
    if (m_isPoleVisible)
        return true;
 
    if (m_vp.useAltAz)
    {
        dX = fabs(p->az().Degrees() - m_vp.focus->az().Degrees());
    }
    else
    {
        dX = fabs(p->ra().Degrees() - m_vp.focus->ra().Degrees());
    }
    if (dX > 180.0)
        dX = 360.0 - dX; // take shorter distance around sky
 
    return dX < m_xrange;
}
 
double Projector::findNorthPA(const SkyPoint *o, float x, float y) const
{
    //Find position angle of North using a test point displaced to the north
    //displace by 100/zoomFactor radians (so distance is always 100 pixels)
    //this is 5730/zoomFactor degrees
 
    // N.B. It is not sufficient to find the angle at the center of
    // the screen, and the angle of the NCP can vary drastically
    // across the screen (e.g. when pointed towards the NCP)
 
    KStarsData *data = KStarsData::Instance();
    double newDec    = o->dec().Degrees() + 5730.0 / m_vp.zoomFactor;
    if (newDec > 90.0)
        newDec = 90.0;
    SkyPoint test(o->ra().Hours(), newDec);
    if (m_vp.useAltAz)
        test.EquatorialToHorizontal(data->lst(), data->geo()->lat());
    Eigen::Vector2f t = toScreenVec(&test);
    float dx    = t.x() - x;
    float dy    = y - t.y(); //backwards because QWidget Y-axis increases to the bottom (FIXME: Check)
    // float dx = dx_ * m_vp.rotationAngle.cos() - dy_ * m_vp.rotationAngle.sin();
    // float dy = dx_ * m_vp.rotationAngle.sin() + dy_ * m_vp.rotationAngle.cos();
    float north;
    if (dy)
    {
        north = atan2f(dx, dy) * 180.0 / dms::PI;
    }
    else
    {
        north = (dx > 0.0 ? -90.0 : 90.0);
    }
 
    return north;
}
 
double Projector::findPA(const SkyObject *o, float x, float y) const
{
    return (findNorthPA(o, x, y) + (m_vp.mirror ? -o->pa() : o->pa()));
}
 
double Projector::findZenithPA(const SkyPoint *o, float x, float y) const
{
    //Find position angle of North using a test point displaced to the north
    //displace by 100/zoomFactor radians (so distance is always 100 pixels)
    //this is 5730/zoomFactor degrees
 
    // N.B. It is not sufficient to find the angle at the center of
    // the screen, and the angle of the zenith can vary drastically
    // across the screen (e.g. when pointed towards the zenith)
 
    KStarsData *data = KStarsData::Instance();
    double newAlt    = o->alt().Degrees() + 5730.0 / m_vp.zoomFactor;
    if (newAlt > 90.0)
        newAlt = 90.0;
    SkyPoint test;
    test.setAlt(newAlt);
    test.setAz(o->az().Degrees());
    if (!m_vp.useAltAz)
        test.HorizontalToEquatorial(data->lst(), data->geo()->lat());
    Eigen::Vector2f t = toScreenVec(&test);
    float dx    = t.x() - x;
    float dy    = y - t.y(); //backwards because QWidget Y-axis increases to the bottom (FIXME: Check)
    // float dx = dx_ * m_vp.rotationAngle.cos() - dy_ * m_vp.rotationAngle.sin();
    // float dy = dx_ * m_vp.rotationAngle.sin() + dy_ * m_vp.rotationAngle.cos();
    float zenith;
    if (dy)
    {
        zenith = atan2f(dx, dy) * 180.0 / dms::PI;
    }
    else
    {
        zenith = (dx > 0.0 ? -90.0 : 90.0);
    }
 
    return zenith;
}
 
QVector<Eigen::Vector2f> Projector::groundPoly(SkyPoint *labelpoint, bool *drawLabel) const
{
    QVector<Eigen::Vector2f> ground;
 
    static const QString horizonLabel = i18n("Horizon");
    float marginLeft, marginRight, marginTop, marginBot;
    SkyLabeler::Instance()->getMargins(horizonLabel, &marginLeft, &marginRight, &marginTop, &marginBot);
 
    //daz is 1/2 the width of the sky in degrees
    double daz = 90.;
    if (m_vp.useAltAz && m_vp.rotationAngle.reduce().Degrees() == 0.0)
    {
        daz = 0.5 * m_vp.width / (dms::DegToRad * m_vp.zoomFactor); //center to edge, in degrees
        if (type() == Projector::Orthographic)
        {
            daz = daz * 1.4;
        }
        daz = qMin(qreal(90.0), daz);
    }
 
    double faz = m_vp.focus->az().Degrees();
    double az1 = faz - daz;
    double az2 = faz + daz;
 
    bool allGround = true;
    bool allSky    = true;
 
    bool inverted = ((m_vp.rotationAngle + 90.0_deg).reduce().Degrees() > 180.);
 
    double inc = 1.0;
    //Add points along horizon
    for (double az = az1; az <= az2 + inc; az += inc)
    {
        SkyPoint p   = pointAt(az);
        bool visible = false;
        Eigen::Vector2f o   = toScreenVec(&p, false, &visible);
        if (visible)
        {
            ground.append(o);
            //Set the label point if this point is onscreen
            if (labelpoint && o.x() < marginRight && o.y() > marginTop && o.y() < marginBot)
                *labelpoint = p;
 
            if (o.y() > 0.)
                (inverted ? allSky : allGround) = false;
            if (o.y() < m_vp.height)
                (inverted ? allGround : allSky) = false;
        }
    }
 
    if (allSky)
    {
        if (drawLabel)
            *drawLabel = false;
        return QVector<Eigen::Vector2f>();
    }
 
    if (allGround)
    {
        ground.clear();
        ground.append(Eigen::Vector2f(-10., -10.));
        ground.append(Eigen::Vector2f(m_vp.width + 10., -10.));
        ground.append(Eigen::Vector2f(m_vp.width + 10., m_vp.height + 10.));
        ground.append(Eigen::Vector2f(-10., m_vp.height + 10.));
        if (drawLabel)
            *drawLabel = false;
        return ground;
    }
 
    if (m_vp.mirror)
        std::reverse(ground.begin(), ground.end());
 
    //In Gnomonic projection, or if sufficiently zoomed in, we can complete
    //the ground polygon by simply adding offscreen points
    //FIXME: not just gnomonic
    if (daz < 25.0 || type() == Projector::Gnomonic)
    {
        const float completion_height = (inverted ?
                                         -10.f : m_vp.height + 10.f);
        ground.append(Eigen::Vector2f(m_vp.width + 10.f, ground.last().y()));
        ground.append(Eigen::Vector2f(m_vp.width + 10.f, completion_height));
        ground.append(Eigen::Vector2f(-10.f, completion_height));
        ground.append(Eigen::Vector2f(-10.f, ground.first().y()));
    }
    else
    {
        double r = m_vp.zoomFactor * radius();
        double t1 =
            atan2(-1. * (ground.last().y() - 0.5 * m_vp.height), ground.last().x() - 0.5 * m_vp.width) / dms::DegToRad;
        double t2 = t1 - 180.;
        for (double t = t1; t >= t2; t -= inc) //step along circumference
        {
            dms a(t);
            double sa(0.), ca(0.);
            a.SinCos(sa, ca);
            ground.append(Eigen::Vector2f(0.5 * m_vp.width + r * ca, 0.5 * m_vp.height - r * sa));
        }
    }
 
    if (drawLabel)
        *drawLabel = true;
    return ground;
}
 
void Projector::updateClipPoly()
{
    m_clipPolygon.clear();
 
    double r   = m_vp.zoomFactor * radius();
    double t1  = 0;
    double t2  = 360;
    double inc = 1.0;
    for (double t = t1; t <= t2; t += inc)
    {
        //step along circumference
        dms a(t);
        double sa(0.), ca(0.);
        a.SinCos(sa, ca);
        m_clipPolygon << QPointF(0.5 * m_vp.width + r * ca, 0.5 * m_vp.height - r * sa);
    }
}
 
QPolygonF Projector::clipPoly() const
{
    return m_clipPolygon;
}
 
bool Projector::unusablePoint(const QPointF &p) const
{
    //r0 is the angular size of the sky horizon, in radians
    double r0 = radius();
    //If the zoom is high enough, all points are usable
    //The center-to-corner distance, in radians
    double r = 0.5 * 1.41421356 * m_vp.width / m_vp.zoomFactor;
    if (r < r0)
        return false;
    //At low zoom, we have to determine whether the point is beyond the sky horizon
    //Convert pixel position to x and y offsets in radians
 
    // N.B. Technically, rotation does not affect the dx² + dy²
    // metric, but we use the derst method for uniformity; this
    // function is not perf critical
    auto p_ = derst(p.x(), p.y());
    double dx = p_[0], dy = p_[1];
    return (dx * dx + dy * dy) > r0 * r0;
}
 
SkyPoint Projector::fromScreen(const QPointF &p, KStarsData* data, bool onlyAltAz) const
{
    dms c;
    double sinc, cosc;
    /** N.B. We don't cache these sin/cos values in the inverse
      * projection because it causes 'shaking' when moving the sky.
      */
    double sinY0, cosY0;
    //Convert pixel position to x and y offsets in radians
    auto p_ = derst(p.x(), p.y());
    double dx = p_[0], dy = p_[1];
 
    double r = sqrt(dx * dx + dy * dy);
    c.setRadians(projectionL(r));
    c.SinCos(sinc, cosc);
 
    if (m_vp.useAltAz)
    {
        dx = -1.0 * dx; //Azimuth goes in opposite direction compared to RA
        SkyPoint::refract(m_vp.focus->alt(), m_vp.useRefraction).SinCos(sinY0, cosY0);
    }
    else
    {
        m_vp.focus->dec().SinCos(sinY0, cosY0);
    }
 
    double Y    = asin(cosc * sinY0 + (r == 0 ? 0 : (dy * sinc * cosY0) / r));
    double atop = dx * sinc;
    double abot = r * cosY0 * cosc - dy * sinY0 * sinc;
    double A    = atan2(atop, abot);
 
    const auto LST = data->lst();
    const auto lat = data->geo()->lat();
    const auto jdf = data->djd();
 
    SkyPoint result;
    if (m_vp.useAltAz)
    {
        dms alt, az;
        alt.setRadians(Y);
        az.setRadians(A + m_vp.focus->az().radians());
        if (m_vp.useRefraction)
            alt = SkyPoint::unrefract(alt);
        result.setAlt(alt);
        result.setAz(az);
        if (onlyAltAz)
            return result;
        result.HorizontalToEquatorialICRS(LST, lat, jdf);
    }
    else
    {
        dms ra, dec;
        dec.setRadians(Y);
        ra.setRadians(A + m_vp.focus->ra().radians());
        SkyPoint p(ra, dec);
        p.catalogueCoord(jdf);
        result.set(ra.reduce(), dec);
        result.setRA0(p.ra0());
        result.setDec0(p.dec0());
        result.EquatorialToHorizontal(LST, lat);
    }
 
    return result;
}
 
Eigen::Vector2f Projector::toScreenVec(const SkyPoint *o, bool oRefract, bool *onVisibleHemisphere) const
{
    double Y, dX;
    double sindX, cosdX, sinY, cosY;
 
    oRefract &= m_vp.useRefraction;
    if (m_vp.useAltAz)
    {
        if (oRefract)
            Y = SkyPoint::refract(o->alt()).radians(); //account for atmospheric refraction
        else
            Y = o->alt().radians();
        dX = m_vp.focus->az().radians() - o->az().radians();
    }
    else
    {
        dX = o->ra().radians() - m_vp.focus->ra().radians();
        Y  = o->dec().radians();
    }
 
    if (!(std::isfinite(Y) && std::isfinite(dX)))
    {
        return Eigen::Vector2f(0, 0);
 
        // JM: Enable this again later when trying to find a solution for it
        //     As it is now creating too much noise in the log file.
        /*
        qDebug() << Q_FUNC_INFO << "Assert in Projector::toScreenVec failed!";
        qDebug() << Q_FUNC_INFO << "using AltAz?" << m_vp.useAltAz << " Refract? " << oRefract;
        const SkyObject *obj;
        qDebug() << Q_FUNC_INFO << "Point supplied has RA0 = " << o->ra0().toHMSString() << " Dec0 = " << o->dec0().toDMSString() << "; alt = " << o->alt().toDMSString() << "; az = " << o->az().toDMSString();
        if ( (obj = dynamic_cast<const SkyObject *>(o) ) ) {
            qDebug() << Q_FUNC_INFO << "Point is object with name = " << obj->name() << " longname = " << obj->longname();
        }
        qDebug() << Q_FUNC_INFO << "dX = " << dX << " and isfinite(dX) is" << std::isfinite(dX);
        qDebug() << Q_FUNC_INFO << "Y = " << Y << " and isfinite(Y) is" << std::isfinite(Y);
 
        //Q_ASSERT( false );
        */
    }
 
    dX = KSUtils::reduceAngle(dX, -dms::PI, dms::PI);
 
    //Convert dX, Y coords to screen pixel coords, using GNU extension if available
#ifdef HAVE_SINCOS
    sincos(dX, &sindX, &cosdX);
    sincos(Y, &sinY, &cosY);
#else
    sindX = sin(dX);
    cosdX = cos(dX);
    sinY  = sin(Y);
    cosY  = cos(Y);
#endif
 
    //c is the cosine of the angular distance from the center
    double c = m_sinY0 * sinY + m_cosY0 * cosY * cosdX;
 
    //If c is less than 0.0, then the "field angle" (angular distance from the focus)
    //is more than 90 degrees.  This is on the "back side" of the celestial sphere
    //and should not be drawn.
    if (onVisibleHemisphere)
        *onVisibleHemisphere =
            (c >
             cosMaxFieldAngle()); // TODO: Isn't it more efficient to bypass the final calculation below if the object is not visible?
 
    double k = projectionK(c);
 
    auto p = rst(k * cosY * sindX, k * (m_cosY0 * sinY - m_sinY0 * cosY * cosdX));
 
#ifdef KSTARS_LITE
    double skyRotation = SkyMapLite::Instance()->getSkyRotation();
    // FIXME: Port above to change the m_vp.rotationAngle, or
    // deprecate it
    Q_ASSERT(false);
#endif
    return p;
}