refine_{radius,prediction}: unite common code (pair_peaks), avoid preliminary integration round

3 files changed, 117 insertions, 177 deletions

diff --git a/src/predict-refine.c b/src/predict-refine.c
index a2a4ce53..52496198 100644
--- a/src/predict-refine.c
+++ b/src/predict-refine.c
@@ -55,28 +55,77 @@ struct reflpeak {
                                * (we assume this never changes) */
 };
 
-static int pair_peaks(ImageFeatureList *flist, UnitCell *cell, RefList *reflist,
-                      struct reflpeak *rps, struct detector *det)
+
+static void twod_mapping(double fs, double ss, double *px, double *py,
+                         struct panel *p)
+{
+	double xs, ys;
+
+	xs = fs*p->fsx + ss*p->ssx;
+	ys = fs*p->fsy + ss*p->ssy;
+
+	*px = (xs + p->cnx) / p->res;
+	*py = (ys + p->cny) / p->res;
+}
+
+
+static double r_dev(struct reflpeak *rp)
+{
+	/* Excitation error term */
+	double rlow, rhigh, p;
+	get_partial(rp->refl, &rlow, &rhigh, &p);
+	return (rlow+rhigh)/2.0;
+}
+
+
+static double x_dev(struct reflpeak *rp, struct detector *det)
+{
+	/* Peak position term */
+	double xpk, ypk, xh, yh;
+	double fsh, ssh;
+	twod_mapping(rp->peak->fs, rp->peak->ss, &xpk, &ypk, rp->panel);
+	get_detector_pos(rp->refl, &fsh, &ssh);
+	twod_mapping(fsh, ssh, &xh, &yh, rp->panel);
+	return xh-xpk;
+}
+
+
+static double y_dev(struct reflpeak *rp, struct detector *det)
+{
+	/* Peak position term */
+	double xpk, ypk, xh, yh;
+	double fsh, ssh;
+	twod_mapping(rp->peak->fs, rp->peak->ss, &xpk, &ypk, rp->panel);
+	get_detector_pos(rp->refl, &fsh, &ssh);
+	twod_mapping(fsh, ssh, &xh, &yh, rp->panel);
+	return yh-ypk;
+}
+
+
+/* Associate a Reflection with each peak in "image" which is close to Bragg.
+ * Reflections will be added to "reflist", which can be NULL if this is not
+ * needed.  "rps" must be an array of sufficient size for all the peaks */
+static int pair_peaks(struct image *image, Crystal *cr,
+                      RefList *reflist, struct reflpeak *rps)
 {
 	int i;
 	const double min_dist = 0.05;
 	int n_acc = 0;
-	int n_notintegrated = 0;
 
-	for ( i=0; i<image_feature_count(flist); i++ ) {
+	for ( i=0; i<image_feature_count(image->features); i++ ) {
 
 		struct imagefeature *f;
 		double h, k, l, hd, kd, ld;
 
 		/* Assume all image "features" are genuine peaks */
-		f = image_get_feature(flist, i);
+		f = image_get_feature(image->features, i);
 		if ( f == NULL ) continue;
 
 		double ax, ay, az;
 		double bx, by, bz;
 		double cx, cy, cz;
 
-		cell_get_cartesian(cell,
+		cell_get_cartesian(crystal_get_cell(cr),
 		                   &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
 
 		/* Decimal and fractional Miller indices of nearest
@@ -94,17 +143,28 @@ static int pair_peaks(ImageFeatureList *flist, UnitCell *cell, RefList *reflist,
 		  && (fabs(l - ld) < min_dist) )
 		{
 			Reflection *refl;
+			struct panel *p;
 
-			/* Dig out the reflection */
-			refl = find_refl(reflist, h, k, l);
+			refl = reflection_new(h, k, l);
 			if ( refl == NULL ) {
-				n_notintegrated++;
-				continue;
+				ERROR("Failed to create reflection\n");
+				return 0;
 			}
 
+			if ( reflist != NULL ) add_refl_to_list(refl, reflist);
+			set_symmetric_indices(refl, h, k, l);
+
+			/* It doesn't matter if the actual predicted location
+			 * doesn't fall on this panel.  We're only interested
+			 * in how far away it is from the peak location.
+			 * The predicted position and excitation errors will be
+			 * filled in by update_partialities(). */
+			p = find_panel(image->det, f->fs, f->ss);
+			set_panel(refl, p);
+
 			rps[n_acc].refl = refl;
 			rps[n_acc].peak = f;
-			rps[n_acc].panel = find_panel(det, f->fs, f->ss);
+			rps[n_acc].panel = p;
 			n_acc++;
 		}
 
@@ -114,19 +174,48 @@ static int pair_peaks(ImageFeatureList *flist, UnitCell *cell, RefList *reflist,
 }
 
 
-static void twod_mapping(double fs, double ss, double *px, double *py,
-                         struct panel *p)
+static int cmpd2(const void *av, const void *bv)
 {
-	double xs, ys;
+	struct reflpeak *a, *b;
 
-	xs = fs*p->fsx + ss*p->ssx;
-	ys = fs*p->fsy + ss*p->ssy;
+	a = (struct reflpeak *)av;
+	b = (struct reflpeak *)bv;
 
-	*px = (xs + p->cnx) / p->res;
-	*py = (ys + p->cny) / p->res;
+	if ( fabs(r_dev(a)) < fabs(r_dev(b)) ) return -1;
+	return 1;
 }
 
 
+void refine_radius(Crystal *cr, struct image *image)
+{
+	int n, n_acc;
+	struct reflpeak *rps;
+	RefList *reflist;
+
+	/* Maximum possible size */
+	rps = malloc(image_feature_count(image->features)
+	                  * sizeof(struct reflpeak));
+	if ( rps == NULL ) return;
+
+	reflist = reflist_new();
+	n_acc = pair_peaks(image, cr, reflist, rps);
+	if ( n_acc < 3 ) {
+		ERROR("Too few paired peaks (%i) to determine radius\n", n_acc);
+		free(rps);
+		return;
+	}
+	crystal_set_reflections(cr, reflist);
+	update_partialities(cr, PMODEL_SCSPHERE);
+	crystal_set_reflections(cr, NULL);
+
+	qsort(rps, n_acc, sizeof(struct reflpeak), cmpd2);
+	n = (n_acc-1) - n_acc/50;
+	if ( n < 2 ) n = 2; /* n_acc is always >= 2 */
+	crystal_set_profile_radius(cr, fabs(r_dev(&rps[n])));
+
+	reflist_free(reflist);
+	free(rps);
+}
 
 
 /* Returns d(xh-xpk)/dP + d(yh-ypk)/dP, where P = any parameter */
@@ -239,39 +328,6 @@ static double y_gradient(int param, struct reflpeak *rp, struct detector *det,
 }
 
 
-static double r_dev(struct reflpeak *rp)
-{
-	/* Excitation error term */
-	double rlow, rhigh, p;
-	get_partial(rp->refl, &rlow, &rhigh, &p);
-	return (rlow+rhigh)/2.0;
-}
-
-
-static double x_dev(struct reflpeak *rp, struct detector *det)
-{
-	/* Peak position term */
-	double xpk, ypk, xh, yh;
-	double fsh, ssh;
-	twod_mapping(rp->peak->fs, rp->peak->ss, &xpk, &ypk, rp->panel);
-	get_detector_pos(rp->refl, &fsh, &ssh);
-	twod_mapping(fsh, ssh, &xh, &yh, rp->panel);
-	return xh-xpk;
-}
-
-
-static double y_dev(struct reflpeak *rp, struct detector *det)
-{
-	/* Peak position term */
-	double xpk, ypk, xh, yh;
-	double fsh, ssh;
-	twod_mapping(rp->peak->fs, rp->peak->ss, &xpk, &ypk, rp->panel);
-	get_detector_pos(rp->refl, &fsh, &ssh);
-	twod_mapping(fsh, ssh, &xh, &yh, rp->panel);
-	return yh-ypk;
-}
-
-
 static int iterate(struct reflpeak *rps, int n, UnitCell *cell,
                    struct image *image)
 {
@@ -459,18 +515,20 @@ int refine_prediction(struct image *image, Crystal *cr)
 	int i;
 	struct reflpeak *rps;
 	double max_I;
+	RefList *reflist;
 
 	rps = malloc(image_feature_count(image->features)
 	                       * sizeof(struct reflpeak));
 	if ( rps == NULL ) return 1;
 
-	n = pair_peaks(image->features, crystal_get_cell(cr),
-	               crystal_get_reflections(cr), rps, image->det);
+	reflist = reflist_new();
+	n = pair_peaks(image, cr, reflist, rps);
 	if ( n < 10 ) {
 		ERROR("Too few paired peaks to refine orientation.\n");
 		free(rps);
 		return 1;
 	}
+	crystal_set_reflections(cr, reflist);
 
 	/* Normalise the intensities to max 1 */
 	max_I = -INFINITY;
@@ -493,6 +551,8 @@ int refine_prediction(struct image *image, Crystal *cr)
 		update_partialities(cr, PMODEL_SCSPHERE);
 	}
 
+	crystal_set_reflections(cr, NULL);
+	reflist_free(reflist);
 	free(rps);
 	return 0;
 }
diff --git a/src/predict-refine.h b/src/predict-refine.h
index 9742f9e2..c763d2ca 100644
--- a/src/predict-refine.h
+++ b/src/predict-refine.h
@@ -39,6 +39,7 @@
 struct image;
 
 extern int refine_prediction(struct image *image, Crystal *cr);
+extern void refine_radius(Crystal *cr, struct image *image);
 
 
 #endif	/* PREDICT_REFINE_H */
diff --git a/src/process_image.c b/src/process_image.c
index f40b8cff..cea685f8 100644
--- a/src/process_image.c
+++ b/src/process_image.c
@@ -53,118 +53,6 @@
 #include "predict-refine.h"
 
 
-static int cmpd2(const void *av, const void *bv)
-{
-	double a, b;
-
-	a = *(double *)av;
-	b = *(double *)bv;
-
-	if ( fabs(a) < fabs(b) ) return -1;
-	return 1;
-}
-
-
-static double *excitation_errors(UnitCell *cell, ImageFeatureList *flist,
-                                 RefList *reflist, int *pnacc)
-{
-	int i;
-	const double min_dist = 0.05;
-	double *acc;
-	int n_acc = 0;
-	int n_notintegrated = 0;
-	int max_acc = 1024;
-
-	acc = malloc(max_acc*sizeof(double));
-	if ( acc == NULL ) {
-		ERROR("Allocation failed when refining radius!\n");
-		return NULL;
-	}
-
-	for ( i=0; i<image_feature_count(flist); i++ ) {
-
-		struct imagefeature *f;
-		double h, k, l, hd, kd, ld;
-
-		/* Assume all image "features" are genuine peaks */
-		f = image_get_feature(flist, i);
-		if ( f == NULL ) continue;
-
-		double ax, ay, az;
-		double bx, by, bz;
-		double cx, cy, cz;
-
-		cell_get_cartesian(cell,
-		                   &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
-
-		/* Decimal and fractional Miller indices of nearest
-		 * reciprocal lattice point */
-		hd = f->rx * ax + f->ry * ay + f->rz * az;
-		kd = f->rx * bx + f->ry * by + f->rz * bz;
-		ld = f->rx * cx + f->ry * cy + f->rz * cz;
-		h = lrint(hd);
-		k = lrint(kd);
-		l = lrint(ld);
-
-		/* Check distance */
-		if ( (fabs(h - hd) < min_dist)
-		  && (fabs(k - kd) < min_dist)
-		  && (fabs(l - ld) < min_dist) )
-		{
-			double rlow, rhigh, p;
-			Reflection *refl;
-
-			/* Dig out the reflection */
-			refl = find_refl(reflist, h, k, l);
-			if ( refl == NULL ) {
-				n_notintegrated++;
-				continue;
-			}
-
-			get_partial(refl, &rlow, &rhigh, &p);
-			acc[n_acc++] = fabs(rlow+rhigh)/2.0;
-			if ( n_acc == max_acc ) {
-				max_acc += 1024;
-				acc = realloc(acc, max_acc*sizeof(double));
-				if ( acc == NULL ) {
-					ERROR("Allocation failed during"
-					      " estimate_resolution!\n");
-					return NULL;
-				}
-			}
-		}
-
-	}
-
-	if ( n_acc < 3 ) {
-		STATUS("WARNING: Too few peaks to estimate profile radius.\n");
-		return NULL;
-	}
-
-	*pnacc = n_acc;
-	return acc;
-}
-
-
-static void refine_radius(Crystal *cr, ImageFeatureList *flist)
-{
-	int n = 0;
-	int n_acc;
-	double *acc;
-
-	acc = excitation_errors(crystal_get_cell(cr), flist,
-	                        crystal_get_reflections(cr), &n_acc);
-	if ( acc == NULL ) return;
-
-	qsort(acc, n_acc, sizeof(double), cmpd2);
-	n = n_acc/50;
-	if ( n < 2 ) n = 2; /* n_acc is always >= 2 */
-	crystal_set_profile_radius(cr, acc[(n_acc-1)-n]);
-
-	free(acc);
-}
-
-
 void process_image(const struct index_args *iargs, struct pattern_args *pargs,
                    Stream *st, int cookie, const char *tmpdir, int results_pipe,
                    int serial)
@@ -311,15 +199,9 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs,
 		}
 	}
 
-	/* Preliminary integration, needed for refinement */
-	integrate_all_4(&image, iargs->int_meth, PMODEL_SCSPHERE,
-	                iargs->push_res,
-	                iargs->ir_inn, iargs->ir_mid, iargs->ir_out,
-	                INTDIAG_NONE, 0, 0, 0, results_pipe);
-
 	/* Measure R before refinement */
 	for ( i=0; i<image.n_crystals; i++ ) {
-		refine_radius(image.crystals[i], image.features);
+		refine_radius(image.crystals[i], &image);
 	}
 
 	/* Integrate all the crystals at once - need all the crystals so that
@@ -341,7 +223,7 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs,
 			/* Reset the profile radius and estimate again with
 			 * better geometry */
 			crystal_set_profile_radius(image.crystals[i], 0.02e9);
-			refine_radius(image.crystals[i], image.features);
+			refine_radius(image.crystals[i], &image);
 
 			new_R = crystal_get_profile_radius(image.crystals[i]);
 
@@ -360,10 +242,7 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs,
 
 	}
 
-	/* The final, definitive, integration */
-	for ( i=0; i<image.n_crystals; i++ ) {
-		reflist_free(crystal_get_reflections(image.crystals[i]));
-	}
+	/* Integrate! */
 	integrate_all_4(&image, iargs->int_meth, PMODEL_SCSPHERE,
 	                iargs->push_res,
 	                iargs->ir_inn, iargs->ir_mid, iargs->ir_out,