Merge branch 'master' of git.bitwiz.org.uk:crystfel

11 files changed, 229 insertions, 78 deletions

diff --git a/src/detector.c b/src/detector.c
index 78fb216f..0748b2ca 100644
--- a/src/detector.c
+++ b/src/detector.c
@@ -61,6 +61,7 @@ void record_image(struct image *image, int do_poisson)
 	double total_energy, energy_density;
 	double ph_per_e;
 	double area;
+	double max_tt = 0.0;
 
 	/* How many photons are scattered per electron? */
 	area = M_PI*pow(BEAM_RADIUS, 2.0);
@@ -127,9 +128,17 @@ void record_image(struct image *image, int do_poisson)
 			ERROR("Processed negative at %i,%i %f\n", x, y, counts);
 		}
 
+		if ( image->twotheta[x + image->width*y] > max_tt ) {
+			max_tt = image->twotheta[x + image->width*y];
+		}
+
 	}
 	progress_bar(x, image->width-1, "Post-processing");
 	}
+
+	STATUS("Max 2theta = %.2f deg, min d = %.2f nm (halve this to get the"
+	       " voxel size for a synthesis)\n",
+	        rad2deg(max_tt), (image->lambda/(2.0*sin(max_tt/2.0)))/1e-9);
 }
 
 
diff --git a/src/geometry-simple.tmp b/src/geometry-simple.tmp
new file mode 100644
index 00000000..3472ff8a
--- /dev/null
+++ b/src/geometry-simple.tmp
@@ -0,0 +1,23 @@
+/* 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/get_hkl.c b/src/get_hkl.c
index ae279a8e..735ae52a 100644
--- a/src/get_hkl.c
+++ b/src/get_hkl.c
@@ -41,7 +41,9 @@ static void show_help(const char *s)
 "      --zone-axis            Generate hk0 intensities only (and add\n"
 "                              Synth2D-style header.\n"
 "  -i, --intensities=<file>   Read intensities from file instead of\n"
-"                              calculating them from scratch.\n"
+"                              calculating them from scratch.  You might use\n"
+"                              this if you need to apply noise or twinning.\n"
+"  -p, --pdb=<file>           PDB file from which to get the structure.\n"
 );
 }
 
@@ -114,6 +116,7 @@ int main(int argc, char *argv[])
 	unsigned int *cts;
 	char *input = NULL;
 	signed int h, k, l;
+	char *filename = NULL;
 
 	/* Long options */
 	const struct option longopts[] = {
@@ -124,11 +127,12 @@ int main(int argc, char *argv[])
 		{"twin",               0, &config_twin,        1},
 		{"zone-axis",          0, &config_za,          1},
 		{"intensities",        1, NULL,               'i'},
+		{"pdb",                1, NULL,               'p'},
 		{0, 0, NULL, 0}
 	};
 
 	/* Short options */
