Skip to content
GitLab
Commit 4be9287d authored by Thomas White's avatar Thomas White
Browse files

indexamajig: Take integration radii on command line

parent 385cd10c
Hide whitespace changes
Inline Side-by-side
doc/man/crystfel_geometry.5
View file @ 4be9287d
......@@ -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
......
doc/man/indexamajig.1
View file @ 4be9287d
......@@ -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
......
libcrystfel/src/peaks.c
View file @ 4be9287d
......@@ -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 ) {
......
libcrystfel/src/peaks.h
View file @ 4be9287d
......@@ -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);
......
src/indexamajig.c
View file @ 4be9287d
......@@ -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;
......
tests/integration_check.c
View file @ 4be9287d
......@@ -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);
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

imprint | privacy | Conditions of Use