aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--Makefile.am14
-rw-r--r--Makefile.in43
-rw-r--r--src/process_hkl.c570
-rw-r--r--src/reflist-utils.c72
-rw-r--r--src/reflist-utils.h30
-rw-r--r--src/reflist.c43
-rw-r--r--src/reflist.h6
-rw-r--r--src/stream.c128
-rw-r--r--src/stream.h4
9 files changed, 389 insertions, 521 deletions
diff --git a/Makefile.am b/Makefile.am
index ce2d06d7..f7437acc 100644
--- a/Makefile.am
+++ b/Makefile.am
@@ -42,10 +42,10 @@ tests_gpu_sim_check_SOURCES = tests/gpu_sim_check.c src/utils.c \
endif
src_process_hkl_SOURCES = src/process_hkl.c src/sfac.c src/statistics.c \
- src/cell.c src/utils.c src/reflections.c \
+ src/cell.c src/utils.c \
src/symmetry.c src/stream.c src/beam-parameters.c \
src/thread-pool.c src/image.c src/detector.c \
- src/hdf5-file.c src/reflist.c
+ src/hdf5-file.c src/reflist.c src/reflist-utils.c
src_indexamajig_SOURCES = src/indexamajig.c src/hdf5-file.c src/utils.c \
src/cell.c src/image.c src/peaks.c src/index.c \
@@ -53,7 +53,8 @@ src_indexamajig_SOURCES = src/indexamajig.c src/hdf5-file.c src/utils.c \
src/sfac.c src/dirax.c src/mosflm.c \
src/reflections.c src/symmetry.c \
src/geometry.c src/thread-pool.c \
- src/beam-parameters.c src/reflist.c src/stream.c
+ src/beam-parameters.c src/reflist.c src/stream.c \
+ src/reflist-utils.c
if HAVE_OPENCL
src_indexamajig_SOURCES += src/diffraction-gpu.c src/cl-utils.c
endif
@@ -99,18 +100,19 @@ src_partialator_SOURCES = src/partialator.c src/cell.c src/hdf5-file.c \
src/geometry.c src/reflections.c src/stream.c \
src/thread-pool.c src/beam-parameters.c \
src/symmetry.c src/post-refinement.c \
- src/hrs-scaling.c src/reflist.c
+ src/hrs-scaling.c src/reflist.c src/reflist-utils.c
if BUILD_CUBEIT
src_cubeit_SOURCES = src/cubeit.c src/cell.c src/hdf5-file.c src/utils.c \
src/detector.c src/render.c src/filters.c src/image.c \
- src/symmetry.c src/stream.c src/thread-pool.c src/reflist.c
+ src/symmetry.c src/stream.c src/thread-pool.c \
+ src/reflist.c src/reflist-utils.c
endif
src_estimate_background_SOURCES = src/estimate_background.c src/stream.c \
src/utils.c src/cell.c src/thread-pool.c \
src/image.c src/detector.c src/hdf5-file.c \
- src/reflist.c
+ src/reflist.c src/reflist-utils.c
tests_list_check_SOURCES = tests/list_check.c src/reflist.c src/thread-pool.c \
src/utils.c
diff --git a/Makefile.in b/Makefile.in
index b65b60a3..402094df 100644
--- a/Makefile.in
+++ b/Makefile.in
@@ -107,7 +107,7 @@ src_compare_hkl_DEPENDENCIES = $(top_builddir)/lib/libgnu.a
am__src_cubeit_SOURCES_DIST = src/cubeit.c src/cell.c src/hdf5-file.c \
src/utils.c src/detector.c src/render.c src/filters.c \
src/image.c src/symmetry.c src/stream.c src/thread-pool.c \
- src/reflist.c
+ src/reflist.c src/reflist-utils.c
@BUILD_CUBEIT_TRUE@am_src_cubeit_OBJECTS = src/cubeit.$(OBJEXT) \
@BUILD_CUBEIT_TRUE@ src/cell.$(OBJEXT) src/hdf5-file.$(OBJEXT) \
@BUILD_CUBEIT_TRUE@ src/utils.$(OBJEXT) src/detector.$(OBJEXT) \
@@ -115,7 +115,8 @@ am__src_cubeit_SOURCES_DIST = src/cubeit.c src/cell.c src/hdf5-file.c \
@BUILD_CUBEIT_TRUE@ src/image.$(OBJEXT) src/symmetry.$(OBJEXT) \
@BUILD_CUBEIT_TRUE@ src/stream.$(OBJEXT) \
@BUILD_CUBEIT_TRUE@ src/thread-pool.$(OBJEXT) \
-@BUILD_CUBEIT_TRUE@ src/reflist.$(OBJEXT)
+@BUILD_CUBEIT_TRUE@ src/reflist.$(OBJEXT) \
+@BUILD_CUBEIT_TRUE@ src/reflist-utils.$(OBJEXT)
src_cubeit_OBJECTS = $(am_src_cubeit_OBJECTS)
src_cubeit_LDADD = $(LDADD)
src_cubeit_DEPENDENCIES = $(top_builddir)/lib/libgnu.a
@@ -124,7 +125,7 @@ am_src_estimate_background_OBJECTS = \
src/utils.$(OBJEXT) src/cell.$(OBJEXT) \
src/thread-pool.$(OBJEXT) src/image.$(OBJEXT) \
src/detector.$(OBJEXT) src/hdf5-file.$(OBJEXT) \
- src/reflist.$(OBJEXT)
+ src/reflist.$(OBJEXT) src/reflist-utils.$(OBJEXT)
src_estimate_background_OBJECTS = \
$(am_src_estimate_background_OBJECTS)
src_estimate_background_LDADD = $(LDADD)
@@ -154,8 +155,8 @@ am__src_indexamajig_SOURCES_DIST = src/indexamajig.c src/hdf5-file.c \
src/filters.c src/diffraction.c src/detector.c src/sfac.c \
src/dirax.c src/mosflm.c src/reflections.c src/symmetry.c \
src/geometry.c src/thread-pool.c src/beam-parameters.c \
- src/reflist.c src/stream.c src/diffraction-gpu.c \
- src/cl-utils.c
+ src/reflist.c src/stream.c src/reflist-utils.c \
+ src/diffraction-gpu.c src/cl-utils.c
@HAVE_OPENCL_TRUE@am__objects_1 = src/diffraction-gpu.$(OBJEXT) \
@HAVE_OPENCL_TRUE@ src/cl-utils.$(OBJEXT)
am_src_indexamajig_OBJECTS = src/indexamajig.$(OBJEXT) \
@@ -166,7 +167,8 @@ am_src_indexamajig_OBJECTS = src/indexamajig.$(OBJEXT) \
src/mosflm.$(OBJEXT) src/reflections.$(OBJEXT) \
src/symmetry.$(OBJEXT) src/geometry.$(OBJEXT) \
src/thread-pool.$(OBJEXT) src/beam-parameters.$(OBJEXT) \
- src/reflist.$(OBJEXT) src/stream.$(OBJEXT) $(am__objects_1)
+ src/reflist.$(OBJEXT) src/stream.$(OBJEXT) \
+ src/reflist-utils.$(OBJEXT) $(am__objects_1)
src_indexamajig_OBJECTS = $(am_src_indexamajig_OBJECTS)
src_indexamajig_LDADD = $(LDADD)
src_indexamajig_DEPENDENCIES = $(top_builddir)/lib/libgnu.a
@@ -177,7 +179,7 @@ am_src_partialator_OBJECTS = src/partialator.$(OBJEXT) \
src/stream.$(OBJEXT) src/thread-pool.$(OBJEXT) \
src/beam-parameters.$(OBJEXT) src/symmetry.$(OBJEXT) \
src/post-refinement.$(OBJEXT) src/hrs-scaling.$(OBJEXT) \
- src/reflist.$(OBJEXT)
+ src/reflist.$(OBJEXT) src/reflist-utils.$(OBJEXT)
src_partialator_OBJECTS = $(am_src_partialator_OBJECTS)
src_partialator_LDADD = $(LDADD)
src_partialator_DEPENDENCIES = $(top_builddir)/lib/libgnu.a
@@ -205,11 +207,11 @@ src_powder_plot_LDADD = $(LDADD)
src_powder_plot_DEPENDENCIES = $(top_builddir)/lib/libgnu.a
am_src_process_hkl_OBJECTS = src/process_hkl.$(OBJEXT) \
src/sfac.$(OBJEXT) src/statistics.$(OBJEXT) src/cell.$(OBJEXT) \
- src/utils.$(OBJEXT) src/reflections.$(OBJEXT) \
- src/symmetry.$(OBJEXT) src/stream.$(OBJEXT) \
- src/beam-parameters.$(OBJEXT) src/thread-pool.$(OBJEXT) \
- src/image.$(OBJEXT) src/detector.$(OBJEXT) \
- src/hdf5-file.$(OBJEXT) src/reflist.$(OBJEXT)
+ src/utils.$(OBJEXT) src/symmetry.$(OBJEXT) \
+ src/stream.$(OBJEXT) src/beam-parameters.$(OBJEXT) \
+ src/thread-pool.$(OBJEXT) src/image.$(OBJEXT) \
+ src/detector.$(OBJEXT) src/hdf5-file.$(OBJEXT) \
+ src/reflist.$(OBJEXT) src/reflist-utils.$(OBJEXT)
src_process_hkl_OBJECTS = $(am_src_process_hkl_OBJECTS)
src_process_hkl_LDADD = $(LDADD)
src_process_hkl_DEPENDENCIES = $(top_builddir)/lib/libgnu.a
@@ -612,17 +614,17 @@ src_pattern_sim_SOURCES = src/pattern_sim.c src/diffraction.c \
@HAVE_OPENCL_TRUE@ src/image.c
src_process_hkl_SOURCES = src/process_hkl.c src/sfac.c src/statistics.c \
- src/cell.c src/utils.c src/reflections.c \
+ src/cell.c src/utils.c \
src/symmetry.c src/stream.c src/beam-parameters.c \
src/thread-pool.c src/image.c src/detector.c \
- src/hdf5-file.c src/reflist.c
+ src/hdf5-file.c src/reflist.c src/reflist-utils.c
src_indexamajig_SOURCES = src/indexamajig.c src/hdf5-file.c \
src/utils.c src/cell.c src/image.c src/peaks.c src/index.c \
src/filters.c src/diffraction.c src/detector.c src/sfac.c \
src/dirax.c src/mosflm.c src/reflections.c src/symmetry.c \
src/geometry.c src/thread-pool.c src/beam-parameters.c \
- src/reflist.c src/stream.c $(am__append_6)
+ src/reflist.c src/stream.c src/reflist-utils.c $(am__append_6)
@BUILD_HDFSEE_TRUE@src_hdfsee_SOURCES = src/hdfsee.c src/dw-hdfsee.c src/render.c \
@BUILD_HDFSEE_TRUE@ src/hdf5-file.c src/utils.c src/image.c src/filters.c \
@BUILD_HDFSEE_TRUE@ src/thread-pool.c src/detector.c
@@ -661,16 +663,17 @@ src_partialator_SOURCES = src/partialator.c src/cell.c src/hdf5-file.c \
src/geometry.c src/reflections.c src/stream.c \
src/thread-pool.c src/beam-parameters.c \
src/symmetry.c src/post-refinement.c \
- src/hrs-scaling.c src/reflist.c
+ src/hrs-scaling.c src/reflist.c src/reflist-utils.c
@BUILD_CUBEIT_TRUE@src_cubeit_SOURCES = src/cubeit.c src/cell.c src/hdf5-file.c src/utils.c \
@BUILD_CUBEIT_TRUE@ src/detector.c src/render.c src/filters.c src/image.c \
-@BUILD_CUBEIT_TRUE@ src/symmetry.c src/stream.c src/thread-pool.c src/reflist.c
+@BUILD_CUBEIT_TRUE@ src/symmetry.c src/stream.c src/thread-pool.c \
+@BUILD_CUBEIT_TRUE@ src/reflist.c src/reflist-utils.c
src_estimate_background_SOURCES = src/estimate_background.c src/stream.c \
src/utils.c src/cell.c src/thread-pool.c \
src/image.c src/detector.c src/hdf5-file.c \
- src/reflist.c
+ src/reflist.c src/reflist-utils.c
tests_list_check_SOURCES = tests/list_check.c src/reflist.c src/thread-pool.c \
src/utils.c
@@ -851,6 +854,8 @@ src/stream.$(OBJEXT): src/$(am__dirstamp) \
src/$(DEPDIR)/$(am__dirstamp)
src/reflist.$(OBJEXT): src/$(am__dirstamp) \
src/$(DEPDIR)/$(am__dirstamp)
+src/reflist-utils.$(OBJEXT): src/$(am__dirstamp) \
+ src/$(DEPDIR)/$(am__dirstamp)
src/cubeit$(EXEEXT): $(src_cubeit_OBJECTS) $(src_cubeit_DEPENDENCIES) src/$(am__dirstamp)
@rm -f src/cubeit$(EXEEXT)
$(AM_V_CCLD)$(LINK) $(src_cubeit_OBJECTS) $(src_cubeit_LDADD) $(LIBS)
@@ -977,6 +982,7 @@ mostlyclean-compile:
-rm -f src/powder_plot.$(OBJEXT)
-rm -f src/process_hkl.$(OBJEXT)
-rm -f src/reflections.$(OBJEXT)
+ -rm -f src/reflist-utils.$(OBJEXT)
-rm -f src/reflist.$(OBJEXT)
-rm -f src/render.$(OBJEXT)
-rm -f src/render_hkl.$(OBJEXT)
@@ -1024,6 +1030,7 @@ distclean-compile:
@AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/powder_plot.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/process_hkl.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/reflections.Po@am__quote@
+@AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/reflist-utils.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/reflist.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/render.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/render_hkl.Po@am__quote@
diff --git a/src/process_hkl.c b/src/process_hkl.c
index 8df2dae1..da3ee3f8 100644
--- a/src/process_hkl.c
+++ b/src/process_hkl.c
@@ -3,7 +3,7 @@
*
* Assemble and process FEL Bragg intensities
*
- * (c) 2006-2010 Thomas White <taw@physics.org>
+ * (c) 2006-2011 Thomas White <taw@physics.org>
*
* Part of CrystFEL - crystallography with a FEL
*
@@ -24,10 +24,12 @@
#include "utils.h"
#include "statistics.h"
#include "sfac.h"
-#include "reflections.h"
+#include "reflist-utils.h"
#include "symmetry.h"
#include "stream.h"
#include "beam-parameters.h"
+#include "reflist.h"
+#include "image.h"
/* Number of divisions for intensity histograms */
@@ -43,7 +45,7 @@ static void show_help(const char *s)
" -h, --help Display this help message.\n"
" -i, --input=<filename> Specify input filename (\"-\" for stdin).\n"
" -o, --output=<filename> Specify output filename for merged intensities\n"
-" (don't specify for no output).\n"
+" Default: processed.hkl).\n"
" -p, --pdb=<filename> PDB file to use (default: molecule.pdb).\n"
" -b, --beam=<file> Get beam parameters from file (used for sigmas).\n"
"\n"
@@ -66,13 +68,10 @@ static void show_help(const char *s)
" --scale Scale each pattern for best fit with the current\n"
" model.\n"
" -y, --symmetry=<sym> Merge according to point group <sym>.\n"
-" -a, --input-symmetry=<a> Specify the apparent (input) symmetry.\n"
" --reference=<file> Compare against intensities from <file> when\n"
" scaling or resolving ambiguities.\n"
" The symmetry of the reference list must be the\n"
" same as that given with '-y'.\n"
-" --outstream=<file> Write an annotated version of the input stream\n"
-" to <file>.\n"
);
}
@@ -119,175 +118,55 @@ static void plot_histogram(double *vals, int n)
}
-/* Note "holo" needn't actually be a holohedral point group, if you want to try
- * something strange like resolving from a low-symmetry group into an even
- * lower symmetry one.
- */
-static ReflItemList *get_twin_possibilities(const char *holo, const char *mero)
-{
- ReflItemList *test_items;
- ReflItemList *twins;
- int np;
-
- np = num_general_equivs(holo) / num_general_equivs(mero);
-
- test_items = new_items();
-
- /* Some arbitrarily chosen reflections which can't be special
- * reflections in any point group, i.e. lots of odd numbers,
- * prime numbers and so on. There's probably an analytical
- * way of working these out, but this will do. */
- add_item(test_items, 1, 2, 3);
- add_item(test_items, 3, 7, 13);
- add_item(test_items, 5, 2, 1);
-
- 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();
- }
-
- return twins;
-}
-
-
-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;
- int best_op = 0;
-
- n = num_items(twins);
-
- for ( i=0; i<n; i++ ) {
-
- int j;
- int op;
- double *trial_ints = new_list_intensity();
- unsigned int *trial_counts = new_list_count();
- double fom;
- ReflItemList *intersection;
-
- op = get_item(twins, i)->op;
-
- for ( j=0; j<num_items(items); j++ ) {
-
- signed int h, k, l;
- struct refl_item *r = get_item(items, j);
-
- get_general_equiv(r->h, r->k, r->l, &h, &k, &l,
- holo, op);
- get_asymm(h, k, l, &h, &k, &l, mero);
-
- set_intensity(trial_ints, h, k, l,
- lookup_intensity(patt, r->h, r->k, r->l));
- set_count(trial_counts, h, k, l, 1);
-
- }
-
- intersection = intersection_items(observed, items);
- fom = stat_pearson_i(trial_ints, model, intersection);
- delete_items(intersection);
-
- free(trial_ints);
- free(trial_counts);
-
- //printf(" %f", fom);
- if ( fom > best_fom ) {
- best_fom = fom;
- best_op = op;
- }
-
- }
- //printf("\n");
-
- return best_op;
-}
-
-
-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,
- ReflItemList *reference, const double *reference_i,
- double *hist_vals,
- signed int hist_h, signed int hist_k,
- signed int hist_l, int *hist_n, double *devs,
- double *means, FILE *outfh)
+static void merge_pattern(RefList *model, RefList *new, int max_only,
+ const char *sym,
+ double *hist_vals, signed int hist_h,
+ signed int hist_k, signed int hist_l, int *hist_n)
{
- int i;
- int twin;
- ReflItemList *sym_items = new_items();
-
- if ( twins != NULL ) {
- if ( reference != NULL ) {
- twin = resolve_twin(reference_i, reference, new, items,
- twins, holo, mero);
- } else {
- twin = resolve_twin(model, observed, new, items,
- twins, holo, mero);
- }
- } else {
- twin = 0;
- }
-
- if ( outfh != NULL ) {
- fprintf(outfh, "twin=%i\n", twin);
- }
+ Reflection *refl;
+ RefListIterator *iter;
- for ( i=0; i<num_items(items); i++ ) {
+ for ( refl = first_refl(new, &iter);
+ refl != NULL;
+ refl = next_refl(refl, iter) ) {
double intensity;
- signed int hs, ks, ls;
signed int h, k, l;
- struct refl_item *item;
-
- item = get_item(items, i);
+ Reflection *model_version;
+ double model_int;
- hs = item->h;
- ks = item->k;
- ls = item->l;
-
- /* Transform into correct side of the twin law.
- * "twin" is always zero if no de-twinning is performed. */
- get_general_equiv(hs, ks, ls, &h, &k, &l, holo, twin);
+ get_indices(refl, &h, &k, &l);
+ model_version = find_refl(model, h, k, l);
+ if ( model_version == NULL ) {
+ model_version = add_refl(model, h, k, l);
+ }
/* Put into the asymmetric unit for the target group */
- get_asymm(h, k, l, &h, &k, &l, mero);
+ get_asymm(h, k, l, &h, &k, &l, sym);
/* Read the intensity from the original location
* (i.e. before screwing around with symmetry) */
- intensity = lookup_intensity(new, hs, ks, ls);
+ intensity = get_intensity(refl);
+
+ /* Get the current model intensity */
+ model_int = get_intensity(model_version);
/* User asked for max only? */
- if ( !mo ) {
- integrate_intensity(model, h, k, l, intensity);
+ if ( !max_only ) {
+ set_int(model_version, model_int + intensity);
} else {
- if ( intensity > lookup_intensity(model, h, k, l) ) {
- set_intensity(model, h, k, l, intensity);
+ if ( intensity > get_intensity(model_version) ) {
+ set_int(model_version, intensity);
}
}
- if ( devs != NULL ) {
- double m;
- m = lookup_intensity(means, h, k, l);
- integrate_intensity(devs, h, k, l,
- pow(intensity-m, 2.0));
- }
-
- if ( !find_item(sym_items, h, k, l) ) {
- add_item(sym_items, h, k, l);
- }
+ double dev = get_sum_squared_dev(model_version);
+ set_sum_squared_dev(model_version,
+ dev + pow(intensity-model_int, 2.0));
- /* Increase count */
- integrate_count(model_counts, h, k, l, 1);
+ /* Increase redundancy */
+ int cur_redundancy = get_redundancy(model_version);
+ set_redundancy(model_version, cur_redundancy+1);
if ( hist_vals != NULL ) {
int p = *hist_n;
@@ -298,11 +177,6 @@ static void merge_pattern(double *model, ReflItemList *observed,
}
}
-
- /* Dump the reflections in this pattern into the overall list */
- union_items(observed, sym_items);
-
- delete_items(sym_items);
}
@@ -314,34 +188,36 @@ enum {
};
-static void scale_intensities(const double *model, ReflItemList *model_items,
- double *new_pattern, ReflItemList *new_items,
- double f0, int f0_valid, const char *sym)
+static void scale_intensities(RefList *model, RefList *new, const char *sym)
{
double s;
double top = 0.0;
double bot = 0.0;
- unsigned int i;
const int scaling = SCALE_INTPERBRAGG;
+ Reflection *refl;
+ RefListIterator *iter;
+ Reflection *model_version;
- for ( i=0; i<num_items(new_items); i++ ) {
+ for ( refl = first_refl(new, &iter);
+ refl != NULL;
+ refl = next_refl(refl, iter) ) {
double i1, i2;
- struct refl_item *it;
signed int hu, ku, lu;
+ signed int h, k, l;
- /* Get the next item in the list of new reflections */
- it = get_item(new_items, i);
+ get_indices(refl, &h, &k, &l);
switch ( scaling ) {
case SCALE_TWOPASS :
- /* Find the (only) partner in the model */
- find_unique_equiv(model_items, it->h, it->k, it->l, sym,
- &hu, &ku, &lu);
+ model_version = find_refl(model, h, k, l);
+ if ( model_version == NULL ) continue;
+
+ get_asymm(h, k, l, &hu, &ku, &lu, sym);
- i1 = lookup_intensity(model, hu, ku, lu);
- i2 = lookup_intensity(new_pattern, it->h, it->k, it->l);
+ i1 = get_intensity(model_version);
+ i2 = get_intensity(refl);
/* Calculate LSQ estimate of scaling factor */
top += i1 * i2;
@@ -352,7 +228,7 @@ static void scale_intensities(const double *model, ReflItemList *model_items,
case SCALE_CONSTINT :
/* Sum up the intensity in the pattern */
- i2 = lookup_intensity(new_pattern, it->h, it->k, it->l);
+ i2 = get_intensity(refl);
top += i2;
break;
@@ -360,7 +236,7 @@ static void scale_intensities(const double *model, ReflItemList *model_items,
case SCALE_INTPERBRAGG :
/* Sum up the intensity in the pattern */
- i2 = lookup_intensity(new_pattern, it->h, it->k, it->l);
+ i2 = get_intensity(refl);
top += i2;
bot += 1.0;
@@ -381,162 +257,108 @@ static void scale_intensities(const double *model, ReflItemList *model_items,
s = 1000.0 / (top/bot);
break;
}
- //if ( f0_valid ) printf("%f %f\n", s, f0);
/* Multiply the new pattern up by "s" */
- for ( i=0; i<LIST_SIZE; i++ ) {
- new_pattern[i] *= s;
+ for ( refl = first_refl(new, &iter);
+ refl != NULL;
+ refl = next_refl(refl, iter) ) {
+
+ double intensity = get_intensity(refl);
+ set_int(refl, intensity*s);
+
}
}
-static void merge_all(FILE *fh, double **pmodel, ReflItemList **pobserved,
- unsigned int **pcounts,
+static void merge_all(FILE *fh, RefList *model,
int config_maxonly, int config_scale, int config_sum,
int config_startafter, int config_stopafter,
- ReflItemList *twins, const char *holo, const char *mero,
+ const char *sym,
int n_total_patterns,
- ReflItemList *reference, double *reference_i,
- double *hist_vals,
- signed int hist_h, signed int hist_k, signed int hist_l,
- int *hist_i, double *devs, double *means, FILE *outfh)
+ double *hist_vals, signed int hist_h,
+ signed int hist_k, signed int hist_l,
+ int *hist_i)
{
- char *rval;
- float f0;
- int n_nof0 = 0;
- int f0_valid = 0;
+ int rval;
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;
-
- if ( config_startafter != 0 ) {
-
- do {
-
- char line[1024];
-
- rval = fgets(line, 1023, fh);
- if ( (strncmp(line, "Reflections from indexing", 25)
- == 0)
- || (strncmp(line, "New pattern", 11) == 0) ) {
- n_patterns++;
- }
-
- if ( n_patterns == config_startafter ) break;
-
- } while ( rval != NULL );
+ int n_used = 0;
+ Reflection *refl;
+ RefListIterator *iter;
+ if ( skip_some_files(fh, config_startafter) ) {
+ ERROR("Failed to skip first %i files.\n", config_startafter);
+ return;
}
do {
- char line[1024];
- signed int h, k, l;
- float intensity;
- int r;
+ struct image image;
- rval = fgets(line, 1023, fh);
- if ( ((strncmp(line, "Reflections from indexing", 25) == 0)
- || (rval == NULL)) && (num_items(items)>0) ) {
+ /* Get data from next chunk */
+ rval = read_chunk(fh, &image);
+ if ( rval ) continue;
- /* Assume a default I0 if we don't have one by now */
- if ( config_scale && !f0_valid ) {
- n_nof0++;
- f0 = 1.0;
- }
+ n_patterns++;
+
+ if ( (image.reflections != NULL) && (image.indexed_cell) ) {
/* Scale if requested */
if ( config_scale ) {
- if ( reference == NULL ) {
- scale_intensities(model, observed,
- new_pattern, items,
- f0, f0_valid, mero);
- } else {
- scale_intensities(reference_i,
- reference,
- new_pattern, items,
- f0, f0_valid, mero);
- }
+ scale_intensities(model, image.reflections,
+ sym);
}
- /* Start of second or later pattern */
- merge_pattern(model, observed, new_pattern, items,
- counts, config_maxonly,
- twins, holo, mero,
- reference, reference_i,
- hist_vals, hist_h, hist_k, hist_l,
- hist_i, devs, means, outfh);
-
- n_patterns++;
- if ( n_patterns == config_stopafter ) {
- progress_bar(n_total_patterns, n_total_patterns,
- "Merging");
- break;
- }
- progress_bar(n_patterns, n_total_patterns, "Merging");
+ merge_pattern(model, image.reflections, config_maxonly,
+ sym, hist_vals, hist_h, hist_k, hist_l,
+ hist_i);
- /* Reset for the next pattern */
- clear_items(items);
-
- f0_valid = 0;
+ n_used++;
}
- if ( outfh != NULL ) {
- fprintf(outfh, "%s", line);
- }
+ reflist_free(image.reflections);
+ image_feature_list_free(image.features);
+ cell_free(image.indexed_cell);
- if ( strncmp(line, "f0 = ", 5) == 0 ) {
- r = sscanf(line, "f0 = %f", &f0);
- if ( r != 1 ) {
- f0 = 1.0;
- f0_valid = 0;
- continue;
- }
- f0_valid = 1;
- }
+ progress_bar(n_patterns, n_total_patterns-config_startafter,
+ "Merging");
- r = sscanf(line, "%i %i %i %f", &h, &k, &l, &intensity);
- if ( r != 4 ) continue;
+ } while ( rval == 0 );
- /* Not interested in the central beam */
- if ( (h==0) && (k==0) && (l==0) ) continue;
+ /* Divide by counts to get mean intensity if necessary */
+ if ( !config_sum && !config_maxonly ) {
- /* The same raw indices (before mapping into the asymmetric
- * unit should not turn up twice in one pattern. */
- if ( find_item(items, h, k, l) != 0 ) {
- ERROR("More than one measurement for %i %i %i in"
- " pattern number %i\n", h, k, l, n_patterns);
- }
- set_intensity(new_pattern, h, k, l, intensity);
+ Reflection *refl;
+ RefListIterator *iter;
- /* NB: This list contains raw indices, before working out
- * where they belong in the asymmetric unit. */
- add_item(items, h, k, l);
+ for ( refl = first_refl(model, &iter);
+ refl != NULL;
+ refl = next_refl(refl, iter) ) {
- } while ( rval != NULL );
+ double intensity = get_intensity(refl);
+ double sum_squared_dev = get_sum_squared_dev(refl);
+ int red = get_redundancy(refl);
- delete_items(items);
- free(new_pattern);
+ set_int(refl, intensity / (double)red);
+ set_sum_squared_dev(refl, sum_squared_dev/(double)red);
- /* Calculate mean intensity if necessary */
- if ( !config_sum && !config_maxonly ) {
- for ( i=0; i<IDIM*IDIM*IDIM; i++ ) {
- if ( counts[i] > 0 ) {
- model[i] /= (double)counts[i];
- }
}
+
}
- *pmodel = model;
- *pcounts = counts;
- *pobserved = observed;
+ /* Calculate ESDs */
+ for ( refl = first_refl(model, &iter);
+ refl != NULL;
+ refl = next_refl(refl, iter) ) {
+
+ double sum_squared_dev = get_sum_squared_dev(refl);
+ int red = get_redundancy(refl);
+
+ set_esd_intensity(refl, sum_squared_dev/(double)red);
+
+ }
- STATUS("%i patterns had no f0 valid value.\n", n_nof0);
+ STATUS("%i of the patterns could be used.\n", n_used);
}
@@ -546,35 +368,22 @@ int main(int argc, char *argv[])
char *filename = NULL;
char *output = NULL;
FILE *fh;
- double *model;
- unsigned int *counts;
+ RefList *model;
UnitCell *cell = NULL;
int config_maxonly = 0;
int config_startafter = 0;
int config_stopafter = 0;
int config_sum = 0;
int config_scale = 0;
- int config_rmerge = 0;
unsigned int n_total_patterns;
char *sym = NULL;
- char *in_sym = NULL;
char *pdb = NULL;
- ReflItemList *twins;
- ReflItemList *observed;
- int i;
char *histo = NULL;
signed int hist_h, hist_k, hist_l;
double *hist_vals = NULL;
int hist_i;
- char *outstream = NULL;
- char *reference = NULL;
- ReflItemList *reference_items;
- double *reference_i;
- FILE *outfh = NULL;
struct beam_params *beam = NULL;
int space_for_hist = 0;
- double *devs;
- double *esds;
/* Long options */
const struct option longopts[] = {
@@ -588,18 +397,14 @@ int main(int argc, char *argv[])
{"sum", 0, &config_sum, 1},
{"scale", 0, &config_scale, 1},
{"symmetry", 1, NULL, 'y'},
- {"input-symmetry", 1, NULL, 'a'},
{"pdb", 1, NULL, 'p'},
{"histogram", 1, NULL, 'g'},
- {"rmerge", 0, &config_rmerge, 1},
- {"outstream", 1, NULL, 2},
- {"reference", 1, NULL, 'r'},
{"beam", 1, NULL, 'b'},
{0, 0, NULL, 0}
};
/* Short options */
- while ((c = getopt_long(argc, argv, "hi:e:r:o:p:y:g:f:a:b:",
+ while ((c = getopt_long(argc, argv, "hi:e:o:p:y:g:f:b:",
longopts, NULL)) != -1) {
switch (c) {
@@ -635,18 +440,6 @@ int main(int argc, char *argv[])
histo = strdup(optarg);
break;
- case 'r' :
- reference = strdup(optarg);
- break;
-
- case 2 :
- outstream = strdup(optarg);
- break;
-
- case 'a' :
- in_sym = strdup(optarg);
- break;
-
case 'b' :
beam = get_beam_parameters(optarg);
if ( beam == NULL ) {
@@ -665,17 +458,15 @@ int main(int argc, char *argv[])
}
- if ( config_sum && config_rmerge ) {
- ERROR("Options --sum and --rmerge do not make sense"
- " together.\n");
- return 1;
- }
-
if ( filename == NULL ) {
ERROR("Please specify filename using the -i option\n");
return 1;
}
+ if ( output == NULL ) {
+ output = strdup("processed.hkl");
+ }
+
if ( pdb == NULL ) {
pdb = strdup("molecule.pdb");
}
@@ -685,35 +476,6 @@ int main(int argc, char *argv[])
if ( sym == NULL ) sym = strdup("1");
- /* Show useful symmetry information */
- if ( (sym != NULL) && (in_sym != NULL) ) {
-
- int np = num_general_equivs(in_sym) / num_general_equivs(sym);
- if ( np > 1 ) {
-
- STATUS("Resolving point group %s into %s "
- "(%i possibilities)\n",
- in_sym, sym, np);
- /* Get the list of twin/Bijvoet possibilities */
- twins = get_twin_possibilities(in_sym, sym);
- STATUS("Twin operation indices from %s are:", in_sym);
- for ( i=0; i<num_items(twins); i++ ) {
- STATUS(" %i", get_item(twins, i)->op);
- }
- STATUS("\n");
-
- } else {
- STATUS("No resolution required to get from %s to %s\n",
- in_sym, sym);
- twins = NULL;
- }
-
- } else {
- STATUS("No twin resolution requested (use -a otherwise).\n");
- twins = NULL;
- in_sym = strdup(sym);
- }
-
/* Open the data stream */
if ( strcmp(filename, "-") == 0 ) {
fh = stdin;
@@ -726,33 +488,13 @@ int main(int argc, char *argv[])
return 1;
}
- /* Read the reference reflections */
- if ( reference != NULL ) {
- reference_i = new_list_intensity();
- reference_items = read_reflections(reference, reference_i,
- NULL, NULL, NULL);
- if ( reference_items == NULL ) {
- ERROR("Couldn't read '%s'\n", reference);
- return 1;
- }
- } else {
- reference_items = NULL;
- reference_i = NULL;
- }
-
- if ( outstream != NULL ) {
- outfh = fopen(outstream, "w");
- if ( outfh == NULL ) {
- ERROR("Couldn't open '%s'\n", outstream);
- return 1;
- }
- }
-
/* Count the number of patterns in the file */
n_total_patterns = count_patterns(fh);
STATUS("There are %i patterns to process\n", n_total_patterns);
rewind(fh);
+ model = reflist_new();
+
if ( histo != NULL ) {
int r;
r = sscanf(histo, "%i,%i,%i", &hist_h, &hist_k, &hist_l);
@@ -771,13 +513,10 @@ int main(int argc, char *argv[])
}
hist_i = 0;
- merge_all(fh, &model, &observed, &counts,
- config_maxonly, config_scale, config_sum,
+ merge_all(fh, model, config_maxonly, config_scale, config_sum,
config_startafter, config_stopafter,
- twins, in_sym, sym, n_total_patterns,
- reference_items, reference_i,
- hist_vals, hist_h, hist_k, hist_l, &hist_i, NULL, NULL,
- outfh);
+ sym, n_total_patterns,
+ hist_vals, hist_h, hist_k, hist_l, &hist_i);
rewind(fh);
if ( space_for_hist && (hist_i >= space_for_hist) ) {
ERROR("Histogram array was too small!\n");
@@ -790,60 +529,17 @@ int main(int argc, char *argv[])
}
STATUS("Extra pass to calculate ESDs...\n");
- devs = new_list_intensity();
- esds = new_list_sigma();
rewind(fh);
- merge_all(fh, &model, &observed, &counts,
- config_maxonly, config_scale, 0,
- config_startafter, config_stopafter,
- twins, in_sym, sym,
- n_total_patterns, reference_items, reference_i,
- NULL, 0, 0, 0, NULL, devs, model, NULL);
+ merge_all(fh, model, config_maxonly, config_scale, 0,
+ config_startafter, config_stopafter, sym, n_total_patterns,
+ NULL, 0, 0, 0, NULL);
- for ( i=0; i<num_items(observed); i++ ) {
-
- struct refl_item *it;
- signed int h, k, l;
- double dev, esd;
- unsigned int count;
-
- it = get_item(observed, i);
- h = it->h; k = it->k, l = it->l;
-
- dev = lookup_intensity(devs, h, k, l);
- count = lookup_count(counts, h, k, l);
-
- if ( count < 2 ) {
- /* If we have only one measurement, the error is 100% */
- esd = lookup_sigma(model, h, k, l);
- } else {
- esd = sqrt(dev) / (double)count;
- }
- set_sigma(esds, h, k, l, esd);
-
- }
-
- if ( output != NULL ) {
-
- double adu_per_photon;
-
- if ( beam == NULL ) {
- adu_per_photon = 167.0;
- STATUS("No beam parameters file provided (use -b), "
- "so I'm assuming 167.0 ADU per photon.\n");
- } else {
- adu_per_photon = beam->adu_per_photon;
- }
-
- write_reflections(output, observed, model, esds,
- NULL, counts, cell);
-
- }
+ write_reflist(output, model, cell);
+ cell_free(cell);
fclose(fh);
free(sym);
- free(model);
- free(counts);
+ reflist_free(model);
free(output);
if ( cell != NULL ) cell_free(cell);
diff --git a/src/reflist-utils.c b/src/reflist-utils.c
new file mode 100644
index 00000000..33d373b8
--- /dev/null
+++ b/src/reflist-utils.c
@@ -0,0 +1,72 @@
+/*
+ * reflist-utils.c
+ *
+ * Utilities to complement the core reflist.c
+ *
+ * (c) 2006-2011 Thomas White <taw@physics.org>
+ *
+ * Part of CrystFEL - crystallography with a FEL
+ *
+ */
+
+
+#include <stdio.h>
+
+
+#include "reflist.h"
+#include "cell.h"
+
+
+void write_reflections_to_file(FILE *fh, RefList *list, UnitCell *cell)
+{
+ Reflection *refl;
+ RefListIterator *iter;
+
+ fprintf(fh, " h k l I phase sigma(I) "
+ " 1/d(nm^-1) counts fs/px ss/px\n");
+
+ for ( refl = first_refl(list, &iter);
+ refl != NULL;
+ refl = next_refl(refl, iter) ) {
+
+ signed int h, k, l;
+ double intensity, esd_i, s;
+ double fs, ss;
+
+ get_indices(refl, &h, &k, &l);
+ get_detector_pos(refl, &fs, &ss);
+ intensity = get_intensity(refl);
+ esd_i = 0.0; /* FIXME! */
+ s = 0.0; /* FIXME! */
+
+ /* h, k, l, I, sigma(I), s */
+ fprintf(fh,
+ "%3i %3i %3i %10.2f %s %10.2f %10.2f %7i %6.1f %6.1f\n",
+ h, k, l, intensity, " -", esd_i, s/1.0e9, 1,
+ fs, ss);
+
+ }
+}
+
+
+int write_reflist(const char *filename, RefList *list, UnitCell *cell)
+{
+ FILE *fh;
+
+ if ( filename == NULL ) {
+ fh = stdout;
+ } else {
+ fh = fopen(filename, "w");
+ }
+
+ if ( fh == NULL ) {
+ ERROR("Couldn't open output file '%s'.\n", filename);
+ return 1;
+ }
+
+ write_reflections_to_file(fh, list, cell);
+
+ fclose(fh);
+
+ return 0;
+}
diff --git a/src/reflist-utils.h b/src/reflist-utils.h
new file mode 100644
index 00000000..d5c1e655
--- /dev/null
+++ b/src/reflist-utils.h
@@ -0,0 +1,30 @@
+/*
+ * reflist-utils.h
+ *
+ * Utilities to complement the core reflist.c
+ *
+ * (c) 2006-2011 Thomas White <taw@physics.org>
+ *
+ * Part of CrystFEL - crystallography with a FEL
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#ifndef REFLIST_UTILS_H
+#define REFLIST_UTILS_H
+
+
+#include "reflist.h"
+#include "cell.h"
+
+
+extern void write_reflections_to_file(FILE *fh, RefList *list, UnitCell *cell);
+
+extern int write_reflist(const char *filename, RefList *list,
+ UnitCell *cell);
+
+
+#endif /* REFLIST_UTILS_H */
diff --git a/src/reflist.c b/src/reflist.c
index 3a53330c..3181b418 100644
--- a/src/reflist.c
+++ b/src/reflist.c
@@ -45,6 +45,13 @@ struct _refldata {
/* Intensity */
double intensity;
double esd_i;
+
+ /* Redundancy */
+ int redundancy;
+
+ /* Total squared difference between all estimates of this reflection
+ * and the estimated mean value */
+ double sum_squared_dev;
};
@@ -225,6 +232,24 @@ int get_scalable(Reflection *refl)
}
+int get_redundancy(Reflection *refl)
+{
+ return refl->data.redundancy;
+}
+
+
+double get_sum_squared_dev(Reflection *refl)
+{
+ return refl->data.sum_squared_dev;
+}
+
+
+double get_esd_intensity(Reflection *refl)
+{
+ return refl->data.esd_i;
+}
+
+
/********************************** Setters ***********************************/
void set_detector_pos(Reflection *refl, double exerr, double x, double y)
@@ -258,6 +283,24 @@ void set_scalable(Reflection *refl, int scalable)
}
+void set_redundancy(Reflection *refl, int red)
+{
+ refl->data.redundancy = red;
+}
+
+
+void set_sum_squared_dev(Reflection *refl, double dev)
+{
+ refl->data.sum_squared_dev = dev;
+}
+
+
+void set_esd_intensity(Reflection *refl, double esd)
+{
+ refl->data.esd_i = esd;
+}
+
+
/********************************* Insertion **********************************/
static void insert_node(Reflection *head, Reflection *new)
diff --git a/src/reflist.h b/src/reflist.h
index e25a16ad..f6d01ecf 100644
--- a/src/reflist.h
+++ b/src/reflist.h
@@ -41,6 +41,9 @@ extern double get_intensity(Reflection *refl);
extern void get_partial(Reflection *refl, double *r1, double *r2, double *p,
int *clamp_low, int *clamp_high);
extern int get_scalable(Reflection *refl);
+extern int get_redundancy(Reflection *refl);
+extern double get_sum_squared_dev(Reflection *refl);
+extern double get_esd_intensity(Reflection *refl);
/* Set */
extern void set_detector_pos(Reflection *refl, double exerr,
@@ -49,6 +52,9 @@ extern void set_partial(Reflection *refl, double r1, double r2, double p,
double clamp_low, double clamp_high);
extern void set_int(Reflection *refl, double intensity);
extern void set_scalable(Reflection *refl, int scalable);
+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);
/* Insertion */
extern Reflection *add_refl(RefList *list, INDICES);
diff --git a/src/stream.c b/src/stream.c
index 2a75216d..4b2bea9f 100644
--- a/src/stream.c
+++ b/src/stream.c
@@ -23,6 +23,7 @@
#include "image.h"
#include "stream.h"
#include "reflist.h"
+#include "reflist-utils.h"
#define CHUNK_START_MARKER "----- Begin chunk -----"
@@ -32,7 +33,6 @@
#define REFLECTION_START_MARKER "Reflections measured after indexing"
#define REFLECTION_END_MARKER "End of reflections"
-
static void exclusive(const char *a, const char *b)
{
ERROR("The stream options '%s' and '%s' are mutually exclusive.\n",
@@ -103,10 +103,9 @@ int count_patterns(FILE *fh)
rval = fgets(line, 1023, fh);
if ( rval == NULL ) continue;
- if ( (strncmp(line, "Reflections from indexing", 25) == 0)
- || (strncmp(line, "New pattern", 11) == 0) ) {
- n_total_patterns++;
- }
+ chomp(line);
+ if ( strcmp(line, CHUNK_END_MARKER) == 0 ) n_total_patterns++;
+
} while ( rval != NULL );
return n_total_patterns;
@@ -147,6 +146,7 @@ static UnitCell *read_orientation_matrix(FILE *fh)
static int read_reflections(FILE *fh, struct image *image)
{
char *rval = NULL;
+ int first = 1;
image->reflections = reflist_new();
@@ -164,17 +164,20 @@ static int read_reflections(FILE *fh, struct image *image)
if ( rval == NULL ) continue;
chomp(line);
- if ( strcmp(line, PEAK_LIST_END_MARKER) == 0 ) return 0;
+ if ( strcmp(line, REFLECTION_END_MARKER) == 0 ) return 0;
r = sscanf(line, "%i %i %i %f %s %f %f %i %f %f",
&h, &k, &l, &intensity, phs, &sigma, &res, &cts,
&fs, &ss);
- if ( r != 10 ) return 1;
-
- refl = add_refl(image->reflections, h, k, l);
- set_int(refl, intensity);
- set_detector_pos(refl, fs, ss, 0.0);
- /* FIXME: Set ESD */
+ if ( (r != 10) && (!first) ) return 1;
+
+ first = 0;
+ if ( r == 10 ) {
+ refl = add_refl(image->reflections, h, k, l);
+ set_int(refl, intensity);
+ set_detector_pos(refl, fs, ss, 0.0);
+ set_esd_intensity(refl, sigma);
+ }
} while ( rval != NULL );
@@ -184,45 +187,10 @@ static int read_reflections(FILE *fh, struct image *image)
}
-static void write_reflections(struct image *image, FILE *ofh)
-{
- Reflection *refl;
- RefListIterator *iter;
-
- fprintf(ofh, REFLECTION_START_MARKER"\n");
- /* FIXME: Unify this with write_reflections() over in reflections.c */
- fprintf(ofh, " h k l I phase sigma(I) "
- " 1/d(nm^-1) counts fs/px ss/px\n");
-
- for ( refl = first_refl(image->reflections, &iter);
- refl != NULL;
- refl = next_refl(refl, iter) ) {
-
- signed int h, k, l;
- double intensity, esd_i, s;
- double fs, ss;
-
- get_indices(refl, &h, &k, &l);
- get_detector_pos(refl, &fs, &ss);
- intensity = get_intensity(refl);
- esd_i = 0.0; /* FIXME! */
- s = 0.0; /* FIXME! */
-
- /* h, k, l, I, sigma(I), s */
- fprintf(ofh,
- "%3i %3i %3i %10.2f %s %10.2f %10.2f %7i %6.1f %6.1f\n",
- h, k, l, intensity, " -", esd_i, s/1.0e9, 1,
- fs, ss);
-
- }
-
- fprintf(ofh, REFLECTION_END_MARKER"\n");
-}
-
-
static int read_peaks(FILE *fh, struct image *image)
{
char *rval = NULL;
+ int first = 1;
image->features = image_feature_list_new();
@@ -239,9 +207,17 @@ static int read_peaks(FILE *fh, struct image *image)
if ( strcmp(line, PEAK_LIST_END_MARKER) == 0 ) return 0;
r = sscanf(line, "%f %f %f %f", &x, &y, &d, &intensity);
- if ( r != 4 ) return 1;
+ if ( (r != 4) && (!first) ) {
+ ERROR("Failed to parse peak list line.\n");
+ ERROR("The failed line was: '%s'\n", line);
+ return 1;
+ }
- image_add_feature(image->features, x, y, image, 1.0, NULL);
+ first = 0;
+ if ( r == 4 ) {
+ image_add_feature(image->features, x, y,
+ image, 1.0, NULL);
+ }
} while ( rval != NULL );
@@ -255,7 +231,7 @@ static void write_peaks(struct image *image, FILE *ofh)
int i;
fprintf(ofh, PEAK_LIST_START_MARKER"\n");
- fprintf(ofh, " fs/px ss/px (1/d)/nm^-1 Intensity\n");
+ fprintf(ofh, " fs/px ss/px (1/d)/nm^-1 Intensity\n");
for ( i=0; i<image_feature_count(image->features); i++ ) {
@@ -269,7 +245,7 @@ static void write_peaks(struct image *image, FILE *ofh)
r = get_q(image, f->fs, f->ss, NULL, 1.0/image->lambda);
q = modulus(r.u, r.v, r.w);
- fprintf(ofh, "%8.3f %8.3f %8.3f %12.3f\n",
+ fprintf(ofh, "%6.1f %6.1f %10.2f %10.2f\n",
f->fs, f->ss, q/1.0e9, f->intensity);
}
@@ -330,8 +306,12 @@ void write_chunk(FILE *ofh, struct image *i, int f)
}
if ( (f & STREAM_PIXELS) || (f & STREAM_INTEGRATED) ) {
+
fprintf(ofh, "\n");
- write_reflections(i, ofh);
+ fprintf(ofh, REFLECTION_START_MARKER"\n");
+ write_reflections_to_file(ofh, i->reflections, i->indexed_cell);
+ fprintf(ofh, REFLECTION_END_MARKER"\n");
+
}
fprintf(ofh, CHUNK_END_MARKER"\n\n");
@@ -374,10 +354,11 @@ int read_chunk(FILE *fh, struct image *image)
if ( find_start_of_chunk(fh) ) return 1;
image->i0_available = 0;
- if ( image->features != NULL ) {
- image_feature_list_free(image->features);
- image->features = NULL;
- }
+ image->i0 = 1.0;
+ image->lambda = -1.0;
+ image->features = NULL;
+ image->reflections = NULL;
+ image->indexed_cell = NULL;
do {
@@ -426,16 +407,22 @@ int read_chunk(FILE *fh, struct image *image)
}
if ( strcmp(line, PEAK_LIST_START_MARKER) == 0 ) {
- if ( read_peaks(fh, image) ) return 1;
+ if ( read_peaks(fh, image) ) {
+ ERROR("Failed while reading peaks\n");
+ return 1;
+ }
}
if ( strcmp(line, REFLECTION_START_MARKER) == 0 ) {
- if ( read_reflections(fh, image) ) return 1;
+ if ( read_reflections(fh, image) ) {
+ ERROR("Failed while reading reflections\n");
+ return 1;
+ }
}
} while ( strcmp(line, CHUNK_END_MARKER) != 0 );
- if ( have_filename && have_cell && have_ev ) return 0;
+ if ( have_filename && have_ev ) return 0;
ERROR("Incomplete chunk found in input file.\n");
return 1;
@@ -496,3 +483,24 @@ int find_chunk(FILE *fh, UnitCell **cell, char **filename, double *ev)
return 1;
}
+
+
+int skip_some_files(FILE *fh, int n)
+{
+ char *rval = NULL;
+ int n_patterns = 0;
+
+ do {
+
+ char line[1024];
+
+ if ( n_patterns == n ) return 0;
+
+ rval = fgets(line, 1023, fh);
+ if ( rval == NULL ) continue;
+ if ( strcmp(line, CHUNK_END_MARKER) == 0 ) n_patterns++;
+
+ } while ( rval != NULL );
+
+ return 1;
+}
diff --git a/src/stream.h b/src/stream.h
index 3869e515..64bd9529 100644
--- a/src/stream.h
+++ b/src/stream.h
@@ -38,4 +38,8 @@ extern void write_chunk(FILE *ofh, struct image *image, int flags);
extern int parse_stream_flags(const char *a);
+extern int read_chunk(FILE *fh, struct image *image);
+
+extern int skip_some_files(FILE *fh, int n);
+
#endif /* STREAM_H */
generated by cgit v1.2.3 (git 2.25.1) at 2024-07-09 05:30:53 +0000