diff --git a/starfit.cpp b/starfit.cpp new file mode 100644 index 0000000..ddf6ac5 --- /dev/null +++ b/starfit.cpp @@ -0,0 +1,179 @@ +#include "starfit.h" +#include +#include +#include + +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 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(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;ysize;y++) + { + for(size_t x=0;xsize;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(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;ysize;y++) + { + for(size_t x=0;xsize;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 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; +} diff --git a/starfit.h b/starfit.h new file mode 100644 index 0000000..7ca22e8 --- /dev/null +++ b/starfit.h @@ -0,0 +1,29 @@ +#ifndef STARFIT_H +#define STARFIT_H + +#include "rawimage.h" +#include + +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 data); + //void fitStars(RawImageAbs *image, const std::vector &peaks); +}; + +#endif // STARFIT_H diff --git a/tenmon.pro b/tenmon.pro index 160e1ee..c2b28a7 100644 --- a/tenmon.pro +++ b/tenmon.pro @@ -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