aboutsummaryrefslogtreecommitdiff
path: root/src/cell_explorer.c
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2014-03-31 16:04:54 +0200
committerThomas White <taw@physics.org>2014-03-31 16:04:54 +0200
commitc1ba58ea57d0104bd5612cd4b9bbd714eba11dd6 (patch)
treea6088ca1978ed36f542571b458a1dfbcc2d8205a /src/cell_explorer.c
parent93291f2fadccf00c6134cf9ebafcd7c6c0986519 (diff)
cell_explorer: Add peak fitting
Diffstat (limited to 'src/cell_explorer.c')
-rw-r--r--src/cell_explorer.c199
1 files changed, 199 insertions, 0 deletions
diff --git a/src/cell_explorer.c b/src/cell_explorer.c
index 6714c92e..3b640f99 100644
--- a/src/cell_explorer.c
+++ b/src/cell_explorer.c
@@ -38,6 +38,7 @@
#include <gtk/gtk.h>
#include <math.h>
#include <gdk/gdkkeysyms.h>
+#include <gsl/gsl_multifit_nlin.h>
#include "stream.h"
#include "image.h"
@@ -92,6 +93,11 @@ typedef struct {
double press_min;
int sel;
+ int have_fit;
+ double fit_a;
+ double fit_b;
+ double fit_c;
+
struct _cellwindow *parent;
} HistoBox;
@@ -381,6 +387,26 @@ static gboolean draw_sig(GtkWidget *da, GdkEventExpose *event, HistoBox *b)
cairo_fill(cr);
}
+ /* Draw fitted curve */
+ if ( b->have_fit ) {
+
+ double A, B, C;
+ A = b->fit_a / max; B = b->fit_b; C = b->fit_c;
+ cairo_new_path(cr);
+
+ /* In the current coordinate system, 0,0 is the bottom left
+ * of the graph (coordinates b->dmin, 0), and 1,1 is the top
+ * right of the graph (coordinate b->dmax, max) */
+ for ( i=0; i<300; i++ ) {
+ double xd = (double)i/300.0;
+ double x = b->dmin + xd*(b->dmax - b->dmin);
+ cairo_line_to(cr, xd, A*exp(-(x-B)*(x-B)/(C*C)));
+ }
+ cairo_set_source_rgba(cr, 1.0, 0.3, 0.0, 1.0);
+ cairo_set_line_width(cr, 0.005);
+ cairo_stroke(cr);
+ }
+
cairo_restore(cr);
draw_axis(cr, b, width, height);
@@ -679,6 +705,172 @@ static gint quit_sig(GtkWidget *widget, CellWindow *w)
}
+struct gaussian_data
+{
+ size_t n;
+ double gstep;
+ double min;
+ int *data;
+};
+
+
+static int gaussian_f(const gsl_vector *p, void *vp, gsl_vector *f)
+{
+ struct gaussian_data *d = vp;
+ int i;
+
+ double a = gsl_vector_get(p, 0);
+ double b = gsl_vector_get(p, 1);
+ double c = gsl_vector_get(p, 2);
+
+ for ( i=0; i<d->n; i++ ) {
+ double x = i*d->gstep + d->min;
+ double ycalc = a*exp(-pow(x-b, 2)/(c*c));
+ gsl_vector_set(f, i, ycalc - d->data[i]);
+ }
+
+ return GSL_SUCCESS;
+}
+
+
+static int gaussian_df(const gsl_vector *p, void *vp, gsl_matrix *J)
+{
+ struct gaussian_data *d = vp;
+ double a = gsl_vector_get(p, 0);
+ double b = gsl_vector_get(p, 1);
+ double c = gsl_vector_get(p, 2);
+ int i;
+
+ for ( i=0; i<d->n; i++ ) {
+ double x = i*d->gstep + d->min;
+ double ycalc = a*exp(-pow(x-b, 2)/(c*c));
+ gsl_matrix_set(J, i, 0, ycalc/a);
+ gsl_matrix_set(J, i, 1, 2.0*ycalc*pow(c, -2)*(x-b));
+ gsl_matrix_set(J, i, 2, 2.0*ycalc*pow(c, -3)*(x-b)*(x-b));
+ }
+
+ return GSL_SUCCESS;
+}
+
+
+static int gaussian_fdf(const gsl_vector *x, void *data, gsl_vector *f,
+ gsl_matrix *J)
+{
+ gaussian_f(x, data, f);
+ gaussian_df(x, data, J);
+ return GSL_SUCCESS;
+}
+
+
+static void fit_param(HistoBox *h)
+{
+ gsl_multifit_fdfsolver *s;
+ gsl_vector *v;
+ gsl_multifit_function_fdf f;
+ struct gaussian_data params;
+ double min, max, gstep;
+ int n, n_iter, status, i, nm;
+ int lowest_bin, highest_bin;
+ int *data_p, *data_a, *data_b, *data_c, *data_i, *data_f;
+ int *data_r, *data_h;
+
+ if ( h->sel1 > h->sel2 ) {
+ min = h->sel2;
+ max = h->sel1;
+ } else {
+ min = h->sel1;
+ max = h->sel2;
+ }
+ gstep = (h->max - h->min)/h->n;
+ n = (max-min)/gstep;
+
+ /* How many bins fit entirely within the selected range? */
+ lowest_bin = floor((min - h->min)/gstep)+1;
+ highest_bin = floor((max - h->min)/gstep)-1;
+ params.n = 1 + highest_bin - lowest_bin;
+ if ( highest_bin-lowest_bin < 3 ) {
+ ERROR("Not enough bins.\n");
+ return;
+ }
+ params.min = h->min + (lowest_bin+0.5)*gstep;
+ params.gstep = gstep;
+
+ data_p = multihistogram_get_data(h->h, CAT_P);
+ data_a = multihistogram_get_data(h->h, CAT_A);
+ data_b = multihistogram_get_data(h->h, CAT_B);
+ data_c = multihistogram_get_data(h->h, CAT_C);
+ data_i = multihistogram_get_data(h->h, CAT_I);
+ data_f = multihistogram_get_data(h->h, CAT_F);
+ data_h = multihistogram_get_data(h->h, CAT_H);
+ data_r = multihistogram_get_data(h->h, CAT_R);
+ params.data = malloc(params.n*sizeof(int));
+ if ( params.data == NULL ) return;
+ nm = 0;
+ for ( i=0; i<params.n; i++ ) {
+ int j = i+lowest_bin;
+ if ( (j < 0) || (j >= h->n) ) {
+ params.data[i] = 0;
+ continue;
+ }
+ params.data[i] = data_p[j] + data_a[j] + data_b[j] + data_c[j]
+ + data_i[j] + data_f[j] + data_h[j] + data_r[j];
+ if ( params.data[i] > nm ) nm = params.data[i];
+ }
+
+ s = gsl_multifit_fdfsolver_alloc(gsl_multifit_fdfsolver_lmsder, n, 3);
+
+ v = gsl_vector_alloc(3);
+ gsl_vector_set(v, 0, nm);
+ gsl_vector_set(v, 1, min+(max-min)/2.0);
+ gsl_vector_set(v, 2, (max-min)/5.0);
+
+ f.f = gaussian_f;
+ f.df = gaussian_df;
+ f.fdf = gaussian_fdf;
+ f.n = n;
+ f.p = 3;
+ f.params = &params;
+ gsl_multifit_fdfsolver_set(s, &f, v);
+
+ n_iter = 0;
+ do {
+ n_iter++;
+ status = gsl_multifit_fdfsolver_iterate(s);
+ if ( status ) break;
+
+ status = gsl_multifit_test_delta(s->dx, s->x, 0.001, 0.001);
+ } while ( (status == GSL_CONTINUE) && (n_iter < 10));
+
+ STATUS("Fitted: %.2f %.2f %.2f after %i iterations\n",
+ gsl_vector_get(s->x, 0), gsl_vector_get(s->x, 1),
+ gsl_vector_get(s->x, 2), n_iter);
+ h->have_fit = 1;
+ h->fit_a = gsl_vector_get(s->x, 0);
+ h->fit_b = gsl_vector_get(s->x, 1);
+ h->fit_c = gsl_vector_get(s->x, 2);
+
+ free(params.data);
+ gsl_multifit_fdfsolver_free(s);
+ gsl_vector_free(v);
+}
+
+
+static gint fit_sig(GtkWidget *widget, CellWindow *w)
+{
+ if ( w->hist_a->show_sel ) fit_param(w->hist_a);
+ if ( w->hist_b->show_sel ) fit_param(w->hist_b);
+ if ( w->hist_c->show_sel ) fit_param(w->hist_c);
+ if ( w->hist_al->show_sel ) fit_param(w->hist_al);
+ if ( w->hist_be->show_sel ) fit_param(w->hist_be);
+ if ( w->hist_ga->show_sel ) fit_param(w->hist_ga);
+
+ /* FIXME: If all six worked, show values in cut+paste window */
+
+ redraw_all(w);
+
+ return TRUE;
+}
+
static gint about_sig(GtkWidget *widget, CellWindow *w)
{
GtkWidget *window;
@@ -720,6 +912,9 @@ static void add_menu_bar(CellWindow *w, GtkWidget *vbox)
"<menu name=\"file\" action=\"FileAction\">"
" <menuitem name=\"quit\" action=\"QuitAction\" />"
"</menu>"
+ "<menu name=\"tools\" action=\"ToolsAction\" >"
+ " <menuitem name=\"fit\" action=\"FitCellAction\" />"
+ "</menu>"
"<menu name=\"help\" action=\"HelpAction\">"
" <menuitem name=\"about\" action=\"AboutAction\" />"
"</menu>"
@@ -731,6 +926,10 @@ static void add_menu_bar(CellWindow *w, GtkWidget *vbox)
{ "QuitAction", GTK_STOCK_QUIT, "_Quit", NULL, NULL,
G_CALLBACK(quit_sig) },
+ { "ToolsAction", NULL, "_Tools", NULL, NULL, NULL },
+ { "FitCellAction", NULL, "_Fit cell", "<Control>F", NULL,
+ G_CALLBACK(fit_sig) },
+
{ "HelpAction", NULL, "_Help", NULL, NULL, NULL },
{ "AboutAction", GTK_STOCK_ABOUT, "_About", NULL, NULL,
G_CALLBACK(about_sig) },