aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--doc/man/crystfel_geometry.56
-rw-r--r--doc/man/indexamajig.126
-rw-r--r--libcrystfel/src/peaks.c36
-rw-r--r--libcrystfel/src/peaks.h12
-rw-r--r--src/indexamajig.c94
-rw-r--r--tests/integration_check.c13
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=<m> Use <m> 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=<a,b,c,angl> 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=<n> Only accept peaks above <n> ADU. Default: 800.\n"
-" --min-gradient=<n> Minimum gradient for Zaefferer peak search.\n"
-" Default: 100,000.\n"
-" --min-snr=<n> Minimum signal-to-noise ratio for peaks.\n"
-" Default: 5.\n"
-" --min-integration-snr=<n> Minimum signal-to-noise ratio for peaks\n"
-" during integration. Default: -infinity.\n"
-" -e, --image=<element> Use this image from the HDF5 file.\n"
-" Example: /data/data0.\n"
-" Default: The first one found.\n"
+" --cell-reduction=<m> Use <m> 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=<tol> 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=<n> Only accept peaks above <n> ADU. Default: 800.\n"
+" --min-gradient=<n> Minimum gradient for Zaefferer peak search.\n"
+" Default: 100,000.\n"
+" --min-snr=<n> Minimum signal-to-noise ratio for peaks.\n"
+" Default: 5.\n"
+" --min-integration-snr=<n> Minimum signal-to-noise ratio for peaks\n"
+" during integration. Default: -infinity.\n"
+" --int-radius=<r> Set the integration radii. Default: 4,5,7.\n"
+"-e, --image=<element> 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 <f> Copy the value of field <f> 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 <stdio.h>
#include <image.h>
-#include <peaks.h>
#include <utils.h>
#include <beam-parameters.h>
+#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);