tenmon/imageinfo.cpp

#include "imageinfo.h"
#include <QSettings>
#include <QTime>
#include <QHeaderView>
#include <wcslib/wcshdr.h>
#include <wcslib/wcsfix.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;
}

QByteArray FITSRecord::valueToByteArray() const
{
    if(value.type() == QVariant::Bool)
        return value.toBool() ? "T" : "F";
    else
        return value.toString().toLatin1();
}

ImageInfo::ImageInfo(QWidget *parent) : QTreeWidget(parent)
{
    setColumnCount(3);
    setHeaderLabels({tr("Property"), tr("Value"), tr("Comment")});
    setIndentation(5);
    QSettings settings;
    header()->restoreState(settings.value("imageinfo/headerstate").toByteArray());
}

ImageInfo::~ImageInfo()
{
    QSettings settings;
    settings.setValue("imageinfo/headerstate", header()->saveState());
}

void ImageInfo::setInfo(const ImageInfoData &info)
{
    clear();
    if(info.fitsHeader.size())
    {
        QTreeWidgetItem *fitsHeader = new QTreeWidgetItem({tr("FITS Header")});
        for(const FITSRecord &record : info.fitsHeader)
        {
            new QTreeWidgetItem(fitsHeader, {record.key, record.value.toString().left(1024), record.comment});
        }
        addTopLevelItem(fitsHeader);
    }
    if(info.info.size())
    {
        QTreeWidgetItem *infoHeader = new QTreeWidgetItem({tr("Image info")});
        for(auto &item : info.info)
        {
            new QTreeWidgetItem(infoHeader, {item.first, item.second});
        }
        addTopLevelItem(infoHeader);
    }
    expandAll();
}

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

WCSData::WCSData(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();
}

WCSData::WCSData(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;

        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();
}

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

bool WCSData::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 WCSData::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;
}

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

    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;
    }
}

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');
}