From 4be9287d7ae69bcfee52efed5a108f31c0f40f7e Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 9 Mar 2012 18:20:36 +0100 Subject: indexamajig: Take integration radii on command line --- doc/man/crystfel_geometry.5 | 6 --- doc/man/indexamajig.1 | 26 ++++++++++++- libcrystfel/src/peaks.c | 36 ++++++++--------- libcrystfel/src/peaks.h | 12 ++---- src/indexamajig.c | 94 ++++++++++++++++++++++++++++++--------------- tests/integration_check.c | 13 ++++--- 6 files changed, 118 insertions(+), 69 deletions(-) diff --git a/doc/man/crystfel_geometry.5 b/doc/man/crystfel_geometry.5 index de7054cb..abf52acc 100644 --- a/doc/man/crystfel_geometry.5 +++ b/doc/man/crystfel_geometry.5 @@ -96,12 +96,6 @@ panel0/res = 9090.91 .br panel0/peak_sep = 6.0 -; You need to specify the peak integration radius, which should be a little -.br -; larger than the actual radii of the peaks in pixels -.br -panel0/integr_radius = 2.0 - ; The camera length (in metres) for this panel .br ; You can also specify the HDF path to a scalar floating point value containing diff --git a/doc/man/indexamajig.1 b/doc/man/indexamajig.1 index 4160d25c..4360bd1d 100644 --- a/doc/man/indexamajig.1 +++ b/doc/man/indexamajig.1 @@ -131,6 +131,30 @@ If the unit cell is centered (i.e. if the space group begins with I, F, H, A, B, .PP The default is \fB--cell-reduction=none\fR. +.SH PEAK INTEGRATION +If the pattern could be successfully indexed and the cell reduction gave a +positive match, peaks will be predicted in the pattern and their intensities +measured. The peak integration relies on knowing an accurate radius to +integrate the peak within, and that the annulus used to estimate the background +level is not so big that that it covers neighbouring peaks. However, +indexamajig cannot (yet) determine these values for you. You need to specify +them using the \fB--int-radius\fR option (see below). +.PP +To determine appropriate values, index some patterns with the default values and +view the results using \fBcheck-near-bragg\fR (in the scripts folder). Set the +binning in \fBhdfsee\fR to 1, and adjust the ring radius until none of the rings +overlap for any of the patterns. This ring radius is the outer radius to use. +Then reduce the radius until the circles match the sizes of the peaks as +closely as possible. This value is the inner radius. The middle radius should +be between the two, ideally between two and three pixels smaller than the outer +radius. +.PP +If it's difficult to do this without setting the middle radius to the +same value as the inner radius, then the peaks are too close together to be +accurately integrated. Perhaps you got greedy with the resolution and put the +detector too close to the interaction region? Improved integration algorithms, +designed to handle such difficult cases, are under development. + .SH OPTIONS .PD 0 .IP "\fB-i\fR \fIfilename\fR" @@ -342,8 +366,6 @@ Don't run more than one indexamajig jobs simultaneously in the same working directory - they'll overwrite each other's DirAx or MOSFLM files, causing subtle problems which can't easily be detected. .PP -The way indexamajig calculates the intensities and errors of peaks sucks, for want of a better word. In particular, I/sigma(I) is not accurately estimated and many of the options don't make sense. Development of improved integration algorithms is in progress, so this issue will be addressed in a future version of CrystFEL. -.PP ReAx indexing is experimental. It works very nicely for some people, and crashes for others. In a future version, it will be improved and fully supported. .SH AUTHOR diff --git a/libcrystfel/src/peaks.c b/libcrystfel/src/peaks.c index a5b0e907..974d5b2e 100644 --- a/libcrystfel/src/peaks.c +++ b/libcrystfel/src/peaks.c @@ -155,10 +155,10 @@ static int cull_peaks(struct image *image) /* 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 *sigma) + double *pfs, double *pss, double *intensity, double *sigma, + double ir_inn, double ir_mid, double ir_out) { signed int fs, ss; - double lim, out_lim, mid_lim; double lim_sq, out_lim_sq, mid_lim_sq; double pk_total; int pk_counts; @@ -177,16 +177,13 @@ int integrate_peak(struct image *image, int cfs, int css, aduph = p->adu_per_eV * ph_lambda_to_eV(image->lambda); - 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); + lim_sq = pow(ir_inn, 2.0); + mid_lim_sq = pow(ir_mid, 2.0); + out_lim_sq = pow(ir_out, 2.0); /* Estimate the background */ - for ( fs=-out_lim; fs<+out_lim; fs++ ) { - for ( ss=-out_lim; ss<+out_lim; ss++ ) { + for ( fs=-ir_out; fs<+ir_out; fs++ ) { + for ( ss=-ir_out; ss<+ir_out; ss++ ) { double val; double tt = 0.0; @@ -236,8 +233,8 @@ int integrate_peak(struct image *image, int cfs, int css, pk_total = 0.0; pk_counts = 0; fsct = 0.0; ssct = 0.0; - for ( fs=-lim; fs<+lim; fs++ ) { - for ( ss=-lim; ss<+lim; ss++ ) { + for ( fs=-ir_inn; fs<+ir_inn; fs++ ) { + for ( ss=-ir_inn; ss<+ir_inn; ss++ ) { double val; double tt = 0.0; @@ -301,7 +298,8 @@ int integrate_peak(struct image *image, int cfs, int css, static void search_peaks_in_panel(struct image *image, float threshold, float min_gradient, float min_snr, - struct panel *p) + struct panel *p, + double ir_inn, double ir_mid, double ir_out) { int fs, ss, stride; float *data; @@ -407,7 +405,8 @@ static void search_peaks_in_panel(struct image *image, float threshold, /* Centroid peak and get better coordinates. */ r = integrate_peak(image, mask_fs, mask_ss, - &f_fs, &f_ss, &intensity, &sigma); + &f_fs, &f_ss, &intensity, &sigma, + ir_inn, ir_mid, ir_out); if ( r ) { /* Bad region - don't detect peak */ @@ -459,7 +458,7 @@ static void search_peaks_in_panel(struct image *image, float threshold, void search_peaks(struct image *image, float threshold, float min_gradient, - float min_snr) + float min_snr, double ir_inn, double ir_mid, double ir_out) { int i; @@ -474,7 +473,7 @@ void search_peaks(struct image *image, float threshold, float min_gradient, if ( p->no_index ) continue; search_peaks_in_panel(image, threshold, min_gradient, - min_snr, p); + min_snr, p, ir_inn, ir_mid, ir_out); } } @@ -611,7 +610,8 @@ 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 min_snr, + double ir_inn, double ir_mid, double ir_out) { struct integr_ind *il; int n, i; @@ -660,7 +660,7 @@ void integrate_reflections(struct image *image, int use_closer, int bgsub, } r = integrate_peak(image, pfs, pss, &fs, &ss, - &intensity, &sigma); + &intensity, &sigma, ir_inn, ir_mid, ir_out); /* Record intensity and set redundancy to 1 on success */ if ( r == 0 ) { diff --git a/libcrystfel/src/peaks.h b/libcrystfel/src/peaks.h index e5e31b91..8940a62e 100644 --- a/libcrystfel/src/peaks.h +++ b/libcrystfel/src/peaks.h @@ -34,22 +34,18 @@ #include "reflist.h" extern void search_peaks(struct image *image, float threshold, - float min_gradient, float min_snr); + float min_gradient, float min_snr, + double ir_inn, double ir_mid, double ir_out); extern void integrate_reflections(struct image *image, - int use_closer, int bgsub, double min_snr); + int use_closer, int bgsub, double min_snr, + double ir_inn, double ir_mid, double ir_out); extern double peak_lattice_agreement(struct image *image, UnitCell *cell, double *pst); extern int peak_sanity_check(struct image *image); -/* Exported so it can be poked by integration_check */ -extern int integrate_peak(struct image *image, - int cfs, int css, - double *pfs, double *pss, - double *intensity, double *sigma); - extern void estimate_resolution(RefList *list, UnitCell *cell, double *min, double *max); diff --git a/src/indexamajig.c b/src/indexamajig.c index ae610063..faee01fc 100644 --- a/src/indexamajig.c +++ b/src/indexamajig.c @@ -91,6 +91,9 @@ struct static_index_args struct beam_params *beam; const char *element; const char *hdf5_peak_path; + double ir_inn; + double ir_mid; + double ir_out; /* Output stream */ pthread_mutex_t *output_mutex; /* Protects the output stream */ @@ -178,34 +181,34 @@ static void show_help(const char *s) "The default is '--record=integrated'.\n" "\n\n" "For more control over the process, you might need:\n\n" -" --cell-reduction= Use as the cell reduction method. Choose from:\n" -" none : no matching, just use the raw cell.\n" -" reduce : full cell reduction.\n" -" compare : match by at most changing the order of\n" -" the indices.\n" -" compare_ab : compare 'a' and 'b' lengths only.\n" -" --tolerance= Set the tolerance for a,b,c axis (in %%)\n" -" and for the angles (in deg) when reducing\n" -" or comparing (default is 5%% and 1.5deg)\n" -" --filter-cm Perform common-mode noise subtraction on images\n" -" before proceeding. Intensities will be extracted\n" -" from the image as it is after this processing.\n" -" --filter-noise Apply an aggressive noise filter which sets all\n" -" pixels in each 3x3 region to zero if any of them\n" -" have negative values. Intensity measurement will\n" -" be performed on the image as it was before this.\n" -" --no-sat-corr Don't correct values of saturated peaks using a\n" -" table included in the HDF5 file.\n" -" --threshold= Only accept peaks above ADU. Default: 800.\n" -" --min-gradient= Minimum gradient for Zaefferer peak search.\n" -" Default: 100,000.\n" -" --min-snr= Minimum signal-to-noise ratio for peaks.\n" -" Default: 5.\n" -" --min-integration-snr= Minimum signal-to-noise ratio for peaks\n" -" during integration. Default: -infinity.\n" -" -e, --image= Use this image from the HDF5 file.\n" -" Example: /data/data0.\n" -" Default: The first one found.\n" +" --cell-reduction= Use as the cell reduction method. Choose from:\n" +" none : no matching, just use the raw cell.\n" +" reduce : full cell reduction.\n" +" compare : match by at most changing the order of\n" +" the indices.\n" +" compare_ab : compare 'a' and 'b' lengths only.\n" +" --tolerance= Set the tolerances for cell reduction.\n" +" Default: 5,5,5,1.5.\n" +" --filter-cm Perform common-mode noise subtraction on images\n" +" before proceeding. Intensities will be extracted\n" +" from the image as it is after this processing.\n" +" --filter-noise Apply an aggressive noise filter which sets all\n" +" pixels in each 3x3 region to zero if any of them\n" +" have negative values. Intensity measurement will\n" +" be performed on the image as it was before this.\n" +" --no-sat-corr Don't correct values of saturated peaks using a\n" +" table included in the HDF5 file.\n" +" --threshold= Only accept peaks above ADU. Default: 800.\n" +" --min-gradient= Minimum gradient for Zaefferer peak search.\n" +" Default: 100,000.\n" +" --min-snr= Minimum signal-to-noise ratio for peaks.\n" +" Default: 5.\n" +" --min-integration-snr= Minimum signal-to-noise ratio for peaks\n" +" during integration. Default: -infinity.\n" +" --int-radius= Set the integration radii. Default: 4,5,7.\n" +"-e, --image= Use this image from the HDF5 file.\n" +" Example: /data/data0.\n" +" Default: The first one found.\n" "\n" "\nFor time-resolved stuff, you might want to use:\n\n" " --copy-hdf5-field Copy the value of field into the stream. You\n" @@ -353,7 +356,10 @@ static void process_image(void *pp, int cookie) case PEAK_ZAEF : search_peaks(&image, pargs->static_args.threshold, pargs->static_args.min_gradient, - pargs->static_args.min_snr); + pargs->static_args.min_snr, + pargs->static_args.ir_inn, + pargs->static_args.ir_mid, + pargs->static_args.ir_out); break; } @@ -382,7 +388,10 @@ static void process_image(void *pp, int cookie) integrate_reflections(&image, pargs->static_args.config_closer, pargs->static_args.config_bgsub, - pargs->static_args.min_int_snr); + pargs->static_args.min_int_snr, + pargs->static_args.ir_inn, + pargs->static_args.ir_mid, + pargs->static_args.ir_out); } @@ -587,6 +596,10 @@ int main(int argc, char *argv[]) char *endptr; char *hdf5_peak_path = NULL; struct copy_hdf5_field *copyme; + char *intrad = NULL; + float ir_inn = 4.0; + float ir_mid = 5.0; + float ir_out = 7.0; copyme = new_copy_hdf5_field_list(); if ( copyme == NULL ) { @@ -631,6 +644,7 @@ int main(int argc, char *argv[]) {"min-snr", 1, NULL, 11}, {"min-integration-snr",1, NULL, 12}, {"tolerance", 1, NULL, 13}, + {"int-radius", 1, NULL, 14}, {0, 0, NULL, 0} }; @@ -762,6 +776,10 @@ int main(int argc, char *argv[]) toler = strdup(optarg); break; + case 14 : + intrad = strdup(optarg); + break; + case 0 : break; @@ -900,6 +918,19 @@ int main(int argc, char *argv[]) } } + if ( intrad != NULL ) { + int r; + r = sscanf(intrad, "%f,%f,%f", &ir_inn, &ir_mid, &ir_out); + if ( r != 3 ) { + ERROR("Invalid parameters for '--int-radius'\n"); + return 1; + } + } else { + STATUS("WARNING: You did not specify --int-radius.\n"); + STATUS("WARNING: I will use the default values, which are" + " probably not appropriate for your patterns.\n"); + } + if ( geometry == NULL ) { ERROR("You need to specify a geometry file with --geometry\n"); return 1; @@ -994,6 +1025,9 @@ int main(int argc, char *argv[]) qargs.static_args.stream_flags = stream_flags; qargs.static_args.hdf5_peak_path = hdf5_peak_path; qargs.static_args.copyme = copyme; + qargs.static_args.ir_inn = ir_inn; + qargs.static_args.ir_mid = ir_mid; + qargs.static_args.ir_out = ir_out; qargs.fh = fh; qargs.prefix = prefix; diff --git a/tests/integration_check.c b/tests/integration_check.c index 242e4076..eeccc6ec 100644 --- a/tests/integration_check.c +++ b/tests/integration_check.c @@ -32,10 +32,11 @@ #include #include -#include #include #include +#include "../libcrystfel/src/peaks.c" + /* The third integration check draws a Poisson background and checks that, on * average, it gets subtracted by the background subtraction. */ @@ -63,7 +64,7 @@ static void third_integration_check(struct image *image, int n_trials, } r = integrate_peak(image, 64, 64, &fsp, &ssp, - &intensity, &sigma); + &intensity, &sigma, 10.0, 15.0, 17.0); if ( r == 0 ) { mean_intensity += intensity; @@ -128,7 +129,7 @@ static void fourth_integration_check(struct image *image, int n_trials, } r = integrate_peak(image, 64, 64, &fsp, &ssp, - &intensity, &sigma); + &intensity, &sigma, 10.0, 15.0, 17.0); if ( r == 0 ) { mean_intensity += intensity; @@ -333,7 +334,8 @@ int main(int argc, char *argv[]) memset(image.data, 0, 128*128*sizeof(float)); /* First check: no intensity -> no peak, or very low intensity */ - r = integrate_peak(&image, 64, 64, &fsp, &ssp, &intensity, &sigma); + r = integrate_peak(&image, 64, 64, &fsp, &ssp, &intensity, &sigma, + 10.0, 15.0, 17.0); STATUS(" First check: integrate_peak() returned %i", r); if ( r == 0 ) { @@ -359,7 +361,8 @@ int main(int argc, char *argv[]) } } - r = integrate_peak(&image, 64, 64, &fsp, &ssp, &intensity, &sigma); + r = integrate_peak(&image, 64, 64, &fsp, &ssp, &intensity, &sigma, + 10.0, 15.0, 17.0); if ( r ) { ERROR(" Second check: integrate_peak() returned %i (wrong).\n", r); -- cgit v1.2.3