From 2e2386192d946988811cb75bac649877ba0fe0c7 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Tue, 19 Nov 2013 15:37:32 +0100 Subject: Break off integrate_rings_once() --- tests/ring_check.c | 314 +++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 314 insertions(+) create mode 100644 tests/ring_check.c (limited to 'tests/ring_check.c') diff --git a/tests/ring_check.c b/tests/ring_check.c new file mode 100644 index 00000000..ac871581 --- /dev/null +++ b/tests/ring_check.c @@ -0,0 +1,314 @@ +/* + * integration_check.c + * + * Check peak integration + * + * Copyright © 2012 Thomas White + * Copyright © 2012 Andrew Martin + * + * 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 . + * + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + + +#include +#include + +#include +#include +#include + +#include "../libcrystfel/src/peaks.c" + + +/* The third integration check draws a Poisson background and checks that, on + * average, it gets subtracted by the background subtraction. */ +static void third_integration_check(struct image *image, int n_trials, + int *fail) +{ + double mean_intensity = 0.0; + double mean_bg = 0.0; + double mean_max = 0.0; + double mean_sigma = 0.0; + int i; + int fs, ss; + int nfail = 0; + + for ( i=0; iwidth; fs++ ) { + for ( ss=0; ssheight; ss++ ) { + image->data[fs+image->width*ss] = poisson_noise(1000.0); + } + } + + r = integrate_peak(image, 64, 64, &fsp, &ssp, + &intensity, &sigma, 10.0, 15.0, 17.0, + NULL, NULL); + + if ( r == 0 ) { + mean_intensity += intensity; + mean_sigma += sigma; + } else { + nfail++; + } + + } + mean_intensity /= n_trials; + mean_bg /= n_trials; + mean_max /= n_trials; + mean_sigma /= n_trials; + + STATUS(" Third check (mean values): intensity = %.2f, sigma = %.2f," + " integration failed %i/%i times\n", + mean_intensity, mean_sigma, nfail, n_trials); + +/* These values are always wrong, because the integration sucks */ +// if ( fabs(mean_intensity) > 5.0 ) { +// ERROR("Mean intensity should be close to zero.\n"); +// *fail = 1; +// } +// if ( fabs(mean_intensity) > mean_sigma/10.0 ) { +// ERROR("Mean intensity should be much less than mean sigma.\n"); +// *fail = 1; +// } +} + + +/* The fourth integration check draws a Poisson background and draws a peak on + * top of it, then checks that the intensity of the peak is correctly recovered + * accounting for the background. */ +static void fourth_integration_check(struct image *image, int n_trials, + int *fail) +{ + double mean_intensity = 0.0; + double mean_sigma = 0.0; + int i; + int fs, ss; + int pcount = 0; + int nfail = 0; + + for ( i=0; iwidth; fs++ ) { + for ( ss=0; ssheight; ss++ ) { + int idx = fs+image->width*ss; + image->data[idx] = poisson_noise(1000.0); + if ( (fs-64)*(fs-64) + (ss-64)*(ss-64) > 9*9 ) continue; + image->data[idx] += 1000.0; + pcount++; + } + } + + r = integrate_peak(image, 64, 64, &fsp, &ssp, + &intensity, &sigma, 10.0, 15.0, 17.0, + NULL, NULL); + + if ( r == 0 ) { + mean_intensity += intensity; + mean_sigma += sigma; + } else { + nfail++; + } + + } + mean_intensity /= n_trials; + mean_sigma /= n_trials; + pcount /= n_trials; + + STATUS(" Fourth check (mean values): intensity = %.2f, sigma = %.2f," + " integration failed %i/%i times\n", + mean_intensity, mean_sigma, nfail, n_trials); + + if ( fabs(mean_intensity - pcount*1000.0) > 4000.0 ) { + ERROR("Mean intensity should be close to %f\n", pcount*1000.0); + *fail = 1; + } + if ( fabs(mean_intensity) < mean_sigma ) { + ERROR("Mean intensity should be greater than mean sigma.\n"); + *fail = 1; + } +} + + +int main(int argc, char *argv[]) +{ + struct image image; + double fsp, ssp, intensity, sigma; + int fs, ss; + FILE *fh; + unsigned int seed; + int fail = 0; + const int n_trials = 100; + int r, npx; + double ex; + + fh = fopen("/dev/urandom", "r"); + fread(&seed, sizeof(seed), 1, fh); + fclose(fh); + srand(seed); + + image.data = malloc(128*128*sizeof(float)); + image.flags = NULL; + image.beam = NULL; + image.lambda = ph_eV_to_lambda(1000.0); + + 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].adu_per_eV = 1.0/1000.0; /* -> 1 adu per photon */ + image.det->panels[0].max_adu = +INFINITY; /* No cutoff */ + + image.width = 128; + image.height = 128; + memset(image.data, 0, 128*128*sizeof(float)); + + image.n_crystals = 0; + image.crystals = NULL; + + /* First check: no intensity -> no peak, or very low intensity */ + r = integrate_peak(&image, 64, 64, &fsp, &ssp, &intensity, &sigma, + 10.0, 15.0, 17.0, NULL, NULL); + STATUS(" First check: integrate_peak() returned %i", r); + if ( r == 0 ) { + + STATUS(", intensity = %.2f, sigma = %.2f\n", intensity, sigma); + + if ( fabs(intensity) > 0.01 ) { + ERROR("Intensity should be very close to zero.\n"); + fail = 1; + } + + } else { + STATUS(" (correct)\n"); + } + + /* Second check: uniform peak gives correct I and low sigma(I) */ + npx = 0; + for ( fs=0; fs 9*9 ) continue; + image.data[fs+image.width*ss] = 1000.0; + npx++; + } + } + + r = integrate_peak(&image, 64, 64, &fsp, &ssp, &intensity, &sigma, + 10.0, 15.0, 17.0, NULL, NULL); + if ( r ) { + ERROR(" Second check: integrate_peak() returned %i (wrong).\n", + r); + fail = 1; + } else { + + STATUS(" Second check: intensity = %.2f, sigma = %.2f\n", + intensity, sigma); + + ex = npx*1000.0; + if ( within_tolerance(ex, intensity, 1.0) == 0 ) { + ERROR("Intensity should be close to %f\n", ex); + fail = 1; + } + + ex = sqrt(npx*1000.0); + if ( within_tolerance(ex, sigma, 1.0) == 0 ) { + ERROR("Sigma should be roughly %f.\n", ex); + fail = 1; + } + + } + + /* Third check: Poisson background should get mostly subtracted */ + third_integration_check(&image, n_trials, &fail); + + /* Fourth check: peak on Poisson background */ + fourth_integration_check(&image, n_trials, &fail); + + /* Fifth check: uniform peak on uniform background */ + npx = 0; + for ( fs=0; fs 9*9 ) continue; + image.data[fs+image.width*ss] += 1000.0; + npx++; + } + } + + r = integrate_peak(&image, 64, 64, &fsp, &ssp, &intensity, &sigma, + 10.0, 15.0, 17.0, NULL, NULL); + if ( r ) { + ERROR(" Fifth check: integrate_peak() returned %i (wrong).\n", + r); + fail = 1; + } else { + + STATUS(" Fifth check: intensity = %.2f, sigma = %.2f\n", + intensity, sigma); + + ex = npx*1000.0; + if ( within_tolerance(ex, intensity, 1.0) == 0 ) { + ERROR("Intensity should be close to %f\n", ex); + fail = 1; + } + + ex = sqrt(npx*1000.0); + if ( within_tolerance(ex, sigma, 1.0) == 0 ) { + ERROR("Sigma should be roughly %f.\n", ex); + fail = 1; + } + + } + + + free(image.beam); + free(image.det->panels); + free(image.det); + free(image.data); + + if ( fail ) return 1; + + return 0; +} -- cgit v1.2.3