Get detector geometry from file

10 files changed, 179 insertions, 66 deletions

diff --git a/src/detector.c b/src/detector.c
index 0748b2ca..988c77fb 100644
--- a/src/detector.c
+++ b/src/detector.c
@@ -33,7 +33,7 @@ int map_position(struct image *image, double dx, double dy,
 	double x = 0.0;
 	double y = 0.0;
 
-	p = find_panel(&image->det, dx, dy);
+	p = find_panel(image->det, dx, dy);
 	if ( p == NULL ) return 1;
 
 	x = ((double)dx - p->cx);
@@ -94,7 +94,7 @@ void record_image(struct image *image, int do_poisson)
 			ERROR("NaN at %i,%i\n", x, y);
 		}
 
-		p = find_panel(&image->det, x, y);
+		p = find_panel(image->det, x, y);
 
 		/* Area of one pixel */
 		pix_area = pow(1.0/p->res, 2.0);
@@ -158,3 +158,98 @@ struct panel *find_panel(struct detector *det, int x, int y)
 
 	return NULL;
 }
+
+
+struct detector *get_detector_geometry(const char *filename)
+{
+	FILE *fh;
+	struct detector *det;
+	char *rval;
+	char **bits;
+	int n_panels = -1;
+	int i;
+
+	fh = fopen(filename, "r");
+	if ( fh == NULL ) return NULL;
+
+	det = malloc(sizeof(struct detector));
+	if ( det == NULL ) {
+		fclose(fh);
+		return NULL;
+	}
+
+	do {
+
+		int n1, n2;
+		char **path;
+		char line[1024];
+		int np;
+
+		rval = fgets(line, 1023, fh);
+		if ( rval == NULL ) break;
+		chomp(line);
+
+		n1 = assplode(line, " \t", &bits, ASSPLODE_NONE);
+		if ( n1 < 3 ) continue;
+
+		if ( bits[1][0] != '=' ) continue;
+
+		if ( strcmp(bits[0], "n_panels") == 0 ) {
+			if ( n_panels != -1 ) {
+				ERROR("Duplicate n_panels statement.\n");
+				fclose(fh);
+				free(det);
+				return NULL;
+			}
+			n_panels = atoi(bits[2]);
+			det->panels = malloc(n_panels * sizeof(struct panel));
+			continue;
+		}
+
+		n2 = assplode(bits[0], "/\\.", &path, ASSPLODE_NONE);
+
+		np = atoi(path[0]);
+
+		if ( strcmp(path[1], "min_x") == 0 ) {
+			det->panels[np].min_x = atof(bits[2]);
+		} else if ( strcmp(path[1], "max_x") == 0 ) {
+			det->panels[np].max_x = atof(bits[2]);
+		} else if ( strcmp(path[1], "min_y") == 0 ) {
+			det->panels[np].min_y = atof(bits[2]);
+		} else if ( strcmp(path[1], "max_y") == 0 ) {
+			det->panels[np].max_y = atof(bits[2]);
+		} else if ( strcmp(path[1], "cx") == 0 ) {
+			det->panels[np].cx = atof(bits[2]);
+		} else if ( strcmp(path[1], "cy") == 0 ) {
+			det->panels[np].cy = atof(bits[2]);
+		} else if ( strcmp(path[1], "clen") == 0 ) {
+			det->panels[np].clen = atof(bits[2]);
+		} else if ( strcmp(path[1], "res") == 0 ) {
+			det->panels[np].res = atof(bits[2]);
+		} else {
+			ERROR("Unrecognised field '%s'\n", path[1]);
+		}
+
+	} while ( rval != NULL );
+
+	if ( n_panels == -1 ) {
+		ERROR("No panel descriptions in geometry file.\n");
+		fclose(fh);
+		free(det->panels);
+		free(det);
+		return NULL;
+	}
+
+	for ( i=0; i<n_panels; i++ ) {
+		STATUS("Panel %i, min_x = %i\n", i, det->panels[i].min_x);
+		STATUS("Panel %i, max_x = %i\n", i, det->panels[i].max_x);
+		STATUS("Panel %i, min_y = %i\n", i, det->panels[i].min_y);
+		STATUS("Panel %i, max_y = %i\n", i, det->panels[i].max_y);
+		STATUS("Panel %i, cx = %f\n", i, det->panels[i].cx);
+		STATUS("Panel %i, cy = %f\n", i, det->panels[i].cy);
+		STATUS("Panel %i, clen = %f\n", i, det->panels[i].clen);
+		STATUS("Panel %i, res = %f\n", i, det->panels[i].res);
+	}
+
+	return det;
+}
diff --git a/src/detector.h b/src/detector.h
index d6e740a0..057b70a8 100644
--- a/src/detector.h
+++ b/src/detector.h
@@ -47,4 +47,6 @@ extern void record_image(struct image *image, int do_poisson);
 
 extern struct panel *find_panel(struct detector *det, int x, int y);
 
