diff options
| -rw-r--r-- | src/hrs-scaling.c | 132 | ||||
| -rw-r--r-- | src/hrs-scaling.h | 26 |
2 files changed, 158 insertions, 0 deletions
diff --git a/src/hrs-scaling.c b/src/hrs-scaling.c new file mode 100644 index 00000000..83c18e17 --- /dev/null +++ b/src/hrs-scaling.c @@ -0,0 +1,132 @@ +/* + * hrs-scaling.c + * + * Intensity scaling using generalised HRS target function + * + * (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 <assert.h> +#include <gsl/gsl_matrix.h> +#include <gsl/gsl_vector.h> +#include <gsl/gsl_linalg.h> + +#include "image.h" +#include "peaks.h" +#include "symmetry.h" +#include "geometry.h" +#include "cell.h" + + +static void show_matrix_eqn(gsl_matrix *M, gsl_vector *v, int r) +{ + int i, j; + + for ( i=0; i<r; i++ ) { + STATUS("[ "); + for ( j=0; j<r; j++ ) { + STATUS("%+9.3e ", gsl_matrix_get(M, i, j)); + } + STATUS("][ a%2i ] = [ %+9.3e ]\n", i, gsl_vector_get(v, i)); + } +} + + +/* Scale the stack of images */ +void scale_intensities(struct image *images, int n, const char *sym) +{ +#if 0 + gsl_matrix *M; + gsl_vector *v; + gsl_vector *shifts; + int j; + + M = gsl_matrix_calloc(n, n); + v = gsl_vector_calloc(n); + + for ( j=0; j<n; j++ ) { + + signed int hind, kind, lind; + signed int ha, ka, la; + double I_full, delta_I; + float I_partial; + float xc, yc; + int h; + struct image *image = &images[j]; + struct cpeak *spots = image->cpeaks; + + for ( h=0; h<image->n_cpeaks; h++ ) { + + int g; + double v_c, gr, I_full; + + 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; + + 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 / image->osf); + + + 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], + image->profile_radius) + * gradient(image, k, spots[h], + image->profile_radius); + M_c *= pow(I_full, 2.0); + + gsl_matrix_set(M, g, k, M_curr + M_c); + + } + + gr = gradient(image, k, spots[h], + image->profile_radius); + v_c = delta_I * I_full * gr; + gsl_vector_set(v, k, v_c); + + } + + } + show_matrix_eqn(M, v, NUM_PARAMS); + + 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); +#endif +} diff --git a/src/hrs-scaling.h b/src/hrs-scaling.h new file mode 100644 index 00000000..28343283 --- /dev/null +++ b/src/hrs-scaling.h @@ -0,0 +1,26 @@ +/* + * hrs-scaling.h + * + * Intensity scaling using generalised HRS target function + * + * (c) 2006-2010 Thomas White <taw@physics.org> + * + * Part of CrystFEL - crystallography with a FEL + * + */ + +#ifndef HRS_SCALING_H +#define HRS_SCALING_H + + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + + +#include "image.h" + +extern void scale_intensities(struct image *image, int n, const char *sym); + + +#endif /* HRS_SCALING_H */ |