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.c593
1 files changed, 593 insertions, 0 deletions
diff --git a/libcrystfel/src/peaks.c b/libcrystfel/src/peaks.c
new file mode 100644
index 00000000..ad524c61
--- /dev/null
+++ b/libcrystfel/src/peaks.c
@@ -0,0 +1,593 @@
+/*
+ * peaks.c
+ *
+ * Peak search and other image analysis
+ *
+ * (c) 2006-2011 Thomas White <taw@physics.org>
+ * 2011 Andrew Martin
+ *
+ * Part of CrystFEL - crystallography with a FEL
+ *
+ */
+
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+#include <string.h>
+#include <assert.h>
+#include <gsl/gsl_statistics_int.h>
+#include <pthread.h>
+#include <fenv.h>
+
+#include "image.h"
+#include "utils.h"
+#include "peaks.h"
+#include "detector.h"
+#include "filters.h"
+#include "diffraction.h"
+#include "reflist-utils.h"
+#include "beam-parameters.h"
+
+
+/* How close a peak must be to an indexed position to be considered "close"
+ * for the purposes of double hit detection and sanity checking. */
+#define PEAK_CLOSE (30.0)
+
+/* How close a peak must be to an indexed position to be considered "close"
+ * for the purposes of integration. */
+#define PEAK_REALLY_CLOSE (10.0)
+
+/* Degree of polarisation of X-ray beam */
+#define POL (1.0)
+
+static int cull_peaks_in_panel(struct image *image, struct panel *p)
+{
+ int i, n;
+ int nelim = 0;
+
+ n = image_feature_count(image->features);
+
+ for ( i=0; i<n; i++ ) {
+
+ struct imagefeature *f;
+ int j, ncol;
+
+ f = image_get_feature(image->features, i);
+ if ( f == NULL ) continue;
+
+ if ( f->fs < p->min_fs ) continue;
+ if ( f->fs > p->max_fs ) continue;
+ if ( f->ss < p->min_ss ) continue;
+ if ( f->ss > p->max_ss ) continue;
+
+ /* How many peaks are in the same column? */
+ ncol = 0;
+ for ( j=0; j<n; j++ ) {
+
+ struct imagefeature *g;
+
+ if ( i==j ) continue;
+
+ g = image_get_feature(image->features, j);
+ if ( g == NULL ) continue;
+
+ if ( p->badrow == 'f' ) {
+ if ( fabs(f->ss - g->ss) < 2.0 ) ncol++;
+ } else if ( p->badrow == 's' ) {
+ if ( fabs(f->fs - g->fs) < 2.0 ) ncol++;
+ } /* else do nothing */
+
+ }
+
+ /* More than three? */
+ if ( ncol <= 3 ) continue;
+
+ /* Yes? Delete them all... */
+ nelim = 0;
+ for ( j=0; j<n; j++ ) {
+ struct imagefeature *g;
+ g = image_get_feature(image->features, j);
+ if ( g == NULL ) continue;
+ if ( p->badrow == 'f' ) {
+ if ( fabs(f->ss - g->ss) < 2.0 ) {
+ image_remove_feature(image->features,
+ j);
+ nelim++;
+ }
+ } else if ( p->badrow == 's' ) {
+ if ( fabs(f->fs - g->ss) < 2.0 ) {
+ image_remove_feature(image->features,
+ j);
+ nelim++;
+ }
+ } else {
+ ERROR("Invalid badrow direction.\n");
+ abort();
+ }
+
+ }
+
+ }
+
+ return nelim;
+}
+
+
+/* Post-processing of the peak list to remove noise */
+static int cull_peaks(struct image *image)
+{
+ int nelim = 0;
+ struct panel *p;
+ int i;
+
+ for ( i=0; i<image->det->n_panels; i++ ) {
+ p = &image->det->panels[i];
+ if ( p->badrow != '-' ) {
+ nelim += cull_peaks_in_panel(image, p);
+ }
+ }
+
+ return nelim;
+}
+
+
+/* Returns non-zero if peak has been vetoed.
+ * i.e. don't use result if return value is not zero. */
+int integrate_peak(struct image *image, int cfs, int css,
+ double *pfs, double *pss, double *intensity,
+ double *pbg, double *pmax, double *sigma,
+ int do_polar, int centroid, int bgsub)
+{
+ signed int fs, ss;
+ double lim, out_lim, mid_lim;
+ double lim_sq, out_lim_sq, mid_lim_sq;
+ double total = 0.0;
+ double fsct = 0.0;
+ double ssct = 0.0;
+ double noise = 0.0;
+ int noise_counts = 0;
+ double max = 0.0;
+ struct panel *p = NULL;
+ int pixel_counts = 0;
+ double noise_mean = 0.0;
+ double noise_meansq = 0.0;
+ struct beam_params *beam;
+ double aduph;
+
+ beam = image->beam;
+ if ( beam != NULL ) {
+ aduph = image->beam->adu_per_photon;
+ } else {
+ aduph = 1.0;
+ }
+
+ p = find_panel(image->det, cfs, css);
+ if ( p == NULL ) return 1;
+ if ( p->no_index ) return 1;
+
+ lim = p->integr_radius;
+ mid_lim = 3.0 + lim;
+ out_lim = 6.0 + lim;
+ lim_sq = pow(lim, 2.0);
+ mid_lim_sq = pow(mid_lim, 2.0);
+ out_lim_sq = pow(out_lim, 2.0);
+
+ for ( fs=-out_lim; fs<+out_lim; fs++ ) {
+ for ( ss=-out_lim; ss<+out_lim; ss++ ) {
+
+ double val;
+ double tt = 0.0;
+ double phi, pa, pb, pol;
+ uint16_t flags;
+ struct panel *p2;
+ int idx;
+
+ /* Outer mask radius */
+ if ( fs*fs + ss*ss > out_lim_sq ) continue;
+
+ if ( ((fs+cfs)>=image->width) || ((fs+cfs)<0) ) continue;
+ if ( ((ss+css)>=image->height) || ((ss+css)<0) ) continue;
+
+ /* Strayed off one panel? */
+ p2 = find_panel(image->det, fs+cfs, ss+css);
+ if ( p2 != p ) return 1;
+
+ idx = fs+cfs+image->width*(ss+css);
+
+ /* Veto this peak if we tried to integrate in a bad region */
+ if ( image->flags != NULL ) {
+
+ flags = image->flags[idx];
+
+ /* It must have all the "good" bits to be valid */
+ if ( !((flags & image->det->mask_good)
+ == image->det->mask_good) ) return 1;
+
+ /* If it has any of the "bad" bits, reject */
+ if ( flags & image->det->mask_bad ) return 1;
+
+ }
+
+ val = image->data[idx];
+
+ if ( do_polar ) {
+
+ tt = get_tt(image, fs+cfs, ss+css);
+
+ phi = atan2(ss+css, fs+cfs);
+ pa = pow(sin(phi)*sin(tt), 2.0);
+ pb = pow(cos(tt), 2.0);
+ pol = 1.0 - 2.0*POL*(1-pa) + POL*(1.0+pb);
+
+ val /= pol;
+
+ }
+
+ if ( val > max ) max = val;
+
+ /* If outside inner mask, estimate noise from this region */
+ if ( fs*fs + ss*ss > mid_lim_sq ) {
+
+ /* Noise
+ * noise and noise_meansq are both in photons (^2) */
+ noise += val / image->beam->adu_per_photon;
+ noise_counts++;
+ noise_meansq += pow(val, 2.0);
+
+ } else if ( fs*fs + ss*ss < lim_sq ) {
+
+ /* Peak */
+ pixel_counts++;
+ total += val;
+ fsct += val*(cfs+fs);
+ ssct += val*(css+ss);
+
+ }
+
+ }
+ }
+
+ noise_mean = noise / noise_counts; /* photons */
+
+ /* The centroid is excitingly undefined if there is no intensity */
+ centroid = 0;
+ if ( centroid && (total != 0) ) {
+ *pfs = ((double)fsct / total) + 0.5;
+ *pss = ((double)ssct / total) + 0.5;
+ } else {
+ *pfs = (double)cfs + 0.5;
+ *pss = (double)css + 0.5;
+ }
+ if ( bgsub ) {
+ *intensity = total - aduph * pixel_counts*noise_mean; /* ADU */
+ } else {
+ *intensity = total; /* ADU */
+ }
+
+ if ( in_bad_region(image->det, *pfs, *pss) ) return 1;
+
+ if ( sigma != NULL ) {
+
+ /* First term is standard deviation of background per pixel
+ * sqrt(pixel_counts) - increase of error for integrated value
+ * sqrt(2) - increase of error for background subtraction */
+ *sigma = sqrt(noise_meansq/noise_counts-(noise_mean*noise_mean))
+ * sqrt(2.0*pixel_counts) * aduph;
+
+ }
+
+ if ( pbg != NULL ) {
+ *pbg = aduph * (noise / noise_counts);
+ }
+ if ( pmax != NULL ) {
+ *pmax = max;
+ }
+
+ return 0;
+}
+
+
+static void search_peaks_in_panel(struct image *image, float threshold,
+ float min_gradient, float min_snr,
+ struct panel *p)
+{
+ int fs, ss, stride;
+ float *data;
+ double d;
+ int idx;
+ double f_fs = 0.0;
+ double f_ss = 0.0;
+ double intensity = 0.0;
+ double sigma = 0.0;
+ double pbg = 0.0;
+ double pmax = 0.0;
+ int nrej_dis = 0;
+ int nrej_pro = 0;
+ int nrej_fra = 0;
+ int nrej_bad = 0;
+ int nrej_snr = 0;
+ int nacc = 0;
+ int ncull;
+ const int pws = p->peak_sep/2;
+
+ data = image->data;
+ stride = image->width;
+
+ for ( fs = p->min_fs+1; fs <= p->max_fs-1; fs++ ) {
+ for ( ss = p->min_ss+1; ss <= p->max_ss-1; ss++ ) {
+
+ double dx1, dx2, dy1, dy2;
+ double dxs, dys;
+ double grad;
+ int mask_fs, mask_ss;
+ int s_fs, s_ss;
+ double max;
+ unsigned int did_something;
+ int r;
+
+ /* Overall threshold */
+ if ( data[fs+stride*ss] < threshold ) continue;
+
+ /* Get gradients */
+ dx1 = data[fs+stride*ss] - data[(fs+1)+stride*ss];
+ dx2 = data[(fs-1)+stride*ss] - data[fs+stride*ss];
+ dy1 = data[fs+stride*ss] - data[(fs+1)+stride*(ss+1)];
+ dy2 = data[fs+stride*(ss-1)] - data[fs+stride*ss];
+
+ /* Average gradient measurements from both sides */
+ dxs = ((dx1*dx1) + (dx2*dx2)) / 2;
+ dys = ((dy1*dy1) + (dy2*dy2)) / 2;
+
+ /* Calculate overall gradient */
+ grad = dxs + dys;
+
+ if ( grad < min_gradient ) continue;
+
+ mask_fs = fs;
+ mask_ss = ss;
+
+ do {
+
+ max = data[mask_fs+stride*mask_ss];
+ did_something = 0;
+
+ for ( s_ss=biggest(mask_ss-pws/2,
+ p->min_ss);
+ s_ss<=smallest(mask_ss+pws/2,
+ p->max_ss);
+ s_ss++ ) {
+ for ( s_fs=biggest(mask_fs-pws/2,
+ p->min_fs);
+ s_fs<=smallest(mask_fs+pws/2,
+ p->max_fs);
+ s_fs++ ) {
+
+ if ( data[s_fs+stride*s_ss] > max ) {
+ max = data[s_fs+stride*s_ss];
+ mask_fs = s_fs;
+ mask_ss = s_ss;
+ did_something = 1;
+ }
+
+ }
+ }
+
+ /* Abort if drifted too far from the foot point */
+ if ( distance(mask_fs, mask_ss, fs, ss) >
+ p->peak_sep/2.0 )
+ {
+ break;
+ }
+
+ } while ( did_something );
+
+ /* Too far from foot point? */
+ if ( distance(mask_fs, mask_ss, fs, ss) > p->peak_sep/2.0 ) {
+ nrej_dis++;
+ continue;
+ }
+
+ /* Should be enforced by bounds used above. Muppet check. */
+ assert(mask_fs <= p->max_fs);
+ assert(mask_ss <= p->max_ss);
+ assert(mask_fs >= p->min_fs);
+ assert(mask_ss >= p->min_ss);
+
+ /* Centroid peak and get better coordinates.
+ * Don't bother doing polarisation/SA correction, because the
+ * intensity of this peak is only an estimate at this stage. */
+ r = integrate_peak(image, mask_fs, mask_ss,
+ &f_fs, &f_ss, &intensity,
+ &pbg, &pmax, &sigma, 0, 1, 1);
+
+ if ( r ) {
+ /* Bad region - don't detect peak */
+ nrej_bad++;
+ continue;
+ }
+
+ /* It is possible for the centroid to fall outside the image */
+ if ( (f_fs < p->min_fs) || (f_fs > p->max_fs)
+ || (f_ss < p->min_ss) || (f_ss > p->max_ss) ) {
+ nrej_fra++;
+ continue;
+ }
+
+ if (intensity/sigma < min_snr) {
+ nrej_snr++;
+ continue;
+ }
+
+ /* Check for a nearby feature */
+ image_feature_closest(image->features, f_fs, f_ss, &d, &idx);
+ if ( d < p->peak_sep/2.0 ) {
+ nrej_pro++;
+ continue;
+ }
+
+ /* Add using "better" coordinates */
+ image_add_feature(image->features, f_fs, f_ss, image, intensity,
+ NULL);
+ nacc++;
+
+ }
+ }
+
+ if ( image->det != NULL ) {
+ ncull = cull_peaks(image);
+ nacc -= ncull;
+ } else {
+ STATUS("Not culling peaks because I don't have a "
+ "detector geometry file.\n");
+ ncull = 0;
+ }
+
+// STATUS("%i accepted, %i box, %i proximity, %i outside panel, "
+// "%i in bad regions, %i with SNR < %g, %i badrow culled.\n",
+// nacc, nrej_dis, nrej_pro, nrej_fra, nrej_bad,
+// nrej_snr, min_snr, ncull);
+}
+
+
+void search_peaks(struct image *image, float threshold, float min_gradient,
+ float min_snr)
+{
+ int i;
+
+ if ( image->features != NULL ) {
+ image_feature_list_free(image->features);
+ }
+ image->features = image_feature_list_new();
+
+ for ( i=0; i<image->det->n_panels; i++ ) {
+
+ struct panel *p = &image->det->panels[i];
+
+ if ( p->no_index ) continue;
+ search_peaks_in_panel(image, threshold, min_gradient, min_snr, p);
+
+ }
+}
+
+
+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;
+
+ /* 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);
+
+ /* 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;
+
+ /* Assume all image "features" are genuine peaks */
+ f = image_get_feature(image->features, i);
+ if ( f == NULL ) continue;
+ n_feat++;
+
+ /* 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) )
+ {
+ n_sane++;
+ 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;
+
+ return 1;
+}
+
+
+/* Integrate the list of predicted reflections in "image" */
+void integrate_reflections(struct image *image, int polar, int use_closer,
+ int bgsub)
+{
+ Reflection *refl;
+ RefListIterator *iter;
+
+ for ( refl = first_refl(image->reflections, &iter);
+ refl != NULL;
+ refl = next_refl(refl, iter) ) {
+
+ double fs, ss, intensity;
+ double d;
+ int idx;
+ double bg, max;
+ double sigma;
+ double pfs, pss;
+ int r;
+
+ get_detector_pos(refl, &pfs, &pss);
+
+ /* Is there a really close feature which was detected? */
+ if ( use_closer ) {
+
+ struct imagefeature *f;
+
+ if ( image->features != NULL ) {
+ f = image_feature_closest(image->features,
+ pfs, pss, &d, &idx);
+ } else {
+ f = NULL;
+ }
+ if ( (f != NULL) && (d < PEAK_REALLY_CLOSE) ) {
+
+ pfs = f->fs;
+ pss = f->ss;
+
+ }
+ }
+
+ r = integrate_peak(image, pfs, pss, &fs, &ss,
+ &intensity, &bg, &max, &sigma, polar, 0,
+ bgsub);
+
+ /* Record intensity and set redundancy to 1 on success */
+ if ( r == 0 ) {
+ set_int(refl, intensity);
+ set_esd_intensity(refl, sigma);
+ set_redundancy(refl, 1);
+ } else {
+ set_redundancy(refl, 0);
+ }
+
+ }
+}
generated by cgit v1.2.3 (git 2.25.1) at 2024-07-29 06:48:19 +0000