aboutsummaryrefslogtreecommitdiff
path: root/src/hrs-scaling.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/hrs-scaling.c')
-rw-r--r--src/hrs-scaling.c132
1 files changed, 132 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
+}