aboutsummaryrefslogtreecommitdiff
path: root/tests/integration_check.c
diff options
context:
space:
mode:
Diffstat (limited to 'tests/integration_check.c')
-rw-r--r--tests/integration_check.c143
1 files changed, 143 insertions, 0 deletions
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;
+}