diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/facetron.c | 201 | ||||
| -rw-r--r-- | src/image.h | 2 |
2 files changed, 200 insertions, 3 deletions
diff --git a/src/facetron.c b/src/facetron.c index f4148ae1..ca58ea2a 100644 --- a/src/facetron.c +++ b/src/facetron.c @@ -22,6 +22,9 @@ #include <getopt.h> #include <assert.h> #include <pthread.h> +#include <gsl/gsl_matrix.h> +#include <gsl/gsl_vector.h> +#include <gsl/gsl_linalg.h> #include "utils.h" #include "hdf5-file.h" @@ -34,6 +37,13 @@ #include "beam-parameters.h" +/* Refineable parameters */ +enum { + REF_SCALE, + NUM_PARAMS +}; + + static void show_help(const char *s) { printf("Syntax: %s [options]\n\n", s); @@ -68,11 +78,197 @@ struct refine_args }; +/* Return the gradient of parameter 'k' given the current status of 'image'. */ +static double gradient(struct image *image, int k, + struct cpeak spot, double I_partial) +{ + switch ( k ) { + + case REF_SCALE : + return I_partial; + + } + + ERROR("No gradient defined for parameter %i\n", k); + abort(); +} + + +/* Apply the given shift to the 'k'th parameter of 'image'. */ +static void apply_shift(struct image *image, int k, double shift) +{ + switch ( k ) { + + case REF_SCALE : + image->osf += shift; + break; + + default : + ERROR("No gradient defined for parameter %i\n", k); + abort(); + + } +} + + +static double mean_partial_dev(struct image *image, struct cpeak *spots, int n, + const char *sym, double *i_full) +{ + int h; + double delta_I; + + for ( h=0; h<n; h++ ) { + + signed int hind, kind, lind; + signed int ha, ka, la; + double I_full; + float I_partial; + float xc, yc; + + hind = spots[h].h; + kind = spots[h].k; + lind = spots[h].l; + + /* Don't attempt to use spots with very small + * partialities, since it won't be accurate. */ + if ( spots[h].p < 0.1 ) continue; + + /* Actual measurement of this reflection from this + * pattern? */ + /* FIXME: Coordinates aren't whole numbers */ + if ( integrate_peak(image, spots[h].x, spots[h].y, + &xc, &yc, &I_partial, NULL, NULL, 1, 1) ) { + continue; + } + + get_asymm(hind, kind, lind, &ha, &ka, &la, sym); + I_full = lookup_intensity(i_full, ha, ka, la); + delta_I += I_partial - spots[h].p * I_full; + + } + + return delta_I / (double)n; +} + + +static double iterate(struct image *image, double *i_full, const char *sym, + struct cpeak *spots, int n) +{ + gsl_matrix *M; + gsl_vector *v; + gsl_vector *shifts; + int h, shift; + + M = gsl_matrix_alloc(NUM_PARAMS, NUM_PARAMS); + v = gsl_vector_alloc(NUM_PARAMS); + + /* Construct the equations, one per reflection in this image */ + for ( h=0; h<n; h++ ) { + + signed int hind, kind, lind; + signed int ha, ka, la; + double I_full, delta_I; + float I_partial; + float xc, yc; + int k; + + hind = spots[h].h; + kind = spots[h].k; + lind = spots[h].l; + + /* Don't attempt to use spots with very small + * partialities, since it won't be accurate. */ + if ( spots[h].p < 0.1 ) continue; + + /* Actual measurement of this reflection from this + * pattern? */ + /* FIXME: Coordinates aren't whole numbers */ + if ( integrate_peak(image, spots[h].x, spots[h].y, + &xc, &yc, &I_partial, NULL, NULL, 1, 1) ) { + continue; + } + + get_asymm(hind, kind, lind, &ha, &ka, &la, sym); + I_full = lookup_intensity(i_full, ha, ka, la); + delta_I = I_partial - spots[h].p * I_full; + + for ( k=0; k<NUM_PARAMS; k++ ) { + + int g; + double v_c; + + for ( g=0; g<NUM_PARAMS; g++ ) { + + double M_curr, M_c; + + M_curr = gsl_matrix_get(M, g, k); + + M_c = gradient(image, g, spots[h], I_partial) + * gradient(image, k, spots[h], I_partial); + M_c *= pow(I_full, 2.0); + + gsl_matrix_set(M, g, k, M_curr + M_c); + + } + + v_c = delta_I * I_full * gradient(image, k, spots[h], + I_partial); + gsl_vector_set(v, k, v_c); + + } + + } + + shifts = gsl_vector_alloc(NUM_PARAMS); + gsl_linalg_HH_solve(M, v, shifts); + for ( shift=0; shift<NUM_PARAMS; shift++ ) { + apply_shift(image, shift, gsl_vector_get(shifts, shift)); + } + + free(spots); + spots = find_intersections(image, image->indexed_cell, &n, 0); + return mean_partial_dev(image, spots, n, sym, i_full); +} + + static void refine_image(int mytask, void *tasks) { struct refine_args *all_args = tasks; struct refine_args *pargs = &all_args[mytask]; - /* Do, er, something. */ + struct image *image = pargs->image; + double nominal_photon_energy = pargs->image->beam->photon_energy; + struct hdfile *hdfile; + struct cpeak *spots; + int n; + double dev; + + 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; + } + + spots = find_intersections(image, image->indexed_cell, &n, 0); + dev = iterate(image, pargs->i_full, pargs->sym, spots, n); + STATUS("Iteration %2i: mean dev = %5.2f\n", 0, dev); + + free(image->data); + if ( image->flags != NULL ) free(image->flags); + hdfile_close(hdfile); + + /* Muppet proofing */ + image->data = NULL; + image->flags = NULL; } @@ -114,7 +310,7 @@ static void integrate_image(int mytask, void *tasks) } /* Figure out which spots should appear in this pattern */ - spots = find_intersections(image, image->indexed_cell, &n, 0); + spots = find_intersections(image, image->indexed_cell, &n, 1); /* For each reflection, estimate the partiality */ for ( j=0; j<n; j++ ) { @@ -414,6 +610,7 @@ int main(int argc, char *argv[]) images[i].bw = beam->bandwidth; images[i].det = det; images[i].beam = beam; + images[i].osf = 1.0; /* Muppet proofing */ images[i].data = NULL; diff --git a/src/image.h b/src/image.h index 6b09dfc4..ae2f3179 100644 --- a/src/image.h +++ b/src/image.h @@ -94,7 +94,6 @@ struct image { /* Information about the crystal */ double m; /* Mosaicity in radians */ - /* Per-shot radiation values */ double lambda; /* Wavelength in m */ double div; /* Divergence in radians */ @@ -102,6 +101,7 @@ struct image { double f0; /* Incident intensity */ int f0_available; /* 0 if f0 wasn't available * from the input. */ + double osf; /* Overall scaling factor */ int width; int height; |