Convert make_pixelmap to DataTemplate

NB The coffset is no longer written to the HDF5 file.  This was a
terrible hack anyway.  The resolution field is just as bad, but as least
easy to get from the detgeom.

3 files changed, 239 insertions, 67 deletions

diff --git a/libcrystfel/src/datatemplate.c b/libcrystfel/src/datatemplate.c
index f996c99b..2593a7f9 100644
--- a/libcrystfel/src/datatemplate.c
+++ b/libcrystfel/src/datatemplate.c
@@ -1309,3 +1309,168 @@ void data_template_add_copy_header(DataTemplate *dt,
 	/* FIXME: Add "header" to list of things to copy */
 	STATUS("Adding %s\n", header);
 }
+
+
+/* If possible, i.e. if there are no references to image headers in
+ * 'dt', generate a detgeom structure from it.
+ *
+ * NB This is probably not the function you're looking for!
+ * The detgeom structure should normally come from loading an image,
+ * reading a stream or similar.  This function should only be used
+ * when there is really no data involved, e.g. in make_pixelmap.
+ */
+struct detgeom *data_template_to_detgeom(const DataTemplate *dt)
+{
+	struct detgeom *detgeom;
+	int i;
+
+	if ( dt == NULL ) return NULL;
+
+	detgeom = malloc(sizeof(struct detgeom));
+	if ( detgeom == NULL ) return;
+
+	detgeom->panels = malloc(dt->n_panels*sizeof(struct detgeom_panel));
+	if ( detgeom->panels == NULL ) return;
+
+	detgeom->n_panels = dt->n_panels;
+
+	for ( i=0; i<dt->n_panels; i++ ) {
+
+		detgeom->panels[i].name = safe_strdup(dt->panels[i].name);
+
+		detgeom->panels[i].pixel_pitch = dt->panels[i].pixel_pitch;
+
+		/* NB cnx,cny are in pixels, cnz is in m */
+		detgeom->panels[i].cnx = dt->panels[i].cnx;
+		detgeom->panels[i].cny = dt->panels[i].cny;
+		detgeom->panels[i].cnz = get_length(image, dt->panels[i].cnz_from);
+
+		/* Apply offset (in m) and then convert cnz from
+		 * m to pixels */
+		detgeom->panels[i].cnz += dt->panels[i].cnz_offset;
+		detgeom->panels[i].cnz /= detgeom->panels[i].pixel_pitch;
+
+		detgeom->panels[i].max_adu = dt->panels[i].max_adu;
+		detgeom->panels[i].adu_per_photon = 1.0;  /* FIXME ! */
+
+		detgeom->panels[i].w = dt->panels[i].orig_max_fs
+		                        - dt->panels[i].orig_min_fs + 1;
+		detgeom->panels[i].h = dt->panels[i].orig_max_ss
+		                        - dt->panels[i].orig_min_ss + 1;
+
+		detgeom->panels[i].fsx = dt->panels[i].fsx;
+		detgeom->panels[i].fsy = dt->panels[i].fsy;
+		detgeom->panels[i].fsz = dt->panels[i].fsz;
+		detgeom->panels[i].ssx = dt->panels[i].ssx;
+		detgeom->panels[i].ssy = dt->panels[i].ssy;
+		detgeom->panels[i].ssz = dt->panels[i].ssz;
+
+	}
+
+	return detgeom;
+}
+
+
+int data_template_get_slab_extents(const DataTemplate *dt,
+                                   int *pw, int *ph)
+{
+	int w, h;
+	char *data_from;
+
+	data_from = dt->panels[0].data;
+
+	w = 0;  h = 0;
+	for ( i=0; i<dt->n_panels; i++ ) {
+
+		struct panel_template *p = &dt->panels[i];
+
+		if ( strcmp(data_from, p->data) != 0 ) {
+			/* Not slabby */
+			return 1;
+		}
+
+		if ( p->dim_structure != NULL ) {
+			int j;
+			for ( j=0; j<p->dim_structure->ndims; j++ ) {
+				if ( p->dim_structure->dims[j] == HYSL_PLACEHOLDER ) {
+					/* Not slabby */
+					return 1;
+				}
+			}
+		}
+
+		if ( p->file_max_fs > w ) {
+			w = p->file_max_fs;
+		}
+		if ( p->file_max_ss > h ) {
+			h = p->file_max_ss;
+		}
+	}
+
+	/* Inclusive -> exclusive */
+	*pw = w + 1;
+	*ph = h + 1;
+	return 0;
+}
+
+
+/* Return non-zero if pixel fs,ss on panel p is in a bad region
+ * as specified in the geometry file (regions only, not including
+ * masks, NaN/inf, no_index etc */
+int data_template_in_bad_region(DataTemplate *dtempl,
+                                int pn, double fs, double ss)
+{
+	double rx, ry;
+	double xs, ys;
+	int i;
+	struct panel_template *p;
+
+	if ( p >= dtempl->n_panels ) {
+		ERROR("Panel index out of range\n");
+		return 0;
+	}
+	p = dtempl->panels[pn];
+
+	/* Convert xs and ys, which are in fast scan/slow scan coordinates,
+	 * to x and y */
+	xs = fs*p->fsx + ss*p->ssx;
+	ys = fs*p->fsy + ss*p->ssy;
+
+	rx = xs + p->cnx;
+	ry = ys + p->cny;
+
+	for ( i=0; i<dtempl->n_bad; i++ ) {
+
+		struct dt_badregion *b = &dtempl->bad[i];
+
+		if ( (b->panel != NULL)
+		  && (strcmp(b->panel, p->name) != 0) ) continue;
+
+		if ( b->is_fsss ) {
+
+			int nfs, nss;
+
+			/* fs/ss bad regions are specified according
+			 * to the original coordinates */
+			nfs = fs + p->orig_min_fs;
+			nss = ss + p->orig_min_ss;
+
+			if ( nfs < b->min_fs ) continue;
+			if ( nfs > b->max_fs ) continue;
+			if ( nss < b->min_ss ) continue;
+			if ( nss > b->max_ss ) continue;
+
+		} else {
+
+			if ( rx < b->min_x ) continue;
+			if ( rx > b->max_x ) continue;
+			if ( ry < b->min_y ) continue;
+			if ( ry > b->max_y ) continue;
+
+		}
+
+		return 1;
+	}
+
+	return 0;
+}
diff --git a/libcrystfel/src/datatemplate.h b/libcrystfel/src/datatemplate.h
index fd99d5f0..9cb60bf1 100644
--- a/libcrystfel/src/datatemplate.h
+++ b/libcrystfel/src/datatemplate.h
@@ -33,6 +33,7 @@
 #include <config.h>
 #endif
 
