Added iterative peaksearch method

git-svn-id: svn://cook.msm.cam.ac.uk:745/diff-tomo/dtr@3 bf6ca9ba-c028-0410-8290-897cf20841d1

7 files changed, 552 insertions, 41 deletions

diff --git a/src/Makefile.am b/src/Makefile.am
index 1e532f1..1b29a7a 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -1,5 +1,5 @@
 bin_PROGRAMS = dtr
-dtr_SOURCES = main.c displaywindow.c trackball.c reflections.c readpng.c mrc.c imagedisplay.c utils.c itrans.c qdrp.c ipr.c
+dtr_SOURCES = main.c displaywindow.c trackball.c reflections.c readpng.c mrc.c imagedisplay.c utils.c itrans.c qdrp.c ipr.c peakdetect.c
 dtr_LDADD = @LIBS@ @GTK_LIBS@ -lm @GTKGLEXT_LIBS@ -lgsl -lgslcblas
 AM_CFLAGS = -Wall -g @CFLAGS@ @GTK_CFLAGS@ @GTKGLEXT_CFLAGS@
 AM_CPPFLAGS = -DDATADIR=\""$(datadir)"\"
diff --git a/src/Makefile.in b/src/Makefile.in
index d07e5be..d6994f2 100644
--- a/src/Makefile.in
+++ b/src/Makefile.in
@@ -1,8 +1,8 @@
-# Makefile.in generated by automake 1.9.6 from Makefile.am.
+# Makefile.in generated by automake 1.10 from Makefile.am.
 # @configure_input@
 
 # Copyright (C) 1994, 1995, 1996, 1997, 1998, 1999, 2000, 2001, 2002,
-# 2003, 2004, 2005  Free Software Foundation, Inc.
+# 2003, 2004, 2005, 2006  Free Software Foundation, Inc.
 # This Makefile.in is free software; the Free Software Foundation
 # gives unlimited permission to copy and/or distribute it,
 # with or without modifications, as long as this notice is preserved.
@@ -14,15 +14,11 @@
 
 @SET_MAKE@
 
-srcdir = @srcdir@
-top_srcdir = @top_srcdir@
 VPATH = @srcdir@
 pkgdatadir = $(datadir)/@PACKAGE@
 pkglibdir = $(libdir)/@PACKAGE@
 pkgincludedir = $(includedir)/@PACKAGE@
-top_builddir = ..
 am__cd = CDPATH="$${ZSH_VERSION+.}$(PATH_SEPARATOR)" && cd
-INSTALL = @INSTALL@
 install_sh_DATA = $(install_sh) -c -m 644
 install_sh_PROGRAM = $(install_sh) -c
 install_sh_SCRIPT = $(install_sh) -c
@@ -50,10 +46,11 @@ PROGRAMS = $(bin_PROGRAMS)
 am_dtr_OBJECTS = main.$(OBJEXT) displaywindow.$(OBJEXT) \
 	trackball.$(OBJEXT) reflections.$(OBJEXT) readpng.$(OBJEXT) \
 	mrc.$(OBJEXT) imagedisplay.$(OBJEXT) utils.$(OBJEXT) \
-	itrans.$(OBJEXT) qdrp.$(OBJEXT) ipr.$(OBJEXT)
+	itrans.$(OBJEXT) qdrp.$(OBJEXT) ipr.$(OBJEXT) \
+	peakdetect.$(OBJEXT)
 dtr_OBJECTS = $(am_dtr_OBJECTS)
 dtr_DEPENDENCIES =
-DEFAULT_INCLUDES = -I. -I$(srcdir) -I$(top_builddir)
+DEFAULT_INCLUDES = -I. -I$(top_builddir)@am__isrc@
 depcomp = $(SHELL) $(top_srcdir)/depcomp
 am__depfiles_maybe = depfiles
 COMPILE = $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) \
@@ -66,8 +63,6 @@ ETAGS = etags
 CTAGS = ctags
 DISTFILES = $(DIST_COMMON) $(DIST_SOURCES) $(TEXINFOS) $(EXTRA_DIST)
 ACLOCAL = @ACLOCAL@
-AMDEP_FALSE = @AMDEP_FALSE@
-AMDEP_TRUE = @AMDEP_TRUE@
 AMTAR = @AMTAR@
 AUTOCONF = @AUTOCONF@
 AUTOHEADER = @AUTOHEADER@
@@ -91,6 +86,7 @@ GTKGLEXT_CFLAGS = @GTKGLEXT_CFLAGS@
 GTKGLEXT_LIBS = @GTKGLEXT_LIBS@
 GTK_CFLAGS = @GTK_CFLAGS@
 GTK_LIBS = @GTK_LIBS@
+INSTALL = @INSTALL@
 INSTALL_DATA = @INSTALL_DATA@
 INSTALL_PROGRAM = @INSTALL_PROGRAM@
 INSTALL_SCRIPT = @INSTALL_SCRIPT@