+extern struct detector *get_detector_geometry(const char *filename);
+
 #endif	/* DETECTOR_H */
diff --git a/src/diffraction.c b/src/diffraction.c
index 4adf5932..fb512e62 100644
--- a/src/diffraction.c
+++ b/src/diffraction.c
@@ -348,7 +348,7 @@ struct rvec get_q(struct image *image, unsigned int xs, unsigned int ys,
 	const unsigned int x = xs / sampling;
 	const unsigned int y = ys / sampling; /* Integer part only */
 
-	p = find_panel(&image->det, x, y);
+	p = find_panel(image->det, x, y);
 	assert(p != NULL);
 
 	rx = ((float)xs - (sampling*p->cx)) / (sampling * p->res);
@@ -377,7 +377,7 @@ double get_tt(struct image *image, unsigned int xs, unsigned int ys)
 	const unsigned int x = xs;
 	const unsigned int y = ys; /* Integer part only */
 
-	p = find_panel(&image->det, x, y);
+	p = find_panel(image->det, x, y);
 
 	rx = ((float)xs - p->cx) / p->res;
 	ry = ((float)ys - p->cy) / p->res;
diff --git a/src/geometry-lcls.tmp b/src/geometry-lcls.tmp
deleted file mode 100644
index 50bf9ef0..00000000
--- a/src/geometry-lcls.tmp
+++ /dev/null
@@ -1,23 +0,0 @@
-/* Set up detector configuration */
-image.det.n_panels = 2;
-image.det.panels = malloc(image.det.n_panels*sizeof(struct panel));
-
-/* Upper panel */
-image.det.panels[0].min_x = 0;
-image.det.panels[0].max_x = 1023;
-image.det.panels[0].min_y = 512;
-image.det.panels[0].max_y = 1023;
-image.det.panels[0].cx = 491.9;
-image.det.panels[0].cy = 440.7;
-image.det.panels[0].clen = 67.8e-3;
-image.det.panels[0].res = 13333.3;  /* 75 micron pixel size */
-
-/* Lower panel */
-image.det.panels[1].min_x = 0;
-image.det.panels[1].max_x = 1023;
-image.det.panels[1].min_y = 0;
-image.det.panels[1].max_y = 511;
-image.det.panels[1].cx = 492.0;
-image.det.panels[1].cy = 779.7;
-image.det.panels[1].clen = 70.8e-3;
-image.det.panels[1].res = 13333.3;  /* 75 micron pixel size */
diff --git a/src/geometry-simple.tmp b/src/geometry-simple.tmp
deleted file mode 100644
index 3472ff8a..00000000
--- a/src/geometry-simple.tmp
+++ /dev/null
@@ -1,23 +0,0 @@
-/* Set up detector configuration */
-image.det.n_panels = 2;
-image.det.panels = malloc(image.det.n_panels*sizeof(struct panel));
-
-/* Upper panel */
-image.det.panels[0].min_x = 0;
-image.det.panels[0].max_x = 1023;
-image.det.panels[0].min_y = 512;
-image.det.panels[0].max_y = 1023;
-image.det.panels[0].cx = 512.0;
-image.det.panels[0].cy = 502.0;
-image.det.panels[0].clen = 50.0e-3;
-image.det.panels[0].res = 13333.3;  /* 75 micron pixel size */
-
-/* Lower panel */
-image.det.panels[1].min_x = 0;
-image.det.panels[1].max_x = 1023;
-image.det.panels[1].min_y = 0;
-image.det.panels[1].max_y = 511;
-image.det.panels[1].cx = 512.0;
-image.det.panels[1].cy = 522.0;
-image.det.panels[1].clen = 50.0e-3;
-image.det.panels[1].res = 13333.3;  /* 75 micron pixel size */
diff --git a/src/image.h b/src/image.h
index 70809d22..3b31b2fa 100644
--- a/src/image.h
+++ b/src/image.h
@@ -86,7 +86,7 @@ struct image {
 	UnitCell                *indexed_cell;
 	UnitCell                *candidate_cells[MAX_CELL_CANDIDATES];
 	int                     ncells;
-	struct detector         det;
+	struct detector         *det;
 	const char              *filename;
 	struct reflhit          *hits;
 	int                     n_hits;
diff --git a/src/indexamajig.c b/src/indexamajig.c
index 0fb97bfe..39ceea2b 100644
--- a/src/indexamajig.c
+++ b/src/indexamajig.c
@@ -59,6 +59,7 @@ struct process_args
 	int config_nomatch;
 	int config_unpolar;
 	int config_sanity;
+	struct detector *det;
 	IndexingMethod indm;
 	const double *intensities;
 	const unsigned int *counts;
@@ -86,6 +87,7 @@ static void show_help(const char *s)
 "      --indexing=<method> Use 'method' for indexing.  Choose from:\n"
 "                           none     : no indexing\n"
 "                           dirax    : invoke DirAx\n"
+"  -g. --geometry=<file>   Get detector geometry from file.\n"
 "\n\nWith just the above options, this program does not do much of practical "
 "use.\nYou should also enable some of the following:\n\n"
 "      --near-bragg        Output a list of reflection intensities to stdout.\n"
@@ -159,7 +161,7 @@ static struct image *get_simage(struct image *template, int alternate)
 	 * - not necessarily the same as the original. */
 	image->width = 1024;
 	image->height = 1024;
-	image->det.n_panels = 2;
+	image->det->n_panels = 2;
 
 	if ( alternate ) {
 
@@ -183,12 +185,12 @@ static struct image *get_simage(struct image *template, int alternate)
 		panels[1].clen = 56.7e-2;  /* 56.7 cm */
 		panels[1].res = 13333.3;   /* 75 microns/pixel */
 
-		image->det.panels = panels;
+		image->det->panels = panels;
 
 	} else {
 
 		/* Copy pointer to old geometry */
-		image->det.panels = template->det.panels;
+		image->det->panels = template->det->panels;
 
 	}
 
@@ -263,6 +265,7 @@ static void *process_image(void *pargsv)
 	image.filename = filename;
 	image.hits = NULL;
 	image.n_hits = 0;
+	image.det = pargs->det;
 
 	/* View head-on (unit cell is tilted) */
 	image.orientation.w = 1.0;
@@ -285,8 +288,6 @@ static void *process_image(void *pargsv)
 		return result;
 	}
 
-	#include "geometry-lcls.tmp"
-
 	hdf5_read(hdfile, &image);
 
 	if ( config_cmfilter ) {
@@ -372,7 +373,6 @@ static void *process_image(void *pargsv)
 done:
 	free(image.data);
 	free(image.flags);
-	free(image.det.panels);
 	image_feature_list_free(image.features);
 	free(image.hits);
 	hdfile_close(hdfile);
@@ -409,6 +409,8 @@ int main(int argc, char *argv[])
 	int config_alternate = 0;
 	int config_unpolar = 0;
 	int config_sanity = 0;
+	struct detector *det;
+	char *geometry = NULL;
 	IndexingMethod indm;
 	char *indm_str = NULL;
 	UnitCell *cell;
@@ -435,6 +437,7 @@ int main(int argc, char *argv[])
 		{"near-bragg",         0, &config_nearbragg,   1},
 		{"write-drx",          0, &config_writedrx,    1},
 		{"indexing",           1, NULL,               'z'},
+		{"geometry",           1, NULL,               'g'},
 		{"simulate",           0, &config_simulate,    1},
 		{"filter-cm",          0, &config_cmfilter,    1},
 		{"filter-noise",       0, &config_noisefilter, 1},
@@ -450,7 +453,8 @@ int main(int argc, char *argv[])
 	};
 
 	/* Short options */
-	while ((c = getopt_long(argc, argv, "hi:wp:j:x:", longopts, NULL)) != -1) {
+	while ((c = getopt_long(argc, argv, "hi:wp:j:x:g:",
+	                        longopts, NULL)) != -1) {
 
 		switch (c) {
 		case 'h' : {
@@ -488,6 +492,11 @@ int main(int argc, char *argv[])
 			break;
 		}
 
+		case 'g' : {
+			geometry = strdup(optarg);
+			break;
+		}
+
 		case 0 : {
 			break;
 		}
@@ -550,6 +559,18 @@ int main(int argc, char *argv[])
 	}
 	free(indm_str);
 
+	if ( geometry == NULL ) {
+		ERROR("You need to specify a geometry file with --geometry\n");
+		return 1;
+	}
+
+	det = get_detector_geometry(geometry);
+	if ( det == NULL ) {
+		ERROR("Failed to read detector geometry from '%s'\n", geometry);
+		return 1;
+	}
+	free(geometry);
+
 	cell = load_cell_from_pdb(pdb);
 	if ( cell == NULL ) {
 		if ( pdb == NULL ) {
@@ -603,6 +624,7 @@ int main(int argc, char *argv[])
 		pargs->config_unpolar = config_unpolar;
 		pargs->config_sanity = config_sanity;
 		pargs->cell = cell;
+		pargs->det = det;
 		pargs->indm = indm;
 		pargs->intensities = intensities;
 		pargs->counts = counts;
@@ -695,6 +717,8 @@ int main(int argc, char *argv[])
 	}
 
 	free(prefix);
+	free(det->panels);
+	free(det);
 	free(cell);
 	fclose(fh);
 
diff --git a/src/pattern_sim.c b/src/pattern_sim.c
index 65999718..44963b9d 100644
--- a/src/pattern_sim.c
+++ b/src/pattern_sim.c
@@ -47,6 +47,7 @@ static void show_help(const char *s)
 "                            intensities file)\n"
 "     --simulation-details  Show technical details of the simulation.\n"
 "     --gpu                 Use the GPU to speed up the calculation.\n"
+"  -g. --geometry=<file>    Get detector geometry from file.\n"
 "\n"
 "     --near-bragg          Output h,k,l,I near Bragg conditions.\n"
 " -n, --number=<N>          Generate N images.  Default 1.\n"
@@ -187,6 +188,7 @@ int main(int argc, char *argv[])
 	char *filename = NULL;
 	char *grad_str = NULL;
 	char *outfile = NULL;
+	char *geometry = NULL;
 	GradientMethod grad;
 	int ndone = 0;    /* Number of simulations done (images or not) */
 	int number = 1;   /* Number used for filename of image */
@@ -208,14 +210,16 @@ int main(int argc, char *argv[])
 		{"no-noise",           0, &config_nonoise,     1},
 		{"intensities",        1, NULL,               'i'},
 		{"powder",             1, NULL,               'w'},
-		{"gradients",          1, NULL,               'g'},
+		{"gradients",          1, NULL,               't'},
 		{"pdb",                1, NULL,               'p'},
 		{"output",             1, NULL,               'o'},
+		{"geometry",           1, NULL,               'g'},
 		{0, 0, NULL, 0}
 	};
 
 	/* Short options */
-	while ((c = getopt_long(argc, argv, "hrn:i:g:p:o:", longopts, NULL)) != -1) {
+	while ((c = getopt_long(argc, argv, "hrn:i:t:p:o:",
+	                        longopts, NULL)) != -1) {
 
 		switch (c) {
 		case 'h' : {
@@ -238,7 +242,7 @@ int main(int argc, char *argv[])
 			break;
 		}
 
-		case 'g' : {
+		case 't' : {
 			grad_str = strdup(optarg);
 			break;
 		}
@@ -258,6 +262,11 @@ int main(int argc, char *argv[])
 			break;
 		}
 
+		case 'g' : {
+			geometry = strdup(optarg);
+			break;
+		}
+
 		case 0 : {
 			break;
 		}
@@ -329,6 +338,11 @@ int main(int argc, char *argv[])
 		return 1;
 	}
 
+	if ( geometry == NULL ) {
+		ERROR("You need to specify a geometry file with --geometry\n");
+		return 1;
+	}
+
 	/* Define image parameters */
 	image.width = 1024;
 	image.height = 1024;
@@ -343,7 +357,12 @@ int main(int argc, char *argv[])
 	image.f0 = 1.0;
 	image.f0_available = 1;
 
-	#include "geometry-lcls.tmp"
+	image.det = get_detector_geometry(geometry);
+	if ( image.det == NULL ) {
+		ERROR("Failed to read detector geometry from '%s'\n", geometry);
+		return 1;
+	}
+	free(geometry);
 
 	powder = calloc(image.width*image.height, sizeof(*powder));
 
@@ -461,7 +480,8 @@ skip:
 		cleanup_gpu(gctx);
 	}
 
-	free(image.det.panels);
+	free(image.det->panels);
+	free(image.det);
 	free(powder);
 	free(cell);
 	free(intensities);
diff --git a/src/peaks.c b/src/peaks.c
index a5dce6c1..bcd1fd12 100644
--- a/src/peaks.c
+++ b/src/peaks.c
@@ -166,7 +166,7 @@ static int integrate_peak(struct image *image, int xp, int yp,
 			if ( !((flags & 0x01) && (flags & 0x04)) ) return 1;
 		}
 
-		p = find_panel(&image->det, x+xp, y+yp);
+		p = find_panel(image->det, x+xp, y+yp);
 		if ( p == NULL ) return 1;
 
 		/* Area of one pixel */
diff --git a/src/powder_plot.c b/src/powder_plot.c
index 3f0c1abd..3c8872f6 100644
--- a/src/powder_plot.c
+++ b/src/powder_plot.c
@@ -35,6 +35,7 @@ static void show_help(const char *s)
 "Compare intensity lists.\n"
 "\n"
 "  -h, --help                 Display this help message.\n"
+"  -g. --geometry=<file>   Get detector geometry from file.\n"
 "\n");
 }
 
@@ -46,16 +47,18 @@ int main(int argc, char *argv[])
 	int x, y;
 	struct hdfile *hdfile;
 	char *filename = NULL;
+	char *geometry = NULL;
 
 	/* Long options */
 	const struct option longopts[] = {
 		{"help",               0, NULL,               'h'},
 		{"input",              1, NULL,               'i'},
+		{"geometry",           1, NULL,               'g'},
 		{0, 0, NULL, 0}
 	};
 
 	/* Short options */
-	while ((c = getopt_long(argc, argv, "hi:", longopts, NULL)) != -1) {
+	while ((c = getopt_long(argc, argv, "hi:g:", longopts, NULL)) != -1) {
 
 		switch (c) {
 		case 'h' : {
@@ -72,6 +75,11 @@ int main(int argc, char *argv[])
 			break;
 		}
 
+		case 'g' : {
+			geometry = strdup(optarg);
+			break;
+		}
+
 		default : {
 			return 1;
 		}
@@ -84,7 +92,17 @@ int main(int argc, char *argv[])
 		return 1;
 	}
 
-	#include "geometry-lcls.tmp"
+	if ( geometry == NULL ) {
+		ERROR("You need to specify a geometry file with --geometry\n");
+		return 1;
+	}
+
+	image.det = get_detector_geometry(geometry);
+	if ( image.det == NULL ) {
+		ERROR("Failed to read detector geometry from '%s'\n", geometry);
+		return 1;
+	}
+	free(geometry);
 
 	hdfile = hdfile_open(filename);
 	hdfile_set_image(hdfile, "/data/data");