+#include "detgeom.h"
 
 /**
  * \file datatemplate.h
@@ -67,6 +68,13 @@ extern int data_template_panel_to_file_coords(const DataTemplate *dt,
 extern void data_template_add_copy_header(DataTemplate *dt,
                                           const char *header);
 
+extern struct detgeom *data_template_to_detgeom(const DataTemplate *dt);
+
+extern int data_template_get_slab_extents(const DataTemplate *dt, int *pw, int *ph);
+
+extern int data_template_in_bad_region(DataTemplate *dtempl,
+                                       int pn, double fs, double ss);
+
 #ifdef __cplusplus
 }
 #endif
diff --git a/src/make_pixelmap.c b/src/make_pixelmap.c
index 1a1ed5f3..eedcb439 100644
--- a/src/make_pixelmap.c
+++ b/src/make_pixelmap.c
@@ -7,7 +7,7 @@
  *                       a research centre of the Helmholtz Association.
  *
  * Authors:
- *   2012-2018 Thomas White <taw@physics.org>
+ *   2012-2020 Thomas White <taw@physics.org>
  *   2016-2016 Omri Mor <omor1@asu.edu>
  *
  * This file is part of CrystFEL.
@@ -41,8 +41,11 @@
 #include <unistd.h>
 #include <getopt.h>
 #include <assert.h>
+#include <hdf5.h>
 
-#include "detector.h"
+#include "utils.h"
+#include "datatemplate.h"
+#include "detgeom.h"
 
 
 static void show_help(const char *s)
@@ -126,7 +129,7 @@ static void create_scalar(hid_t gh, const char *name, float val)
 
 static void write_pixelmap_hdf5(const char *filename,
                                 float *x, float *y, float *z,
-                                int width, int height, float res, float coffset)
+                                int width, int height, float res)
 {
 	hid_t fh;
 
@@ -141,7 +144,6 @@ static void write_pixelmap_hdf5(const char *filename,
 	create_array(fh, "z", z, H5T_NATIVE_FLOAT, width, height);
 
 	create_scalar(fh, "res", res);
-	create_scalar(fh, "coffset", coffset);
 
 	H5Fclose(fh);
 }
@@ -171,12 +173,12 @@ int main(int argc, char *argv[])
 	int c;
 	char *input_file = NULL;
 	char *output_file = NULL;
-	struct detector *det = NULL;
+	DataTemplate *dtempl;
+	struct detgeom *detgeom;
 	int fs, ss, w, h;
 	float *x, *y, *z;
 	uint16_t *b;
-	int i;
-	float res, coffset;
+	float res;
 	int badmap = 0;
 	int good_pixel_val = 1;
 	int bad_pixel_val = 0;
@@ -248,24 +250,28 @@ int main(int argc, char *argv[])
 	}
 
 	/* Load geometry */
