tenmon/solver.cpp

#include "solver.h"
#include <QJsonObject>
#include <QJsonDocument>
#include <fitsio.h>
#include <QStandardPaths>
#include <wcslib/wcshdr.h>
#include <wcslib/wcsutil.h>
#include "rawimage.h"
#include "loadrunable.h"
#include "scriptengine.h"

Solver::Solver(QObject *parent) : QObject(parent)
{
    _solver = new StellarSolver(this);

    QStringList indexFolder = getIndexPaths();

    _solver->setProperty("ProcessType", SSolver::SOLVE);
    _solver->setIndexFolderPaths(indexFolder);
    _solver->setParameterProfile(SSolver::Parameters::ALL_STARS);

    connect(_solver, &StellarSolver::finished, this, &Solver::finished);
}

Solver::~Solver()
{
}

void Solver::setIndexFolders(const QStringList &indexPaths)
{
    _solver->setIndexFolderPaths(indexPaths);
}

bool Solver::loadImage(const QString &path)
{
    _loaded = false;
    std::shared_ptr<RawImage> image;
    ImageInfoData info;
    if(::loadImage(path, info, image, true))
    {
        loadImage(image, path);
    }
    return false;
}

bool Solver::loadImage(std::shared_ptr<RawImage> &image, const QString &path)
{
    _rawImage = image;
    switch(_rawImage->type())
    {
    case RawImage::UINT8:
        _stats.dataType = TBYTE;
        break;
    case RawImage::UINT16:
        _stats.dataType = TUSHORT;
        break;
    case RawImage::UINT32:
        _stats.dataType = TUINT;
        break;
    case RawImage::FLOAT32:
        _stats.dataType = TFLOAT;
        break;
    case RawImage::FLOAT64:
        _stats.dataType = TDOUBLE;
        break;
    default:
        _error = tr("Unsupported image data type");
        return false;
        break;
    }
    _stats.bytesPerPixel = _rawImage->typeSize(_rawImage->type());
    _stats.channels = _rawImage->channels();
    _stats.width = _rawImage->width();
    _stats.height = _rawImage->height();
    _stats.samples_per_channel = _stats.width * _stats.height;

    _loaded = _solver->loadNewImageBuffer(_stats, (const uint8_t*)_rawImage->data());
    _path = path;
    return _loaded;
}

bool Solver::solveImage()
{
    if(_loaded && !_solver->isRunning())
    {
        _process = SSolver::ProcessType::SOLVE;
        _solver->setProperty("ProcessType", _process);
        _solver->start();
        return true;
    }
    return false;
}

bool Solver::extractSources(bool hfr)
{
    if(_loaded && !_solver->isRunning())
    {
        _process = hfr ? SSolver::ProcessType::EXTRACT_WITH_HFR : SSolver::ProcessType::EXTRACT;
        _solver->setProperty("ProcessType", _process);
        _solver->start();
        return true;
    }
    return false;
}

const FITSImage::Solution& Solver::getSolution() const
{
    return _solver->getSolution();
}

const QList<FITSImage::Star>& Solver::getStars() const
{
    return _solver->getStarList();
}

double Solver::getHFR() const
{
    double hfr = 0.0;
    auto stars = getStars();
    if(stars.empty())return -1.0;

    for(auto &star : stars)
    {
        hfr += star.HFR;
    }
    return hfr / stars.size();
}

QString Solver::errorMessage() const
{
    return _error;
}

void Solver::updateHeader()
{
    if(!_solver->solvingDone())return;

    FITSImage::Solution solution = getSolution();
    qDebug() << "RA" << solution.ra << "DEC" << solution.dec << "Orient" << solution.orientation << "field wxh" << solution.fieldWidth << solution.fieldHeight << solution.pixscale;
    qDebug() << "error" << solution.raError << solution.decError;

    double rotationDeg = 360.0 - solution.orientation;
    if(rotationDeg > 360)rotationDeg -= 360;
    double rotationRad = rotationDeg / 180.0 * M_PI;

    double cdeltx = (solution.parity == FITSImage::NEGATIVE ? solution.pixscale : -solution.pixscale) / 3600.0;
    double cdelty = solution.pixscale / 3600.0;

    Script::File file(_path, nullptr);
    Script::FITSRecordModify modify;
    modify.updateKeyword("CRPIX1", _stats.width / 2.0, QByteArray("x pixel coordinate of the reference point"));
    modify.updateKeyword("CRPIX2", _stats.height / 2.0, QByteArray("y pixel coordinate of the reference point"));
    modify.updateKeyword("CDELT1", cdeltx, QByteArray("X pixel size (deg)"));
    modify.updateKeyword("CDELT2", cdelty, QByteArray("Y pixel size (deg)"));
    modify.updateKeyword("CRVAL1", solution.ra, QByteArray("RA of reference pixel (deg)"));
    modify.updateKeyword("CRVAL2", solution.dec, QByteArray("DEC of reference pixel (deg)"));
    modify.updateKeyword("CD1_1", std::cos(rotationRad) * cdeltx, QByteArray("CD matrix to convert (x,y) to (RA, DEC)"));
    modify.updateKeyword("CD1_2",-std::sin(rotationRad) * cdelty, QByteArray("CD matrix to convert (x,y) to (RA, DEC)"));
    modify.updateKeyword("CD2_1", std::sin(rotationRad) * cdeltx, QByteArray("CD matrix to convert (x,y) to (RA, DEC)"));
    modify.updateKeyword("CD2_2", std::cos(rotationRad) * cdelty, QByteArray("CD matrix to convert (x,y) to (RA, DEC)"));
    modify.updateKeyword("CROTA1", rotationDeg, QByteArray("Image twist X axis (deg)"));
    modify.updateKeyword("CROTA2", rotationDeg, QByteArray("Image twist Y axis (deg)"));
    modify.updateKeyword("CTYPE1", "RA---TAN", QByteArray("first parameter RA, projection TANgential"));
    modify.updateKeyword("CTYPE2", "DEC--TAN", QByteArray("first parameter DEC, projection TANgential"));
    modify.updateKeyword("RADESYS", "ICRS", QByteArray("International Celestial Reference System"));
    modify.updateKeyword("EQUINOX", 2000, QByteArray("Equinox of coordinates"));
    file.modifyFITSRecords(&modify);
}

void Solver::setParameters(Parameters::ParametersProfile profile)
{
    _solver->setParameterProfile(profile);
}

void Solver::setParameters(const Parameters &parameters)
{
    _solver->setParameters(parameters);
}

QStringList Solver::getIndexPaths()
{
    QStringList paths;// = StellarSolver::getDefaultIndexFolderPaths();
    paths.prepend(getTenmonIndexPath());
    return paths;
}

QString Solver::getTenmonIndexPath()
{
    return QStandardPaths::writableLocation(QStandardPaths::AppDataLocation) + "/astrometry";
}

void Solver::finished()
{
    switch(_process)
    {
    case SSolver::ProcessType::SOLVE:
        emit solvingDone();
        break;
    case SSolver::ProcessType::EXTRACT_WITH_HFR:
    case SSolver::ProcessType::EXTRACT:
        emit extractionDone();
        break;
    }

}