#include "imageinfodata.h"
#include <QTime>
#include <QRectF>
#include <QRegularExpression>
#include <wcslib/wcshdr.h>
#include <wcslib/wcsfix.h>
#include "database.h"
#include "libxisf.h"

static const QVector<QByteArray> noEditableKey = {"SIMPLE", "BITPIX", "NAXIS", "NAXIS1", "NAXIS2", "NAXIS3", "EXTEND", "BZERO", "BSCALE"};

bool FITSRecord::editable() const
{
    return noEditableKey.count(key);
}

FITSRecord::FITSRecord(const QByteArray &key, const QVariant &value, const QByteArray &comment) :
    key(key), value(value), comment(comment)
{
}

FITSRecord::FITSRecord(const LibXISF::FITSKeyword &record)
{
    key = record.name.c_str();
    comment = record.comment.c_str();

    QString string = record.value.c_str();
    if(string.startsWith('\'') && string.endsWith('\''))
    {
        string.chop(1);
        string.remove(0, 1);
    }
    bool isint;
    bool isdouble;
    double vald = string.toDouble(&isdouble);
    long long vall = string.toLongLong(&isint);
    if(isint)
        value = vall;
    else if(isdouble)
        value = vald;
    else if(string == "T" || string == "F")
        value = string == "T";
    else
        value = string;
}

FITSRecord::FITSRecord(const LibXISF::Property &property)
{
    key = property.id.c_str();
    value = QString::fromStdString(property.value.toString());
    comment = property.comment.c_str();
    xisf = true;
}

void WCSDataT::freeWCS()
{
    wcsvfree(&nwcs, &wcs);
    nwcs = 0;
    wcs = nullptr;
}

WCSDataT::WCSDataT(int width, int height, char *header, int nrec) :
    width(width),
    height(height)
{
    int nreject = 0;
    int status = wcspih(header, nrec, 1, 0, &nreject, &nwcs, &wcs);
    if(status != 0)
    {
        freeWCS();
        return;
    }
    status = cdfix(wcs);
    if(status > 0 || wcs->crpix[0] == 0)
        freeWCS();
}

WCSDataT::WCSDataT(int width, int height, const QVector<FITSRecord> &header) :
    width(width),
    height(height)
{
    int status = 0;

    QByteArray str;
    int nrec = 1;
    for(const FITSRecord &record : header)
    {
        if(record.key.startsWith("PV"))continue;
        if(record.xisf)continue;

        QByteArray rec;
        rec.append(record.key.leftJustified(8, ' '));
        rec.append("= ");
        rec.append(record.value.toString().toLatin1());
        rec.append(" / ");
        rec.append(record.comment);
        str.append(rec.leftJustified(80, ' ', true));
        nrec++;
    }
    str.append(QByteArray("END").leftJustified(80));

    int nreject = 0;
    status = wcspih(str.data(), nrec, 1, 0, &nreject, &nwcs, &wcs);
    if(status != 0)
    {
        freeWCS();
        return;
    }
    status = cdfix(wcs);
    if(status > 0 || wcs->crpix[0] == 0)
        freeWCS();
}

WCSDataT::~WCSDataT()
{
    if(wcs)
        freeWCS();
}

bool WCSDataT::pixelToWorld(const QPointF &pixel, SkyPoint &point) const
{
    if(!valid())return false;

    double pixcrd[2] = {pixel.x(), pixel.y()};
    double imgcrd[8] = {0};
    double phi = 0;
    double theta = 0;
    double world[8] = {0};
    int stat[NWCSFIX] = {0};
    int status = wcsp2s(wcs, 1, 2, pixcrd, imgcrd, &phi, &theta, world, stat);
    if(status == 0)
    {
        point = SkyPoint(world[0], world[1]);
        return true;
    }
    return false;
}

bool WCSDataT::worldToPixel(const SkyPoint &point, QPointF &pixel) const
{
    if(!valid())return false;

    double world[2] = {point.RA(), point.DEC()};
    double phi = 0;
    double theta = 0;
    double imgcrd[8] = {0};
    double pixcrd[8] = {0};
    int stat[NWCSFIX] = {0};
    int status = wcss2p(wcs, 1, 2, world, &phi, &theta, imgcrd, pixcrd, stat);
    if(status == 0)
    {
        pixel = QPointF(pixcrd[0], pixcrd[1]);
        return true;
    }
    return false;
}