@@ -101,6 +97,7 @@ LIBS = @LIBS@
 LN_S = @LN_S@
 LTLIBOBJS = @LTLIBOBJS@
 MAKEINFO = @MAKEINFO@
+MKDIR_P = @MKDIR_P@
 OBJEXT = @OBJEXT@
 PACKAGE = @PACKAGE@
 PACKAGE_BUGREPORT = @PACKAGE_BUGREPORT@
@@ -114,9 +111,11 @@ SET_MAKE = @SET_MAKE@
 SHELL = @SHELL@
 STRIP = @STRIP@
 VERSION = @VERSION@
+abs_builddir = @abs_builddir@
+abs_srcdir = @abs_srcdir@
+abs_top_builddir = @abs_top_builddir@
+abs_top_srcdir = @abs_top_srcdir@
 ac_ct_CC = @ac_ct_CC@
-am__fastdepCC_FALSE = @am__fastdepCC_FALSE@
-am__fastdepCC_TRUE = @am__fastdepCC_TRUE@
 am__include = @am__include@
 am__leading_dot = @am__leading_dot@
 am__quote = @am__quote@
@@ -124,6 +123,7 @@ am__tar = @am__tar@
 am__untar = @am__untar@
 bindir = @bindir@
 build_alias = @build_alias@
+builddir = @builddir@
 datadir = @datadir@
 datarootdir = @datarootdir@
 docdir = @docdir@
@@ -147,9 +147,12 @@ program_transform_name = @program_transform_name@
 psdir = @psdir@
 sbindir = @sbindir@
 sharedstatedir = @sharedstatedir@
+srcdir = @srcdir@
 sysconfdir = @sysconfdir@
 target_alias = @target_alias@
-dtr_SOURCES = main.c displaywindow.c trackball.c reflections.c readpng.c mrc.c imagedisplay.c utils.c itrans.c qdrp.c ipr.c
+top_builddir = @top_builddir@
+top_srcdir = @top_srcdir@
+dtr_SOURCES = main.c displaywindow.c trackball.c reflections.c readpng.c mrc.c imagedisplay.c utils.c itrans.c qdrp.c ipr.c peakdetect.c
 dtr_LDADD = @LIBS@ @GTK_LIBS@ -lm @GTKGLEXT_LIBS@ -lgsl -lgslcblas
 AM_CFLAGS = -Wall -g @CFLAGS@ @GTK_CFLAGS@ @GTKGLEXT_CFLAGS@
 AM_CPPFLAGS = -DDATADIR=\""$(datadir)"\"
@@ -188,7 +191,7 @@ $(ACLOCAL_M4):  $(am__aclocal_m4_deps)
 	cd $(top_builddir) && $(MAKE) $(AM_MAKEFLAGS) am--refresh
 install-binPROGRAMS: $(bin_PROGRAMS)
 	@$(NORMAL_INSTALL)
-	test -z "$(bindir)" || $(mkdir_p) "$(DESTDIR)$(bindir)"
+	test -z "$(bindir)" || $(MKDIR_P) "$(DESTDIR)$(bindir)"
 	@list='$(bin_PROGRAMS)'; for p in $$list; do \
 	  p1=`echo $$p|sed 's/$(EXEEXT)$$//'`; \
 	  if test -f $$p \
@@ -211,7 +214,7 @@ clean-binPROGRAMS:
 	-test -z "$(bin_PROGRAMS)" || rm -f $(bin_PROGRAMS)
 dtr$(EXEEXT): $(dtr_OBJECTS) $(dtr_DEPENDENCIES) 
 	@rm -f dtr$(EXEEXT)
-	$(LINK) $(dtr_LDFLAGS) $(dtr_OBJECTS) $(dtr_LDADD) $(LIBS)
+	$(LINK) $(dtr_OBJECTS) $(dtr_LDADD) $(LIBS)
 
 mostlyclean-compile:
 	-rm -f *.$(OBJEXT)
@@ -225,6 +228,7 @@ distclean-compile:
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/itrans.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/main.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/mrc.Po@am__quote@
+@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/peakdetect.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/qdrp.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/readpng.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/reflections.Po@am__quote@
@@ -232,19 +236,18 @@ distclean-compile:
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/utils.Po@am__quote@
 
 .c.o:
-@am__fastdepCC_TRUE@	if $(COMPILE) -MT $@ -MD -MP -MF "$(DEPDIR)/$*.Tpo" -c -o $@ $<; \
-@am__fastdepCC_TRUE@	then mv -f "$(DEPDIR)/$*.Tpo" "$(DEPDIR)/$*.Po"; else rm -f "$(DEPDIR)/$*.Tpo"; exit 1; fi
+@am__fastdepCC_TRUE@	$(COMPILE) -MT $@ -MD -MP -MF $(DEPDIR)/$*.Tpo -c -o $@ $<
+@am__fastdepCC_TRUE@	mv -f $(DEPDIR)/$*.Tpo $(DEPDIR)/$*.Po
 @AMDEP_TRUE@@am__fastdepCC_FALSE@	source='$<' object='$@' libtool=no @AMDEPBACKSLASH@
 @AMDEP_TRUE@@am__fastdepCC_FALSE@	DEPDIR=$(DEPDIR) $(CCDEPMODE) $(depcomp) @AMDEPBACKSLASH@
 @am__fastdepCC_FALSE@	$(COMPILE) -c $<
 
 .c.obj:
