Commit 136d9cc9 authored by Thomas White's avatar Thomas White
Browse files

Work on post refinement

parent e0f621b0
......@@ -23,9 +23,7 @@
#include "peaks.h"
#include "beam-parameters.h"
#include "reflist.h"
#define MAX_CPEAKS (256 * 256)
#include "reflist-utils.h"
static signed int locate_peak(double x, double y, double z, double k,
......@@ -292,3 +290,72 @@ double integrate_all(struct image *image, RefList *reflections)
return itot;
}
/* Calculate partialities and apply them to the image's raw_reflections */
void update_partialities(struct image *image, const char *sym,
int *n_expected, int *n_found, int *n_notfound)
{
Reflection *refl;
RefListIterator *iter;
RefList *predicted;
predicted = find_intersections(image, image->indexed_cell, 0);
for ( refl = first_refl(predicted, &iter);
refl != NULL;
refl = next_refl(refl, iter) )
{
Reflection *peak_in_pattern;
double r1, r2, p, x, y;
signed int h, k, l;
int clamp1, clamp2, scalable;
/* Get predicted indices and location */
get_indices(refl, &h, &k, &l);
get_detector_pos(refl, &x, &y);
if ( n_expected != NULL ) (*n_expected)++;
/* Look for this reflection in the pattern */
peak_in_pattern = find_refl(image->raw_reflections, h, k, l);
if ( peak_in_pattern == NULL ) {
if ( n_notfound != NULL ) (*n_notfound)++;
continue;
}
if ( n_found != NULL ) (*n_found)++;
/* Transfer partiality stuff */
get_partial(refl, &r1, &r2, &p, &clamp1, &clamp2);
set_partial(peak_in_pattern, r1, r2, p, clamp1, clamp2);
/* Transfer detector location */
get_detector_pos(refl, &x, &y);
set_detector_pos(peak_in_pattern, 0.0, x, y);
}
reflist_free(predicted);
}
void update_partialities_and_asymm(struct image *image, const char *sym,
ReflItemList *obs,
int *n_expected, int *n_found,
int *n_notfound)
{
/* Get rid of the old list, about to be replaced */
reflist_free(image->reflections);
image->reflections = NULL;
/* Fill in partialities */
update_partialities(image, sym, n_expected, n_found, n_notfound);
/* Rewrite the reflections with the asymmetric indices
* to get the list used for scaling and post refinement */
image->reflections = asymmetric_indices(image->raw_reflections,
sym, obs);
/* Need these lists to work fast */
optimise_reflist(image->reflections);
}
......@@ -24,5 +24,8 @@ extern RefList *find_intersections(struct image *image, UnitCell *cell,
extern double integrate_all(struct image *image, RefList *reflections);
extern void update_partialities_and_asymm(struct image *image, const char *sym,
ReflItemList *obs,
int *n_expected, int *n_found,
int *n_notfound);
#endif /* GEOMETRY_H */
......@@ -302,6 +302,7 @@ 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);
......
......@@ -87,6 +87,7 @@ typedef struct _imagefeaturelist ImageFeatureList;
* int height;
*
* RefList *reflections;
* RefList *raw_reflections;
*
* ImageFeatureList *features;
* };
......@@ -150,6 +151,9 @@ struct image {
/* Integrated (or about-to-be-integrated) reflections */
RefList *reflections;
/* Raw version of "reflections", e.g. in point group 1 */
RefList *raw_reflections;
/* Detected peaks */
ImageFeatureList *features;
......
......@@ -197,6 +197,7 @@ int main(int argc, char *argv[])
int n_expected = 0;
int n_notfound = 0;
char *cref;
int n_usable_patterns = 0;
/* Long options */
const struct option longopts[] = {
......@@ -315,93 +316,47 @@ int main(int argc, char *argv[])
obs = new_items();
for ( i=0; i<n_total_patterns; i++ ) {
RefList *predicted;
RefList *measured;
Reflection *refl;
RefListIterator *iter;
images[n_usable_patterns].det = NULL;
if ( read_chunk(fh, &images[i]) != 0 ) {
if ( read_chunk(fh, &images[n_usable_patterns]) != 0 ) {
/* Should not happen, because we counted the patterns
* earlier. */
ERROR("Failed to read chunk from the input stream.\n");
return 1;
}
/* FIXME: This leaves gaps in the array */
if ( images[i].indexed_cell == NULL ) continue;
/* Won't be needing this, if it exists */
image_feature_list_free(images[i].features);
images[i].features = NULL;
images[i].div = beam->divergence;
images[i].bw = beam->bandwidth;
images[i].det = det;
images[i].width = det->max_fs;
images[i].height = det->max_ss;
images[i].osf = 1.0;
images[i].profile_radius = 0.005e9;
/* Muppet proofing */
images[i].data = NULL;
images[i].flags = NULL;
images[i].beam = NULL;
/* Calculate initial partialities and fill in intensities from
* the stream */
predicted = find_intersections(&images[i],
images[i].indexed_cell, 0);
/* We start again with a new reflection list, this time with
* the asymmetric indices */
measured = images[i].reflections;
images[i].reflections = reflist_new();
for ( refl = first_refl(predicted, &iter);
refl != NULL;
refl = next_refl(refl, iter) ) {
image_feature_list_free(images[n_usable_patterns].features);
images[n_usable_patterns].features = NULL;
Reflection *peak_in_pattern;
Reflection *new;
signed int h, k, l, ha, ka, la;
double r1, r2, p, x, y;
int clamp1, clamp2;
/* Get predicted indices and location */
get_indices(refl, &h, &k, &l);
get_detector_pos(refl, &x, &y);
n_expected++;
/* Look for this reflection in the pattern */
peak_in_pattern = find_refl(measured, h, k, l);
if ( peak_in_pattern == NULL ) {
n_notfound++;
continue;
}
n_found++;
/* "n_usable_patterns" will not be incremented in this case */
if ( images[n_usable_patterns].indexed_cell == NULL ) continue;
/* Put it into the asymmetric cell */
get_asymm(h, k, l, &ha, &ka, &la, sym);
if ( find_item(obs, ha, ka, la) == 0 ) {
add_item(obs, ha, ka, la);
}
/* Fill in initial estimates of stuff */
images[n_usable_patterns].div = beam->divergence;
images[n_usable_patterns].bw = beam->bandwidth;
images[n_usable_patterns].det = det;
images[n_usable_patterns].width = det->max_fs;
images[n_usable_patterns].height = det->max_ss;
images[n_usable_patterns].osf = 1.0;
images[n_usable_patterns].profile_radius = 0.005e9;
/* Create new reflection and copy data across */
new = add_refl(images[i].reflections, ha, ka, la);
get_partial(refl, &r1, &r2, &p, &clamp1, &clamp2);
get_detector_pos(refl, &x, &y);
set_int(new, get_intensity(peak_in_pattern));
set_partial(new, r1, r2, p, clamp1, clamp2);
set_detector_pos(new, 0.0, x, y);
/* Muppet proofing */
images[n_usable_patterns].data = NULL;
images[n_usable_patterns].flags = NULL;
images[n_usable_patterns].beam = NULL;
}
reflist_free(measured);
reflist_free(predicted);
/* This is the raw list of reflections */
images[n_usable_patterns].raw_reflections =
images[n_usable_patterns].reflections;
images[n_usable_patterns].reflections = NULL;
/* Do magic on the reflection list to make things go faster */
optimise_reflist(images[i].reflections);
update_partialities_and_asymm(&images[n_usable_patterns], sym,
obs, &n_expected, &n_found,
&n_notfound);
progress_bar(i, n_total_patterns-1, "Loading pattern data");
n_usable_patterns++;
}
fclose(fh);
......@@ -410,12 +365,12 @@ int main(int argc, char *argv[])
STATUS("Mean measurements per unique reflection: %5.2f\n",
(double)n_found / num_items(obs));
cref = find_common_reflections(images, n_total_patterns);
cref = find_common_reflections(images, n_usable_patterns);
/* Make initial estimates */
STATUS("Performing initial scaling.\n");
select_scalable_reflections(images, n_total_patterns);
full = scale_intensities(images, n_total_patterns, sym, obs, cref);
full = scale_intensities(images, n_usable_patterns, sym, obs, cref);
/* Iterate */
for ( i=0; i<n_iter; i++ ) {
......@@ -446,8 +401,8 @@ int main(int argc, char *argv[])
/* Re-estimate all the full intensities */
reflist_free(full);
select_scalable_reflections(images, n_total_patterns);
full = scale_intensities(images, n_total_patterns,
select_scalable_reflections(images, n_usable_patterns);
full = scale_intensities(images, n_usable_patterns,
sym, obs, cref);
fclose(fhg);
......@@ -456,7 +411,7 @@ int main(int argc, char *argv[])
}
STATUS("Final scale factors:\n");
for ( i=0; i<n_total_patterns; i++ ) {
for ( i=0; i<n_usable_patterns; i++ ) {
STATUS("%4i : %5.2f\n", i, images[i].osf);
}
......@@ -464,18 +419,18 @@ int main(int argc, char *argv[])
write_reflist(outfile, full, images[0].indexed_cell);
/* Clean up */
for ( i=0; i<n_total_patterns; i++ ) {
for ( i=0; i<n_usable_patterns; i++ ) {
reflist_free(images[i].reflections);
reflist_free(images[i].raw_reflections);
}
reflist_free(full);
delete_items(obs);
free(sym);
free(outfile);
free(det->panels);
free(det);
free_detector_geometry(det);
free(beam);
free(cref);
for ( i=0; i<n_total_patterns; i++ ) {
for ( i=0; i<n_usable_patterns; i++ ) {
cell_free(images[i].indexed_cell);
free(images[i].filename);
}
......
......@@ -11,6 +11,7 @@
#include <stdio.h>
#include <assert.h>
#include "reflist.h"
......@@ -359,3 +360,40 @@ RefList *read_reflections(const char *filename)
return out;
}
RefList *asymmetric_indices(RefList *in, const char *sym, ReflItemList *obs)
{
Reflection *refl;
RefListIterator *iter;
RefList *new;
new = reflist_new();
for ( refl = first_refl(in, &iter);
refl != NULL;
refl = next_refl(refl, iter) ) {
signed int h, k, l;
signed int ha, ka, la;
Reflection *cr;
get_indices(refl, &h, &k, &l);
get_asymm(h, k, l, &ha, &ka, &la, sym);
cr = add_refl(new, ha, ka, la);
assert(cr != NULL);
copy_data(cr, refl);
if ( obs != NULL ) {
if ( !find_item(obs, ha, ka, la) ) {
add_item(obs, ha, ka, la);
}
}
}
return new;
}
......@@ -41,4 +41,7 @@ extern int find_equiv_in_list(RefList *list, signed int h, signed int k,
signed int l, const char *sym, signed int *hu,
signed int *ku, signed int *lu);
extern RefList *asymmetric_indices(RefList *in, const char *sym,
ReflItemList *obs);
#endif /* REFLIST_UTILS_H */
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