bool WCSDataT::calculateBounds(double &minRa, double &maxRa, double &minDec, double &maxDec, double &crVal1, double &crVal2) const
{
    if(wcs == nullptr)return false;

    minRa = 1000;
    maxRa = -1000;
    minDec = 1000;
    maxDec = -1000;

    if(wcs->crval)
    {
        crVal1 = wcs->crval[0];
        crVal2 = wcs->crval[1];
    }
    else
    {
        crVal1 = crVal2 = NAN;
    }

    auto update = [&](const QPointF &pixel)
    {
        SkyPoint point;
        pixelToWorld(pixel, point);
        minRa = std::min(minRa, point.RA());
        maxRa = std::max(maxRa, point.RA());
        minDec = std::min(minDec, point.DEC());
        maxDec = std::max(maxDec, point.DEC());
    };

    for(int x=0; x<width; x++)
    {
        update(QPointF(x, 0));
        update(QPointF(x, height - 1));
    }

    for(int y=0; y<height; y++)
    {
        update(QPointF(0, y));
        update(QPointF(width - 1, y));
    }

    QPointF ncp;
    QPointF scp;
    QRectF s(0, 0, width - 1, height - 1);
    if(worldToPixel(SkyPoint(0, 90), ncp))
    {
        if(s.contains(ncp))
            maxDec = 90;
    }

    if(worldToPixel(SkyPoint(0, -90), scp))
    {
        if(s.contains(scp))
            minDec = -90;
    }
    return true;
}

double hav(double x)
{
    return (1.0 - std::cos(x)) * 0.5;
}

double haverSine(const SkyPoint &a, SkyPoint &b)
{
    const double ToRAD = M_PI / 180.0;
    double d = hav((a.DEC() - b.DEC()) * ToRAD) + std::cos(a.DEC() * ToRAD) * std::cos(b.DEC() * ToRAD) * hav((a.RA() - b.RA()) * ToRAD);
    return std::acos(1.0 - 2.0 * d) * (180.0 / M_PI);
}

SkyPointScale WCSDataT::getRaDecScale() const
{
    SkyPointScale ret;
    pixelToWorld(QPointF(width/2.0, height/2.0), ret.point);
    SkyPoint pointX;
    SkyPoint pointY;
    pixelToWorld(QPointF(width/2.0+1, height/2.0), pointX);
    pixelToWorld(QPointF(width/2.0, height/2.0+1), pointY);
    double scaleX = haverSine(ret.point, pointX) * 3600.0;
    double scaleY = haverSine(ret.point, pointY) * 3600.0;
    ret.scaleLow = std::min(scaleX, scaleY);
    ret.scaleHigh = std::max(scaleX, scaleY);
    ret.scaleValid = true;
    return ret;
}

SkyPoint::SkyPoint() : ra(NAN), dec(NAN)
{
}

SkyPoint::SkyPoint(double ra, double dec) : ra(ra), dec(dec)
{
}

void SkyPoint::set(double ra, double dec)
{
    this->ra = ra;
    this->dec = dec;
}

QString SkyPoint::toString() const
{
    if(std::isnan(ra) || std::isnan(dec))
        return QString();

    QTime t(0, 0);
    t = t.addSecs(ra * 240);

    double deg, min, sec;
    min = std::abs(std::modf(dec, &deg) * 60);
    sec = std::modf(min, &min) * 60;
    return QString("RA: %1 DEC: %2° %3' %4\"").arg(t.toString("HH'h' mm'm' ss's'")).arg(deg, 2, 'f', 0, '0').arg(min, 2, 'f', 0, '0').arg(sec, 2, 'f', 0, '0');
}

QString SkyPoint::RAString() const
{
    return toHMS(ra / 15);
}

QString SkyPoint::DECString() const
{
    return toDMS(dec);
}

double SkyPoint::fromHMS(const QString &hms)
{
    double deg = fromDMS(hms);
    if(std::isnan(deg))return deg;
    return deg * 15.0;
}

double SkyPoint::fromDMS(const QString &dms)
{
    double deg = 0.0;
    QString str = dms.trimmed();
    str.remove(QRegularExpression("[hdms°'\"]"));
    str.replace(':', ' ');
    str.replace(QRegularExpression("\\s+"), " ");
    QStringList fields = str.split(' ');
    double sign = 1.0;

    bool ok = false;
    if(fields.size() >= 1)
        deg = fields.at(0).toDouble(&ok);
    if(!ok)return NAN;
    if(deg < 0.0)
        sign = -1.0;
    if(fields.size() >= 2)
        deg += sign * fields.at(1).toDouble() / 60.0;
    if(fields.size() >= 3)
        deg += sign * fields.at(2).toDouble() / 3600.0;

    return deg;
}

QString SkyPoint::toHMS(double decHour)
{
    double h,m,s,md;
    md = std::modf(decHour, &h) * 60.0;
    s = std::modf(md, &m) * 60.0;

    return QString("%1h %2m %3s").arg((int)h, 2, 10, QChar('0')).arg((int)m, 2, 10, QChar('0')).arg((int)s, 2, 10, QChar('0'));
}

