From 31851dfe017e89ebcf00c89ea4c3a0f12f7e86a6 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Mon, 2 May 2011 16:09:12 +0200 Subject: Get rid of "image.raw_reflections" --- src/geometry.c | 80 ++++++++++++++++++++++++++++----------------------- src/geometry.h | 8 +++--- src/image.h | 4 --- src/partialator.c | 33 +++++++++++---------- src/post-refinement.c | 41 ++++++++++++-------------- src/reflist-utils.c | 9 ++---- src/reflist-utils.h | 3 +- src/reflist.c | 37 ++++++++++++++++++++++++ src/reflist.h | 5 ++++ 9 files changed, 130 insertions(+), 90 deletions(-) diff --git a/src/geometry.c b/src/geometry.c index db261e6b..d5845e5d 100644 --- a/src/geometry.c +++ b/src/geometry.c @@ -24,6 +24,7 @@ #include "beam-parameters.h" #include "reflist.h" #include "reflist-utils.h" +#include "symmetry.h" static signed int locate_peak(double x, double y, double z, double k, @@ -293,10 +294,8 @@ double integrate_all(struct image *image, RefList *reflections) /* Decide which reflections can be scaled */ -static int select_scalable_reflections(RefList *list) +static void select_scalable_reflections(RefList *list, ReflItemList *sc_l) { - int n_scalable = 0; - Reflection *refl; RefListIterator *iter; @@ -312,16 +311,27 @@ static int select_scalable_reflections(RefList *list) if ( v < 0.1 ) scalable = 0; set_scalable(refl, scalable); - if ( scalable ) n_scalable++; + if ( scalable && (sc_l != NULL) ) { - } + signed int h, k, l; + + get_indices(refl, &h, &k, &l); /* Should already be + * asymmetric */ + if ( !find_item(sc_l, h, k, l) ) { + add_item(sc_l, h, k, l); + } + + } - return n_scalable; + } } -/* Calculate partialities and apply them to the image's raw_reflections */ +/* Calculate partialities and apply them to the image's raw_reflections, + * returning a ReflItemList of the currentl scalable (asymmetric) reflections. + */ void update_partialities(struct image *image, const char *sym, + ReflItemList *scalable, int *n_expected, int *n_found, int *n_notfound) { Reflection *refl; @@ -335,56 +345,54 @@ void update_partialities(struct image *image, const char *sym, refl = next_refl(refl, iter) ) { - Reflection *peak_in_pattern; + Reflection *p_peak; double r1, r2, p, x, y; signed int h, k, l; + signed int ha, ka, la; int clamp1, clamp2; + int found = 0; /* Get predicted indices and location */ get_indices(refl, &h, &k, &l); get_detector_pos(refl, &x, &y); if ( n_expected != NULL ) (*n_expected)++; + /* Get the asymmetric indices, with which the reflection in the + * image's list will be indexed. */ + get_asymm(h, k, l, &ha, &ka, &la, sym); + /* 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)++; + p_peak = find_refl(image->reflections, ha, ka, la); + do { + + signed int hs, ks, ls; + + if ( p_peak != NULL ) { + get_symmetric_indices(p_peak, &hs, &ks, &ls); + if ( (hs==h) && (ks==k) && (ls==l) ) found = 1; + } + + if ( !found && (p_peak != NULL ) ) { + p_peak = next_found_refl(p_peak); + } + + } while ( !found && (p_peak != NULL) ); + if ( !found ) { + 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); + set_partial(p_peak, r1, r2, p, clamp1, clamp2); /* Transfer detector location */ get_detector_pos(refl, &x, &y); - set_detector_pos(peak_in_pattern, 0.0, x, y); + set_detector_pos(p_peak, 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); - select_scalable_reflections(image->reflections); - - /* Need these lists to work fast */ - optimise_reflist(image->reflections); + select_scalable_reflections(image->reflections, scalable); } diff --git a/src/geometry.h b/src/geometry.h index f8c007fa..171368a7 100644 --- a/src/geometry.h +++ b/src/geometry.h @@ -24,8 +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); +extern void update_partialities(struct image *image, const char *sym, + ReflItemList *scalable, + int *n_expected, int *n_found, int *n_notfound); + #endif /* GEOMETRY_H */ diff --git a/src/image.h b/src/image.h index 22e0ddae..d6e22aee 100644 --- a/src/image.h +++ b/src/image.h @@ -87,7 +87,6 @@ typedef struct _imagefeaturelist ImageFeatureList; * int height; * * RefList *reflections; - * RefList *raw_reflections; * * ImageFeatureList *features; * }; @@ -151,9 +150,6 @@ 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; diff --git a/src/partialator.c b/src/partialator.c index 9f931496..dc29d0d6 100644 --- a/src/partialator.c +++ b/src/partialator.c @@ -155,7 +155,7 @@ int main(int argc, char *argv[]) FILE *fh; int nthreads = 1; struct detector *det; - ReflItemList *obs; + ReflItemList *scalable; int i; int n_total_patterns; struct image *images; @@ -282,9 +282,11 @@ int main(int argc, char *argv[]) /* Fill in what we know about the images so far */ rewind(fh); - obs = new_items(); + scalable = new_items(); for ( i=0; iindexed_cell, &asx, &asy, &asz, &bsx, &bsy, &bsz, &csx, &csy, &csz); - xl = hi*asx + ki*bsx + li*csx; - yl = hi*asy + ki*bsy + li*csy; - zl = hi*asz + ki*bsz + li*csz; + xl = hs*asx + ks*bsx + ls*csx; + yl = hs*asy + ks*bsy + ls*csy; + zl = hs*asz + ks*bsz + ls*csz; - ds = 2.0 * resolution(image->indexed_cell, hi, ki, li); - tt = angle_between(0.0, 0.0, 1.0, xl, yl, zl+k); + ds = 2.0 * resolution(image->indexed_cell, hs, ks, ls); + tt = angle_between(0.0, 0.0, 1.0, xl, yl, zl+1.0/image->lambda); azi = angle_between(1.0, 0.0, 0.0, xl, yl, 0.0); get_partial(refl, &r1, &r2, &p, &clamp_low, &clamp_high); @@ -142,23 +142,23 @@ static double gradient(struct image *image, int k, Reflection *refl, double r) /* Cell parameters and orientation */ case REF_ASX : - return hi * sin(tt) * cos(azi) * g; + return hs * sin(tt) * cos(azi) * g; case REF_BSX : - return ki * sin(tt) * g; + return ks * sin(tt) * g; case REF_CSX : - return li * sin(tt) * g; + return ls * sin(tt) * g; case REF_ASY : - return hi * sin(tt) * g; + return hs * sin(tt) * g; case REF_BSY : - return ki * sin(tt) * g; + return ks * sin(tt) * g; case REF_CSY : - return li * sin(tt) * g; + return ls * sin(tt) * g; case REF_ASZ : - return hi * cos(tt) * g; + return hs * cos(tt) * g; case REF_BSZ : - return ki * cos(tt) * g; + return ks * cos(tt) * g; case REF_CSZ : - return li * cos(tt) * g; + return ls * cos(tt) * g; } @@ -285,8 +285,7 @@ static double pr_iterate(struct image *image, const RefList *full, /* Calculate all gradients for this reflection */ for ( k=0; kprofile_radius); + gr = gradient(image, k, refl, image->profile_radius); gradients[k] = gr; } @@ -403,8 +402,7 @@ static void plot_curve(struct image *image, const RefList *full, ax = origval + (double)i*shval; cell_set_reciprocal(cell, ax, ay, az, bx, by, bz, cx, cy, cz); - update_partialities_and_asymm(image, sym, - NULL, NULL, NULL, NULL); + update_partialities(image, sym, NULL, NULL, NULL, NULL); dev = mean_partial_dev(image, full, sym); STATUS("%i %e %e\n", i, ax, dev); @@ -431,8 +429,7 @@ void pr_refine(struct image *image, const RefList *full, const char *sym) double dev; max_shift = pr_iterate(image, full, sym); - update_partialities_and_asymm(image, sym, - NULL, NULL, NULL, NULL); + update_partialities(image, sym, NULL, NULL, NULL, NULL); dev = mean_partial_dev(image, full, sym); STATUS("PR Iteration %2i: max shift = %5.2f dev = %5.2f\n", diff --git a/src/reflist-utils.c b/src/reflist-utils.c index 0a5c4863..016fa317 100644 --- a/src/reflist-utils.c +++ b/src/reflist-utils.c @@ -362,7 +362,7 @@ RefList *read_reflections(const char *filename) } -RefList *asymmetric_indices(RefList *in, const char *sym, ReflItemList *obs) +RefList *asymmetric_indices(RefList *in, const char *sym) { Reflection *refl; RefListIterator *iter; @@ -386,12 +386,7 @@ RefList *asymmetric_indices(RefList *in, const char *sym, ReflItemList *obs) assert(cr != NULL); copy_data(cr, refl); - - if ( obs != NULL ) { - if ( !find_item(obs, ha, ka, la) ) { - add_item(obs, ha, ka, la); - } - } + set_symmetric_indices(cr, h, k, l); } diff --git a/src/reflist-utils.h b/src/reflist-utils.h index 2efa3a76..775e375a 100644 --- a/src/reflist-utils.h +++ b/src/reflist-utils.h @@ -41,7 +41,6 @@ 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); +extern RefList *asymmetric_indices(RefList *in, const char *sym); #endif /* REFLIST_UTILS_H */ diff --git a/src/reflist.c b/src/reflist.c index e420b736..d23b952c 100644 --- a/src/reflist.c +++ b/src/reflist.c @@ -45,6 +45,11 @@ struct _refldata { + /* Symmetric indices (i.e. the "real" indices) */ + signed int hs; + signed int ks; + signed int ls; + /* Partiality and related geometrical stuff */ double r1; /* First excitation error */ double r2; /* Second excitation error */ @@ -287,6 +292,29 @@ void get_indices(const Reflection *refl, } +/** + * get_symmetric_indices: + * @refl: A %Reflection + * @h: Location at which to store the 'h' index of the reflection + * @k: Location at which to store the 'k' index of the reflection + * @l: Location at which to store the 'l' index of the reflection + * + * This function gives the symmetric indices, that is, the "real" indices before + * squashing down to the asymmetric reciprocal unit. This may be useful if the + * list is indexed according to the asymmetric indices, but you still need + * access to the symmetric version. This happens during post-refinement. + * + **/ +void get_symmetric_indices(const Reflection *refl, + signed int *hs, signed int *ks, + signed int *ls) +{ + *hs = refl->data.hs; + *ks = refl->data.ks; + *ls = refl->data.ls; +} + + /** * get_partiality: * @refl: A %Reflection @@ -489,6 +517,15 @@ void set_ph(Reflection *refl, double phase) } +void set_symmetric_indices(Reflection *refl, + signed int hs, signed int ks, signed int ls) +{ + refl->data.hs = hs; + refl->data.ks = ks; + refl->data.ls = ls; +} + + /********************************* Insertion **********************************/ static void insert_node(Reflection *head, Reflection *new) diff --git a/src/reflist.h b/src/reflist.h index f380ba85..580ea33c 100644 --- a/src/reflist.h +++ b/src/reflist.h @@ -47,6 +47,9 @@ extern void get_detector_pos(const Reflection *refl, double *fs, double *ss); extern double get_partiality(const Reflection *refl); extern void get_indices(const Reflection *refl, signed int *h, signed int *k, signed int *l); +extern void get_symmetric_indices(const Reflection *refl, + signed int *hs, signed int *ks, + signed int *ls); extern double get_intensity(const Reflection *refl); extern void get_partial(const Reflection *refl, double *r1, double *r2, double *p, int *clamp_low, int *clamp_high); @@ -68,6 +71,8 @@ extern void set_redundancy(Reflection *refl, int red); extern void set_sum_squared_dev(Reflection *refl, double dev); extern void set_esd_intensity(Reflection *refl, double esd); extern void set_ph(Reflection *refl, double phase); +extern void set_symmetric_indices(Reflection *refl, + signed int hs, signed int ks, signed int ls); /* Insertion */ extern Reflection *add_refl(RefList *list, -- cgit v1.2.3