aboutsummaryrefslogtreecommitdiff
path: root/src/predict-refine.c
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2015-04-16 17:17:42 +0200
committerThomas White <taw@physics.org>2015-04-20 15:50:41 +0200
commit55e593d3402861c127cbd96de9e36c85b45bf96f (patch)
tree092928ae04641f3e2664d95e2209fe33e8945f15 /src/predict-refine.c
parente502ab129889a8849807406cd5f80680f6eb259c (diff)
New pairing procedure for prediction refinement (and auto-R)
Diffstat (limited to 'src/predict-refine.c')
-rw-r--r--src/predict-refine.c119
1 files changed, 80 insertions, 39 deletions
diff --git a/src/predict-refine.c b/src/predict-refine.c
index 34805436..62e9c128 100644
--- a/src/predict-refine.c
+++ b/src/predict-refine.c
@@ -162,27 +162,32 @@ 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 = 0;
+ double ax, ay, az;
+ double bx, by, bz;
+ double cx, cy, cz;
+ RefList *all_reflist;
+
+ all_reflist = reflist_new();
+ cell_get_cartesian(crystal_get_cell(cr),
+ &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
+
+ /* First, create a RefList containing the most likely indices for each
+ * peak, with no exclusion criteria */
for ( i=0; i<image_feature_count(image->features); i++ ) {
struct imagefeature *f;
double h, k, l, hd, kd, ld;
+ Reflection *refl;
+ struct panel *p;
/* Assume all image "features" are genuine peaks */
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(crystal_get_cell(cr),
- &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
-
- /* Decimal and fractional Miller indices of nearest
- * reciprocal lattice point */
+ /* 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;
@@ -190,38 +195,74 @@ static int pair_peaks(struct image *image, Crystal *cr,
k = lrint(kd);
l = lrint(ld);
- /* Check distance */
- if ( (fabs(h - hd) < min_dist)
- && (fabs(k - kd) < min_dist)
- && (fabs(l - ld) < min_dist) )
- {
- Reflection *refl;
- struct panel *p;
-
- refl = reflection_new(h, k, l);
- if ( refl == NULL ) {
- ERROR("Failed to create reflection\n");
- return 0;
- }
+ refl = reflection_new(h, k, l);
+ if ( refl == NULL ) {
+ ERROR("Failed to create reflection\n");
+ return 0;
+ }
+
+ add_refl_to_list(refl, all_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].refl = refl;
+ rps[n].peak = f;
+ rps[n].panel = p;
+ n++;
+
+ }
+
+ /* Get the excitation errors and detector positions for the candidate
+ * reflections */
+ crystal_set_reflections(cr, all_reflist);
+ update_partialities(cr, PMODEL_SCSPHERE);
+
+ /* Pass over the peaks again, keeping only the ones which look like
+ * good pairings */
+ for ( i=0; i<n; i++ ) {
+
+ double p, rlow, rhigh;
+ double fs, ss, pd;
+ signed int h, k, l;
+ Reflection *refl = rps[i].refl;
+
+ get_indices(refl, &h, &k, &l);
+
+ /* Is the supposed reflection anywhere near Bragg? */
+ get_partial(refl, &rlow, &rhigh, &p);
+ if ( (rlow-rhigh)/2.0 > 0.02e9 ) {
+ STATUS("rejecting %i %i %i because exerr=%e nm^-1\n",
+ h, k, l, 1e9*(rlow-rhigh)/2.0);
+ continue;
+ }
+
+ /* Is the supposed reflection anywhere near the peak? */
+ get_detector_pos(refl, &fs, &ss);
+ pd = pow(fs - rps[i].peak->fs, 2.0)
+ + pow(ss - rps[i].peak->ss, 2.0);
+ if ( pd > 100.0 * 100.0 ) {
+ STATUS("rejecting %i %i %i because %f %f -> %f %f\n",
+ h, k, l, fs, ss, rps[i].peak->fs, rps[i].peak->ss);
+ continue;
+ }
- 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 = p;
- n_acc++;
+ rps[n_acc] = rps[i];
+ rps[n_acc].refl = reflection_new(h, k, l);
+ copy_data(rps[n_acc].refl, refl);
+ if ( reflist != NULL ) {
+ add_refl_to_list(rps[n_acc].refl, reflist);
}
+ n_acc++;
}
+ reflist_free(all_reflist);
return n_acc;
}