aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2015-04-23 17:51:50 +0200
committerThomas White <taw@physics.org>2015-04-28 14:47:19 +0200
commit9d5afa09933ad3660d6b342c5828e5893e045eed (patch)
treea395be00a0da05f1089329fbf8153f401f344a14
parent51fc27d06e7cde3e0e57415d06624a7be712aed2 (diff)
Add outlier check for paired peaks
-rw-r--r--src/predict-refine.c62
1 files changed, 50 insertions, 12 deletions
diff --git a/src/predict-refine.c b/src/predict-refine.c
index c0e2aad5..50fc7f8a 100644
--- a/src/predict-refine.c
+++ b/src/predict-refine.c
@@ -160,6 +160,52 @@ static void write_pairs(const char *filename, struct reflpeak *rps, int n,
}
+static int cmpd2(const void *av, const void *bv)
+{
+ struct reflpeak *a, *b;
+
+ a = (struct reflpeak *)av;
+ b = (struct reflpeak *)bv;
+
+ if ( fabs(r_dev(a)) < fabs(r_dev(b)) ) return -1;
+ return 1;
+}
+
+
+static int check_outlier_transition(struct reflpeak *rps, int n,
+ struct detector *det)
+{
+ int i;
+
+ if ( n < 3 ) return n;
+
+ qsort(rps, n, sizeof(struct reflpeak), cmpd2);
+ write_pairs("pairs-before-outlier.lst", rps, n, det);
+
+ for ( i=1; i<n-1; i++ ) {
+
+ int j;
+ double grad = fabs(r_dev(&rps[i])) / i;
+
+ STATUS("Trying gradient %i\n", i);
+ for ( j=i+1; j<n; j++ ) {
+ if ( fabs(r_dev(&rps[j])) < 0.001e9+grad*j ) {
+ STATUS("Exit %i\n", j);
+ break;
+ }
+ }
+ if ( j == n ) {
+ STATUS("Outlier transition found at position %i / %i\n",
+ i, n);
+ return i;
+ }
+ }
+
+ STATUS("No outlier transition found.\n");
+ return n;
+}
+
+
/* 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 */
@@ -269,19 +315,11 @@ static int pair_peaks(struct image *image, Crystal *cr,
}
reflist_free(all_reflist);
- return n_acc;
-}
-
+ /* Sort the pairings by excitation error and look for a transition
+ * between good pairings and outliers */
+ n_acc = check_outlier_transition(rps, n_acc, image->det);
-static int cmpd2(const void *av, const void *bv)
-{
- struct reflpeak *a, *b;
-
- a = (struct reflpeak *)av;
- b = (struct reflpeak *)bv;
-
- if ( fabs(r_dev(a)) < fabs(r_dev(b)) ) return -1;
- return 1;
+ return n_acc;
}