Initial stellarsolver implementation

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)

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<FITSRecord> &header) :
+WCSDataT::WCSDataT(int width, int height, const QVector<FITSRecord> &header) :
     width(width),
     height(height)
 {
@@ -156,13 +156,13 @@ WCSData::WCSData(int width, int height, const QVector<FITSRecord> &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;
 

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);
+    WCSDataT(int width, int height, char *header, int nrec);
-    WCSData(int width, int height, const QVector<FITSRecord> &header);
+    WCSDataT(int width, int height, const QVector<FITSRecord> &header);
-    WCSData(const WCSData &) = delete;
+    WCSDataT(const WCSDataT &) = delete;
-    ~WCSData();
+    ~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<FITSRecord> fitsHeader;
     QVector<QPair<QString, QString>> info;
-    std::shared_ptr<WCSData> wcs;
+    std::shared_ptr<WCSDataT> wcs;
 };
 
 Q_DECLARE_METATYPE(ImageInfoData);

imagewidget.cpp

@@ -168,7 +168,7 @@ void ImageWidgetGL::setImage(std::shared_ptr<RawImage> image, int index)
     else setOffset(m_dx, m_dy);
 }
 
-void ImageWidgetGL::setWCS(std::shared_ptr<WCSData> wcs)
+void ImageWidgetGL::setWCS(std::shared_ptr<WCSDataT> wcs)
 {
     m_wcs = wcs;
 }

imagewidget.h

@@ -20,7 +20,7 @@ public:
     virtual ~ImageWidget(){}
 
     virtual void setImage(std::shared_ptr<RawImage> image, int index) = 0;
-    virtual void setWCS(std::shared_ptr<WCSData> wcs) = 0;
+    virtual void setWCS(std::shared_ptr<WCSDataT> 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<QOpenGLTexture> m_lut;
     GLuint m_debayerTex = 0;
     std::shared_ptr<RawImage> m_rawImage;
-    std::shared_ptr<WCSData> m_wcs;
+    std::shared_ptr<WCSDataT> 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<RawImage> image, int index) override;
-    void setWCS(std::shared_ptr<WCSData> wcs) override;
+    void setWCS(std::shared_ptr<WCSDataT> wcs) override;
     void zoom(int zoom, const QPointF &mousePos = QPointF()) override;
     void bestFit() override;
     void allocateThumbnails(const QStringList &paths) override;

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<WCSData>(naxes[0], naxes[1], header, nrec);
+        info.wcs = std::make_shared<WCSDataT>(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<RawImage> &image)
+bool loadFITS(const QString path, ImageInfoData &info, std::shared_ptr<RawImage> &image, bool planar)
 {
     fitsfile *file;
     int status = 0;
@@ -261,7 +261,7 @@ bool loadFITS(const QString path, ImageInfoData &info, std::shared_ptr<RawImage>
                     s[i] -= INT16_MIN;
             }
 
-            if(img.channels() == 1)
+            if(img.channels() == 1 || planar)
                 image = std::make_shared<RawImage>(std::move(img));
             else
                 image = RawImage::fromPlanar(img);
@@ -294,7 +294,7 @@ bool loadFITS(const QString path, ImageInfoData &info, std::shared_ptr<RawImage>
     return true;
 }
 
-bool loadXISF(const QString &path, ImageInfoData &info, std::shared_ptr<RawImage> &image)
+bool loadXISF(const QString &path, ImageInfoData &info, std::shared_ptr<RawImage> &image, bool planar)
 {
     try
     {
@@ -314,7 +314,7 @@ bool loadXISF(const QString &path, ImageInfoData &info, std::shared_ptr<RawImage
             info.fitsHeader.append(prop);
         }
 
-        info.wcs = std::make_shared<WCSData>(xisfImage.width(), xisfImage.height(), info.fitsHeader);
+        info.wcs = std::make_shared<WCSDataT>(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<RawImage
         }
         else if(tmpImage.channelCount() == 3 || tmpImage.channelCount() == 4)
         {
-            image = RawImage::fromPlanar(tmpImage.imageData(), tmpImage.width(), tmpImage.height(), tmpImage.channelCount(), type);
+            if(planar)
+            {
+                image = std::make_shared<RawImage>(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<WCSData>(image.width(), image.height(), info.fitsHeader);
+        info.wcs = std::make_shared<WCSDataT>(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> &rawImage)
+bool loadImage(const QString &path, ImageInfoData &info, std::shared_ptr<RawImage> &rawImage, bool planar)
 {
     bool ret = false;
     QElapsedTimer timer;
@@ -494,12 +503,12 @@ bool loadImage(const QString &path, ImageInfoData &info, std::shared_ptr<RawImag
     }
     else if(path.endsWith(".FIT", Qt::CaseInsensitive) || path.endsWith(".FITS", Qt::CaseInsensitive))
     {
-        ret = loadFITS(path, info, rawImage);
+        ret = loadFITS(path, info, rawImage, planar);
         qDebug() << "LoadFITS" << timer.elapsed();
     }
     else if(path.endsWith(".XISF", Qt::CaseInsensitive))
     {
-        ret = loadXISF(path, info, rawImage);
+        ret = loadXISF(path, info, rawImage, planar);
         qDebug() << "LoadXISF" << timer.elapsed();
     }
     else

loadrunable.h

@@ -10,7 +10,7 @@ class RawImage;
 
 bool readFITSHeader(const QString &path, ImageInfoData &info);
 bool readXISFHeader(const QString &path, ImageInfoData &info);
-bool loadImage(const QString &path, ImageInfoData &info, std::shared_ptr<RawImage> &rawImage);
+bool loadImage(const QString &path, ImageInfoData &info, std::shared_ptr<RawImage> &rawImage, bool planar = false);
 
 class Image;
 

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

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;

solver.cpp

@@ -0,0 +1,107 @@
+#include "solver.h"
+#include <fitsio.h>
+#include <wcslib/wcshdr.h>
+#include <wcslib/wcsutil.h>
+#include "rawimage.h"
+#include "loadrunable.h"
+#include "scriptengine.h"
+
+Solver::Solver(QObject *parent) : QObject(parent)
+{
+    _solver = std::make_unique<StellarSolver>();
+
+    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> 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);
+}

solver.h

@@ -0,0 +1,22 @@
+#ifndef SOLVER_H
+#define SOLVER_H
+
+#include <stellarsolver.h>
+
+class Solver : public QObject
+{
+    Q_OBJECT
+    std::unique_ptr<StellarSolver> _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