-@am__fastdepCC_TRUE@	if $(COMPILE) -MT $@ -MD -MP -MF "$(DEPDIR)/$*.Tpo" -c -o $@ `$(CYGPATH_W) '$<'`; \
-@am__fastdepCC_TRUE@	then mv -f "$(DEPDIR)/$*.Tpo" "$(DEPDIR)/$*.Po"; else rm -f "$(DEPDIR)/$*.Tpo"; exit 1; fi
+@am__fastdepCC_TRUE@	$(COMPILE) -MT $@ -MD -MP -MF $(DEPDIR)/$*.Tpo -c -o $@ `$(CYGPATH_W) '$<'`
+@am__fastdepCC_TRUE@	mv -f $(DEPDIR)/$*.Tpo $(DEPDIR)/$*.Po
 @AMDEP_TRUE@@am__fastdepCC_FALSE@	source='$<' object='$@' libtool=no @AMDEPBACKSLASH@
 @AMDEP_TRUE@@am__fastdepCC_FALSE@	DEPDIR=$(DEPDIR) $(CCDEPMODE) $(depcomp) @AMDEPBACKSLASH@
 @am__fastdepCC_FALSE@	$(COMPILE) -c `$(CYGPATH_W) '$<'`
-uninstall-info-am:
 
 ID: $(HEADERS) $(SOURCES) $(LISP) $(TAGS_FILES)
 	list='$(SOURCES) $(HEADERS) $(LISP) $(TAGS_FILES)'; \
@@ -295,22 +298,21 @@ distclean-tags:
 	-rm -f TAGS ID GTAGS GRTAGS GSYMS GPATH tags
 
 distdir: $(DISTFILES)
-	@srcdirstrip=`echo "$(srcdir)" | sed 's|.|.|g'`; \
-	topsrcdirstrip=`echo "$(top_srcdir)" | sed 's|.|.|g'`; \
-	list='$(DISTFILES)'; for file in $$list; do \
-	  case $$file in \
-	    $(srcdir)/*) file=`echo "$$file" | sed "s|^$$srcdirstrip/||"`;; \
-	    $(top_srcdir)/*) file=`echo "$$file" | sed "s|^$$topsrcdirstrip/|$(top_builddir)/|"`;; \
-	  esac; \
+	@srcdirstrip=`echo "$(srcdir)" | sed 's/[].[^$$\\*]/\\\\&/g'`; \
+	topsrcdirstrip=`echo "$(top_srcdir)" | sed 's/[].[^$$\\*]/\\\\&/g'`; \
+	list='$(DISTFILES)'; \
+	  dist_files=`for file in $$list; do echo $$file; done | \
+	  sed -e "s|^$$srcdirstrip/||;t" \
+	      -e "s|^$$topsrcdirstrip/|$(top_builddir)/|;t"`; \
+	case $$dist_files in \
+	  */*) $(MKDIR_P) `echo "$$dist_files" | \
+			   sed '/\//!d;s|^|$(distdir)/|;s,/[^/]*$$,,' | \
+			   sort -u` ;; \
+	esac; \
+	for file in $$dist_files; do \
 	  if test -f $$file || test -d $$file; then d=.; else d=$(srcdir); fi; \
-	  dir=`echo "$$file" | sed -e 's,/[^/]*$$,,'`; \
-	  if test "$$dir" != "$$file" && test "$$dir" != "."; then \
-	    dir="/$$dir"; \
-	    $(mkdir_p) "$(distdir)$$dir"; \
-	  else \
-	    dir=''; \
-	  fi; \
 	  if test -d $$d/$$file; then \
+	    dir=`echo "/$$file" | sed -e 's,/[^/]*$$,,'`; \
 	    if test -d $(srcdir)/$$file && test $$d != $(srcdir); then \
 	      cp -pR $(srcdir)/$$file $(distdir)$$dir || exit 1; \
 	    fi; \
@@ -326,7 +328,7 @@ check: check-am
 all-am: Makefile $(PROGRAMS)
 installdirs:
 	for dir in "$(DESTDIR)$(bindir)"; do \
-	  test -z "$$dir" || $(mkdir_p) "$$dir"; \
+	  test -z "$$dir" || $(MKDIR_P) "$$dir"; \
 	done
 install: install-am
 install-exec: install-exec-am
@@ -374,12 +376,20 @@ info-am:
 
 install-data-am:
 
+install-dvi: install-dvi-am
+
 install-exec-am: install-binPROGRAMS
 
+install-html: install-html-am
+
 install-info: install-info-am
 
 install-man:
 
+install-pdf: install-pdf-am
+
+install-ps: install-ps-am
+
 installcheck-am:
 
 maintainer-clean: maintainer-clean-am
