1 files changed, 103 insertions, 88 deletions

diff --git a/src/process_image.c b/src/process_image.c
index 6ccd0c9a..783ca464 100644
--- a/src/process_image.c
+++ b/src/process_image.c
@@ -3,11 +3,11 @@
  *
  * The processing pipeline for one image
  *
- * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY,
  *                       a research centre of the Helmholtz Association.
  *
  * Authors:
- *   2010-2014 Thomas White <taw@physics.org>
+ *   2010-2015 Thomas White <taw@physics.org>
  *   2014      Valerio Mariani
  *
  * This file is part of CrystFEL.
@@ -50,79 +50,35 @@
 #include "reflist-utils.h"
 #include "process_image.h"
 #include "integration.h"
+#include "predict-refine.h"
 
 
-static int cmpd2(const void *av, const void *bv)
+static void try_refine_autoR(struct image *image, Crystal *cr)
 {
-	double *ap, *bp;
-	double a, b;
+	double old_R, new_R;
+	char notes[1024];
 
-	ap = (double *)av;
-	bp = (double *)bv;
+	refine_radius(cr, image);
+	old_R = crystal_get_profile_radius(cr);
 
-	a = ap[1];
-	b = bp[1];
-
-	if ( fabs(a) < fabs(b) ) return -1;
-	return 1;
-}
-
-
-static void refine_radius(Crystal *cr)
-{
-	Reflection *refl;
-	RefListIterator *iter;
-	double vals[num_reflections(crystal_get_reflections(cr))*2];
-	int n = 0;
-	int i;
-	double ti = 0.0;  /* Total intensity */
-
-	for ( refl = first_refl(crystal_get_reflections(cr), &iter);
-	      refl != NULL;
-	      refl = next_refl(refl, iter) )
-	{
-		double rlow, rhigh, p;
-		int val = 0;
-
-		if ( get_intensity(refl) > 9.0*get_esd_intensity(refl) ) {
-			val = 1;
-		}
-
-		get_partial(refl, &rlow, &rhigh, &p);
-
-		vals[(2*n)+0] = val;
-		vals[(2*n)+1] = fabs((rhigh+rlow)/2.0);
-		n++;
-
-	}
-
-	/* Sort in ascending order of absolute "deviation from Bragg" */
-	qsort(vals, n, sizeof(double)*2, cmpd2);
-
-	/* Calculate cumulative number of very strong reflections as a function
-	 * of absolute deviation from Bragg */
-	for ( i=0; i<n-1; i++ ) {
-		ti += vals[2*i];
-		vals[2*i] = ti;
-	}
-
-	if ( ti < 10 ) {
-		ERROR("WARNING: Not enough strong reflections (%.0f) to estimate "
-		      "crystal parameters (trying anyway).\n", ti);
+	if ( refine_prediction(image, cr) ) {
+		crystal_set_user_flag(cr, 1);
+		return;
 	}
 
-	/* Find the cutoff where we get 90% of the strong spots */
-	for ( i=0; i<n-1; i++ ) {
-		if ( vals[2*i] > 0.90*ti ) break;
-	}
+	/* Estimate radius again with better geometry */
+	refine_radius(cr, image);
+	new_R = crystal_get_profile_radius(cr);
 
-	crystal_set_profile_radius(cr, fabs(vals[2*i+1]));
+	snprintf(notes, 1024, "predict_refine/R old = %.5f new = %.5f nm^-1",
+	                      old_R/1e9, new_R/1e9);
+	crystal_add_notes(cr, notes);
 }
 
 
 void process_image(const struct index_args *iargs, struct pattern_args *pargs,
                    Stream *st, int cookie, const char *tmpdir, int results_pipe,
-                   int serial)
+                   int serial, pthread_mutex_t *term_lock)
 {
 	float *data_for_measurement;
 	size_t data_size;
@@ -133,6 +89,7 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs,
 	int r;
 	int ret;
 	char *rn;
+	int n_crystals_left;
 
 	image.features = NULL;
 	image.data = NULL;
@@ -142,10 +99,11 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs,
 	image.filename = pargs->filename_p_e->filename;
 	image.event = pargs->filename_p_e->ev;
 	image.beam = iargs->beam;
-	image.det = iargs->det;
+	image.det = copy_geom(iargs->det);
 	image.crystals = NULL;
 	image.n_crystals = 0;
 	image.serial = serial;
+	image.indexed_by = INDEXING_NONE;
 
 	hdfile = hdfile_open(image.filename);
 	if ( hdfile == NULL ) {
@@ -177,10 +135,7 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs,
 	switch ( iargs->peaks ) {
 
 		case PEAK_HDF5:
-		/* Get peaks from HDF5 */
-
-		if ( get_peaks(&image, hdfile, iargs->hdf5_peak_path,
-		               iargs->cxi_hdf5_peaks, pargs->filename_p_e) ) {
+		if ( get_peaks(&image, hdfile, iargs->hdf5_peak_path) ) {
 			ERROR("Failed to get peaks from HDF5 file.\n");
 		}
 		if ( !iargs->no_revalidate ) {
@@ -191,6 +146,19 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs,
 		}
 		break;
 
+		case PEAK_CXI:
+		if ( get_peaks_cxi(&image, hdfile, iargs->hdf5_peak_path,
+		                   pargs->filename_p_e) ) {
+			ERROR("Failed to get peaks from CXI file.\n");
+		}
+		if ( !iargs->no_revalidate ) {
+			validate_peaks(&image, iargs->min_snr,
+				       iargs->pk_inn, iargs->pk_mid,
+		                       iargs->pk_out, iargs->use_saturated,
+				       iargs->check_hdf5_snr);
+		}
+		break;
+
 		case PEAK_ZAEF:
 		search_peaks(&image, iargs->threshold,
 		             iargs->min_gradient, iargs->min_snr,
@@ -226,36 +194,73 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs,
 	}
 	free(rn);
 
-	pargs->n_crystals = image.n_crystals;
 	for ( i=0; i<image.n_crystals; i++ ) {
 		crystal_set_image(image.crystals[i], &image);
+		crystal_set_user_flag(image.crystals[i], 0);
 	}
 
-	/* Default parameters */
-	image.div = 0.0;
-	image.bw = 0.00000001;
-	for ( i=0; i<image.n_crystals; i++ ) {
-		crystal_set_profile_radius(image.crystals[i], 0.01e9);
-		crystal_set_mosaicity(image.crystals[i], 0.0);  /* radians */
+	/* Set beam/crystal parameters */
+	if ( iargs->fix_divergence >= 0.0 ) {
+		image.div = iargs->fix_divergence;
+	} else {
+		image.div = 0.0;
+	}
+	if ( iargs->fix_bandwidth >= 0.0 ) {
+		image.bw = iargs->fix_bandwidth;
+	} else {
+		image.bw = 0.00000001;
+	}
+	if ( iargs->fix_profile_r >= 0.0 ) {
+		for ( i=0; i<image.n_crystals; i++ ) {
+			crystal_set_profile_radius(image.crystals[i],
+			                           iargs->fix_profile_r);
+			crystal_set_mosaicity(image.crystals[i], 0.0);
+		}
+	} else {
+		for ( i=0; i<image.n_crystals; i++ ) {
+			crystal_set_profile_radius(image.crystals[i], 0.02e9);
+			crystal_set_mosaicity(image.crystals[i], 0.0);
+		}
 	}
 
-	/* Integrate all the crystals at once - need all the crystals so that
-	 * overlaps can be detected. */
-	integrate_all_4(&image, iargs->int_meth, PMODEL_SCSPHERE, iargs->push_res,
-	                iargs->ir_inn, iargs->ir_mid, iargs->ir_out,
-	                INTDIAG_NONE, 0, 0, 0, results_pipe);
+	if ( iargs->fix_profile_r < 0.0 ) {
+
+		for ( i=0; i<image.n_crystals; i++ ) {
+			if ( iargs->predict_refine ) {
+				try_refine_autoR(&image, image.crystals[i]);
+			} else {
+				refine_radius(image.crystals[i], &image);
+			}
+		}
+
+	} else {
 
+		for ( i=0; i<image.n_crystals; i++ ) {
+			if ( iargs->predict_refine ) {
+				refine_prediction(&image, image.crystals[i]);
+			}
+		}
+
+	}
+
+	/* If there are no crystals left, set the indexing flag back to zero */
+	n_crystals_left = 0;
 	for ( i=0; i<image.n_crystals; i++ ) {
-		refine_radius(image.crystals[i]);
-		reflist_free(crystal_get_reflections(image.crystals[i]));
+		if ( crystal_get_user_flag(image.crystals[i]) == 0 ) {
+			n_crystals_left++;
+		}
+	}
+	if ( n_crystals_left == 0 ) {
+		image.indexed_by = INDEXING_NONE;
 	}
 
+	/* Integrate! */
 	integrate_all_4(&image, iargs->int_meth, PMODEL_SCSPHERE,
-		                iargs->push_res,
-			        iargs->ir_inn, iargs->ir_mid, iargs->ir_out,
-			        iargs->int_diag, iargs->int_diag_h,
-			        iargs->int_diag_k, iargs->int_diag_l,
-			        results_pipe);
+	                iargs->push_res,
+	                iargs->ir_inn, iargs->ir_mid, iargs->ir_out,
+	                iargs->int_diag, iargs->int_diag_h,
+	                iargs->int_diag_k, iargs->int_diag_l,
+	                term_lock);
 
 	ret = write_chunk(st, &image, hdfile,
 	                  iargs->stream_peaks, iargs->stream_refls,
@@ -269,8 +274,17 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs,
 		n += crystal_get_num_implausible_reflections(image.crystals[i]);
 	}
 	if ( n > 0 ) {
-		STATUS("WARNING: %i implausibly negative reflection%s in %s.\n",
-		       n, n>1?"s":"", image.filename);
+		STATUS("WARNING: %i implausibly negative reflection%s in %s "
+		       "%s\n", n, n>1?"s":"", image.filename,
+		       get_event_string(image.event));
+	}
+
+	/* Count crystals which are still good */
+	pargs->n_crystals = 0;
+	for ( i=0; i<image.n_crystals; i++ ) {
+		if ( crystal_get_user_flag(image.crystals[i]) == 0 ) {
+			pargs->n_crystals++;
+		}
 	}
 
 	for ( i=0; i<image.n_crystals; i++ ) {
@@ -290,5 +304,6 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs,
 	free(image.data);
 	if ( image.flags != NULL ) free(image.flags);
 	image_feature_list_free(image.features);
+	free_detector_geometry(image.det);
 	hdfile_close(hdfile);
 }