Commit f69a7404 authored by Andrew Martin's avatar Andrew Martin Committed by Thomas White
Browse files

Added background subtraction and sigma

parent 45177822
......@@ -285,7 +285,7 @@ double integrate_all(struct image *image, RefList *reflections)
get_detector_pos(refl, &xp, &yp);
if ( integrate_peak(image, xp, yp, &x, &y,
&intensity, NULL, NULL, 0, 0) ) continue;
&intensity, NULL, NULL, NULL, 0, 0) ) continue;
itot += intensity;
}
......
......@@ -4,6 +4,7 @@
* Peak search and other image analysis
*
* (c) 2006-2011 Thomas White <taw@physics.org>
* 2011 Andrew Martin
*
* Part of CrystFEL - crystallography with a FEL
*
......@@ -142,7 +143,7 @@ static int cull_peaks(struct image *image)
* 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 *pbg, double *pmax, double *sigma,
int do_polar, int centroid)
{
signed int fs, ss;
......@@ -155,6 +156,10 @@ int integrate_peak(struct image *image, int cfs, int css,
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;
p = find_panel(image->det, cfs, css);
if ( p == NULL ) return 1;
......@@ -203,16 +208,6 @@ int integrate_peak(struct image *image, int cfs, int css,
val = image->data[idx];
/* Inner mask */
if ( fs*fs + ss*ss > lim_sq ) {
/* Estimate noise from this region */
noise += fabs(val);
noise_counts++;
continue;
}
if ( val > max ) max = val;
if ( do_polar ) {
tt = get_tt(image, fs+cfs, ss+css);
......@@ -226,27 +221,54 @@ int integrate_peak(struct image *image, int cfs, int css,
}
total += val;
if ( val > max ) max = val;
/* Inner mask */
if ( fs*fs + ss*ss > lim_sq ) {
/* Estimate noise from this region */
noise += val; //fabs(val);
noise_counts++;
/* Estimate the standard deviation of the noise */
noise_meansq += fabs(val)*fabs(val) ;
continue;
}
pixel_counts ++;
total += val;
fsct += val*(cfs+fs);
ssct += val*(css+ss);
}
}
noise_mean = noise / noise_counts;
/* The centroid is excitingly undefined if there is no intensity */
if ( centroid && (total != 0) ) {
*pfs = (double)fsct / total;
*pss = (double)ssct / total;
*intensity = total;
*intensity = total - pixel_counts*noise_mean;
} else {
*pfs = (double)cfs;
*pss = (double)css;
*intensity = total;
*intensity = total - pixel_counts*noise_mean;
}
if ( in_bad_region(image->det, *pfs, *pss) ) return 1;
if ( sigma != NULL ) {
/*
* 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);
/* printf(" counts %d %d \n", noise_counts, pixel_counts);
printf(" intensity, bg, diff, %f, %f, %f \n", total, pixel_counts*noise_mean, total - pixel_counts*noise_mean);
printf(" sigma = %f \n", *sigma); */
}
if ( pbg != NULL ) {
*pbg = (noise / noise_counts);
}
......@@ -362,7 +384,7 @@ static void search_peaks_in_panel(struct image *image, float threshold,
* Don't bother doing polarisation/SA correction, because the
* intensity of this peak is only an estimate at this stage. */
r = integrate_peak(image, mask_fs, mask_ss,
&f_fs, &f_ss, &intensity, NULL, NULL, 0, 1);
&f_fs, &f_ss, &intensity, NULL, NULL, NULL, 0, 1);
if ( r ) {
/* Bad region - don't detect peak */
nrej_bad++;
......@@ -587,6 +609,7 @@ void integrate_reflections(struct image *image, int polar, int use_closer)
double d;
int idx;
double bg, max;
double sigma;
struct panel *p;
double pfs, pss;
int r;
......@@ -616,11 +639,12 @@ void integrate_reflections(struct image *image, int polar, int use_closer)
}
r = integrate_peak(image, pfs, pss, &fs, &ss,
&intensity, &bg, &max, polar, 0);
&intensity, &bg, &max, &sigma, polar, 0);
/* Record intensity and set redundancy to 1 on success */
if ( r == 0 ) {
set_int(refl, intensity);
set_esd_intensity(refl,sigma);
set_redundancy(refl, 1);
}
......
......@@ -35,7 +35,7 @@ extern RefList *find_projected_peaks(struct image *image, UnitCell *cell,
extern int integrate_peak(struct image *image, int xp, int yp,
double *xc, double *yc, double *intensity,
double *pbg, double *pmax,
double *pbg, double *pmax, double *sigma,
int do_polar, int centroid);
#endif /* PEAKS_H */
......@@ -4,7 +4,7 @@
* Assemble and process FEL Bragg intensities
*
* (c) 2006-2011 Thomas White <taw@physics.org>
*
* 2011 Andrew Martin
* Part of CrystFEL - crystallography with a FEL
*
*/
......@@ -171,8 +171,12 @@ static void merge_pattern(RefList *model, RefList *new, int max_only,
} else if ( pass == 2 ) {
double dev = get_sum_squared_dev(model_version);
set_sum_squared_dev(model_version,
dev + pow(intensity-model_int, 2.0));
double refl_esd = get_esd_intensity(refl);
set_sum_squared_dev(model_version,
// dev + pow(refl_esd,2.0) );
dev + pow(intensity, 2.0) );
// dev + pow(intensity-model_int, 2.0));
if ( hist_vals != NULL ) {
int p = *hist_n;
......@@ -188,6 +192,8 @@ static void merge_pattern(RefList *model, RefList *new, int max_only,
}
}
}
......@@ -365,11 +371,19 @@ static void merge_all(FILE *fh, RefList *model,
refl = next_refl(refl, iter) ) {
double sum_squared_dev = get_sum_squared_dev(refl);
double intensity = get_intensity(refl);
int red = get_redundancy(refl);
set_esd_intensity(refl,
sqrt(sum_squared_dev)/(double)red);
int h, k, l;
get_indices(refl,&h,&k,&l);
/*
set_esd_intensity(refl,
sqrt(sum_squared_dev)/(double)red);
*/
set_esd_intensity(refl,
sqrt((sum_squared_dev/(double)red) - pow(intensity,2.0) )/(double)red);
}
}
......
......@@ -27,7 +27,7 @@
int main(int argc, char *argv[])
{
struct image image;
double fsp, ssp, intensity, bg, max;
double fsp, ssp, intensity, bg, max, sigma;
int fs, ss;
image.data = malloc(128*128*sizeof(float));
......@@ -63,7 +63,7 @@ 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, 0, 1);
integrate_peak(&image, 64, 64, &fsp, &ssp, &intensity, &bg, &max, &sigma, 0, 1);
STATUS(" First check: intensity = %.2f bg = %.2f max = %.2f\n",
intensity, bg, max);
if ( intensity != 0.0 ) {
......@@ -82,7 +82,7 @@ int main(int argc, char *argv[])
image.data[fs+image.width*ss] = 1000.0;
}
}
integrate_peak(&image, 64, 64, &fsp, &ssp, &intensity, &bg, &max, 0, 1);
integrate_peak(&image, 64, 64, &fsp, &ssp, &intensity, &bg, &max, &sigma, 0, 1);
STATUS("Second check: intensity = %.2f, bg = %.2f, max = %.2f\n",
intensity, bg, max);
if ( fabs(intensity - M_PI*9.0*9.0*1000.0) > 4000.0 ) {
......@@ -104,7 +104,7 @@ int main(int argc, char *argv[])
image.data[fs+image.width*ss] = poisson_noise(10.0) - 10.0;
}
}
integrate_peak(&image, 64, 64, &fsp, &ssp, &intensity, &bg, &max, 0, 1);
integrate_peak(&image, 64, 64, &fsp, &ssp, &intensity, &bg, &max, &sigma, 0, 1);
STATUS(" Third check: intensity = %.2f, bg = %.2f, max = %.2f\n",
intensity, bg, max);
if ( fabs(intensity) > 100.0 ) {
......@@ -122,7 +122,7 @@ int main(int argc, char *argv[])
image.data[fs+image.width*ss] = 1000.0;
}
}
integrate_peak(&image, 64, 64, &fsp, &ssp, &intensity, &bg, &max, 0, 1);
integrate_peak(&image, 64, 64, &fsp, &ssp, &intensity, &bg, &max, &sigma, 0, 1);
STATUS("Fourth check: intensity = %.2f, bg = %.2f, max = %.2f\n",
intensity, bg, max);
if ( fabs(intensity - M_PI*9.0*9.0*1000.0) > 4000.0 ) {
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment