Initial re-parameterisation of prediction refinement with fs/ss coordinates

The gradients need to be updated.

5 files changed, 48 insertions, 70 deletions

diff --git a/libcrystfel/src/geometry.c b/libcrystfel/src/geometry.c
index 8e8a445d..5289856b 100644
--- a/libcrystfel/src/geometry.c
+++ b/libcrystfel/src/geometry.c
@@ -1086,9 +1086,11 @@ void polarisation_correction(RefList *list, UnitCell *cell,
 }
 
 
-/* Returns dx_h/dP, where P = any parameter */
-double x_gradient(int param, Reflection *refl, UnitCell *cell,
-                  struct detgeom_panel *p)
+/* Returns dfs_refl/dP, where P = any parameter *
+ * fs_refl is the fast scan position of refl in panel 'p',
+ * in metres (not pixels) */
+double fs_gradient(int param, Reflection *refl, UnitCell *cell,
+                   struct detgeom_panel *p)
 {
 	signed int h, k, l;
 	double xl, zl, kpred;
@@ -1140,16 +1142,20 @@ double x_gradient(int param, Reflection *refl, UnitCell *cell,
 		case GPARAM_CLEN :
 		return -xl / (kpred+zl);
 
+		default:
+		ERROR("Positional gradient requested for parameter %i?\n", param);
+		abort();
+
 	}
 
-	ERROR("Positional gradient requested for parameter %i?\n", param);
-	abort();
 }
 
 
-/* Returns dy_h/dP, where P = any parameter */
-double y_gradient(int param, Reflection *refl, UnitCell *cell,
-                  struct detgeom_panel *p)
+/* Returns dss_refl/dP, where P = any parameter
+ * ss_refl is the slow scan position of refl in panel 'p',
+ * in metres (not pixels) */
+double ss_gradient(int param, Reflection *refl, UnitCell *cell,
+                   struct detgeom_panel *p)
 {
 	signed int h, k, l;
 	double yl, zl, kpred;
@@ -1201,8 +1207,9 @@ double y_gradient(int param, Reflection *refl, UnitCell *cell,
 		case GPARAM_CLEN :
 		return -yl / (kpred+zl);
 
-	}
+		default :
+		ERROR("Positional gradient requested for parameter %i?\n", param);
+		abort();
 
-	ERROR("Positional gradient requested for parameter %i?\n", param);
-	abort();
+	}
 }
diff --git a/libcrystfel/src/geometry.h b/libcrystfel/src/geometry.h
index 21bd4d39..a7057f12 100644
--- a/libcrystfel/src/geometry.h
+++ b/libcrystfel/src/geometry.h
@@ -8,7 +8,7 @@
  * Copyright © 2012 Richard Kirian
  *
  * Authors:
- *   2010-2020 Thomas White <taw@physics.org>
+ *   2010-2021 Thomas White <taw@physics.org>
  *   2012      Richard Kirian
  *
  * This file is part of CrystFEL.
@@ -120,10 +120,10 @@ extern void polarisation_correction(RefList *list, UnitCell *cell,
 extern double sphere_fraction(double rlow, double rhigh, double pr);
 extern double gaussian_fraction(double rlow, double rhigh, double pr);
 
-extern double x_gradient(int param, Reflection *refl, UnitCell *cell,
-                         struct detgeom_panel *p);
-extern double y_gradient(int param, Reflection *refl, UnitCell *cell,
-                         struct detgeom_panel *p);
+extern double fs_gradient(int param, Reflection *refl, UnitCell *cell,
+                          struct detgeom_panel *p);
+extern double ss_gradient(int param, Reflection *refl, UnitCell *cell,
+                          struct detgeom_panel *p);
 
 #ifdef __cplusplus
 }
diff --git a/libcrystfel/src/predict-refine.c b/libcrystfel/src/predict-refine.c
index 5210846a..1ed3bfed 100644
--- a/libcrystfel/src/predict-refine.c
+++ b/libcrystfel/src/predict-refine.c
@@ -68,19 +68,6 @@ static const enum gparam rv[] =
 };
 
 
