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() --- Makefile.am | 6 +- libcrystfel/src/integration.c | 350 +++++++++++++++++++++++++++++++----------- tests/.gitignore | 1 + tests/integration_check.c | 185 +--------------------- tests/ring_check.c | 314 +++++++++++++++++++++++++++++++++++++ 5 files changed, 584 insertions(+), 272 deletions(-) create mode 100644 tests/ring_check.c diff --git a/Makefile.am b/Makefile.am index 69baeaee..270ddf21 100644 --- a/Makefile.am +++ b/Makefile.am @@ -10,7 +10,7 @@ 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/pr_pl_gradient_check tests/ring_check MERGE_CHECKS = tests/first_merge_check tests/second_merge_check \ tests/third_merge_check tests/fourth_merge_check @@ -25,7 +25,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/cell_check tests/ring_check EXTRA_DIST += $(MERGE_CHECKS) $(PARTIAL_CHECKS) @@ -103,6 +103,8 @@ tests_centering_check_SOURCES = tests/centering_check.c tests_transformation_check_SOURCES = tests/transformation_check.c +tests_ring_check_SOURCES = tests/ring_check.c + tests_cell_check_SOURCES = tests/cell_check.c INCLUDES = -I$(top_srcdir)/libcrystfel/src -I$(top_srcdir)/data diff --git a/libcrystfel/src/integration.c b/libcrystfel/src/integration.c index 1b68cf0b..5bd520e4 100644 --- a/libcrystfel/src/integration.c +++ b/libcrystfel/src/integration.c @@ -43,6 +43,7 @@ #endif #include "reflist.h" +#include "reflist-utils.h" #include "cell.h" #include "crystal.h" #include "cell-utils.h" @@ -155,6 +156,8 @@ enum boxmask_val struct intcontext { + IntegrationMethod meth; + int halfw; int w; enum boxmask_val *bm; /* Box mask */ @@ -175,6 +178,10 @@ struct intcontext int ir_inn; int ir_mid; int ir_out; + + int n_saturated; + int n_implausible; + double limit; }; @@ -1099,8 +1106,8 @@ static double calc_sigma(struct intcontext *ic, struct peak_box *bx) } -static void mean_var_background(struct intcontext *ic, struct peak_box *bx, - double *pmean, double *pvar) +static void mean_var_area(struct intcontext *ic, struct peak_box *bx, + enum boxmask_val v, double *pmean, double *pvar) { int p, q; double sum = 0.0; @@ -1110,7 +1117,7 @@ static void mean_var_background(struct intcontext *ic, struct peak_box *bx, for ( p=0; pw; p++ ) { for ( q=0; qw; q++ ) { - if ( bx->bm[p + ic->w*q] != BM_BG ) continue; + if ( bx->bm[p + ic->w*q] != v ) continue; sum += bx->a*p + bx->b*q + bx->c; n++; } @@ -1120,7 +1127,7 @@ static void mean_var_background(struct intcontext *ic, struct peak_box *bx, n = 0; for ( p=0; pw; p++ ) { for ( q=0; qw; q++ ) { - if ( bx->bm[p + ic->w*q] != BM_BG ) continue; + if ( bx->bm[p + ic->w*q] != v ) continue; var += pow(bx->a*p + bx->b*q + bx->c - mean, 2.0); n++; } @@ -1336,6 +1343,9 @@ static void integrate_prof2d(IntegrationMethod meth, Crystal *cr, ic.halfw = ir_out; ic.image = image; ic.k = 1.0/image->lambda; + ic.meth = meth; + ic.n_saturated = 0; + ic.n_implausible = 0; ic.cell = cell; if ( init_intcontext(&ic) ) { ERROR("Failed to initialise integration.\n"); @@ -1483,6 +1493,132 @@ static void integrate_prof2d(IntegrationMethod meth, Crystal *cr, } +static void integrate_rings_once(Reflection *refl, struct image *image, + struct intcontext *ic, UnitCell *cell) +{ + double pfs, pss; + signed int h, k, l; + struct peak_box *bx; + int pn; + struct panel *p; + int fid_fs, fid_ss; /* Center coordinates, rounded, + * in overall data block */ + int cfs, css; /* Corner coordinates */ + double intensity; + double sigma; + int saturated; + double one_over_d; + int r; + double bgmean, sig2_bg, sig2_poisson, aduph; + double pkmean, sig2_pk; + + set_redundancy(refl, 0); + + get_detector_pos(refl, &pfs, &pss); + + /* Explicit truncation of digits after the decimal point. + * This is actually the correct thing to do here, not + * e.g. lrint(). pfs/pss is the position of the spot, measured + * in numbers of pixels, from the panel corner (not the center + * of the first pixel). So any coordinate from 2.0 to 2.9999 + * belongs to pixel index 2. */ + fid_fs = pfs; + fid_ss = pss; + pn = find_panel_number(image->det, fid_fs, fid_ss); + p = &image->det->panels[pn]; + + cfs = (fid_fs-p->min_fs) - ic->halfw; + css = (fid_ss-p->min_ss) - ic->halfw; + + bx = add_box(ic); + bx->refl = refl; + bx->cfs = cfs; + bx->css = css; + bx->p = p; + bx->pn = pn; + get_indices(refl, &h, &k, &l); + bx->verbose = VERBOSITY; + + if ( ic->meth & INTEGRATION_CENTER ) { + r = center_and_check_box(ic, bx, &saturated); + } else { + r = check_box(ic, bx, &saturated); + if ( !r ) { + fit_bg(ic, bx); + if ( bx->verbose ) show_peak_box(ic, bx); + } + bx->offs_fs = 0.0; + bx->offs_ss = 0.0; + } + if ( r ) { + delete_box(ic, bx); + return; + } + + if ( saturated ) { + ic->n_saturated++; + if ( !(ic->meth & INTEGRATION_SATURATED) ) { + delete_box(ic, bx); + return; + } + } + + intensity = tentative_intensity(ic, bx); + mean_var_area(ic, bx, BM_BG, &bgmean, &sig2_bg); + + mean_var_area(ic, bx, BM_PK, &pkmean, &sig2_pk); + if ( bx->verbose ) { + STATUS("bg mean, var = %.2f, %.2f\n", bgmean, sig2_bg); + STATUS("pk mean, var = %.2f, %.2f\n", pkmean, sig2_pk); + } + + aduph = bx->p->adu_per_eV * ph_lambda_to_eV(ic->image->lambda); + sig2_poisson = aduph * intensity; + + /* If intensity is within one photon of nothing, set the Poisson + * error to be one photon */ + if ( fabs(intensity / aduph) < 1.0 ) { + sig2_poisson = aduph; + if ( bx->verbose ) { + STATUS("number of photons less than 1\n"); + } + } + + /* If intensity is negative by more than one photon, assume that + * the peak is in the background and add a Poisson error of + * appropriate size */ + if ( intensity < -aduph ) { + sig2_poisson = -aduph*intensity; + if ( bx->verbose ) { + STATUS("very negative (%.2f photons)\n", + intensity/aduph); + } + } + + sigma = sqrt(sig2_poisson + bx->m*sig2_bg); + + if ( intensity < -5.0*sigma ) { + delete_box(ic, bx); + ic->n_implausible++; + return; + } + + /* Record intensity and set redundancy to 1 */ + bx->intensity = intensity; + set_intensity(refl, intensity); + set_esd_intensity(refl, sigma); + set_redundancy(refl, 1); + + one_over_d = resolution(cell, h, k, l); + if ( one_over_d > ic->limit ) ic->limit = one_over_d; + + /* Update position */ + pfs += bx->offs_fs; + pss += bx->offs_ss; + set_detector_pos(refl, 0.0, pfs, pss); +} + + static void integrate_rings(IntegrationMethod meth, Crystal *cr, struct image *image, double ir_inn, double ir_mid, double ir_out) @@ -1492,8 +1628,6 @@ static void integrate_rings(IntegrationMethod meth, Crystal *cr, RefListIterator *iter; UnitCell *cell; struct intcontext ic; - int n_saturated = 0; - double limit = 0.0; list = find_intersections(image, cr); if ( list == NULL ) return; @@ -1505,10 +1639,13 @@ static void integrate_rings(IntegrationMethod meth, Crystal *cr, ic.halfw = ir_out; ic.image = image; ic.k = 1.0/image->lambda; + ic.n_saturated = 0; + ic.n_implausible = 0; ic.cell = cell; ic.ir_inn = ir_inn; ic.ir_mid = ir_mid; ic.ir_out = ir_out; + ic.limit = 0.0; if ( init_intcontext(&ic) ) { ERROR("Failed to initialise integration.\n"); return; @@ -1519,101 +1656,131 @@ static void integrate_rings(IntegrationMethod meth, Crystal *cr, refl != NULL; refl = next_refl(refl, iter) ) { - double pfs, pss; - signed int h, k, l; - struct peak_box *bx; - int pn; - struct panel *p; - int fid_fs, fid_ss; /* Center coordinates, rounded, - * in overall data block */ - int cfs, css; /* Corner coordinates */ - double intensity; - double sigma; - int saturated; - double one_over_d; - int r; - double bgmean, sig2_bg, sig2_poisson, aduph; + integrate_rings_once(refl, image, &ic, cell); + } - set_redundancy(refl, 0); + refine_rigid_groups(&ic); - get_detector_pos(refl, &pfs, &pss); + free_intcontext(&ic); - /* Explicit truncation of digits after the decimal point. - * This is actually the correct thing to do here, not - * e.g. lrint(). pfs/pss is the position of the spot, measured - * in numbers of pixels, from the panel corner (not the center - * of the first pixel). So any coordinate from 2.0 to 2.9999 - * belongs to pixel index 2. */ - fid_fs = pfs; - fid_ss = pss; - pn = find_panel_number(image->det, fid_fs, fid_ss); - p = &image->det->panels[pn]; + crystal_set_num_saturated_reflections(cr, ic.n_saturated); + crystal_set_resolution_limit(cr, ic.limit); + crystal_set_reflections(cr, list); - cfs = (fid_fs-p->min_fs) - ic.halfw; - css = (fid_ss-p->min_ss) - ic.halfw; + if ( ic.n_implausible ) { + STATUS("Warning: %i implausibly negative reflections.\n", + ic.n_implausible); + } +} - bx = add_box(&ic); - bx->refl = refl; - bx->cfs = cfs; - bx->css = css; - bx->p = p; - bx->pn = pn; - get_indices(refl, &h, &k, &l); - bx->verbose = VERBOSITY; - if ( meth & INTEGRATION_CENTER ) { - r = center_and_check_box(&ic, bx, &saturated); - } else { - r = check_box(&ic, bx, &saturated); - if ( !r ) { - fit_bg(&ic, bx); - if ( bx->verbose ) show_peak_box(&ic, bx); - } - bx->offs_fs = 0.0; - bx->offs_ss = 0.0; - } - if ( r ) { - delete_box(&ic, bx); - continue; - } +static void resolution_cutoff(RefList *list, UnitCell *cell) +{ + double *rmins; + double *rmaxs; + double rmin, rmax; + double total_vol, vol_per_shell; + int i; + int *sh_nref; + double *sh_snr; + Reflection **highest_snr; + Reflection *refl; + RefListIterator *iter; - if ( saturated ) { - n_saturated++; - if ( !(meth & INTEGRATION_SATURATED) ) { - delete_box(&ic, bx); - continue; - } - } + const int nshells = 100; + const double min_snr = 10.0; + + rmins = malloc(nshells*sizeof(double)); + rmaxs = malloc(nshells*sizeof(double)); + sh_nref = malloc(nshells*sizeof(int)); + sh_snr = malloc(nshells*sizeof(double)); + highest_snr = malloc(nshells*sizeof(Reflection *)); + + if ( (rmins==NULL) || (rmaxs==NULL) || (sh_nref==NULL) + || (sh_snr==NULL) || (highest_snr==NULL) ) + { + ERROR("Couldn't allocate memory for resolution shells.\n"); + return; + } - intensity = tentative_intensity(&ic, bx); - mean_var_background(&ic, bx, &bgmean, &sig2_bg); - aduph = bx->p->adu_per_eV * ph_lambda_to_eV(ic.image->lambda); - sig2_poisson = aduph * intensity; - sigma = sqrt(sig2_poisson + bx->m*sig2_bg); + resolution_limits(list, cell, &rmin, &rmax); - /* Record intensity and set redundancy to 1 */ - bx->intensity = intensity; - set_intensity(refl, intensity); - set_esd_intensity(refl, sigma); - set_redundancy(refl, 1); + total_vol = pow(rmax, 3.0) - pow(rmin, 3.0); + vol_per_shell = total_vol / nshells; + rmins[0] = rmin; + for ( i=1; i limit ) limit = one_over_d; + double r; - /* Update position */ - pfs += bx->offs_fs; - pss += bx->offs_ss; - set_detector_pos(refl, 0.0, pfs, pss); + r = vol_per_shell + pow(rmins[i-1], 3.0); + r = pow(r, 1.0/3.0); + + /* Shells of constant volume */ + rmaxs[i-1] = r; + rmins[i] = r; } + rmaxs[nshells-1] = rmax; - refine_rigid_groups(&ic); + for ( i=0; irmins[j]) && (res<=rmaxs[j]) ) { + bin = j; + break; + } + } + + /* Allow for slight rounding errors */ + if ( (bin == -1) && (res <= rmins[0]) ) bin = 0; + if ( (bin == -1) && (res >= rmaxs[nshells-1]) ) bin = 0; + assert(bin != -1); + + sh_nref[bin]++; + sh_snr[bin] += snr; + } + + for ( i=nshells-1; i>=0; i-- ) { + if ( sh_nref[i] > 2 ) break; + } + double cen = (rmins[i]+rmaxs[i])/2.0; + STATUS("The highest shell is %i\n", i); + printf("%10.3f nm^-1 or %.2f A\n", cen/1e9, 10e9/cen); + + FILE *fh; + fh = fopen("cutoff.dat", "w"); + for ( i=0; icrystals[i], image, ir_inn, ir_mid, ir_out); - return; + break; case INTEGRATION_PROF2D : integrate_prof2d(meth, image->crystals[i], image, ir_inn, ir_mid, ir_out); - return; + break; default : ERROR("Unrecognised integration method %i\n", meth); - return; + break; } + /* Set resolution limit */ + resolution_cutoff(crystal_get_reflections(image->crystals[i]), + crystal_get_cell(image->crystals[i])); + + } } 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 - * Copyright © 2012 Andrew Martin + * Copyright © 2013 Thomas White * * This file is part of CrystFEL. * @@ -35,126 +34,9 @@ #include #include -#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; 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[]) { @@ -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 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 + * 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