facetron: Add pre-refinement

7 files changed, 239 insertions, 61 deletions

diff --git a/src/Makefile.am b/src/Makefile.am
index be423094..34f8a9a0 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -55,7 +55,8 @@ calibrate_detector_SOURCES = calibrate_detector.c utils.c hdf5-file.c image.c \
                              sfac.c
 calibrate_detector_LDADD = @LIBS@
 
-facetron_SOURCES = facetron.c cell.c hdf5-file.c utils.c detector.c
+facetron_SOURCES = facetron.c cell.c hdf5-file.c utils.c detector.c peaks.c \
+                   image.c diffraction.c sfac.c
 facetron_LDADD = @LIBS@
 
 cubeit_SOURCES = cubeit.c cell.c hdf5-file.c utils.c
diff --git a/src/Makefile.in b/src/Makefile.in
index 4fae1697..4716d60e 100644
--- a/src/Makefile.in
+++ b/src/Makefile.in
@@ -69,7 +69,9 @@ am_cubeit_OBJECTS = cubeit.$(OBJEXT) cell.$(OBJEXT) \
 cubeit_OBJECTS = $(am_cubeit_OBJECTS)
 cubeit_DEPENDENCIES =
 am_facetron_OBJECTS = facetron.$(OBJEXT) cell.$(OBJEXT) \
-	hdf5-file.$(OBJEXT) utils.$(OBJEXT) detector.$(OBJEXT)
+	hdf5-file.$(OBJEXT) utils.$(OBJEXT) detector.$(OBJEXT) \
+	peaks.$(OBJEXT) image.$(OBJEXT) diffraction.$(OBJEXT) \
+	sfac.$(OBJEXT)
 facetron_OBJECTS = $(am_facetron_OBJECTS)
 facetron_DEPENDENCIES =
 am_get_hkl_OBJECTS = get_hkl.$(OBJEXT) sfac.$(OBJEXT) cell.$(OBJEXT) \
@@ -282,7 +284,9 @@ calibrate_detector_SOURCES = calibrate_detector.c utils.c hdf5-file.c image.c \
                              sfac.c
 
 calibrate_detector_LDADD = @LIBS@
-facetron_SOURCES = facetron.c cell.c hdf5-file.c utils.c detector.c
+facetron_SOURCES = facetron.c cell.c hdf5-file.c utils.c detector.c peaks.c \
+                   image.c diffraction.c sfac.c
+
 facetron_LDADD = @LIBS@
 cubeit_SOURCES = cubeit.c cell.c hdf5-file.c utils.c
 cubeit_LDADD = @LIBS@
diff --git a/src/cell.c b/src/cell.c
index 79261249..cc252b70 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -179,6 +179,21 @@ UnitCell *cell_new_from_cell(UnitCell *orig)
 }
 
 
+void cell_set_reciprocal(UnitCell *cell,
+                        double asx, double asy, double asz,
+                        double bsx, double bsy, double bsz,
+                        double csx, double csy, double csz)
+{
+	if ( cell == NULL ) return;
+
+	cell->axs = asx;  cell->ays = asy;  cell->azs = asz;
+	cell->bxs = bsx;  cell->bys = bsy;  cell->bzs = bsz;
+	cell->cxs = csx;  cell->cys = csy;  cell->czs = csz;
+
+	cell->rep = CELL_REP_RECIP;
+}
+
+
 /************************* Getter helper functions ****************************/
 
 static int cell_crystallographic_to_cartesian(UnitCell *cell,
diff --git a/src/cell.h b/src/cell.h
index 113e6512..42ef2eb6 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -62,6 +62,10 @@ extern int cell_get_reciprocal(UnitCell *cell,
                                double *bsx, double *bsy, double *bsz,
                                double *csx, double *csy, double *csz);
 
+extern void cell_set_reciprocal(UnitCell *cell,
+                               double asx, double asy, double asz,
+                               double bsx, double bsy, double bsz,
+                               double csx, double csy, double csz);
 
 extern double resolution(UnitCell *cell,
                          signed int h, signed int k, signed int l);
diff --git a/src/facetron.c b/src/facetron.c
index 7a1eec53..0fa6b643 100644
--- a/src/facetron.c
+++ b/src/facetron.c
@@ -26,6 +26,7 @@
 #include "cell.h"
 #include "hdf5-file.h"
 #include "detector.h"
+#include "peaks.h"
 
 
 #define MAX_HITS (1024)
@@ -46,9 +47,78 @@ static void show_help(const char *s)
 }
 
 
