Add make_pixelmap to core CrystFEL

5 files changed, 417 insertions, 1 deletions

diff --git a/AUTHORS b/AUTHORS
index e27035d5..b00a2a08 100644
--- a/AUTHORS
+++ b/AUTHORS
@@ -110,3 +110,6 @@ Contributors
 * Yaroslav Gevorkov <yaroslav.gevorkov@desy.de>
   peakfinder9
   xgandalf
+
+* Omri Mor <omor1@asu.edu>
+  Some work on make_pixelmap
diff --git a/CMakeLists.txt b/CMakeLists.txt
index a3f15912..e3238209 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -344,6 +344,14 @@ target_include_directories(whirligig PRIVATE ${COMMON_INCLUDES})
 target_link_libraries(whirligig ${COMMON_LIBRARIES})
 list(APPEND CRYSTFEL_EXECUTABLES whirligig)
 
+# ----------------------------------------------------------------------
+# make_pixelmap
+
+set(MAKE_PIXELMAP_SOURCES src/make_pixelmap.c)
+add_executable(make_pixelmap ${MAKE_PIXELMAP_SOURCES})
+target_include_directories(make_pixelmap PRIVATE ${COMMON_INCLUDES})
+target_link_libraries(make_pixelmap ${COMMON_LIBRARIES} ${HDF5_C_LIBRARIES})
+list(APPEND CRYSTFEL_EXECUTABLES make_pixelmap)
 
 # ----------------------------------------------------------------------
 # Install targets
@@ -376,6 +384,7 @@ install(FILES
 	doc/man/process_hkl.1
 	doc/man/render_hkl.1
 	doc/man/whirligig.1
+	doc/man/make_pixelmap.1
 	DESTINATION ${CMAKE_INSTALL_MANDIR}/man1
 )
 
diff --git a/README b/README
index a744f30d..42835723 100644
--- a/README
+++ b/README
@@ -31,6 +31,7 @@ Authors:
   Helen Ginn  <helen@strubi.ox.ac.uk>
   Nicolas Riebesel <nicolas.riebesel@tuhh.de>
   Yaroslav Gevorkov <yaroslav.gevorkov@desy.de>
+  Omri Mor <omor1@asu.edu>
 
 CrystFEL is free software: you can redistribute it and/or modify it under the
 terms of the GNU General Public License as published by the Free Software
@@ -100,7 +101,9 @@ In addition, there is also:
 
 	- render_hkl, for turning reflection lists into pretty graphics.
 
-	- list_events, for creating event lists from multi-event files
+	- list_events, for creating event lists from multi-event files.
+
+	- make_pixelmap, for creating pixel maps for Cheetah, OnDA etc.
 
 There is also a folder full of scripts for achieving many related tasks.
 
