aboutsummaryrefslogtreecommitdiff
path: root/src/process_hkl.c
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2010-07-14 17:58:55 +0200
committerThomas White <taw@physics.org>2012-02-22 15:26:53 +0100
commit6a5422356c15962726df2261aa53354b0ff12662 (patch)
treeb6c5ae80b837bda38957cef07816b511515ffdcb /src/process_hkl.c
parent6a476e010468f27e02df6bb90a1ea197bd9d039d (diff)
Reduce the scope of "count"
Lists of counts had pervaded every corner of CrystFEL, being used as markers for the presence of reflections. Now we have a better way of doing this, the ReflItemList, and few parts of the suite apart from process_hkl have any business knowing how many observations were made of a particular reflection.
Diffstat (limited to 'src/process_hkl.c')
-rw-r--r--src/process_hkl.c80
1 files changed, 48 insertions, 32 deletions
diff --git a/src/process_hkl.c b/src/process_hkl.c
index b5b6df9a..986477fb 100644
--- a/src/process_hkl.c
+++ b/src/process_hkl.c
@@ -89,6 +89,8 @@ static ReflItemList *get_twin_possibilities(const char *holo, const char *mero)
twins = get_twins(test_items, holo, mero);
delete_items(test_items);
+ /* Idiot check. Wouldn't be necessary if I could prove that the above
+ * set of arbitrarily chosen reflections were always general. */
if ( num_items(twins) != np ) {
ERROR("Whoops! Couldn't find all the twinning possiblities.\n");
abort();
@@ -98,10 +100,9 @@ static ReflItemList *get_twin_possibilities(const char *holo, const char *mero)
}
-static int resolve_twin(ReflItemList *items, ReflItemList *twins,
- const double *patt, const double *model,
- const unsigned int *model_counts,
- const char *holo, const char *mero)
+static int resolve_twin(const double *model, ReflItemList *observed,
+ const double *patt, ReflItemList *items,
+ ReflItemList *twins, const char *holo, const char *mero)
{
int n, i;
double best_fom = 0.0;
@@ -116,6 +117,7 @@ static int resolve_twin(ReflItemList *items, ReflItemList *twins,
double *trial_ints = new_list_intensity();
unsigned int *trial_counts = new_list_count();
double fom;
+ ReflItemList *intersection;
op = get_item(twins, i)->op;
@@ -134,8 +136,9 @@ static int resolve_twin(ReflItemList *items, ReflItemList *twins,
}
- fom = stat_pearson(trial_ints, trial_counts,
- model, model_counts);
+ intersection = intersection_items(observed, items);
+ fom = stat_pearson(trial_ints, model, intersection);
+ delete_items(intersection);
free(trial_ints);
free(trial_counts);
@@ -153,9 +156,9 @@ static int resolve_twin(ReflItemList *items, ReflItemList *twins,
}
-static void merge_pattern(double *model, const double *new,
- unsigned int *model_counts,
- ReflItemList *items, int mo, int sum,
+static void merge_pattern(double *model, ReflItemList *observed,
+ const double *new, ReflItemList *items,
+ unsigned int *model_counts, int mo,
ReflItemList *twins,
const char *holo, const char *mero)
{
@@ -163,8 +166,8 @@ static void merge_pattern(double *model, const double *new,
int twin;
if ( twins != NULL ) {
- twin = resolve_twin(items, twins, new, model,
- model_counts, holo, mero);
+ twin = resolve_twin(model, observed, new, items,
+ twins, holo, mero);
} else {
twin = 0;
}
@@ -187,25 +190,23 @@ static void merge_pattern(double *model, const double *new,
intensity = lookup_intensity(new, h, k, l);
+ /* User asked for max only? */
if ( !mo ) {
integrate_intensity(model, h, k, l, intensity);
- if ( sum ) {
- set_count(model_counts, h, k, l, 1);
- } else {
- integrate_count(model_counts, h, k, l, 1);
- }
} else {
if ( intensity > lookup_intensity(model, h, k, l) ) {
set_intensity(model, h, k, l, intensity);
}
- set_count(model_counts, h, k, l, 1);
}
+ integrate_count(model_counts, h, k, l, 1);
+
}
}
-static void merge_all(FILE *fh, double *model, unsigned int *model_counts,
+static void merge_all(FILE *fh, double **pmodel, ReflItemList **pobserved,
+ unsigned int **pcounts,
int config_maxonly, int config_scale, int config_sum,
int config_stopafter,
ReflItemList *twins, const char *holo, const char *mero,
@@ -218,6 +219,10 @@ static void merge_all(FILE *fh, double *model, unsigned int *model_counts,
int n_patterns = 0;
double *new_pattern = new_list_intensity();
ReflItemList *items = new_items();
+ ReflItemList *observed = new_items();
+ double *model = new_list_intensity();
+ unsigned int *counts = new_list_count();
+ int i;
do {
@@ -244,19 +249,20 @@ static void merge_all(FILE *fh, double *model, unsigned int *model_counts,
/* Scale if requested */
if ( config_scale ) {
- scale_intensities(model, new_pattern,
- model_counts, items, f0,
- f0_valid);
+ scale_intensities(model, observed,
+ new_pattern, items,
+ f0, f0_valid);
}
/* Start of second or later pattern */
- merge_pattern(model, new_pattern, model_counts,
- items, config_maxonly, config_sum, twins,
- holo, mero);
+ merge_pattern(model, observed, new_pattern, items,
+ counts, config_maxonly,
+ twins, holo, mero);
if ( n_patterns == config_stopafter ) break;
n_patterns++;
+ union_items(observed, items);
clear_items(items);
progress_bar(n_patterns, n_total_patterns, "Merging");
@@ -292,6 +298,19 @@ static void merge_all(FILE *fh, double *model, unsigned int *model_counts,
delete_items(items);
free(new_pattern);
+ if ( !config_sum ) {
+ for ( i=0; i<IDIM*IDIM*IDIM; i++ ) {
+ if ( counts[i] > 0 ) {
+ model[i] /= (double)counts[i];
+ counts[i] = 1;
+ }
+ }
+ }
+
+ *pmodel = model;
+ *pcounts = counts;
+ *pobserved = observed;
+
STATUS("%i patterns had no f0 valid value.\n", n_nof0);
}
@@ -322,7 +341,7 @@ int main(int argc, char *argv[])
char *output = NULL;
FILE *fh;
double *model;
- unsigned int *model_counts;
+ unsigned int *counts;
UnitCell *cell;
int config_maxonly = 0;
int config_stopafter = 0;
@@ -332,6 +351,7 @@ int main(int argc, char *argv[])
char *sym = NULL;
char *pdb = NULL;
ReflItemList *twins;
+ ReflItemList *observed;
int i;
/* Long options */
@@ -400,9 +420,6 @@ int main(int argc, char *argv[])
sym = strdup("1");
}
- model = new_list_intensity();
- model_counts = new_list_count();
-
cell = load_cell_from_pdb(pdb);
free(pdb);
@@ -443,7 +460,7 @@ int main(int argc, char *argv[])
STATUS("There are %i patterns to process\n", n_total_patterns);
rewind(fh);
- merge_all(fh, model, model_counts,
+ merge_all(fh, &model, &observed, &counts,
config_maxonly, config_scale, config_sum, config_stopafter,
twins, holo, sym, n_total_patterns);
rewind(fh);
@@ -451,13 +468,12 @@ int main(int argc, char *argv[])
fclose(fh);
if ( output != NULL ) {
- write_reflections(output, model_counts, model, NULL,
- 0, cell, 1);
+ write_reflections(output, observed, model, NULL, counts, cell);
}
free(sym);
free(model);
- free(model_counts);
+ free(counts);
free(output);
free(cell);