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.

1 files changed, 86 insertions, 168 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);
 }