-void twod_mapping(double fs, double ss, double *px, double *py,
-                  struct detgeom_panel *p, double dx, double dy)
-{
-	double xs, ys;
-
-	xs = fs*p->fsx + ss*p->ssx;  /* pixels */
-	ys = fs*p->fsy + ss*p->ssy;  /* pixels */
-
-	*px = (xs + p->cnx) * p->pixel_pitch + dx;  /* metres */
-	*py = (ys + p->cny) * p->pixel_pitch + dy;  /* metres */
-}
-
-
 static double r_dev(struct reflpeak *rp)
 {
 	/* Excitation error term */
@@ -88,29 +75,19 @@ static double r_dev(struct reflpeak *rp)
 }
 
 
-double x_dev(struct reflpeak *rp, struct detgeom *det,
-             double dx, double dy)
+double fs_dev(struct reflpeak *rp, struct detgeom *det)
 {
-	/* Peak position term */
-	double xpk, ypk, xh, yh;
 	double fsh, ssh;
-	twod_mapping(rp->peak->fs, rp->peak->ss, &xpk, &ypk, rp->panel, dx, dy);
 	get_detector_pos(rp->refl, &fsh, &ssh);
-	twod_mapping(fsh, ssh, &xh, &yh, rp->panel, dx, dy);
-	return xh-xpk;
+	return fsh - rp->peak->fs;
 }
 
 
-double y_dev(struct reflpeak *rp, struct detgeom *det,
-             double dx, double dy)
+double ss_dev(struct reflpeak *rp, struct detgeom *det)
 {
-	/* Peak position term */
-	double xpk, ypk, xh, yh;
 	double fsh, ssh;
-	twod_mapping(rp->peak->fs, rp->peak->ss, &xpk, &ypk, rp->panel, dx, dy);
 	get_detector_pos(rp->refl, &fsh, &ssh);
-	twod_mapping(fsh, ssh, &xh, &yh, rp->panel, dx, dy);
-	return yh-ypk;
+	return ssh - rp->peak->ss;
 }
 
 
@@ -397,10 +374,10 @@ static int iterate(struct reflpeak *rps, int n, UnitCell *cell,
 
 		}
 
-		/* Positional x terms */
+		/* Positional fs terms */
 		for ( k=0; k<num_params; k++ ) {
-			gradients[k] = x_gradient(rv[k], rps[i].refl, cell,
-			                          rps[i].panel);
+			gradients[k] = fs_gradient(rv[k], rps[i].refl, cell,
+			                           rps[i].panel);
 		}
 
 		for ( k=0; k<num_params; k++ ) {
@@ -422,17 +399,17 @@ static int iterate(struct reflpeak *rps, int n, UnitCell *cell,
 
 			}
 
-			v_c = x_dev(&rps[i], image->detgeom, *total_x, *total_y);
+			v_c = fs_dev(&rps[i], image->detgeom, *total_x, *total_y);
 			v_c *= -gradients[k];
 			v_curr = gsl_vector_get(v, k);
 			gsl_vector_set(v, k, v_curr + v_c);
 
 		}
 
