Move fine peak search to new file

4 files changed, 180 insertions, 123 deletions

diff --git a/src/Makefile.am b/src/Makefile.am
index 9010373c..948d2fbf 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -17,7 +17,7 @@ process_hkl_SOURCES = process_hkl.c sfac.c statistics.c cell.c utils.c \
 process_hkl_LDADD = @LIBS@
 
 indexamajig_SOURCES = indexamajig.c hdf5-file.c utils.c dirax.c cell.c image.c \
-                      intensities.c ewald.c
+                      intensities.c ewald.c peaks.c
 indexamajig_LDADD = @LIBS@
 
 if HAVE_GTK
diff --git a/src/dirax.c b/src/dirax.c
index 070e714d..9930ca03 100644
--- a/src/dirax.c
+++ b/src/dirax.c
@@ -34,6 +34,7 @@
 #include "dirax.h"
 #include "utils.h"
 #include "sfac.h"
+#include "peaks.h"
 
 
 typedef enum {
@@ -377,13 +378,15 @@ static int map_position(struct image *image, double x, double y,
 }
 
 
-#define PEAK_WINDOW_SIZE (10)
-
-static void search_peaks(struct image *image, int dump_peaks)
+void index_pattern(struct image *image, int no_index, int dump_peaks)
 {
+	unsigned int opts;
+	int saved_stderr;
 	FILE *fh;
-	int x, y, width, height;
-	int16_t *data;
+	int i;
+
+	/* Do peak search and splurge out 'xfel.drx' */
+	search_peaks(image, dump_peaks);
 
 	fh = fopen("xfel.drx", "w");
 	if ( !fh ) {
@@ -392,128 +395,20 @@ static void search_peaks(struct image *image, int dump_peaks)
 	}
 	fprintf(fh, "%f\n", 0.5);  /* Lie about the wavelength.  */
 
-	data = image->data;
-	width = image->width;
-	height = image->height;
-
-	if ( image->features != NULL ) {
-		image_feature_list_free(image->features);
-	}
-	image->features = image_feature_list_new();
-
-	for ( x=1; x<image->width-1; x++ ) {
-	for ( y=1; y<image->height-1; y++ ) {
-
-		double dx1, dx2, dy1, dy2;
-		double dxs, dys;
-		double grad;
-		int mask_x, mask_y;
-		int sx, sy;
-		double max;
-		unsigned int did_something = 1;
-
-		/* Overall threshold */
-		if ( data[x+width*y] < 800 ) continue;
-
-		/* Ignore streak */
-		if ( abs(x-image->x_centre) < 15 ) continue;
+	for ( i=0; i<image_feature_count(image->features); i++ ) {
 
-		/* Get gradients */
-		dx1 = data[x+width*y] - data[(x+1)+width*y];
-		dx2 = data[(x-1)+width*y] - data[x+width*y];
-		dy1 = data[x+width*y] - data[(x+1)+width*(y+1)];
-		dy2 = data[x+width*(y-1)] - data[x+width*y];
+		struct imagefeature *f;
+		double rx = 0.0;
+		double ry = 0.0;
+		double rz = 0.0;
 
-		/* Average gradient measurements from both sides */
-		dxs = ((dx1*dx1) + (dx2*dx2)) / 2;
-		dys = ((dy1*dy1) + (dy2*dy2)) / 2;
-
-		/* Calculate overall gradient */
-		grad = dxs + dys;
-
-		if ( grad < 2000000 ) continue;
-
-		mask_x = x;
-		mask_y = y;
-
-		while ( (did_something)
-		     && (distance(mask_x, mask_y, x, y)<50) ) {
-
-			max = data[mask_x+width*mask_y];
-			did_something = 0;
-
-			for ( sy=biggest(mask_y-PEAK_WINDOW_SIZE/2, 0);
-			      sy<smallest(mask_y+PEAK_WINDOW_SIZE/2, height);
-			      sy++ ) {
-			for ( sx=biggest(mask_x-PEAK_WINDOW_SIZE/2, 0);
-			      sx<smallest(mask_x+PEAK_WINDOW_SIZE/2, width);
-			      sx++ ) {
-
-				if ( data[sx+width*sy] > max ) {
-					max = data[sx+width*sy];
-					mask_x = sx;
-					mask_y = sy;
-					did_something = 1;
-				}
-
-			}
-			}
-
-		}
+		f = image_get_feature(image->features, i);
+		map_position(image, f->x, f->y, &rx, &ry, &rz);
+		fprintf(fh, "%10f %10f %10f %8f\n",
+		        rx/1e10, ry/1e10, rz/1e10, 1.0);
 
-		if ( !did_something ) {
-
-			double d;
-			int idx;
-
-			assert(mask_x<image->width);
-			assert(mask_y<image->height);
-			assert(mask_x>=0);
-			assert(mask_y>=0);
-
-			/* Too far from foot point? */
-			if ( distance(mask_x, mask_y, x, y) > 50.0 )  continue;
-
-			/* Check for a feature at exactly the
-			 * same coordinates */
-			image_feature_closest(image->features, mask_x, mask_y,
-			                      &d, &idx);
-
-			if ( d > 1.0 ) {
-
-				double rx = 0.0;
-				double ry = 0.0;
-				double rz = 0.0;
-
-				/* Map and record reflection */
-				if ( dump_peaks ) {
-					printf("%i %i\n", x, y);
-				}
-
-				image_add_feature(image->features,
-				                  mask_x, mask_y, image, 1.0);
-				map_position(image, x, y, &rx, &ry, &rz);
-				fprintf(fh, "%10f %10f %10f %8f\n",
-				        rx/1e10, ry/1e10, rz/1e10, 1.0);
-			}
-
-		}
-
-
-	}
 	}
-
 	fclose(fh);
-}
-
-
-void index_pattern(struct image *image, int no_index, int dump_peaks)
-{
-	unsigned int opts;
-	int saved_stderr;
-
-	/* Do peak search and splurge out 'xfel.drx' */
-	search_peaks(image, dump_peaks);
 
 	if ( no_index ) return;
 
diff --git a/src/peaks.c b/src/peaks.c
new file mode 100644
index 00000000..afc05401
--- /dev/null
+++ b/src/peaks.c
@@ -0,0 +1,138 @@
+/*
+ * peaks.c
+ *
+ * Peak search and other image analysis
+ *
+ * (c) 2006-2009 Thomas White <taw@physics.org>
+ *
+ * Part of CrystFEL - crystallography with a FEL
+ *
+ */
+
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+#include <string.h>
+#include <assert.h>
+
+#include "image.h"
+#include "utils.h"
+
+
+#define PEAK_WINDOW_SIZE (10)
+
+void search_peaks(struct image *image, int dump_peaks)
+{
+	int x, y, width, height;
+	int16_t *data;
+
+	data = image->data;
+	width = image->width;
+	height = image->height;
+
+	if ( image->features != NULL ) {
+		image_feature_list_free(image->features);
+	}
+	image->features = image_feature_list_new();
+
+	for ( x=1; x<image->width-1; x++ ) {
+	for ( y=1; y<image->height-1; y++ ) {
+
+		double dx1, dx2, dy1, dy2;
+		double dxs, dys;
+		double grad;
+		int mask_x, mask_y;
+		int sx, sy;
+		double max;
+		unsigned int did_something = 1;
+
+		/* Overall threshold */
+		if ( data[x+width*y] < 800 ) continue;
+
+		/* Ignore streak */
+		if ( abs(x-image->x_centre) < 15 ) continue;
+
+		/* Get gradients */
+		dx1 = data[x+width*y] - data[(x+1)+width*y];
+		dx2 = data[(x-1)+width*y] - data[x+width*y];
+		dy1 = data[x+width*y] - data[(x+1)+width*(y+1)];
+		dy2 = data[x+width*(y-1)] - data[x+width*y];
+
+		/* Average gradient measurements from both sides */
+		dxs = ((dx1*dx1) + (dx2*dx2)) / 2;
+		dys = ((dy1*dy1) + (dy2*dy2)) / 2;
+
+		/* Calculate overall gradient */
+		grad = dxs + dys;
+
+		if ( grad < 2000000 ) continue;
+
+		mask_x = x;
+		mask_y = y;
+
+		while ( (did_something)
+		     && (distance(mask_x, mask_y, x, y)<50) ) {
+
+			max = data[mask_x+width*mask_y];
+			did_something = 0;
+
+			for ( sy=biggest(mask_y-PEAK_WINDOW_SIZE/2, 0);
+			      sy<smallest(mask_y+PEAK_WINDOW_SIZE/2, height);
+			      sy++ ) {
+			for ( sx=biggest(mask_x-PEAK_WINDOW_SIZE/2, 0);
+			      sx<smallest(mask_x+PEAK_WINDOW_SIZE/2, width);
+			      sx++ ) {
+
+				if ( data[sx+width*sy] > max ) {
+					max = data[sx+width*sy];
+					mask_x = sx;
+					mask_y = sy;
+					did_something = 1;
+				}
+
+			}
+			}
+
+		}
+
+		if ( !did_something ) {
+
+			double d;
+			int idx;
+
+			assert(mask_x<image->width);
+			assert(mask_y<image->height);
+			assert(mask_x>=0);
+			assert(mask_y>=0);
+
+			/* Too far from foot point? */
+			if ( distance(mask_x, mask_y, x, y) > 50.0 )  continue;
+
+			/* Check for a feature at exactly the
+			 * same coordinates */
+			image_feature_closest(image->features, mask_x, mask_y,
+			                      &d, &idx);
+
+			if ( d > 1.0 ) {
+
+				/* Map and record reflection */
+				if ( dump_peaks ) {
+					printf("%i %i\n", x, y);
+				}
+
+				image_add_feature(image->features,
+				                  mask_x, mask_y, image, 1.0);
+
+			}
+
+		}
+
+
+	}
+	}
+}
diff --git a/src/peaks.h b/src/peaks.h
new file mode 100644
index 00000000..d7bbca4d
--- /dev/null
+++ b/src/peaks.h
@@ -0,0 +1,24 @@
+/*
+ * peaks.h
+ *
+ * Peak search and other image analysis
+ *
+ * (c) 2006-2009 Thomas White <taw@physics.org>
+ *
+ * Part of CrystFEL - crystallography with a FEL
+ *
+ */
+
+
+#ifndef PEAKS_H
+#define PEAKS_H
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+
+extern void search_peaks(struct image *image, int dump_peaks);
+
+
+#endif	/* PEAKS_H */