diff --git a/doc/man/make_pixelmap.1 b/doc/man/make_pixelmap.1
new file mode 100644
index 00000000..22fba821
--- /dev/null
+++ b/doc/man/make_pixelmap.1
@@ -0,0 +1,63 @@
+.\"
+.\" make_pixelmap man page
+.\"
+.\" Copyright © 2015-2018 Deutsches Elektronen-Synchrotron DESY,
+.\"                       a research centre of the Helmholtz Association.
+.\"
+.\" Part of CrystFEL - crystallography with a FEL
+.\"
+
+.TH MAKE_PIXELMAP 1
+.SH NAME
+make_pixelmap \- create pixelmap from CrystFEL geometry file
+.SH SYNOPSIS
+.PP
+\fBmake_pixelmap \fR[options] \fIinput.geom
+.PP
+\fBmake_pixelmap --help\fI
+
+.SH DESCRIPTION
+\fBmake_pixelmap\fR takes a CrystFEL geometry file as input, and produces an HDF5 "pixel map" as used by Cheetah and OnDA.  The HDF5 file will contain three arrays, "/x", "/y" and "/z", containing the x, y and z coordinates for each pixel, respectively, in metres, in real space.  The file will also contain "/res" and "/coffset", containing the resolution of the detector (one over the pixel size in metres) and the default coffset value for the detector (or that for the first panel, if there is no default value set).
+.P
+\fBmake_pixelmap\fR currently assumes that the file layout is "slabby", i.e. all of the image data appears in one "slab", with the fs and ss coordinates of each pixel being unique across the entire detector.
+.P
+With \fB--badmap\fR, make_pixelmap will instead write a map of bad pixels according to the bad regions specified in the geometry file.
+
+.SH OPTIONS
+
+.IP \fB--badmap
+.PD
+Create a bad pixel mask instead of a pixel geometry map.  The values to use for good and bad pixels in this map are given by \fB--good-pixel\fR and \fB--bad-pixel\fR respectively.  Note that CrystFEL also has other ways of specifying bad regions, e.g. by providing a separate mask.  This program is only concerned with the regions specified using "bad_" in the geometry file.
+
+.PD 0
+.IP \fB--good-pixel=\fIn
+.IP \fB--bad-pixel=\fIn
+.PD
+Specify the values to use in the mask for good and bad pixels, respectively.  The defaults are \fB--good-pixel=1\fR and \fB--bad-pixel=0\fR.
+
+.PD 0
+.IP "\fB-o \fIoutput.h5\fR"
+.IP \fB--output=\fIoutput.h5\fR
+.PD
+Set the output filename.  The default is \fB--output=pixelmap.h5\fR.
+
+.SH AUTHOR
+This page was written by Thomas White.
+
+.SH REPORTING BUGS
+Report bugs to <taw@physics.org>, or visit <http://www.desy.de/~twhite/crystfel>.
+
+.SH COPYRIGHT AND DISCLAIMER
+Copyright © 2015-2018 Deutsches Elektronen-Synchrotron DESY, a research centre of the Helmholtz Association.
+.P
+make_pixelmap, and this manual, are part of CrystFEL.
+.P
+CrystFEL is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
+.P
+CrystFEL is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
+.P
+You should have received a copy of the GNU General Public License along with CrystFEL.  If not, see <http://www.gnu.org/licenses/>.
+
+.SH SEE ALSO
+.BR crystfel (7),
+.BR crystfel_geometry (5)
diff --git a/src/make_pixelmap.c b/src/make_pixelmap.c
new file mode 100644
index 00000000..dfd91907
--- /dev/null
+++ b/src/make_pixelmap.c
@@ -0,0 +1,338 @@
+/*
+ * make_pixelmap.c
+ *
+ * Create a pixel map file from a CrystFEL geometry description
+ *
+ * Copyright © 2012-2018 Deutsches Elektronen-Synchrotron DESY,
+ *                       a research centre of the Helmholtz Association.
+ *
+ * Authors:
+ *   2012-2018 Thomas White <taw@physics.org>
+ *   2016-2016 Omri Mor <omor1@asu.edu>
+ *
+ * This file is part of CrystFEL.
+ *
+ * CrystFEL is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * CrystFEL is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with CrystFEL.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#define _ISOC99_SOURCE
+#include <math.h>
+#include <stdarg.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <unistd.h>
+#include <getopt.h>
+#include <assert.h>
+
+#include "crystfel/detector.h"
+
+
+static void show_help(const char *s)
+{
+	printf("Syntax: %s [options] <input.geom>\n\n", s);
+	printf(
+"Create a pixel map file.\n"
+"\n"
+" -h, --help                Display this help message.\n"
+"\n"
+" -o, --output=<filename>   Output filename.  Default: pixelmap.h5.\n"
+"     --badmap              Generate bad pixel map instead of geometry\n"
+"     --good-pixel=<n>      Value for good pixels in bad map.  Default 1.\n"
+"     --bad-pixel=<n>       Value for bad pixels in bad map.  Default 0.\n"
+);
+}
+
+
+static void create_array(hid_t gh, const char *name, void *vals,
+                         hid_t type, int width, int height)
+{
+	hid_t dh, sh;
+	herr_t r;
+	hsize_t size[2];
+	hsize_t max_size[2];
+
+	/* Note the "swap" here, according to section 3.2.5,
+	 * "C versus Fortran Dataspaces", of the HDF5 user's guide. */
+	size[0] = height;
+	size[1] = width;
+	max_size[0] = height;
+	max_size[1] = width;
+	sh = H5Screate_simple(2, size, max_size);
+
+	dh = H5Dcreate2(gh, name, type, sh,
+	                H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+	if ( dh < 0 ) {
+		ERROR("Couldn't create dataset\n");
+		return;
+	}
+
+	r = H5Dwrite(dh, type, H5S_ALL, H5S_ALL, H5P_DEFAULT, vals);
+	if ( r < 0 ) {
+		ERROR("Couldn't write data\n");
+		H5Dclose(dh);
+		return;
+	}
+
+	H5Sclose(sh);
+	H5Dclose(dh);
+}
+
+
+static void create_scalar(hid_t gh, const char *name, float val)
+{
+	hid_t dh, sh;
+	herr_t r;
+
+	sh = H5Screate(H5S_SCALAR);
+
+	STATUS("%s: %e\n", name, val);
+
+	dh = H5Dcreate2(gh, name, H5T_NATIVE_FLOAT, sh,
+	                H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+	if ( dh < 0 ) {
+		ERROR("Couldn't create dataset\n");
+		return;
+	}
+
+	r = H5Dwrite(dh, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT, &val);
+	if ( r < 0 ) {
+		ERROR("Couldn't write data\n");
+		H5Dclose(dh);
+		return;
+	}
+
+	H5Sclose(sh);
+	H5Dclose(dh);
+}
+
+
+static void write_pixelmap_hdf5(const char *filename,
+                                float *x, float *y, float *z,
+                                int width, int height, float res, float coffset)
+{
+	hid_t fh;
+
+	fh = H5Fcreate(filename, H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
+	if ( fh < 0 ) {
+		ERROR("Couldn't create file: %s\n", filename);
+		return;
+	}
+
+	create_array(fh, "x", x, H5T_NATIVE_FLOAT, width, height);
+	create_array(fh, "y", y, H5T_NATIVE_FLOAT, width, height);
+	create_array(fh, "z", z, H5T_NATIVE_FLOAT, width, height);
+
+	create_scalar(fh, "res", res);
+	create_scalar(fh, "coffset", coffset);
+
+	H5Fclose(fh);
+}
+
+
+static void write_badmap_hdf5(const char *filename, uint16_t *b,
+                              int width, int height)
+{
+	hid_t fh, gh;
+
+	fh = H5Fcreate(filename, H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
+	if ( fh < 0 ) {
+		ERROR("Couldn't create file: %s\n", filename);
+		return;
+	}
+
+	gh = H5Gcreate(fh, "data", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+
+	create_array(gh, "/data/data", b, H5T_STD_U16LE, width, height);
+
+	H5Fclose(fh);
+}
+
+
+int main(int argc, char *argv[])
+{
+	int c;
+	char *input_file = NULL;
+	char *output_file = NULL;
+	struct detector *det = NULL;
+	int fs, ss, w, h;
+	float *x, *y, *z;
+	uint16_t *b;
+	int i;
+	float res, coffset;
+	int badmap = 0;
+	int good_pixel_val = 1;
+	int bad_pixel_val = 0;
+
+	/* Long options */
+	const struct option longopts[] = {
+		{"help",               0, NULL,               'h'},
+		{"output",             1, NULL,               'o'},
+		{"badmap",             0, &badmap,             1},
+		{"good-pixel",         1, NULL,              301},
+		{"bad-pixel",          1, NULL,              302},
+		{0, 0, NULL, 0}
+	};
+
+	/* Short options */
+	while ((c = getopt_long(argc, argv, "ho:",
+	                        longopts, NULL)) != -1) {
+
+		switch (c) {
+
+			case 'h' :
+			show_help(argv[0]);
+			return 0;
+
+			case 'o' :
+			output_file = strdup(optarg);
+			break;
+
+			case 301:
+			if (sscanf(optarg, "%d", &good_pixel_val) != 1)
+			{
+				ERROR("Invalid value for --good-pixel\n");
+				return 1;
+			}
+			break;
+
+			case 302:
+			if (sscanf(optarg, "%d", &bad_pixel_val) != 1)
+			{
+				ERROR("Invalid value for --bad-pixel\n");
+				return 1;
+			}
+			break;
+
+			case 0 :
+			break;
+
+			default :
+			ERROR("Unhandled option '%c'\n", c);
+			break;
+
+		}
+
+	}
+
+	if ( argc != (optind+1) ) {
+		ERROR("Please provide exactly one geometry file name.\n");
+		return 1;
+	}
+
+	input_file = strdup(argv[optind++]);
+
+	if ( output_file == NULL ) {
+		output_file = strdup("pixelmap.h5");
+	}
+
+	/* Load geometry */
+	det = get_detector_geometry(input_file, NULL);
+	if ( det == NULL ) {
+		ERROR("Failed to read geometry from '%s'\n", input_file);
+		return 1;
+	}
+	free(input_file);
+
+	/* 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;
+		}
+	}
+	w += 1;  h += 1;  /* Inclusive -> Exclusive */
+	STATUS("Data slab size: %i x %i\n", w, h);
+
+	x = malloc(w*h*sizeof(float));
+	y = malloc(w*h*sizeof(float));
+	z = malloc(w*h*sizeof(float));
+	b = malloc(w*h*sizeof(uint16_t));
+	if ( (x==NULL) || (y==NULL) || (z==NULL) || (b==NULL) ) {
+		ERROR("Failed to allocate memory.\n");
+		return 1;
+	}
+
+	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;
+			}
+
+		}
+
+		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;
+	}
+
+	if ( badmap ) {
+		write_badmap_hdf5(output_file, b, w, h);
+	} else {
+		write_pixelmap_hdf5(output_file, x, y, z, w, h, res, coffset);
+	}
+
+	return 0;
+}