diff options
| -rw-r--r-- | src/peaks.c | 24 | ||||
| -rw-r--r-- | tests/integration_check.c | 327 |
2 files changed, 283 insertions, 68 deletions
diff --git a/src/peaks.c b/src/peaks.c index 4b87ba08..936601a3 100644 --- a/src/peaks.c +++ b/src/peaks.c @@ -222,20 +222,22 @@ int integrate_peak(struct image *image, int cfs, int css, if ( val > max ) max = val; - /* Inner mask */ + /* If outside inner mask, estimate noise from this region */ if ( fs*fs + ss*ss > lim_sq ) { - /* Estimate noise from this region */ - noise += val; //fabs(val); + + /* Noise */ + noise += val; noise_counts++; - /* Estimate the standard deviation of the noise */ - noise_meansq += fabs(val)*fabs(val) ; - continue; - } + noise_meansq += pow(val, 2.0); - pixel_counts ++; - total += val; - fsct += val*(cfs+fs); - ssct += val*(css+ss); + } else { + + /* Peak */ + pixel_counts++; + total += val; + fsct += val*(cfs+fs); + ssct += val*(css+ss); + } } } diff --git a/tests/integration_check.c b/tests/integration_check.c index 82278fdc..d6577032 100644 --- a/tests/integration_check.c +++ b/tests/integration_check.c @@ -25,11 +25,261 @@ #include "../src/beam-parameters.h" +/* 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; + + for ( i=0; i<n_trials; i++ ) { + + double intensity, bg, max, sigma; + double fsp, ssp; + + for ( fs=0; fs<image->width; fs++ ) { + for ( ss=0; ss<image->height; ss++ ) { + image->data[fs+image->width*ss] = poisson_noise(10.0); + } + } + integrate_peak(image, 64, 64, &fsp, &ssp, &intensity, + &bg, &max, &sigma, 0, 1, 1); + + mean_intensity += intensity; + mean_bg += bg; + mean_max += max; + mean_sigma += sigma; + + } + mean_intensity /= n_trials; + mean_bg /= n_trials; + mean_max /= n_trials; + mean_sigma /= n_trials; + + STATUS(" Third check (mean values): intensity = %.2f, bg = %.2f," + " max = %.2f, sigma = %.2f\n", + mean_intensity, mean_bg, mean_max, mean_sigma); + if ( fabs(mean_intensity) > 5.0 ) { + ERROR("Mean intensity should be close to zero.\n"); + *fail = 1; + } + if ( fabs(mean_bg-10.0) > 0.3 ) { + ERROR("Mean background should be close to ten.\n"); + *fail = 1; + } + if ( fabs(mean_intensity) > mean_sigma ) { + ERROR("Mean intensity should be 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_bg = 0.0; + double mean_max = 0.0; + double mean_sigma = 0.0; + int i; + int fs, ss; + int pcount = 0; + + for ( i=0; i<n_trials; i++ ) { + + double intensity, bg, max, sigma; + double fsp, ssp; + + for ( fs=0; fs<image->width; fs++ ) { + for ( ss=0; ss<image->height; ss++ ) { + int idx = fs+image->width*ss; + image->data[idx] = 0.0; + if ( (fs-64)*(fs-64) + (ss-64)*(ss-64) > 9*9 ) { + image->data[idx] = poisson_noise(10.0); + } else { + image->data[idx] += 1000.0; + pcount++; + } + } + } + integrate_peak(image, 64, 64, &fsp, &ssp, &intensity, + &bg, &max, &sigma, 0, 1, 1); + + mean_intensity += intensity; + mean_bg += bg; + mean_max += max; + mean_sigma += sigma; + + } + mean_intensity /= n_trials; + mean_bg /= n_trials; + mean_max /= n_trials; + mean_sigma /= n_trials; + pcount /= n_trials; + + STATUS(" Fourth check (mean values): intensity = %.2f, bg = %.2f," + " max = %.2f, sigma = %.2f\n", + mean_intensity, mean_bg, mean_max, mean_sigma); + 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_bg-10.0) > 0.3 ) { + ERROR("Mean background should be close to ten.\n"); + *fail = 1; + } + if ( fabs(mean_intensity) < mean_sigma ) { + ERROR("Mean intensity should be greater than mean sigma.\n"); + *fail = 1; + } +} + + +/* The fifth check integrates a Poisson background with background subtraction + * switched off, and checks that the result is what would be expected. */ +static void fifth_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 pcount = 0; + + for ( i=0; i<n_trials; i++ ) { + + double intensity, bg, max, sigma; + double fsp, ssp; + + 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(10.0); + if ( (fs-64)*(fs-64) + (ss-64)*(ss-64) <= 10*10 ) { + pcount++; + } + } + } + integrate_peak(image, 64, 64, &fsp, &ssp, &intensity, + &bg, &max, &sigma, 0, 1, 0); + + mean_intensity += intensity; + mean_bg += bg; + mean_max += max; + mean_sigma += sigma; + + } + mean_intensity /= n_trials; + mean_bg /= n_trials; + mean_max /= n_trials; + mean_sigma /= n_trials; + pcount /= n_trials; + + STATUS(" Fifth check (mean values): intensity = %.2f, bg = %.2f," + " max = %.2f, sigma = %.2f\n", + mean_intensity, mean_bg, mean_max, mean_sigma); + double s = pcount*10.0; + if ( fabs(mean_intensity - s) > 1.0 ) { + ERROR("Mean intensity should be close to %f.\n", pcount*10.0); + *fail = 1; + } + if ( fabs(mean_bg-10.0) > 0.3 ) { + ERROR("Mean background should be close to ten.\n"); + *fail = 1; + } +} + + +/* The sixth check is like the fourth check, except that the background + * subtraction is switched off */ +static void sixth_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 pcount = 0; + int npcount = 0; + + for ( i=0; i<n_trials; i++ ) { + + double intensity, bg, max, sigma; + double fsp, ssp; + + for ( fs=0; fs<image->width; fs++ ) { + for ( ss=0; ss<image->height; ss++ ) { + int idx = fs+image->width*ss; + double r = (fs-64)*(fs-64) + (ss-64)*(ss-64); + image->data[idx] = poisson_noise(10.0); + if ( r < 9*9 ) { + image->data[idx] += 1000.0; + pcount++; + } else if ( r <= 10*10 ) { + npcount++; + } + } + } + integrate_peak(image, 64, 64, &fsp, &ssp, &intensity, + &bg, &max, &sigma, 0, 1, 0); + + mean_intensity += intensity; + mean_bg += bg; + mean_max += max; + mean_sigma += sigma; + + } + mean_intensity /= n_trials; + mean_bg /= n_trials; + mean_max /= n_trials; + mean_sigma /= n_trials; + pcount /= n_trials; + npcount /= n_trials; + + STATUS(" Sixth check (mean values): intensity = %.2f, bg = %.2f," + " max = %.2f, sigma = %.2f\n", + mean_intensity, mean_bg, mean_max, mean_sigma); + + double s = pcount*1010.0 + npcount*10.0; + if ( fabs(mean_intensity - s) > 4000.0 ) { + ERROR("Mean intensity should be close to %f.\n", s); + *fail = 1; + } + if ( fabs(mean_bg-10.0) > 0.3 ) { + ERROR("Mean background should be close to ten.\n"); + *fail = 1; + } + + STATUS(" (Absolute value of mean sigma not (yet) tested)\n"); +} + + int main(int argc, char *argv[]) { struct image image; double fsp, ssp, intensity, bg, max, sigma; int fs, ss; + FILE *fh; + unsigned int seed; + int fail = 0; + const int n_trials = 1000; + + 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; @@ -65,16 +315,16 @@ int main(int argc, char *argv[]) /* First check: no intensity -> zero intensity and bg */ integrate_peak(&image, 64, 64, &fsp, &ssp, &intensity, - &bg, &max, &sigma, 0, 1); - STATUS(" First check: intensity = %.2f, bg = %.2f, max = %.2f," + &bg, &max, &sigma, 0, 1, 1); + STATUS(" First check: intensity = %.2f, bg = %.2f, max = %.2f," " sigma = %.2f\n", intensity, bg, max, sigma); if ( intensity != 0.0 ) { ERROR("Intensity should be zero.\n"); - return 1; + fail = 1; } if ( bg != 0.0 ) { ERROR("Background should be zero.\n"); - return 1; + fail = 1; } /* Second check: uniform peak gives correct value and no bg */ @@ -85,81 +335,44 @@ int main(int argc, char *argv[]) } } integrate_peak(&image, 64, 64, &fsp, &ssp, &intensity, - &bg, &max, &sigma, 0, 1); - STATUS("Second check: intensity = %.2f, bg = %.2f, max = %.2f," + &bg, &max, &sigma, 0, 1, 1); + STATUS(" Second check: intensity = %.2f, bg = %.2f, max = %.2f," " sigma = %.2f\n", intensity, bg, max, sigma); if ( fabs(intensity - M_PI*9.0*9.0*1000.0) > 4000.0 ) { ERROR("Intensity should be close to 1000*pi*integr_r^2\n"); - return 1; + fail = 1; } if ( bg != 0.0 ) { ERROR("Background should be zero.\n"); - return 1; + fail = 1; } if ( max != 1000.0 ) { ERROR("Max should be 1000.\n"); - return 1; + fail = 1; } if ( sigma != 0.0 ) { ERROR("Sigma should be zero.\n"); - return 1; + fail = 1; } /* Third check: Poisson background should get mostly subtracted */ - for ( fs=0; fs<image.width; fs++ ) { - for ( ss=0; ss<image.height; ss++ ) { - image.data[fs+image.width*ss] = poisson_noise(10.0); - } - } - integrate_peak(&image, 64, 64, &fsp, &ssp, &intensity, - &bg, &max, &sigma, 0, 1); - STATUS(" Third check: intensity = %.2f, bg = %.2f, max = %.2f," - " sigma = %.2f\n", intensity, bg, max, sigma); - if ( fabs(intensity) > 100.0 ) { - ERROR("Intensity should be close to zero.\n"); - return 1; - } - if ( fabs(bg-10.0) > 10.0 ) { - ERROR("Background should be close to ten.\n"); - return 1; - } - if ( fabs(intensity) > sigma ) { - ERROR("Intensity should be less than sigma.\n"); - return 1; - } - if ( fabs(sigma-71.0) > 10.0 ) { - ERROR("Sigma should be close to 71.\n"); - return 1; - } + third_integration_check(&image, n_trials, &fail); /* Fourth check: peak on Poisson background */ - 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; - } - } - integrate_peak(&image, 64, 64, &fsp, &ssp, &intensity, - &bg, &max, &sigma, 0, 1); - STATUS("Fourth check: intensity = %.2f, bg = %.2f, max = %.2f," - " sigma = %.2f\n", intensity, bg, max, sigma); - if ( fabs(intensity - M_PI*9.0*9.0*1000.0) > 10000.0 ) { - ERROR("Intensity should be close to 1000*pi*integr_r^2.\n"); - return 1; - } - if ( fabs(bg - 10.0) > 1.0 ) { - ERROR("Background should be close to 10.\n"); - return 1; - } - if ( fabs(sigma-71.0) > 10.0 ) { - ERROR("Sigma should be close to 71.\n"); - return 1; - } + fourth_integration_check(&image, n_trials, &fail); + + /* Fifth check: Like third check but without background subtraction */ + fifth_integration_check(&image, n_trials, &fail); + + /* Sixth check: Like fourth check but without background subtraction */ + sixth_integration_check(&image, n_trials, &fail); free(image.beam); free(image.det->panels); free(image.det); free(image.data); + if ( fail ) return 1; + return 0; } |