Merge branch 'master' of git.bitwiz.org.uk:crystfel

13 files changed, 358 insertions, 43 deletions

diff --git a/scripts/double-hit b/scripts/double-hit
index 831e2df6..55a95ee1 100755
--- a/scripts/double-hit
+++ b/scripts/double-hit
@@ -10,6 +10,16 @@ my $line;
 my $filename;
 my $full_filename;
 
+my $hsteps = 128;
+my @hist;
+my $hmin = 0;
+my $hmax = 5;
+my $hstep = ($hmax - $hmin)/$hsteps;
+
+for ( my $i=0; $i<$hsteps; $i++ ) {
+	$hist[$i] = 0;
+}
+
 while ( $line = <FH> ) {
 
 	if ( $line =~ /^Reflections\ from\ indexing\ in\ (.+)$/ ) {
@@ -50,6 +60,9 @@ while ( $line = <FH> ) {
 		if ( $div == 1.0 ) {
 			printf("Single hit: %s\n", $full_filename);
 		}
+
+		my $bin = int(($div-$hmin)/$hstep);
+		$hist[$bin]++;
 	}
 
 }
@@ -58,6 +71,7 @@ close(FH);
 close(OFH);
 
 open(GP, "| gnuplot");
+
 print(GP "set term postscript enhanced font \"Helvetica,20\"\n");
 print(GP "set output \"doublehit.ps\"\n");
 print(GP "set xtics nomirror out rotate by -60\n");
@@ -66,6 +80,19 @@ print(GP "set timefmt \"%Y/%b/%d-%H:%M:%S\"\n");
 print(GP "set format x \"%d/%b %H:%M\"\n");
 print(GP "unset key\n");
 print(GP "plot [] [0:5] \"doublehit.dat\" u 1:2 w points\n");
+
+open(OFH, "> doublehit-hist.dat");
+for ( my $i=0; $i<$hsteps; $i++ ) {
+	printf(OFH "%f %f\n", $hmin+$hstep*$i+($hstep/2), $hist[$i]);
+}
+close(OFH);
+
+print(GP "set xtics nomirror out rotate by 0\n");
+print(GP "unset xdata\n");
+print(GP "set format x \"% g\"\n");
+print(GP "unset key\n");
+print(GP "plot [] [0:1000] \"doublehit-hist.dat\" u 1:2 w histeps\n");
+
 close(GP);
 
 system("ps2pdf doublehit.ps");
diff --git a/scripts/i0-analysis b/scripts/i0-analysis
new file mode 100755
index 00000000..1c797139
--- /dev/null
+++ b/scripts/i0-analysis
@@ -0,0 +1,106 @@
+#!/usr/bin/perl -w
+
+use strict;
+use File::Basename;
+
+open(FH, $ARGV[0]);
+
+my $line;
+my $filename;
+my $full_filename;
+
+my $hsteps = 128;
+my @hist;
+
+open(GP, "| gnuplot");
+print(GP "set term postscript enhanced font \"Helvetica,20\"\n");
+print(GP "set output \"i0.ps\"\n");
+print(GP "unset key\n");
+print(GP "set xtics nomirror out rotate by 0\n");
+print(GP "unset xdata\n");
+print(GP "set format x \"% g\"\n");
+print(GP "unset key\n");
+
+
+# First graph: distribution of scaling factors
+my $hmin = -0.1;
+my $hmax = 1.5;
+my $hstep = ($hmax - $hmin)/$hsteps;
+for ( my $i=0; $i<$hsteps; $i++ ) {
+	$hist[$i] = 0;
+}
+while ( $line = <FH> ) {
+	if ( $line =~ /^([\d\.]+)\ ([\d\.]+)$/ ) {
+		my $scale = $1;
+		my $i0 = $2;
+		my $bin = int(($scale-$hmin)/$hstep);
+		$hist[$bin]++;
+	}
+}
+close(FH);
+open(OFH, "> i0-scale-hist.dat");
+for ( my $i=0; $i<$hsteps; $i++ ) {
+	printf(OFH "%f %f\n", $hmin+$hstep*$i+($hstep/2), $hist[$i]);
+}
+close(OFH);
+print(GP "set title \"Distribution of Scaling Factors\"\n");
+print(GP "plot [] [] \"i0-scale-hist.dat\" u 1:2 w histeps\n");
+
+
+# Second graph: distribution of gas detector values
+$hmin = -0.5;
+$hmax = 3.0;
+$hstep = ($hmax - $hmin)/$hsteps;
+open(FH, $ARGV[0]);
+for ( my $i=0; $i<$hsteps; $i++ ) {
+	$hist[$i] = 0;
+}
+while ( $line = <FH> ) {
+	if ( $line =~ /^([\d\.]+)\ ([\d\.]+)$/ ) {
+		my $scale = $1;
+		my $i0 = $2;
+		my $bin = int(($i0-$hmin)/$hstep);
+		$hist[$bin]++;
+	}
+}
+close(FH);
+open(OFH, "> i0-hist.dat");
+for ( my $i=0; $i<$hsteps; $i++ ) {
+	printf(OFH "%f %f\n", $hmin+$hstep*$i+($hstep/2), $hist[$i]);
+}
+close(OFH);
+print(GP "set title \"Distribution of Gas Detector Values\"\n");
+print(GP "plot [] [] \"i0-hist.dat\" u 1:2 w histeps\n");
+
+
+# Third graph: distribution of scale factors after scaling i0=1
+$hmin = -0.1;
+$hmax = 1.5;
+$hstep = ($hmax - $hmin)/$hsteps;
+open(FH, $ARGV[0]);
+for ( my $i=0; $i<$hsteps; $i++ ) {
+	$hist[$i] = 0;
+}
+while ( $line = <FH> ) {
+	if ( $line =~ /^([\d\.]+)\ ([\d\.]+)$/ ) {
+		my $scale = $1;
+		my $i0 = $2;
+		my $bin = int((($scale/$i0)-$hmin)/$hstep);
+		$hist[$bin]++;
+	}
+}
+close(FH);
+open(OFH, "> i0-prior-hist.dat");
+for ( my $i=0; $i<$hsteps; $i++ ) {
+	printf(OFH "%f %f\n", $hmin+$hstep*$i+($hstep/2), $hist[$i]);
+}
+close(OFH);
+print(GP "set title \"Distribution of Scaling Factors with I0=1\"\n");
+print(GP "plot [] [] \"i0-prior-hist.dat\" u 1:2 w histeps\n");
+
+
+close(GP);
+
+system("ps2pdf i0.ps");
+unlink("i0.dat");
+unlink("i0.ps");
diff --git a/src/diffraction.c b/src/diffraction.c
index 1960f786..b3e80e8f 100644
--- a/src/diffraction.c
+++ b/src/diffraction.c
@@ -240,8 +240,8 @@ struct rvec get_q(struct image *image, unsigned int xs, unsigned int ys,
 {
 	struct rvec q;
 	float twotheta, r, az;
-	float rx = 0.0;
-	float ry = 0.0;
+	float rx;
+	float ry;
 	struct panel *p;
 
 	const unsigned int x = xs / sampling;
@@ -267,6 +267,26 @@ struct rvec get_q(struct image *image, unsigned int xs, unsigned int ys,
 }
 
 
+double get_tt(struct image *image, unsigned int xs, unsigned int ys)
+{
+	float r, rx, ry;
+	struct panel *p;
+
+	const unsigned int x = xs;
+	const unsigned int y = ys; /* Integer part only */
+
+	p = find_panel(&image->det, x, y);
+
+	rx = ((float)xs - p->cx) / p->res;
+	ry = ((float)ys - p->cy) / p->res;
+
+	/* Calculate q-vector for this sub-pixel */
+	r = sqrt(pow(rx, 2.0) + pow(ry, 2.0));
+
+	return atan2(r, p->clen);
+}
+
+
 void get_diffraction(struct image *image, int na, int nb, int nc,
                      const double *intensities, const unsigned int *counts,
                      UnitCell *cell, int do_water, GradientMethod m)
diff --git a/src/diffraction.h b/src/diffraction.h
index 85d06cc2..ed583feb 100644
--- a/src/diffraction.h
+++ b/src/diffraction.h
@@ -32,4 +32,6 @@ extern void get_diffraction(struct image *image, int na, int nb, int nc,
 extern struct rvec get_q(struct image *image, unsigned int xs, unsigned int ys,
                          unsigned int sampling, float *ttp, float k);
 
+extern double get_tt(struct image *image, unsigned int xs, unsigned int ys);
+
 #endif	/* DIFFRACTION_H */
diff --git a/src/hdf5-file.c b/src/hdf5-file.c
index d76f2bd4..5573267b 100644
--- a/src/hdf5-file.c
+++ b/src/hdf5-file.c
@@ -204,6 +204,24 @@ static double get_wavelength(struct hdfile *f)
 }
 
 
+static double get_f0(struct hdfile *f)
+{
+	herr_t r;
+	hid_t dh;
+	double f0;
+
+	dh = H5Dopen(f->fh, "/LCLS/f_11_ENRC", H5P_DEFAULT);
+	if ( dh < 0 ) return -1.0;
+
+	r = H5Dread(dh, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL,
+		            H5P_DEFAULT, &f0);
+	H5Dclose(dh);
+	if ( r < 0 ) return -1.0;
+
+	return f0;
+}
+
+
 static void debodge_saturation(struct hdfile *f, struct image *image)
 {
 	hid_t dh, sh;
@@ -321,6 +339,15 @@ int hdf5_read(struct hdfile *f, struct image *image)
 		image->lambda = ph_en_to_lambda(eV_to_J(2000.0));
 	}
 
+	image->f0 = get_f0(f);
+	if ( image->f0 < 0.0 ) {
+		ERROR("Couldn't read incident intensity - using 1.0.\n");
+		image->f0 = 1.0;
+		image->f0_available = 0;
+	} else {
+		image->f0_available = 1;
+	}
+
 	debodge_saturation(f, image);
 
 	return 0;
diff --git a/src/image.h b/src/image.h
index 4ba9e5f9..db1b491e 100644
--- a/src/image.h
+++ b/src/image.h
@@ -67,6 +67,16 @@ struct rvec
 };
 
 
+struct reflhit {
+	signed int h;
+	signed int k;
+	signed int l;
+	double min_distance;
+	int x;
+	int y;
+};
+
+
 /* Structure describing an image */
 struct image {
 
@@ -77,6 +87,8 @@ struct image {
 	int			ncells;
 	struct detector		det;
 	const char		*filename;
+	struct reflhit		*hits;
+	int			n_hits;
 
 	int			id;
 
@@ -85,6 +97,10 @@ struct image {
 	/* Wavelength must always be given */
 	double			lambda;		/* Wavelength in m */
 
+	/* Incident intensity (if unknown, put 1.0) */
+	double			f0;
+	int			f0_available;
+
 	int			width;
 	int			height;
 
diff --git a/src/indexamajig.c b/src/indexamajig.c
index bd3b53e2..6cf9cd15 100644
--- a/src/indexamajig.c
+++ b/src/indexamajig.c
@@ -57,6 +57,7 @@ struct process_args
 	int config_gpu;
 	int config_simulate;
 	int config_nomatch;
+	int config_unpolar;
 	IndexingMethod indm;
 	const double *intensities;
 	const unsigned int *counts;
@@ -66,6 +67,7 @@ struct process_args
 struct process_result
 {
 	int hit;
+	int peaks_sane;
 };
 
 
@@ -118,6 +120,7 @@ static void show_help(const char *s)
 "      --no-match          Don't attempt to match the indexed cell to the\n"
 "                           model, just proceed with the one generated by the\n"
 "                           auto-indexing procedure.\n"
+"      --unpolarized       Don't correct for the polarisation of the X-rays.\n"
 "\n\nOptions for greater performance or verbosity:\n\n"
 "      --verbose           Be verbose about indexing.\n"
 "      --gpu               Use the GPU to speed up the simulation.\n"
@@ -190,6 +193,11 @@ static struct image *get_simage(struct image *template, int alternate)
 	image->features = template->features;
 	image->filename = template->filename;
 	image->indexed_cell = template->indexed_cell;
+	image->f0 = template->f0;
+
+	/* Prevent muppetry */
+	image->hits = NULL;
+	image->n_hits = 0;
 
 	return image;
 }
@@ -237,6 +245,7 @@ static void *process_image(void *pargsv)
 	int config_gpu = pargs->config_gpu;
 	int config_simulate = pargs->config_simulate;
 	int config_nomatch = pargs->config_nomatch;
+	int config_unpolar = pargs->config_unpolar;
 	IndexingMethod indm = pargs->indm;
 	const double *intensities = pargs->intensities;
 	const unsigned int *counts = pargs->counts;
@@ -248,6 +257,8 @@ static void *process_image(void *pargsv)
 	image.indexed_cell = NULL;
 	image.id = pargs->id;
 	image.filename = filename;
+	image.hits = NULL;
+	image.n_hits = 0;
 
 	STATUS("Processing '%s'\n", image.filename);
 
@@ -295,6 +306,12 @@ static void *process_image(void *pargsv)
 
 	/* Perform 'fine' peak search */
 	search_peaks(&image);
+
+	/* Get rid of noise-filtered version at this point
+	 * - it was strictly for the purposes of peak detection. */
+	free(image.data);
+	image.data = data_for_measurement;
+
 	if ( image_feature_count(image.features) < 5 ) goto done;
 
 	if ( config_dumpfound ) dump_peaks(&image, pargs->output_mutex);
@@ -312,19 +329,26 @@ static void *process_image(void *pargsv)
 	}
 
 	/* No cell at this point?  Then we're done. */
+	result->peaks_sane = 0;
 	if ( image.indexed_cell == NULL ) goto done;
 
-	simage = get_simage(&image, config_alternate);
+	/* Sanity check */
+	if ( !peak_sanity_check(&image, image.indexed_cell) ) {
+		STATUS("Failed peak sanity check.\n");
+		result->peaks_sane = 0;
+		goto done;
+	} else {
+		result->peaks_sane = 1;
+	}
 
 	/* Measure intensities if requested */
 	if ( config_nearbragg ) {
-		/* Use original data (temporarily) */
-		simage->data = data_for_measurement;
-		output_intensities(simage, image.indexed_cell,
-		                   pargs->output_mutex);
-		simage->data = NULL;
+		output_intensities(&image, image.indexed_cell,
+		                   pargs->output_mutex, config_unpolar);
 	}
 
+	simage = get_simage(&image, config_alternate);
+
 	/* Simulate if requested */
 	if ( config_simulate ) {
 		if ( config_gpu ) {
@@ -350,7 +374,7 @@ done:
 	free(image.data);
 	free(image.det.panels);
 	image_feature_list_free(image.features);
-	free(data_for_measurement);
+	free(image.hits);
 	hdfile_close(hdfile);
 
 	if ( image.indexed_cell == NULL ) {
@@ -371,6 +395,7 @@ int main(int argc, char *argv[])
 	char *rval = NULL;
 	int n_images;
 	int n_hits;
+	int n_sane;
 	int config_noindex = 0;
 	int config_dumpfound = 0;
 	int config_nearbragg = 0;
@@ -382,6 +407,7 @@ int main(int argc, char *argv[])
 	int config_gpu = 0;
 	int config_verbose = 0;
 	int config_alternate = 0;
+	int config_unpolar = 0;
 	IndexingMethod indm;
 	char *indm_str = NULL;
 	UnitCell *cell;
@@ -417,6 +443,7 @@ int main(int argc, char *argv[])
 		{"intensities",        1, NULL,               'q'},
 		{"pdb",                1, NULL,               'p'},
 		{"prefix",             1, NULL,               'x'},
+		{"unpolarized",        0, &config_unpolar,     1},
 		{0, 0, NULL, 0}
 	};
 
@@ -535,6 +562,7 @@ int main(int argc, char *argv[])
 	gsl_set_error_handler_off();
 	n_images = 0;
 	n_hits = 0;
+	n_sane = 0;
 
 	for ( i=0; i<nthreads; i++ ) {
 		worker_args[i] = malloc(sizeof(struct process_args));
@@ -570,6 +598,7 @@ int main(int argc, char *argv[])
 		pargs->config_gpu = config_gpu;
 		pargs->config_simulate = config_simulate;
 		pargs->config_nomatch = config_nomatch;
+		pargs->config_unpolar = config_unpolar;
 		pargs->cell = cell;
 		pargs->indm = indm;
 		pargs->intensities = intensities;
@@ -616,6 +645,7 @@ int main(int argc, char *argv[])
 
 			if ( result != NULL ) {
 				n_hits += result->hit;
+				n_sane += result->peaks_sane;
 				free(result);
 			}
 
@@ -666,7 +696,7 @@ int main(int argc, char *argv[])
 	fclose(fh);
 
 	STATUS("There were %i images.\n", n_images);
-	STATUS("%i hits were found.\n", n_hits);
+	STATUS("%i hits were found, of which %i were sane.\n", n_hits, n_sane);
 
 	if ( gctx != NULL ) {
 		cleanup_gpu(gctx);
diff --git a/src/likelihood.c b/src/likelihood.c
index e8c62010..95f005c7 100644
--- a/src/likelihood.c
+++ b/src/likelihood.c
@@ -26,13 +26,13 @@ void detwin_intensities(const double *model, double *new_pattern,
 
 void scale_intensities(const double *model, double *new_pattern,
                        const unsigned int *model_counts,
-                       unsigned int *new_counts)
+                       unsigned int *new_counts, double f0, int f0_valid)
 {
 	double s;
 	unsigned int i;
 
 	s = stat_scale_intensity(model, model_counts, new_pattern, new_counts);
-	printf("%f\n", s);
+	if ( f0_valid ) printf("%f %f\n", s, f0);
 
 	/* NaN -> abort */
 	if ( isnan(s) ) return;
diff --git a/src/likelihood.h b/src/likelihood.h
index 9a773829..0a2bb407 100644
--- a/src/likelihood.h
+++ b/src/likelihood.h
@@ -23,7 +23,8 @@ extern void detwin_intensities(const double *model, double *new_pattern,
 
 extern void scale_intensities(const double *model, double *new_pattern,
 	                      const unsigned int *model_counts,
-                              unsigned int *new_counts);
+                              unsigned int *new_counts, double f0,
+                              int f0_valid);
 
 
 #endif	/* LIKELIHOOD_H */
diff --git a/src/pattern_sim.c b/src/pattern_sim.c
index a0d81a9d..e09ab72a 100644
--- a/src/pattern_sim.c
+++ b/src/pattern_sim.c
@@ -346,7 +346,7 @@ int main(int argc, char *argv[])
 		record_image(&image, !config_nonoise);
 
 		if ( config_nearbragg ) {
-			output_intensities(&image, cell, NULL);
+			output_intensities(&image, cell, NULL, 1);
 		}
 
 		if ( config_powder ) {
diff --git a/src/peaks.c b/src/peaks.c
index cf4eacf3..22e12949 100644
--- a/src/peaks.c
+++ b/src/peaks.c
@@ -34,21 +34,15 @@
 #define MAX_HITS (1024)
 
 /* How close a peak must be to an indexed position to be considered "close"
- * for the purposes of double hit detection etc. */
+ * for the purposes of double hit detection and sanity checking. */
 #define PEAK_CLOSE (30.0)
 
 /* How close a peak must be to an indexed position to be considered "close"
  * for the purposes of integration. */
 #define PEAK_REALLY_CLOSE (10.0)
 
-struct reflhit {
-	signed int h;
-	signed int k;
-	signed int l;
-	double min_distance;
-	int x;
-	int y;
-};
+/* Degree of polarisation of X-ray beam */
+#define POL (1.0)
 
 
 #define PEAK_WINDOW_SIZE (10)
@@ -140,7 +134,8 @@ static void cull_peaks(struct image *image)
 
 
 static void integrate_peak(struct image *image, int xp, int yp,
-                           float *xc, float *yc, float *intensity)
+                           float *xc, float *yc, float *intensity,
+                           int do_polar)
 {
 	signed int x, y;
 	const int lim = INTEGRATION_RADIUS * INTEGRATION_RADIUS;
@@ -154,6 +149,7 @@ static void integrate_peak(struct image *image, int xp, int yp,
 		struct panel *p;
 		double val, sa, pix_area, Lsq, dsq, proj_area;
 		float tt;
+		double phi, pa, pb, pol;
 
 		/* Circular mask */
 		if ( x*x + y*y > lim ) continue;
@@ -168,7 +164,7 @@ static void integrate_peak(struct image *image, int xp, int yp,
 		Lsq = pow(p->clen, 2.0);
 
 		/* Area of pixel as seen from crystal (approximate) */
-		get_q(image, x+xp, y+yp, 1, &tt, 1.0 / image->lambda);
+		tt = get_tt(image, x+xp, y+yp);
 		proj_area = pix_area * cos(tt);
 
 		/* Calculate distance from crystal to pixel */
@@ -178,7 +174,16 @@ static void integrate_peak(struct image *image, int xp, int yp,
 		/* Projected area of pixel divided by distance squared */
 		sa = 1.0e7 * proj_area / (dsq + Lsq);
 
-		val = image->data[(x+xp)+image->width*(y+yp)] / sa;
+		if ( do_polar ) {
+			phi = atan2(y+yp, x+xp);
+			pa = pow(sin(phi)*sin(tt), 2.0);
+			pb = pow(cos(tt), 2.0);
+			pol = 1.0 - 2.0*POL*(1-pa) + POL*(1.0+pb);
+		} else {
+			pol = 1.0;
+		}
+
+		val = image->data[(x+xp)+image->width*(y+yp)] / (sa*pol);
 
 		total += val;
 
@@ -305,8 +310,10 @@ void search_peaks(struct image *image)
 			continue;
 		}
 
-		/* Centroid peak and get better coordinates */
-		integrate_peak(image, mask_x, mask_y, &fx, &fy, &intensity);
+		/* Centroid peak and get better coordinates.
+		 * Don't bother doing polarisation correction, because the
+		 * intensity of this peak is only an estimate at this stage. */
+		integrate_peak(image, mask_x, mask_y, &fx, &fy, &intensity, 0);
 
 		/* It is possible for the centroid to fall outside the image */
 		if ( (fx < 0.0) || (fx > image->width)
@@ -369,23 +376,17 @@ void dump_peaks(struct image *image, pthread_mutex_t *mutex)
 }
 
 
-void output_intensities(struct image *image, UnitCell *cell,
-                        pthread_mutex_t *mutex)
+static int find_projected_peaks(struct image *image, UnitCell *cell)
 {
 	int x, y;
 	double ax, ay, az;
 	double bx, by, bz;
 	double cx, cy, cz;
-	double a, b, c, al, be, ga;
-	double asx, asy, asz;
-	double bsx, bsy, bsz;
-	double csx, csy, csz;
-	struct reflhit hits[MAX_HITS];
+	struct reflhit *hits;
 	int n_hits = 0;
-	int i;
-	int n_found;
-	int n_indclose = 0;
-	int n_foundclose = 0;
+
+	hits = malloc(sizeof(struct reflhit)*MAX_HITS);
+	if ( hits == NULL ) return 0;
 
 	cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
 
@@ -446,6 +447,59 @@ void output_intensities(struct image *image, UnitCell *cell,
 	}
 
 	STATUS("Found %i reflections\n", n_hits);
+	image->hits = hits;
+	image->n_hits = n_hits;
+
+	return n_hits;
+}
+
+
+int peak_sanity_check(struct image *image, UnitCell *cell)
+{
+	int i;
+	int n_sane = 0;
+
+	find_projected_peaks(image, cell);
+	if ( image->n_hits == 0 ) return 0;  /* Failed sanity check: no peaks */
+
+	for ( i=0; i<image->n_hits; i++ ) {
+
+		double d;
+		int idx;
+		struct imagefeature *f;
+
+		f = image_feature_closest(image->features,
+                                          image->hits[i].x, image->hits[i].y,
+		                          &d, &idx);
+		if ( (f != NULL) && (d < PEAK_CLOSE) ) {
+			n_sane++;
+		}
+
+	}
+
+	STATUS("Sanity factor: %f / %f = %f\n", (float)n_sane, (float)image->n_hits,
+                                                (float)n_sane / (float)image->n_hits);
+	if ( (float)n_sane / (float)image->n_hits < 0.8 ) return 0;
+
+	return 1;
+}
+
+
+void output_intensities(struct image *image, UnitCell *cell,
+                        pthread_mutex_t *mutex, int unpolar)
+{
+	int i;
+	int n_found;
+	int n_indclose = 0;
+	int n_foundclose = 0;
+	double asx, asy, asz;
+	double bsx, bsy, bsz;
+	double csx, csy, csz;
+	double a, b, c, al, be, ga;
+	int n_hits = image->n_hits;
+	struct reflhit *hits = image->hits;
+
+	if ( image->n_hits == 0 ) return;
 
 	/* Get exclusive access to the output stream if necessary */
 	if ( mutex != NULL ) pthread_mutex_lock(mutex);
@@ -470,6 +524,13 @@ void output_intensities(struct image *image, UnitCell *cell,
 	printf("cstar = %+9.7f %+9.7f %+9.7f nm^-1\n",
 	       csx/1e9, csy/1e9, csz/1e9);
 
+	if ( image->f0_available ) {
+		printf("f0 = %7.5f (arbitrary gas detector units)\n",
+		       image->f0);
+	} else {
+		printf("f0 = invalid\n");
+	}
+
 	for ( i=0; i<n_hits; i++ ) {
 
 		float x, y, intensity;
@@ -485,13 +546,14 @@ void output_intensities(struct image *image, UnitCell *cell,
 			/* f->intensity was measured on the filtered pattern,
 			 * so instead re-integrate using old coordinates.
 			 * This will produce further revised coordinates. */
-			integrate_peak(image, f->x, f->y, &x, &y, &intensity);
+			integrate_peak(image, f->x, f->y, &x, &y, &intensity,
+			               !unpolar);
 			intensity = f->intensity;
 
 		} else {
 
 			integrate_peak(image, hits[i].x, hits[i].y,
-			               &x, &y, &intensity);
+			               &x, &y, &intensity, !unpolar);
 
 		}
 
diff --git a/src/peaks.h b/src/peaks.h
index a5c771b7..10cf5b0a 100644
--- a/src/peaks.h
+++ b/src/peaks.h
@@ -22,6 +22,7 @@
 extern void search_peaks(struct image *image);
 extern void dump_peaks(struct image *image, pthread_mutex_t *mutex);
 extern void output_intensities(struct image *image, UnitCell *cell,
-                               pthread_mutex_t *mutex);
+                               pthread_mutex_t *mutex, int unpolar);
+extern int peak_sanity_check(struct image *image, UnitCell *cell);
 
 #endif	/* PEAKS_H */
diff --git a/src/process_hkl.c b/src/process_hkl.c
index 53638e7b..7c42cd3e 100644
--- a/src/process_hkl.c
+++ b/src/process_hkl.c
@@ -245,6 +245,8 @@ int main(int argc, char *argv[])
 	unsigned int *new_counts = NULL;
 	unsigned int n_total_patterns;
 	unsigned int *truecounts = NULL;
+	float f0;
+	int f0_valid;
 
 	/* Long options */
 	const struct option longopts[] = {
@@ -359,6 +361,7 @@ int main(int argc, char *argv[])
 	STATUS("There are %i patterns to process\n", n_total_patterns);
 
 	n_patterns = 0;
+	f0_valid = 0;
 	do {
 
 		char line[1024];
@@ -382,10 +385,17 @@ int main(int argc, char *argv[])
 				                   model_counts, new_counts);
 			}
 
+			/* Assume a default I0 if we don't have one by now */
+			if ( config_scale && !f0_valid ) {
+				ERROR("No f0 value.\n");
+				f0 = 1.0;
+			}
+
 			/* Scale if requested */
 			if ( config_scale ) {
 				scale_intensities(model, new_pattern,
-				                  model_counts, new_counts);
+				                  model_counts, new_counts, f0,
+				                  f0_valid);
 			}
 
 			/* Start of second or later pattern */
@@ -409,6 +419,19 @@ int main(int argc, char *argv[])
 
 			progress_bar(n_patterns, n_total_patterns, "Merging");
 
+			f0_valid = 0;
+
+		}
+
+		if ( strncmp(line, "f0 = ", 5) == 0 ) {
+			r = sscanf(line, "f0 = %f", &f0);
+			if ( r != 1 ) {
+				ERROR("Couldn't understand f0 line.\n");
+				f0 = 1.0;
+				f0_valid = 0;
+				continue;
+			}
+			f0_valid = 1;
 		}
 
 		r = sscanf(line, "%i %i %i %f", &h, &k, &l, &intensity);