@@ -399,18 +409,22 @@ ps: ps-am
 
 ps-am:
 
-uninstall-am: uninstall-binPROGRAMS uninstall-info-am
+uninstall-am: uninstall-binPROGRAMS
+
+.MAKE: install-am install-strip
 
 .PHONY: CTAGS GTAGS all all-am check check-am clean clean-binPROGRAMS \
 	clean-generic ctags distclean distclean-compile \
 	distclean-generic distclean-tags distdir dvi dvi-am html \
 	html-am info info-am install install-am install-binPROGRAMS \
-	install-data install-data-am install-exec install-exec-am \
-	install-info install-info-am install-man install-strip \
+	install-data install-data-am install-dvi install-dvi-am \
+	install-exec install-exec-am install-html install-html-am \
+	install-info install-info-am install-man install-pdf \
+	install-pdf-am install-ps install-ps-am install-strip \
 	installcheck installcheck-am installdirs maintainer-clean \
 	maintainer-clean-generic mostlyclean mostlyclean-compile \
 	mostlyclean-generic pdf pdf-am ps ps-am tags uninstall \
-	uninstall-am uninstall-binPROGRAMS uninstall-info-am
+	uninstall-am uninstall-binPROGRAMS
 
 # Tell versions [3.59,3.63) of GNU make to not export all variables.
 # Otherwise a system limit (for SysV at least) may be exceeded.
