Make panel assignments invariant during prediction- and post-refinement

Reflections appearing and disappearing are problematic when trying to do a
least-squares refinement.  Therefore, assume that reflections stay on
panel and keep them under consideration even if their partialities go to
zero (i.e. they drift off Bragg).  This should stabilise both
refinements, and simplifies quite a lot of code.

Collateral "damage": the old "select_intersection()" is now gone.

7 files changed, 129 insertions, 206 deletions

diff --git a/libcrystfel/src/geometry.c b/libcrystfel/src/geometry.c
index b55b696d..35492fde 100644
--- a/libcrystfel/src/geometry.c
+++ b/libcrystfel/src/geometry.c
@@ -47,58 +47,66 @@
 #include "geometry.h"
 
 
-static signed int locate_peak(double x, double y, double z, double k,
-                              struct detector *det, double *xdap, double *ydap)
+static int locate_peak_on_panel(double x, double y, double z, double k,
+                                struct panel *p,
+                                double *pfs, double *pss)
 {
-	int i;
-	signed int found = -1;
 	const double den = k + z;
+	double fs, ss, plx, ply, xd, yd;
 
-	*xdap = -1;  *ydap = -1;
+	/* Coordinates of peak relative to central beam, in m */
+	xd = p->clen * x / den;
+	yd = p->clen * y / den;
 
-	for ( i=0; i<det->n_panels; i++ ) {
+	/* Convert to pixels */
+	xd *= p->res;
+	yd *= p->res;
 
-		double xd, yd;
-		double fs, ss, plx, ply;
-		struct panel *p;
+	/* Convert to relative to the panel corner */
+	plx = xd - p->cnx;
+	ply = yd - p->cny;
+
+	fs = p->xfs*plx + p->yfs*ply;
+	ss = p->xss*plx + p->yss*ply;
+
+	fs += p->min_fs;
+	ss += p->min_ss;
+
+	*pfs = fs;  *pss = ss;
+
+	/* Now, is this on this panel? */
+	if ( fs < p->min_fs ) return 0;
+	if ( fs > p->max_fs ) return 0;
+	if ( ss < p->min_ss ) return 0;
+	if ( ss > p->max_ss ) return 0;
 
-		p = &det->panels[i];
 
-		/* Coordinates of peak relative to central beam, in m */
-		xd = p->clen * x / den;
-		yd = p->clen * y / den;
+	return 1;
+}
 
-		/* Convert to pixels */
-		xd *= p->res;
-		yd *= p->res;
+static signed int locate_peak(double x, double y, double z, double k,
+                              struct detector *det, double *pfs, double *pss)
+{
+	int i;
 
-		/* Convert to relative to the panel corner */
-		plx = xd - p->cnx;
-		ply = yd - p->cny;
+	*pfs = -1;  *pss = -1;
 
-		fs = p->xfs*plx + p->yfs*ply;
-		ss = p->xss*plx + p->yss*ply;
+	for ( i=0; i<det->n_panels; i++ ) {
 
-		fs += p->min_fs;
-		ss += p->min_ss;
+		struct panel *p;
+
+		p = &det->panels[i];
 
-		/* Now, is this on this panel? */
-		if ( fs < p->min_fs ) continue;
-		if ( fs > p->max_fs ) continue;
-		if ( ss < p->min_ss ) continue;
-		if ( ss > p->max_ss ) continue;
+		if ( locate_peak_on_panel(x, y, z, k, p, pfs, pss) ) {
 
-		/* If peak appears on multiple panels, reject it */
-		if ( found != -1 ) return -1;
+			/* Woohoo! */
+			return i;
 
-		/* Woohoo! */
-		found = i;
-		*xdap = fs;
-		*ydap = ss;
+		}
 
 	}
 
-	return found;
+	return -1;
 }
 
 
@@ -182,7 +190,8 @@ static double partiality(PartialityModel pmodel,
 static Reflection *check_reflection(struct image *image, Crystal *cryst,
                                     PartialityModel pmodel,
                                     signed int h, signed int k, signed int l,
-                                    double xl, double yl, double zl)
+                                    double xl, double yl, double zl,
+                                    Reflection *updateme)
 {
 	const int output = 0;
 	double tl;
@@ -195,7 +204,7 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst,
 	double del;
 
 	/* Don't predict 000 */
-	if ( abs(h)+abs(k)+abs(l) == 0 ) return NULL;
+	if ( (updateme == NULL) && (abs(h)+abs(k)+abs(l) == 0) ) return NULL;
 
 	pr = crystal_get_profile_radius(cryst);
 	del = image->div + crystal_get_mosaicity(cryst);
@@ -207,7 +216,7 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst,
 	khigh = 1.0/(image->lambda + image->lambda*image->bw/2.0);
 
 	/* If the point is looking "backscattery", reject it straight away */
-	if ( zl < -khigh/2.0 ) return NULL;
+	if ( (updateme == NULL) && (zl < -khigh/2.0) ) return NULL;
 
 	tl = sqrt(xl*xl + yl*yl);
 
@@ -220,28 +229,47 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst,
 	rlow = klow - distance(cet, cez, tl, zl);  /* Loss of precision */
 
 	/* Condition for reflection to be excited at all */
-	if ( (signbit(rlow) == signbit(rhigh))
+	if ( (updateme == NULL)
+	     && (signbit(rlow) == signbit(rhigh))
 	     && (fabs(rlow) > pr)
 	     && (fabs(rhigh) > pr) ) return NULL;
 
 	/* Calculate partiality */
 	part = partiality(pmodel, rlow, rhigh, pr);
 
-	/* Add peak to list */
-	refl = reflection_new(h, k, l);
+	if ( updateme == NULL ) {
+		refl = reflection_new(h, k, l);
+	} else {
+		refl = updateme;
+	}
+
+	/* If we are updating a previous reflection, assume it stays
+	 * on the same panel and calculate the new position even if it's
+	 * fallen off the edge of the panel. */
+	if ( updateme != NULL ) {
+
+		double fs, ss;
+		locate_peak_on_panel(xl, yl, zl, 1.0/image->lambda,
+		                     get_panel(updateme), &fs, &ss);
+		set_detector_pos(refl, fs, ss);
+
+	}
 
-	/* If we have detector information, check the spot is measured.
-	 * Otherwise, we make do with calculating the partialiaty etc. */
-	if ( image->det != NULL ) {
-		double xda, yda;        /* Position on detector */
-		signed int p;           /* Panel number */
-		p = locate_peak(xl, yl, zl, 1.0/image->lambda, image->det,
-		                &xda, &yda);
+	/* Otherwise, calculate position if we have a detector structure, and
+	 * if we don't then just make do with partiality calculation */
+	if ( (image->det != NULL) && (updateme == NULL) ) {
+
+		double fs, ss;        /* Position on detector */
+		signed int p;         /* Panel number */
+		p = locate_peak(xl, yl, zl, 1.0/image->lambda,
+		                image->det, &fs, &ss);
 		if ( p == -1 ) {
 			reflection_free(refl);
 			return NULL;
 		}
-		set_detector_pos(refl, 0.0, xda, yda);
+		set_detector_pos(refl, fs, ss);
+		set_panel(refl, &image->det->panels[p]);
+
 	}
 
 	if ( unlikely(rlow < rhigh) ) {
@@ -249,8 +277,11 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst,
 		ERROR("%3i %3i %3i  rlow = %e, rhigh = %e\n",
 		      h, k, l, rlow, rhigh);
 		ERROR("div + m = %e, R = %e, bw = %e\n", del, pr, image->bw);
-		reflection_free(refl);
-		return NULL;
+		/* If we are updating, this is (kind of) OK */
+		if ( updateme == NULL ) {
+			reflection_free(refl);
+			return NULL;
+		}
 	}
 
 	set_partial(refl, rlow, rhigh, part);
@@ -339,7 +370,7 @@ RefList *find_intersections_to_res(struct image *image, Crystal *cryst,
 		zl = h*asz + k*bsz + l*csz;
 
 		refl = check_reflection(image, cryst, pmodel,
-		                        h, k, l, xl, yl, zl);
+		                        h, k, l, xl, yl, zl, NULL);
 
 		if ( refl != NULL ) {
 			add_refl_to_list(refl, reflections);
@@ -353,67 +384,6 @@ RefList *find_intersections_to_res(struct image *image, Crystal *cryst,
 }
 
 
-/* Deprecated: select reflections using Kirian-style pixel proximity */
-RefList *select_intersections(struct image *image, Crystal *cryst)
-{
-	double ax, ay, az;
-	double bx, by, bz;
-	double cx, cy, cz;
-	const double min_dist = 0.25;
-	RefList *list;
-	int i;
-
-	/* Round towards nearest */
-	fesetround(1);
-
-	/* Cell basis vectors for this image */
-	cell_get_cartesian(crystal_get_cell(cryst), &ax, &ay, &az,
-	                   &bx, &by, &bz, &cx, &cy, &cz);
-
-	list = reflist_new();
-	if ( list == NULL ) return NULL;
-
-	/* 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;
-		double dsq;
-
-		f = image_get_feature(image->features, i);
-		if ( f == NULL ) continue;
-
-		/* 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 */
-		dsq = pow(h-hd, 2.0) + pow(k-kd, 2.0) + pow(l-ld, 2.0);
-
-		if ( sqrt(dsq) < min_dist ) {
-
-			Reflection *refl;
-
-			refl = add_refl(list, h, k, l);
-			set_detector_pos(refl, sqrt(dsq), f->fs, f->ss);
-
-		}
-
-	}
-
-	return list;
-}
-
-
 static void set_unity_partialities(Crystal *cryst)
 {
 	Reflection *refl;
@@ -430,8 +400,7 @@ static void set_unity_partialities(Crystal *cryst)
 
 
 /* Calculate partialities and apply them to the image's reflections */
-void update_partialities_2(Crystal *cryst, PartialityModel pmodel,
-                           int *n_gained, int *n_lost, double *mean_p_change)
+void update_partialities(Crystal *cryst, PartialityModel pmodel)
 {
 	Reflection *refl;
 	RefListIterator *iter;
@@ -439,8 +408,6 @@ void update_partialities_2(Crystal *cryst, PartialityModel pmodel,
 	double bsx, bsy, bsz;
 	double csx, csy, csz;
 	struct image *image = crystal_get_image(cryst);
-	double total_p_change = 0.0;
-	int n = 0;
 
 	if ( pmodel == PMODEL_UNITY ) {
 		set_unity_partialities(cryst);
@@ -454,69 +421,20 @@ void update_partialities_2(Crystal *cryst, PartialityModel pmodel,
 	      refl != NULL;
 	      refl = next_refl(refl, iter) )
 	{
-		Reflection *vals;
-		double r1, r2, L, p, x, y;
 		double xl, yl, zl;
 		signed int h, k, l;
-		double old_p;
 
 		get_symmetric_indices(refl, &h, &k, &l);
-		old_p = get_partiality(refl);
 
 		/* Get the coordinates of the reciprocal lattice point */
 		xl = h*asx + k*bsx + l*csx;
 		yl = h*asy + k*bsy + l*csy;
 		zl = h*asz + k*bsz + l*csz;
 
-		vals = check_reflection(image, cryst, pmodel,
-		                        h, k, l, xl, yl, zl);
-
-		if ( vals == NULL ) {
-
-			if ( get_redundancy(refl) != 0 ) {
-				(*n_lost)++;
-				set_partiality(refl, 0.0);
-				set_redundancy(refl, 0);
-			}
-
-		} else {
-
-			if ( get_redundancy(refl) == 0 ) {
-				(*n_gained)++;
-				set_redundancy(refl, 1);
-			}
-
-			/* Transfer partiality stuff */
-			get_partial(vals, &r1, &r2, &p);
-			set_partial(refl, r1, r2, p);
-			L = get_lorentz(vals);
-			set_lorentz(refl, L);
-
-			/* Transfer detector location */
-			get_detector_pos(vals, &x, &y);
-			set_detector_pos(refl, 0.0, x, y);
-
-			reflection_free(vals);
-
-			total_p_change += fabs(p - old_p);
-			n++;
-
-		}
+		check_reflection(image, cryst, pmodel,
+		                 h, k, l, xl, yl, zl, refl);
 
 	}
-
-	*mean_p_change = total_p_change / n;
-}
-
-
-/* Wrapper to maintain API compatibility */
-void update_partialities(Crystal *cryst, PartialityModel pmodel)
-{
-	int n_gained = 0;
-	int n_lost = 0;
-	double mean_p_change = 0.0;
-	update_partialities_2(cryst, pmodel, &n_gained, &n_lost,
-	                      &mean_p_change);
 }
 
 
diff --git a/libcrystfel/src/geometry.h b/libcrystfel/src/geometry.h
index f39b6f8b..162c0355 100644
--- a/libcrystfel/src/geometry.h
+++ b/libcrystfel/src/geometry.h
@@ -3,12 +3,12 @@
  *
  * Geometry of diffraction
  *
- * Copyright © 2013-2014 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2013-2015 Deutsches Elektronen-Synchrotron DESY,
  *                       a research centre of the Helmholtz Association.
  * Copyright © 2012 Richard Kirian
  *
  * Authors:
- *   2010-2014 Thomas White <taw@physics.org>
+ *   2010-2015 Thomas White <taw@physics.org>
  *   2012      Richard Kirian
  *
  * This file is part of CrystFEL.
@@ -67,13 +67,7 @@ extern RefList *find_intersections_to_res(struct image *image, Crystal *cryst,
                                           PartialityModel pmodel,
 					  double max_res);
 
-/* Deprecated: select reflections using Kirian-style pixel proximity */
-extern RefList *select_intersections(struct image *image, Crystal *cryst);
-
 extern void update_partialities(Crystal *cryst, PartialityModel pmodel);
-extern void update_partialities_2(Crystal *cryst, PartialityModel pmodel,
-                                  int *n_gained, int *n_lost,
-                                  double *mean_p_change);
 extern void polarisation_correction(RefList *list, UnitCell *cell,
                                     struct image *image);
 
diff --git a/libcrystfel/src/integration.c b/libcrystfel/src/integration.c
index a05f64f1..0e35037d 100644
--- a/libcrystfel/src/integration.c
+++ b/libcrystfel/src/integration.c
@@ -1193,7 +1193,7 @@ static void integrate_prof2d_once(struct intcontext *ic, struct peak_box *bx,
 		get_detector_pos(bx->refl, &pfs, &pss);
 		pfs += bx->offs_fs;
 		pss += bx->offs_ss;
-		set_detector_pos(bx->refl, 0.0, pfs, pss);
+		set_detector_pos(bx->refl, pfs, pss);
 
 		if ( bx->intensity < -5.0*bx->sigma ) {
 			ic->n_implausible++;
@@ -1447,7 +1447,7 @@ static void integrate_rings_once(Reflection *refl, struct image *image,
 	/* Update position */
 	pfs += bx->offs_fs;
 	pss += bx->offs_ss;
-	set_detector_pos(refl, 0.0, pfs, pss);
+	set_detector_pos(refl, pfs, pss);
 
 	if ( get_int_diag(ic, refl) ) show_peak_box(ic, bx, results_pipe);
 
diff --git a/libcrystfel/src/reflist-utils.c b/libcrystfel/src/reflist-utils.c
index 722f800c..82943158 100644
--- a/libcrystfel/src/reflist-utils.c
+++ b/libcrystfel/src/reflist-utils.c
@@ -390,7 +390,7 @@ RefList *read_reflections_from_file(FILE *fh)
 
 			refl = add_refl(out, h, k, l);
 			set_intensity(refl, intensity);
-			set_detector_pos(refl, 0.0, fs, ss);
+			set_detector_pos(refl, fs, ss);
 
 			set_esd_intensity(refl, sigma);
 			set_redundancy(refl, cts);
diff --git a/libcrystfel/src/reflist.c b/libcrystfel/src/reflist.c
index 33c1f948..795894a1 100644
--- a/libcrystfel/src/reflist.c
+++ b/libcrystfel/src/reflist.c
@@ -75,10 +75,7 @@ struct _refldata {
 	/* Location in image */
 	double fs;
 	double ss;
-
-	/* The distance from the exact Bragg position to the coordinates
-	 * given above. */
-	double excitation_error;
+	struct panel *panel;
 
 	/* Non-zero if this reflection can be used for scaling */
 	int scalable;
@@ -316,29 +313,31 @@ Reflection *next_found_refl(Reflection *refl)
 
 /********************************** Getters ***********************************/
 
+
 /**
- * get_excitation_error:
+ * get_detector_pos:
  * @refl: A %Reflection
+ * @fs: Location at which to store the fast scan offset of the reflection
+ * @ss: Location at which to store the slow scan offset of the reflection
  *
- * Returns: The excitation error for the reflection.
  **/
-double get_excitation_error(const Reflection *refl)
+void get_detector_pos(const Reflection *refl, double *fs, double *ss)
 {
-	return refl->data.excitation_error;
+	*fs = refl->data.fs;
+	*ss = refl->data.ss;
 }
 
 
 /**
- * get_detector_pos:
+ * get_panel:
  * @refl: A %Reflection
- * @fs: Location at which to store the fast scan offset of the reflection
- * @ss: Location at which to store the slow scan offset of the reflection
+ *
+ * Returns: the panel which the reflection appears on
  *
  **/
-void get_detector_pos(const Reflection *refl, double *fs, double *ss)
+struct panel *get_panel(const Reflection *refl)
 {
-	*fs = refl->data.fs;
-	*ss = refl->data.ss;
+	return refl->data.panel;
 }
 
 
@@ -570,20 +569,32 @@ void copy_data(Reflection *to, const Reflection *from)
 /**
  * set_detector_pos:
  * @refl: A %Reflection
- * @exerr: The excitation error for this reflection
  * @fs: The fast scan offset of the reflection
  * @ss: The slow scan offset of the reflection
  *
  **/
-void set_detector_pos(Reflection *refl, double exerr, double fs, double ss)
+void set_detector_pos(Reflection *refl, double fs, double ss)
 {
-	refl->data.excitation_error = exerr;
 	refl->data.fs = fs;
 	refl->data.ss = ss;
 }
 
 
 /**
+ * set_panel:
+ * @refl: A %Reflection
+ * @panel: Pointer to the panel structure on which the reflection appears
+ *
+ * Note that the pointer will be stored, not the contents of the structure.
+ *
+ **/
+void set_panel(Reflection *refl, struct panel *p)
+{
+	refl->data.panel = p;
+}
+
+
+/**
  * set_partial:
  * @refl: A %Reflection
  * @rlow: The "low" excitation error
diff --git a/libcrystfel/src/reflist.h b/libcrystfel/src/reflist.h
index 85a87c54..8b0d8e04 100644
--- a/libcrystfel/src/reflist.h
+++ b/libcrystfel/src/reflist.h
@@ -3,11 +3,11 @@
  *
  * Fast reflection/peak list
  *
- * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY,
  *                       a research centre of the Helmholtz Association.
  *
  * Authors:
- *   2011-2014 Thomas White <taw@physics.org>
+ *   2011-2015 Thomas White <taw@physics.org>
  *
  * This file is part of CrystFEL.
  *
@@ -84,8 +84,8 @@ extern Reflection *find_refl(const RefList *list, signed int h, signed int k, si
 extern Reflection *next_found_refl(Reflection *refl);
 
 /* Get */
-extern double get_excitation_error(const Reflection *refl);
 extern void get_detector_pos(const Reflection *refl, double *fs, double *ss);
+extern struct panel *get_panel(const Reflection *refl);
 extern double get_partiality(const Reflection *refl);
 extern double get_lorentz(const Reflection *refl);
 extern void get_indices(const Reflection *refl,
@@ -106,8 +106,8 @@ extern double get_mean_bg(const Reflection *refl);
 
 /* Set */
 extern void copy_data(Reflection *to, const Reflection *from);
-extern void set_detector_pos(Reflection *refl, double exerr,
-                             double fs, double ss);
+extern void set_detector_pos(Reflection *refl, double fs, double ss);
+extern void set_panel(Reflection *refl, struct panel *p);
 extern void set_partial(Reflection *refl, double rlow, double rhigh, double p);
 extern void set_partiality(Reflection *refl, double p);
 extern void set_lorentz(Reflection *refl, double L);
diff --git a/libcrystfel/src/stream.c b/libcrystfel/src/stream.c
index 8586c850..78ac9a40 100644
--- a/libcrystfel/src/stream.c
+++ b/libcrystfel/src/stream.c
@@ -319,7 +319,7 @@ static RefList *read_stream_reflections_2_3(FILE *fh, struct detector *det)
 				p = find_panel_by_name(det,pn);
 				write_fs = fs - p->orig_min_fs + p->min_fs;
 				write_ss = ss - p->orig_min_ss + p->min_ss;
-				set_detector_pos(refl, 0.0, write_fs, write_ss);
+				set_detector_pos(refl, write_fs, write_ss);
 			}
 			set_esd_intensity(refl, sigma);
 			set_peak(refl, pk);
@@ -383,11 +383,11 @@ static RefList *read_stream_reflections_2_1(FILE *fh, struct detector *det)
 				p = find_orig_panel(det, fs, ss);
 				write_fs = fs - p->orig_min_fs + p->min_fs;
 				write_ss = ss - p->orig_min_ss + p->min_ss;
-				set_detector_pos(refl, 0.0, write_fs, write_ss);
+				set_detector_pos(refl, write_fs, write_ss);
 
 			} else {
 
-				set_detector_pos(refl, 0.0, fs, ss);
+				set_detector_pos(refl, fs, ss);
 
 			}
 			set_esd_intensity(refl, sigma);
@@ -448,11 +448,11 @@ static RefList *read_stream_reflections_2_2(FILE *fh, struct detector *det)
 				p = find_orig_panel(det, fs, ss);
 				write_fs = fs - p->orig_min_fs + p->min_fs;
 				write_ss = ss - p->orig_min_ss + p->min_ss;
-				set_detector_pos(refl, 0.0, write_fs, write_ss);
+				set_detector_pos(refl, write_fs, write_ss);
 
 			} else {
 
-				set_detector_pos(refl, 0.0, fs, ss);
+				set_detector_pos(refl, fs, ss);
 
 			}