-/* Predict peaks */
-static int find_intersections(struct image *image, UnitCell *cell,
-                              double divergence, double bandwidth)
+static UnitCell *read_orientation_matrix(FILE *fh)
+{
+	float u, v, w;
+	struct rvec as, bs, cs;
+	UnitCell *cell;
+	char line[1024];
+
+	if ( fgets(line, 1023, fh) == NULL ) return NULL;
+	if ( sscanf(line, "astar = %f %f %f", &u, &v, &w) != 3 ) {
+		ERROR("Couldn't read a-star\n");
+		return NULL;
+	}
+	as.u = u*1e9;  as.v = v*1e9;  as.w = w*1e9;
+	if ( fgets(line, 1023, fh) == NULL ) return NULL;
+	if ( sscanf(line, "bstar = %f %f %f", &u, &v, &w) != 3 ) {
+		ERROR("Couldn't read b-star\n");
+		return NULL;
+	}
+	bs.u = u*1e9;  bs.v = v*1e9;  bs.w = w*1e9;
+	if ( fgets(line, 1023, fh) == NULL ) return NULL;
+	if ( sscanf(line, "cstar = %f %f %f", &u, &v, &w) != 3 ) {
+		ERROR("Couldn't read c-star\n");
+		return NULL;
+	}
+	cs.u = u*1e9;  cs.v = v*1e9;  cs.w = w*1e9;
+	cell = cell_new_from_axes(as, bs, cs);
+
+	return cell;
+}
+
+
+static int find_chunk(FILE *fh, UnitCell **cell, char **filename)
+{
+	char line[1024];
+	char *rval = NULL;
+
+	do {
+
+		rval = fgets(line, 1023, fh);
+		if ( rval == NULL ) continue;
+
+		chomp(line);
+
+		if ( strncmp(line, "Reflections from indexing", 25) != 0 ) {
+			continue;
+		}
+
+		*filename = strdup(line+29);
+
+		/* Skip two lines (while checking for errors) */
+		rval = fgets(line, 1023, fh);
+		if ( rval == NULL ) continue;
+		rval = fgets(line, 1023, fh);
+		if ( rval == NULL ) continue;
+
+		*cell = read_orientation_matrix(fh);
+		if ( *cell == NULL ) {
+			STATUS("Got filename but no cell for %s\n", *filename);
+			continue;
+		}
+
+		return 0;
+
+	} while ( rval != NULL );
+
+	return 1;
+}
+
+
+static struct reflhit *find_intersections(struct image *image, UnitCell *cell,
+                                          double divergence, double bandwidth,
+                                          int *n, int output)
 {
 	double asx, asy, asz;
 	double bsx, bsy, bsz;
@@ -60,7 +130,10 @@ static int find_intersections(struct image *image, UnitCell *cell,
 	signed int h, k, l;
 
 	hits = malloc(sizeof(struct reflhit)*MAX_HITS);
-	if ( hits == NULL ) return 0;
+	if ( hits == NULL ) {
+		*n = 0;
+		return NULL;
+	}
 
 	cell_get_reciprocal(cell, &asx, &asy, &asz,
 	                          &bsx, &bsy, &bsz,
@@ -140,84 +213,155 @@ static int find_intersections(struct image *image, UnitCell *cell,
 		hits[np].k = k;
 		hits[np].l = l;
 		np++;
-		printf("%i %i %i 0.0 (at %f,%f) %i\n", h, k, l, xda, yda, p);
+
+		if ( output ) {
+			printf("%i %i %i 0.0 (at %f,%f)\n", h, k, l, xda, yda);
+		}
 
 	}
 	}
 	}
 
-	free(hits);  /* FIXME! */
-
-	return np;
+	*n = np;
+	return hits;
 }
 
 
