More ReAx improvements

3 files changed, 162 insertions, 32 deletions

diff --git a/libcrystfel/src/peaks.c b/libcrystfel/src/peaks.c
index 0cbce55f..713dd0a6 100644
--- a/libcrystfel/src/peaks.c
+++ b/libcrystfel/src/peaks.c
@@ -473,9 +473,8 @@ void search_peaks(struct image *image, float threshold, float min_gradient,
 }
 
 
-int peak_sanity_check(struct image *image)
+double peak_lattice_agreement(struct image *image, UnitCell *cell, double *pst)
 {
-
 	int i;
 	int n_feat = 0;
 	int n_sane = 0;
@@ -483,14 +482,13 @@ int peak_sanity_check(struct image *image)
 	double bx, by, bz;
 	double cx, cy, cz;
 	double min_dist = 0.25;
+	double stot = 0.0;
 
 	/* Round towards nearest */
 	fesetround(1);
 
 	/* Cell basis vectors for this image */
-	cell_get_cartesian(image->indexed_cell, &ax, &ay, &az,
-	                                        &bx, &by, &bz,
-	                                        &cx, &cy, &cz);
+	cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
 
 	/* Loop over peaks, checking proximity to nearest reflection */
 	for ( i=0; i<image_feature_count(image->features); i++ ) {
@@ -520,17 +518,25 @@ int peak_sanity_check(struct image *image)
 		if ( (fabs(h - hd) < min_dist) && (fabs(k - kd) < min_dist)
 		  && (fabs(l - ld) < min_dist) )
 		{
+			double sval;
 			n_sane++;
+			sval = pow(h-hd, 2.0) + pow(k-kd, 2.0) + pow(l-ld, 2.0);
+			stot += 1.0 - sval;
 			continue;
 		}
 
 	}
 
-	/* return 0 means fail test, return 1 means pass test */
-	// printf("%d out of %d peaks are \"sane\"\n",n_sane,n_feat);
-	if ( (float)n_sane / (float)n_feat < 0.5 ) return 0;
+	*pst = stot;
+	return (double)n_sane / (float)n_feat;
+}
+
 
-	return 1;
+int peak_sanity_check(struct image *image)
+{
+	double stot;
+	/* 0 means failed test, 1 means passed test */
+	return peak_lattice_agreement(image, image->indexed_cell, &stot) >= 0.5;
 }
 
 
diff --git a/libcrystfel/src/peaks.h b/libcrystfel/src/peaks.h
index d52c75e4..ffeb896b 100644
--- a/libcrystfel/src/peaks.h
+++ b/libcrystfel/src/peaks.h
@@ -28,7 +28,10 @@ extern void integrate_reflections(struct image *image,
                                   int polar, int use_closer, int bgsub,
                                   double min_snr);
 
-extern int peak_sanity_check(struct image * image);
+extern double peak_lattice_agreement(struct image *image, UnitCell *cell,
+                                     double *pst);
+
+extern int peak_sanity_check(struct image *image);
 
 /* Exported so it can be poked by integration_check */
 extern int integrate_peak(struct image *image, int cfs, int css,
diff --git a/libcrystfel/src/reax.c b/libcrystfel/src/reax.c
index 80353a2c..ea8f143a 100644
--- a/libcrystfel/src/reax.c
+++ b/libcrystfel/src/reax.c
@@ -838,35 +838,97 @@ static int right_handed(struct rvec a, struct rvec b, struct rvec c)
 }
 
 
