Peak integration changes

1 files changed, 113 insertions, 86 deletions

diff --git a/libcrystfel/src/peaks.c b/libcrystfel/src/peaks.c
index dabc3b1a..c350388b 100644
--- a/libcrystfel/src/peaks.c
+++ b/libcrystfel/src/peaks.c
@@ -3,12 +3,12 @@
  *
  * Peak search and other image analysis
  *
- * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY,
  *                  a research centre of the Helmholtz Association.
  * Copyright © 2012 Richard Kirian
  *
  * Authors:
- *   2010-2012 Thomas White <taw@physics.org>
+ *   2010-2013 Thomas White <taw@physics.org>
  *   2012      Kenneth Beyerlein <kenneth.beyerlein@desy.de>
  *   2011      Andrew Martin <andrew.martin@desy.de>
  *   2011      Richard Kirian
@@ -147,24 +147,15 @@ static int cull_peaks(struct image *image)
 }
 
 
-/* cfs, css relative to panel origin */
-static int *make_BgMask(struct image *image, struct panel *p,
-                        double ir_out, double ir_inn)
+static void add_crystal_to_mask(struct image *image, struct panel *p,
+                                double ir_out, double ir_inn, int w, int h,
+                                int *mask, Crystal *cr)
 {
 	Reflection *refl;
 	RefListIterator *iter;
-	int *mask;
-	int w, h;
-
-	w = p->max_fs - p->min_fs + 1;
-	h = p->max_ss - p->min_ss + 1;
-	mask = calloc(w*h, sizeof(int));
-	if ( mask == NULL ) return NULL;
-
-	if ( image->reflections == NULL ) return mask;
 
 	/* Loop over all reflections */
-	for ( refl = first_refl(image->reflections, &iter);
+	for ( refl = first_refl(crystal_get_reflections(cr), &iter);
 	      refl != NULL;
 	      refl = next_refl(refl, iter) )
 	{
@@ -205,6 +196,28 @@ static int *make_BgMask(struct image *image, struct panel *p,
 		}
 
 	}
+}
+
+
+/* cfs, css relative to panel origin */
+static int *make_BgMask(struct image *image, struct panel *p,
+                        double ir_out, double ir_inn)
+{
+	int *mask;
+	int w, h;
+	int i;
+
+	w = p->max_fs - p->min_fs + 1;
+	h = p->max_ss - p->min_ss + 1;
+	mask = calloc(w*h, sizeof(int));
+	if ( mask == NULL ) return NULL;
+
+	if ( image->crystals == NULL ) return mask;
+
+	for ( i=0; i<image->n_crystals; i++ ) {
+		add_crystal_to_mask(image, p, ir_inn, ir_out,
+		                    w, h, mask, image->crystals[i]);
+	}
 
 	return mask;
 }
@@ -581,29 +594,19 @@ void search_peaks(struct image *image, float threshold, float min_gradient,
 }
 
 
-double peak_lattice_agreement(struct image *image, UnitCell *cell, double *pst)
+int peak_sanity_check(struct image *image)
 {
-	int i;
 	int n_feat = 0;
 	int n_sane = 0;
-	double ax, ay, az;
-	double bx, by, bz;
-	double cx, cy, cz;
-	double min_dist = 0.25;
-	double stot = 0.0;
-
-	/* Round towards nearest */
-	fesetround(1);
-
-	/* Cell basis vectors for this image */
-	cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
+	int i;
+	const double min_dist = 0.25;
 
-	/* Loop over peaks, checking proximity to nearest reflection */
 	for ( i=0; i<image_feature_count(image->features); i++ ) {
 
 		struct imagefeature *f;
 		struct rvec q;
 		double h,k,l,hd,kd,ld;
+		int j;
 
 		/* Assume all image "features" are genuine peaks */
 		f = image_get_feature(image->features, i);
@@ -613,38 +616,42 @@ double peak_lattice_agreement(struct image *image, UnitCell *cell, double *pst)
 		/* Reciprocal space position of found peak */
 		q = get_q(image, f->fs, f->ss, NULL, 1.0/image->lambda);
 
-		/* Decimal and fractional Miller indices of nearest
-		 * reciprocal lattice point */
-		hd = q.u * ax + q.v * ay + q.w * az;
-		kd = q.u * bx + q.v * by + q.w * bz;
-		ld = q.u * cx + q.v * cy + q.w * cz;
-		h = lrint(hd);
-		k = lrint(kd);
-		l = lrint(ld);
-
-		/* Check distance */
-		if ( (fabs(h - hd) < min_dist) && (fabs(k - kd) < min_dist)
-		  && (fabs(l - ld) < min_dist) )
-		{
-			double sval;
-			n_sane++;
-			sval = pow(h-hd, 2.0) + pow(k-kd, 2.0) + pow(l-ld, 2.0);
-			stot += 1.0 - sval;
-			continue;
-		}
+		for ( j=0; j<image->n_crystals; j++ ) {
+
+			double ax, ay, az;
+			double bx, by, bz;
+			double cx, cy, cz;
+
+			cell_get_cartesian(crystal_get_cell(image->crystals[j]),
+			                   &ax, &ay, &az,
+			                   &bx, &by, &bz,
+			                   &cx, &cy, &cz);
+
+			/* Decimal and fractional Miller indices of nearest
+			 * reciprocal lattice point */
+			hd = q.u * ax + q.v * ay + q.w * az;
+			kd = q.u * bx + q.v * by + q.w * bz;
+			ld = q.u * cx + q.v * cy + q.w * cz;
+			h = lrint(hd);
+			k = lrint(kd);
+			l = lrint(ld);
+
+			/* Check distance */
+			if ( (fabs(h - hd) < min_dist)
+			  && (fabs(k - kd) < min_dist)
+			  && (fabs(l - ld) < min_dist) )
+			{
+				n_sane++;
+				continue;
+			}
 
-	}
+		}
 
-	*pst = stot;
-	return (double)n_sane / (float)n_feat;
-}
 
+	}
 
