Commit 2f5b4463 authored by Thomas White's avatar Thomas White
Browse files

Don't constantly re-integrate peaks

parent f20d9a29
......@@ -37,10 +37,6 @@
#include "reflist.h"
/* Maximum number of iterations of NLSq to do for each image per macrocycle. */
#define MAX_CYCLES (100)
static void show_help(const char *s)
{
printf("Syntax: %s [options]\n\n", s);
......@@ -57,9 +53,6 @@ static void show_help(const char *s)
" initial values for parameters, and nominal\n"
" wavelengths if no per-shot value is found in \n"
" an HDF5 file.\n"
" -x, --prefix=<p> Prefix filenames from input file with <p>.\n"
" --basename Remove the directory parts of the filenames.\n"
" --no-check-prefix Don't attempt to correct the --prefix.\n"
" -y, --symmetry=<sym> Merge according to symmetry <sym>.\n"
" -n, --iterations=<n> Run <n> cycles of scaling and post-refinement.\n"
"\n"
......@@ -83,58 +76,8 @@ static void refine_image(int mytask, void *tasks)
struct refine_args *all_args = tasks;
struct refine_args *pargs = &all_args[mytask];
struct image *image = pargs->image;
double nominal_photon_energy = pargs->image->beam->photon_energy;
struct hdfile *hdfile;
int i;
double dev, last_dev;
RefList *reflections;
hdfile = hdfile_open(image->filename);
if ( hdfile == NULL ) {
ERROR("Couldn't open '%s'\n", image->filename);
return;
} else if ( hdfile_set_image(hdfile, "/data/data0") ) {
ERROR("Couldn't select path\n");
hdfile_close(hdfile);
return;
}
if ( hdf5_read(hdfile, pargs->image, 0, nominal_photon_energy) ) {
ERROR("Couldn't read '%s'\n", image->filename);
hdfile_close(hdfile);
return;
}
double a, b, c, al, be, ga;
cell_get_parameters(image->indexed_cell, &a, &b, &c, &al, &be, &ga);
STATUS("Initial cell: %5.2f %5.2f %5.2f nm %5.2f %5.2f %5.2f deg\n",
a/1.0e-9, b/1.0e-9, c/1.0e-9,
rad2deg(al), rad2deg(be), rad2deg(ga));
/* FIXME: Don't do this each time */
reflections = find_intersections(image, image->indexed_cell, 0);
dev = +INFINITY;
i = 0;
do {
last_dev = dev;
dev = pr_iterate(image, pargs->i_full, pargs->sym, reflections);
STATUS("Iteration %2i: mean dev = %5.2f\n", i, dev);
i++;
} while ( (fabs(last_dev - dev) > 1.0) && (i < MAX_CYCLES) );
mean_partial_dev(image, reflections, pargs->sym,
pargs->i_full, pargs->graph);
if ( pargs->pgraph ) {
fprintf(pargs->pgraph, "%5i %5.2f\n", mytask, dev);
}
free(image->data);
if ( image->flags != NULL ) free(image->flags);
hdfile_close(hdfile);
reflist_free(reflections);
/* Muppet proofing */
image->data = NULL;
image->flags = NULL;
pr_refine(image, pargs->i_full, pargs->sym);
}
......@@ -165,87 +108,6 @@ static void refine_all(struct image *images, int n_total_patterns,
}
static void uniquify(Reflection *refl, const char *sym)
{
signed int h, k, l;
signed int ha, ka, la;
get_indices(refl, &h, &k, &l);
get_asymm(h, k, l, &ha, &ka, &la, sym);
set_indices(refl, h, k, l);
}
/* FIXME: Get rid of this */
static void integrate_image(struct image *image, ReflItemList *obs,
const char *sym)
{
RefList *reflections;
Reflection *refl;
struct hdfile *hdfile;
double nominal_photon_energy = image->beam->photon_energy;
hdfile = hdfile_open(image->filename);
if ( hdfile == NULL ) {
ERROR("Couldn't open '%s'\n", image->filename);
return;
} else if ( hdfile_set_image(hdfile, "/data/data0") ) {
ERROR("Couldn't select path\n");
hdfile_close(hdfile);
return;
}
if ( hdf5_read(hdfile, image, 0, nominal_photon_energy) ) {
ERROR("Couldn't read '%s'\n", image->filename);
hdfile_close(hdfile);
return;
}
/* Figure out which spots should appear in this pattern */
reflections = find_intersections(image, image->indexed_cell, 0);
/* For each reflection, estimate the partiality */
for ( refl = first_refl(reflections);
refl != NULL;
refl = next_refl(refl) ) {
signed int h, k, l;
float i_partial;
float xc, yc;
double x, y;
uniquify(refl, sym);
get_indices(refl, &h, &k, &l);
/* Don't attempt to use spots with very small
* partialities, since it won't be accurate. */
if ( get_partiality(refl) < 0.1 ) continue;
/* Actual measurement of this reflection from this pattern? */
get_detector_pos(refl, &x, &y);
if ( integrate_peak(image, x, y,
&xc, &yc, &i_partial, NULL, NULL, 1, 0) ) {
delete_refl(refl);
continue;
}
set_int(refl, i_partial);
if ( !find_item(obs, h, k, l) ) add_item(obs, h, k, l);
}
image->reflections = reflections;
free(image->data);
if ( image->flags != NULL ) free(image->flags);
hdfile_close(hdfile);
/* Muppet proofing */
image->data = NULL;
image->flags = NULL;
}
/* Decide which reflections can be scaled */
static void select_scalable_reflections(struct image *images, int n)
{
......@@ -280,12 +142,9 @@ int main(int argc, char *argv[])
char *infile = NULL;
char *outfile = NULL;
char *geomfile = NULL;
char *prefix = NULL;
char *sym = NULL;
FILE *fh;
int nthreads = 1;
int config_basename = 0;
int config_checkprefix = 1;
struct detector *det;
unsigned int *cts;
ReflItemList *obs;
......@@ -303,9 +162,6 @@ int main(int argc, char *argv[])
{"output", 1, NULL, 'o'},
{"geometry", 1, NULL, 'g'},
{"beam", 1, NULL, 'b'},
{"prefix", 1, NULL, 'x'},
{"basename", 0, &config_basename, 1},
{"no-check-prefix", 0, &config_checkprefix, 0},
{"symmetry", 1, NULL, 'y'},
{"iterations", 1, NULL, 'n'},
{0, 0, NULL, 0}
......@@ -329,10 +185,6 @@ int main(int argc, char *argv[])
geomfile = strdup(optarg);
break;
case 'x' :
prefix = strdup(optarg);
break;
case 'j' :
nthreads = atoi(optarg);
break;
......@@ -387,15 +239,6 @@ int main(int argc, char *argv[])
outfile = strdup("facetron.hkl");
}
/* Sanitise prefix */
if ( prefix == NULL ) {
prefix = strdup("");
} else {
if ( config_checkprefix ) {
prefix = check_prefix(prefix);
}
}
if ( sym == NULL ) sym = strdup("1");
/* Get detector geometry */
......@@ -427,7 +270,11 @@ int main(int argc, char *argv[])
UnitCell *cell;
char *filename;
char *fnamereal;
char *rval;
char line[1024];
RefList *peaks;
RefList *transfer;
Reflection *refl;
if ( find_chunk(fh, &cell, &filename) == 1 ) {
ERROR("Couldn't get all of the filenames and cells"
......@@ -437,38 +284,83 @@ int main(int argc, char *argv[])
images[i].indexed_cell = cell;
/* Mangle the filename now */
if ( config_basename ) {
char *tmp;
tmp = safe_basename(filename);
free(filename);
filename = tmp;
}
fnamereal = malloc(1024);
snprintf(fnamereal, 1023, "%s%s", prefix, filename);
free(filename);
images[i].filename = fnamereal;
images[i].filename = filename;
images[i].div = beam->divergence;
images[i].bw = beam->bandwidth;
images[i].det = det;
images[i].beam = beam;
images[i].osf = 1.0;
images[i].profile_radius = 0.005e9;
images[i].reflections = reflist_new();
/* Muppet proofing */
images[i].data = NULL;
images[i].flags = NULL;
/* Get reflections from this image.
* FIXME: Use the ones from the stream */
integrate_image(&images[i], obs, sym);
/* Read integrated intensities from pattern */
peaks = reflist_new();
do {
Reflection *refl;
signed int h, k, l;
float intensity;
int r;
rval = fgets(line, 1023, fh);
chomp(line);
if ( (strlen(line) == 0) || (rval == NULL) ) break;
r = sscanf(line, "%i %i %i %f", &h, &k, &l, &intensity);
if ( r != 4 ) continue;
/* Not interested in the central beam */
if ( (h==0) && (k==0) && (l==0) ) continue;
refl = add_refl(peaks, h, k, l);
set_int(refl, intensity);
} while ( (strlen(line) != 0) && (rval != NULL) );
/* Calculate initial partialities and fill in intensities from
* the stream */
transfer = find_intersections(&images[i], cell, 0);
images[i].reflections = reflist_new();
for ( refl = first_refl(transfer);
refl != NULL;
refl = next_refl(refl) ) {
Reflection *peak;
Reflection *new;
signed int h, k, l, ha, ka, la;
double r1, r2, p, x, y;
int clamp1, clamp2;
get_indices(refl, &h, &k, &l);
peak = find_refl(peaks, h, k, l);
if ( peak == NULL ) {
if ( (h==0) && (k==0) && (l==0) ) continue;
STATUS("%3i %3i %3i not found\n", h, k, l);
continue;
}
get_asymm(h, k, l, &ha, &ka, &la, sym);
add_item(obs, ha, ka, la);
new = add_refl(images[i].reflections, ha, ka, la);
get_partial(refl, &r1, &r2, &p, &clamp1, &clamp2);
get_detector_pos(refl, &x, &y);
set_partial(new, r1, r2, p, clamp1, clamp2);
set_detector_pos(new, 0.0, x, y);
}
reflist_free(peaks);
reflist_free(transfer);
progress_bar(i, n_total_patterns-1, "Loading pattern data");
}
fclose(fh);
free(prefix);
cts = new_list_count();
......
......@@ -28,6 +28,28 @@
#include "cell.h"
/* Maximum number of iterations of NLSq to do for each image per macrocycle. */
#define MAX_CYCLES (100)
/* Refineable parameters */
enum {
REF_SCALE,
REF_DIV,
REF_R,
REF_ASX,
REF_BSX,
REF_CSX,
REF_ASY,
REF_BSY,
REF_CSY,
REF_ASZ,
REF_BSZ,
REF_CSZ,
NUM_PARAMS
};
/* Returns dp/dr at "r" */
static double partiality_gradient(double r, double profile_radius)
{
......@@ -65,7 +87,7 @@ static double partiality_rgradient(double r, double profile_radius)
/* Return the gradient of parameter 'k' given the current status of 'image'. */
double gradient(struct image *image, int k, Reflection *refl, double r)
static double gradient(struct image *image, int k, Reflection *refl, double r)
{
double ds, tt, azi;
double nom, den;
......@@ -187,7 +209,7 @@ static void apply_cell_shift(UnitCell *cell, int k, double shift)
/* Apply the given shift to the 'k'th parameter of 'image'. */
void apply_shift(struct image *image, int k, double shift)
static void apply_shift(struct image *image, int k, double shift)
{
switch ( k ) {
......@@ -226,67 +248,19 @@ void apply_shift(struct image *image, int k, double shift)
}
double mean_partial_dev(struct image *image, RefList *reflections,
const char *sym, double *i_full, FILE *graph)
{
int n_used;
double delta_I = 0.0;
Reflection *refl;
n_used = 0;
for ( refl = first_refl(reflections);
refl != NULL;
refl = next_refl(refl) ) {
signed int hind, kind, lind;
signed int ha, ka, la;
double I_full;
float I_partial;
float xc, yc;
double x, y;
double p;
get_indices(refl, &hind, &kind, &lind);
if ( !get_scalable(refl) ) continue;
/* Actual measurement of this reflection from this pattern? */
get_detector_pos(refl, &x, &y);
/* FIXME: Get rid of this */
if ( integrate_peak(image, x, y,
&xc, &yc, &I_partial, NULL, NULL,
1, 0) ) {
continue;
}
get_asymm(hind, kind, lind, &ha, &ka, &la, sym);
I_full = lookup_intensity(i_full, ha, ka, la);
p = get_partiality(refl);
delta_I += fabs(I_partial - (p * I_full / image->osf));
n_used++;
if ( graph != NULL ) {
fprintf(graph, "%3i %3i %3i %5.2f (at %5.2f,%5.2f)"
" %5.2f %5.2f\n",
hind, kind, lind, I_partial/p, xc, yc,
p, I_partial / I_full);
}
}
return delta_I / (double)n_used;
}
/* Perform one cycle of post refinement on 'image' against 'i_full' */
double pr_iterate(struct image *image, double *i_full, const char *sym,
RefList *reflections)
static double pr_iterate(struct image *image, const double *i_full,
const char *sym)
{
gsl_matrix *M;
gsl_vector *v;
gsl_vector *shifts;
int param;
Reflection *refl;
RefList *reflections;
double max_shift;
reflections = image->reflections;
M = gsl_matrix_calloc(NUM_PARAMS, NUM_PARAMS);
v = gsl_vector_calloc(NUM_PARAMS);
......@@ -355,14 +329,31 @@ double pr_iterate(struct image *image, double *i_full, const char *sym,
shifts = gsl_vector_alloc(NUM_PARAMS);
gsl_linalg_HH_solve(M, v, shifts);
max_shift = 0.0;
for ( param=0; param<NUM_PARAMS; param++ ) {
double shift = gsl_vector_get(shifts, param);
apply_shift(image, param, shift);
if ( fabs(shift) > max_shift ) max_shift = fabs(shift);
}
gsl_matrix_free(M);
gsl_vector_free(v);
gsl_vector_free(shifts);
return mean_partial_dev(image, reflections, sym, i_full, NULL);
return max_shift;
}
void pr_refine(struct image *image, const double *i_full, const char *sym)
{
double max_shift;
int i;
i = 0;
do {
max_shift = pr_iterate(image, i_full, sym);
STATUS("Iteration %2i: max shift = %5.2f\n", i, max_shift);
i++;
} while ( (max_shift > 0.01) && (i < MAX_CYCLES) );
}
......@@ -3,7 +3,7 @@
*
* Post refinement
*
* (c) 2006-2010 Thomas White <taw@physics.org>
* (c) 2006-2011 Thomas White <taw@physics.org>
*
* Part of CrystFEL - crystallography with a FEL
*
......@@ -21,36 +21,11 @@
#include <stdio.h>
#include "image.h"
#include "reflist.h"
#include "utils.h"
/* Refineable parameters */
enum {
REF_SCALE,
REF_DIV,
REF_R,
REF_ASX,
REF_BSX,
REF_CSX,
REF_ASY,
REF_BSY,
REF_CSY,
REF_ASZ,
REF_BSZ,
REF_CSZ,
NUM_PARAMS
};
/* Apply the given shift to the 'k'th parameter of 'image'. */
void apply_shift(struct image *image, int k, double shift);
double mean_partial_dev(struct image *image, RefList *reflections,
const char *sym, double *i_full, FILE *graph);
double pr_iterate(struct image *image, double *i_full, const char *sym,
RefList *reflections);
extern void pr_refine(struct image *image, const double *i_full,
const char *sym);
#endif /* POST_REFINEMENT_H */
......@@ -19,7 +19,6 @@
struct _reflection {
/* Listy stuff */
RefList *list;
unsigned int serial; /* Unique serial number, key */
struct _reflection *child[2]; /* Child nodes */
struct _reflection *parent; /* Parent node */
......@@ -224,32 +223,6 @@ void set_partial(Reflection *refl, double r1, double r2, double p,
}
void set_indices(Reflection *refl, signed int h, signed int k, signed int l)
{
/* Tricky, because the indices determine the position in the tree */
Reflection copy;
Reflection *new;
Reflection transfer;
/* Copy all data */
copy = *refl;
/* Delete and re-add with new indices */
delete_refl(refl);
new = add_refl(copy.list, h, k, l);
/* Transfer data back */
transfer = *new;
*new = copy;
new->list = transfer.list;
new->parent = transfer.parent;
new->child[0] = transfer.child[0];
new->child[1] = transfer.child[1];
new->h = transfer.h; new->k = transfer.k; new->l = transfer.l;
new->serial = transfer.serial;
}
void set_int(Reflection *refl, double intensity)
{
refl->intensity = intensity;
......
......@@ -46,8 +46,6 @@ extern void set_detector_pos(Reflection *refl, double exerr,
double x, double y);
extern void set_partial(Reflection *refl, double r1, double r2, double p,
double clamp_low, double clamp_high);
extern void set_indices(Reflection *refl,
signed int h, signed int k, signed int l);
extern void set_int(Reflection *refl, double intensity);
extern void set_scalable(Reflection *refl, int scalable);
......
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