aboutsummaryrefslogtreecommitdiff
path: root/libcrystfel/src/peaks.c
diff options
context:
space:
mode:
Diffstat (limited to 'libcrystfel/src/peaks.c')
-rw-r--r--libcrystfel/src/peaks.c199
1 files changed, 113 insertions, 86 deletions
diff --git a/libcrystfel/src/peaks.c b/libcrystfel/src/peaks.c
index dabc3b1a..c350388b 100644
--- a/libcrystfel/src/peaks.c
+++ b/libcrystfel/src/peaks.c
@@ -3,12 +3,12 @@
*
* Peak search and other image analysis
*
- * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY,
* a research centre of the Helmholtz Association.
* Copyright © 2012 Richard Kirian
*
* Authors:
- * 2010-2012 Thomas White <taw@physics.org>
+ * 2010-2013 Thomas White <taw@physics.org>
* 2012 Kenneth Beyerlein <kenneth.beyerlein@desy.de>
* 2011 Andrew Martin <andrew.martin@desy.de>
* 2011 Richard Kirian
@@ -147,24 +147,15 @@ static int cull_peaks(struct image *image)
}
-/* cfs, css relative to panel origin */
-static int *make_BgMask(struct image *image, struct panel *p,
- double ir_out, double ir_inn)
+static void add_crystal_to_mask(struct image *image, struct panel *p,
+ double ir_out, double ir_inn, int w, int h,
+ int *mask, Crystal *cr)
{
Reflection *refl;
RefListIterator *iter;
- int *mask;
- int w, h;
-
- w = p->max_fs - p->min_fs + 1;
- h = p->max_ss - p->min_ss + 1;
- mask = calloc(w*h, sizeof(int));
- if ( mask == NULL ) return NULL;
-
- if ( image->reflections == NULL ) return mask;
/* Loop over all reflections */
- for ( refl = first_refl(image->reflections, &iter);
+ for ( refl = first_refl(crystal_get_reflections(cr), &iter);
refl != NULL;
refl = next_refl(refl, iter) )
{
@@ -205,6 +196,28 @@ static int *make_BgMask(struct image *image, struct panel *p,
}
}
+}
+
+
+/* cfs, css relative to panel origin */
+static int *make_BgMask(struct image *image, struct panel *p,
+ double ir_out, double ir_inn)
+{
+ int *mask;
+ int w, h;
+ int i;
+
+ w = p->max_fs - p->min_fs + 1;
+ h = p->max_ss - p->min_ss + 1;
+ mask = calloc(w*h, sizeof(int));
+ if ( mask == NULL ) return NULL;
+
+ if ( image->crystals == NULL ) return mask;
+
+ for ( i=0; i<image->n_crystals; i++ ) {
+ add_crystal_to_mask(image, p, ir_inn, ir_out,
+ w, h, mask, image->crystals[i]);
+ }
return mask;
}
@@ -581,29 +594,19 @@ void search_peaks(struct image *image, float threshold, float min_gradient,
}
-double peak_lattice_agreement(struct image *image, UnitCell *cell, double *pst)
+int peak_sanity_check(struct image *image)
{
- int i;
int n_feat = 0;
int n_sane = 0;
- double ax, ay, az;
- 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(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
+ int i;
+ const double min_dist = 0.25;
- /* Loop over peaks, checking proximity to nearest reflection */
for ( i=0; i<image_feature_count(image->features); i++ ) {
struct imagefeature *f;
struct rvec q;
double h,k,l,hd,kd,ld;
+ int j;
/* Assume all image "features" are genuine peaks */
f = image_get_feature(image->features, i);
@@ -613,38 +616,42 @@ double peak_lattice_agreement(struct image *image, UnitCell *cell, double *pst)
/* Reciprocal space position of found peak */
q = get_q(image, f->fs, f->ss, NULL, 1.0/image->lambda);
- /* Decimal and fractional Miller indices of nearest
- * reciprocal lattice point */
- hd = q.u * ax + q.v * ay + q.w * az;
- kd = q.u * bx + q.v * by + q.w * bz;
- ld = q.u * cx + q.v * cy + q.w * 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 sval;
- n_sane++;
- sval = pow(h-hd, 2.0) + pow(k-kd, 2.0) + pow(l-ld, 2.0);
- stot += 1.0 - sval;
- continue;
- }
+ for ( j=0; j<image->n_crystals; j++ ) {
+
+ double ax, ay, az;
+ double bx, by, bz;
+ double cx, cy, cz;
+
+ cell_get_cartesian(crystal_get_cell(image->crystals[j]),
+ &ax, &ay, &az,
+ &bx, &by, &bz,
+ &cx, &cy, &cz);
+
+ /* Decimal and fractional Miller indices of nearest
+ * reciprocal lattice point */
+ hd = q.u * ax + q.v * ay + q.w * az;
+ kd = q.u * bx + q.v * by + q.w * bz;
+ ld = q.u * cx + q.v * cy + q.w * 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) )
+ {
+ n_sane++;
+ continue;
+ }
- }
+ }
- *pst = stot;
- return (double)n_sane / (float)n_feat;
-}
+ }
-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;
+ return ((double)n_sane / n_feat) >= 0.5;
}
@@ -705,43 +712,28 @@ static struct integr_ind *sort_reflections(RefList *list, UnitCell *cell,
}
-/* Integrate the list of predicted reflections in "image" */
-void integrate_reflections(struct image *image, int use_closer, int bgsub,
- double min_snr,
- double ir_inn, double ir_mid, double ir_out,
- int integrate_saturated)
+static void integrate_crystal(Crystal *cr, struct image *image, int use_closer,
+ int bgsub, double min_snr,
+ double ir_inn, double ir_mid, double ir_out,
+ int integrate_saturated, int **bgMasks)
{
+ RefList *reflections;
struct integr_ind *il;
int n, i;
double av = 0.0;
int first = 1;
- int **bgMasks;
+ int n_saturated = 0;
+
+ reflections = crystal_get_reflections(cr);
- if ( num_reflections(image->reflections) == 0 ) return;
+ if ( num_reflections(reflections) == 0 ) return;
- il = sort_reflections(image->reflections, image->indexed_cell, &n);
+ il = sort_reflections(reflections, crystal_get_cell(cr), &n);
if ( il == NULL ) {
ERROR("Couldn't sort reflections\n");
return;
}
- /* Make background masks for all panels */
- bgMasks = calloc(image->det->n_panels, sizeof(int *));
- if ( bgMasks == NULL ) {
- ERROR("Couldn't create list of background masks.\n");
- return;
- }
- for ( i=0; i<image->det->n_panels; i++ ) {
- int *mask;
- mask = make_BgMask(image, &image->det->panels[i],
- ir_out, ir_inn);
- if ( mask == NULL ) {
- ERROR("Couldn't create background mask.\n");
- return;
- }
- bgMasks[i] = mask;
- }
-
for ( i=0; i<n; i++ ) {
double fs, ss, intensity;
@@ -809,7 +801,7 @@ void integrate_reflections(struct image *image, int use_closer, int bgsub,
1, bgMasks[pnum], &saturated);
if ( saturated ) {
- image->n_saturated++;
+ n_saturated++;
if ( !integrate_saturated ) r = 1;
}
@@ -844,14 +836,49 @@ void integrate_reflections(struct image *image, int use_closer, int bgsub,
}
}
+ crystal_set_num_saturated_reflections(cr, n_saturated);
+ crystal_set_resolution_limit(cr, 0.0);
+
+ free(il);
+}
+
+
+/* Integrate the list of predicted reflections in "image" */
+void integrate_reflections(struct image *image, int use_closer, int bgsub,
+ double min_snr,
+ double ir_inn, double ir_mid, double ir_out,
+ int integrate_saturated)
+{
+ int i;
+ int **bgMasks;
+
+ /* Make background masks for all panels */
+ bgMasks = calloc(image->det->n_panels, sizeof(int *));
+ if ( bgMasks == NULL ) {
+ ERROR("Couldn't create list of background masks.\n");
+ return;
+ }
+ for ( i=0; i<image->det->n_panels; i++ ) {
+ int *mask;
+ mask = make_BgMask(image, &image->det->panels[i],
+ ir_out, ir_inn);
+ if ( mask == NULL ) {
+ ERROR("Couldn't create background mask.\n");
+ return;
+ }
+ bgMasks[i] = mask;
+ }
+
+ for ( i=0; i<image->n_crystals; i++ ) {
+ integrate_crystal(image->crystals[i], image, use_closer,
+ bgsub, min_snr, ir_inn, ir_mid, ir_out,
+ integrate_saturated, bgMasks);
+ }
+
for ( i=0; i<image->det->n_panels; i++ ) {
free(bgMasks[i]);
}
free(bgMasks);
-
- image->diffracting_resolution = 0.0;
-
- free(il);
}
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-15 08:26:19 +0000