QString SkyPoint::toDMS(double deg)
{
    int sign = deg < 0.0 ? -1 : 1;
    deg *= sign;
    double d,m,s,md;
    md = std::modf(deg, &d) * 60.0;
    s = std::modf(md, &m) * 60.0;

    return QString("%1˚ %2' %3\"").arg((int)d * sign, 2, 10, QChar('0')).arg((int)m, 2, 10, QChar('0')).arg((int)s, 2, 10, QChar('0'));
}

SkyPoint SkyPoint::operator+(const SkyPoint &p)
{
    SkyPoint ret;
    ret.ra = ra + p.ra;
    ret.dec = dec + p.dec;
    return ret;
}

SkyPointScale ImageInfoData::getCenterRaDec() const
{
    SkyPointScale ret;
    if(wcs && wcs->valid())
    {
        ret = wcs->getRaDecScale();
    }
    else
    {
        double ra,dec,focalLen,scale,pixSizeX,pixSizeY;
        int binX = 1;
        int binY = 1;
        ra = dec = focalLen = scale = pixSizeX = pixSizeY = NAN;
        bool ok;
        for(const FITSRecord &header : fitsHeader)
        {
            if(header.key == "OBJCTRA")
            {
                double tmp = SkyPoint::fromHMS(header.value.toString());
                if(!std::isnan(tmp))ra = tmp;
            }
            else if(header.key == "RA" && std::isnan(ra))
            {
                double tmp = header.value.toDouble(&ok);
                if(ok)ra = tmp;
            }
            else if(header.key == "OBJCTDEC")
            {
                double tmp = SkyPoint::fromDMS(header.value.toString());
                if(!std::isnan(tmp))dec = tmp;
            }
            else if(header.key == "DEC" && std::isnan(dec))
            {
                double tmp = SkyPoint::fromDMS(header.value.toString());
                if(!std::isnan(tmp))dec = tmp;
            }
            else if(header.key == "SCALE")
            {
                double tmp = header.value.toDouble(&ok);
                if(ok)scale = tmp;
            }
            else if(header.key == "FOCALLEN")
            {
                double tmp = header.value.toDouble(&ok);
                if(ok)focalLen = tmp;
            }
            else if(header.key == "PIXSIZE1" || header.key == "XPIXSZ")
            {
                pixSizeX = header.value.toDouble();
            }
            else if(header.key == "PIXSIZE2" || header.key == "YPIXSZ")
            {
                pixSizeY = header.value.toDouble();
            }
            else if(header.key == "XBINNING")
            {
                int tmp = header.value.toInt(&ok);
                if(ok)binX = tmp;
            }
            else if(header.key == "YBINNING")
            {
                int tmp = header.value.toInt(&ok);
                if(ok)binY = tmp;
            }
        }

        ret.point.set(ra, dec);
        if(!std::isnan(scale))
        {
            ret.scaleLow = ret.scaleHigh = scale;
            ret.scaleValid = true;
        }
        else if(!(std::isnan(focalLen) || std::isnan(pixSizeX) || std::isnan(pixSizeY)))
        {
            const double r = 206.2648097656; // (180 * 3600) / (1000 * pi) magic number to convert pixel size to focal length ratio to arcsec.
            ret.scaleLow = std::min(pixSizeX * binX / focalLen * r, pixSizeY * binY / focalLen * r);
            ret.scaleHigh = std::max(pixSizeX * binX / focalLen * r, pixSizeY * binY / focalLen * r);
            ret.scaleValid = true;
        }
    }

    if(ret.scaleValid)
    {
        ret.scaleLow *= 0.8;
        ret.scaleHigh *= 1.2;
    }
    return ret;
}

SkyPoint greatCircle(SkyPoint &p, double dist, double azm)
{
    dist = dist * M_PI / 180;
    azm = azm * M_PI / 180;
    double dec0 = p.DEC() * M_PI / 180;
    double ra0 = p.RA() * M_PI / 180;
    double dec1 = std::asin(std::sin(dec0) * std::cos(dist) + std::cos(dec0) * std::sin(dist) * std::cos(azm));
    double ra1 = ra0 + std::atan2(std::sin(azm) * std::sin(dist) * std::cos(dec0), std::cos(dist) - std::sin(dec0) * std::sin(dec1));
    return SkyPoint(ra1 * 180 / M_PI, dec1 * 180 / M_PI);
}

