Move processing pipeline to separate file

4 files changed, 329 insertions, 250 deletions

diff --git a/src/im-sandbox.c b/src/im-sandbox.c
index 504c52c4..290f811f 100644
--- a/src/im-sandbox.c
+++ b/src/im-sandbox.c
@@ -41,8 +41,6 @@
 #include <string.h>
 #include <unistd.h>
 #include <getopt.h>
-#include <hdf5.h>
-#include <gsl/gsl_errno.h>
 #include <pthread.h>
 #include <sys/wait.h>
 #include <fcntl.h>
@@ -55,36 +53,14 @@
 #include <sys/time.h>
 #endif
 
-#include "utils.h"
-#include "hdf5-file.h"
-#include "index.h"
-#include "peaks.h"
-#include "detector.h"
-#include "filters.h"
-#include "thread-pool.h"
-#include "beam-parameters.h"
-#include "geometry.h"
-#include "stream.h"
-#include "reflist-utils.h"
-
 #include "im-sandbox.h"
+#include "process_image.h"
 
 
 /* Write statistics at APPROXIMATELY this interval */
 #define STATS_EVERY_N_SECONDS (5)
 
 
-/* Information about the indexing process for one pattern */
-struct pattern_args
-{
-	/* "Input" */
-	char *filename;
-
-	/* "Output" */
-	int n_crystals;
-};
-
-
 struct sb_reader
 {
 	pthread_mutex_t lock;
@@ -194,190 +170,6 @@ static char *get_pattern(FILE *fh, char **use_this_one_instead,
 }
 
 
-static void process_image(const struct index_args *iargs,
-                          struct pattern_args *pargs, Stream *st,
-                          int cookie)
-{
-	float *data_for_measurement;
-	size_t data_size;
-	int check;
-	struct hdfile *hdfile;
-	struct image image;
-	int i;
-	char filename[1024];
-
-	/* Prefix to jump out of temporary folder */
-	if ( pargs->filename[0] != '/' ) {
-		snprintf(filename, 1023, "../../%s", pargs->filename);
-	} else {
-		snprintf(filename, 1023, "%s", pargs->filename);
-	}
-
-	image.features = NULL;
-	image.data = NULL;
-	image.flags = NULL;
-	image.copyme = iargs->copyme;
-	image.id = cookie;
-	image.filename = pargs->filename;  /* Relative to top level */
-	image.beam = iargs->beam;
-	image.det = iargs->det;
-	image.crystals = NULL;
-	image.n_crystals = 0;
-
-	hdfile = hdfile_open(filename);  /* Relative to temporary folder */
-	if ( hdfile == NULL ) return;
-
-	if ( iargs->element != NULL ) {
-
-		int r;
-		r = hdfile_set_image(hdfile, iargs->element);
-		if ( r ) {
-			ERROR("Couldn't select path '%s'\n", iargs->element);
-			hdfile_close(hdfile);
-			return;
-		}
-
-	} else {
-
-		int r;
-		r = hdfile_set_first_image(hdfile, "/");
-		if ( r ) {
-			ERROR("Couldn't select first path\n");
-			hdfile_close(hdfile);
-			return;
-		}
-
-	}
-
-	check = hdf5_read(hdfile, &image, 1);
-	if ( check ) {
-		hdfile_close(hdfile);
-		return;
-	}
-
-	if ( (image.width != image.det->max_fs + 1 )
-	  || (image.height != image.det->max_ss + 1))
-	{
-		ERROR("Image size doesn't match geometry size"
-			" - rejecting image.\n");
-		ERROR("Image size: %i,%i.  Geometry size: %i,%i\n",
-		      image.width, image.height,
-		      image.det->max_fs + 1, image.det->max_ss + 1);
-		hdfile_close(hdfile);
-		return;
-	}
-
-	fill_in_values(image.det, hdfile);
-	fill_in_beam_parameters(image.beam, hdfile);
-
-	image.lambda = ph_en_to_lambda(eV_to_J(image.beam->photon_energy));
-
-	if ( (image.beam->photon_energy < 0.0) || (image.lambda > 1000) ) {
-		/* Error message covers a silly value in the beam file or in
-		 * the HDF5 file. */
-		ERROR("Nonsensical wavelength (%e m or %e eV) value for %s.\n",
-		      image.lambda, image.beam->photon_energy, image.filename);
-		hdfile_close(hdfile);
-		return;
-	}
-
-	if ( iargs->cmfilter ) filter_cm(&image);
-
-	/* Take snapshot of image after CM subtraction but before
-	 * the aggressive noise filter. */
-	data_size = image.width * image.height * sizeof(float);
-	data_for_measurement = malloc(data_size);
-
-	if ( iargs->noisefilter ) {
-		filter_noise(&image, data_for_measurement);
-	} else {
-		memcpy(data_for_measurement, image.data, data_size);
-	}
-
-	switch ( iargs->peaks ) {
-
-		case PEAK_HDF5:
-		// Get peaks from HDF5
-		if (get_peaks(&image, hdfile,
-			iargs->hdf5_peak_path)) {
-			ERROR("Failed to get peaks from HDF5 file.\n");
-		}
-		if ( !iargs->no_revalidate ) {
-			validate_peaks(&image, iargs->min_int_snr,
-				       iargs->ir_inn, iargs->ir_mid,
-				       iargs->ir_out, iargs->use_saturated);
-		}
-		break;
-
-		case PEAK_ZAEF:
-		search_peaks(&image, iargs->threshold,
-		             iargs->min_gradient, iargs->min_snr,
-		             iargs->ir_inn, iargs->ir_mid, iargs->ir_out,
-		             iargs->use_saturated);
-		break;
-
-	}
-
-	/* Get rid of noise-filtered version at this point
-	 * - it was strictly for the purposes of peak detection. */
-	free(image.data);
-	image.data = data_for_measurement;
-
-	/* Index the pattern */
-	index_pattern(&image, iargs->indm, iargs->ipriv);
-
-	pargs->n_crystals = image.n_crystals;
-
-	/* Default beam parameters */
-	image.div = image.beam->divergence;
-	image.bw = image.beam->bandwidth;
-
-	/* Integrate each crystal's diffraction spots */
-	for ( i=0; i<image.n_crystals; i++ ) {
-
-		RefList *reflections;
-
-		/* Set default crystal parameter(s) */
-		crystal_set_profile_radius(image.crystals[i],
-		                           image.beam->profile_radius);
-
-		if ( iargs->integrate_found ) {
-			reflections = select_intersections(&image,
-			                                   image.crystals[i]);
-		} else {
-			reflections = find_intersections(&image,
-			                                 image.crystals[i]);
-		}
-
-		crystal_set_reflections(image.crystals[i], reflections);
-
-	}
-
-	/* Integrate all the crystals at once - need all the crystals so that
-	 * overlaps can be detected. */
-	integrate_reflections(&image, iargs->closer,
-	                              iargs->bgsub,
-	                              iargs->min_int_snr,
-	                              iargs->ir_inn,
-	                              iargs->ir_mid,
-	                              iargs->ir_out,
-	                              iargs->integrate_saturated,
-	                              iargs->res_cutoff);
-
-	write_chunk(st, &image, hdfile,
-	            iargs->stream_peaks, iargs->stream_refls);
-
-	for ( i=0; i<image.n_crystals; i++ ) {
-		crystal_free(image.crystals[i]);
-	}
-
-	free(image.data);
-	if ( image.flags != NULL ) free(image.flags);
-	image_feature_list_free(image.features);
-	hdfile_close(hdfile);
-}
-
-
 static void run_work(const struct index_args *iargs,
                      int filename_pipe, int results_pipe, Stream *st,
                      int cookie)
diff --git a/src/im-sandbox.h b/src/im-sandbox.h
index dfddabe6..7df41738 100644
--- a/src/im-sandbox.h
+++ b/src/im-sandbox.h
@@ -31,48 +31,10 @@
  *
  */
 
+#include "index.h"
 #include "stream.h"
-
-enum {
-	PEAK_ZAEF,
-	PEAK_HDF5,
-};
-
-
-/* Information about the indexing process which is common to all patterns */
-struct index_args
-{
-	UnitCell *cell;
-	int cmfilter;
-	int noisefilter;
-	int satcorr;
-	int closer;
-	int bgsub;
-	float threshold;
-	float min_gradient;
-	float min_snr;
-	double min_int_snr;
-	struct detector *det;
-	IndexingMethod *indm;
-	IndexingPrivate **ipriv;
-	int peaks;                /* Peak detection method */
-	float tols[4];
-	struct beam_params *beam;
-	char *element;
-	char *hdf5_peak_path;
-	float ir_inn;
-	float ir_mid;
-	float ir_out;
-	struct copy_hdf5_field *copyme;
-	int integrate_saturated;
-	int use_saturated;
-	int no_revalidate;
-	int integrate_found;
-	int stream_peaks;
-	int stream_refls;
-	int res_cutoff;
-};
-
+#include "cell.h"
+#include "process_image.h"
 
 extern void create_sandbox(struct index_args *iargs, int n_proc, char *prefix,
                            int config_basename, FILE *fh,
diff --git a/src/process_image.c b/src/process_image.c
new file mode 100644
index 00000000..010054c8
--- /dev/null
+++ b/src/process_image.c
@@ -0,0 +1,231 @@
+/*
+ * process_image.c
+ *
+ * The processing pipeline for one image
+ *
+ * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY,
+ *                       a research centre of the Helmholtz Association.
+ *
+ * Authors:
+ *   2010-2013 Thomas White <taw@physics.org>
+ *
+ * This file is part of CrystFEL.
+ *
+ * CrystFEL is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * CrystFEL is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with CrystFEL.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdlib.h>
+#include <hdf5.h>
+#include <gsl/gsl_errno.h>
+
+#include "utils.h"
+#include "hdf5-file.h"
+#include "index.h"
+#include "peaks.h"
+#include "detector.h"
+#include "filters.h"
+#include "thread-pool.h"
+#include "beam-parameters.h"
+#include "geometry.h"
+#include "stream.h"
+#include "reflist-utils.h"
+#include "process_image.h"
+
+
+void process_image(const struct index_args *iargs, struct pattern_args *pargs,
+                   Stream *st, int cookie)
+{
+	float *data_for_measurement;
+	size_t data_size;
+	int check;
+	struct hdfile *hdfile;
+	struct image image;
+	int i;
+	char filename[1024];
+
+	/* Prefix to jump out of temporary folder */
+	if ( pargs->filename[0] != '/' ) {
+		snprintf(filename, 1023, "../../%s", pargs->filename);
+	} else {
+		snprintf(filename, 1023, "%s", pargs->filename);
+	}
+
+	image.features = NULL;
+	image.data = NULL;
+	image.flags = NULL;
+	image.copyme = iargs->copyme;
+	image.id = cookie;
+	image.filename = pargs->filename;  /* Relative to top level */
+	image.beam = iargs->beam;
+	image.det = iargs->det;
+	image.crystals = NULL;
+	image.n_crystals = 0;
+
+	hdfile = hdfile_open(filename);  /* Relative to temporary folder */
+	if ( hdfile == NULL ) return;
+
+	if ( iargs->element != NULL ) {
+
+		int r;
+		r = hdfile_set_image(hdfile, iargs->element);
+		if ( r ) {
+			ERROR("Couldn't select path '%s'\n", iargs->element);
+			hdfile_close(hdfile);
+			return;
+		}
+
+	} else {
+
+		int r;
+		r = hdfile_set_first_image(hdfile, "/");
+		if ( r ) {
+			ERROR("Couldn't select first path\n");
+			hdfile_close(hdfile);
+			return;
+		}
+
+	}
+
+	check = hdf5_read(hdfile, &image, 1);
+	if ( check ) {
+		hdfile_close(hdfile);
+		return;
+	}
+
+	if ( (image.width != image.det->max_fs + 1 )
+	  || (image.height != image.det->max_ss + 1))
+	{
+		ERROR("Image size doesn't match geometry size"
+			" - rejecting image.\n");
+		ERROR("Image size: %i,%i.  Geometry size: %i,%i\n",
+		      image.width, image.height,
+		      image.det->max_fs + 1, image.det->max_ss + 1);
+		hdfile_close(hdfile);
+		return;
+	}
+
+	fill_in_values(image.det, hdfile);
+	fill_in_beam_parameters(image.beam, hdfile);
+
+	image.lambda = ph_en_to_lambda(eV_to_J(image.beam->photon_energy));
+
+	if ( (image.beam->photon_energy < 0.0) || (image.lambda > 1000) ) {
+		/* Error message covers a silly value in the beam file or in
+		 * the HDF5 file. */
+		ERROR("Nonsensical wavelength (%e m or %e eV) value for %s.\n",
+		      image.lambda, image.beam->photon_energy, image.filename);
+		hdfile_close(hdfile);
+		return;
+	}
+
+	if ( iargs->cmfilter ) filter_cm(&image);
+
+	/* Take snapshot of image after CM subtraction but before
+	 * the aggressive noise filter. */
+	data_size = image.width * image.height * sizeof(float);
+	data_for_measurement = malloc(data_size);
+
+	if ( iargs->noisefilter ) {
+		filter_noise(&image, data_for_measurement);
+	} else {
+		memcpy(data_for_measurement, image.data, data_size);
+	}
+
+	switch ( iargs->peaks ) {
+
+		case PEAK_HDF5:
+		// Get peaks from HDF5
+		if (get_peaks(&image, hdfile,
+			iargs->hdf5_peak_path)) {
+			ERROR("Failed to get peaks from HDF5 file.\n");
+		}
+		if ( !iargs->no_revalidate ) {
+			validate_peaks(&image, iargs->min_int_snr,
+				       iargs->ir_inn, iargs->ir_mid,
+				       iargs->ir_out, iargs->use_saturated);
+		}
+		break;
+
+		case PEAK_ZAEF:
+		search_peaks(&image, iargs->threshold,
+		             iargs->min_gradient, iargs->min_snr,
+		             iargs->ir_inn, iargs->ir_mid, iargs->ir_out,
+		             iargs->use_saturated);
+		break;
+
+	}
+
+	/* Get rid of noise-filtered version at this point
+	 * - it was strictly for the purposes of peak detection. */
+	free(image.data);
+	image.data = data_for_measurement;
+
+	/* Index the pattern */
+	index_pattern(&image, iargs->indm, iargs->ipriv);
+
+	pargs->n_crystals = image.n_crystals;
+
+	/* Default beam parameters */
+	image.div = image.beam->divergence;
+	image.bw = image.beam->bandwidth;
+
+	/* Integrate each crystal's diffraction spots */
+	for ( i=0; i<image.n_crystals; i++ ) {
+
+		RefList *reflections;
+
+		/* Set default crystal parameter(s) */
+		crystal_set_profile_radius(image.crystals[i],
+		                           image.beam->profile_radius);
+
+		if ( iargs->integrate_found ) {
+			reflections = select_intersections(&image,
+			                                   image.crystals[i]);
+		} else {
+			reflections = find_intersections(&image,
+			                                 image.crystals[i]);
+		}
+
+		crystal_set_reflections(image.crystals[i], reflections);
+
+	}
+
+	/* Integrate all the crystals at once - need all the crystals so that
+	 * overlaps can be detected. */
+	integrate_reflections(&image, iargs->closer,
+	                              iargs->bgsub,
+	                              iargs->min_int_snr,
+	                              iargs->ir_inn,
+	                              iargs->ir_mid,
+	                              iargs->ir_out,
+	                              iargs->integrate_saturated,
+	                              iargs->res_cutoff);
+
+	write_chunk(st, &image, hdfile,
+	            iargs->stream_peaks, iargs->stream_refls);
+
+	for ( i=0; i<image.n_crystals; i++ ) {
+		crystal_free(image.crystals[i]);
+	}
+
+	free(image.data);
+	if ( image.flags != NULL ) free(image.flags);
+	image_feature_list_free(image.features);
+	hdfile_close(hdfile);
+}
diff --git a/src/process_image.h b/src/process_image.h
new file mode 100644
index 00000000..f2f034d9
--- /dev/null
+++ b/src/process_image.h
@@ -0,0 +1,94 @@
+/*
+ * process_image.h
+ *
+ * The processing pipeline for one image
+ *
+ * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY,
+ *                       a research centre of the Helmholtz Association.
+ *
+ * Authors:
+ *   2010-2013 Thomas White <taw@physics.org>
+ *
+ * This file is part of CrystFEL.
+ *
+ * CrystFEL is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * CrystFEL is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with CrystFEL.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+#ifndef PROCESS_IMAGE_H
+#define PROCESS_IMAGE_H
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+
+enum {
+	PEAK_ZAEF,
+	PEAK_HDF5,
+};
+
+
+/* Information about the indexing process which is common to all patterns */
+struct index_args
+{
+	UnitCell *cell;
+	int cmfilter;
+	int noisefilter;
+	int satcorr;
+	int closer;
+	int bgsub;
+	float threshold;
+	float min_gradient;
+	float min_snr;
+	double min_int_snr;
+	struct detector *det;
+	IndexingMethod *indm;
+	IndexingPrivate **ipriv;
+	int peaks;                /* Peak detection method */
+	float tols[4];
+	struct beam_params *beam;
+	char *element;
+	char *hdf5_peak_path;
+	float ir_inn;
+	float ir_mid;
+	float ir_out;
+	struct copy_hdf5_field *copyme;
+	int integrate_saturated;
+	int use_saturated;
+	int no_revalidate;
+	int integrate_found;
+	int stream_peaks;
+	int stream_refls;
+	int res_cutoff;
+};
+
+
+/* Information about the indexing process for one pattern */
+struct pattern_args
+{
+	/* "Input" */
+	char *filename;
+
+	/* "Output" */
+	int n_crystals;
+};
+
+
+extern void process_image(const struct index_args *iargs,
+                          struct pattern_args *pargs, Stream *st,
+                          int cookie);
+
+
+#endif	/* PROCESS_IMAGEs_H */