Fix error in starfit model, preparation for angled version

loadrunable.cpp

@@ -6,9 +6,11 @@
 #include <QPainter>
 #include <QElapsedTimer>
 #include <QDebug>
+#include <iostream>
 #include <libexif/exif-data.h>
 #include <fitsio2.h>
 #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<Peak> &peaks)
     img = pix.toImage();
 }
 
+void drawStars(QImage &img, const std::vector<Star> &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<double> &data, const Star &star)
+{
+    QString d = "d=[";
+    QString m = "m=[";
+    for(int y=0; y<radius; y++)
+    {
+        for(int x=0; x<radius; x++)
+        {
+            d += QString::number(data[y*radius+x]) + ",";
+            m += QString::number(gauss_model(star.m_am, star.m_x, star.m_y, star.m_sx, star.m_sy, x, y)) + ",";
+        }
+        d += ";";
+        m += ";";
+    }
+    d += "];";
+    m += "];";
+    //std::cout << star.m_am << " " << star.m_sx << star.m_sy << std::endl;
+    std::cout << d.toStdString() << std::endl;
+    std::cout << m.toStdString() << std::endl << std::endl;
+}
+
 bool loadRAW(QString path,  ImageInfoData &info, RawImageAbs **image, QImage *qimage)
 {
     if(!image && !qimage)
@@ -271,14 +306,41 @@ void LoadRunable::run()
                 rawImage->quarter();
                 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();
+            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<Star> stars;
+                for(uint i=0; i<peaks.size(); i++)
+                {
+                    Peak p = peaks[i];
+                    std::vector<double> 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();
         }
     }
 

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<double> &r) const = 0;
     virtual int findPeaks(uint64_t background, double distance, std::vector<Peak> &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<double> &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<y+h;i++)
+            for(int o=x;o<x+w;o++)
+                r[d++] = pixel(o, i);
+    }
     int findPeaks(uint64_t background, double distance, std::vector<Peak> &peaks) const
     {
         std::vector<Peak> tmpPeaks;
@@ -165,8 +180,9 @@ public:
         }
         return num;
     }
-    void medianFilter()
+    RawImageAbs* medianFilter() const
     {
+        RawImage<T> *ret = new RawImage<T>;
         std::vector<T> 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()
     {

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<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);
+    double th = gsl_vector_get(X, PARAM_TH);
+
+    int i = 0;
+    double a = sin(th);
+    double b = cos(th);
+    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*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<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);
+    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<double> data)
+Star StarFit::fitStar(const std::vector<double> &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);
+    double min = *std::min_element(data.begin(), data.end());
-    gsl_vector_set(start, PARAM_AM, 1000.0);
+    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<double> 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<double> 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);

starfit.h

@@ -4,25 +4,30 @@
 #include "rawimage.h"
 #include <gsl/gsl_multifit_nlinear.h>
 
+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<double> data);
+    Star fitStar(const std::vector<double> &data, bool angle);
     //void fitStars(RawImageAbs *image, const std::vector<Peak> &peaks);
 };