From 9dacf62226cfd1af957681374f7a230d971ab50f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Du=C5=A1an=20Poizl?= Date: Fri, 11 Oct 2019 13:16:25 +0200 Subject: [PATCH] Fix error in starfit model, preparation for angled version --- loadrunable.cpp | 68 +++++++++++++++++++++++++-- rawimage.h | 25 ++++++++-- starfit.cpp | 121 +++++++++++++++++++++++++++++++++++++++++------- starfit.h | 7 ++- 4 files changed, 198 insertions(+), 23 deletions(-) diff --git a/loadrunable.cpp b/loadrunable.cpp index 8c5fef3..76ad78c 100644 --- a/loadrunable.cpp +++ b/loadrunable.cpp @@ -6,9 +6,11 @@ #include #include #include +#include #include #include #include "rawimage.h" +#include "starfit.h" LoadRunable::LoadRunable(const QString &file, Image *receiver, AnalyzeLevel level) : m_file(file), @@ -40,6 +42,39 @@ void drawPeaks(QImage &img, const std::vector &peaks) img = pix.toImage(); } +void drawStars(QImage &img, const std::vector &stars) +{ + QPixmap pix = QPixmap::fromImage(img); + QPainter painter(&pix); + painter.setPen(Qt::red); + for(auto star : stars) + { + painter.drawEllipse(QPointF(star.m_x, star.m_y), star.m_sx, star.m_sy); + } + img = pix.toImage(); +} + +void printStarModel(int radius, const std::vector &data, const Star &star) +{ + QString d = "d=["; + QString m = "m=["; + for(int y=0; yquarter(); qDebug() << "quarter" << timer.restart(); } - rawImage->medianFilter(); + RawImageAbs *medianImage = rawImage->medianFilter(); qDebug() << "median" << timer.restart(); - int numPeaks = rawImage->findPeaks(median+stdDev, 20, peaks); + int numPeaks = medianImage->findPeaks(median+stdDev, 20, peaks); + delete medianImage; qDebug() << "peaks" << timer.restart(); - drawPeaks(img, peaks); + if(m_analyzeLevel == Peaks) + drawPeaks(img, peaks); qDebug() << "draw peaks" << timer.restart(); info.append(StringPair(QObject::tr("Peaks"), QString::number(numPeaks))); info.append(StringPair(QObject::tr("Peaks draw"), QString::number(peaks.size()))); + + if(m_analyzeLevel>= Stars) + { + const int radius = 13; + StarFit starFit(radius); + std::vector stars; + for(uint i=0; i r; + int x = p.x(); + int y = p.y(); + rawImage->rect(x, y, radius, radius, r); + Star star = starFit.fitStar(r, false); + if(star.valid()) + { + //printStarModel(radius, r, star); + star.m_x += x; + star.m_y += y; + stars.push_back(star); + } + } + drawStars(img, stars); + } + qDebug() << "Star fit" << timer.restart(); } } diff --git a/rawimage.h b/rawimage.h index 76d616b..097a100 100644 --- a/rawimage.h +++ b/rawimage.h @@ -34,8 +34,9 @@ class RawImageAbs public: virtual ~RawImageAbs(){} virtual bool imageStats(uint64_t *mean, double *stdDev, uint64_t *median, uint64_t *min, uint64_t *max) const = 0; + virtual void rect(int &x, int &y, int w, int h, std::vector &r) const = 0; virtual int findPeaks(uint64_t background, double distance, std::vector &peaks) const = 0; - virtual void medianFilter() = 0; + virtual RawImageAbs* medianFilter() const = 0; virtual void quarter() = 0; }; @@ -125,6 +126,20 @@ public: if(y>=m_height)y=0; return m_img[y*m_width+x]; } + void rect(int &x, int &y, int w, int h, std::vector &r) const + { + r.resize(w*h); + x -= w/2; + y -= h/2; + if(x<0)x = 0; + if(y<0)y = 0; + if(x+w >= m_width)x = m_width-w; + if(y+h >= m_height)y = m_height-h; + uint32_t d = 0; + for(int i=y;i &peaks) const { std::vector tmpPeaks; @@ -165,8 +180,9 @@ public: } return num; } - void medianFilter() + RawImageAbs* medianFilter() const { + RawImage *ret = new RawImage; std::vector tmp; tmp.resize(m_width*m_height); #pragma omp parallel for @@ -181,7 +197,10 @@ public: tmp[y*m_width+x] = array[4]; } } - m_img = std::move(tmp); + ret->m_width = m_width; + ret->m_height = m_height; + ret->m_img = std::move(tmp); + return ret; } void quarter() { diff --git a/starfit.cpp b/starfit.cpp index ddf6ac5..aa8c505 100644 --- a/starfit.cpp +++ b/starfit.cpp @@ -8,6 +8,7 @@ const int PARAM_X0 = 1; const int PARAM_Y0 = 2; const int PARAM_SX = 3; const int PARAM_SY = 4; +const int PARAM_TH = 5; const int MAX_ITER = 20; const double TOL = 1.0e-3; @@ -79,12 +80,84 @@ int func_df(const gsl_vector *X, void *params, gsl_matrix *J) return GSL_SUCCESS; } +int func_f_an(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); + double th = gsl_vector_get(X, PARAM_TH); + + int i = 0; + double a = sin(th); + double b = cos(th); + for(size_t y=0;ysize;y++) + { + for(size_t x=0;xsize;x++) + { + double v = gauss_model(am, x0, y0, sx, sy, x*b-y*a, x*a+y*b); + gsl_vector_set(f, i, d->val[i] - v); + i++; + } + } + + return GSL_SUCCESS; +} + +int func_df_af(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); + QString r = "Iter: " + QString::number(iter) + + " Am: " + QString::number(gsl_vector_get(x, PARAM_AM)) + + " X0: " + QString::number(gsl_vector_get(x, PARAM_X0)) + + " Y0: " + QString::number(gsl_vector_get(x, PARAM_Y0)) + + " SX: " + QString::number(gsl_vector_get(x, PARAM_SX)) + + " SY: " + QString::number(gsl_vector_get(x, PARAM_SY)) + + " J(X) :" + QString::number(1.0/rcond) + + " av: " + QString::number(gsl_multifit_nlinear_avratio(w)); + std::cout << r.toStdString() << std::endl; +} + +void callback_an(const size_t iter, void *, const gsl_multifit_nlinear_workspace *w) { double rcond; gsl_vector *x = gsl_multifit_nlinear_position(w); @@ -94,15 +167,9 @@ void callback(const size_t iter, void *, const gsl_multifit_nlinear_workspace *w << "Y0:" << gsl_vector_get(x, PARAM_Y0) << "SX:" << gsl_vector_get(x, PARAM_SX) << "SY:" << gsl_vector_get(x, PARAM_SY) + << "TH:" << gsl_vector_get(x, PARAM_TH) << "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() @@ -110,6 +177,11 @@ Star::Star() m_am = m_x = m_y = m_sx = m_sy = NAN; } +bool Star::valid() const +{ + return !isnan(m_am); +} + StarFit::StarFit(int size) { m_size = size; @@ -121,6 +193,12 @@ StarFit::StarFit(int size) m_fdf.fvv = nullptr; m_fdf.n = size*size; m_fdf.p = 5;//number of model parameters amplitude, x, y, fwhm_x, fwhm_y + + m_fdf_an.f = func_f_an; + m_fdf_an.df = nullptr; + m_fdf_an.fvv = nullptr; + m_fdf_an.n = size*size; + m_fdf_an.p = 6;//number of model parameters amplitude, x, y, sigma_x, sigma_y, angle } StarFit::~StarFit() @@ -132,25 +210,36 @@ Star StarFit::fitStar(RawImageAbs *image, const Peak &peak) } -Star StarFit::fitStar(std::vector data) +Star StarFit::fitStar(const std::vector &data, bool angle) { + gsl_multifit_nlinear_fdf *fdf = angle ? &m_fdf_an : &m_fdf; Star star; StarData d; + d.val = data; d.size = m_size; d.val = data; - m_fdf.params = &d; + fdf->params = &d; int info; - gsl_vector *start = gsl_vector_alloc(m_fdf.p); - gsl_vector_set(start, PARAM_AM, 1000.0); + double min = *std::min_element(data.begin(), data.end()); + double max = *std::max_element(data.begin(), data.end()) - min; + for(double &v : d.val) + { + v -= min; + } + + gsl_vector *start = gsl_vector_alloc(fdf->p); + gsl_vector_set(start, PARAM_AM, max); 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); + if(angle) + gsl_vector_set(start, PARAM_TH, 0.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_workspace *workspace = gsl_multifit_nlinear_alloc(gsl_multifit_nlinear_trust, &m_fdf_params, fdf->n, fdf->p); - gsl_multifit_nlinear_init(start, &m_fdf, workspace); + gsl_multifit_nlinear_init(start, fdf, workspace); gsl_vector *f = gsl_multifit_nlinear_residual(workspace); double cost, cost0; @@ -159,8 +248,6 @@ Star StarFit::fitStar(std::vector data) gsl_blas_ddot(f, f, &cost); - qDebug() << "finished" << ret << info << cost0 << cost; - if(ret==0) { gsl_vector *y = gsl_multifit_nlinear_position(workspace); @@ -169,7 +256,9 @@ Star StarFit::fitStar(std::vector data) 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); - + if(angle) + star.m_theta = gsl_vector_get(y, PARAM_TH); + //qDebug() << "finished" << star.m_am << star.m_sx << star.m_sy; } gsl_vector_free(start); diff --git a/starfit.h b/starfit.h index 7ca22e8..9292fc5 100644 --- a/starfit.h +++ b/starfit.h @@ -4,25 +4,30 @@ #include "rawimage.h" #include +double gauss_model(double a, double x0, double y0, double sx, double sy, double x, double y); + struct Star { double m_am; double m_x,m_y; double m_sx,m_sy; + double m_theta; Star(); + bool valid() const; }; class StarFit { int m_size; gsl_multifit_nlinear_fdf m_fdf; + gsl_multifit_nlinear_fdf m_fdf_an; 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); + Star fitStar(const std::vector &data, bool angle); //void fitStars(RawImageAbs *image, const std::vector &peaks); };