Commit 77c97f68 authored by Thomas White's avatar Thomas White Committed by Thomas White
Browse files

More work on scaling

parent 20b8402e
......@@ -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;
}
......@@ -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 */
......@@ -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;
};
......
......@@ -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);
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment