Loads of scaling report stuff

3 files changed, 387 insertions, 28 deletions

diff --git a/src/partialator.c b/src/partialator.c
index 7b8eefa1..19c1b058 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -282,6 +282,8 @@ int main(int argc, char *argv[])
 	RefList *reference = NULL;
 	int n_dud;
 	int have_reference = 0;
+	char cmdline[1024];
+	SRContext *sr;
 
 	/* Long options */
 	const struct option longopts[] = {
@@ -296,6 +298,12 @@ int main(int argc, char *argv[])
 		{0, 0, NULL, 0}
 	};
 
+	cmdline[0] = '\0';
+	for ( i=1; i<argc; i++ ) {
+		strncat(cmdline, argv[i], 1023-strlen(cmdline));
+		strncat(cmdline, " ", 1023-strlen(cmdline));
+	}
+
 	/* Short options */
 	while ((c = getopt_long(argc, argv, "hi:g:x:j:y:o:b:",
 	                        longopts, NULL)) != -1)
@@ -414,6 +422,8 @@ int main(int argc, char *argv[])
 		return 1;
 	}
 
+	sr = sr_header("scaling-report.pdf", infile, cmdline);
+
 	/* Fill in what we know about the images so far */
 	rewind(fh);
 	nobs = 0;
@@ -480,6 +490,8 @@ int main(int argc, char *argv[])
 	STATUS("Performing initial scaling.\n");
 	full = scale_intensities(images, n_usable_patterns, reference);
 
+	sr_before(sr, images, n_usable_patterns, full);
+
 	/* Iterate */
 	for ( i=0; i<n_iter; i++ ) {
 
@@ -537,18 +549,16 @@ int main(int argc, char *argv[])
 
 	}
 
-	STATUS("Final scale factors:\n");
 	n_dud = 0;
 	for ( i=0; i<n_usable_patterns; i++ ) {
 		if ( images[i].pr_dud ) n_dud++;
-		STATUS("%4i : %5.2f\n", i, images[i].osf);
 	}
 	STATUS("%i images could not be refined on the last cycle.\n", n_dud);
 
 	/* Output results */
 	write_reflist(outfile, full, images[0].indexed_cell);
 