-	det = get_detector_geometry(input_file, NULL);
-	if ( det == NULL ) {
-		ERROR("Failed to read geometry from '%s'\n", input_file);
+	dtempl = data_template_new_from_file(input_file);
+	if ( dtempl == NULL ) {
+		ERROR("Failed to read geometry from '%s'\n",
+		      input_file);
 		return 1;
 	}
 	free(input_file);
 
+	detgeom = data_template_to_detgeom(dtempl);
+	if ( detgeom == NULL ) {
+		ERROR("Could not make detector structure.\n");
+		ERROR("Geometry file must not contain references to "
+		      "image header values\n");
+		return 1;
+	}
+
 	/* Determine max orig fs and ss */
-	w = 0;  h = 0;
-	for ( i=0; i<det->n_panels; i++ ) {
-		if ( det->panels[i].orig_max_fs > w ) {
-			w = det->panels[i].orig_max_fs;
-		}
-		if ( det->panels[i].orig_max_ss > h ) {
-			h = det->panels[i].orig_max_ss;
-		}
+	if ( data_template_get_slab_extents(dtempl, &w, &h) ) {
+		ERROR("Pixel map can only be created if all panels "
+		      "have their data in one \"slab\".\n");
+		return 1;
 	}
-	w += 1;  h += 1;  /* Inclusive -> Exclusive */
 	STATUS("Data slab size: %i x %i\n", w, h);
 
 	x = malloc(w*h*sizeof(float));
@@ -277,66 +283,59 @@ int main(int argc, char *argv[])
 		return 1;
 	}
 
+	/* For every pixel in the slab ... */
 	for ( ss=0; ss<h; ss++ ) {
-		for ( fs=0; fs<w; fs++ ) {
-
-			double rx, ry;
-			struct panel *p;
-			double xs, ys;
-			double cfs, css;
-			int nfs, nss;
-
-			p = find_orig_panel(det, fs, ss);
-
-			/* Add half a pixel to fs and ss to get the fs,ss
-			 * coordinates of the CENTRE of the pixel */
-			cfs = fs + 0.5;
-			css = ss + 0.5;
-			xs = (cfs - p->orig_min_fs)*p->fsx
-			      + (css - p->orig_min_ss)*p->ssx;
-			ys = (cfs - p->orig_min_fs)*p->fsy
-			     + (css - p->orig_min_ss)*p->ssy;
-
-			rx = (xs + p->cnx) / p->res;
-			ry = (ys + p->cny) / p->res;
-
-			x[fs + w*ss] = rx;
-			y[fs + w*ss] = ry;
-			z[fs + w*ss] = 0.0;  /* FIXME */
-
-			nfs = fs - p->orig_min_fs;
-			nss = ss - p->orig_min_ss;
-			if ( in_bad_region(det, p, nfs, nss) ) {
-				b[fs + w*ss] = bad_pixel_val;
-			} else {
-				b[fs + w*ss] = good_pixel_val;
-			}
+	for ( fs=0; fs<w; fs++ ) {
+
+		double rx, ry;
+		double xs, ys;
+		float cfs, css;
+		int pn;
+		struct detgeom_panel *p;
+
+		/* Add half a pixel to fs and ss to get the fs,ss
+		 * coordinates of the CENTRE of the pixel */
+		cfs = fs + 0.5;
+		css = ss + 0.5;
+
+		if ( data_template_file_to_panel_coords(dtempl,
+		                                        &cfs, &css,
+		                                        &pn) )
+		{
+			ERROR("Couldn't convert coordinates\n");
+			return 1;
+		}
 
+		p = &detgeom->panels[pn];
+
+		xs = cfs*p->fsx + css*p->ssx;
+		ys = cfs*p->fsy + css*p->ssy;
+
+		rx = (xs + p->cnx) * p->pixel_pitch;
+		ry = (ys + p->cny) * p->pixel_pitch;
+
+		x[fs + w*ss] = rx;
+		y[fs + w*ss] = ry;
+		z[fs + w*ss] = 0.0;  /* 2D part only */
+
+		if ( data_template_in_bad_region(dtempl, pn, cfs, css) ) {
+			b[fs + w*ss] = bad_pixel_val;
+		} else {
+			b[fs + w*ss] = good_pixel_val;
 		}
 
-		progress_bar(ss, h, "Converting");
 	}
-	STATUS("\n");
-
-	res = det->defaults.res;
-	/* use the res of the first panel if not in global definitions
-	 * assumes one of these exist */
-	if ( res == -1.0 ) {
-		res = det->panels[0].res;
 	}
 
-	coffset = det->defaults.coffset;
-	/* use the coffset of the first panel if not in global definitions
-	 * assumes one of these exist*/
-	if ( coffset == 0.0 ) {
-		coffset = det->panels[0].coffset;
-	}
+	res = 1.0 / detgeom->panels[0].pixel_pitch;
 
 	if ( badmap ) {
 		write_badmap_hdf5(output_file, b, w, h);
 	} else {
-		write_pixelmap_hdf5(output_file, x, y, z, w, h, res, coffset);
+		write_pixelmap_hdf5(output_file, x, y, z, w, h, res);
 	}
 
+	data_template_free(dtempl);
+
 	return 0;
 }