StarFit class that can fit Gauss model over star image

2019-10-05 12:24:52 +02:00
parent 1eb70fe257
commit 702d2111cc
3 changed files with 213 additions and 3 deletions
@@ -0,0 +1,179 @@
+#include "starfit.h"
+#include <gsl/gsl_blas.h>
+#include <QDebug>
+#include <iostream>
+
+const int PARAM_AM = 0;
+const int PARAM_X0 = 1;
+const int PARAM_Y0 = 2;
+const int PARAM_SX = 3;
+const int PARAM_SY = 4;
+
+const int MAX_ITER = 20;
+const double TOL = 1.0e-3;
+
+struct StarData
+{
+    size_t size;
+    std::vector<double> val;
+};
+
+// a * exp(-0.5*((x-x0)/sx)^2 + ((y-y0)/sy)^2)
+double gauss_model(double a, double x0, double y0, double sx, double sy, double x, double y)
+{
+    double _x = (x-x0)/sx;
+    double _y = (y-y0)/sy;
+    return a*exp(-0.5*(_x*_x + _y*_y));
+}
+
+int func_f(const gsl_vector *X, void *params, gsl_vector *f)
+{
+    StarData *d = static_cast<StarData*>(params);
+    double am = gsl_vector_get(X, PARAM_AM);
+    double x0 = gsl_vector_get(X, PARAM_X0);
+    double y0 = gsl_vector_get(X, PARAM_Y0);
+    double sx = gsl_vector_get(X, PARAM_SX);
+    double sy = gsl_vector_get(X, PARAM_SY);
+
+    int i = 0;
+    for(size_t y=0;y<d->size;y++)
+    {
+        for(size_t x=0;x<d->size;x++)
+        {
+            double v = gauss_model(am, x0, y0, sx, sy, x, y);
+            gsl_vector_set(f, i, d->val[i] - v);
+            i++;
+        }
+    }
+
+    return GSL_SUCCESS;
+}
+
+int func_df(const gsl_vector *X, void *params, gsl_matrix *J)
+{
+    StarData *d = static_cast<StarData*>(params);
+    double am = gsl_vector_get(X, PARAM_AM);
+    double x0 = gsl_vector_get(X, PARAM_X0);
+    double y0 = gsl_vector_get(X, PARAM_Y0);
+    double sx = gsl_vector_get(X, PARAM_SX);
+    double sy = gsl_vector_get(X, PARAM_SY);
+
+    int i = 0;
+    for(size_t y=0;y<d->size;y++)
+    {
+        for(size_t x=0;x<d->size;x++)
+        {
+            double tx = x-x0;
+            double ty = y-y0;
+            double e = gauss_model(am, x0, y0, sx, sy, x, y);
+
+            gsl_matrix_set(J, i, PARAM_AM, -e/am);
+            gsl_matrix_set(J, i, PARAM_X0, -e*(tx/(sx*sx)));
+            gsl_matrix_set(J, i, PARAM_Y0, -e*(ty/(sy*sy)));
+            gsl_matrix_set(J, i, PARAM_SX, -e*(tx*tx/(sx*sx*sx)));
+            gsl_matrix_set(J, i, PARAM_SY, -e*(ty*ty/(sy*sy*sy)));
+            i++;
+        }
+    }
+
+    return GSL_SUCCESS;
+}
+
+//int func_fvv(const gsl_vector *x, const gsl_vector * v, void *params, gsl_vector *fvv)
+//{
+//    return GSL_SUCCESS;
+//}
+
+void callback(const size_t iter, void *, const gsl_multifit_nlinear_workspace *w)
+{
+    double rcond;
+    gsl_vector *x = gsl_multifit_nlinear_position(w);
+    gsl_multifit_nlinear_rcond(&rcond, w);
+    qDebug() << "Iter:" << iter << "Am:" << gsl_vector_get(x, PARAM_AM)
+                << "X0:" << gsl_vector_get(x, PARAM_X0)
+                << "Y0:" << gsl_vector_get(x, PARAM_Y0)
+                << "SX:" << gsl_vector_get(x, PARAM_SX)
+                << "SY:" << gsl_vector_get(x, PARAM_SY)
+                << "J(X):" << 1.0/rcond
+                << "av:" << gsl_multifit_nlinear_avratio(w);
+
+    gsl_matrix *j = gsl_multifit_nlinear_jac(w);
+    QString r = "=[";
+    for(int i=5*13;i<6*13;i++)
+        r += QString("%1, ").arg(gsl_matrix_get(j, i, PARAM_X0));
+    r += "];";
+    qDebug() << "J:" << r;
+}
+
+Star::Star()
+{
+    m_am = m_x = m_y = m_sx = m_sy = NAN;
+}
+
+StarFit::StarFit(int size)
+{
+    m_size = size;
+    m_fdf_params = gsl_multifit_nlinear_default_parameters();
+    m_fdf_params.trs = gsl_multifit_nlinear_trs_lmaccel;
+
+    m_fdf.f = func_f;
+    m_fdf.df = func_df;
+    m_fdf.fvv = nullptr;
+    m_fdf.n = size*size;
+    m_fdf.p = 5;//number of model parameters amplitude, x, y, fwhm_x, fwhm_y
+}
+
+StarFit::~StarFit()
+{
+}
+
+Star StarFit::fitStar(RawImageAbs *image, const Peak &peak)
+{
+
+}
+
+Star StarFit::fitStar(std::vector<double> data)
+{
+    Star star;
+    StarData d;
+    d.size = m_size;
+    d.val = data;
+    m_fdf.params = &d;
+    int info;
+
+    gsl_vector *start = gsl_vector_alloc(m_fdf.p);
+    gsl_vector_set(start, PARAM_AM, 1000.0);
+    gsl_vector_set(start, PARAM_X0, m_size/2);
+    gsl_vector_set(start, PARAM_Y0, m_size/2);
+    gsl_vector_set(start, PARAM_SX, 1.0);
+    gsl_vector_set(start, PARAM_SY, 1.0);
+
+    gsl_multifit_nlinear_workspace *workspace = gsl_multifit_nlinear_alloc(gsl_multifit_nlinear_trust, &m_fdf_params, m_fdf.n, m_fdf.p);
+
+    gsl_multifit_nlinear_init(start, &m_fdf, workspace);
+    gsl_vector *f = gsl_multifit_nlinear_residual(workspace);
+
+    double cost, cost0;
+    gsl_blas_ddot(f, f, &cost0);
+    int ret = gsl_multifit_nlinear_driver(MAX_ITER, TOL, TOL, TOL, nullptr, nullptr, &info, workspace);
+
+    gsl_blas_ddot(f, f, &cost);
+
+    qDebug() << "finished" << ret << info << cost0 << cost;
+
+    if(ret==0)
+    {
+        gsl_vector *y = gsl_multifit_nlinear_position(workspace);
+        star.m_am = gsl_vector_get(y, PARAM_AM);
+        star.m_x = gsl_vector_get(y, PARAM_X0);
+        star.m_y = gsl_vector_get(y, PARAM_Y0);
+        star.m_sx = gsl_vector_get(y, PARAM_SX);
+        star.m_sy = gsl_vector_get(y, PARAM_SY);
+
+    }
+
+    gsl_vector_free(start);
+    gsl_multifit_nlinear_free(workspace);
+
+    return star;
+}
@@ -0,0 +1,29 @@
+#ifndef STARFIT_H
+#define STARFIT_H
+
+#include "rawimage.h"
+#include <gsl/gsl_multifit_nlinear.h>
+
+struct Star
+{
+    double m_am;
+    double m_x,m_y;
+    double m_sx,m_sy;
+    Star();
+};
+
+class StarFit
+{
+    int m_size;
+    gsl_multifit_nlinear_fdf m_fdf;
+    gsl_multifit_nlinear_parameters m_fdf_params;
+    gsl_vector *m_vector;
+public:
+    StarFit(int size);
+    ~StarFit();
+    Star fitStar(RawImageAbs *image, const Peak &peak);
+    Star fitStar(std::vector<double> data);
+    //void fitStars(RawImageAbs *image, const std::vector<Peak> &peaks);
+};
+
+#endif // STARFIT_H
@@ -15,7 +15,7 @@ CONFIG += c++11

 QMAKE_CXXFLAGS += -fopenmp

-LIBS += -lraw -lexif -lcfitsio -fopenmp
+LIBS += -lraw -lexif -lcfitsio -lgsl -lgslcblas -fopenmp

 win32:LIBS += -LC:\msys64\mingw64\lib -LC:\msys64\mingw64\bin
 win32:INCLUDEPATH += C:\msys64\mingw64\include\ C:\msys64\mingw64\include\cfitsio
@@ -26,7 +26,8 @@ SOURCES += main.cpp\
    imageringlist.cpp \
    database.cpp \
    loadrunable.cpp \
-    imageinfo.cpp
+    imageinfo.cpp \
+    starfit.cpp

 HEADERS  += mainwindow.h \
    imagescrollarea.h \
@@ -34,4 +35,5 @@ HEADERS  += mainwindow.h \
    database.h \
    loadrunable.h \
    imageinfo.h \
-    rawimage.h
+    rawimage.h \
+    starfit.h