-		/* Positional y terms */
+		/* Positional ss terms */
 		for ( k=0; k<num_params; k++ ) {
-			gradients[k] = y_gradient(rv[k], rps[i].refl, cell,
-			                          rps[i].panel);
+			gradients[k] = ss_gradient(rv[k], rps[i].refl, cell,
+			                           rps[i].panel);
 		}
 
 		for ( k=0; k<num_params; k++ ) {
@@ -454,7 +431,7 @@ static int iterate(struct reflpeak *rps, int n, UnitCell *cell,
 
 			}
 
-			v_c = y_dev(&rps[i], image->detgeom, *total_x, *total_y);
+			v_c = ss_dev(&rps[i], image->detgeom, *total_x, *total_y);
 			v_c *= -gradients[k];
 			v_curr = gsl_vector_get(v, k);
 			gsl_vector_set(v, k, v_curr + v_c);
@@ -535,13 +512,13 @@ static double pred_residual(struct reflpeak *rps, int n, struct detgeom *det,
 
 	r = 0.0;
 	for ( i=0; i<n; i++ ) {
-		r += pow(x_dev(&rps[i], det, dx, dy), 2.0);
+		r += pow(fs_dev(&rps[i], det, dx, dy), 2.0);
 	}
 	res += r;
 
 	r = 0.0;
 	for ( i=0; i<n; i++ ) {
-		r += pow(y_dev(&rps[i], det, dx, dy), 2.0);
+		r += pow(ss_dev(&rps[i], det, dx, dy), 2.0);
 	}
 	res += r;
 
diff --git a/libcrystfel/src/predict-refine.h b/libcrystfel/src/predict-refine.h
index d1e5fbe7..e74e00a8 100644
--- a/libcrystfel/src/predict-refine.h
+++ b/libcrystfel/src/predict-refine.h
@@ -46,17 +46,12 @@ struct reflpeak {
  * Prediction refinement: refinement of indexing solutions before integration.
  */
 
-extern void twod_mapping(double fs, double ss, double *px, double *py,
-                         struct detgeom_panel *p, double dx, double dy);
-
 extern int refine_prediction(struct image *image, Crystal *cr, Mille *mille);
 
 extern int refine_radius(Crystal *cr, struct image *image);
 
-extern double x_dev(struct reflpeak *rp, struct detgeom *det,
-                    double dx, double dy);
+extern double fs_dev(struct reflpeak *rp, struct detgeom *det);
 
-extern double y_dev(struct reflpeak *rp, struct detgeom *det,
-                    double dx, double dy);
+extern double ss_dev(struct reflpeak *rp, struct detgeom *det);
 
 #endif	/* PREDICT_REFINE_H */
diff --git a/tests/prediction_gradient_check.c b/tests/prediction_gradient_check.c
index 61ea49af..9a5d6ec7 100644
--- a/tests/prediction_gradient_check.c
+++ b/tests/prediction_gradient_check.c
@@ -56,7 +56,7 @@ static void scan(RefList *reflections, RefList *compare,
 	{
 		signed int h, k, l;
 		Reflection *refl2;
-		double fs, ss, xh, yh;
+		double fs, ss;
 		int pn;
 
 		get_indices(refl, &h, &k, &l);
@@ -69,7 +69,6 @@ static void scan(RefList *reflections, RefList *compare,
 
 		get_detector_pos(refl2, &fs, &ss);
 		pn = get_panel_number(refl2);
-		twod_mapping(fs, ss, &xh, &yh, &det->panels[pn], 0, 0);
 
 		switch ( checkrxy ) {
 
@@ -78,11 +77,11 @@ static void scan(RefList *reflections, RefList *compare,
 			break;
 
 			case 1 :
-			vals[idx][i] = xh;
+			vals[idx][i] = fs * det->panels[pn].pixel_pitch;
 			break;
 
 			case 2 :
-			vals[idx][i] = yh;
+			vals[idx][i] = ss * det->panels[pn].pixel_pitch;
 			break;
 		}
 
@@ -286,11 +285,11 @@ static double test_gradients(Crystal *cr, double incr_val, int refine,
 				image = crystal_get_image(cr);
 
 				if ( checkrxy == 1 ) {
-					cgrad = x_gradient(refine, refl,
+					cgrad = fs_gradient(refine, refl,
 					         crystal_get_cell(cr),
 					         &image->detgeom->panels[0]);
 				} else {
-					cgrad = y_gradient(refine, refl,
+					cgrad = ss_gradient(refine, refl,
 					         crystal_get_cell(cr),
 					         &image->detgeom->panels[0]);
 				}
@@ -408,10 +407,10 @@ int main(int argc, char *argv[])
 	image.detgeom->panels[0].max_adu = INFINITY;
 	image.detgeom->panels[0].fsx = 1.0;
 	image.detgeom->panels[0].fsy = 0.0;
-	image.detgeom->panels[0].fsz = 0.0;
+	image.detgeom->panels[0].fsz = 0.5;
 	image.detgeom->panels[0].ssx = 0.0;
 	image.detgeom->panels[0].ssy = 1.0;
-	image.detgeom->panels[0].ssz = 0.0;
+	image.detgeom->panels[0].ssz = 0.5;
 	image.detgeom->panels[0].cnx = -500.0;
 	image.detgeom->panels[0].cny = -500.0;
 	image.detgeom->panels[0].cnz = 1000.0; /* pixels */
@@ -449,11 +448,11 @@ int main(int argc, char *argv[])
 			STATUS("Excitation error:\n");
 			break;
 			case 1:
-			STATUS("x coordinate:\n");
+			STATUS("fs coordinate:\n");
 			break;
 			default:
 			case 2:
-			STATUS("y coordinate:\n");
+			STATUS("ss coordinate:\n");
 			break;
 			STATUS("WTF??\n");
 			break;