Convert prediction and integration to detgeom

5 files changed, 39 insertions, 49 deletions

diff --git a/libcrystfel/src/geometry.c b/libcrystfel/src/geometry.c
index 8f8a9672..a1b76a87 100644
--- a/libcrystfel/src/geometry.c
+++ b/libcrystfel/src/geometry.c
@@ -52,7 +52,7 @@
 /** \file geometry.h */
 
 static int locate_peak_on_panel(double x, double y, double z, double k,
-                                struct panel *p,
+                                struct detgeom_panel *p,
                                 double *pfs, double *pss)
 {
 	double ctt, tta, phi;
@@ -85,7 +85,7 @@ static int locate_peak_on_panel(double x, double y, double z, double k,
 	gsl_matrix_set(M, 1, 0, p->cny);
 	gsl_matrix_set(M, 1, 1, p->fsy);
 	gsl_matrix_set(M, 1, 2, p->ssy);
-	gsl_matrix_set(M, 2, 0, p->clen*p->res);
+	gsl_matrix_set(M, 2, 0, p->cnz);
 	gsl_matrix_set(M, 2, 1, p->fsz);
 	gsl_matrix_set(M, 2, 2, p->ssz);
 
@@ -113,7 +113,8 @@ static int locate_peak_on_panel(double x, double y, double z, double k,
 }
 
 static signed int locate_peak(double x, double y, double z, double k,
-                              struct detector *det, double *pfs, double *pss)
+                              struct detgeom *det,
+                              double *pfs, double *pss)
 {
 	int i;
 
@@ -121,7 +122,7 @@ static signed int locate_peak(double x, double y, double z, double k,
 
 	for ( i=0; i<det->n_panels; i++ ) {
 
-		struct panel *p;
+		struct detgeom_panel *p;
 
 		p = &det->panels[i];
 
@@ -379,29 +380,31 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst,
 	/* 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 ( (image->det != NULL) && (updateme != NULL) ) {
+	if ( (image->detgeom != NULL) && (updateme != NULL) ) {
 
 		double fs, ss;
+		assert(get_panel_number(updateme) <= image->detgeom->n_panels);
 		locate_peak_on_panel(xl, yl, zl, mean_kpred,
-		                     get_panel(updateme), &fs, &ss);
+		                     &image->detgeom->panels[get_panel_number(updateme)],
+		                     &fs, &ss);
 		set_detector_pos(refl, fs, ss);
 
 	}
 
-	/* Otherwise, calculate position if we have a detector structure, and
+	/* 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 */
+		double fs, ss;        /* position on detector */
+		signed int p;         /* panel number */
 		p = locate_peak(xl, yl, zl, mean_kpred,
-		                image->det, &fs, &ss);
+		                image->detgeom, &fs, &ss);
 		if ( p == -1 ) {
 			reflection_free(refl);
 			return NULL;
 		}
 		set_detector_pos(refl, fs, ss);
-		set_panel(refl, &image->det->panels[p]);
+		set_panel_number(refl, p);
 
 	}
 
diff --git a/libcrystfel/src/integration.c b/libcrystfel/src/integration.c
index 3d8258c1..681aef02 100644
--- a/libcrystfel/src/integration.c
+++ b/libcrystfel/src/integration.c
@@ -52,6 +52,7 @@
 #include "image.h"
 #include "peaks.h"
 #include "integration.h"
+#include "detgeom.h"
 
 
 /** \file integration.h */
@@ -110,7 +111,7 @@ struct peak_box
 	enum boxmask_val *bm;  /* Box mask */
 
 	int pn;           /* Panel number */
-	struct panel *p;  /* The panel itself */
+	struct detgeom_panel *p;  /* The panel itself */
 
 	/* Fitted background parameters */
 	double a;
@@ -273,9 +274,6 @@ static void show_peak_box(struct intcontext *ic, struct peak_box *bx,
 
 	get_indices(bx->refl, &h, &k, &l);
 	get_detector_pos(bx->refl, &fs, &ss);
-	/* Convert coordinates to match arrangement of panels in HDF5 file */
-	fs = fs + bx->p->orig_min_fs;
-	ss = ss + bx->p->orig_min_ss;
 	printw("Indices %i %i %i\nPanel %s\nPosition fs = %.1f, ss = %.1f\n\n",
 	       h, k, l, bx->p->name, fs, ss);
 
@@ -806,7 +804,7 @@ static int check_box(struct intcontext *ic, struct peak_box *bx, int *sat)
 		int fs, ss;
 		double hd, kd, ld;
 		signed int h, k, l;
-		struct rvec dv;
+		double dv[3];
 		float lsat;
 
 		fs = bx->cfs + p;
@@ -864,10 +862,11 @@ static int check_box(struct intcontext *ic, struct peak_box *bx, int *sat)
 
 		/* Ignore if this pixel is closer to the next reciprocal lattice
 		 * point */
-		dv = get_q_for_panel(bx->p, fs, ss, NULL, ic->k);
-		hd = dv.u * adx + dv.v * ady + dv.w * adz;
-		kd = dv.u * bdx + dv.v * bdy + dv.w * bdz;
-		ld = dv.u * cdx + dv.v * cdy + dv.w * cdz;
+		detgeom_transform_coords(bx->p, fs, ss,
+		                         1.0/ic->k, dv);
+		hd = dv[0] * adx + dv[1] * ady + dv[2] * adz;
+		kd = dv[0] * bdx + dv[1] * bdy + dv[2] * bdz;
+		ld = dv[0] * cdx + dv[1] * cdy + dv[2] * cdz;
 		h = lrint(hd);
 		k = lrint(kd);
 		l = lrint(ld);
@@ -1273,7 +1272,6 @@ static void setup_profile_boxes(struct intcontext *ic, RefList *list)
 		double pfs, pss;
 		struct peak_box *bx;
 		int pn;
-		struct panel *p;
 		int fid_fs, fid_ss;  /* Center coordinates, rounded,
 		                      * in overall data block */
 		int cfs, css;  /* Corner coordinates */
@@ -1283,12 +1281,7 @@ static void setup_profile_boxes(struct intcontext *ic, RefList *list)
 		set_redundancy(refl, 0);
 
 		get_detector_pos(refl, &pfs, &pss);
-		p = get_panel(refl);
-		pn = panel_number(ic->image->det, p);
-		if ( pn == ic->image->det->n_panels ) {
-			ERROR("Couldn't find panel %p\n", p);
-			continue;
-		}
+		pn = get_panel_number(refl);
 
 		/* Explicit truncation of digits after the decimal point.
 		 * This is actually the correct thing to do here, not
@@ -1307,7 +1300,7 @@ static void setup_profile_boxes(struct intcontext *ic, RefList *list)
 		bx->refl = refl;
 		bx->cfs = cfs;
 		bx->css = css;
-		bx->p = p;
+		bx->p = &ic->image->detgeom->panels[pn];
 		bx->pn = pn;
 
 		/* Which reference profile? */
@@ -1406,7 +1399,6 @@ static int integrate_rings_once(Reflection *refl, struct image *image,
 	double pfs, pss;
 	struct peak_box *bx;
 	int pn;
-	struct panel *p;
 	int fid_fs, fid_ss;  /* Center coordinates, rounded,
 	                      * in overall data block */
 	int cfs, css;  /* Corner coordinates */
@@ -1419,12 +1411,7 @@ static int integrate_rings_once(Reflection *refl, struct image *image,
 	set_redundancy(refl, 0);
 
 	get_detector_pos(refl, &pfs, &pss);
-	p = get_panel(refl);
-	pn = panel_number(image->det, p);
-	if ( pn == image->det->n_panels ) {
-		ERROR("Couldn't find panel %p\n", p);
-		return 1;
-	}
+	pn = get_panel_number(refl);
 
 	/* Explicit truncation of digits after the decimal point.
 	 * This is actually the correct thing to do here, not
@@ -1442,7 +1429,7 @@ static int integrate_rings_once(Reflection *refl, struct image *image,
 	bx->refl = refl;
 	bx->cfs = cfs;
 	bx->css = css;
-	bx->p = p;
+	bx->p = &image->detgeom->panels[pn];
 	bx->pn = pn;
 
 	if ( ic->meth & INTEGRATION_CENTER ) {
diff --git a/libcrystfel/src/reflist.c b/libcrystfel/src/reflist.c
index 9cedfc86..df8f24fc 100644
--- a/libcrystfel/src/reflist.c
+++ b/libcrystfel/src/reflist.c
@@ -53,7 +53,7 @@ struct _refldata {
 	/* Location in image */
 	double fs;
 	double ss;
-	struct panel *panel;
+	int panel_number;
 
 	/* Non-zero if this reflection can be used for scaling */
 	int scalable;
@@ -308,12 +308,13 @@ void get_detector_pos(const Reflection *refl, double *fs, double *ss)
 /**
  * \param refl: Reflection
  *
- * \returns the panel which the reflection appears on
+ * \returns panel number (index in detgeom/DataTemplate structure)
+ *           which the reflection appears on
  *
  **/
-struct panel *get_panel(const Reflection *refl)
+int get_panel_number(const Reflection *refl)
 {
-	return refl->data.panel;
+	return refl->data.panel_number;
 }
 
 
@@ -589,14 +590,13 @@ void set_detector_pos(Reflection *refl, double fs, double ss)
 
 /**
  * \param refl: Reflection
- * \param p: Pointer to the panel structure on which the reflection appears
- *
- * Note that the pointer will be stored, not the contents of the structure.
+ * \param pn: Panel number (index in detgeom/DataTemplate structure) of
+ *      the panel on which the reflection appears.
  *
  **/
-void set_panel(Reflection *refl, struct panel *p)
+void set_panel_number(Reflection *refl, int pn)
 {
-	refl->data.panel = p;
+	refl->data.panel_number = pn;
 }
 
 
diff --git a/libcrystfel/src/reflist.h b/libcrystfel/src/reflist.h
index 3085a2aa..91f740a5 100644
--- a/libcrystfel/src/reflist.h
+++ b/libcrystfel/src/reflist.h
@@ -7,7 +7,7 @@
  *                       a research centre of the Helmholtz Association.
  *
  * Authors:
- *   2011-2018 Thomas White <taw@physics.org>
+ *   2011-2020 Thomas White <taw@physics.org>
  *
  * This file is part of CrystFEL.
  *
@@ -109,7 +109,7 @@ extern Reflection *next_found_refl(Reflection *refl);
 
 /* Get */
 extern void get_detector_pos(const Reflection *refl, double *fs, double *ss);
-extern struct panel *get_panel(const Reflection *refl);
+extern int get_panel_number(const Reflection *refl);
 extern double get_partiality(const Reflection *refl);
 extern double get_khalf(const Reflection *refl);
 extern double get_kpred(const Reflection *refl);
@@ -134,7 +134,7 @@ extern struct reflection_contributions *get_contributions(const Reflection *refl
 /* Set */
 extern void copy_data(Reflection *to, const Reflection *from);
 extern void set_detector_pos(Reflection *refl, double fs, double ss);
-extern void set_panel(Reflection *refl, struct panel *p);
+extern void set_panel_number(Reflection *refl, int pn);
 extern void set_kpred(Reflection *refl, double kpred);
 extern void set_khalf(Reflection *refl, double khalf);
 extern void set_exerr(Reflection *refl, double exerr);
diff --git a/tests/integration_check.c b/tests/integration_check.c
index da808cce..50e99d65 100644
--- a/tests/integration_check.c
+++ b/tests/integration_check.c
@@ -124,7 +124,7 @@ int main(int argc, char *argv[])
 		list = reflist_new();
 		refl = add_refl(list, 0, 0, 0);
 		set_detector_pos(refl, 64, 64);
-		set_panel(refl, &image.det->panels[0]);
+		set_panel_number(refl, 0);
 		cell = cell_new();
 		cell_set_lattice_type(cell, L_CUBIC);
 		cell_set_centering(cell, 'P');