diff --git a/src/control.h b/src/control.h
index 0493757..d557f41 100644
--- a/src/control.h
+++ b/src/control.h
@@ -32,7 +32,8 @@ typedef enum {
 	PEAKSEARCH_THRESHOLD,
 	PEAKSEARCH_ADAPTIVE_THRESHOLD,
 	PEAKSEARCH_LSQ,
-	PEAKSEARCH_ZAEFFERER
+	PEAKSEARCH_ZAEFFERER,
+	PEAKSEARCH_ITERATIVE
 } PeakSearchMode;
 
 typedef enum {
diff --git a/src/itrans.c b/src/itrans.c
index 0480f3f..1b66915 100644
--- a/src/itrans.c
+++ b/src/itrans.c
@@ -17,6 +17,7 @@
 #include "control.h"
 #include "imagedisplay.h"
 #include "utils.h"
+#include "peakdetect.h"
 
 #define MAX_LSQ_ITERATIONS 100
 #define PEAK_WINDOW_SIZE 20
@@ -499,6 +500,30 @@ static unsigned int itrans_peaksearch_lsq(int16_t *image, ControlContext *ctx, d
 	return n_reflections;
 }
 
+static unsigned int itrans_peaksearch_iterative(int16_t *image, ControlContext *ctx, double tilt_degrees, ImageDisplay *imagedisplay) {
+	printf("Starting itrans_peaksearch_iterative\n");
+	int n_reflections;	
+	gsl_matrix *m;
+	gsl_matrix *p;
+	printf("Calling createImageMatrix\n");
+	m = createImageMatrix(ctx,image);
+	printf("Calling iterate\n");
+	p = iterate(m, &n_reflections);
+	int i;
+	double px,py;
+	printf("Found %d reflections\n",n_reflections);
+	for (i=0;i<n_reflections;i++) {
+		px = gsl_matrix_get(p,0,i);
+		py = gsl_matrix_get(p,1,i);
+		printf("Reflection %d (%lf,%lf)\n",i,px,py);
+		reflection_add_from_dp(ctx, (px-ctx->x_centre), (py-ctx->y_centre), tilt_degrees, 1.);
+		if (ctx->first_image) imagedisplay_mark_point(imagedisplay, (unsigned int)px, (unsigned int)py);
+	}
+	gsl_matrix_free(m);
+	gsl_matrix_free(p);
+	return (unsigned)n_reflections;
+}
+
 void itrans_process_image(int16_t *image, ControlContext *ctx, double tilt_degrees) {
 	
 	unsigned int n_reflections;
@@ -514,6 +539,7 @@ void itrans_process_image(int16_t *image, ControlContext *ctx, double tilt_degre
 		case PEAKSEARCH_ADAPTIVE_THRESHOLD : n_reflections = itrans_peaksearch_adaptive_threshold(image, ctx, tilt_degrees, imagedisplay); break;
 		case PEAKSEARCH_LSQ : n_reflections = itrans_peaksearch_lsq(image, ctx, tilt_degrees, imagedisplay); break;
 		case PEAKSEARCH_ZAEFFERER : n_reflections = itrans_peaksearch_zaefferer(image, ctx, tilt_degrees, imagedisplay); break;
+		case PEAKSEARCH_ITERATIVE : n_reflections = itrans_peaksearch_iterative(image, ctx, tilt_degrees, imagedisplay); break;
 		default: n_reflections = 0;
 	}
 	
diff --git a/src/main.c b/src/main.c
index 272c9a8..a89fdf9 100644
--- a/src/main.c
+++ b/src/main.c
@@ -45,6 +45,7 @@ static gint main_method_window_response(GtkWidget *method_window, gint response,
 			case 1 : ctx->psmode = PEAKSEARCH_ADAPTIVE_THRESHOLD; break;
 			case 2 : ctx->psmode = PEAKSEARCH_LSQ; break;
 			case 3 : ctx->psmode = PEAKSEARCH_ZAEFFERER; break;
+			case 4 : ctx->psmode = PEAKSEARCH_ITERATIVE; break;
 			default: abort();
 		}
 		
@@ -113,6 +114,7 @@ void main_method_dialog_open(ControlContext *ctx) {
 	gtk_combo_box_append_text(GTK_COMBO_BOX(ctx->combo_peaksearch), "Adaptive Thresholding");
 	gtk_combo_box_append_text(GTK_COMBO_BOX(ctx->combo_peaksearch), "Least-Squares Fit");
 	gtk_combo_box_append_text(GTK_COMBO_BOX(ctx->combo_peaksearch), "Zaefferer Gradient Search");
+	gtk_combo_box_append_text(GTK_COMBO_BOX(ctx->combo_peaksearch), "Iterative Statistical");
 	gtk_combo_box_set_active(GTK_COMBO_BOX(ctx->combo_peaksearch), 3);
 	gtk_table_attach_defaults(GTK_TABLE(table), ctx->combo_peaksearch, 2, 3, 2, 3);
 
diff --git a/src/peakdetect.c b/src/peakdetect.c
new file mode 100644
index 0000000..666e23f
--- /dev/null
+++ b/src/peakdetect.c
@@ -0,0 +1,449 @@
+/*
+ * peakdetect.c
+ *
+ * Peakdetection routines, utilities and automation
+ *
+ * (c) 2007 Gordon Ball <gfb21@cam.ac.uk>
+ *  dtr - Diffraction Tomography Reconstruction
+ *
+ */
+
+#include <math.h>
+#include "reflections.h"
+#include "control.h"
+#include "imagedisplay.h"
+#include "utils.h"
+#include "peakdetect.h"
+
+/*
+ * Renormalise a gsl_matrix to 0->1
+ * Re-written to use gsl_matrix native functions
+ */
+void renormalise(gsl_matrix *m) {
+	double max,min;
+	gsl_matrix_minmax(m,&min,&max);
+	printf("Renormalising min=%lf max=%lf\n",min,max);
+	gsl_matrix_add_constant(m,0.-min);
+	gsl_matrix_scale(m,1./(max-min));
+}
+
+/*
+ * Create a gsl_matrix for performing image operations on
+ * from a raw image and control context
+ * Converting to gsl_matrix because many of the required operations
+ * can be done as matrices to save effort
+ * Renormalises matrix to 0->1
+ */
+gsl_matrix* createImageMatrix(ControlContext *ctx, int16_t *image) {
+	gsl_matrix *raw;
+	raw = gsl_matrix_alloc(ctx->width,ctx->height);
+	
+	int i,j;
+	for ( i=0; i<raw->size1;i++) {
+		for (j=0; j<raw->size2;j++) {
+			gsl_matrix_set(raw,i,j,(double)image[i+j*ctx->width]);
+		}
+	}
+	printf("Created %dx%d matrix\n",ctx->width,ctx->height);
+	renormalise(raw);
+	return raw;
+}
+
+
+
+/*
+ * Find and return the mean value of the matrix
+ */
+double image_mean(gsl_matrix *m) {
+	double mean=0.;
+	int i,j;
+	for (i=0;i<m->size1;i++) {
+		for (j=0;j<m->size2;j++) {
+			mean += gsl_matrix_get(m,i,j);
+		}
+	}
+	return mean / (m->size1 * m->size2);
+}
+
+/*
+ * Return the standard deviation, sigma, of the matrix values
+ * \sqrt(\sum((x-mean)^2)/n)
+ */
+double image_sigma(gsl_matrix *m, double mean) {
+	double diff2=0;
+	int i,j;
+	for (i=0;i<m->size1;i++) {
+		for (j=0;j<m->size2;j++) {
+			diff2 += (gsl_matrix_get(m,i,j)-mean)*(gsl_matrix_get(m,i,j)-mean);
+		}
+	}
+	return sqrt(diff2/(m->size1 * m->size2));
+}
+
+/*
+ * Filter all pixels with value < mean + k*sigma to 0
+ * Set all matching pixels to 1
+ */
+void sigma_filter(gsl_matrix *m, double k) {
+	double mean,sigma;
+	int i,j;
+	mean = image_mean(m);
+	sigma = image_sigma(m,mean);
+	for (i=0;i<m->size1;i++) {
+		for (j=0;j<m->size2;j++) {
+			if (gsl_matrix_get(m,i,j) >= mean+k*sigma) {
+				gsl_matrix_set(m,i,j,1.);
+			} else {
+				gsl_matrix_set(m,i,j,0.);
+			}
+		}
+	}
+}
+
+/*
+ * Calculate the mean within a circle centred (x,y) of radius r
+ *
+ * TODO: Use a mask instead of calculating valid points
+ */
+double circle_mean(gsl_matrix *m, int x, int y, int r, gsl_matrix *mask) {
+	double mean=0.;
+	int i,j,n=0;
+	for (i=x-r;i<=x+r;i++) {
+		for (j=y-r;j<=y+r;j++) {
+			//printf("cm: ij=(%d,%d) mask=(%d,%d)\n",i,j,i-x+r,j-y+r);
+			if (gsl_matrix_get(mask,i-x+r,j-y+r)>0.) {
+				mean += gsl_matrix_get(m,i,j);
+				//printf("cm: (%d,%d) mean=%lf val=%lf\n",i,j,mean,gsl_matrix_get(m,i,j));
+				n++;
+			}
+		}
+	}
+	//printf("cm: (%d,%d) summean=%lf n=%d\n",x,y,mean,n);
+	return mean/n;
+}
+
+/*
+ * Calculate sigma within a circle centred (x,y) of radius r
+ *
+ * TODO: Use a mask instead of calculating valid points
+ */
+double circle_sigma(gsl_matrix *m, int x, int y, int r, gsl_matrix *mask, double mean) {
+	double diff2=0.;
+	int i,j,n=0;
+	for (i=x-r;i<=x+r;i++) {
+		for (j=y-r;j<=y+r;j++) {
+			if (gsl_matrix_get(mask,i-x+r,j-y+r)>0) {
+				diff2 += (gsl_matrix_get(m,i,j)-mean)*(gsl_matrix_get(m,i,j)-mean);
+				n++;
+			}
+		}
+	}
+	return sqrt(diff2/n);
+}
+
+/*
+ * Calculate a circular mask to save recalculating it every time
+ */
+gsl_matrix* circle_mask(int r) {
+	gsl_matrix *m;
+	m = gsl_matrix_calloc(2*r+1,2*r+1);
+	int i,j;
+	for (i=0;i<2*r+1;i++) {
+		for (j=0;j<2*r+1;j++) {
+			if (sqrt((r-i)*(r-i)+(r-j)*(r-j))<=(double)r) {
+				 gsl_matrix_set(m,i,j,1.);
+			}
+			
+		}
+	}
+	return m;
+}
+
+/*
+ * Variation on the above filter where instead of using
+ * sigma and mean for the whole image, it is found for a local
+ * circle of radius r pixels.
+ * The central point is also calculated as an approximately gaussian
+ * average over local pixels rather than just a single pixel.
+ * This takes a long time - 20-30 seconds for a 1024^2 image and r=10,
+ * obviously large values of r will make it even slower.
+ * The appropriate r value needs to be considered carefully - it needs to
+ * be somewhat smaller than the average inter-reflection seperation.
+ * 10 pixels worked well for the sapphire.mrc images, but this might
+ * not be the case for other images. Images for which this would be very
+ * large probably need to be resampled to a lower resolution.
+ *
+ * TODO: Pass a mask to the ancilliary functions instead of having
+ * them calculate several hundred million sqrts
+ * 
+ * self-referencing problem being dealt with - output being written onto the array before the next point it computed
+ * problem carried over from the OO version where a new object was created by each operation
+ */
+
+void local_sigma_filter(gsl_matrix *m, int r, double k) {
+	//printf("lsf: starting\n");
+	double mean,sigma;
+	double local;
+	//printf("lsf: generating circle mask\n");
+	gsl_matrix *mask = circle_mask(r);
+	gsl_matrix *new;
+	new = gsl_matrix_alloc(m->size1,m->size2);
+	int i,j;
+	int interval = (m->size1-r)/20;
+	//printf("lsf: starting loop\n");
+	printf("lsf: ");
+	//for (i=r;i<m->size1-r;i++) {
+	//	for (j=r;j<m->size2-r;j++) {
+	for (i=0;i<m->size1;i++) {
+		for (j=0;j<m->size2;j++) {
+			if ((i>=r && i<m->size1-r) && (j>=r && j<m->size2-r)) {
+				//printf("lsf: evaluating (%d,%d)\n",i,j);
+				mean = circle_mean(m,i,j,r,mask);
+				//printf("lsf: mean=%lf",mean);
+				sigma = circle_sigma(m,i,j,r,mask,mean);
+				//printf(" sigma=%lf",sigma);
+				local = gsl_matrix_get(m,i,j);
+				local += gsl_matrix_get(m,i+1,j) + gsl_matrix_get(m,i-1,j) + gsl_matrix_get(m,i,j+1) + gsl_matrix_get(m,i,j-1);
+				local += .5*(gsl_matrix_get(m,i+1,j+1) + gsl_matrix_get(m,i-1,j+1) + gsl_matrix_get(m,i+1,j-1) + gsl_matrix_get(m,i-1,j-1));
+				local /= 7.;
+				//printf(" local=%lf\n",local);
+				if (local > mean+k*sigma) {
+					gsl_matrix_set(new,i,j,1.);
+				} else {
+					gsl_matrix_set(new,i,j,0.);
+				}
+			} else {
+				gsl_matrix_set(new,i,j,0.);
+			}
+		}
+		if (i % interval == 0) printf(".");
+	}
+	printf("done\n");
+	gsl_matrix_memcpy(m,new);
+	gsl_matrix_free(new);
+}
+
+
+
+/*
+ * Apply an arbitary kernel to the image - each point takes the value
+ * of the sum of the kernel convolved with the surrounding pixels.
+ * The kernel needs to be a (2n+1)^2 array (centred on (n,n)). It needs
+ * to be square and should probably be normalised.
+ * Don't do daft things like rectangular kernels or kernels larger
+ * than the image - this doesn't check such things and will cry.
+ * 
+ * Also suffers from self-reference problem
+ */
+void apply_kernel(gsl_matrix *m, gsl_matrix *kernel) {
+	int size = kernel->size1;
+	int half = (size-1)/2;
+	gsl_matrix *l;
+	gsl_matrix_view lv;
+	gsl_matrix *new;
+	new = gsl_matrix_calloc(m->size1,m->size2);
+	double val;
+	int i,j,x,y;
+	for (i=0;i<m->size1;i++) {
+		for (j=0;j<m->size2;j++) {
+			if ((i>=half && i<m->size1-half) && (j>=half && j<m->size2-half)) {
+				lv = gsl_matrix_submatrix(m,i-half,j-half,size,size);
+				l = &lv.matrix;
+				val = 0.;
+				for (x=0;x<size;x++) {
+					for (y=0;y<size;y++) {
+						val += gsl_matrix_get(l,x,y)*gsl_matrix_get(kernel,x,y);
+					}
+				}
+				//gsl_matrix_free(l);
+				gsl_matrix_set(new,i,j,val);
+			}
+		}
+	}
+	gsl_matrix_memcpy(m,new);
+	gsl_matrix_free(new);
+}
+
+/*
+ * Generate the simplist possible kernel - a flat one
+ */
+gsl_matrix* generate_flat_kernel(int half) {
+	gsl_matrix *k;
+	k = gsl_matrix_alloc(2*half+1,2*half+1);
+	gsl_matrix_set_all(k,1./((2*half+1)*(2*half+1)));
+	return k;
+}
+
+/*
+ * expands or contracts a gsl_matrix by copying the columns to a new one
+ */
+gsl_matrix* matrix_expand(gsl_matrix *m, int oldsize, int newsize) {
+	gsl_matrix *new;
+	
+	printf("me: %d->%d\n",oldsize,newsize);
+	
+	new = gsl_matrix_calloc(2,newsize);
+	//printf("me: calloc(2,%d)\n",newsize);
+	int j;
+	for (j = 0; j < oldsize; j++) {
+		if (j < newsize) {
+			//printf("me: copying col %d\n",j);
+			gsl_matrix_set(new,0,j,gsl_matrix_get(m,0,j));
+			gsl_matrix_set(new,1,j,gsl_matrix_get(m,1,j));
+		}
+	}
+	//printf("me: freeing old matrix\n");
+	gsl_matrix_free(m);
+	//printf("me: new s1=%d s2=%d\n",new->size1,new->size2);
+	return new;
+	//printf("me: m s1=%d s2=%d\n",m->size1,m->size2);
+}
+
+/*
+ * Stack-based flood-fill iteration routine 
+ * have to return a pointer to com each time because if the matrix size has to be changed then we need to keep
+ * track of the location of the resized instance
+ */
+ 
+gsl_matrix* do_ff(int i, int j, int* mask, double threshold, gsl_matrix *m, gsl_matrix *com, int *com_n, int *com_size) {
+	if (i>=0 && i<m->size1) {
+		if (j>=0 && j<m->size2) {
+			if (mask[i+j*m->size1]==0) {
+				if (gsl_matrix_get(m,i,j)>threshold) {
+					//printf("dff: found valid point (%d,%d)\n",i,j);
+					gsl_matrix_set(com,0,*com_n,(double)i);
+					gsl_matrix_set(com,1,*com_n,(double)j);
+					*com_n=*com_n+1;
+					if (*com_n == *com_size) {
+						com = matrix_expand(com,*com_size,*com_size*2);
+						*com_size *= 2;
+					}
+					mask[i+j*m->size1]=1;
+					com = do_ff(i+1,j,mask,threshold,m,com,com_n,com_size);
+					com = do_ff(i-1,j,mask,threshold,m,com,com_n,com_size);
+					com = do_ff(i,j+1,mask,threshold,m,com,com_n,com_size);
+					com = do_ff(i,j-1,mask,threshold,m,com,com_n,com_size);
+				}
+			}
+		}
+	}
+	return com;
+}
+
+/*
+ * Find points by floodfilling all pixels above threshold
+ * Implements a stack-based flood-filling method which may
+ * cause stack overflows if the blocks are extremely large -
+ * dependent on implementation (default 4MiB stack for unix?)
+ * Implements a crude variable-sized array method which hopefully works
+ * Returns a gsl_matrix with x coordinates in row 0 and y
+ * coordinates in row 1, which should be of the right length
+ * Variable count is set to the number of points found
+ */
+
+gsl_matrix* floodfill(gsl_matrix *m, double threshold, int *count) {
+	//printf("ff: starting\n");
+	int *mask;
+	mask = calloc(m->size1*m->size2,sizeof(int));
+	
+	int size=32,com_size;
+	int i,j,k,n=0;
+	int com_n;
+	gsl_matrix *p;
+	gsl_matrix *com;
+	p = gsl_matrix_calloc(2,size);
+	
+	double com_x,com_y;
+	printf("ff: starting loop\n");
+	for (i=0;i<m->size1;i++) {
+		for (j=0;j<m->size2;j++) {
+			if (gsl_matrix_get(m,i,j)>threshold) {
+				if (mask[i+j*m->size1]==0) {
+					//printf("ff: found starting point (%d,%d)\n",i,j);
+					com_size=32;
+					com_n=0;
+					com_x = com_y = 0.;					
+					com = gsl_matrix_calloc(2,com_size); //this is going to hold the points found for this location
+					//printf("ff: starting floodfill stack\n");
+					com = do_ff(i, j, mask, threshold, m, com, &com_n, &com_size);
+					//printf("ff: ended floodfill stack\n");
+					for (k=0;k<com_n;k++) {
+						com_x += gsl_matrix_get(com,0,k);
+						com_y += gsl_matrix_get(com,1,k);
+					}
+					com_x /= com_n;
+					com_y /= com_n;
+					printf("ff: point CoM (%lf,%lf)\n",com_x,com_y);
+					gsl_matrix_set(p,0,n,com_x);
+					gsl_matrix_set(p,1,n,com_y);
+					n++;
+					if (n==size) {
+						p = matrix_expand(p,size,size*2);
+						size *= 2;
+					}
+				}
+			}
+		}
+	}
+	printf("ff: ending loop, found %d\n",n);
+	*count = n;
+	printf("pcheck s1=%d s2=%d\n",p->size1,p->size2);
+	p = matrix_expand(p,size,n);
+	printf("pcheck s1=%d s2=%d\n",p->size1,p->size2);
+	return p;
+}
+
+
+
+/*
+ * implements the iteration based automatic method
+ * returns a gsl_matrix formatted as described in flood-fill
+ */
+gsl_matrix* iterate(gsl_matrix *m, int *count) {
+	printf("Iterate: starting\n");
+	int old = m->size1*m->size2;
+	int cur;
+	double k;
+	double mean,sigma;
+	gsl_matrix *p;
+	gsl_matrix *kernel;
+	
+	mean = image_mean(m);
+	sigma = image_sigma(m,mean);
+	printf("Iterate: mean=%lf sigma=%lf\n",mean,sigma);
+	
+	
+	printf("Iterate: generating kernel\n");
+	kernel = generate_flat_kernel(1);
+	printf("Iterate: performing local_sigma_filter\n");
+	local_sigma_filter(m,10,1.);
+	
+	mean = image_mean(m);
+	sigma = image_sigma(m,mean);
+	printf("Iterate: mean=%lf sigma=%lf\n",mean,sigma);
+	
+	
+	printf("Iterate: starting loop\n");
+	while (1) {
+		printf("Iterate: smoothing");
+		apply_kernel(m,kernel);
+		printf(" (1)");
+		apply_kernel(m,kernel);
+		printf(" (2)\n");
+		mean = image_mean(m);
+		sigma = image_sigma(m,mean);
+		printf("Iterate: mean=%lf sigma=%lf\n",mean,sigma);
+		k = (0.5-mean)/sigma;
+		printf("Iterate: applying sigma_filter(%lf)\n",k);
+		sigma_filter(m,k);
+		printf("Iterate: floodfilling\n");
+		p = floodfill(m,0,&cur);
+		printf("Iterate: %d points found\n",cur);
+		if (old < 1.05*cur) break;
+		old = cur;
+	}
+	gsl_matrix_free(kernel); 
+	printf("Iterate: finished\n");
+	*count = cur;
+	return p;
+}
diff --git a/src/peakdetect.h b/src/peakdetect.h
new file mode 100644
index 0000000..13837d1
--- /dev/null
+++ b/src/peakdetect.h
@@ -0,0 +1,19 @@
+/*
+ * peakdetect.h
+ */
+
+#include <gsl/gsl_matrix.h>
+
+extern gsl_matrix* createImageMatrix(ControlContext *ctx, int16_t *image);
+
+extern void sigma_filter(gsl_matrix *m, double k);
+
+extern void local_sigma_filter(gsl_matrix *m, int r, double k);
+
+extern void apply_kernel(gsl_matrix *m, gsl_matrix *kernel);
+
+extern gsl_matrix* generate_flat_kernel(int half);
+
+extern gsl_matrix* floodfill(gsl_matrix *m, double threshold, int* count);
+
+extern gsl_matrix* iterate(gsl_matrix *m, int* count);