-static UnitCell *read_orientation_matrix(FILE *fh)
+static double integrate_all(struct image *image, struct reflhit *hits, int n)
 {
-	float u, v, w;
-	struct rvec as, bs, cs;
-	UnitCell *cell;
-	char line[1024];
+	double itot = 0.0;
+	int i;
 
-	if ( fgets(line, 1023, fh) == NULL ) return NULL;
-	if ( sscanf(line, "astar = %f %f %f", &u, &v, &w) != 3 ) {
-		ERROR("Couldn't read a-star\n");
-		return NULL;
-	}
-	as.u = u*1e9;  as.v = v*1e9;  as.w = w*1e9;
-	if ( fgets(line, 1023, fh) == NULL ) return NULL;
-	if ( sscanf(line, "bstar = %f %f %f", &u, &v, &w) != 3 ) {
-		ERROR("Couldn't read b-star\n");
-		return NULL;
+	for ( i=0; i<n; i++ ) {
+
+		float x, y, intensity;
+
+		if ( integrate_peak(image, hits[i].x, hits[i].y, &x, &y,
+                                    &intensity, 0, 0) ) continue;
+
+		itot += intensity;
 	}
-	bs.u = u*1e9;  bs.v = v*1e9;  bs.w = w*1e9;
-	if ( fgets(line, 1023, fh) == NULL ) return NULL;
-	if ( sscanf(line, "cstar = %f %f %f", &u, &v, &w) != 3 ) {
-		ERROR("Couldn't read c-star\n");
-		return NULL;
+
+	return itot;
+}
+
+
+static void get_trial_cell(UnitCell *out, UnitCell *in, int i, double step)
+{
+	double asx, asy, asz;
+	double bsx, bsy, bsz;
+	double csx, csy, csz;
+
+	cell_get_reciprocal(in, &asx, &asy, &asz, &bsx, &bsy, &bsz,
+	                    &csx, &csy, &csz);
+
+	switch ( i ) {
+	case  0 : asx += step; break;
+	case  1 : asx -= step; break;
+	case  2 : asy += step; break;
+	case  3 : asy -= step; break;
+	case  4 : asz += step; break;
+	case  5 : asz -= step; break;
+	case  6 : bsx += step; break;
+	case  7 : bsx -= step; break;
+	case  8 : bsy += step; break;
+	case  9 : bsy -= step; break;
+	case 10 : bsz += step; break;
+	case 11 : bsz -= step; break;
+	case 12 : csx += step; break;
+	case 13 : csx -= step; break;
+	case 14 : csy += step; break;
+	case 15 : csy -= step; break;
+	case 16 : csz += step; break;
+	case 17 : csz -= step; break;
+	default : break;
 	}
-	cs.u = u*1e9;  cs.v = v*1e9;  cs.w = w*1e9;
-	cell = cell_new_from_axes(as, bs, cs);
 
-	return cell;
+	cell_set_reciprocal(out, asx, asy, asz, bsx, bsy, bsz, csx, csy, csz);
 }
 
 
-static int find_chunk(FILE *fh, UnitCell **cell, char **filename)
+static void try_refine(struct image *image, UnitCell *cell,
+                       double da, double dw, double step)
 {
-	char line[1024];
-	char *rval = NULL;
+	struct reflhit *hits;
+	int np;
+	double itot;
+	UnitCell *trial_cell;
+	int did_something;
+
+	trial_cell = cell_new();
+
+	hits = find_intersections(image, cell, da, dw, &np, 0);
+	itot = integrate_all(image, hits, np);
+	STATUS("%f\n", itot);
 
 	do {
 
-		rval = fgets(line, 1023, fh);
-		if ( rval == NULL ) continue;
+		int i;
 
-		chomp(line);
+		did_something = 0;
 
-		if ( strncmp(line, "Reflections from indexing", 25) != 0 ) {
-			continue;
-		}
+		for ( i=0; i<18; i++ ) {
 
-		*filename = strdup(line+29);
+			struct reflhit *trial_hits;
+			double trial_itot;
 
-		/* Skip two lines (while checking for errors) */
-		rval = fgets(line, 1023, fh);
-		if ( rval == NULL ) continue;
-		rval = fgets(line, 1023, fh);
-		if ( rval == NULL ) continue;
+			get_trial_cell(trial_cell, cell, i, step);
 
-		*cell = read_orientation_matrix(fh);
-		if ( *cell == NULL ) {
-			STATUS("Got filename but no cell for %s\n", *filename);
-			continue;
+			trial_hits = find_intersections(image, trial_cell,
+				                        da, dw, &np, 0);
+			trial_itot = integrate_all(image, hits, np);
+
+			if ( trial_itot > itot ) {
+
+				double asx, asy, asz;
+				double bsx, bsy, bsz;
+				double csx, csy, csz;
+
+				cell_get_reciprocal(trial_cell, &asx, &asy,
+				                    &asz, &bsx, &bsy, &bsz,
+				                    &csx, &csy, &csz);
+				cell_set_reciprocal(cell, asx, asy, asz, bsx,
+				                    bsy, bsz, csx, csy, csz);
+
+				itot = trial_itot;
+				free(hits);
+				hits = trial_hits;
+
+				did_something = 1;
+
+			} else {
+
+				free(trial_hits);
+
+			}
 		}
 
-		return 0;
+	} while ( did_something );
 
-	} while ( rval != NULL );
+	free(trial_cell);
+}
 
-	return 1;
+
+/* Predict peaks */
+static void pre_refine(struct image *image, UnitCell *cell,
+                       double *da, double *dw)
+{
+	double asx, asy, asz;
+	double bsx, bsy, bsz;
+	double csx, csy, csz;
+	double istep, step;
+
+	/* Start by changing parameters by 1% */
+	cell_get_reciprocal(cell, &asx, &asy,
+	                    &asz, &bsx, &bsy, &bsz,
+	                    &csx, &csy, &csz);
+	istep = (asx+asy+asz+bsx+bsy+bsz+csx+csy+csz)/900.0;
+
+	for ( step=istep; step>istep/100.0; step-=istep/100.0 ) {
+		try_refine(image, cell, *da, *dw, step);
+	}
 }
 
 
@@ -288,6 +432,8 @@ int main(int argc, char *argv[])
 
 		struct image image;
 		struct hdfile *hdfile;
+		double da, dw;
+		int np;
 
 		STATUS("Integrating intensities from '%s'\n", filename);
 
@@ -302,7 +448,13 @@ int main(int argc, char *argv[])
 		}
 
 		hdf5_read(hdfile, &image, 0);
-		find_intersections(&image, cell, 5.0e-3, 3.0/100.0);
+
+		da = 5.0e-3;     /* Initial divergence */
+		dw = 3.0/100.0;  /* Initial bandwidth */
+
+		pre_refine(&image, cell, &da, &dw);
+
+		find_intersections(&image, cell, da, dw, &np, 1);
 
 		hdfile_close(hdfile);
 		free(cell);
diff --git a/src/peaks.c b/src/peaks.c
index 5f0a8a27..383f4663 100644
--- a/src/peaks.c
+++ b/src/peaks.c
@@ -169,9 +169,9 @@ static int cull_peaks(struct image *image)
 
 /* 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, int do_sa)
+int integrate_peak(struct image *image, int xp, int yp,
+                   float *xc, float *yc, float *intensity,
+                   int do_polar, int do_sa)
 {
 	signed int x, y;
 	const int lim = INTEGRATION_RADIUS * INTEGRATION_RADIUS;
diff --git a/src/peaks.h b/src/peaks.h
index 2a24c5fa..f262a54e 100644
--- a/src/peaks.h
+++ b/src/peaks.h
@@ -25,5 +25,7 @@ extern void output_intensities(struct image *image, UnitCell *cell,
                                pthread_mutex_t *mutex, int polar, int sa);
 extern int peak_sanity_check(struct image *image, UnitCell *cell);
 extern int find_projected_peaks(struct image *image, UnitCell *cell);
-
+extern int integrate_peak(struct image *image, int xp, int yp,
+                          float *xc, float *yc, float *intensity,
+                          int do_polar, int do_sa);
 #endif	/* PEAKS_H */