SkyGrid WCSDataT::prepareGrid(uint32_t w, uint32_t h, Database *database)
{
    SkyGrid skyGrid;
    if(!wcs)return skyGrid;

    double minRa, maxRa, minDec, maxDec, crVal1, crVal2;
    calculateBounds(minRa, maxRa, minDec, maxDec, crVal1, crVal2);

    QPointF a,b;
    worldToPixel(SkyPoint(crVal1, crVal2), a);
    worldToPixel(SkyPoint(crVal1 + 0.01, crVal2), b);
    skyGrid.rot_ang = std::atan2(b.y() - a.y(), b.x() - a.x()) / M_PI * -180.0;

    if(database)
    {
        skyGrid.objects = database->getObjects(minRa, maxRa, minDec, maxDec);
        for(auto &object : skyGrid.objects)
        {
            QPointF p;
            if(worldToPixel(object.skyPoint, p))
                object.pixel = p;

            QPointF majax;
            worldToPixel(greatCircle(object.skyPoint, (object.min_ax + object.maj_ax) / 120.0, object.pos_ang), majax);
            majax -= p;
            object.maj_ax = std::sqrt(QPointF::dotProduct(majax, majax));
        }
    }

    double raStep = 15;
    double decStep = 15;
    double raRange = maxRa - minRa;
    double decRange = maxDec - minDec;
    const QVector<double> raSteps = {15, 5, 2.5, 1.25, 0.25, 20/240.0, 10/240.0, 5/240.0, 1/240.0};
    const QVector<double> decSteps = {20, 10, 5, 2, 1, 20/60.0, 10/60.0, 5/60.0, 2/60.0, 1/60.0, 20/3600.0, 10/3600.0, 5/3600.0, 2/3600.0, 1/3600.0};

    for(double ra : raSteps)
    {
        if(ra * 5 <= raRange)
        {
            raStep = ra;
            break;
        }
    }

    for(double dec : decSteps)
    {
        if(dec * 5 <= decRange)
        {
            decStep = dec;
            break;
        }
    }

    minRa -= std::fmod(minRa, raStep);
    minDec -= std::fmod(minDec, decStep);
    if(minRa < 0)minRa -= raStep;
    if(minDec < 0)minDec -= decStep;

    QRectF clip(0, 0, w, h);
    const double step = 0.2;

    maxRa += raStep;
    maxDec += decStep;
    for(double ra = minRa; ra <= maxRa; ra += raStep)
    {
        QPointF p;
        worldToPixel(SkyPoint(ra, minDec), p);
        skyGrid.grid.moveTo(p);
        for(double dec = minDec + decStep * step; dec <= maxDec; dec += decStep * step)
        {
            worldToPixel(SkyPoint(ra, dec), p);
            skyGrid.grid.lineTo(p);
        }
    }

    for(double dec = minDec; dec <= maxDec; dec += decStep)
    {
        QPointF p;
        worldToPixel(SkyPoint(minRa, dec), p);
        skyGrid.grid.moveTo(p);
        for(double ra = minRa + raStep * step; ra <= maxRa; ra += raStep * step)
        {
            worldToPixel(SkyPoint(ra, dec), p);
            skyGrid.grid.lineTo(p);
        }
    }

    SkyPoint sp1, sp2,orig;
    pixelToWorld(QPointF(-1, -1), orig);
    sp1 = orig;
    for(uint32_t x = 0; x < w; x++)
    {
        QPointF p(x, 0);
        if(!pixelToWorld(p, sp2))
            break;

        if(static_cast<int>(sp1.RA() / raStep) != static_cast<int>(sp2.RA() / raStep))
            skyGrid.text.append({p, std::abs(sp1.RA()) > std::abs(sp2.RA()) ? sp1.RAString() : sp2.RAString()});

        if(static_cast<int>(sp1.DEC() / decStep) != static_cast<int>(sp2.DEC() / decStep))
            skyGrid.text.append({p, std::abs(sp1.DEC()) > std::abs(sp2.DEC()) ? sp1.DECString() : sp2.DECString()});

        sp1 = sp2;
    }

    sp1 = orig;
    for(uint32_t y = 0; y < h; y++)
    {
        QPointF p(0, y);
        if(!pixelToWorld(p, sp2))
            break;

        if(static_cast<int>(sp1.RA() / raStep) != static_cast<int>(sp2.RA() / raStep))
            skyGrid.text.append({p, std::abs(sp1.RA()) > std::abs(sp2.RA()) ? sp1.RAString() : sp2.RAString()});

        if(static_cast<int>(sp1.DEC() / decStep) != static_cast<int>(sp2.DEC() / decStep))
            skyGrid.text.append({p, std::abs(sp1.DEC()) > std::abs(sp2.DEC()) ? sp1.DECString() : sp2.DECString()});

        sp1 = sp2;
    }

    skyGrid.empty = false;
    return skyGrid;
}

void SkyGrid::clear()
{
    empty = true;
    grid.clear();
    text.clear();
    objects.clear();
}