Commit 8eab85fe authored by Thomas White's avatar Thomas White
Browse files

Use reference reflections for scaling

parent 29db694f
......@@ -135,8 +135,62 @@ static double s_vh(struct image *images, int n,
}
static gsl_vector *solve_by_eigenvalue_filtration(gsl_vector *v, gsl_matrix *M)
{
gsl_eigen_symmv_workspace *work;
gsl_vector *e_val;
gsl_matrix *e_vec;
int val, n, frame;
gsl_vector *shifts;
n = v->size;
if ( v->size != M->size1 ) return NULL;
if ( v->size != M->size2 ) return NULL;
/* Diagonalise */
work = gsl_eigen_symmv_alloc(n);
e_val = gsl_vector_alloc(n);
e_vec = gsl_matrix_alloc(n, n);
val = gsl_eigen_symmv(M, e_val, e_vec, work);
if ( val ) {
ERROR("Couldn't diagonalise matrix.\n");
return NULL;
}
gsl_eigen_symmv_free(work);
val = gsl_eigen_symmv_sort(e_val, e_vec, GSL_EIGEN_SORT_ABS_DESC);
/* Rotate the "solution vector" */
gsl_vector *rprime;
rprime = gsl_vector_alloc(n);
val = gsl_blas_dgemv(CblasTrans, 1.0, e_vec, v, 0.0, rprime);
/* Solve (now easy) */
gsl_vector *sprime;
sprime = gsl_vector_alloc(n);
for ( frame=0; frame<n-1; frame++ ) {
double num, den;
num = gsl_vector_get(rprime, frame);
den = gsl_vector_get(e_val, frame);
gsl_vector_set(sprime, frame, num/den);
}
gsl_vector_set(sprime, n-1, 0.0); /* Set last shift to zero */
/* Rotate back */
shifts = gsl_vector_alloc(n);
val = gsl_blas_dgemv(CblasNoTrans, 1.0, e_vec, sprime, 0.0, shifts);
gsl_vector_free(sprime);
gsl_vector_free(rprime);
gsl_matrix_free(e_vec);
gsl_vector_free(e_val);
return shifts;
}
static double iterate_scale(struct image *images, int n,
ReflItemList *obs, const char *sym, char *cref)
ReflItemList *obs, const char *sym, char *cref,
double *reference)
{
gsl_matrix *M;
gsl_vector *v;
......@@ -197,16 +251,26 @@ static double iterate_scale(struct image *images, int n,
int b; /* Frame (scale factor) number */
struct image *image_a = &images[a];
double vc, Ih, uh, vh, rha, vha, uha;
double vc, Ih, uh, rha, vha, uha;
double vval;
/* Determine the "solution" vector component */
uha = uha_arr[a];
vha = vha_arr[a];
uh = lookup_intensity(uh_arr, h, k, l);
vh = lookup_intensity(vh_arr, h, k, l);
Ih = vh / uh;
if ( isnan(Ih) ) Ih = 0.0; /* 0 / 0 = 0, not NaN */
if ( !reference ) {
double vh;
vh = lookup_intensity(vh_arr, h, k, l);
Ih = vh / uh;
/* 0 / 0 = 0, not NaN */
if ( isnan(Ih) ) Ih = 0.0;
} else {
/* Look up by asymmetric indices */
Ih = lookup_intensity(reference, h, k, l);
}
rha = vha - image_a->osf * uha * Ih;
vc = Ih * rha;
......@@ -220,6 +284,9 @@ static double iterate_scale(struct image *images, int n,
/* Matrix is symmetric */
if ( b > a ) continue;
/* Off-diagonals zero if reference available */
if ( (reference != NULL) && (a!=b) ) continue;
/* Zero if no common reflections */
if ( cref[a+n*b] != 0 ) continue;
......@@ -252,43 +319,24 @@ static double iterate_scale(struct image *images, int n,
free(uha_arr);
free(vha_arr);
/* Fox and Holmes method */
gsl_eigen_symmv_workspace *work;
gsl_vector *e_val;
gsl_matrix *e_vec;
int val;
/* Diagonalise */
work = gsl_eigen_symmv_alloc(n);
e_val = gsl_vector_alloc(n);
e_vec = gsl_matrix_alloc(n, n);
val = gsl_eigen_symmv(M, e_val, e_vec, work);
if ( val ) {
ERROR("Couldn't diagonalise matrix.\n");
return 0.0;
show_matrix_eqn(M, v, n);
if ( reference == NULL ) {
shifts = solve_by_eigenvalue_filtration(v, M);
} else {
shifts = gsl_vector_alloc(n);
for ( frame=0; frame<n-1; frame++ ) {
double num, den;
num = gsl_vector_get(v, frame);
den = gsl_matrix_get(M, frame, frame);
gsl_vector_set(shifts, frame, num/den);
}
}
gsl_eigen_symmv_free(work);
val = gsl_eigen_symmv_sort(e_val, e_vec, GSL_EIGEN_SORT_ABS_DESC);
/* Rotate the "solution vector" */
gsl_vector *rprime;
rprime = gsl_vector_alloc(n);
val = gsl_blas_dgemv(CblasTrans, 1.0, e_vec, v, 0.0, rprime);
/* Solve (now easy) */
gsl_vector *sprime;
sprime = gsl_vector_alloc(n);
for ( frame=0; frame<n-1; frame++ ) {
double num, den;
num = gsl_vector_get(rprime, frame);
den = gsl_vector_get(e_val, frame);
gsl_vector_set(sprime, frame, num/den);
}
gsl_vector_set(sprime, n-1, 0.0); /* Set last shift to zero */
gsl_matrix_free(M);
gsl_vector_free(v);
/* Rotate back */
shifts = gsl_vector_alloc(n);
val = gsl_blas_dgemv(CblasNoTrans, 1.0, e_vec, sprime, 0.0, shifts);
if ( shifts == NULL ) return 0.0;
/* Apply shifts */
max_shift = 0.0;
......@@ -305,12 +353,6 @@ static double iterate_scale(struct image *images, int n,
}
gsl_vector_free(shifts);
gsl_vector_free(sprime);
gsl_vector_free(rprime);
gsl_matrix_free(e_vec);
gsl_vector_free(e_val);
gsl_matrix_free(M);
gsl_vector_free(v);
return max_shift;
}
......@@ -394,7 +436,7 @@ static void normalise_osfs(struct image *images, int n)
/* Scale the stack of images */
RefList *scale_intensities(struct image *images, int n, const char *sym,
ReflItemList *obs, char *cref)
ReflItemList *obs, char *cref, double *reference)
{
int m;
RefList *full;
......@@ -408,7 +450,7 @@ RefList *scale_intensities(struct image *images, int n, const char *sym,
i = 0;
do {
max_shift = iterate_scale(images, n, obs, sym, cref);
max_shift = iterate_scale(images, n, obs, sym, cref, reference);
STATUS("Scaling iteration %2i: max shift = %5.2f\n",
i, max_shift);
i++;
......
......@@ -21,7 +21,8 @@
#include "image.h"
extern RefList *scale_intensities(struct image *images, int n, const char *sym,
ReflItemList *obs, char *cref);
ReflItemList *obs, char *cref,
double *reference);
extern char *find_common_reflections(struct image *images, int n);
......
......@@ -171,7 +171,7 @@ int main(int argc, char *argv[])
char *cref;
int n_usable_patterns = 0;
char *reference_file = NULL;
RefList *reference;
double *reference = NULL;
/* Long options */
const struct option longopts[] = {
......@@ -230,7 +230,7 @@ int main(int argc, char *argv[])
break;
case 1 :
reference = strdup(optarg);
reference_file = strdup(optarg);
break;
case 0 :
......@@ -278,9 +278,15 @@ int main(int argc, char *argv[])
}
if ( reference_file != NULL ) {
reference = read_reflections(reference_file);
RefList *list;
RefList *symmed;
list = read_reflections(reference_file);
free(reference_file);
if ( reference == NULL ) return 1;
if ( list == NULL ) return 1;
symmed = asymmetric_indices(list, sym);
reflist_free(list);
reference = intensities_from_list(symmed);
reflist_free(symmed);
}
n_total_patterns = count_patterns(fh);
......@@ -360,7 +366,7 @@ int main(int argc, char *argv[])
/* Make initial estimates */
STATUS("Performing initial scaling.\n");
full = scale_intensities(images, n_usable_patterns, sym,
scalable, cref);
scalable, cref, reference);
for ( i=0; i<num_items(scalable); i++ ) {
Reflection *f;
......@@ -425,7 +431,7 @@ int main(int argc, char *argv[])
/* Re-estimate all the full intensities */
reflist_free(full);
full = scale_intensities(images, n_usable_patterns,
sym, scalable, cref);
sym, scalable, cref, reference);
fclose(fhg);
fclose(fhp);
......@@ -451,6 +457,7 @@ int main(int argc, char *argv[])
free_detector_geometry(det);
free(beam);
free(cref);
free(reference);
for ( i=0; i<n_usable_patterns; i++ ) {
cell_free(images[i].indexed_cell);
free(images[i].filename);
......
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