diff options
| -rw-r--r-- | doc/man/indexamajig.1 | 5 | ||||
| -rw-r--r-- | libcrystfel/Makefile.am | 4 | ||||
| -rw-r--r-- | libcrystfel/src/geometry.c | 12 | ||||
| -rw-r--r-- | libcrystfel/src/integration.c | 672 | ||||
| -rw-r--r-- | libcrystfel/src/integration.h | 75 | ||||
| -rw-r--r-- | libcrystfel/src/peaks.c | 240 | ||||
| -rw-r--r-- | libcrystfel/src/peaks.h | 13 | ||||
| -rw-r--r-- | src/indexamajig.c | 33 | ||||
| -rw-r--r-- | src/process_image.c | 28 | ||||
| -rw-r--r-- | src/process_image.h | 7 |
10 files changed, 804 insertions, 285 deletions
diff --git a/doc/man/indexamajig.1 b/doc/man/indexamajig.1 index 738614ef..550d04e8 100644 --- a/doc/man/indexamajig.1 +++ b/doc/man/indexamajig.1 @@ -300,11 +300,6 @@ with old scripts because this option selects the behaviour which is now the default. .PD 0 -.IP \fB--no-bg-sub\fR -.PD -Don't subtract local background estimates from integrated intensities. This is almost never useful, but might help to detect problems from troublesome background noise. - -.PD 0 .IP \fB--use-saturated\fR .PD Normally, peaks which contain one or more pixels above max_adu (defined in the detector geometry file) will not be used for indexing. Using this option skips this check, possibly improving the indexing rate if there is a large proportion of saturated peaks. diff --git a/libcrystfel/Makefile.am b/libcrystfel/Makefile.am index a9bda55b..2812fa52 100644 --- a/libcrystfel/Makefile.am +++ b/libcrystfel/Makefile.am @@ -9,7 +9,7 @@ libcrystfel_la_SOURCES = src/reflist.c src/utils.c src/cell.c src/detector.c \ src/reflist-utils.c src/filters.c \ src/render.c src/index.c src/dirax.c src/mosflm.c \ src/cell-utils.c src/integer_matrix.c src/crystal.c \ - src/grainspotter.c src/xds.c + src/grainspotter.c src/xds.c src/integration.c if HAVE_FFTW libcrystfel_la_SOURCES += src/reax.c @@ -26,7 +26,7 @@ libcrystfel_la_include_HEADERS = src/beam-parameters.h src/hdf5-file.h \ src/filters.h src/dirax.h src/mosflm.h \ src/reax.h src/cell-utils.h \ src/integer_matrix.h src/crystal.h \ - src/grainspotter.h src/xds.h + src/grainspotter.h src/xds.h src/integration.h AM_CPPFLAGS = -DDATADIR=\""$(datadir)"\" -I$(top_builddir)/lib -Wall AM_CPPFLAGS += -I$(top_srcdir)/lib @LIBCRYSTFEL_CFLAGS@ diff --git a/libcrystfel/src/geometry.c b/libcrystfel/src/geometry.c index efb3f2ce..273ac7f0 100644 --- a/libcrystfel/src/geometry.c +++ b/libcrystfel/src/geometry.c @@ -135,11 +135,13 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst, double cet, cez; double pr; double L; + double del; /* Don't predict 000 */ if ( abs(h)+abs(k)+abs(l) == 0 ) return NULL; pr = crystal_get_profile_radius(cryst); + del = image->div + crystal_get_mosaicity(cryst); /* "low" gives the largest Ewald sphere (wavelength short => k large) * "high" gives the smallest Ewald sphere (wavelength long => k small) @@ -152,12 +154,12 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst, tl = sqrt(xl*xl + yl*yl); - cet = -sin(image->div/2.0) * khigh; - cez = -cos(image->div/2.0) * khigh; + cet = -sin(del/2.0) * khigh; + cez = -cos(del/2.0) * khigh; rhigh = khigh - distance(cet, cez, tl, zl); /* Loss of precision */ - cet = sin(image->div/2.0) * klow; - cez = -cos(image->div/2.0) * klow; + cet = sin(del/2.0) * klow; + cez = -cos(del/2.0) * klow; rlow = klow - distance(cet, cez, tl, zl); /* Loss of precision */ if ( (signbit(rlow) == signbit(rhigh)) @@ -168,7 +170,7 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst, ERROR("Reflection with rlow < rhigh!\n"); ERROR("%3i %3i %3i rlow = %e, rhigh = %e\n", h, k, l, rlow, rhigh); - ERROR("div = %e\n", image->div); + ERROR("div + m = %e\n", del); return NULL; } diff --git a/libcrystfel/src/integration.c b/libcrystfel/src/integration.c new file mode 100644 index 00000000..d6860fec --- /dev/null +++ b/libcrystfel/src/integration.c @@ -0,0 +1,672 @@ +/* + * integration.c + * + * Integration of intensities + * + * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. + * + * Authors: + * 2010-2013 Thomas White <taw@physics.org> + * + * This file is part of CrystFEL. + * + * CrystFEL is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * CrystFEL is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with CrystFEL. If not, see <http://www.gnu.org/licenses/>. + * + */ + + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include <stdlib.h> +#include <assert.h> +#include <gsl/gsl_matrix.h> +#include <gsl/gsl_vector.h> +#include <gsl/gsl_linalg.h> +#include <gsl/gsl_eigen.h> + +#include "reflist.h" +#include "cell.h" +#include "crystal.h" +#include "cell-utils.h" +#include "geometry.h" +#include "image.h" +#include "peaks.h" +#include "integration.h" + + +struct integr_ind +{ + double res; + Reflection *refl; +}; + + +static int compare_resolution(const void *av, const void *bv) +{ + const struct integr_ind *a = av; + const struct integr_ind *b = bv; + + return a->res > b->res; +} + + +static struct integr_ind *sort_reflections(RefList *list, UnitCell *cell, + int *np) +{ + struct integr_ind *il; + Reflection *refl; + RefListIterator *iter; + int i, n; + + n = num_reflections(list); + *np = 0; /* For now */ + + if ( n == 0 ) return NULL; + + il = calloc(n, sizeof(struct integr_ind)); + if ( il == NULL ) return NULL; + + i = 0; + for ( refl = first_refl(list, &iter); + refl != NULL; + refl = next_refl(refl, iter) ) + { + signed int h, k, l; + double res; + + if ( get_redundancy(refl) == 0 ) continue; + + get_indices(refl, &h, &k, &l); + res = resolution(cell, h, k, l); + + il[i].res = res; + il[i].refl = refl; + + i++; + assert(i <= n); + } + + qsort(il, n, sizeof(struct integr_ind), compare_resolution); + + *np = n; + return il; +} + + +static void check_eigen(gsl_vector *e_val) +{ + int i; + double vmax, vmin; + const int n = e_val->size; + const double max_condition = 1e6; + const int verbose = 0; + + if ( verbose ) STATUS("Eigenvalues:\n"); + vmin = +INFINITY; + vmax = 0.0; + for ( i=0; i<n; i++ ) { + double val = gsl_vector_get(e_val, i); + if ( verbose ) STATUS("%i: %e\n", i, val); + if ( val > vmax ) vmax = val; + if ( val < vmin ) vmin = val; + } + + for ( i=0; i<n; i++ ) { + double val = gsl_vector_get(e_val, i); + if ( val < vmax/max_condition ) { + gsl_vector_set(e_val, i, 0.0); + } + } + + vmin = +INFINITY; + vmax = 0.0; + for ( i=0; i<n; i++ ) { + double val = gsl_vector_get(e_val, i); + if ( val == 0.0 ) continue; + if ( val > vmax ) vmax = val; + if ( val < vmin ) vmin = val; + } + if ( verbose ) { + STATUS("Condition number: %e / %e = %5.2f\n", + vmax, vmin, vmax/vmin); + } +} + + +static gsl_vector *solve_svd(gsl_vector *v, gsl_matrix *M) +{ + gsl_matrix *s_vec; + gsl_vector *s_val; + int err, n; + gsl_vector *shifts; + + n = v->size; + if ( v->size != M->size1 ) return NULL; + if ( v->size != M->size2 ) return NULL; + + s_val = gsl_vector_calloc(n); + s_vec = gsl_matrix_calloc(n, n); + + err = gsl_linalg_SV_decomp_jacobi(M, s_vec, s_val); + if ( err ) { + ERROR("SVD failed: %s\n", gsl_strerror(err)); + gsl_matrix_free(s_vec); + gsl_vector_free(s_val); + return NULL; + } + /* "M" is now "U" */ + + check_eigen(s_val); + + shifts = gsl_vector_calloc(n); + err = gsl_linalg_SV_solve(M, s_vec, s_val, v, shifts); + if ( err ) { + ERROR("Matrix solution failed: %s\n", gsl_strerror(err)); + gsl_matrix_free(s_vec); + gsl_vector_free(s_val); + gsl_vector_free(shifts); + return NULL; + } + + gsl_matrix_free(s_vec); + gsl_vector_free(s_val); + + return shifts; +} + + +struct intcontext +{ + int halfw; + int w; + struct image *image; + gsl_matrix *bgm; +}; + + +static void addm(gsl_matrix *m, int i, int j, double val) +{ + double v = gsl_matrix_get(m, i, j); + gsl_matrix_set(m, i, j, v+val); +} + + +static void addv(gsl_vector *v, int i, double val) +{ + double k = gsl_vector_get(v, i); + gsl_vector_set(v, i, k+val); +} + + +static double boxi(struct intcontext *ic, + int fid_fs, int fid_ss, int p, int q) +{ + int fs, ss; + + fs = fid_fs + p; + ss = fid_ss + q; + + /* FIXME: Wrong panel? */ + return ic->image->data[fs + ic->image->width*ss]; +} + + +static void fit_bg(struct intcontext *ic, + int fid_fs, int fid_ss, + double *pa, double *pb, double *pc) +{ + int p, q; + gsl_vector *v; + gsl_vector *ans; + + v = gsl_vector_calloc(3); + + for ( p=0; p<ic->w; p++ ) { + for ( q=0; q<ic->w; q++ ) { + double bi; + bi = boxi(ic, fid_fs, fid_ss, p, q); + + addv(v, 0, bi*p); + addv(v, 1, bi*q); + addv(v, 2, bi); + } + } + + ans = solve_svd(v, ic->bgm); + gsl_vector_free(v); + + *pa = gsl_vector_get(ans, 0); + *pb = gsl_vector_get(ans, 1); + *pc = gsl_vector_get(ans, 2); + + gsl_vector_free(ans); +} + + +static void init_intcontext(struct intcontext *ic) +{ + int p, q; + + ic->w = 2*ic->halfw + 1; + + ic->bgm = gsl_matrix_calloc(3, 3); + if ( ic->bgm == NULL ) { + ERROR("Failed to initialise matrix.\n"); + return; + } + + for ( p=0; p<ic->w; p++ ) { + for ( q=0; q<ic->w; q++ ) { + addm(ic->bgm, 0, 0, p*p); + addm(ic->bgm, 0, 1, p*q); + addm(ic->bgm, 0, 2, p); + addm(ic->bgm, 1, 0, p*q); + addm(ic->bgm, 1, 1, q*q); + addm(ic->bgm, 1, 2, q); + addm(ic->bgm, 2, 0, p); + addm(ic->bgm, 2, 1, q); + addm(ic->bgm, 2, 2, 1); + } + } +} + + +static void measure_all_intensities(RefList *list, struct image *image) +{ + Reflection *refl; + RefListIterator *iter; + struct intcontext ic; + + ic.halfw = 4; + ic.image = image; + init_intcontext(&ic); + + for ( refl = first_refl(list, &iter); + refl != NULL; + refl = next_refl(refl, iter) ) + { + double pfs, pss; + int fid_fs, fid_ss; + double a, b, c; + + get_detector_pos(refl, &pfs, &pss); + fid_fs = lrint(pfs); + fid_ss = lrint(pss); + fit_bg(&ic, fid_fs, fid_ss, &a, &b, &c); + + //set_intensity(refl, intensity); + } +} + + +static void estimate_mosaicity(Crystal *cr, struct image *image) +{ + int msteps = 50; + int i; + const double mest = crystal_get_mosaicity(cr); + const double mmax = 2.0 * mest; + RefList *list; + + STATUS("Initial estimate: m = %f\n", mest); + + crystal_set_mosaicity(cr, mmax); + list = find_intersections(image, cr); + crystal_set_reflections(cr, list); + measure_all_intensities(list, image); + + for ( i=1; i<=msteps; i++ ) { + + /* "m" varies from just over zero up to 2x the given estimate */ + Reflection *refl; + RefListIterator *iter; + const double m = mmax*((double)i/msteps); + int n_gained = 0; + int n_lost = 0; + double i_gained = 0.0; + double i_lost = 0.0; + + crystal_set_mosaicity(cr, m); + update_partialities(cr, PMODEL_SPHERE); + + for ( refl = first_refl(list, &iter); + refl != NULL; + refl = next_refl(refl, iter) ) + { + if ( get_redundancy(refl) == 0 ) { + if ( get_temp1(refl) > 0.0 ) { + i_lost += get_intensity(refl); + n_lost++; + } + set_temp1(refl, -1.0); + } else if ( get_temp1(refl) < 0.0 ) { + i_gained += get_intensity(refl); + n_gained++; + set_temp1(refl, 1.0); + } + } + + if ( i > 1 ) { + STATUS("%.2e %10.2f %4i %10.2f %4i %10.2f\n", m, + i_gained, n_gained, i_lost, n_lost, + i_gained - i_lost); + } + + } +} + + +static void estimate_resolution(RefList *reflections, Crystal *cr, + struct image *image) +{ + struct integr_ind *il; + int n, i; + int score = 1000; /* FIXME */ + int cutoff = 0; + double limit = 0.0; + + if ( num_reflections(reflections) == 0 ) return; + + il = sort_reflections(reflections, crystal_get_cell(cr), &n); + if ( il == NULL ) { + ERROR("Couldn't sort reflections\n"); + return; + } + + for ( i=0; i<n; i++ ) { + + double intensity, sigma, snr; + Reflection *refl; + + refl = il[i].refl; + intensity = get_intensity(refl); + sigma = get_esd_intensity(refl); + + /* I/sigma(I) cutoff + * Rejects reflections below --min-integration-snr, or if the + * SNR is clearly silly. Silly indicates that the intensity + * was zero. */ + snr = fabs(intensity)/sigma; + + /* Record intensity and set redundancy to 1 on success */ + if ( !cutoff ) { + set_redundancy(refl, 1); + } else { + set_redundancy(refl, 0); + } + + if ( snr > 3.0 ) { + score++; + } else { + score--; + } + + //STATUS("%5.2f A, %5.2f, %i\n", 1e10/il[i].res, snr, score); + if ( score == 0 ) { + limit = il[i].res; + cutoff = 1; + } + + } + + crystal_set_resolution_limit(cr, limit); + + free(il); +} + + +static void integrate_refine(Crystal *cr, struct image *image, int use_closer, + double min_snr, + double ir_inn, double ir_mid, double ir_out, + int integrate_saturated, int **bgMasks) +{ + RefList *reflections; + + /* Create initial list of reflections with nominal parameters */ + reflections = find_intersections(image, cr); + measure_all_intensities(reflections, image); + + /* Find resolution limit of pattern using this list */ + estimate_resolution(reflections, cr, image); + + reflist_free(reflections); + + STATUS("Initial resolution estimate = %.2f nm^-1 or %.2f A\n", + crystal_get_resolution_limit(cr)/1e9, + 1e9 / crystal_get_resolution_limit(cr)); + + /* Estimate the mosaicity of the crystal using this resolution limit */ + estimate_mosaicity(cr, image); + + /* Create new list of reflections with refined mosaicity */ + reflections = find_intersections(image, cr); + measure_all_intensities(reflections, image); + + estimate_resolution(reflections, cr, image); +} + + +static void integrate_rings(Crystal *cr, struct image *image, int use_closer, + double min_snr, + double ir_inn, double ir_mid, double ir_out, + int integrate_saturated, int **bgMasks) +{ + RefList *list; + Reflection *refl; + RefListIterator *iter; + UnitCell *cell; + int n_saturated = 0; + double limit = 0.0; + + list = find_intersections(image, cr); + if ( list == NULL ) return; + + if ( num_reflections(list) == 0 ) return; + + cell = crystal_get_cell(cr); + + for ( refl = first_refl(list, &iter); + refl != NULL; + refl = next_refl(refl, iter) ) + { + double fs, ss, intensity; + double d; + int idx; + double sigma, snr; + double pfs, pss; + int r; + signed int h, k, l; + struct panel *p; + int pnum, j, found; + int saturated; + double one_over_d; + + get_detector_pos(refl, &pfs, &pss); + get_indices(refl, &h, &k, &l); + + /* 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; + } + + /* FIXME: Horrible hardcoded value */ + if ( (f != NULL) && (d < 10.0) ) { + + double exe; + + exe = get_excitation_error(refl); + + pfs = f->fs; + pss = f->ss; + + set_detector_pos(refl, exe, pfs, pss); + + } + + } + + p = find_panel(image->det, pfs, pss); + if ( p == NULL ) continue; /* Next peak */ + found = 0; + for ( j=0; j<image->det->n_panels; j++ ) { + if ( &image->det->panels[j] == p ) { + pnum = j; + found = 1; + break; + } + } + if ( !found ) { + ERROR("Couldn't find panel %p in list.\n", p); + return; + } + + r = integrate_peak(image, pfs, pss, &fs, &ss, + &intensity, &sigma, ir_inn, ir_mid, ir_out, + bgMasks[pnum], &saturated); + + if ( !r && saturated ) { + n_saturated++; + if ( !integrate_saturated ) r = 1; + } + + /* I/sigma(I) cutoff + * Rejects reflections below --min-integration-snr, or if the + * SNR is clearly silly. Silly indicates that the intensity + * was zero. */ + snr = fabs(intensity)/sigma; + if ( !r && (isnan(snr) || (snr < min_snr)) ) { + r = 1; + } + + /* Record intensity and set redundancy to 1 on success */ + if ( !r ) { + set_intensity(refl, intensity); + set_esd_intensity(refl, sigma); + set_redundancy(refl, 1); + } else { + set_redundancy(refl, 0); + } + + one_over_d = resolution(cell, h, k, l); + if ( one_over_d > limit ) limit = one_over_d; + + } + + crystal_set_num_saturated_reflections(cr, n_saturated); + crystal_set_resolution_limit(cr, limit); + crystal_set_reflections(cr, list); +} + + +void integrate_all(struct image *image, IntegrationMethod meth, + int use_closer, 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_inn); + if ( mask == NULL ) { + ERROR("Couldn't create background mask.\n"); + return; + } + bgMasks[i] = mask; + } + + for ( i=0; i<image->n_crystals; i++ ) { + + switch ( meth & INTEGRATION_METHOD_MASK ) { + + case INTEGRATION_NONE : + return; + + case INTEGRATION_RINGS : + integrate_rings(image->crystals[i], image, use_closer, + min_snr, ir_inn, ir_mid, ir_out, + integrate_saturated, bgMasks); + return; + + case INTEGRATION_REFINE : + integrate_refine(image->crystals[i], image, use_closer, + min_snr, ir_inn, ir_mid, ir_out, + integrate_saturated, bgMasks); + return; + + default : + ERROR("Unrecognised integration method %i\n", meth); + return; + + } + + } + + for ( i=0; i<image->det->n_panels; i++ ) { + free(bgMasks[i]); + } + free(bgMasks); +} + + +IntegrationMethod integration_method(const char *str, int *err) +{ + int n, i; + char **methods; + IntegrationMethod meth = INTEGRATION_NONE; + + if ( err != NULL ) *err = 0; + n = assplode(str, ",-", &methods, ASSPLODE_NONE); + + for ( i=0; i<n; i++ ) { + + if ( strcmp(methods[i], "rings") == 0) { + meth = INTEGRATION_DEFAULTS_RINGS; + + } else if ( strcmp(methods[i], "refine") == 0) { + meth = INTEGRATION_DEFAULTS_REFINE; + + } else if ( strcmp(methods[i], "none") == 0) { + return INTEGRATION_NONE; + + } else { + ERROR("Bad integration method: '%s'\n", str); + if ( err != NULL ) *err = 1; + return INTEGRATION_NONE; + } + + free(methods[i]); + + } + free(methods); + + return meth; + +} diff --git a/libcrystfel/src/integration.h b/libcrystfel/src/integration.h new file mode 100644 index 00000000..12e47ff5 --- /dev/null +++ b/libcrystfel/src/integration.h @@ -0,0 +1,75 @@ +/* + * integration.h + * + * Integration of intensities + * + * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. + * + * Authors: + * 2010-2013 Thomas White <taw@physics.org> + * + * This file is part of CrystFEL. + * + * CrystFEL is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * CrystFEL is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with CrystFEL. If not, see <http://www.gnu.org/licenses/>. + * + */ + +#ifndef INTEGRATION_H +#define INTEGRATION_H + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#define INTEGRATION_DEFAULTS_RINGS (INTEGRATION_RINGS | INTEGRATION_SATURATED) +#define INTEGRATION_DEFAULTS_REFINE (INTEGRATION_REFINE | INTEGRATION_SATURATED) + +/** + * IntegrationMethod: + * @INTEGRATION_NONE: No integration at all + * @INTEGRATION_RINGS: Predict reflections and integrate them all + * @INTEGRATION_REFINE: As @INTEGRATION_RINGS, but + * + * @INTEGRATION_SATURATED: Integrate saturated reflections + * + * An enumeration of all the available integration methods. + **/ +typedef enum { + + INTEGRATION_NONE = 0, + + /* The core integration methods themselves */ + INTEGRATION_RINGS = 1, + INTEGRATION_REFINE = 2, + + /* Bits at the top of the IntegrationMethod are flags which modify the + * behaviour of the integration. */ + INTEGRATION_SATURATED = 256, + INTEGRATION_CLOSER = 512, + +} IntegrationMethod; + +/* This defines the bits in "IntegrationMethod" which are used to represent the + * core of the integration method */ +#define INTEGRATION_METHOD_MASK (0xff) + +extern IntegrationMethod integration_method(const char *t, int *err); + +extern void integrate_all(struct image *image, IntegrationMethod meth, + int use_closer, double min_snr, + double ir_inn, double ir_mid, double ir_out, + int integrate_saturated); + +#endif /* INTEGRATION_H */ diff --git a/libcrystfel/src/peaks.c b/libcrystfel/src/peaks.c index 4b8a8ee8..05c564ba 100644 --- a/libcrystfel/src/peaks.c +++ b/libcrystfel/src/peaks.c @@ -52,6 +52,7 @@ #include "reflist-utils.h" #include "beam-parameters.h" #include "cell-utils.h" +#include "geometry.h" /* Degree of polarisation of X-ray beam */ @@ -200,7 +201,7 @@ static void add_crystal_to_mask(struct image *image, struct panel *p, /* cfs, css relative to panel origin */ -static int *make_BgMask(struct image *image, struct panel *p, double ir_inn) +int *make_BgMask(struct image *image, struct panel *p, double ir_inn) { int *mask; int w, h; @@ -224,11 +225,11 @@ static int *make_BgMask(struct image *image, struct panel *p, double ir_inn) /* Returns non-zero if peak has been vetoed. * i.e. don't use result if return value is not zero. */ -static int integrate_peak(struct image *image, int cfs, int css, - double *pfs, double *pss, - double *intensity, double *sigma, - double ir_inn, double ir_mid, double ir_out, - int *bgPkMask, int *saturated) +int integrate_peak(struct image *image, int cfs, int css, + double *pfs, double *pss, + double *intensity, double *sigma, + double ir_inn, double ir_mid, double ir_out, + int *bgPkMask, int *saturated) { signed int dfs, dss; double lim_sq, out_lim_sq, mid_lim_sq; @@ -653,233 +654,6 @@ int peak_sanity_check(struct image *image, Crystal **crystals, int n_cryst) } -struct integr_ind -{ - double res; - Reflection *refl; -}; - - -static int compare_resolution(const void *av, const void *bv) -{ - const struct integr_ind *a = av; - const struct integr_ind *b = bv; - - return a->res > b->res; -} - - -static struct integr_ind *sort_reflections(RefList *list, UnitCell *cell, - int *np) -{ - struct integr_ind *il; - Reflection *refl; - RefListIterator *iter; - int i, n; - - n = num_reflections(list); - *np = 0; /* For now */ - - if ( n == 0 ) return NULL; - - il = calloc(n, sizeof(struct integr_ind)); - if ( il == NULL ) return NULL; - - i = 0; - for ( refl = first_refl(list, &iter); - refl != NULL; - refl = next_refl(refl, iter) ) - { - signed int h, k, l; - double res; - - get_indices(refl, &h, &k, &l); - res = resolution(cell, h, k, l); - - il[i].res = res; - il[i].refl = refl; - - i++; - assert(i <= n); - } - - qsort(il, n, sizeof(struct integr_ind), compare_resolution); - - *np = n; - return il; -} - - -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, - int res_cutoff) -{ - RefList *reflections; - struct integr_ind *il; - int n, i; - double av = 0.0; - int first = 1; - int n_saturated = 0; - - reflections = crystal_get_reflections(cr); - - if ( num_reflections(reflections) == 0 ) return; - - il = sort_reflections(reflections, crystal_get_cell(cr), &n); - if ( il == NULL ) { - ERROR("Couldn't sort reflections\n"); - return; - } - - for ( i=0; i<n; i++ ) { - - double fs, ss, intensity; - double d; - int idx; - double sigma, snr; - double pfs, pss; - int r; - Reflection *refl; - signed int h, k, l; - struct panel *p; - int pnum, j, found; - int saturated; - - refl = il[i].refl; - - get_detector_pos(refl, &pfs, &pss); - get_indices(refl, &h, &k, &l); - - /* 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; - } - - /* FIXME: Horrible hardcoded value */ - if ( (f != NULL) && (d < 10.0) ) { - - double exe; - - exe = get_excitation_error(refl); - - pfs = f->fs; - pss = f->ss; - - set_detector_pos(refl, exe, pfs, pss); - - } - - } - - p = find_panel(image->det, pfs, pss); - if ( p == NULL ) continue; /* Next peak */ - found = 0; - for ( j=0; j<image->det->n_panels; j++ ) { - if ( &image->det->panels[j] == p ) { - pnum = j; - found = 1; - break; - } - } - if ( !found ) { - ERROR("Couldn't find panel %p in list.\n", p); - return; - } - - r = integrate_peak(image, pfs, pss, &fs, &ss, - &intensity, &sigma, ir_inn, ir_mid, ir_out, - bgMasks[pnum], &saturated); - - if ( !r && saturated ) { - n_saturated++; - if ( !integrate_saturated ) r = 1; - } - - /* I/sigma(I) cutoff - * Rejects reflections below --min-integration-snr, or if the - * SNR is clearly silly. Silly indicates that the intensity - * was zero. */ - snr = fabs(intensity)/sigma; - if ( isnan(snr) || (snr < min_snr) ) { - r = 1; - } - - /* Record intensity and set redundancy to 1 on success */ - if ( r == 0 ) { - set_intensity(refl, intensity); - set_esd_intensity(refl, sigma); - set_redundancy(refl, 1); - } else { - set_redundancy(refl, 0); - } - - if ( snr > 1.0 ) { - if ( first ) { - av = snr; - first = 0; - } else { - av = av + 0.1*(snr - av); - } - //STATUS("%5.2f A, %5.2f, av %5.2f\n", - // 1e10/il[i].res, snr, av); - if ( res_cutoff && (av < 1.0) ) break; - } - } - - 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 res_cutoff) -{ - 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_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, res_cutoff); - } - - for ( i=0; i<image->det->n_panels; i++ ) { - free(bgMasks[i]); - } - free(bgMasks); -} - - void validate_peaks(struct image *image, double min_snr, int ir_inn, int ir_mid, int ir_out, int use_saturated) { diff --git a/libcrystfel/src/peaks.h b/libcrystfel/src/peaks.h index f788bae5..b94ebee5 100644 --- a/libcrystfel/src/peaks.h +++ b/libcrystfel/src/peaks.h @@ -37,16 +37,13 @@ #include "reflist.h" +extern int *make_BgMask(struct image *image, struct panel *p, double ir_inn); + extern void search_peaks(struct image *image, float threshold, float min_gradient, float min_snr, double ir_inn, double ir_mid, double ir_out, int use_saturated); -extern 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 res_cutoff); - extern int peak_sanity_check(struct image *image, Crystal **crystals, int n_cryst); @@ -54,4 +51,10 @@ extern void validate_peaks(struct image *image, double min_snr, int ir_inn, int ir_mid, int ir_out, int use_saturated); +extern int integrate_peak(struct image *image, int cfs, int css, + double *pfs, double *pss, + double *intensity, double *sigma, + double ir_inn, double ir_mid, double ir_out, + int *bgPkMask, int *saturated); + #endif /* PEAKS_H */ diff --git a/src/indexamajig.c b/src/indexamajig.c index 208c0ed0..7d26c15d 100644 --- a/src/indexamajig.c +++ b/src/indexamajig.c @@ -66,6 +66,7 @@ #include "stream.h" #include "reflist-utils.h" #include "cell-utils.h" +#include "integration.h" #include "im-sandbox.h" @@ -100,6 +101,7 @@ static void show_help(const char *s) " hdf5 : Get from a table in HDF5 file.\n" " --hdf5-peaks=<p> Find peaks table in HDF5 file here.\n" " Default: /processing/hitfinder/peakinfo\n" +" --integration=<meth> Perform final pattern integration using <meth>.\n" "\n\n" "For more control over the process, you might need:\n\n" " --tolerance=<tol> Set the tolerances for cell comparison.\n" @@ -139,17 +141,12 @@ static void show_help(const char *s) " --no-check-prefix Don't attempt to correct the --prefix.\n" " --closer-peak Don't integrate from the location of a nearby peak\n" " instead of the predicted spot. Don't use.\n" -" --no-bg-sub Don't subtract local background estimates from\n" -" integrated intensities.\n" " --use-saturated During the initial peak search, don't reject\n" " peaks which contain pixels above max_adu.\n" " --integrate-saturated During the final integration stage, don't reject\n" " peaks which contain pixels above max_adu.\n" " --no-revalidate Don't re-integrate and check HDF5 peaks for\n" " validity.\n" -" --integrate-found Skip the spot prediction step, and just integrate\n" -" the intensities of the spots found by the initial\n" -" peak search.\n" " --no-peaks-in-stream Do not record peak search results in the stream.\n" " --no-refls-in-stream Do not record integrated reflections in the stream.\n" ); @@ -179,6 +176,7 @@ int main(int argc, char *argv[]) char *use_this_one_instead; struct index_args iargs; char *intrad = NULL; + char *int_str = NULL; /* Defaults */ iargs.cell = NULL; @@ -186,7 +184,6 @@ int main(int argc, char *argv[]) iargs.median_filter = 0; iargs.satcorr = 1; iargs.closer = 0; - iargs.bgsub = 1; iargs.tols[0] = 5.0; iargs.tols[1] = 5.0; iargs.tols[2] = 5.0; @@ -207,10 +204,8 @@ int main(int argc, char *argv[]) iargs.use_saturated = 0; iargs.integrate_saturated = 0; iargs.no_revalidate = 0; - iargs.integrate_found = 0; iargs.stream_peaks = 1; iargs.stream_refls = 1; - iargs.res_cutoff = 0; iargs.copyme = new_copy_hdf5_field_list(); if ( iargs.copyme == NULL ) { ERROR("Couldn't allocate HDF5 field list.\n"); @@ -218,6 +213,7 @@ int main(int argc, char *argv[]) } iargs.indm = NULL; /* No default */ iargs.ipriv = NULL; /* No default */ + iargs.int_meth = integration_method("rings", NULL); /* Long options */ const struct option longopts[] = { @@ -242,15 +238,11 @@ int main(int argc, char *argv[]) {"no-closer-peak", 0, &iargs.closer, 0}, {"closer-peak", 0, &iargs.closer, 1}, {"basename", 0, &config_basename, 1}, - {"bg-sub", 0, &iargs.bgsub, 1}, - {"no-bg-sub", 0, &iargs.bgsub, 0}, {"no-peaks-in-stream", 0, &iargs.stream_peaks, 0}, {"no-refls-in-stream", 0, &iargs.stream_refls, 0}, - {"res-cutoff", 0, &iargs.res_cutoff, 1}, {"integrate-saturated",0, &iargs.integrate_saturated,1}, {"use-saturated", 0, &iargs.use_saturated, 1}, {"no-revalidate", 0, &iargs.no_revalidate, 1}, - {"integrate-found", 0, &iargs.integrate_found, 1}, /* Long-only options with arguments */ {"peaks", 1, NULL, 2}, @@ -267,6 +259,7 @@ int main(int argc, char *argv[]) {"tolerance", 1, NULL, 13}, {"int-radius", 1, NULL, 14}, {"median-filter", 1, NULL, 15}, + {"integration", 1, NULL, 16}, {0, 0, NULL, 0} }; @@ -391,6 +384,10 @@ int main(int argc, char *argv[]) iargs.median_filter = atoi(optarg); break; + case 16 : + int_str = strdup(optarg); + break; + case 0 : break; @@ -465,6 +462,18 @@ int main(int argc, char *argv[]) free(indm_str); } + if ( int_str != NULL ) { + + int err; + + iargs.int_meth = integration_method(int_str, &err); + if ( err ) { + ERROR("Invalid integration method '%s'\n", int_str); + return 1; + } + free(int_str); + } + if ( toler != NULL ) { int ttt; ttt = sscanf(toler, "%f,%f,%f,%f", diff --git a/src/process_image.c b/src/process_image.c index ef673c84..bf694623 100644 --- a/src/process_image.c +++ b/src/process_image.c @@ -46,6 +46,7 @@ #include "stream.h" #include "reflist-utils.h" #include "process_image.h" +#include "integration.h" void process_image(const struct index_args *iargs, struct pattern_args *pargs, @@ -188,34 +189,21 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs, /* Integrate each crystal's diffraction spots */ for ( i=0; i<image.n_crystals; i++ ) { - RefList *reflections; - /* Set default crystal parameter(s) */ crystal_set_profile_radius(image.crystals[i], image.beam->profile_radius); - - if ( iargs->integrate_found ) { - reflections = select_intersections(&image, - image.crystals[i]); - } else { - reflections = find_intersections(&image, - image.crystals[i]); - } - - crystal_set_reflections(image.crystals[i], reflections); + crystal_set_mosaicity(image.crystals[i], 2e-3); /* radians */ + crystal_set_image(image.crystals[i], &image); } /* Integrate all the crystals at once - need all the crystals so that * overlaps can be detected. */ - integrate_reflections(&image, iargs->closer, - iargs->bgsub, - iargs->min_int_snr, - iargs->ir_inn, - iargs->ir_mid, - iargs->ir_out, - iargs->integrate_saturated, - iargs->res_cutoff); + integrate_all(&image, iargs->int_meth, + iargs->closer, + iargs->min_int_snr, + iargs->ir_inn, iargs->ir_mid, iargs->ir_out, + iargs->integrate_saturated); write_chunk(st, &image, hdfile, iargs->stream_peaks, iargs->stream_refls); diff --git a/src/process_image.h b/src/process_image.h index 5fe11e33..03041562 100644 --- a/src/process_image.h +++ b/src/process_image.h @@ -34,6 +34,9 @@ #endif +#include "integration.h" + + enum { PEAK_ZAEF, PEAK_HDF5, @@ -49,7 +52,6 @@ struct index_args int median_filter; int satcorr; int closer; - int bgsub; float threshold; float min_gradient; float min_snr; @@ -69,10 +71,9 @@ struct index_args int integrate_saturated; int use_saturated; int no_revalidate; - int integrate_found; int stream_peaks; int stream_refls; - int res_cutoff; + IntegrationMethod int_meth; }; |