From ee2a9703c723d807c4e590906ab6a45868651a81 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Mon, 3 Mar 2014 15:44:26 +0100 Subject: Include background and peak height in stream --- libcrystfel/src/integration.c | 18 +++++++ libcrystfel/src/reflist.c | 55 ++++++++++++++++++++ libcrystfel/src/reflist.h | 4 ++ libcrystfel/src/stream.c | 115 +++++++++++++++++++++++++++++++++++++++--- 4 files changed, 185 insertions(+), 7 deletions(-) diff --git a/libcrystfel/src/integration.c b/libcrystfel/src/integration.c index b86399e9..bcd0a88b 100644 --- a/libcrystfel/src/integration.c +++ b/libcrystfel/src/integration.c @@ -216,6 +216,9 @@ struct peak_box double sigma; double J; /* Profile scaling factor */ + /* Highest (non-ignored) value in peak */ + double peak; + /* Offsets to final observed position */ double offs_fs; double offs_ss; @@ -841,6 +844,7 @@ static int check_box(struct intcontext *ic, struct peak_box *bx, int *sat) &cdx, &cdy, &cdz); get_indices(bx->refl, &hr, &kr, &lr); + bx->peak = -INFINITY; for ( p=0; pw; p++ ) { for ( q=0; qw; q++ ) { @@ -874,6 +878,12 @@ static int check_box(struct intcontext *ic, struct peak_box *bx, int *sat) if ( sat != NULL ) *sat = 1; } + if ( (bx->bm[p+ic->w*q] != BM_IG) + && (bx->bm[p+ic->w*q] != BM_BH) + && (boxi(ic, bx, p, q) > bx->peak) ) { + bx->peak = boxi(ic, bx, p, q); + } + /* Ignore if this pixel is closer to the next reciprocal lattice * point */ dv = get_q_for_panel(bx->p, fs, ss, NULL, ic->k); @@ -1342,9 +1352,15 @@ static void integrate_prof2d_once(struct intcontext *ic, struct peak_box *bx) if ( bg_ok(bx) ) { double pfs, pss; + double bgmean; + double sig2_bg; /* unused */ + + mean_var_area(ic, bx, BM_BG, &bgmean, &sig2_bg); set_intensity(bx->refl, bx->intensity); set_esd_intensity(bx->refl, bx->sigma); + set_peak(bx->refl, bx->peak); + set_mean_bg(bx->refl, bgmean); set_redundancy(bx->refl, 1); /* Update position */ @@ -1600,6 +1616,8 @@ static void integrate_rings_once(Reflection *refl, struct image *image, set_intensity(refl, intensity); set_esd_intensity(refl, sigma); set_redundancy(refl, 1); + set_mean_bg(refl, bgmean); + set_peak(refl, bx->peak); get_indices(refl, &h, &k, &l); one_over_d = resolution(cell, h, k, l); diff --git a/libcrystfel/src/reflist.c b/libcrystfel/src/reflist.c index f6fad33d..7bc98788 100644 --- a/libcrystfel/src/reflist.c +++ b/libcrystfel/src/reflist.c @@ -100,6 +100,10 @@ struct _refldata { /* Redundancy */ int redundancy; + /* Peak height and mean background */ + double peak; + double mean_bg; + /* User-specified temporary values */ double temp1; double temp2; @@ -521,6 +525,33 @@ double get_phase(const Reflection *refl, int *have_phase) } +/** + * get_peak: + * @refl: A %Reflection + * + * Returns: the peak height (value of the highest pixel, before background + * subtraction) for this reflection. + * + **/ +double get_peak(const Reflection *refl) +{ + return refl->data.peak; +} + + +/** + * get_mean_bg: + * @refl: A %Reflection + * + * Returns: the mean background level for this reflection. + * + **/ +double get_mean_bg(const Reflection *refl) +{ + return refl->data.mean_bg; +} + + /** * get_temp1: * @refl: A %Reflection @@ -720,6 +751,30 @@ void set_phase(Reflection *refl, double phase) } +/** + * set_peak: + * @refl: A %Reflection + * @peak: New peak height for the reflection + * + **/ +void set_peak(Reflection *refl, double peak) +{ + refl->data.peak = peak; +} + + +/** + * set_mean_bg: + * @refl: A %Reflection + * @mean_bg: New peak height for the reflection + * + **/ +void set_mean_bg(Reflection *refl, double mean_bg) +{ + refl->data.mean_bg = mean_bg; +} + + /** * set_symmetric_indices: * @refl: A %Reflection diff --git a/libcrystfel/src/reflist.h b/libcrystfel/src/reflist.h index 7969235c..61cac9d5 100644 --- a/libcrystfel/src/reflist.h +++ b/libcrystfel/src/reflist.h @@ -94,6 +94,8 @@ extern double get_temp1(const Reflection *refl); extern double get_temp2(const Reflection *refl); extern double get_esd_intensity(const Reflection *refl); extern double get_phase(const Reflection *refl, int *have_phase); +extern double get_peak(const Reflection *refl); +extern double get_mean_bg(const Reflection *refl); /* Set */ extern void copy_data(Reflection *to, const Reflection *from); @@ -111,6 +113,8 @@ extern void set_temp1(Reflection *refl, double temp); extern void set_temp2(Reflection *refl, double temp); extern void set_esd_intensity(Reflection *refl, double esd); extern void set_phase(Reflection *refl, double phase); +extern void set_peak(Reflection *refl, double pk); +extern void set_mean_bg(Reflection *refl, double bg); extern void set_symmetric_indices(Reflection *refl, signed int hs, signed int ks, signed int ls); diff --git a/libcrystfel/src/stream.c b/libcrystfel/src/stream.c index 4cf02414..b204f9a9 100644 --- a/libcrystfel/src/stream.c +++ b/libcrystfel/src/stream.c @@ -51,7 +51,10 @@ #include "reflist-utils.h" #define LATEST_MAJOR_VERSION (2) -#define LATEST_MINOR_VERSION (1) +#define LATEST_MINOR_VERSION (2) + +#define AT_LEAST_VERSION(st, a, b) ((st->major_version>=(a)) \ + && (st->minor_version>=(b))) struct _stream @@ -130,6 +133,91 @@ static void write_peaks(struct image *image, FILE *ofh) } +static RefList *read_stream_reflections(FILE *fh) +{ + char *rval = NULL; + int first = 1; + RefList *out; + + out = reflist_new(); + + do { + + char line[1024]; + signed int h, k, l; + float intensity, sigma, fs, ss, pk, bg; + int r; + Reflection *refl; + + rval = fgets(line, 1023, fh); + if ( rval == NULL ) continue; + chomp(line); + + if ( strcmp(line, REFLECTION_END_MARKER) == 0 ) return out; + + r = sscanf(line, "%i %i %i %f %f %f %f %f %f", + &h, &k, &l, &intensity, &sigma, &fs, &ss, &pk, &bg); + if ( (r != 9) && (!first) ) { + reflist_free(out); + return NULL; + } + + first = 0; + if ( r == 9 ) { + + refl = add_refl(out, h, k, l); + set_intensity(refl, intensity); + set_detector_pos(refl, 0.0, fs, ss); + set_esd_intensity(refl, sigma); + set_redundancy(refl, 1); + set_peak(refl, pk); + set_mean_bg(refl, bg); + + } + + } while ( rval != NULL ); + + /* Got read error of some kind before finding PEAK_LIST_END_MARKER */ + return NULL; +} + + +static void write_stream_reflections(FILE *fh, RefList *list) +{ + Reflection *refl; + RefListIterator *iter; + + fprintf(fh, " h k l I sigma(I) fs/px ss/px" + " peak background\n"); + + for ( refl = first_refl(list, &iter); + refl != NULL; + refl = next_refl(refl, iter) ) + { + + signed int h, k, l; + double intensity, esd_i, bg, pk; + double fs, ss; + + get_indices(refl, &h, &k, &l); + get_detector_pos(refl, &fs, &ss); + intensity = get_intensity(refl); + esd_i = get_esd_intensity(refl); + pk = get_peak(refl); + bg = get_mean_bg(refl); + + /* Reflections with redundancy = 0 are not written */ + if ( get_redundancy(refl) == 0 ) continue; + + fprintf(fh, + "%4i %4i %4i %10.2f %10.2f %6.1f %6.1f %10.2f %10.2f\n", + h, k, l, intensity, esd_i, fs, ss, pk, bg); + + } + +} + + static void write_crystal(Stream *st, Crystal *cr, int include_reflections) { UnitCell *cell; @@ -185,7 +273,11 @@ static void write_crystal(Stream *st, Crystal *cr, int include_reflections) if ( reflist != NULL ) { fprintf(st->fh, REFLECTION_START_MARKER"\n"); - write_reflections_to_file(st->fh, reflist); + if ( AT_LEAST_VERSION(st, 2, 2) ) { + write_stream_reflections(st->fh, reflist); + } else { + write_reflections_to_file(st->fh, reflist); + } fprintf(st->fh, REFLECTION_END_MARKER"\n"); } else { @@ -366,16 +458,21 @@ void read_crystal(Stream *st, struct image *image) if ( strcmp(line, REFLECTION_START_MARKER) == 0 ) { - RefList *reflections; - - reflections = read_reflections_from_file(st->fh); + RefList *reflist; - if ( reflections == NULL ) { + /* The reflection list format in the stream diverges + * after 2.2 */ + if ( AT_LEAST_VERSION(st, 2, 2) ) { + reflist = read_stream_reflections(st->fh); + } else { + reflist = read_reflections_from_file(st->fh); + } + if ( reflist == NULL ) { ERROR("Failed while reading reflections\n"); break; } - crystal_set_reflections(cr, reflections); + crystal_set_reflections(cr, reflist); } @@ -564,6 +661,9 @@ Stream *open_stream_for_read(const char *filename) } else if ( strncmp(line, "CrystFEL stream format 2.1", 26) == 0 ) { st->major_version = 2; st->minor_version = 1; + } else if ( strncmp(line, "CrystFEL stream format 2.2", 26) == 0 ) { + st->major_version = 2; + st->minor_version = 2; } else { ERROR("Invalid stream, or stream format is too new.\n"); close_stream(st); @@ -634,6 +734,7 @@ int is_stream(const char *filename) if ( strncmp(line, "CrystFEL stream format 2.0", 26) == 0 ) return 1; if ( strncmp(line, "CrystFEL stream format 2.1", 26) == 0 ) return 1; + if ( strncmp(line, "CrystFEL stream format 2.2", 26) == 0 ) return 1; return 0; } -- cgit v1.2.3