diff options
Diffstat (limited to 'tests')
| -rw-r--r-- | tests/.gitignore | 1 | ||||
| -rw-r--r-- | tests/integration_check.c | 185 | ||||
| -rw-r--r-- | tests/ring_check.c | 314 |
3 files changed, 319 insertions, 181 deletions
diff --git a/tests/.gitignore b/tests/.gitignore index 31ce158f..fabb48f3 100644 --- a/tests/.gitignore +++ b/tests/.gitignore @@ -10,4 +10,5 @@ symmetry_check centering_check transformation_check cell_check +ring_check .dirstamp diff --git a/tests/integration_check.c b/tests/integration_check.c index ac871581..4b302955 100644 --- a/tests/integration_check.c +++ b/tests/integration_check.c @@ -1,10 +1,9 @@ /* * integration_check.c * - * Check peak integration + * Check reflection integration * - * Copyright © 2012 Thomas White <taw@physics.org> - * Copyright © 2012 Andrew Martin <andrew.martin@desy.de> + * Copyright © 2013 Thomas White <taw@physics.org> * * This file is part of CrystFEL. * @@ -35,126 +34,9 @@ #include <utils.h> #include <beam-parameters.h> -#include "../libcrystfel/src/peaks.c" +#include "../libcrystfel/src/integration.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; i<n_trials; i++ ) { - - double intensity, sigma; - double fsp, ssp; - int r; - - for ( fs=0; fs<image->width; fs++ ) { - for ( ss=0; ss<image->height; 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; i<n_trials; i++ ) { - - double intensity, sigma; - double fsp, ssp; - int r; - - for ( fs=0; fs<image->width; fs++ ) { - for ( ss=0; ss<image->height; 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[]) { @@ -164,7 +46,6 @@ int main(int argc, char *argv[]) FILE *fh; unsigned int seed; int fail = 0; - const int n_trials = 100; int r, npx; double ex; @@ -208,65 +89,7 @@ int main(int argc, char *argv[]) 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<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; - 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 */ + /* Uniform peak on uniform background */ npx = 0; for ( fs=0; fs<image.width; fs++ ) { for ( ss=0; ss<image.height; ss++ ) { 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 <taw@physics.org> + * Copyright © 2012 Andrew Martin <andrew.martin@desy.de> + * + * 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 "../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; i<n_trials; i++ ) { + + double intensity, sigma; + double fsp, ssp; + int r; + + for ( fs=0; fs<image->width; fs++ ) { + for ( ss=0; ss<image->height; 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; i<n_trials; i++ ) { + + double intensity, sigma; + double fsp, ssp; + int r; + + for ( fs=0; fs<image->width; fs++ ) { + for ( ss=0; ss<image->height; 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<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; + 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<image.width; fs++ ) { + for ( ss=0; ss<image.height; ss++ ) { + image.data[fs+image.width*ss] = 1000.0; + if ( (fs-64)*(fs-64) + (ss-64)*(ss-64) > 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; +} |