facetron: Work on geometry

6 files changed, 178 insertions, 80 deletions

diff --git a/src/Makefile.am b/src/Makefile.am
index 29949f7c..e2062d15 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -59,7 +59,8 @@ calibrate_detector_SOURCES = calibrate_detector.c utils.c hdf5-file.c image.c \
 calibrate_detector_LDADD = @LIBS@
 
 facetron_SOURCES = facetron.c cell.c hdf5-file.c utils.c detector.c peaks.c \
-                   image.c geometry.c reflections.c stream.c thread-pool.c
+                   image.c geometry.c reflections.c stream.c thread-pool.c \
+                   beam-parameters.c
 facetron_LDADD = @LIBS@
 
 cubeit_SOURCES = cubeit.c cell.c hdf5-file.c utils.c detector.c render.c \
diff --git a/src/Makefile.in b/src/Makefile.in
index aa0519e8..54dff463 100644
--- a/src/Makefile.in
+++ b/src/Makefile.in
@@ -312,7 +312,8 @@ calibrate_detector_SOURCES = calibrate_detector.c utils.c hdf5-file.c image.c \
 
 calibrate_detector_LDADD = @LIBS@
 facetron_SOURCES = facetron.c cell.c hdf5-file.c utils.c detector.c peaks.c \
-                   image.c geometry.c reflections.c stream.c thread-pool.c
+                   image.c geometry.c reflections.c stream.c thread-pool.c \
+                   beam-parameters.c
 
 facetron_LDADD = @LIBS@
 cubeit_SOURCES = cubeit.c cell.c hdf5-file.c utils.c detector.c render.c \
diff --git a/src/facetron.c b/src/facetron.c
index 8d66d549..806ce54f 100644
--- a/src/facetron.c
+++ b/src/facetron.c
@@ -31,6 +31,7 @@
 #include "geometry.h"
 #include "peaks.h"
 #include "thread-pool.h"
+#include "beam-parameters.h"
 
 
 static void show_help(const char *s)
@@ -45,6 +46,9 @@ static void show_help(const char *s)
 "                              (must be a file, not e.g. stdin)\n"
 "  -o, --output=<filename>    Output filename.  Default: facetron.hkl.\n"
 "  -g. --geometry=<file>      Get detector geometry from file.\n"
+"  -b, --beam=<file>          Get beam parameters from file (provides initial\n"
+"                              values for parameters, and nominal wavelengths\n"
+"                              if no per-shot value is found in an HDF5 file.\n"
 "  -x, --prefix=<p>           Prefix filenames from input file with <p>.\n"
 "      --basename             Remove the directory parts of the filenames.\n"
 "      --no-check-prefix      Don't attempt to correct the --prefix.\n"
@@ -91,6 +95,7 @@ static void integrate_image(int mytask, void *tasks)
 	int j, n;
 	struct hdfile *hdfile;
 	struct image *image = pargs->image;
+	double nominal_photon_energy = pargs->image->beam->photon_energy;
 
 	hdfile = hdfile_open(image->filename);
 	if ( hdfile == NULL ) {
@@ -102,34 +107,31 @@ static void integrate_image(int mytask, void *tasks)
 		return;
 	}
 
