Add unit test for peak integration

2 files changed, 144 insertions, 0 deletions

diff --git a/tests/.gitignore b/tests/.gitignore
index b855d2bd..938b5a4a 100644
--- a/tests/.gitignore
+++ b/tests/.gitignore
@@ -2,4 +2,5 @@
 .deps
 list_check
 gpu_sim_check
+integration_check
 .dirstamp
diff --git a/tests/integration_check.c b/tests/integration_check.c
new file mode 100644
index 00000000..f80181ee
--- /dev/null
+++ b/tests/integration_check.c
@@ -0,0 +1,143 @@
+/*
+ * integration_check.c
+ *
+ * Check peak integration
+ *
+ * (c) 2011 Thomas White <taw@physics.org>
+ *
+ * Part of CrystFEL - crystallography with a FEL
+ *
+ */
+
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+
+#include <stdlib.h>
+#include <stdio.h>
+
+#include "../src/image.h"
+#include "../src/peaks.h"
+#include "../src/utils.h"
+#include "../src/beam-parameters.h"
+
+
+int main(int argc, char *argv[])
+{
+	struct image image;
+	double fsp, ssp, intensity, bg, max;
+	int fs, ss;
+
+	image.data = malloc(128*128*sizeof(float));
+	image.flags = NULL;
+
+	image.det = calloc(1, sizeof(struct detector));
+	image.det->n_panels = 1;
+	image.det->panels = calloc(1, sizeof(struct panel));
+
+	image.det->panels[0].min_fs = 0;
+	image.det->panels[0].max_fs = 128;
+	image.det->panels[0].min_ss = 0;
+	image.det->panels[0].max_ss = 128;
+	image.det->panels[0].fsx = 1.0;
+	image.det->panels[0].fsy = 0.0;
+	image.det->panels[0].ssx = 0.0;
+	image.det->panels[0].ssy = 1.0;
+	image.det->panels[0].xfs = 1.0;
+	image.det->panels[0].yfs = 0.0;
+	image.det->panels[0].xss = 0.0;
+	image.det->panels[0].yss = 1.0;
+	image.det->panels[0].cnx = -64.0;
+	image.det->panels[0].cny = -64.0;
+	image.det->panels[0].clen = 1.0;
+	image.det->panels[0].res = 1.0;
+	image.det->panels[0].integr_radius = 10.0;
+
+	image.width = 128;
+	image.height = 128;
+	memset(image.data, 0, 128*128*sizeof(float));
+
+	image.beam = calloc(1, sizeof(struct beam_params));
+	image.beam->adu_per_photon = 100.0;
+
+	/* First check: no intensity -> zero intensity and bg */
+	integrate_peak(&image, 64, 64, &fsp, &ssp, &intensity, &bg, &max, 0, 1);
+	STATUS(" First check: intensity = %.2f bg = %.2f max = %.2f\n",
+	       intensity, bg, max);
+	if ( intensity != 0.0 ) {
+		ERROR("Intensity should be zero.\n");
+		return 1;
+	}
+	if ( bg != 0.0 ) {
+		ERROR("Background should be zero\n");
+		return 1;
+	}
+
+	/* Second check: uniform peak gives correct value and no bg */
+	for ( fs=0; fs<image.width; fs++ ) {
+	for ( ss=0; ss<image.height; ss++ ) {
+		if ( (fs-64)*(fs-64) + (ss-64)*(ss-64) > 9*9 ) continue;
+		image.data[fs+image.width*ss] = 1000.0;
+	}
+	}
+	integrate_peak(&image, 64, 64, &fsp, &ssp, &intensity, &bg, &max, 0, 1);
+	STATUS("Second check: intensity = %.2f, bg = %.2f, max = %.2f\n",
+	       intensity, bg, max);
+	if ( fabs(intensity - M_PI*9.0*9.0*1000.0) > 4000.0 ) {
+		ERROR("Intensity should be close to 1000*pi*integr_r^2\n");
+		return 1;
+	}
+	if ( bg != 0.0 ) {
+		ERROR("Background should be zero\n");
+		return 1;
+	}
+	if ( max != 1000.0 ) {
+		ERROR("Max should be 1000\n");
+		return 1;
+	}
+
+	/* Third check: flat Poisson background should hoover up bg only */
+	for ( fs=0; fs<image.width; fs++ ) {
+	for ( ss=0; ss<image.height; ss++ ) {
+		image.data[fs+image.width*ss] = poisson_noise(10.0) - 10.0;
+	}
+	}
+	integrate_peak(&image, 64, 64, &fsp, &ssp, &intensity, &bg, &max, 0, 1);
+	STATUS(" Third check: intensity = %.2f, bg = %.2f, max = %.2f\n",
+	       intensity, bg, max);
+	if ( fabs(intensity) > 100.0 ) {
+		ERROR("Intensity should be close to zero\n");
+		return 1;
+	}
+	if ( fabs(bg - sqrt(10.0)) > 1.0 ) {
+		ERROR("Background should be close to sqrt(10)\n");
+		return 1;
+	}
+
+	for ( fs=0; fs<image.width; fs++ ) {
+	for ( ss=0; ss<image.height; ss++ ) {
+		if ( (fs-64)*(fs-64) + (ss-64)*(ss-64) > 9*9 ) continue;
+		image.data[fs+image.width*ss] = 1000.0;
+	}
+	}
+	integrate_peak(&image, 64, 64, &fsp, &ssp, &intensity, &bg, &max, 0, 1);
+	STATUS("Fourth check: intensity = %.2f, bg = %.2f, max = %.2f\n",
+	       intensity, bg, max);
+	if ( fabs(intensity - M_PI*9.0*9.0*1000.0) > 4000.0 ) {
+		ERROR("Intensity should be close to 1000*pi*peak_r^2\n");
+		return 1;
+	}
+	if ( fabs(bg - sqrt(10.0)) > 1.0 ) {
+		ERROR("Background should be close to sqrt(10)\n");
+		return 1;
+	}
+
+	free(image.beam);
+	free(image.det->panels);
+	free(image.det);
+	free(image.data);
+
+	return 0;
+}