aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--src/hrs-scaling.c292
-rw-r--r--src/hrs-scaling.h3
-rw-r--r--src/image.h4
-rw-r--r--src/partialator.c40
4 files changed, 273 insertions, 66 deletions
diff --git a/src/hrs-scaling.c b/src/hrs-scaling.c
index 83c18e17..11306029 100644
--- a/src/hrs-scaling.c
+++ b/src/hrs-scaling.c
@@ -27,6 +27,10 @@
#include "cell.h"
+/* Maximum number of iterations of NLSq scaling per macrocycle. */
+#define MAX_CYCLES (10)
+
+
static void show_matrix_eqn(gsl_matrix *M, gsl_vector *v, int r)
{
int i, j;
@@ -41,92 +45,276 @@ static void show_matrix_eqn(gsl_matrix *M, gsl_vector *v, int r)
}
-/* Scale the stack of images */
-void scale_intensities(struct image *images, int n, const char *sym)
+static double gradient(int j, signed int h, signed int k, signed int l,
+ struct image *images, int n, const char *sym)
+{
+ struct image *image;
+ struct cpeak *spots;
+ double num1 = 0.0;
+ double num2 = 0.0;
+ double den = 0.0;
+ double num1_this = 0.0;
+ double num2_this = 0.0;
+ int m, i;
+
+ /* "Everything" terms */
+ for ( m=0; m<n; m++ ) {
+
+ image = &images[m];
+ spots = image->cpeaks;
+
+ for ( i=0; i<image->n_cpeaks; i++ ) {
+
+ signed int ha, ka, la;
+
+ if ( !spots[i].valid ) continue;
+ if ( spots[i].p < 0.1 ) continue;
+
+ get_asymm(spots[i].h, spots[i].k, spots[i].l,
+ &ha, &ka, &la, sym);
+
+ if ( ha != h ) continue;
+ if ( ka != k ) continue;
+ if ( la != l ) continue;
+
+ num1 += pow(image->osf, 2.0) * pow(spots[i].p, 2.0);
+ num2 += image->osf * spots[i].intensity * spots[i].p;
+ den += pow(image->osf, 2.0) * pow(spots[i].p, 2.0);
+
+ }
+
+ }
+
+ /* "This frame" terms */
+ image = &images[j];
+ spots = image->cpeaks;
+ for ( i=0; i<image->n_cpeaks; i++ ) {
+
+ signed int ha, ka, la;
+
+ if ( !spots[i].valid ) continue;
+ if ( spots[i].p < 0.1 ) continue;
+
+ get_asymm(spots[i].h, spots[i].k, spots[i].l,
+ &ha, &ka, &la, sym);
+
+ if ( ha != h ) continue;
+ if ( ka != k ) continue;
+ if ( la != l ) continue;
+
+ num1_this += spots[i].intensity * spots[i].p;
+ num2_this += pow(spots[i].p, 2.0);
+
+ }
+
+ return (num1*num1_this - num2*num2_this) / den;
+}
+
+
+static double iterate_scale(struct image *images, int n, const char *sym,
+ double *I_full_old)
{
-#if 0
gsl_matrix *M;
gsl_vector *v;
gsl_vector *shifts;
- int j;
+ int m;
+ double max_shift;
+ int crossed = 0;
- M = gsl_matrix_calloc(n, n);
- v = gsl_vector_calloc(n);
+ M = gsl_matrix_calloc(n-1, n-1);
+ v = gsl_vector_calloc(n-1);
- for ( j=0; j<n; j++ ) {
+ for ( m=0; m<n; m++ ) { /* "Equation number": one equation per frame */
- signed int hind, kind, lind;
- signed int ha, ka, la;
- double I_full, delta_I;
- float I_partial;
- float xc, yc;
+ int k; /* Frame (scale factor) number */
int h;
- struct image *image = &images[j];
+ int mcomp;
+ double vc_tot = 0.0;
+ struct image *image = &images[m];
struct cpeak *spots = image->cpeaks;
- for ( h=0; h<image->n_cpeaks; h++ ) {
+ if ( m == crossed ) continue;
- int g;
- double v_c, gr, I_full;
+ /* Determine the "solution" vector component */
+ for ( h=0; h<image->n_cpeaks; h++ ) {
- hind = spots[h].h;
- kind = spots[h].k;
- lind = spots[h].l;
+ double v_c;
+ double ic, ip;
+ signed int ha, ka, la;
- /* Don't attempt to use spots with very small
- * partialities, since it won't be accurate. */
+ if ( !spots[h].valid ) continue;
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(spots[h].h, spots[h].k, spots[h].l,
+ &ha, &ka, &la, sym);
+ ic = lookup_intensity(I_full_old, ha, ka, la);
+
+ v_c = ip - (spots[h].p * image->osf * ic);
+ v_c *= pow(spots[h].p, 2.0);
+ v_c *= pow(image->osf, 2.0);
+ v_c *= ic;
+ v_c *= gradient(m, ha, ka, la, images, n, sym);
+ vc_tot += v_c;
+
+ }
+
+ mcomp = m;
+ if ( m > crossed ) mcomp--;
+ gsl_vector_set(v, mcomp, vc_tot);
- 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);
+ /* Now fill in the matrix components */
+ for ( k=0; k<n; k++ ) {
+ double val = 0.0;
+ int kcomp;
- for ( g=0; g<NUM_PARAMS; g++ ) {
+ if ( k == crossed ) continue;
- double M_curr, M_c;
+ for ( h=0; h<image->n_cpeaks; h++ ) {
- M_curr = gsl_matrix_get(M, g, k);
+ signed int ha, ka, la;
+ double con;
+ double ic;
- 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);
+ if ( !spots[h].valid ) continue;
+ if ( spots[h].p < 0.1 ) continue;
- gsl_matrix_set(M, g, k, M_curr + M_c);
+ get_asymm(spots[h].h, spots[h].k, spots[h].l,
+ &ha, &ka, &la, sym);
+ ic = lookup_intensity(I_full_old, ha, ka, la);
+
+ con = -pow(spots[h].p, 3.0);
+ con *= pow(image->osf, 3.0);
+ con *= ic;
+ con *= gradient(m, ha, ka, la, images, n, sym);
+ con *= gradient(k, ha, ka, la, images, n, sym);
+ val += con;
}
- gr = gradient(image, k, spots[h],
- image->profile_radius);
- v_c = delta_I * I_full * gr;
- gsl_vector_set(v, k, v_c);
+ kcomp = k;
+ if ( k > crossed ) kcomp--;
+ gsl_matrix_set(M, mcomp, kcomp, val);
}
+
}
- show_matrix_eqn(M, v, NUM_PARAMS);
+ show_matrix_eqn(M, v, n-1);
- shifts = gsl_vector_alloc(NUM_PARAMS);
+ shifts = gsl_vector_alloc(n-1);
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);
+ max_shift = 0.0;
+ for ( m=0; m<n-1; m++ ) {
+
+ double shift = gsl_vector_get(shifts, m);
+ int mimg;
+
+ mimg = m;
+ if ( mimg >= crossed ) mimg++;
+
+ images[mimg].osf += shift;
+
+ if ( shift > max_shift ) max_shift = shift;
+
}
+ gsl_vector_free(shifts);
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
+ return max_shift;
+}
+
+
+static double *lsq_intensities(struct image *images, int n, const char *sym,
+ ReflItemList *obs)
+{
+ double *I_full;
+ int i;
+
+ I_full = new_list_intensity();
+ for ( i=0; i<num_items(obs); i++ ) {
+
+ signed int h, k, l;
+ struct refl_item *it = get_item(obs, i);
+ double num = 0.0;
+ double den = 0.0;
+ int m;
+
+ get_asymm(it->h, it->k, it->l, &h, &k, &l, sym);
+
+ /* For each frame */
+ for ( m=0; m<n; m++ ) {
+
+ double G;
+ int a;
+
+ G = images[m].osf;
+
+ /* For each peak */
+ for ( a=0; a<images[m].n_cpeaks; a++ ) {
+
+ signed int ha, ka, la;
+
+ if ( !images[m].cpeaks[a].valid ) continue;
+ if ( images[m].cpeaks[a].p < 0.1 ) continue;
+
+ /* Correct reflection? */
+ get_asymm(images[m].cpeaks[a].h,
+ images[m].cpeaks[a].k,
+ images[m].cpeaks[a].l,
+ &ha, &ka, &la, sym);
+ if ( ha != h ) continue;
+ if ( ka != k ) continue;
+ if ( la != l ) continue;
+
+ num += images[m].cpeaks[a].intensity
+ * images[m].cpeaks[a].p * G;
+
+ den += pow(images[m].cpeaks[a].p, 2.0)
+ * pow(G, 2.0);
+
+ }
+
+ }
+
+ set_intensity(I_full, h, k, l, num/den);
+
+ }
+
+ return I_full;
+}
+
+
+/* Scale the stack of images */
+double *scale_intensities(struct image *images, int n, const char *sym,
+ ReflItemList *obs)
+{
+ int m;
+ double *I_full;
+ int i;
+ double max_shift;
+
+ /* Start with all scale factors equal */
+ for ( m=0; m<n; m++ ) {
+ images[m].osf = 1.0;
+ }
+
+ /* Calculate LSQ intensities using these scale factors */
+ I_full = lsq_intensities(images, n, sym, obs);
+
+ /* Iterate */
+ i = 0;
+ do {
+
+ max_shift = iterate_scale(images, n, sym, I_full);
+ STATUS("Iteration %2i: max shift = %5.2f\n", i, max_shift);
+ free(I_full);
+ I_full = lsq_intensities(images, n, sym, obs);
+ i++;
+
+ } while ( (max_shift > 0.1) && (i < MAX_CYCLES) );
+
+ return I_full;
}
diff --git a/src/hrs-scaling.h b/src/hrs-scaling.h
index 28343283..ae384aa8 100644
--- a/src/hrs-scaling.h
+++ b/src/hrs-scaling.h
@@ -20,7 +20,8 @@
#include "image.h"
-extern void scale_intensities(struct image *image, int n, const char *sym);
+extern double *scale_intensities(struct image *image, int n, const char *sym,
+ ReflItemList *obs);
#endif /* HRS_SCALING_H */
diff --git a/src/image.h b/src/image.h
index 0d2f6c2e..078231a0 100644
--- a/src/image.h
+++ b/src/image.h
@@ -61,6 +61,7 @@ struct cpeak
signed int l;
double min_distance;
+ int valid;
/* Partiality */
double r1; /* First excitation error */
@@ -72,6 +73,9 @@ struct cpeak
/* Location in image */
int x;
int y;
+
+ /* Intensity */
+ double intensity;
};
diff --git a/src/partialator.c b/src/partialator.c
index c023fcde..b90c8c27 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -163,7 +163,7 @@ static void refine_all(struct image *images, int n_total_patterns,
}
-static void integrate_image(struct image *image)
+static void integrate_image(struct image *image, ReflItemList *obs)
{
struct cpeak *spots;
int j, n;
@@ -210,15 +210,22 @@ static void integrate_image(struct image *image)
if ( integrate_peak(image, spots[j].x, spots[j].y,
&xc, &yc, &i_partial, NULL, NULL,
1, 1, 0) ) {
+ spots[j].valid = 0;
continue;
}
+ spots[j].intensity = i_partial;
+ spots[j].valid = 1;
+
+ if ( !find_item(obs, h, k, l) ) add_item(obs, h, k, l);
+
}
+ image->cpeaks = spots;
+ image->n_cpeaks = n;
free(image->data);
if ( image->flags != NULL ) free(image->flags);
hdfile_close(hdfile);
- image->cpeaks = spots;
/* Muppet proofing */
image->data = NULL;
@@ -239,7 +246,6 @@ int main(int argc, char *argv[])
int config_basename = 0;
int config_checkprefix = 1;
struct detector *det;
- double *i_full;
unsigned int *cts;
ReflItemList *obs;
int i;
@@ -247,6 +253,7 @@ int main(int argc, char *argv[])
struct image *images;
int n_iter = 10;
struct beam_params *beam = NULL;
+ double *I_full;
/* Long options */
const struct option longopts[] = {
@@ -363,10 +370,6 @@ int main(int argc, char *argv[])
return 1;
}
- /* Prepare for iteration */
- i_full = new_list_intensity();
- obs = new_items();
-
n_total_patterns = count_patterns(fh);
STATUS("There are %i patterns to process\n", n_total_patterns);
@@ -378,6 +381,7 @@ int main(int argc, char *argv[])
/* Fill in what we know about the images so far */
rewind(fh);
+ obs = new_items();
for ( i=0; i<n_total_patterns; i++ ) {
UnitCell *cell;
@@ -416,7 +420,7 @@ int main(int argc, char *argv[])
images[i].flags = NULL;
/* Get reflections from this image */
- integrate_image(&images[i]);
+ integrate_image(&images[i], obs);
progress_bar(i, n_total_patterns-1, "Loading pattern data");
@@ -428,7 +432,7 @@ int main(int argc, char *argv[])
/* Make initial estimates */
STATUS("Performing initial scaling.\n");
- scale_intensities(images, n_total_patterns, sym);
+ I_full = scale_intensities(images, n_total_patterns, sym, obs);
/* Iterate */
for ( i=0; i<n_iter; i++ ) {
@@ -454,21 +458,31 @@ int main(int argc, char *argv[])
}
/* Refine the geometry of all patterns to get the best fit */
- refine_all(images, n_total_patterns, det, sym, obs, i_full,
+ refine_all(images, n_total_patterns, det, sym, obs, I_full,
nthreads, fhg, fhp);
/* Re-estimate all the full intensities */
- scale_intensities(images, n_total_patterns, sym);
+ free(I_full);
+ I_full = scale_intensities(images, n_total_patterns, sym, obs);
fclose(fhg);
fclose(fhp);
+
+ }
+
+ STATUS("Final scale factors:\n");
+ for ( i=0; i<n_total_patterns; i++ ) {
+ STATUS("%4i : %e\n", i, images[i].osf);
}
/* Output results */
- write_reflections(outfile, obs, i_full, NULL, NULL, cts, NULL);
+ write_reflections(outfile, obs, I_full, NULL, NULL, cts, NULL);
/* Clean up */
- free(i_full);
+ for ( i=0; i<n_total_patterns; i++ ) {
+ free(images[i].cpeaks);
+ }
+ free(I_full);
delete_items(obs);
free(sym);
free(outfile);
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-18 05:31:42 +0000