aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--Makefile.am6
-rw-r--r--libcrystfel/src/integration.c350
-rw-r--r--tests/.gitignore1
-rw-r--r--tests/integration_check.c185
-rw-r--r--tests/ring_check.c314
5 files changed, 584 insertions, 272 deletions
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; p<ic->w; p++ ) {
for ( q=0; q<ic->w; 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; p<ic->w; p++ ) {
for ( q=0; q<ic->w; 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<nshells; i++ ) {
- one_over_d = resolution(cell, h, k, l);
- if ( one_over_d > 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; i<nshells; i++ ) {
+ sh_nref[i] = 0;
+ sh_snr[i] = 0.0;
+ }
- free_intcontext(&ic);
+ for ( refl = first_refl(list, &iter);
+ refl != NULL;
+ refl = next_refl(refl, iter) )
+ {
+ signed int h, k, l;
+ double res, intensity, sigma, snr;
+ int bin, j;
- crystal_set_num_saturated_reflections(cr, n_saturated);
- crystal_set_resolution_limit(cr, limit);
- crystal_set_reflections(cr, list);
+ intensity = get_intensity(refl);
+ sigma = get_esd_intensity(refl);
+ snr = intensity / sigma;
+
+ if ( snr < min_snr ) continue;
+
+ get_indices(refl, &h, &k, &l);
+ res = 2.0*resolution(cell, h, k, l);
+
+ bin = -1;
+ for ( j=0; j<nshells; j++ ) {
+ if ( (res>rmins[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; i<nshells; i++ ) {
+ double cen = (rmins[i]+rmaxs[i])/2.0;
+ double val = sh_nref[i];
+ fprintf(fh, "%10.3f %.2f %.4f\n", cen/1e9, 10e9/cen, val);
+ }
+ fclose(fh);
+
+ free(rmins);
+ free(rmaxs);
+ free(sh_nref);
+ free(sh_snr);
}
@@ -1627,24 +1794,29 @@ void integrate_all(struct image *image, IntegrationMethod meth,
switch ( meth & INTEGRATION_METHOD_MASK ) {
case INTEGRATION_NONE :
- return;
+ break;
case INTEGRATION_RINGS :
integrate_rings(meth, image->crystals[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 <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;
+}
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-18 12:02:40 +0000