Include background and peak height in stream

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; p<ic->w; p++ ) {
 	for ( q=0; q<ic->w; 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;
@@ -522,6 +526,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
  *
@@ -721,6 +752,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
  * @hs: The 'h' index of the 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;
 }