Add unit test for peak integration

4 files changed, 180 insertions, 6 deletions

diff --git a/Makefile.am b/Makefile.am
index db6d1286..3d6949ed 100644
--- a/Makefile.am
+++ b/Makefile.am
@@ -7,10 +7,11 @@ bin_PROGRAMS = src/pattern_sim src/process_hkl src/get_hkl src/indexamajig \
                src/calibrate_detector src/partialator \
                src/check_hkl src/sum_stack
 
-noinst_PROGRAMS = tests/list_check
+noinst_PROGRAMS = tests/list_check tests/integration_check
 
 TESTS = tests/list_check tests/first_merge_check tests/second_merge_check \
-        tests/third_merge_check tests/fourth_merge_check
+        tests/third_merge_check tests/fourth_merge_check \
+        tests/integration_check
 
 EXTRA_DIST += tests/first_merge_check tests/second_merge_check \
               tests/third_merge_check tests/fourth_merge_check \
@@ -118,6 +119,11 @@ endif
 tests_list_check_SOURCES = tests/list_check.c src/reflist.c src/thread-pool.c \
                            src/utils.c
 
+tests_integration_check_SOURCES = tests/integration_check.c src/peaks.c \
+                                  src/reflist.c src/cell.c src/detector.c \
+                                  src/thread-pool.c src/utils.c src/image.c \
+                                  src/hdf5-file.c
+
 INCLUDES = "-I$(top_srcdir)/data"
 
 EXTRA_DIST += src/cell.h src/hdf5-file.h src/image.h src/utils.h \
diff --git a/Makefile.in b/Makefile.in
index ab99bcbd..9571f064 100644
--- a/Makefile.in
+++ b/Makefile.in
@@ -41,10 +41,12 @@ bin_PROGRAMS = src/pattern_sim$(EXEEXT) src/process_hkl$(EXEEXT) \
 	src/render_hkl$(EXEEXT) src/calibrate_detector$(EXEEXT) \
 	src/partialator$(EXEEXT) src/check_hkl$(EXEEXT) \
 	src/sum_stack$(EXEEXT) $(am__EXEEXT_1) $(am__EXEEXT_2)
-noinst_PROGRAMS = tests/list_check$(EXEEXT) $(am__EXEEXT_3)
+noinst_PROGRAMS = tests/list_check$(EXEEXT) \
+	tests/integration_check$(EXEEXT) $(am__EXEEXT_3)
 TESTS = tests/list_check$(EXEEXT) tests/first_merge_check \
 	tests/second_merge_check tests/third_merge_check \
-	tests/fourth_merge_check $(am__EXEEXT_3)
+	tests/fourth_merge_check tests/integration_check$(EXEEXT) \
+	$(am__EXEEXT_3)
 @BUILD_HDFSEE_TRUE@am__append_1 = src/hdfsee
 @BUILD_CUBEIT_TRUE@am__append_2 = src/cubeit
 @HAVE_OPENCL_TRUE@am__append_3 = src/diffraction-gpu.c src/cl-utils.c
@@ -240,6 +242,16 @@ am__tests_gpu_sim_check_SOURCES_DIST = tests/gpu_sim_check.c \
 tests_gpu_sim_check_OBJECTS = $(am_tests_gpu_sim_check_OBJECTS)
 tests_gpu_sim_check_LDADD = $(LDADD)
 tests_gpu_sim_check_DEPENDENCIES = $(top_builddir)/lib/libgnu.a
+am_tests_integration_check_OBJECTS =  \
+	tests/integration_check.$(OBJEXT) src/peaks.$(OBJEXT) \
+	src/reflist.$(OBJEXT) src/cell.$(OBJEXT) \
+	src/detector.$(OBJEXT) src/thread-pool.$(OBJEXT) \
+	src/utils.$(OBJEXT) src/image.$(OBJEXT) \
+	src/hdf5-file.$(OBJEXT)
+tests_integration_check_OBJECTS =  \
+	$(am_tests_integration_check_OBJECTS)
+tests_integration_check_LDADD = $(LDADD)
+tests_integration_check_DEPENDENCIES = $(top_builddir)/lib/libgnu.a
 am_tests_list_check_OBJECTS = tests/list_check.$(OBJEXT) \
 	src/reflist.$(OBJEXT) src/thread-pool.$(OBJEXT) \
 	src/utils.$(OBJEXT)
