tenmon/solver.cpp

#include "solver.h"
#include <QJsonObject>
#include <QJsonDocument>
#include <fitsio.h>
#include <QStandardPaths>
#include <QSettings>
#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);
    connect(_solver, &StellarSolver::logOutput, this, &Solver::logOutput);

    _solver->setProperty("ProcessType", SSolver::SOLVE);
    QSettings settings;
    setIndexFolder(settings.value("platesolving/indexPath", Solver::getTenmonIndexPath()).toString());
    int profileIdx = settings.value("platesolving/profile", 0).toInt();
    auto profiles = _solver->getBuiltInProfiles();
    _solver->setParameters(profiles[profileIdx]);

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

Solver::~Solver()
{
}

void Solver::setIndexFolder(const QString &indexPath)
{
    _solver->setIndexFolderPaths(QStringList(indexPath));
}

bool Solver::loadImage(const QString &path)
{
    if(path == _path)return true;

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

bool Solver::loadImage(std::shared_ptr<RawImage> &image, const QString &path)
{
    _rawImage = image;
    if(_rawImage->channels() > 1)
        _rawImagePlanar = _rawImage->toPlanar();
    else
        _rawImagePlanar = _rawImage;

    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(_rawImagePlanar->type());
    _stats.channels = _rawImagePlanar->channels();
    _stats.width = _rawImagePlanar->width();
    _stats.height = _rawImagePlanar->height();
    _stats.samples_per_channel = _stats.width * _stats.height;

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

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

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

void Solver::abort()
{
    _solver->abort();
}

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

bool Solver::updateHeader(QString &error)
{
    if(!_solver->solvingDone())
    {
        error = tr("Solving is not finished");
        return false;
    }

    FITSImage::Solution solution = getSolution();

    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.removeKeyword("RADECSYS");
    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"));
    bool ret = file.modifyFITSRecords(&modify);
    if(!ret)error = tr("Failed to update file header");
    return ret;
}

void Solver::setParameters(Parameters::ParametersProfile profile)
{
    auto profileParam = _solver->getBuiltInProfiles().at(profile);
    profileParam.partition = false;
    _solver->setParameters(profileParam);
}

void Solver::setParameters(const Parameters &parameters)
{
    auto profile = parameters;
    profile.partition = false;
    _solver->setParameters(profile);
}

void Solver::setSearchScale(double fovLow, double fowHigh, SSolver::ScaleUnits units)
{
    _solver->setSearchScale(fovLow, fowHigh, units);
}

void Solver::setSearchPosition(double ra, double dec)
{
    _solver->setSearchPositionRaDec(ra, dec);
}

void Solver::clearStartingPositionAndScale()
{
    _solver->clearSearchPosition();
    _solver->clearSearchScale();
}

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

}