aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--libcrystfel/src/peaks.c140
-rw-r--r--libcrystfel/src/peaks.h8
-rw-r--r--tests/integration_check.c40
3 files changed, 89 insertions, 99 deletions
diff --git a/libcrystfel/src/peaks.c b/libcrystfel/src/peaks.c
index 66907298..e1a828e5 100644
--- a/libcrystfel/src/peaks.c
+++ b/libcrystfel/src/peaks.c
@@ -151,23 +151,20 @@ static int cull_peaks(struct image *image)
/* Returns non-zero if peak has been vetoed.
* i.e. don't use result if return value is not zero. */
int integrate_peak(struct image *image, int cfs, int css,
- double *pfs, double *pss, double *intensity,
- double *pbg, double *pmax, double *sigma,
- int centroid, int bgsub)
+ double *pfs, double *pss, double *intensity, double *sigma)
{
signed int fs, ss;
double lim, out_lim, mid_lim;
double lim_sq, out_lim_sq, mid_lim_sq;
- double total = 0.0;
- double fsct = 0.0;
- double ssct = 0.0;
- double noise = 0.0;
- int noise_counts = 0;
- double max = 0.0;
- struct panel *p = NULL;
- int pixel_counts = 0;
- double noise_mean = 0.0;
- double noise_meansq = 0.0;
+ double pk_total;
+ int pk_counts;
+ double fsct, ssct;
+ double bg_tot = 0.0;
+ int bg_counts = 0;
+ struct panel *p;
+ double bg_mean, bg_var;
+ double bg_tot_sq = 0.0;
+ double var;
struct beam_params *beam;
double aduph;
@@ -189,6 +186,7 @@ int integrate_peak(struct image *image, int cfs, int css,
mid_lim_sq = pow(mid_lim, 2.0);
out_lim_sq = pow(out_lim, 2.0);
+ /* Estimate the background */
for ( fs=-out_lim; fs<+out_lim; fs++ ) {
for ( ss=-out_lim; ss<+out_lim; ss++ ) {
@@ -199,11 +197,9 @@ int integrate_peak(struct image *image, int cfs, int css,
struct panel *p2;
int idx;
- /* Outer mask radius */
+ /* Restrict to annulus */
if ( fs*fs + ss*ss > out_lim_sq ) continue;
-
- if ( ((fs+cfs)>=image->width) || ((fs+cfs)<0) ) continue;
- if ( ((ss+css)>=image->height) || ((ss+css)<0) ) continue;
+ if ( fs*fs + ss*ss < mid_lim_sq ) continue;
/* Strayed off one panel? */
p2 = find_panel(image->det, fs+cfs, ss+css);
@@ -227,65 +223,80 @@ int integrate_peak(struct image *image, int cfs, int css,
val = image->data[idx];
- if ( val > max ) max = val;
+ bg_tot += val;
+ bg_tot_sq += pow(val, 2.0);
+ bg_counts++;
- /* If outside inner mask, estimate noise from this region */
- if ( fs*fs + ss*ss > mid_lim_sq ) {
+ }
+ }
- /* Noise
- * noise and noise_meansq are both in photons (^2) */
- noise += val / aduph;
- noise_counts++;
- noise_meansq += pow(val/aduph, 2.0);
+ if ( bg_counts == 0 ) return 1;
+ bg_mean = bg_tot / bg_counts;
+ bg_var = (bg_tot_sq/bg_counts) - pow(bg_mean, 2.0);
- } else if ( fs*fs + ss*ss < lim_sq ) {
+ /* Measure the peak */
+ pk_total = 0.0;
+ pk_counts = 0;
+ fsct = 0.0; ssct = 0.0;
+ for ( fs=-lim; fs<+lim; fs++ ) {
+ for ( ss=-lim; ss<+lim; ss++ ) {
- /* Peak */
- pixel_counts++;
- total += val;
- fsct += val*(cfs+fs);
- ssct += val*(css+ss);
+ double val;
+ double tt = 0.0;
+ double phi, pa, pb, pol;
+ uint16_t flags;
+ struct panel *p2;
+ int idx;
- }
+ /* Inner mask radius */
+ if ( fs*fs + ss*ss > lim_sq ) continue;
- }
- }
+ /* Strayed off one panel? */
+ p2 = find_panel(image->det, fs+cfs, ss+css);
+ if ( p2 != p ) return 1;
- noise_mean = noise / noise_counts; /* photons */
+ idx = fs+cfs+image->width*(ss+css);
- /* The centroid is excitingly undefined if there is no intensity */
- if ( centroid && (total != 0) ) {
- *pfs = ((double)fsct / total) + 0.5;
- *pss = ((double)ssct / total) + 0.5;
- } else {
- *pfs = (double)cfs + 0.5;
- *pss = (double)css + 0.5;
- }
- if ( bgsub ) {
- *intensity = total - aduph * pixel_counts*noise_mean; /* ADU */
- } else {
- *intensity = total; /* ADU */
- }
+ /* Veto this peak if we tried to integrate in a bad region */
+ if ( image->flags != NULL ) {
- if ( in_bad_region(image->det, *pfs, *pss) ) return 1;
+ flags = image->flags[idx];
- if ( sigma != NULL ) {
+ /* It must have all the "good" bits to be valid */
+ if ( !((flags & image->det->mask_good)
+ == image->det->mask_good) ) return 1;
- /* First term is standard deviation of background per pixel
- * sqrt(pixel_counts) - increase of error for integrated value
- * sqrt(2) - increase of error for background subtraction */
- *sigma = sqrt(noise_meansq/noise_counts-(noise_mean*noise_mean))
- * sqrt(2.0*pixel_counts) * aduph;
+ /* If it has any of the "bad" bits, reject */
+ if ( flags & image->det->mask_bad ) return 1;
- }
+ }
+
+ val = image->data[idx] - bg_mean;
+
+ /* At least half a photon? */
+ if ( val < aduph / 2.0 ) continue;
+
+ pk_counts++;
+ pk_total += val;
+
+ fsct += val*(cfs+fs);
+ ssct += val*(css+ss);
- if ( pbg != NULL ) {
- *pbg = aduph * (noise / noise_counts);
}
- if ( pmax != NULL ) {
- *pmax = max;
}
+ if ( pk_counts == 0 ) return 1;
+
+ *pfs = ((double)fsct / pk_total) + 0.5;
+ *pss = ((double)ssct / pk_total) + 0.5;
+
+ var = pk_counts * bg_var;
+ var += aduph * pk_total;
+ if ( var < 0.0 ) return 1;
+
+ if ( intensity != NULL ) *intensity = pk_total;
+ if ( sigma != NULL ) *sigma = sqrt(var);
+
return 0;
}
@@ -396,12 +407,9 @@ static void search_peaks_in_panel(struct image *image, float threshold,
assert(mask_fs >= p->min_fs);
assert(mask_ss >= p->min_ss);
- /* Centroid peak and get better coordinates.
- * Don't bother doing polarisation/SA correction, because the
- * intensity of this peak is only an estimate at this stage. */
+ /* Centroid peak and get better coordinates. */
r = integrate_peak(image, mask_fs, mask_ss,
- &f_fs, &f_ss, &intensity,
- &pbg, &pmax, &sigma, 1, 1);
+ &f_fs, &f_ss, &intensity, &sigma);
if ( r ) {
/* Bad region - don't detect peak */
@@ -654,7 +662,7 @@ void integrate_reflections(struct image *image, int use_closer, int bgsub,
}
r = integrate_peak(image, pfs, pss, &fs, &ss,
- &intensity, &bg, &max, &sigma, 0, bgsub);
+ &intensity, &sigma);
/* Record intensity and set redundancy to 1 on success */
if ( r == 0 ) {
diff --git a/libcrystfel/src/peaks.h b/libcrystfel/src/peaks.h
index 8b64a081..e5e31b91 100644
--- a/libcrystfel/src/peaks.h
+++ b/libcrystfel/src/peaks.h
@@ -45,10 +45,10 @@ extern double peak_lattice_agreement(struct image *image, UnitCell *cell,
extern int peak_sanity_check(struct image *image);
/* Exported so it can be poked by integration_check */
-extern int integrate_peak(struct image *image, int cfs, int css,
- double *pfs, double *pss, double *intensity,
- double *pbg, double *pmax, double *sigma,
- int centroid, int bgsub);
+extern int integrate_peak(struct image *image,
+ int cfs, int css,
+ double *pfs, double *pss,
+ double *intensity, double *sigma);
extern void estimate_resolution(RefList *list, UnitCell *cell,
double *min, double *max);
diff --git a/tests/integration_check.c b/tests/integration_check.c
index aa9e2705..1e439611 100644
--- a/tests/integration_check.c
+++ b/tests/integration_check.c
@@ -59,8 +59,7 @@ static void third_integration_check(struct image *image, int n_trials,
image->data[fs+image->width*ss] = poisson_noise(10.0);
}
}
- integrate_peak(image, 64, 64, &fsp, &ssp, &intensity,
- &bg, &max, &sigma, 1, 1);
+ integrate_peak(image, 64, 64, &fsp, &ssp, &intensity, &sigma);
mean_intensity += intensity;
mean_bg += bg;
@@ -122,8 +121,7 @@ static void fourth_integration_check(struct image *image, int n_trials,
}
}
}
- integrate_peak(image, 64, 64, &fsp, &ssp, &intensity,
- &bg, &max, &sigma, 1, 1);
+ integrate_peak(image, 64, 64, &fsp, &ssp, &intensity, &sigma);
mean_intensity += intensity;
mean_bg += bg;
@@ -182,8 +180,7 @@ static void fifth_integration_check(struct image *image, int n_trials,
}
}
}
- integrate_peak(image, 64, 64, &fsp, &ssp, &intensity,
- &bg, &max, &sigma, 1, 0);
+ integrate_peak(image, 64, 64, &fsp, &ssp, &intensity, &sigma);
mean_intensity += intensity;
mean_bg += bg;
@@ -244,8 +241,7 @@ static void sixth_integration_check(struct image *image, int n_trials,
}
}
}
- integrate_peak(image, 64, 64, &fsp, &ssp, &intensity,
- &bg, &max, &sigma, 1, 0);
+ integrate_peak(image, 64, 64, &fsp, &ssp, &intensity, &sigma);
mean_intensity += intensity;
mean_bg += bg;
@@ -281,7 +277,7 @@ static void sixth_integration_check(struct image *image, int n_trials,
int main(int argc, char *argv[])
{
struct image image;
- double fsp, ssp, intensity, bg, max, sigma;
+ double fsp, ssp, intensity, sigma;
int fs, ss;
FILE *fh;
unsigned int seed;
@@ -326,18 +322,13 @@ int main(int argc, char *argv[])
image.beam->adu_per_photon = 100.0;
/* First check: no intensity -> zero intensity and bg */
- integrate_peak(&image, 64, 64, &fsp, &ssp, &intensity,
- &bg, &max, &sigma, 1, 1);
- STATUS(" First check: intensity = %.2f, bg = %.2f, max = %.2f,"
- " sigma = %.2f\n", intensity, bg, max, sigma);
+ integrate_peak(&image, 64, 64, &fsp, &ssp, &intensity, &sigma);
+ STATUS(" First check: intensity = %.2f, sigma = %.2f\n",
+ intensity, sigma);
if ( intensity != 0.0 ) {
ERROR("Intensity should be zero.\n");
fail = 1;
}
- if ( bg != 0.0 ) {
- ERROR("Background should be zero.\n");
- fail = 1;
- }
/* Second check: uniform peak gives correct value and no bg */
for ( fs=0; fs<image.width; fs++ ) {
@@ -346,22 +337,13 @@ int main(int argc, char *argv[])
image.data[fs+image.width*ss] = 1000.0;
}
}
- integrate_peak(&image, 64, 64, &fsp, &ssp, &intensity,
- &bg, &max, &sigma, 1, 1);
- STATUS(" Second check: intensity = %.2f, bg = %.2f, max = %.2f,"
- " sigma = %.2f\n", intensity, bg, max, sigma);
+ integrate_peak(&image, 64, 64, &fsp, &ssp, &intensity, &sigma);
+ STATUS(" Second check: intensity = %.2f, sigma = %.2f\n",
+ intensity, 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");
fail = 1;
}
- if ( bg != 0.0 ) {
- ERROR("Background should be zero.\n");
- fail = 1;
- }
- if ( max != 1000.0 ) {
- ERROR("Max should be 1000.\n");
- fail = 1;
- }
if ( sigma != 0.0 ) {
ERROR("Sigma should be zero.\n");
fail = 1;
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-05 23:34:02 +0000