-static int check_twinning(UnitCell *c1, UnitCell *c2)
+struct cell_candidate
+{
+	UnitCell *cell;
+	double fom;
+};
+
+
+struct cell_candidate_list
+{
+	struct cell_candidate *cand;
+	int n_cand;
+};
+
+
+static int check_twinning(UnitCell *c1, UnitCell *c2, int verbose)
 {
 	int i;
+	int n_dup;
+	const int n_trials = 40;
 
 	double asx, asy, asz;
 	double bsx, bsy, bsz;
 	double csx, csy, csz;
+	double ax, ay, az;
+	double bx, by, bz;
+	double cx, cy, cz;
+
+	cell_get_reciprocal(c1, &asx, &asy, &asz,
+	                        &bsx, &bsy, &bsz,
+	                        &csx, &csy, &csz);
+	cell_get_cartesian(c2, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
+
+	n_dup = 0;
+	for ( i=0; i<n_trials; i++ ) {
 
-	for ( i=0; i<100; i++ ) {
 		signed int h, k, l;
+		double h2, k2, l2;
+		double rx, ry, rz;
+		double dev;
+		signed int h2i, k2i, l2i;
 
 		h = flat_noise(0, 10);
 		k = flat_noise(0, 10);
 		l = flat_noise(0, 10);
 
+		/* Position of this (randomly selected)
+		 * reciprocal lattice point */
+		rx = h*asx + k*bsx + l*csx;
+		ry = h*asy + k*bsy + l*csy;
+		rz = h*asz + k*bsz + l*csz;
+
+		/* Indices of this point in the basis of the other cell */
+		h2 = rx*ax + ry*ay + rz*az;
+		k2 = rx*bx + ry*by + rz*bz;
+		l2 = rx*cx + ry*cy + rz*cz;
+
+		h2i = lrint(h2);
+		k2i = lrint(k2);
+		l2i = lrint(l2);
+
+		dev = pow(h2i-h2, 2.0) + pow(k2i-k2, 2.0) + pow(l2i-l2, 2.0);
+
+		if ( verbose ) {
+			STATUS("%3i %3i %3i -> %5.2f %5.2f %5.2f -> "
+			       "%3i %3i %3i -> %5.2f\n", h, k, l,
+			       h2, k2, l2, h2i, k2i, l2i, dev);
+		}
+
+		if ( dev < 0.1 ) {
+			n_dup++;
+		}
+
+	}
+
+	if ( verbose ) {
+		STATUS("%i duplicates.\n", n_dup);
 	}
 
+	if ( n_dup > 10 ) return 1;
 	return 0;
 }
 
 
 /* Return true if "cnew" accounts for more than 25% of the peaks predicted by
  * any of the "ncells" cells in "cells". */
-static int twinned(UnitCell *cnew, UnitCell **cells, int ncells)
+static int twinned(UnitCell *cnew, struct cell_candidate_list *cl)
 {
 	int i;
 
-	for ( i=0; i<ncells; i++ ) {
-		if ( check_twinning(cnew, cells[i]) ) return 1;
+	for ( i=0; i<cl->n_cand; i++ ) {
+		if ( check_twinning(cnew, cl->cand[i].cell, 0) ) return 1;
 	}
 
 	return 0;
@@ -885,26 +947,64 @@ static int check_vector_combination(struct dvec *vi, struct dvec *vj,
 	ang = angle_between(vi->x, vi->y, vi->z, vj->x, vj->y, vj->z);
 	if ( fabs(ang-ga) > angtol ) return 0;
 
-	ang = angle_between(vi->x, vi->y, vi->z,
-	                    vk->x, vk->y, vk->z);
+	ang = angle_between(vi->x, vi->y, vi->z, vk->x, vk->y, vk->z);
 	if ( fabs(ang-be) > angtol ) return 0;
 
-	ang = angle_between(vj->x, vj->y, vj->z,
-	                    vk->x, vk->y, vk->z);
+	ang = angle_between(vj->x, vj->y, vj->z, vk->x, vk->y, vk->z);
 	if ( fabs(ang-al) > angtol ) return 0;
 
 	return 1;
 }
 
 
+static void add_cell_candidate(struct cell_candidate_list *cl, UnitCell *cnew,
+                               double fom)
+{
+	struct cell_candidate cshift;
+	int i, cpos;
+
+	cpos = cl->n_cand;
+	for ( i=0; i<cl->n_cand; i++ ) {
+		if ( fom > cl->cand[i].fom ) {
+			cpos = i;
+			break;
+		}
+	}
+
+	cshift.cell = cnew;
+	cshift.fom = fom;
+
+	for ( i=cpos; i<cl->n_cand; i++ ) {
+
+		struct cell_candidate cshift2;
+		cshift2 = cl->cand[i];
+		cl->cand[i] = cshift;
+		cshift = cshift2;
+
+	}
+
+	if ( cl->n_cand >= MAX_CELL_CANDIDATES ) {
+		/* "cshift" just fell off the end of the list */
+	} else {
+		cl->cand[cl->n_cand++] = cshift;
+	}
+}
+
+
 static void assemble_cells_from_candidates(struct image *image,
                                            struct reax_search *s,
                                            UnitCell *cell)
 {
 	int i, j, k;
 	signed int ti, tj, tk;
-	UnitCell *cells[MAX_CELL_CANDIDATES];
-	int ncells = 0;
+	struct cell_candidate_list cl;
+
+	cl.cand = calloc(MAX_CELL_CANDIDATES, sizeof(struct cell_candidate));
+	if ( cl.cand == NULL ) {
+		ERROR("Failed to allocate cell candidate list.\n");
+		return;
+	}
+	cl.n_cand = 0;
 
 	/* Find candidates for axes 0 and 1 which have the right angle */
 	for ( i=0; i<s->search[0].n_cand; i++ ) {
@@ -917,6 +1017,7 @@ static void assemble_cells_from_candidates(struct image *image,
 		struct dvec vi, vj, vk;
 		struct rvec ai, bi, ci;
 		UnitCell *cnew;
+		double fom;
 
 		vi = s->search[0].cand[i].v;
 		vj = s->search[1].cand[j].v;
@@ -937,17 +1038,15 @@ static void assemble_cells_from_candidates(struct image *image,
 		/* We have three vectors with the right angles */
 		cnew = cell_new_from_direct_axes(ai, bi, ci);
 
-		if ( twinned(cnew, cells, ncells) ) {
+		refine_cell(image, cnew, image->features);
+
+		if ( twinned(cnew, &cl) ) {
 			cell_free(cnew);
 			continue;
 		}
 
-		refine_cell(image, cnew, image->features);
-
-		/* FIXME: Rank according to quality of prediction */
-		if ( ncells < MAX_CELL_CANDIDATES ) {
-			cells[ncells++] = cnew;
-		}
+		peak_lattice_agreement(image, cnew, &fom);
+		add_cell_candidate(&cl, cnew, fom);
 
 	}
 	}
@@ -956,11 +1055,33 @@ static void assemble_cells_from_candidates(struct image *image,
 	}
 	}
 
-	image->ncells = ncells;
-	assert(ncells <= MAX_CELL_CANDIDATES);
-	for ( i=0; i<ncells; i++ ) {
-		image->candidate_cells[i] = cells[i];
+	for ( i=0; i<cl.n_cand; i++ ) {
+		double a, b, c, al, be, ga;
+		double aA, bA, cA, alA, beA, gaA;
+		int w = 0;
+//		STATUS("%i: %f\n", i, cl.cand[i].fom);
+		cell_get_parameters(cl.cand[i].cell, &a, &b, &c, &al, &be, &ga);
+		cell_get_parameters(cl.cand[i].cell, &aA, &bA, &cA,
+		                                     &alA, &beA, &gaA);
+		if ( (a - aA) > aA/10.0 ) w = 1;
+		if ( (b - bA) > bA/10.0 ) w = 1;
+		if ( (c - cA) > cA/10.0 ) w = 1;
+		if ( (al - alA) > deg2rad(5.0) ) w = 1;
+		if ( (be - beA) > deg2rad(5.0) ) w = 1;
+		if ( (ga - gaA) > deg2rad(5.0) ) w = 1;
+		if ( w ) {
+			STATUS("This cell is a long way from that sought:\n");
+			cell_print(cl.cand[i].cell);
+		}
 	}
+
+	image->ncells = cl.n_cand;
+	assert(image->ncells <= MAX_CELL_CANDIDATES);
+	for ( i=0; i<cl.n_cand; i++ ) {
+		image->candidate_cells[i] = cl.cand[i].cell;
+	}
+
+	free(cl.cand);
 }