diff options
| -rw-r--r-- | Makefile.am | 7 | ||||
| -rw-r--r-- | tests/.gitignore | 1 | ||||
| -rw-r--r-- | tests/prof2d_check.c | 192 |
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; +} |