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() --- libcrystfel/src/integration.c | 350 +++++++++++++++++++++++++++++++----------- 1 file changed, 261 insertions(+), 89 deletions(-) (limited to 'libcrystfel/src/integration.c') 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])); + + } } -- cgit v1.2.3