-	if ( hdf5_read(hdfile, pargs->image, 0) ) {
+	/* FIXME: Nominal photon energy */
+	if ( hdf5_read(hdfile, pargs->image, 0, nominal_photon_energy) ) {
 		ERROR("Couldn't read '%s'\n", image->filename);
 		hdfile_close(hdfile);
 		return;
 	}
 
 	/* Figure out which spots should appear in this pattern */
-	spots = find_intersections(image, image->indexed_cell,
-	                           image->div, image->bw, &n, 0);
+	spots = find_intersections(image, image->indexed_cell, &n, 1);
 
 	/* For each reflection, estimate the partiality */
 	for ( j=0; j<n; j++ ) {
 
 		signed int h, k, l;
 		float i_partial;
-		double p;
 		float xc, yc;
+		float i_full_est;
 
 		h = spots[j].h;
 		k = spots[j].k;
 		l = spots[j].l;
 
-		/* Calculated partiality of this spot in this pattern */
-		p = partiality(image, h, k, l);
-
 		/* Don't attempt to use spots with very small
 		 * partialities, since it won't be accurate. */
-		if ( p < 0.1 ) continue;
+		if ( spots[j].p < 0.1 ) continue;
 
 		/* Actual measurement of this reflection from this
 		 * pattern? */
@@ -139,8 +141,10 @@ static void integrate_image(int mytask, void *tasks)
 			continue;
 		}
 
+		i_full_est = i_partial * spots[j].p;
+
 		pthread_mutex_lock(pargs->list_lock);
-		integrate_intensity(pargs->i_full, h, k, l, i_partial);
+		integrate_intensity(pargs->i_full, h, k, l, i_full_est);
 		integrate_count(pargs->cts, h, k, l, 1);
 		if ( !find_item(pargs->obs, h, k, l) ) {
 			add_item(pargs->obs, h, k, l);
@@ -249,6 +253,7 @@ int main(int argc, char *argv[])
 	int n_total_patterns;
 	struct image *images;
 	int n_iter = 10;
+	struct beam_params *beam = NULL;
 
 	/* Long options */
 	const struct option longopts[] = {
@@ -256,6 +261,7 @@ int main(int argc, char *argv[])
 		{"input",              1, NULL,               'i'},
 		{"output",             1, NULL,               'o'},
 		{"geometry",           1, NULL,               'g'},
+		{"beam",               1, NULL,               'b'},
 		{"prefix",             1, NULL,               'x'},
 		{"basename",           0, &config_basename,    1},
 		{"no-check-prefix",    0, &config_checkprefix, 0},
@@ -265,7 +271,7 @@ int main(int argc, char *argv[])
 	};
 
 	/* Short options */
-	while ((c = getopt_long(argc, argv, "hi:g:x:j:y:o:",
+	while ((c = getopt_long(argc, argv, "hi:g:x:j:y:o:b:",
 	                        longopts, NULL)) != -1)
 	{
 
@@ -302,6 +308,15 @@ int main(int argc, char *argv[])
 			n_iter = atoi(optarg);
 			break;
 
+		case 'b' :
+			beam = get_beam_parameters(optarg);
+			if ( beam == NULL ) {
+				ERROR("Failed to load beam parameters"
+				      " from '%s'\n", optarg);
+				return 1;
+			}
+			break;
+
 		case 0 :
 			break;
 
@@ -348,6 +363,11 @@ int main(int argc, char *argv[])
 	}
 	free(geomfile);
 
+	if ( beam == NULL ) {
+		ERROR("You must provide a beam parameters file.\n");
+		return 1;
+	}
+
 	/* Prepare for iteration */
 	i_full = new_list_intensity();
 	obs = new_items();
@@ -388,13 +408,14 @@ int main(int argc, char *argv[])
 		snprintf(fnamereal, 1023, "%s%s", prefix, filename);
 
 		images[i].filename = fnamereal;
-		images[i].div = 0.5e-3;
-		images[i].bw = 0.001;
+		images[i].div = beam->divergence;
+		images[i].bw = beam->bandwidth;
 		images[i].orientation.w = 1.0;
 		images[i].orientation.x = 0.0;
 		images[i].orientation.y = 0.0;
 		images[i].orientation.z = 0.0;
 		images[i].det = det;
+		images[i].beam = beam;
 
 		/* Muppet proofing */
 		images[i].data = NULL;
diff --git a/src/geometry.c b/src/geometry.c
index cca2d050..2b94e7ba 100644
--- a/src/geometry.c
+++ b/src/geometry.c
@@ -15,18 +15,93 @@
 
 
 #include <stdlib.h>
+#include <gsl/gsl_poly.h>
+#include <assert.h>
 
 #include "utils.h"
 #include "cell.h"
 #include "image.h"
 #include "peaks.h"
+#include "beam-parameters.h"
 
 
-#define MAX_CPEAKS (1024)
+#define MAX_CPEAKS (256 * 256)
+
+
+static signed int locate_peak(double x, double y, double z, double lambda,
+                              struct detector *det, double *xdap, double *ydap)
+{
+	int p;
+	signed int found = -1;
+	const double den = 1.0/lambda + z;
+
+	for ( p=0; p<det->n_panels; p++ ) {
+
+		double xd, yd, cl;
+		double xda, yda;
+
+		/* Camera length for this panel */
+		cl = det->panels[p].clen;
+
+		/* Coordinates of peak relative to central beam, in m */
+		xd = cl * x / den;
+		yd = cl * y / den;
+
+		/* Convert to pixels */
+		xd *= det->panels[p].res;
+		yd *= det->panels[p].res;
+
+		/* Add the coordinates of the central beam */
+		xda = xd + det->panels[p].cx;
+		yda = yd + det->panels[p].cy;
+
+		/* Now, is this on this panel? */
+		if ( xda < det->panels[p].min_x ) continue;
+		if ( xda > det->panels[p].max_x ) continue;
+		if ( yda < det->panels[p].min_y ) continue;
+		if ( yda > det->panels[p].max_y ) continue;
+
+		/* If peak appears on multiple panels, reject it */
+		if ( found != -1 ) return -1;
+
+		/* Woohoo! */
+		found = p;
+		*xdap = xda;
+		*ydap = yda;
+
+	}
+
+	return found;
+}
+
+
+static double excitation_error(double xl, double yl, double zl,
+                               double ds_sq, double lambda)
+{
+	double tt;
+	double a, b, c;
+	double r1, r2;
+	int n;
+
+	tt = angle_between(0.0, 0.0, 1.0,  xl, yl, zl+1.0/lambda);
+
+	a = 1.0;
+	b = - cos(tt) * 2.0/lambda;
+	c = pow(lambda, -2.0) - ds_sq;
+
+	/* FIXME: I don't trust GSL's quadratic solver */
+	n = gsl_poly_solve_quadratic(a, b, c, &r1, &r2);
+	assert(n == 2);
+
+	r1 -= 1.0/lambda;
+	r2 -= 1.0/lambda;
+
+	if ( r1 > r2 ) return r2;
+	return r1;
+}
 
 
 struct cpeak *find_intersections(struct image *image, UnitCell *cell,
-                                 double divergence, double bandwidth,
                                  int *n, int output)
 {
 	double asx, asy, asz;
@@ -37,6 +112,11 @@ struct cpeak *find_intersections(struct image *image, UnitCell *cell,
 	int hmax, kmax, lmax;
 	double mres;
 	signed int h, k, l;
+	double bandwidth = image->bw;
+	double divergence = image->div;
+	double lambda = image->lambda;
+	const double profile_cutoff = 0.05e9;  /* 0.1 nm^-1 */
+	double llow, lhigh;    /* Wavelength */
 
 	cpeaks = malloc(sizeof(struct cpeak)*MAX_CPEAKS);
 	if ( cpeaks == NULL ) {
@@ -48,103 +128,91 @@ struct cpeak *find_intersections(struct image *image, UnitCell *cell,
 	                          &bsx, &bsy, &bsz,
 	                          &csx, &csy, &csz);
 
+	/* FIXME: Account for left-handed indexing */
+	asz = -asz;  bsz = -bsz;  csz = -csz;
+
 	mres = 1.0 / 8.0e-10;  /* 8 Angstroms */
 	hmax = mres / modulus(asx, asy, asz);
 	kmax = mres / modulus(bsx, bsy, bsz);
 	lmax = mres / modulus(csx, csy, csz);
 
+	/* "low" gives the largest Ewald sphere,
+	 * "high" gives the smallest Ewald sphere. */
+	llow = lambda - lambda*bandwidth/2.0;
+	lhigh = lambda + lambda*bandwidth/2.0;
+
 	for ( h=-hmax; h<hmax; h++ ) {
 	for ( k=-kmax; k<kmax; k++ ) {
 	for ( l=-lmax; l<lmax; l++ ) {
 
 		double xl, yl, zl;
-		double ds_sq, dps_sq;
-		double delta, divfact;
-		double llow, lhigh;
-		double xd, yd, cl;
+		double ds, ds_sq;
+		double rlow, rhigh;    /* "Excitation error" */
+		signed int plow, phigh;
+		double xdalow, xdahigh;
+		double ydalow, ydahigh;
 		double xda, yda;
-		int p;
-		int found = 0;
 
+		/* Ignore central beam */
 		if ( (h==0) && (k==0) && (l==0) ) continue;
 
-		llow = image->lambda - image->lambda*bandwidth/2.0;
-		lhigh = image->lambda + image->lambda*bandwidth/2.0;
-
 		/* Get the coordinates of the reciprocal lattice point */
 		zl = h*asz + k*bsz + l*csz;
-		if ( zl < 0.0 ) continue;  /* Do this check very early */
+		if ( zl > 0.0 ) continue;  /* Throw out if it's "in front" */
 		xl = h*asx + k*bsx + l*csx;
 		yl = h*asy + k*bsy + l*csy;
 
+		/* Calculate reciprocal lattice point modulus (and square) */
 		ds_sq = modulus_squared(xl, yl, zl);  /* d*^2 */
-		delta = divergence/image->lambda;
-		dps_sq = ds_sq + pow(delta, 2.0);  /* d'*^2 */
-
-		/* In range? */
-		divfact = 2.0 * delta * sqrt(xl*xl + yl*yl);
-		if ( ds_sq - 2.0*zl/llow > 0.0 ) continue;
-		if ( ds_sq - 2.0*zl/lhigh < 0.0 ) continue;
-
-		/* Work out which panel this peak would fall on */
-		for ( p=0; p<image->det->n_panels; p++ ) {
-
-			/* Camera length for this panel */
-			cl = image->det->panels[p].clen;
-
-			/* Coordinates of peak relative to central beam, in m */
-			xd = cl*xl / (ds_sq/(2.0*zl) - zl);
-			yd = cl*yl / (ds_sq/(2.0*zl) - zl);
-
-			/* Convert to pixels */
-			xd *= image->det->panels[p].res;
-			yd *= image->det->panels[p].res;
-
-			/* Add the coordinates of the central beam */
-			xda = xd + image->det->panels[p].cx;
-			yda = yd + image->det->panels[p].cy;
-
-			/* Now, is this on this panel? */
-			if ( xda < image->det->panels[p].min_x ) continue;
-			if ( xda > image->det->panels[p].max_x ) continue;
-			if ( yda < image->det->panels[p].min_y ) continue;
-			if ( yda > image->det->panels[p].max_y ) continue;
-
-			/* Woohoo! */
-			found = 1;
-			break;
-
+		ds = sqrt(ds_sq);
+		if ( ds > mres ) continue;  /* Outside resolution range */
+
+		/* Calculate excitation errors */
+		rlow = excitation_error(xl, yl, zl, ds_sq, llow);
+		rhigh = excitation_error(xl, yl, zl, ds_sq, lhigh);
+
+		if ( (fabs(rlow) > profile_cutoff)
+		  && (fabs(rhigh) > profile_cutoff) ) {
+			/* Reflection is not close to Bragg at either of
+			 * the two extremes of wavelength, so skip it. */
+			continue;
 		}
 
-		if ( !found ) continue;
-
+		/* Locate peak on detector, and check it doesn't span panels */
+		plow = locate_peak(xl, yl, zl, llow, image->det,
+		                   &xdalow, &ydalow);
+		if ( plow == -1 ) continue;
+		phigh = locate_peak(xl, yl, zl, lhigh, image->det,
+		                    &xdahigh, &ydahigh);
+		if ( phigh == -1 ) continue;
+		if ( plow != phigh ) continue;
+
+		xda = xdahigh;
+		yda = ydahigh;
 		cpeaks[np].h = h;
 		cpeaks[np].k = k;
 		cpeaks[np].l = l;
-		cpeaks[np].x = xda;
-		cpeaks[np].y = yda;
+		cpeaks[np].x = xdahigh;
+		cpeaks[np].y = ydahigh;
 		np++;
 
 		if ( output ) {
-			printf("%i %i %i 0.0 (at %f,%f)\n", h, k, l, xda, yda);
+			printf("%3i %3i %3i %6f (at %5.2f,%5.2f) %9e %9e\n",
+			       h, k, l, 0.0, xda, yda, rlow, rhigh);
 		}
 
+		if ( np == MAX_CPEAKS ) goto out;
+
 	}
 	}
 	}
 
+out:
 	*n = np;
 	return cpeaks;
 }
 
 
-/* Return the partiality of this reflection in this image */
-double partiality(struct image *image, signed int h, signed int k, signed int l)
-{
-	return 1.0;
-}
-
-
 double integrate_all(struct image *image, struct cpeak *cpeaks, int n)
 {
 	double itot = 0.0;
diff --git a/src/geometry.h b/src/geometry.h
index 8ebfc1b2..0782a3e1 100644
--- a/src/geometry.h
+++ b/src/geometry.h
@@ -18,11 +18,8 @@
 #endif
 
 extern struct cpeak *find_intersections(struct image *image, UnitCell *cell,
-                                        double divergence, double bandwidth,
                                         int *n, int output);
 
-extern double partiality(struct image *image, signed int h,
-                         signed int k, signed int l);
 extern double integrate_all(struct image *image, struct cpeak *cpeaks, int n);
 
 
diff --git a/src/image.h b/src/image.h
index 0f6db034..f6cb8a58 100644
--- a/src/image.h
+++ b/src/image.h
@@ -53,11 +53,21 @@ typedef struct _imagefeaturelist ImageFeatureList;
 
 
 /* This structure represents a predicted peak in an image */
-struct cpeak {
+struct cpeak
+{
+	/* Indices */
 	signed int h;
 	signed int k;
 	signed int l;
+
 	double min_distance;
+
+	/* Partiality */
+	double r1;
+	double r2;
+	double p;
+
+	/* Location in image */
 	int x;
 	int y;
 };
@@ -73,7 +83,7 @@ struct image {
 	UnitCell                *candidate_cells[MAX_CELL_CANDIDATES];
 	int                     ncells;
 	struct detector         *det;
-	struct beam_params      *beam;
+	struct beam_params      *beam;  /* The nominal beam parameters */
 	char                    *filename;
 	struct cpeak            *cpeaks;
 	int                     n_cpeaks;
@@ -86,7 +96,7 @@ struct image {
 	double                  m;  /* Mosaicity in radians */
 
 
-	/* Information about the radiation */
+	/* Per-shot radiation values */
 	double                  lambda;        /* Wavelength in m */
 	double                  div;           /* Divergence in radians */
 	double                  bw;            /* Bandwidth as a fraction */