aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--src/geometry.c80
-rw-r--r--src/geometry.h8
-rw-r--r--src/image.h4
-rw-r--r--src/partialator.c33
-rw-r--r--src/post-refinement.c41
-rw-r--r--src/reflist-utils.c9
-rw-r--r--src/reflist-utils.h3
-rw-r--r--src/reflist.c37
-rw-r--r--src/reflist.h5
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; i<n_total_patterns; i++ ) {
+ RefList *as;
+
images[n_usable_patterns].det = NULL;
if ( read_chunk(fh, &images[n_usable_patterns]) != 0 ) {
@@ -316,13 +318,14 @@ int main(int argc, char *argv[])
images[n_usable_patterns].beam = NULL;
/* This is the raw list of reflections */
- images[n_usable_patterns].raw_reflections =
- images[n_usable_patterns].reflections;
- images[n_usable_patterns].reflections = NULL;
+ as = asymmetric_indices(images[n_usable_patterns].reflections,
+ sym);
+ optimise_reflist(as);
+ reflist_free(images[n_usable_patterns].reflections);
+ images[n_usable_patterns].reflections = as;
- update_partialities_and_asymm(&images[n_usable_patterns], sym,
- obs, &n_expected, &n_found,
- &n_notfound);
+ update_partialities(&images[n_usable_patterns], sym, scalable,
+ &n_expected, &n_found, &n_notfound);
progress_bar(i, n_total_patterns-1, "Loading pattern data");
n_usable_patterns++;
@@ -331,14 +334,15 @@ int main(int argc, char *argv[])
fclose(fh);
STATUS("Found %5.2f%% of the expected peaks (missed %i of %i).\n",
100.0 * (double)n_found / n_expected, n_notfound, n_expected);
- STATUS("Mean measurements per unique reflection: %5.2f\n",
- (double)n_found / num_items(obs));
+ STATUS("Mean measurements per scalable unique reflection: %5.2f\n",
+ (double)n_found / num_items(scalable));
cref = find_common_reflections(images, n_usable_patterns);
/* Make initial estimates */
STATUS("Performing initial scaling.\n");
- full = scale_intensities(images, n_usable_patterns, sym, obs, cref);
+ full = scale_intensities(images, n_usable_patterns, sym,
+ scalable, cref);
/* Iterate */
for ( i=0; i<n_iter; i++ ) {
@@ -364,13 +368,13 @@ 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, full,
+ refine_all(images, n_total_patterns, det, sym, scalable, full,
nthreads, fhg, fhp);
/* Re-estimate all the full intensities */
reflist_free(full);
full = scale_intensities(images, n_usable_patterns,
- sym, obs, cref);
+ sym, scalable, cref);
fclose(fhg);
fclose(fhp);
@@ -388,10 +392,9 @@ int main(int argc, char *argv[])
/* Clean up */
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);
+ delete_items(scalable);
free(sym);
free(outfile);
free_detector_geometry(det);
diff --git a/src/post-refinement.c b/src/post-refinement.c
index d804dbc2..1fd4d13b 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -95,21 +95,21 @@ static double gradient(struct image *image, int k, Reflection *refl, double r)
double bsx, bsy, bsz;
double csx, csy, csz;
double xl, yl, zl;
- signed int hi, ki, li;
+ signed int hs, ks, ls;
double r1, r2, p;
int clamp_low, clamp_high;
- get_indices(refl, &hi, &ki, &li);
+ get_symmetric_indices(refl, &hs, &ks, &ls);
cell_get_reciprocal(image->indexed_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; k<NUM_PARAMS; k++ ) {
double gr;
- gr = gradient(image, k, refl,
- image->profile_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 */
@@ -288,6 +293,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,
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-08 18:19:48 +0000