@@ -273,7 +285,7 @@ SOURCES = $(src_calibrate_detector_SOURCES) $(src_check_hkl_SOURCES) \
 	$(src_pattern_sim_SOURCES) $(src_powder_plot_SOURCES) \
 	$(src_process_hkl_SOURCES) $(src_render_hkl_SOURCES) \
 	$(src_sum_stack_SOURCES) $(tests_gpu_sim_check_SOURCES) \
-	$(tests_list_check_SOURCES)
+	$(tests_integration_check_SOURCES) $(tests_list_check_SOURCES)
 DIST_SOURCES = $(src_calibrate_detector_SOURCES) \
 	$(src_check_hkl_SOURCES) $(src_compare_hkl_SOURCES) \
 	$(am__src_cubeit_SOURCES_DIST) $(src_get_hkl_SOURCES) \
@@ -283,7 +295,7 @@ DIST_SOURCES = $(src_calibrate_detector_SOURCES) \
 	$(src_process_hkl_SOURCES) $(src_render_hkl_SOURCES) \
 	$(src_sum_stack_SOURCES) \
 	$(am__tests_gpu_sim_check_SOURCES_DIST) \
-	$(tests_list_check_SOURCES)
+	$(tests_integration_check_SOURCES) $(tests_list_check_SOURCES)
 RECURSIVE_TARGETS = all-recursive check-recursive dvi-recursive \
 	html-recursive info-recursive install-data-recursive \
 	install-dvi-recursive install-exec-recursive \
@@ -665,6 +677,11 @@ src_partialator_SOURCES = src/partialator.c src/cell.c src/hdf5-file.c \
 tests_list_check_SOURCES = tests/list_check.c src/reflist.c src/thread-pool.c \
                            src/utils.c
 
+tests_integration_check_SOURCES = tests/integration_check.c src/peaks.c \
+                                  src/reflist.c src/cell.c src/detector.c \
+                                  src/thread-pool.c src/utils.c src/image.c \
+                                  src/hdf5-file.c
+
 INCLUDES = "-I$(top_srcdir)/data"
 crystfeldir = $(datadir)/crystfel
 crystfel_DATA = data/diffraction.cl data/defs.h data/hdfsee.ui
@@ -915,6 +932,11 @@ tests/gpu_sim_check.$(OBJEXT): tests/$(am__dirstamp) \
 tests/gpu_sim_check$(EXEEXT): $(tests_gpu_sim_check_OBJECTS) $(tests_gpu_sim_check_DEPENDENCIES) tests/$(am__dirstamp)
 	@rm -f tests/gpu_sim_check$(EXEEXT)
 	$(AM_V_CCLD)$(LINK) $(tests_gpu_sim_check_OBJECTS) $(tests_gpu_sim_check_LDADD) $(LIBS)
+tests/integration_check.$(OBJEXT): tests/$(am__dirstamp) \
+	tests/$(DEPDIR)/$(am__dirstamp)
+tests/integration_check$(EXEEXT): $(tests_integration_check_OBJECTS) $(tests_integration_check_DEPENDENCIES) tests/$(am__dirstamp)
+	@rm -f tests/integration_check$(EXEEXT)
+	$(AM_V_CCLD)$(LINK) $(tests_integration_check_OBJECTS) $(tests_integration_check_LDADD) $(LIBS)
 tests/list_check.$(OBJEXT): tests/$(am__dirstamp) \
 	tests/$(DEPDIR)/$(am__dirstamp)
 tests/list_check$(EXEEXT): $(tests_list_check_OBJECTS) $(tests_list_check_DEPENDENCIES) tests/$(am__dirstamp)
@@ -963,6 +985,7 @@ mostlyclean-compile:
 	-rm -f src/thread-pool.$(OBJEXT)
 	-rm -f src/utils.$(OBJEXT)
 	-rm -f tests/gpu_sim_check.$(OBJEXT)
+	-rm -f tests/integration_check.$(OBJEXT)
 	-rm -f tests/list_check.$(OBJEXT)
 
 distclean-compile:
@@ -1008,6 +1031,7 @@ distclean-compile:
 @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/thread-pool.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/utils.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@tests/$(DEPDIR)/gpu_sim_check.Po@am__quote@
+@AMDEP_TRUE@@am__include@ @am__quote@tests/$(DEPDIR)/integration_check.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@tests/$(DEPDIR)/list_check.Po@am__quote@
 
 .c.o:
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;
+}