From da79197376d17a42f141a0bfa8ea89e135665376 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Du=C5=A1an=20Poizl?= Date: Tue, 17 Sep 2024 23:05:27 +0200 Subject: [PATCH] Initial stellarsolver implementation --- CMakeLists.txt | 9 ++++ imageinfo.cpp | 14 +++---- imageinfo.h | 12 +++--- imagewidget.cpp | 2 +- imagewidget.h | 6 +-- loadrunable.cpp | 29 ++++++++----- loadrunable.h | 2 +- scriptengine.cpp | 12 +++--- scriptengine.h | 2 +- solver.cpp | 107 +++++++++++++++++++++++++++++++++++++++++++++++ solver.h | 22 ++++++++++ 11 files changed, 182 insertions(+), 35 deletions(-) create mode 100644 solver.cpp create mode 100644 solver.h diff --git a/CMakeLists.txt b/CMakeLists.txt index 41cf88f..5c61ab4 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -25,6 +25,7 @@ find_library(FITS_LIB cfitsio REQUIRED) find_library(RAW_LIB NAMES raw_r REQUIRED) find_library(WCS_LIB wcs wcslib REQUIRED) find_library(LCMS2_LIB lcms2 REQUIRED) +find_library(STELLARSOLVER_LIB stellarsolver) add_subdirectory(libXISF) @@ -92,6 +93,14 @@ if(LIBRAW_STATIC) target_link_libraries(tenmon PRIVATE jasper) endif() +find_path(STELLARSOLVER_INCLUDE stellarsolver.h PATH_SUFFIXES libstellarsolver) +if(STELLARSOLVER_INCLUDE AND STELLARSOLVER_LIB) + target_include_directories(tenmon PRIVATE ${STELLARSOLVER_INCLUDE}) + target_link_libraries(tenmon PRIVATE ${STELLARSOLVER_LIB}) + target_compile_definitions(tenmon PRIVATE "PLATESOLVER") + target_sources(tenmon PRIVATE solver.cpp solver.h) +endif(STELLARSOLVER_INCLUDE AND STELLARSOLVER_LIB) + option(FLATPAK "Flatpak build" OFF) if(FLATPAK) target_compile_definitions(tenmon PRIVATE FLATPAK) diff --git a/imageinfo.cpp b/imageinfo.cpp index 1f0d8a8..acacc9b 100644 --- a/imageinfo.cpp +++ b/imageinfo.cpp @@ -98,14 +98,14 @@ void ImageInfo::setInfo(const ImageInfoData &info) expandAll(); } -void WCSData::freeWCS() +void WCSDataT::freeWCS() { wcsvfree(&nwcs, &wcs); nwcs = 0; wcs = nullptr; } -WCSData::WCSData(int width, int height, char *header, int nrec) : +WCSDataT::WCSDataT(int width, int height, char *header, int nrec) : width(width), height(height) { @@ -121,7 +121,7 @@ WCSData::WCSData(int width, int height, char *header, int nrec) : freeWCS(); } -WCSData::WCSData(int width, int height, const QVector &header) : +WCSDataT::WCSDataT(int width, int height, const QVector &header) : width(width), height(height) { @@ -156,13 +156,13 @@ WCSData::WCSData(int width, int height, const QVector &header) : freeWCS(); } -WCSData::~WCSData() +WCSDataT::~WCSDataT() { if(wcs) freeWCS(); } -bool WCSData::pixelToWorld(const QPointF &pixel, SkyPoint &point) const +bool WCSDataT::pixelToWorld(const QPointF &pixel, SkyPoint &point) const { if(!valid())return false; @@ -181,7 +181,7 @@ bool WCSData::pixelToWorld(const QPointF &pixel, SkyPoint &point) const return false; } -bool WCSData::worldToPixel(const SkyPoint &point, QPointF &pixel) const +bool WCSDataT::worldToPixel(const SkyPoint &point, QPointF &pixel) const { if(!valid())return false; @@ -200,7 +200,7 @@ bool WCSData::worldToPixel(const SkyPoint &point, QPointF &pixel) const return false; } -void WCSData::calculateBounds(double &minRa, double &maxRa, double &minDec, double &maxDec, double &crVal1, double &crVal2) const +void WCSDataT::calculateBounds(double &minRa, double &maxRa, double &minDec, double &maxDec, double &crVal1, double &crVal2) const { if(wcs == nullptr)return; diff --git a/imageinfo.h b/imageinfo.h index e6f3ac6..a422055 100644 --- a/imageinfo.h +++ b/imageinfo.h @@ -35,7 +35,7 @@ public: QString toString() const; }; -class WCSData +class WCSDataT { int nwcs = 0; struct wcsprm *wcs = nullptr; @@ -43,10 +43,10 @@ class WCSData int height; void freeWCS(); public: - WCSData(int width, int height, char *header, int nrec); - WCSData(int width, int height, const QVector &header); - WCSData(const WCSData &) = delete; - ~WCSData(); + WCSDataT(int width, int height, char *header, int nrec); + WCSDataT(int width, int height, const QVector &header); + WCSDataT(const WCSDataT &) = delete; + ~WCSDataT(); bool pixelToWorld(const QPointF &pixel, SkyPoint &point) const; bool worldToPixel(const SkyPoint &point, QPointF &pixel) const; void calculateBounds(double &minRa, double &maxRa, double &minDec, double &maxDec, double &crVal1, double &crVal2) const; @@ -57,7 +57,7 @@ struct ImageInfoData { QVector fitsHeader; QVector> info; - std::shared_ptr wcs; + std::shared_ptr wcs; }; Q_DECLARE_METATYPE(ImageInfoData); diff --git a/imagewidget.cpp b/imagewidget.cpp index 1b656f4..e39f94c 100644 --- a/imagewidget.cpp +++ b/imagewidget.cpp @@ -168,7 +168,7 @@ void ImageWidgetGL::setImage(std::shared_ptr image, int index) else setOffset(m_dx, m_dy); } -void ImageWidgetGL::setWCS(std::shared_ptr wcs) +void ImageWidgetGL::setWCS(std::shared_ptr wcs) { m_wcs = wcs; } diff --git a/imagewidget.h b/imagewidget.h index fbda655..f7faec4 100644 --- a/imagewidget.h +++ b/imagewidget.h @@ -20,7 +20,7 @@ public: virtual ~ImageWidget(){} virtual void setImage(std::shared_ptr image, int index) = 0; - virtual void setWCS(std::shared_ptr wcs) = 0; + virtual void setWCS(std::shared_ptr wcs) = 0; virtual void zoom(int zoom, const QPointF &mousePos = QPointF()) = 0; virtual void bestFit() = 0; @@ -65,7 +65,7 @@ class ImageWidgetGL : public QOpenGLWidget, public ImageWidget std::unique_ptr m_lut; GLuint m_debayerTex = 0; std::shared_ptr m_rawImage; - std::shared_ptr m_wcs; + std::shared_ptr m_wcs; int m_width, m_height; int m_imgWidth = -1, m_imgHeight = -1; int m_currentImg = 0; @@ -95,7 +95,7 @@ public: explicit ImageWidgetGL(Database *database, QWidget *parent = nullptr); ~ImageWidgetGL() override; void setImage(std::shared_ptr image, int index) override; - void setWCS(std::shared_ptr wcs) override; + void setWCS(std::shared_ptr wcs) override; void zoom(int zoom, const QPointF &mousePos = QPointF()) override; void bestFit() override; void allocateThumbnails(const QStringList &paths) override; diff --git a/loadrunable.cpp b/loadrunable.cpp index c80427d..67f8b6d 100644 --- a/loadrunable.cpp +++ b/loadrunable.cpp @@ -176,14 +176,14 @@ int loadFITSHeader(fitsfile *file, ImageInfoData &info) fits_hdr2str(file, TRUE, (char**)exclist, 2, &header, &nrec, &status); if(status == 0) { - info.wcs = std::make_shared(naxes[0], naxes[1], header, nrec); + info.wcs = std::make_shared(naxes[0], naxes[1], header, nrec); if(!info.wcs->valid())info.wcs.reset(); } fits_free_memory(header, &status); return status; } -bool loadFITS(const QString path, ImageInfoData &info, std::shared_ptr &image) +bool loadFITS(const QString path, ImageInfoData &info, std::shared_ptr &image, bool planar) { fitsfile *file; int status = 0; @@ -261,7 +261,7 @@ bool loadFITS(const QString path, ImageInfoData &info, std::shared_ptr s[i] -= INT16_MIN; } - if(img.channels() == 1) + if(img.channels() == 1 || planar) image = std::make_shared(std::move(img)); else image = RawImage::fromPlanar(img); @@ -294,7 +294,7 @@ bool loadFITS(const QString path, ImageInfoData &info, std::shared_ptr return true; } -bool loadXISF(const QString &path, ImageInfoData &info, std::shared_ptr &image) +bool loadXISF(const QString &path, ImageInfoData &info, std::shared_ptr &image, bool planar) { try { @@ -314,7 +314,7 @@ bool loadXISF(const QString &path, ImageInfoData &info, std::shared_ptr(xisfImage.width(), xisfImage.height(), info.fitsHeader); + info.wcs = std::make_shared(xisfImage.width(), xisfImage.height(), info.fitsHeader); info.info.append({QObject::tr("Width"), QString::number(xisfImage.width())}); info.info.append({QObject::tr("Height"), QString::number(xisfImage.height())}); if(!info.wcs->valid())info.wcs.reset(); @@ -341,7 +341,16 @@ bool loadXISF(const QString &path, ImageInfoData &info, std::shared_ptr(tmpImage.width(), tmpImage.height(), tmpImage.channelCount(), type); + std::memcpy(image->data(), tmpImage.imageData(), tmpImage.imageDataSize()); + } + else + { + image = RawImage::fromPlanar(tmpImage.imageData(), tmpImage.width(), tmpImage.height(), tmpImage.channelCount(), type); + } + image->setICCProfile(tmpImage.iccProfile()); return true; } @@ -471,7 +480,7 @@ bool readXISFHeader(const QString &path, ImageInfoData &info) { info.fitsHeader.append(prop); } - info.wcs = std::make_shared(image.width(), image.height(), info.fitsHeader); + info.wcs = std::make_shared(image.width(), image.height(), info.fitsHeader); if(!info.wcs->valid())info.wcs.reset(); } catch (LibXISF::Error &err) @@ -482,7 +491,7 @@ bool readXISFHeader(const QString &path, ImageInfoData &info) return true; } -bool loadImage(const QString &path, ImageInfoData &info, std::shared_ptr &rawImage) +bool loadImage(const QString &path, ImageInfoData &info, std::shared_ptr &rawImage, bool planar) { bool ret = false; QElapsedTimer timer; @@ -494,12 +503,12 @@ bool loadImage(const QString &path, ImageInfoData &info, std::shared_ptr &rawImage); +bool loadImage(const QString &path, ImageInfoData &info, std::shared_ptr &rawImage, bool planar = false); class Image; diff --git a/scriptengine.cpp b/scriptengine.cpp index 14c3f0f..670dcd6 100644 --- a/scriptengine.cpp +++ b/scriptengine.cpp @@ -359,7 +359,7 @@ bool File::modifyFITSRecords(const FITSRecordModify *modify) fits_get_num_hdus(file, &num, &status); if(status) { - _engine->newMessage("Failed to open FITS file", true); + if(_engine)_engine->newMessage("Failed to open FITS file", true); return false; } int imgtype; @@ -412,7 +412,7 @@ bool File::modifyFITSRecords(const FITSRecordModify *modify) break; } default: - _engine->newMessage("Unknown type for KEY " + record.key, true); + if(_engine)_engine->newMessage("Unknown type for KEY " + record.key, true); return false; break; } @@ -420,7 +420,7 @@ bool File::modifyFITSRecords(const FITSRecordModify *modify) { char error[100]; fits_get_errstatus(status, error); - _engine->newMessage(QString("Error when updating KEY %1 %2").arg(record.key).arg(error), true); + if(_engine)_engine->newMessage(QString("Error when updating KEY %1 %2").arg(record.key).arg(error), true); return false; } } @@ -456,7 +456,7 @@ bool File::modifyFITSRecords(const FITSRecordModify *modify) break; } default: - _engine->newMessage("Unknown type for KEY " + record.key, true); + if(_engine)_engine->newMessage("Unknown type for KEY " + record.key, true); return false; break; } @@ -464,7 +464,7 @@ bool File::modifyFITSRecords(const FITSRecordModify *modify) { char error[100]; fits_get_errstatus(status, error); - _engine->newMessage(QString("Error when adding KEY {} {}").arg(record.key).arg(error), true); + if(_engine)_engine->newMessage(QString("Error when adding KEY {} {}").arg(record.key).arg(error), true); return false; } } @@ -496,7 +496,7 @@ bool File::modifyFITSRecords(const FITSRecordModify *modify) } catch(LibXISF::Error &err) { - _engine->newMessage("Failed to modify file " + _path + " " + err.what(), true); + if(_engine)_engine->newMessage("Failed to modify file " + _path + " " + err.what(), true); return false; } } diff --git a/scriptengine.h b/scriptengine.h index 8c1b497..36f4f9d 100644 --- a/scriptengine.h +++ b/scriptengine.h @@ -76,7 +76,7 @@ class FITSRecordModify; class File : public QObject { Q_OBJECT - ScriptEngine *_engine; + ScriptEngine *_engine = nullptr; QString _path; QString _root; QFileInfo _info; diff --git a/solver.cpp b/solver.cpp new file mode 100644 index 0000000..f82de3a --- /dev/null +++ b/solver.cpp @@ -0,0 +1,107 @@ +#include "solver.h" +#include +#include +#include +#include "rawimage.h" +#include "loadrunable.h" +#include "scriptengine.h" + +Solver::Solver(QObject *parent) : QObject(parent) +{ + _solver = std::make_unique(); + + QStringList indexFolder = StellarSolver::getDefaultIndexFolderPaths(); + + _solver->setProperty("ProcessType", SSolver::SOLVE); + _solver->setIndexFolderPaths(indexFolder); + _solver->setParameterProfile(SSolver::Parameters::ALL_STARS); +} + +Solver::~Solver() +{ +} + +bool Solver::solveImage(const QString &path) +{ + ImageInfoData info; + std::shared_ptr rawImage; + if(loadImage(path, info, rawImage, true)) + { + _path = path; + 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; + + _solver->loadNewImageBuffer(_stats, (const uint8_t*)rawImage->data()); + return _solver->solve(); + } + return false; +} + +FITSImage::Solution Solver::getSolution() const +{ + return _solver->getSolution(); +} + +QString Solver::errorMessage() const +{ + return _error; +} + +void Solver::updateHeader() +{ + 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); +} diff --git a/solver.h b/solver.h new file mode 100644 index 0000000..422e878 --- /dev/null +++ b/solver.h @@ -0,0 +1,22 @@ +#ifndef SOLVER_H +#define SOLVER_H + +#include + +class Solver : public QObject +{ + Q_OBJECT + std::unique_ptr _solver; + FITSImage::Statistic _stats; + QString _path; + QString _error; +public: + explicit Solver(QObject *parent = nullptr); + ~Solver(); + bool solveImage(const QString &path); + FITSImage::Solution getSolution() const; + QString errorMessage() const; + void updateHeader(); +}; + +#endif // SOLVER_H