-	while ((c = getopt_long(argc, argv, "ht:o:i:", longopts, NULL)) != -1) {
+	while ((c = getopt_long(argc, argv, "ht:o:i:p:", longopts, NULL)) != -1) {
 
 		switch (c) {
 		case 'h' : {
@@ -151,6 +155,11 @@ int main(int argc, char *argv[])
 			break;
 		}
 
+		case 'p' : {
+			filename = strdup(optarg);
+			break;
+		}
+
 		case 0 : {
 			break;
 		}
@@ -162,10 +171,14 @@ int main(int argc, char *argv[])
 
 	}
 
-	mol = load_molecule();
+	if ( filename == NULL ) {
+		filename = strdup("molecule.pdb");
+	}
+
+	mol = load_molecule(filename);
 	cts = new_list_count();
 	if ( input == NULL ) {
-		ideal_ref = get_reflections(mol, eV_to_J(1790.0), 1/(0.6e-9), cts);
+		ideal_ref = get_reflections(mol, eV_to_J(1790.0), 1/(0.05e-9), cts);
 	} else {
 		ideal_ref = read_reflections(input, cts);
 		free(input);
diff --git a/src/hdf5-file.c b/src/hdf5-file.c
index 5573267b..4272cbe6 100644
--- a/src/hdf5-file.c
+++ b/src/hdf5-file.c
@@ -113,6 +113,10 @@ static void cleanup(hid_t fh)
 		type = H5Iget_type(id);
 
 		if ( type == H5I_GROUP ) H5Gclose(id);
+		if ( type == H5I_DATASET ) H5Dclose(id);
+		if ( type == H5I_DATATYPE ) H5Tclose(id);
+		if ( type == H5I_DATASPACE ) H5Sclose(id);
+		if ( type == H5I_ATTR ) H5Aclose(id);
 
 	}
 }
@@ -318,6 +322,11 @@ int hdf5_read(struct hdfile *f, struct image *image)
 {
 	herr_t r;
 	float *buf;
+	uint16_t *flags;
+	hid_t mask_dh;
+
+	image->height = f->nx;
+	image->width = f->ny;
 
 	buf = malloc(sizeof(float)*f->nx*f->ny);
 
@@ -327,10 +336,23 @@ int hdf5_read(struct hdfile *f, struct image *image)
 		ERROR("Couldn't read data\n");
 		return 1;
 	}
-
 	image->data = buf;
-	image->height = f->nx;
-	image->width = f->ny;
+
+	mask_dh = H5Dopen(f->fh, "/processing/hitfinder/masks", H5P_DEFAULT);
+	if ( mask_dh <= 0 ) {
+		ERROR("Couldn't open flags\n");
+	} else {
+		flags = malloc(sizeof(uint16_t)*f->nx*f->ny);
+		r = H5Dread(mask_dh, H5T_NATIVE_UINT16, H5S_ALL, H5S_ALL,
+		            H5P_DEFAULT, flags);
+		if ( r < 0 ) {
+			ERROR("Couldn't read flags\n");
+			image->flags = NULL;
+		} else {
+			image->flags = flags;
+		}
+	}
+	H5Dclose(mask_dh);
 
 	/* Read wavelength from file */
 	image->lambda = get_wavelength(f);
diff --git a/src/image.h b/src/image.h
index db1b491e..70809d22 100644
--- a/src/image.h
+++ b/src/image.h
@@ -32,26 +32,26 @@
 /* Structure describing a feature in an image */
 struct imagefeature {
 
-	struct image			*parent;
-	double				x;
-	double				y;
-	double				intensity;
+	struct image                    *parent;
+	double                          x;
+	double                          y;
+	double                          intensity;
 
 	/* Partner for this feature (in another feature list) or NULL */
-	struct imagefeature_struct	*partner;
+	struct imagefeature_struct      *partner;
 
 	/* Distance between this feature and its partner, if any. */
-	double				partner_d;
+	double                          partner_d;
 
 	/* Reciprocal space coordinates (m^-1 of course) of this feature */
-	double				rx;
-	double				ry;
-	double				rz;
+	double                          rx;
+	double                          ry;
+	double                          rz;
 
 	/* Internal use only */
-	int				valid;
+	int                             valid;
 
-	const char			*name;
+	const char                      *name;
 };
 
 /* An opaque type representing a list of image features */
@@ -80,46 +80,47 @@ struct reflhit {
 /* Structure describing an image */
 struct image {
 
-	float			*data;
-	double			*twotheta;
-	UnitCell		*indexed_cell;
-	UnitCell		*candidate_cells[MAX_CELL_CANDIDATES];
-	int			ncells;
-	struct detector		det;
-	const char		*filename;
-	struct reflhit		*hits;
-	int			n_hits;
+	float                   *data;
+	uint16_t                *flags;
+	double                  *twotheta;
+	UnitCell                *indexed_cell;
+	UnitCell                *candidate_cells[MAX_CELL_CANDIDATES];
+	int                     ncells;
+	struct detector         det;
+	const char              *filename;
+	struct reflhit          *hits;
+	int                     n_hits;
 
-	int			id;
+	int                     id;
 
-	struct quaternion	orientation;
+	struct quaternion       orientation;
 
 	/* Wavelength must always be given */
-	double			lambda;		/* Wavelength in m */
+	double                  lambda;   /* Wavelength in m */
 
 	/* Incident intensity (if unknown, put 1.0) */
-	double			f0;
-	int			f0_available;
+	double                  f0;
+	int                     f0_available;
 
-	int			width;
-	int			height;
+	int                     width;
+	int                     height;
 
-	ImageFeatureList	*features;	/* "Experimental" features */
+	ImageFeatureList        *features;   /* "Experimental" features */
 
 	/* DirAx auto-indexing low-level stuff */
-	int			dirax_pty;
-	pid_t			dirax_pid;
-	char			*dirax_rbuffer;
-	int			dirax_rbufpos;
-	int			dirax_rbuflen;
+	int                     dirax_pty;
+	pid_t                   dirax_pid;
+	char                    *dirax_rbuffer;
+	int                     dirax_rbufpos;
+	int                     dirax_rbuflen;
 
 	/* DirAx auto-indexing high-level stuff */
-	int			dirax_step;
-	int			dirax_read_cell;
-	int			best_acl;
-	int			best_acl_nh;
-	int			acls_tried[MAX_CELL_CANDIDATES];
-	int			n_acls_tried;
+	int                     dirax_step;
+	int                     dirax_read_cell;
+	int                     best_acl;
+	int                     best_acl_nh;
+	int                     acls_tried[MAX_CELL_CANDIDATES];
+	int                     n_acls_tried;
 
 };
 
diff --git a/src/indexamajig.c b/src/indexamajig.c
index 6cf9cd15..c365048a 100644
--- a/src/indexamajig.c
+++ b/src/indexamajig.c
@@ -260,18 +260,24 @@ static void *process_image(void *pargsv)
 	image.hits = NULL;
 	image.n_hits = 0;
 
+	/* View head-on (unit cell is tilted) */
+	image.orientation.w = 1.0;
+	image.orientation.x = 0.0;
+	image.orientation.y = 0.0;
+	image.orientation.z = 0.0;
+
 	STATUS("Processing '%s'\n", image.filename);
 
 	result = malloc(sizeof(*result));
 	if ( result == NULL ) return NULL;
+	result->peaks_sane = 0;
+	result->hit = 0;
 
 	hdfile = hdfile_open(filename);
 	if ( hdfile == NULL ) {
-		result->hit = 0;
 		return result;
 	} else if ( hdfile_set_first_image(hdfile, "/") ) {
 		ERROR("Couldn't select path\n");
-		result->hit = 0;
 		return result;
 	}
 
@@ -291,17 +297,7 @@ static void *process_image(void *pargsv)
 	if ( config_noisefilter ) {
 		filter_noise(&image, data_for_measurement);
 	} else {
-
-		int x, y;
-
-		for ( x=0; x<image.width; x++ ) {
-		for ( y=0; y<image.height; y++ ) {
-			float val;
-			val = image.data[x+image.width*y];
-			data_for_measurement[x+image.width*y] = val;
-		}
-		}
-
+		memcpy(data_for_measurement, image.data, data_size);
 	}
 
 	/* Perform 'fine' peak search */
@@ -329,13 +325,11 @@ static void *process_image(void *pargsv)
 	}
 
 	/* No cell at this point?  Then we're done. */
-	result->peaks_sane = 0;
 	if ( image.indexed_cell == NULL ) goto done;
 
 	/* Sanity check */
 	if ( !peak_sanity_check(&image, image.indexed_cell) ) {
 		STATUS("Failed peak sanity check.\n");
-		result->peaks_sane = 0;
 		goto done;
 	} else {
 		result->peaks_sane = 1;
@@ -372,6 +366,7 @@ 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);
diff --git a/src/parameters-lcls.tmp b/src/parameters-lcls.tmp
index 7ae2af48..1fb35065 100644
--- a/src/parameters-lcls.tmp
+++ b/src/parameters-lcls.tmp
@@ -12,3 +12,6 @@
 
 /* Radius of X-ray beam */
 #define BEAM_RADIUS (3.0e-6 / 2.0)
+
+/* Photon energy in eV */
+#define PHOTON_ENERGY (8000.0)
diff --git a/src/pattern_sim.c b/src/pattern_sim.c
index e09ab72a..130b6681 100644
--- a/src/pattern_sim.c
+++ b/src/pattern_sim.c
@@ -31,6 +31,7 @@
 #include "peaks.h"
 #include "sfac.h"
 #include "reflections.h"
+#include "parameters-lcls.tmp"
 
 
 static void show_help(const char *s)
@@ -41,12 +42,21 @@ static void show_help(const char *s)
 "pulses of X-rays from a free electron laser.\n"
 "\n"
 " -h, --help                Display this help message.\n"
+" -p, --pdb=<file>          PDB file from which to get the unit cell.\n"
+"                            (The actual Bragg intensities come from the\n"
+"                            intensities file)\n"
 "     --simulation-details  Show technical details of the simulation.\n"
 "     --gpu                 Use the GPU to speed up the calculation.\n"
 "\n"
 "     --near-bragg          Output h,k,l,I near Bragg conditions.\n"
 " -n, --number=<N>          Generate N images.  Default 1.\n"
 "     --no-images           Do not output any HDF5 files.\n"
+" -o, --output=<filename>   Output HDF5 filename.  Default: sim.h5.\n"
+"                            If you choose to simulate more than one pattern,\n"
+"                            the filename given will be postfixed with a\n"
+"                            hyphen, the image number and '.h5'.  In this\n"
+"                            case, the default value is 'sim', such that the\n"
+"                            files produced are sim-1.h5, sim-2.h5 and so on.\n"
 " -r, --random-orientation  Use a randomly generated orientation\n"
 "                            (a new orientation will be used for each image).\n"
 "     --powder              Write a summed pattern of all images simulated by\n"
@@ -166,7 +176,9 @@ int main(int argc, char *argv[])
 	int config_nosfac = 0;
 	int config_gpu = 0;
 	int config_powder = 0;
+	char *filename = NULL;
 	char *grad_str = NULL;
+	char *outfile = NULL;
 	GradientMethod grad;
 	int ndone = 0;    /* Number of simulations done (images or not) */
 	int number = 1;   /* Number used for filename of image */
@@ -189,11 +201,13 @@ int main(int argc, char *argv[])
 		{"intensities",        1, NULL,               'i'},
 		{"powder",             0, &config_powder,      1},
 		{"gradients",          1, NULL,               'g'},
+		{"pdb",                1, NULL,               'p'},
+		{"output",             1, NULL,               'o'},
 		{0, 0, NULL, 0}
 	};
 
 	/* Short options */
