Commit a1e0854a authored by Thomas White's avatar Thomas White
Browse files

indexamajig: Add --integrate-found

parent 10a8969c
......@@ -334,6 +334,11 @@ Normally, reflections which contain one or more pixels above max_adu (defined in
.PD
When using \fB--peaks=hdf5\fR, the peaks will be put through the same checks as if you were using \fB--peaks=zaef\fR. These checks reject peaks which are too close to panel edges, are saturated (unless you use \fB--use-saturated\fR), fall short of the minimum SNR value given by \fB--min-snr\fR, have other nearby peaks (closer than twice the inner integration radius, see \fB--int-radius\fR), or have any part in a bad region. Using this option skips this validation step, and uses the peaks directly.
.PD 0
.IP \fB--integrate-found\fR
.PD
Without this option, peaks will be predicted in each pattern based on the crystal orientation from autoindexing and the parameters in the beam description file. With this option, indices will be assigned to the peaks found by the peak search, and the positions of those peaks used for the final integration stage.
.SH BUGS
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
......
......@@ -33,6 +33,7 @@
#include <stdlib.h>
#include <assert.h>
#include <fenv.h>
#include "utils.h"
#include "cell.h"
......@@ -276,6 +277,65 @@ RefList *find_intersections(struct image *image, UnitCell *cell)
}
RefList *select_intersections(struct image *image, UnitCell *cell)
{
double ax, ay, az;
double bx, by, bz;
double cx, cy, cz;
const double min_dist = 0.25;
RefList *list;
int i;
/* Round towards nearest */
fesetround(1);
/* Cell basis vectors for this image */
cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
list = reflist_new();
if ( list == NULL ) return NULL;
/* 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;
double dsq;
f = image_get_feature(image->features, i);
if ( f == NULL ) continue;
/* 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 */
dsq = pow(h-hd, 2.0) + pow(k-kd, 2.0) + pow(l-ld, 2.0);
if ( sqrt(dsq) < min_dist ) {
Reflection *refl;
refl = add_refl(list, h, k, l);
set_detector_pos(refl, sqrt(dsq), f->fs, f->ss);
}
}
return list;
}
/* Calculate partialities and apply them to the image's reflections */
void update_partialities(struct image *image)
{
......
......@@ -43,9 +43,10 @@
extern "C" {
#endif
RefList *find_intersections(struct image *image, UnitCell *cell);
extern RefList *find_intersections(struct image *image, UnitCell *cell);
extern RefList *select_intersections(struct image *image, UnitCell *cell);
void update_partialities(struct image *image);
extern void update_partialities(struct image *image);
#ifdef __cplusplus
}
......
......@@ -301,9 +301,17 @@ static void process_image(const struct index_args *iargs,
iargs->config_insane, iargs->tols);
if ( image.indexed_cell != NULL ) {
pargs->indexable = 1;
image.reflections = find_intersections(&image,
image.indexed_cell);
if ( iargs->integrate_found ) {
image.reflections = select_intersections(&image,
image.indexed_cell);
} else {
image.reflections = find_intersections(&image,
image.indexed_cell);
}
if (image.reflections != NULL) {
integrate_reflections(&image,
iargs->config_closer,
......@@ -314,6 +322,7 @@ static void process_image(const struct index_args *iargs,
iargs->ir_out,
iargs->integrate_saturated);
}
} else {
image.reflections = NULL;
}
......
......@@ -69,6 +69,7 @@ struct index_args
int integrate_saturated;
int use_saturated;
int no_revalidate;
int integrate_found;
};
......
......@@ -167,6 +167,9 @@ static void show_help(const char *s)
" 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"
);
}
......@@ -251,6 +254,7 @@ int main(int argc, char *argv[])
int integrate_saturated = 0;
int use_saturated = 0;
int no_revalidate = 0;
int integrate_found = 0;
copyme = new_copy_hdf5_field_list();
if ( copyme == NULL ) {
......@@ -300,6 +304,7 @@ int main(int argc, char *argv[])
{"integrate-saturated",0, &integrate_saturated,1},
{"use-saturated",0, &use_saturated, 1},
{"no-revalidate", 0, &no_revalidate, 1},
{"integrate-found", 0, &integrate_found, 1},
{0, 0, NULL, 0}
};
......@@ -649,6 +654,7 @@ int main(int argc, char *argv[])
iargs.use_saturated = use_saturated;
iargs.integrate_saturated = integrate_saturated;
iargs.no_revalidate = no_revalidate;
iargs.integrate_found = integrate_found;
create_sandbox(&iargs, n_proc, prefix, config_basename, fh,
use_this_one_instead, ofh);
......
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