-int peak_sanity_check(struct image *image)
-{
-	double stot;
 	/* 0 means failed test, 1 means passed test */
-	return peak_lattice_agreement(image, image->indexed_cell, &stot) >= 0.5;
+	return ((double)n_sane / n_feat) >= 0.5;
 }
 
 
@@ -705,43 +712,28 @@ static struct integr_ind *sort_reflections(RefList *list, UnitCell *cell,
 }
 
 
-/* Integrate the list of predicted reflections in "image" */
-void integrate_reflections(struct image *image, int use_closer, int bgsub,
-                           double min_snr,
-                           double ir_inn, double ir_mid, double ir_out,
-                           int integrate_saturated)
+static void integrate_crystal(Crystal *cr, struct image *image, int use_closer,
+                              int bgsub, double min_snr,
+                              double ir_inn, double ir_mid, double ir_out,
+                              int integrate_saturated, int **bgMasks)
 {
+	RefList *reflections;
 	struct integr_ind *il;
 	int n, i;
 	double av = 0.0;
 	int first = 1;
-	int **bgMasks;
+	int n_saturated = 0;
+
+	reflections = crystal_get_reflections(cr);
 
-	if ( num_reflections(image->reflections) == 0 ) return;
+	if ( num_reflections(reflections) == 0 ) return;
 
-	il = sort_reflections(image->reflections, image->indexed_cell, &n);
+	il = sort_reflections(reflections, crystal_get_cell(cr), &n);
 	if ( il == NULL ) {
 		ERROR("Couldn't sort reflections\n");
 		return;
 	}
 
-	/* Make background masks for all panels */
-	bgMasks = calloc(image->det->n_panels, sizeof(int *));
-	if ( bgMasks == NULL ) {
-		ERROR("Couldn't create list of background masks.\n");
-		return;
-	}
-	for ( i=0; i<image->det->n_panels; i++ ) {
-		int *mask;
-		mask = make_BgMask(image, &image->det->panels[i],
-		                   ir_out, ir_inn);
-		if ( mask == NULL ) {
-			ERROR("Couldn't create background mask.\n");
-			return;
-		}
-		bgMasks[i] = mask;
-	}
-
 	for ( i=0; i<n; i++ ) {
 
 		double fs, ss, intensity;
@@ -809,7 +801,7 @@ void integrate_reflections(struct image *image, int use_closer, int bgsub,
 		                   1, bgMasks[pnum], &saturated);
 
 		if ( saturated ) {
-			image->n_saturated++;
+			n_saturated++;
 			if ( !integrate_saturated ) r = 1;
 		}
 
@@ -844,14 +836,49 @@ void integrate_reflections(struct image *image, int use_closer, int bgsub,
 		}
 	}
 
+	crystal_set_num_saturated_reflections(cr, n_saturated);
+	crystal_set_resolution_limit(cr, 0.0);
+
+	free(il);
+}
+
+
+/* Integrate the list of predicted reflections in "image" */
+void integrate_reflections(struct image *image, int use_closer, int bgsub,
+                           double min_snr,
+                           double ir_inn, double ir_mid, double ir_out,
+                           int integrate_saturated)
+{
+	int i;
+	int **bgMasks;
+
+	/* Make background masks for all panels */
+	bgMasks = calloc(image->det->n_panels, sizeof(int *));
+	if ( bgMasks == NULL ) {
+		ERROR("Couldn't create list of background masks.\n");
+		return;
+	}
+	for ( i=0; i<image->det->n_panels; i++ ) {
+		int *mask;
+		mask = make_BgMask(image, &image->det->panels[i],
+		                   ir_out, ir_inn);
+		if ( mask == NULL ) {
+			ERROR("Couldn't create background mask.\n");
+			return;
+		}
+		bgMasks[i] = mask;
+	}
+
+	for ( i=0; i<image->n_crystals; i++ ) {
+		integrate_crystal(image->crystals[i], image, use_closer,
+		                  bgsub, min_snr, ir_inn, ir_mid, ir_out,
+		                  integrate_saturated, bgMasks);
+	}
+
 	for ( i=0; i<image->det->n_panels; i++ ) {
 		free(bgMasks[i]);
 	}
 	free(bgMasks);
-
-	image->diffracting_resolution = 0.0;
-
-	free(il);
 }