-	while ((c = getopt_long(argc, argv, "hrn:i:g:", longopts, NULL)) != -1) {
+	while ((c = getopt_long(argc, argv, "hrn:i:g:p:o:", longopts, NULL)) != -1) {
 
 		switch (c) {
 		case 'h' : {
@@ -221,6 +235,16 @@ int main(int argc, char *argv[])
 			break;
 		}
 
+		case 'p' : {
+			filename = strdup(optarg);
+			break;
+		}
+
+		case 'o' : {
+			outfile = strdup(optarg);
+			break;
+		}
+
 		case 0 : {
 			break;
 		}
@@ -232,6 +256,18 @@ int main(int argc, char *argv[])
 
 	}
 
+	if ( filename == NULL ) {
+		filename = strdup("molecule.pdb");
+	}
+
+	if ( outfile == NULL ) {
+		if ( n_images == 1 ) {
+			outfile = strdup("sim.h5");
+		} else {
+			outfile = strdup("sim");
+		}
+	}
+
 	if ( config_simdetails ) {
 		show_details();
 		return 0;
@@ -281,8 +317,8 @@ int main(int argc, char *argv[])
 	/* Define image parameters */
 	image.width = 1024;
 	image.height = 1024;
-	image.lambda = ph_en_to_lambda(eV_to_J(1790.0));  /* Wavelength */
-	cell = load_cell_from_pdb("molecule.pdb");
+	image.lambda = ph_en_to_lambda(eV_to_J(PHOTON_ENERGY)); /* Wavelength */
+	cell = load_cell_from_pdb(filename);
 	image.filename = NULL;
 
 	#include "geometry-lcls.tmp"
@@ -372,7 +408,13 @@ int main(int argc, char *argv[])
 
 			char filename[1024];
 
-			snprintf(filename, 1023, "results/sim-%i.h5", number);
+			if ( n_images != 1 ) {
+				snprintf(filename, 1023, "%s-%i.h5",
+				         outfile, number);
+			} else {
+				strncpy(filename, outfile, 1023);
+			}
+
 			number++;
 
 			/* Write the output file */
diff --git a/src/peaks.c b/src/peaks.c
index 22e12949..c6b79d25 100644
--- a/src/peaks.c
+++ b/src/peaks.c
@@ -133,7 +133,9 @@ static void cull_peaks(struct image *image)
 }
 
 
-static void integrate_peak(struct image *image, int xp, int yp,
+/* Returns non-zero if peak has been vetoed.
+ * i.e. don't use result if return value is not zero. */
+static int integrate_peak(struct image *image, int xp, int yp,
                            float *xc, float *yc, float *intensity,
                            int do_polar)
 {
@@ -150,6 +152,7 @@ static void integrate_peak(struct image *image, int xp, int yp,
 		double val, sa, pix_area, Lsq, dsq, proj_area;
 		float tt;
 		double phi, pa, pb, pol;
+		uint16_t flags;
 
 		/* Circular mask */
 		if ( x*x + y*y > lim ) continue;
@@ -157,6 +160,12 @@ static void integrate_peak(struct image *image, int xp, int yp,
 		if ( ((x+xp)>=image->width) || ((x+xp)<0) ) continue;
 		if ( ((y+yp)>=image->height) || ((y+yp)<0) ) continue;
 
+		/* Veto this peak if we tried to integrate in a bad region */
+		if ( image->flags != NULL ) {
+			flags = image->flags[(x+xp)+image->width*(y+yp)];
+			if ( !((flags & 0x01) && (flags & 0x04)) ) return 1;
+		}
+
 		p = find_panel(&image->det, x+xp, y+yp);
 
 		/* Area of one pixel */
@@ -203,6 +212,8 @@ static void integrate_peak(struct image *image, int xp, int yp,
 		*yc = (float)yp;
 		*intensity = 0;
 	}
+
+	return 0;
 }
 
 
@@ -219,6 +230,7 @@ void search_peaks(struct image *image)
 	int nrej_hot = 0;
 	int nrej_pro = 0;
 	int nrej_fra = 0;
+	int nrej_bad = 0;
 	int nacc = 0;
 
 	data = image->data;
@@ -240,6 +252,7 @@ void search_peaks(struct image *image)
 		int sx, sy;
 		double max;
 		unsigned int did_something;
+		int r;
 
 		/* Overall threshold */
 		if ( data[x+width*y] < 800 ) continue;
@@ -313,7 +326,13 @@ void search_peaks(struct image *image)
 		/* Centroid peak and get better coordinates.
 		 * Don't bother doing polarisation correction, because the
 		 * intensity of this peak is only an estimate at this stage. */
-		integrate_peak(image, mask_x, mask_y, &fx, &fy, &intensity, 0);
+		r = integrate_peak(image, mask_x, mask_y,
+		                   &fx, &fy, &intensity, 0);
+		if ( r ) {
+			/* Bad region - don't detect peak */
+			nrej_bad++;
+			continue;
+		}
 
 		/* It is possible for the centroid to fall outside the image */
 		if ( (fx < 0.0) || (fx > image->width)
@@ -337,8 +356,9 @@ void search_peaks(struct image *image)
 	}
 	}
 
-	STATUS("%i accepted, %i box, %i hot, %i proximity, %i outside frame\n",
-	       nacc, nrej_dis, nrej_hot, nrej_pro, nrej_fra);
+	STATUS("%i accepted, %i box, %i hot, %i proximity, %i outside frame, "
+	       "%i in bad regions.\n",
+	       nacc, nrej_dis, nrej_hot, nrej_pro, nrej_fra, nrej_bad);
 
 	cull_peaks(image);
 }
@@ -492,6 +512,8 @@ void output_intensities(struct image *image, UnitCell *cell,
 	int n_found;
 	int n_indclose = 0;
 	int n_foundclose = 0;
+	int n_veto = 0;
+	int n_veto_second = 0;
 	double asx, asy, asz;
 	double bsx, bsy, bsz;
 	double csx, csy, csz;
@@ -543,17 +565,35 @@ void output_intensities(struct image *image, UnitCell *cell,
 		                          &d, &idx);
 		if ( (f != NULL) && (d < PEAK_REALLY_CLOSE) ) {
 
+			int r;
+
 			/* f->intensity was measured on the filtered pattern,
 			 * so instead re-integrate using old coordinates.
 			 * This will produce further revised coordinates. */
-			integrate_peak(image, f->x, f->y, &x, &y, &intensity,
-			               !unpolar);
+			r = integrate_peak(image, f->x, f->y, &x, &y,
+			                   &intensity, !unpolar);
+			if ( r ) {
+				/* The original peak (which also went through
+				 * integrate_peak(), but with the mangled
+				 * image data) would have been rejected if it
+				 * was in a bad region.  Integration of the same
+				 * peak included a bad region this time. */
+				n_veto_second++;
+				continue;
+			}
 			intensity = f->intensity;
 
 		} else {
 
-			integrate_peak(image, hits[i].x, hits[i].y,
-			               &x, &y, &intensity, !unpolar);
+			int r;
+
+			r = integrate_peak(image, hits[i].x, hits[i].y,
+			                   &x, &y, &intensity, !unpolar);
+			if ( r ) {
+				/* Plain old ordinary peak veto */
+				n_veto++;
+				continue;
+			}
 
 		}
 
@@ -592,9 +632,12 @@ void output_intensities(struct image *image, UnitCell *cell,
 
 	printf("Peak statistics: %i peaks found by the peak search out of "
 	       "%i were close to indexed positions. "
-	       "%i indexed positions out of %i were close to detected peaks\n",
+	       "%i indexed positions out of %i were close to detected peaks.\n",
 	       n_foundclose, n_found, n_indclose, n_hits);
-
+	printf("%i integrations using indexed locations were aborted because "
+	       "they hit one or more bad pixels.\n", n_veto);
+	printf("%i integrations using peak search locations were aborted "
+	       "because they hit one or more bad pixels.\n", n_veto_second);
 	/* Blank line at end */
 	printf("\n");
 
diff --git a/src/sfac.c b/src/sfac.c
index 62c1e76e..aa7b3bf6 100644
--- a/src/sfac.c
+++ b/src/sfac.c
@@ -311,7 +311,7 @@ static void centre_molecule(struct molecule *mol)
 
 /* Load PDB file into a memory format suitable for efficient(ish) structure
  * factor calculation */
-struct molecule *load_molecule()
+struct molecule *load_molecule(const char *filename)
 {
 	struct molecule *mol;
 	FILE *fh;
@@ -327,7 +327,7 @@ struct molecule *load_molecule()
 	mol->reflections = NULL;
 	mol->cell = NULL;
 
-	fh = fopen("molecule.pdb", "r");
+	fh = fopen(filename, "r");
 	if ( fh == NULL ) {
                 ERROR("Couldn't open PDB file\n");
                 return NULL;
diff --git a/src/sfac.h b/src/sfac.h
index f6b4874f..54104700 100644
--- a/src/sfac.h
+++ b/src/sfac.h
@@ -58,7 +58,7 @@ struct molecule
 /* This is so that the water background calculation can use it */
 extern double complex get_sfac(const char *n, double s, double en);
 
-extern struct molecule *load_molecule(void);
+extern struct molecule *load_molecule(const char *filename);
 extern void free_molecule(struct molecule *mol);
 extern double *get_reflections(struct molecule *mol, double en, double res,
                                unsigned int *counts);