diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/Makefile.am | 4 | ||||
| -rw-r--r-- | src/Makefile.in | 8 | ||||
| -rw-r--r-- | src/facetron.c | 197 | ||||
| -rw-r--r-- | src/post-refinement.c | 213 | ||||
| -rw-r--r-- | src/post-refinement.h | 50 |
5 files changed, 272 insertions, 200 deletions
diff --git a/src/Makefile.am b/src/Makefile.am index dbbdc9a2..5599da98 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -64,7 +64,7 @@ calibrate_detector_LDADD = @LIBS@ facetron_SOURCES = facetron.c cell.c hdf5-file.c utils.c detector.c peaks.c \ image.c geometry.c reflections.c stream.c thread-pool.c \ - beam-parameters.c symmetry.c + beam-parameters.c symmetry.c post-refinement.c facetron_LDADD = @LIBS@ cubeit_SOURCES = cubeit.c cell.c hdf5-file.c utils.c detector.c render.c \ @@ -87,4 +87,4 @@ EXTRA_DIST = cell.h hdf5-file.h image.h utils.h diffraction.h detector.h \ render.h hdfsee.h dirax.h peaks.h index.h filters.h \ diffraction-gpu.h cl-utils.h symmetry.h \ povray.h index-priv.h geometry.h templates.h render_hkl.h \ - stream.h thread-pool.h beam-parameters.h + stream.h thread-pool.h beam-parameters.h post-refinement.h diff --git a/src/Makefile.in b/src/Makefile.in index 00f9af02..fbdcd53c 100644 --- a/src/Makefile.in +++ b/src/Makefile.in @@ -84,7 +84,8 @@ am_facetron_OBJECTS = facetron.$(OBJEXT) cell.$(OBJEXT) \ hdf5-file.$(OBJEXT) utils.$(OBJEXT) detector.$(OBJEXT) \ peaks.$(OBJEXT) image.$(OBJEXT) geometry.$(OBJEXT) \ reflections.$(OBJEXT) stream.$(OBJEXT) thread-pool.$(OBJEXT) \ - beam-parameters.$(OBJEXT) symmetry.$(OBJEXT) + beam-parameters.$(OBJEXT) symmetry.$(OBJEXT) \ + post-refinement.$(OBJEXT) facetron_OBJECTS = $(am_facetron_OBJECTS) facetron_DEPENDENCIES = am_get_hkl_OBJECTS = get_hkl.$(OBJEXT) sfac.$(OBJEXT) cell.$(OBJEXT) \ @@ -324,7 +325,7 @@ calibrate_detector_SOURCES = calibrate_detector.c utils.c hdf5-file.c image.c \ calibrate_detector_LDADD = @LIBS@ facetron_SOURCES = facetron.c cell.c hdf5-file.c utils.c detector.c peaks.c \ image.c geometry.c reflections.c stream.c thread-pool.c \ - beam-parameters.c symmetry.c + beam-parameters.c symmetry.c post-refinement.c facetron_LDADD = @LIBS@ cubeit_SOURCES = cubeit.c cell.c hdf5-file.c utils.c detector.c render.c \ @@ -345,7 +346,7 @@ EXTRA_DIST = cell.h hdf5-file.h image.h utils.h diffraction.h detector.h \ render.h hdfsee.h dirax.h peaks.h index.h filters.h \ diffraction-gpu.h cl-utils.h symmetry.h \ povray.h index-priv.h geometry.h templates.h render_hkl.h \ - stream.h thread-pool.h beam-parameters.h + stream.h thread-pool.h beam-parameters.h post-refinement.h all: all-am @@ -491,6 +492,7 @@ distclean-compile: @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/indexamajig.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/pattern_sim.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/peaks.Po@am__quote@ +@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/post-refinement.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/povray.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/powder_plot.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/process_hkl.Po@am__quote@ diff --git a/src/facetron.c b/src/facetron.c index cf4f099e..9ea7bcea 100644 --- a/src/facetron.c +++ b/src/facetron.c @@ -22,9 +22,6 @@ #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" @@ -35,18 +32,12 @@ #include "peaks.h" #include "thread-pool.h" #include "beam-parameters.h" +#include "post-refinement.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, - NUM_PARAMS -}; - static void show_help(const char *s) { @@ -83,190 +74,6 @@ 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) -{ - double ds; - double nom, den; - - ds = 2.0 * resolution(image->indexed_cell, spot.h, spot.k, spot.l); - - switch ( k ) { - - case REF_SCALE : - return I_partial; - - case REF_DIV : - nom = sqrt(2.0) * ds * sin(image->div); - den = sqrt(1.0 - cos(image->div)); - return nom/den; - - } - - 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; - - case REF_DIV : - STATUS("Shifting div by %e\n", shift); - image->div += shift; - break; - - default : - ERROR("No shift 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, FILE *graph) -{ - int h; - double delta_I = 0.0; - - 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; - } - I_partial *= image->osf; - - get_asymm(hind, kind, lind, &ha, &ka, &la, sym); - I_full = lookup_intensity(i_full, ha, ka, la); - delta_I += fabs(I_partial - spots[h].p * I_full); - - 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/spots[h].p, xc, yc, - spots[h].p, I_partial / I_full); - } - - } - - return delta_I / (double)n; -} - - -static double iterate(struct image *image, double *i_full, const char *sym, - struct cpeak **pspots, int *n) -{ - gsl_matrix *M; - gsl_vector *v; - gsl_vector *shifts; - int h, param; - struct cpeak *spots = *pspots; - - M = gsl_matrix_calloc(NUM_PARAMS, NUM_PARAMS); - v = gsl_vector_calloc(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; - } - I_partial *= image->osf; - - 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 ( param=0; param<NUM_PARAMS; param++ ) { - double shift = gsl_vector_get(shifts, param); - apply_shift(image, param, shift); - } - - gsl_matrix_free(M); - gsl_vector_free(v); - gsl_vector_free(shifts); - - free(spots); - spots = find_intersections(image, image->indexed_cell, n, 0); - *pspots = spots; - return mean_partial_dev(image, spots, *n, sym, i_full, NULL); -} - - static void refine_image(int mytask, void *tasks) { struct refine_args *all_args = tasks; @@ -299,7 +106,7 @@ static void refine_image(int mytask, void *tasks) i = 0; do { last_dev = dev; - dev = iterate(image, pargs->i_full, pargs->sym, &spots, &n); + dev = pr_iterate(image, pargs->i_full, pargs->sym, &spots, &n); STATUS("Iteration %2i: mean dev = %5.2f\n", i, dev); i++; } while ( (fabs(last_dev - dev) > 1.0) && (i < MAX_CYCLES) ); diff --git a/src/post-refinement.c b/src/post-refinement.c new file mode 100644 index 00000000..5ca0dc98 --- /dev/null +++ b/src/post-refinement.c @@ -0,0 +1,213 @@ +/* + * post-refinement.c + * + * Post refinement + * + * (c) 2006-2010 Thomas White <taw@physics.org> + * + * Part of CrystFEL - crystallography with a FEL + * + */ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + + +#include <stdlib.h> +#include <gsl/gsl_poly.h> +#include <assert.h> +#include <gsl/gsl_matrix.h> +#include <gsl/gsl_vector.h> +#include <gsl/gsl_linalg.h> + +#include "image.h" +#include "post-refinement.h" +#include "peaks.h" +#include "symmetry.h" +#include "geometry.h" + + +/* Return the gradient of parameter 'k' given the current status of 'image'. */ +double gradient(struct image *image, int k, + struct cpeak spot, double I_partial) +{ + double ds; + double nom, den; + + ds = 2.0 * resolution(image->indexed_cell, spot.h, spot.k, spot.l); + + switch ( k ) { + + case REF_SCALE : + return I_partial; + + case REF_DIV : + nom = sqrt(2.0) * ds * sin(image->div); + den = sqrt(1.0 - cos(image->div)); + return nom/den; + + } + + ERROR("No gradient defined for parameter %i\n", k); + abort(); +} + + +/* Apply the given shift to the 'k'th parameter of 'image'. */ +void apply_shift(struct image *image, int k, double shift) +{ + switch ( k ) { + + case REF_SCALE : + image->osf += shift; + break; + + case REF_DIV : + STATUS("Shifting div by %e\n", shift); + image->div += shift; + break; + + default : + ERROR("No shift defined for parameter %i\n", k); + abort(); + + } +} + + +double mean_partial_dev(struct image *image, struct cpeak *spots, int n, + const char *sym, double *i_full, FILE *graph) +{ + int h; + double delta_I = 0.0; + + 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; + } + I_partial *= image->osf; + + get_asymm(hind, kind, lind, &ha, &ka, &la, sym); + I_full = lookup_intensity(i_full, ha, ka, la); + delta_I += fabs(I_partial - spots[h].p * I_full); + + 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/spots[h].p, xc, yc, + spots[h].p, I_partial / I_full); + } + + } + + return delta_I / (double)n; +} + + +/* Perform one cycle of post refinement on 'image' against 'i_full' */ +double pr_iterate(struct image *image, double *i_full, const char *sym, + struct cpeak **pspots, int *n) +{ + gsl_matrix *M; + gsl_vector *v; + gsl_vector *shifts; + int h, param; + struct cpeak *spots = *pspots; + + M = gsl_matrix_calloc(NUM_PARAMS, NUM_PARAMS); + v = gsl_vector_calloc(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; + } + I_partial *= image->osf; + + 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 ( param=0; param<NUM_PARAMS; param++ ) { + double shift = gsl_vector_get(shifts, param); + apply_shift(image, param, shift); + } + + gsl_matrix_free(M); + gsl_vector_free(v); + gsl_vector_free(shifts); + + free(spots); + spots = find_intersections(image, image->indexed_cell, n, 0); + *pspots = spots; + return mean_partial_dev(image, spots, *n, sym, i_full, NULL); +} diff --git a/src/post-refinement.h b/src/post-refinement.h new file mode 100644 index 00000000..1ef7a1fd --- /dev/null +++ b/src/post-refinement.h @@ -0,0 +1,50 @@ +/* + * post-refinement.h + * + * Post refinement + * + * (c) 2006-2010 Thomas White <taw@physics.org> + * + * Part of CrystFEL - crystallography with a FEL + * + */ + +#ifndef POST_REFINEMENT_H +#define POST_REFINEMENT_H + + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + + +#include <stdio.h> + +#include "image.h" + + +/* Refineable parameters */ +enum { + REF_SCALE, + REF_DIV, + NUM_PARAMS +}; + + +/* Return the gradient of parameter 'k' given the current status of 'image'. */ +double gradient(struct image *image, int k, struct cpeak spot, + double I_partial); + +/* 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, struct cpeak *spots, int n, + const char *sym, double *i_full, FILE *graph); + + +double pr_iterate(struct image *image, double *i_full, const char *sym, + struct cpeak **pspots, int *n); + + +#endif /* POST_REFINEMENT_H */ |