aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--Makefile.am7
-rw-r--r--tests/.gitignore1
-rw-r--r--tests/prof2d_check.c192
3 files changed, 198 insertions, 2 deletions
diff --git a/Makefile.am b/Makefile.am
index 270ddf21..f117d67d 100644
--- a/Makefile.am
+++ b/Makefile.am
@@ -10,7 +10,8 @@ noinst_PROGRAMS = tests/list_check tests/integration_check \
tests/pr_p_gradient_check tests/symmetry_check \
tests/centering_check tests/transformation_check \
tests/cell_check tests/pr_l_gradient_check \
- tests/pr_pl_gradient_check tests/ring_check
+ tests/pr_pl_gradient_check tests/ring_check \
+ tests/prof2d_check
MERGE_CHECKS = tests/first_merge_check tests/second_merge_check \
tests/third_merge_check tests/fourth_merge_check
@@ -25,7 +26,7 @@ PARTIAL_CHECKS = tests/partialator_merge_check_1 \
TESTS = tests/list_check $(MERGE_CHECKS) $(PARTIAL_CHECKS) \
tests/integration_check \
tests/symmetry_check tests/centering_check tests/transformation_check \
- tests/cell_check tests/ring_check
+ tests/cell_check tests/ring_check tests/prof2d_check
EXTRA_DIST += $(MERGE_CHECKS) $(PARTIAL_CHECKS)
@@ -88,6 +89,8 @@ tests_list_check_SOURCES = tests/list_check.c
tests_integration_check_SOURCES = tests/integration_check.c
+tests_prof2d_check_SOURCES = tests/prof2d_check.c
+
tests_symmetry_check_SOURCES = tests/symmetry_check.c
tests_pr_p_gradient_check_SOURCES = tests/pr_p_gradient_check.c \
diff --git a/tests/.gitignore b/tests/.gitignore
index fabb48f3..1f76a5d5 100644
--- a/tests/.gitignore
+++ b/tests/.gitignore
@@ -12,3 +12,4 @@ transformation_check
cell_check
ring_check
.dirstamp
+prof2d_check
diff --git a/tests/prof2d_check.c b/tests/prof2d_check.c
new file mode 100644
index 00000000..24cf4be5
--- /dev/null
+++ b/tests/prof2d_check.c
@@ -0,0 +1,192 @@
+/*
+ * prof2d_check.c
+ *
+ * Check 2D profile fitting
+ *
+ * Copyright © 2013 Thomas White <taw@physics.org>
+ *
+ * This file is part of CrystFEL.
+ *
+ * CrystFEL is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * CrystFEL is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with CrystFEL. If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+
+#include <stdlib.h>
+#include <stdio.h>
+
+#include <image.h>
+#include <utils.h>
+#include <beam-parameters.h>
+#include <histogram.h>
+
+#include "../libcrystfel/src/integration.c"
+
+
+#define ADD_PX(fs, ss, val) \
+ if ( ((fs)>0) && ((ss)>0) && ((fs)<w) && ((ss)<h) ) { \
+ image.dp[0][(fs)+w*(ss)] += (val); \
+ }
+
+int main(int argc, char *argv[])
+{
+ struct image image;
+ int fs, ss;
+ FILE *fh;
+ unsigned int seed;
+ int fail = 0;
+ const int w = 1024;
+ const int h = 1024;
+ RefList *list;
+ RefListIterator *iter;
+ Reflection *refl;
+ UnitCell *cell;
+ Crystal *cr;
+ const int ir_inn = 2;
+ const int ir_mid = 4;
+ const int ir_out = 6;
+ Histogram *hi;
+ double esd_sum = 0.0;
+ int n = 0;
+
+ fh = fopen("/dev/urandom", "r");
+ fread(&seed, sizeof(seed), 1, fh);
+ fclose(fh);
+ srand(seed);
+
+ image.flags = NULL;
+ image.beam = NULL;
+ image.lambda = ph_eV_to_lambda(9000.0);
+ image.bw = 0.000001;
+ image.div = 0.0;
+
+ image.det = calloc(1, sizeof(struct detector));
+ image.det->n_panels = 1;
+ image.det->panels = calloc(1, sizeof(struct panel));
+
+ image.dp = calloc(1, sizeof(float *));
+ image.bad = calloc(1, sizeof(int *));
+
+ image.width = w;
+ image.height = h;
+ image.det->panels[0].min_fs = 0;
+ image.det->panels[0].max_fs = w;
+ image.det->panels[0].min_ss = 0;
+ image.det->panels[0].max_ss = h;
+ image.det->panels[0].w = w;
+ image.det->panels[0].h = h;
+ 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 = -w/2;
+ image.det->panels[0].cny = -h/2;
+ image.det->panels[0].clen = 60.0e-3;
+ image.det->panels[0].res = 100000; /* 10 px per mm */
+ image.det->panels[0].adu_per_eV = 10.0/9000.0; /* 10 adu/ph */
+ image.det->panels[0].max_adu = +INFINITY; /* No cutoff */
+
+ image.det->furthest_out_panel = &image.det->panels[0];
+ image.det->furthest_out_fs = 0;
+ image.det->furthest_out_ss = 0;
+
+ image.dp[0] = malloc(w*h*sizeof(float));
+ memset(image.dp[0], 0, w*h*sizeof(float));
+ image.bad[0] = malloc(w*h*sizeof(int));
+ memset(image.bad[0], 0, w*h*sizeof(int));
+
+ cell = cell_new();
+ cell_set_lattice_type(cell, L_CUBIC);
+ cell_set_centering(cell, 'P');
+ cell_set_parameters(cell, 300.0e-10, 300.0e-10, 300.0e-10,
+ deg2rad(90.0), deg2rad(90.0), deg2rad(90.0));
+ cell = cell_rotate(cell, random_quaternion());
+
+ cr = crystal_new();
+ crystal_set_profile_radius(cr, 0.001e9);
+ crystal_set_mosaicity(cr, 0.0); /* radians */
+ crystal_set_image(cr, &image);
+ crystal_set_cell(cr, cell);
+
+ image.n_crystals = 1;
+ image.crystals = &cr;
+
+ list = find_intersections(&image, cr);
+
+ for ( fs=0; fs<w; fs++ ) {
+ for ( ss=0; ss<h; ss++ ) {
+ image.dp[0][fs+w*ss] = 10.0*poisson_noise(40);
+ }
+ }
+
+ for ( refl = first_refl(list, &iter);
+ refl != NULL;
+ refl = next_refl(refl, iter) )
+ {
+ double pfs, pss;
+ int fs, ss;
+ const int pk_ph = 1000;
+
+ get_detector_pos(refl, &pfs, &pss);
+ fs = pfs; ss = pss;
+
+ ADD_PX(fs, ss, 10.0*poisson_noise(pk_ph));
+ ADD_PX(fs-1, ss, 10.0*poisson_noise(pk_ph));
+ ADD_PX(fs+1, ss, 10.0*poisson_noise(pk_ph));
+ ADD_PX(fs, ss-1, 10.0*poisson_noise(pk_ph));
+ ADD_PX(fs, ss+1, 10.0*poisson_noise(pk_ph));
+
+ }
+
+ reflist_free(list); /* integrate_prof2d() will predict again */
+ integrate_prof2d(INTEGRATION_PROF2D, cr, &image, INTDIAG_NONE, 0, 0, 0,
+ ir_inn, ir_mid, ir_out);
+
+ hi = histogram_init();
+
+ list = crystal_get_reflections(cr);
+ for ( refl = first_refl(list, &iter);
+ refl != NULL;
+ refl = next_refl(refl, iter) )
+ {
+ if ( get_redundancy(refl) == 0 ) continue;
+ histogram_add_value(hi, get_intensity(refl));
+ esd_sum += get_esd_intensity(refl);
+ n++;
+ }
+
+ printf("Mean calculated sigma(I) = %.2f (%i measurements)\n",
+ esd_sum / n, n);
+
+ histogram_show(hi);
+
+ histogram_free(hi);
+ free(image.beam);
+ free(image.det->panels);
+ free(image.det);
+ free(image.dp[0]);
+ free(image.dp);
+
+ if ( fail ) return 1;
+
+ return 0;
+}