aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--src/geometry.c73
-rw-r--r--src/geometry.h5
-rw-r--r--src/hrs-scaling.c1
-rw-r--r--src/image.h4
-rw-r--r--src/partialator.c117
-rw-r--r--src/reflist-utils.c38
-rw-r--r--src/reflist-utils.h3
7 files changed, 156 insertions, 85 deletions
diff --git a/src/geometry.c b/src/geometry.c
index 749c194e..3ff39f31 100644
--- a/src/geometry.c
+++ b/src/geometry.c
@@ -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);
+}
diff --git a/src/geometry.h b/src/geometry.h
index a0761456..f8c007fa 100644
--- a/src/geometry.h
+++ b/src/geometry.h
@@ -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 */
diff --git a/src/hrs-scaling.c b/src/hrs-scaling.c
index 03458670..4b9b932d 100644
--- a/src/hrs-scaling.c
+++ b/src/hrs-scaling.c
@@ -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);
diff --git a/src/image.h b/src/image.h
index d6e22aee..22e0ddae 100644
--- a/src/image.h
+++ b/src/image.h
@@ -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;
diff --git a/src/partialator.c b/src/partialator.c
index 33dcad68..54022e9e 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -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);
}
diff --git a/src/reflist-utils.c b/src/reflist-utils.c
index ea590aaf..0a5c4863 100644
--- a/src/reflist-utils.c
+++ b/src/reflist-utils.c
@@ -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;
+}
diff --git a/src/reflist-utils.h b/src/reflist-utils.h
index b7ecd3ea..2efa3a76 100644
--- a/src/reflist-utils.h
+++ b/src/reflist-utils.h
@@ -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 */
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-15 02:38:18 +0000