-	scaling_report("scaling-report.pdf", images, n_usable_patterns, infile);
+	sr_after(sr, images, n_usable_patterns, full);
 
 	/* Clean up */
 	for ( i=0; i<n_usable_patterns; i++ ) {
diff --git a/src/scaling-report.c b/src/scaling-report.c
index de0e396e..eb4429e4 100644
--- a/src/scaling-report.c
+++ b/src/scaling-report.c
@@ -16,54 +16,384 @@
 #include <cairo.h>
 #include <cairo-pdf.h>
 #include <pango/pangocairo.h>
+#include <math.h>
 
 #include "image.h"
+#include "scaling-report.h"
 
 
-static void write_title(cairo_t *cr, const char *filename, double w, double h)
+#define PAGE_WIDTH (842.0)
+
+enum justification
+{
+	J_CENTER,
+	J_LEFT,
+	J_RIGHT,
+};
+
+
+struct _srcontext
+{
+	cairo_surface_t *surf;
+	cairo_t *cr;
+	double w;
+	double h;
+};
+
+
+static void show_text(cairo_t *cr, const char *text, double y,
+                      enum justification j, char *font)
 {
-	char text[1024];
 	PangoLayout *layout;
 	PangoFontDescription *fontdesc;
 	int width, height;
+	PangoAlignment just;
 
-	snprintf(text, 1023, "Scaling report: %s", filename);
+	if ( font == NULL ) font = "Sans 10";
 
 	layout = pango_cairo_create_layout(cr);
+	pango_layout_set_ellipsize(layout, PANGO_ELLIPSIZE_NONE);
+	pango_layout_set_width(layout, PANGO_SCALE*(PAGE_WIDTH-20.0));
+
+	switch ( j )
+	{
+		case J_CENTER : just = PANGO_ALIGN_CENTER; break;
+		case J_LEFT : just = PANGO_ALIGN_LEFT; break;
+		case J_RIGHT : just = PANGO_ALIGN_RIGHT; break;
+		default: just = PANGO_ALIGN_LEFT; break;
+	}
+
+	pango_layout_set_alignment(layout, just);
+	pango_layout_set_wrap(layout, PANGO_WRAP_CHAR);
+	pango_layout_set_spacing(layout, 4.0*PANGO_SCALE);
+
 	pango_layout_set_text(layout, text, -1);
-	fontdesc = pango_font_description_from_string("Sans 14 Bold");
+
+	fontdesc = pango_font_description_from_string(font);
 	pango_layout_set_font_description(layout, fontdesc);
 
 	pango_cairo_update_layout(cr, layout);
 	pango_layout_get_size(layout, &width, &height);
 
-	cairo_move_to(cr, 0.5-width/PANGO_SCALE, 10.0);
+	cairo_move_to(cr, 10.0, y);
+	cairo_set_source_rgb(cr, 0.0, 0.0, 0.0);
+	pango_cairo_show_layout(cr, layout);
+}
 
+
+static void show_text_simple(cairo_t *cr, const char *text, double x, double y,
+                             char *font, double rot, enum justification j)
+{
+	PangoLayout *layout;
+	PangoFontDescription *fontdesc;
+	int width, height;
+
+	cairo_save(cr);
+
+	if ( font == NULL ) font = "Sans 10";
+
+	layout = pango_cairo_create_layout(cr);
+	pango_layout_set_ellipsize(layout, PANGO_ELLIPSIZE_NONE);
+	pango_layout_set_alignment(layout, PANGO_ALIGN_LEFT);
+
+	pango_layout_set_text(layout, text, -1);
+
+	fontdesc = pango_font_description_from_string(font);
+	pango_layout_set_font_description(layout, fontdesc);
+
+	pango_cairo_update_layout(cr, layout);
+	pango_layout_get_size(layout, &width, &height);
+
+	cairo_new_path(cr);
+	cairo_translate(cr, x, y);
+	cairo_rotate(cr, rot);
+	if ( j == J_CENTER ) {
+		cairo_translate(cr, -(width/2.0)/PANGO_SCALE,
+		                    -(height/2.0)/PANGO_SCALE);
+	} else if ( j == J_RIGHT ) {
+		cairo_translate(cr, -width/PANGO_SCALE,
+		                    -(height/2.0)/PANGO_SCALE);
+	} else {
+		cairo_translate(cr, 0.0, -(height/2.0)/PANGO_SCALE);
+	}
 	cairo_set_source_rgb(cr, 0.0, 0.0, 0.0);
 	pango_cairo_show_layout(cr, layout);
+
+	cairo_restore(cr);
 }
 
 
-void scaling_report(const char *filename, const struct image *images, int n,
-                    const char *stream_filename)
+static void plot_point(cairo_t *cr, double g_width, double g_height,
+                       double pcalc, double pobs)
 {
-	cairo_surface_t *surface;
-	cairo_t *cr;
-	const double w = 842.0;
-	const double h = 595.0;
+	int bad = 0;
+
+	if ( pobs > 1.0 ) {
+		pobs = 1.01;
+		bad = 1;
+
+	}
+	if ( pcalc > 1.0 ) {
+		pcalc = 1.01;
+		bad = 1;
+	}
+
+	if ( bad ) {
+		cairo_set_source_rgb(cr, 1.0, 0.0, 0.0);
+	}
+
+	cairo_arc(cr, g_width*pcalc, g_height*(1.0-pobs),
+		      1.0, 0.0, 2.0*M_PI);
+	cairo_fill(cr);
+
+	if ( bad ) {
+		cairo_set_source_rgb(cr, 0.0, 0.7, 0.0);
+	}
+}
+
+static void partiality_graph(cairo_t *cr, const struct image *images, int n,
+                             RefList *full)
+{
+	const double g_width = 200.0;
+	const double g_height = 200.0;
+	int i;
+	const int nbins = 50;
+	double totals[nbins];
+	int counts[nbins];
+	double prob;
+
+	show_text_simple(cr, "Observed partiality", -20.0, g_height/2.0,
+	                      NULL, -M_PI_2, J_CENTER);
+	show_text_simple(cr, "Calculated partiality", g_width/2.0, g_height+20.0,
+	                      NULL, 0.0, J_CENTER);
+
+	for ( i=0; i<nbins; i++ ) {
+		totals[i] = 0.0;
+		counts[i] = 0;
+	}
+
+	cairo_set_source_rgb(cr, 0.0, 0.7, 0.0);
+	prob = 1.0 / n;
+	for ( i=0; i<n; i++ ) {
+
+		Reflection *refl;
+		RefListIterator *iter;
+
+		if ( images[i].pr_dud ) continue;
+
+		for ( refl = first_refl(images[i].reflections, &iter);
+		      refl != NULL;
+		      refl = next_refl(refl, iter) )
+		{
+			double Ipart, Ifull, pobs, pcalc;
+			signed int h, k, l;
+			Reflection *f;
+			int bin;
+
+			if ( !get_scalable(refl) ) continue;
+
+			get_indices(refl, &h, &k, &l);
+			f = find_refl(full, h, k, l);
+			if ( f == NULL ) continue;
+
+			Ipart = get_intensity(refl);
+			Ifull = get_intensity(f);
+
+			pobs = Ipart/Ifull;
+			pcalc = get_partiality(refl);
+
+			bin = nbins * pcalc;
+			totals[bin] += pobs;
+			counts[bin]++;
+
+			if ( (double)random()/RAND_MAX < prob ) {
+				plot_point(cr, g_width, g_height, pcalc, pobs);
+			}
 
-	surface = cairo_pdf_surface_create(filename, w, h);
+		}
 
-	if ( cairo_surface_status(surface) != CAIRO_STATUS_SUCCESS ) {
+	}
+
+	cairo_new_path(cr);
+	cairo_rectangle(cr, 0.0, 0.0, g_width, g_height);
+	cairo_clip(cr);
+
+	cairo_new_path(cr);
+	cairo_move_to(cr, 0.0, g_height);
+	for ( i=0; i<nbins; i++ ) {
+		if ( counts[i] == 0 ) continue;
+		cairo_line_to(cr, g_width*(double)i/nbins,
+		                  g_height - g_height*(totals[i]/counts[i]));
+	}
+	cairo_set_source_rgb(cr, 0.0, 0.0, 0.0);
+	cairo_stroke(cr);
+
+	cairo_reset_clip(cr);
+
+	cairo_new_path(cr);
+	cairo_rectangle(cr, 0.0, 0.0, g_width, g_height);
+	cairo_set_source_rgb(cr, 0.0, 0.0, 0.0);
+	cairo_set_line_width(cr, 1.5);
+	cairo_stroke(cr);
+}
+
+
+static void scale_factor_histogram(cairo_t *cr, const struct image *images,
+                                   int n, const char *title)
+{
+	int f_max;
+	int i, b;
+	const int nbins = 100;
+	double osf_max, osf_inc;
+	double osf_low[nbins];
+	double osf_high[nbins];
+	int counts[nbins];
+	const double g_width = 320.0;
+	const double g_height = 200.0;
+	char tmp[32];
+
+	show_text_simple(cr, title, g_width/2.0, -18.0,
+	                      "Sans Bold 10", 0.0, J_CENTER);
+
+	show_text_simple(cr, "Frequency", -50.0, g_height/2.0,
+	                      NULL, -M_PI_2, J_CENTER);
+
+	osf_max = 0.0;
+	for ( i=0; i<n; i++ ) {
+		double osf = images[i].osf;
+		if ( osf > osf_max ) osf_max = osf;
+	}
+	osf_max = ceil(osf_max);
+	osf_inc = osf_max / nbins;
+
+	for ( b=0; b<nbins; b++ ) {
+		osf_low[b] = b*osf_inc;
+		osf_high[b] = (b+1)*osf_inc;
+		counts[b] = 0;
+	}
+
+	for ( i=0; i<n; i++ ) {
+
+		double osf = images[i].osf;
+
+		for ( b=0; b<nbins; b++ ) {
+			if ( (osf > osf_low[b]) && (osf < osf_high[b]) ) {
+				counts[b]++;
+				break;
+			}
+		}
+	}
+
+	f_max = 0;
+	for ( b=0; b<nbins; b++ ) {
+		if ( counts[b] > f_max ) f_max = counts[b];
+	}
+	f_max = (f_max/10)*10 + 10;
+
+	show_text_simple(cr, "0", -10.0, g_height, NULL, 0.0, J_RIGHT);
+	snprintf(tmp, 31, "%i", f_max);
+	show_text_simple(cr, tmp, -10.0, 0.0, NULL, 0.0, J_RIGHT);
+
+	show_text_simple(cr, "0.00", 0.0, g_height+10.0,
+	                     NULL, -M_PI/3.0, J_RIGHT);
+	snprintf(tmp, 32, "%5.2f", osf_max);
+	show_text_simple(cr, tmp, g_width, g_height+10.0,
+	                     NULL, -M_PI/3.0, J_RIGHT);
+
+	for ( b=0; b<nbins; b++ ) {
+
+		double bar_height;
+
+		bar_height = ((double)counts[b]/f_max)*g_height;
+
+		cairo_new_path(cr);
+		cairo_rectangle(cr, (g_width/nbins)*b, g_height,
+		                    g_width/nbins, -bar_height);
+		cairo_set_source_rgb(cr, 0.0, 0.0, 1.0);
+		cairo_set_line_width(cr, 1.0);
+		cairo_stroke(cr);
+
+	}
+
+	cairo_new_path(cr);
+	cairo_rectangle(cr, 0.0, 0.0, g_width, g_height);
+	cairo_set_source_rgb(cr, 0.0, 0.0, 0.0);
+	cairo_set_line_width(cr, 1.5);
+	cairo_stroke(cr);
+}
+
+
+static void watermark(struct _srcontext *sr)
+{
+	show_text(sr->cr, "Written by partialator from CrystFEL"
+	                 " version "PACKAGE_VERSION, sr->h-15.0, J_RIGHT,
+	                 "Sans 7");
+}
+
+
+SRContext *sr_header(const char *filename, const char *stream_filename,
+                     const char *cmdline)
+{
+	char tmp[1024];
+	struct _srcontext *sr;
+
+	sr = malloc(sizeof(*sr));
+	if ( sr == NULL ) return NULL;
+
+	sr->w = PAGE_WIDTH;
+	sr->h = 595.0;
+
+	sr->surf = cairo_pdf_surface_create(filename, sr->w, sr->h);
+
+	if ( cairo_surface_status(sr->surf) != CAIRO_STATUS_SUCCESS ) {
 		fprintf(stderr, "Couldn't create Cairo surface\n");
-		cairo_surface_destroy(surface);
-		return;
+		cairo_surface_destroy(sr->surf);
+		free(sr);
+		return NULL;
 	}
 
-	cr = cairo_create(surface);
+	sr->cr = cairo_create(sr->surf);
+
+	snprintf(tmp, 1023, "Scaling report: %s", stream_filename);
+	show_text(sr->cr, tmp, 10.0, J_CENTER, "Sans Bold 16");
+	snprintf(tmp, 1023, "partialator %s", cmdline);
+	show_text(sr->cr, tmp, 45.0, J_LEFT, "Mono 7");
+	watermark(sr);
+
+	return sr;
+}
+
+
+void sr_before(SRContext *sr, struct image *images, int n, RefList *full)
+{
+	if ( sr == NULL ) return;
+
+	cairo_save(sr->cr);
+	cairo_translate(sr->cr, 75.0, 100.0);
+	scale_factor_histogram(sr->cr, images, n, "Before refinement");
+	cairo_translate(sr->cr, 60.0, 235.0);
+	partiality_graph(sr->cr, images, n, full);
+	cairo_restore(sr->cr);
+}
+
+
+void sr_after(SRContext *sr, struct image *images, int n, RefList *full)
+{
+	if ( sr == NULL ) return;
+
+	cairo_save(sr->cr);
+	cairo_translate(sr->cr, 475.0, 100.0);
+	scale_factor_histogram(sr->cr, images, n, "After refinement");
+	cairo_translate(sr->cr, 60.0, 235.0);
+	partiality_graph(sr->cr, images, n, full);
+	cairo_restore(sr->cr);
+
+	cairo_surface_show_page(sr->surf);
+	watermark(sr);
+
 
-	write_title(cr, stream_filename, w, h);
+	cairo_surface_finish(sr->surf);
+	cairo_destroy(sr->cr);
 
-	cairo_surface_finish(surface);
-	cairo_destroy(cr);
+	free(sr);
 }
diff --git a/src/scaling-report.h b/src/scaling-report.h
index 3186c5f7..063a0f39 100644
--- a/src/scaling-report.h
+++ b/src/scaling-report.h
@@ -19,16 +19,35 @@
 
 #include "utils.h"
 
+typedef struct _srcontext SRContext;  /* Opaque */
+
 #ifdef HAVE_CAIRO
-extern void scaling_report(const char *filename, const struct image *images,
-                           int n, const char *stream_filename);
+
+extern SRContext *sr_header(const char *filename, const char *stream_filename,
+                            const char *cmdline);
+
+extern void sr_before(SRContext *sr, struct image *images, int n,
+                      RefList *full);
+
+extern void sr_after(SRContext *sr, struct image *images, int n,
+                     RefList *full);
+
 #else
-static inline void scaling_report(const char *filename,
-                                  const struct image *images, int n,
-                                  const char *stream_filename)
+
+SRContext *sr_header(const char *filename, const char *stream_filename,
+                     const char *cmdline)
 {
-	ERROR("Not writing scaling report - no Cairo support.\n");
+	return NULL;
 }
+
+void sr_before(SRContext *sr, struct image *images, int n, RefList *full)
+{
+}
+
+void sr_after(SRContext *sr, struct image *images, int n, RefList *full)
+{
+}
+
 #endif
 
 #endif	/* SCALING_REPORT_H */