From 53cf97428ebfc7f53fee5970342a836bf45017ec Mon Sep 17 00:00:00 2001 From: Thomas White Date: Thu, 7 Jun 2012 14:49:28 +0200 Subject: Use max_adu only for final integration Saturated pixels seem to be OK during initial peak search --- libcrystfel/src/peaks.c | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) (limited to 'libcrystfel/src/peaks.c') diff --git a/libcrystfel/src/peaks.c b/libcrystfel/src/peaks.c index 477d9345..670e02bb 100644 --- a/libcrystfel/src/peaks.c +++ b/libcrystfel/src/peaks.c @@ -151,7 +151,8 @@ static int cull_peaks(struct image *image) static int integrate_peak(struct image *image, int cfs, int css, double *pfs, double *pss, double *intensity, double *sigma, - double ir_inn, double ir_mid, double ir_out) + double ir_inn, double ir_mid, double ir_out, + int use_max_adu) { signed int fs, ss; double lim_sq, out_lim_sq, mid_lim_sq; @@ -212,7 +213,7 @@ static int integrate_peak(struct image *image, int cfs, int css, val = image->data[idx]; /* Veto peak if it contains saturation in bg region */ - if ( val > p->max_adu ) return 1; + if ( use_max_adu && (val > p->max_adu) ) return 1; bg_tot += val; bg_tot_sq += pow(val, 2.0); @@ -263,7 +264,7 @@ static int integrate_peak(struct image *image, int cfs, int css, val = image->data[idx] - bg_mean; /* Veto peak if it contains saturation */ - if ( image->data[idx] > p->max_adu ) return 1; + if ( use_max_adu && (image->data[idx] > p->max_adu) ) return 1; pk_counts++; pk_total += val; @@ -398,7 +399,7 @@ static void search_peaks_in_panel(struct image *image, float threshold, /* Centroid peak and get better coordinates. */ r = integrate_peak(image, mask_fs, mask_ss, &f_fs, &f_ss, &intensity, &sigma, - ir_inn, ir_mid, ir_out); + ir_inn, ir_mid, ir_out, 0); if ( r ) { /* Bad region - don't detect peak */ @@ -656,7 +657,8 @@ void integrate_reflections(struct image *image, int use_closer, int bgsub, } r = integrate_peak(image, pfs, pss, &fs, &ss, - &intensity, &sigma, ir_inn, ir_mid, ir_out); + &intensity, &sigma, ir_inn, ir_mid, ir_out, + 1); /* Record intensity and set redundancy to 1 on success */ if ( r == 0 ) { -- cgit v1.2.3 From b90ddc13d5638aa5e9e4094228f464ce48989eb7 Mon Sep 17 00:00:00 2001 From: Kenneth Beyerlein Date: Wed, 13 Jun 2012 17:32:07 +0200 Subject: Mask peaks in the background calculation of integrate_peak. --- libcrystfel/src/peaks.c | 57 +++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 57 insertions(+) (limited to 'libcrystfel/src/peaks.c') diff --git a/libcrystfel/src/peaks.c b/libcrystfel/src/peaks.c index 477d9345..5987f8ff 100644 --- a/libcrystfel/src/peaks.c +++ b/libcrystfel/src/peaks.c @@ -145,6 +145,52 @@ static int cull_peaks(struct image *image) return nelim; } +static void make_BgMask(struct image *image, unsigned char *mask, + double ir_out, int cfs, int css, double ir_in, + struct panel *p) +{ + signed int fs, ss; + RefList *list = image->reflections; + Reflection *refl; + RefListIterator *iter; + double pk2_fs, pk2_ss; + double d_fs, d_ss, distSq; + double limSq = ir_in*ir_in; + struct panel *p2; + + if (num_reflections(list)==0) return; + + /* Loop over all reflections */ + for ( refl = first_refl(list, &iter); + refl != NULL; + refl = next_refl(refl, iter) ) + { + + get_detector_pos(refl, &pk2_fs, &pk2_ss); + + /*Determine if reflection is in the same panel */ + p2 = find_panel(image->det, pk2_fs, pk2_ss); + if ( p2 != p ) continue; + + /* If other peak area overlaps larger bg area, color in the mask */ + if (fabs(pk2_fs-cfs)det, cfs, css); if ( p == NULL ) return 1; if ( p->no_index ) return 1; + /* Call function to block regions where there is expected to be a peak */ + make_BgMask(image, bgPkMask, ir_out, cfs, css, ir_inn, p); + aduph = p->adu_per_eV * ph_lambda_to_eV(image->lambda); lim_sq = pow(ir_inn, 2.0); @@ -189,6 +241,9 @@ static int integrate_peak(struct image *image, int cfs, int css, if ( fs*fs + ss*ss > out_lim_sq ) continue; if ( fs*fs + ss*ss < mid_lim_sq ) continue; + /* Check if there is a peak in the background region */ + if (bgPkMask[(int)(fs+ir_out+2*ir_out*(ss+ir_out))]) continue; + /* Strayed off one panel? */ p2 = find_panel(image->det, fs+cfs, ss+css); if ( p2 != p ) return 1; @@ -286,6 +341,8 @@ static int integrate_peak(struct image *image, int cfs, int css, if ( intensity != NULL ) *intensity = pk_total; if ( sigma != NULL ) *sigma = sqrt(var); + free(bgPkMask); + return 0; } -- cgit v1.2.3 From 469518a4cc84348e859e226a8b4d9ef7d7a1479b Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 20 Jun 2012 16:00:17 +0200 Subject: Fussiness --- libcrystfel/src/peaks.c | 75 +++++++++++++++++++++++++++---------------------- 1 file changed, 42 insertions(+), 33 deletions(-) (limited to 'libcrystfel/src/peaks.c') diff --git a/libcrystfel/src/peaks.c b/libcrystfel/src/peaks.c index 35cc1f96..e541ec47 100644 --- a/libcrystfel/src/peaks.c +++ b/libcrystfel/src/peaks.c @@ -145,53 +145,62 @@ static int cull_peaks(struct image *image) return nelim; } -static void make_BgMask(struct image *image, unsigned char *mask, - double ir_out, int cfs, int css, double ir_in, - struct panel *p) +static unsigned char *make_BgMask(struct image *image, struct panel *p, + double ir_out, int cfs, int css, double ir_in) { - signed int fs, ss; - RefList *list = image->reflections; Reflection *refl; RefListIterator *iter; - double pk2_fs, pk2_ss; - double d_fs, d_ss, distSq; - double limSq = ir_in*ir_in; - struct panel *p2; + unsigned char *mask; - if (num_reflections(list)==0) return; + mask = calloc(4*ir_out*ir_out, sizeof(unsigned char)); + if ( mask == NULL ) return NULL; /* Loop over all reflections */ - for ( refl = first_refl(list, &iter); - refl != NULL; - refl = next_refl(refl, iter) ) + for ( refl = first_refl(image->reflections, &iter); + refl != NULL; + refl = next_refl(refl, iter) ) { + struct panel *p2; + double pk2_fs, pk2_ss; + + get_detector_pos(refl, &pk2_fs, &pk2_ss); + + /* Determine if reflection is in the same panel */ + p2 = find_panel(image->det, pk2_fs, pk2_ss); + if ( p2 != p ) continue; + + /* If other peak area overlaps larger bg area, set mask */ + if ( (fabs(pk2_fs-cfs)det, pk2_fs, pk2_ss); - if ( p2 != p ) continue; + double d_fs, d_ss, distSq; - /* If other peak area overlaps larger bg area, color in the mask */ - if (fabs(pk2_fs-cfs)det, cfs, css); if ( p == NULL ) return 1; if ( p->no_index ) return 1; - /* Call function to block regions where there is expected to be a peak */ - make_BgMask(image, bgPkMask, ir_out, cfs, css, ir_inn, p); + /* Determine regions where there is expected to be a peak */ + bgPkMask = make_BgMask(image, p, ir_out, cfs, css, ir_inn); + if ( bgPkMask == NULL ) return 1; aduph = p->adu_per_eV * ph_lambda_to_eV(image->lambda); -- cgit v1.2.3 From 64a6e175e924f57aa10bbd67ec45df1403a73276 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 20 Jun 2012 16:00:25 +0200 Subject: Revert "Remove deprecated SNR cutoff" This reverts commit ae51e75490daf47e2deefe83e72a1f5c092fa023. --- libcrystfel/src/peaks.c | 3 +++ 1 file changed, 3 insertions(+) (limited to 'libcrystfel/src/peaks.c') diff --git a/libcrystfel/src/peaks.c b/libcrystfel/src/peaks.c index e541ec47..6a7cca34 100644 --- a/libcrystfel/src/peaks.c +++ b/libcrystfel/src/peaks.c @@ -726,6 +726,9 @@ void integrate_reflections(struct image *image, int use_closer, int bgsub, &intensity, &sigma, ir_inn, ir_mid, ir_out, 1); + /* I/sigma(I) cutoff */ + if ( !r && (intensity/sigma < min_snr) ) r = 1; + /* Record intensity and set redundancy to 1 on success */ if ( r == 0 ) { set_intensity(refl, intensity); -- cgit v1.2.3 From eedb9e2a4a4305e89344453ff890c3b2a844b2ee Mon Sep 17 00:00:00 2001 From: Thomas White Date: Thu, 21 Jun 2012 11:25:06 +0200 Subject: Update authorship --- libcrystfel/src/peaks.c | 1 + 1 file changed, 1 insertion(+) (limited to 'libcrystfel/src/peaks.c') diff --git a/libcrystfel/src/peaks.c b/libcrystfel/src/peaks.c index 6a7cca34..3ae7de8e 100644 --- a/libcrystfel/src/peaks.c +++ b/libcrystfel/src/peaks.c @@ -9,6 +9,7 @@ * * Authors: * 2010-2012 Thomas White + * 2012 Kenneth Beyerlein * 2011 Andrew Martin * 2011 Richard Kirian * -- cgit v1.2.3 From 0ce14ba49cd6150e09a9dd32ea3451e3e47f16bf Mon Sep 17 00:00:00 2001 From: Thomas White Date: Thu, 21 Jun 2012 12:20:11 +0200 Subject: More fussiness --- libcrystfel/src/peaks.c | 36 +++++++++++++++++------------------- 1 file changed, 17 insertions(+), 19 deletions(-) (limited to 'libcrystfel/src/peaks.c') diff --git a/libcrystfel/src/peaks.c b/libcrystfel/src/peaks.c index 3ae7de8e..99a12ae5 100644 --- a/libcrystfel/src/peaks.c +++ b/libcrystfel/src/peaks.c @@ -146,14 +146,18 @@ static int cull_peaks(struct image *image) return nelim; } + static unsigned char *make_BgMask(struct image *image, struct panel *p, double ir_out, int cfs, int css, double ir_in) { Reflection *refl; RefListIterator *iter; unsigned char *mask; + int w, h; - mask = calloc(4*ir_out*ir_out, sizeof(unsigned char)); + w = p->max_fs - p->min_fs + 1; + h = p->max_ss - p->min_ss + 1; + mask = calloc(w*h, sizeof(unsigned char)); if ( mask == NULL ) return NULL; /* Loop over all reflections */ @@ -163,6 +167,7 @@ static unsigned char *make_BgMask(struct image *image, struct panel *p, { struct panel *p2; double pk2_fs, pk2_ss; + signed int fs, ss; get_detector_pos(refl, &pk2_fs, &pk2_ss); @@ -170,31 +175,25 @@ static unsigned char *make_BgMask(struct image *image, struct panel *p, p2 = find_panel(image->det, pk2_fs, pk2_ss); if ( p2 != p ) continue; - /* If other peak area overlaps larger bg area, set mask */ - if ( (fabs(pk2_fs-cfs)det, cfs, css); if ( p == NULL ) return 1; if ( p->no_index ) return 1; -- cgit v1.2.3 From f668e3b3beec19cfe0162860a9aea9d934121212 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Thu, 21 Jun 2012 15:28:18 +0200 Subject: More fussiness --- libcrystfel/src/peaks.c | 94 ++++++++++++++++++++++++++++--------------------- 1 file changed, 54 insertions(+), 40 deletions(-) (limited to 'libcrystfel/src/peaks.c') diff --git a/libcrystfel/src/peaks.c b/libcrystfel/src/peaks.c index 99a12ae5..fe836c3f 100644 --- a/libcrystfel/src/peaks.c +++ b/libcrystfel/src/peaks.c @@ -147,17 +147,18 @@ static int cull_peaks(struct image *image) } -static unsigned char *make_BgMask(struct image *image, struct panel *p, - double ir_out, int cfs, int css, double ir_in) +/* cfs, css relative to panel origin */ +static int *make_BgMask(struct image *image, struct panel *p, + double ir_out, int cfs, int css, double ir_in) { Reflection *refl; RefListIterator *iter; - unsigned char *mask; + int *mask; int w, h; w = p->max_fs - p->min_fs + 1; h = p->max_ss - p->min_ss + 1; - mask = calloc(w*h, sizeof(unsigned char)); + mask = calloc(w*h, sizeof(int)); if ( mask == NULL ) return NULL; /* Loop over all reflections */ @@ -167,7 +168,8 @@ static unsigned char *make_BgMask(struct image *image, struct panel *p, { struct panel *p2; double pk2_fs, pk2_ss; - signed int fs, ss; + signed int dfs, dss; + double pk2_cfs, pk2_css; get_detector_pos(refl, &pk2_fs, &pk2_ss); @@ -175,22 +177,27 @@ static unsigned char *make_BgMask(struct image *image, struct panel *p, p2 = find_panel(image->det, pk2_fs, pk2_ss); if ( p2 != p ) continue; - for ( fs=-ir_in; fs<=ir_in; fs++ ) { - for ( ss=-ir_in; ss<=ir_in; ss++ ) { + pk2_cfs = pk2_fs - p->min_fs; + pk2_css = pk2_ss - p->min_ss; - double d_fs, d_ss, distSq; + for ( dfs=-ir_in; dfs<=ir_in; dfs++ ) { + for ( dss=-ir_in; dss<=ir_in; dss++ ) { - d_fs = cfs + fs; - d_ss = css + ss; - distSq = d_fs*d_fs + d_ss*d_ss; + signed int fs, ss; - if ( distSq < ir_in*ir_in ) { + /* In peak region for this peak? */ + if ( dfs*dfs + dss*dss > ir_in*ir_in ) continue; - int idx; - idx = fs+ir_out+2*ir_out*(ss+ir_out); - mask[idx] = 1; + fs = pk2_cfs + dfs; + ss = pk2_css + dss; - } + /* On panel? */ + if ( fs >= w ) continue; + if ( ss >= h ) continue; + if ( fs < 0 ) continue; + if ( ss < 0 ) continue; + + mask[fs + ss*w] = 1; } } @@ -209,7 +216,7 @@ static int integrate_peak(struct image *image, int cfs, int css, double ir_inn, double ir_mid, double ir_out, int use_max_adu) { - signed int fs, ss; + signed int dfs, dss; double lim_sq, out_lim_sq, mid_lim_sq; double pk_total; int pk_counts; @@ -221,14 +228,19 @@ static int integrate_peak(struct image *image, int cfs, int css, double bg_tot_sq = 0.0; double var; double aduph; - unsigned char *bgPkMask; + int *bgPkMask; + int p_cfs, p_css, p_w, p_h; p = find_panel(image->det, cfs, css); if ( p == NULL ) return 1; if ( p->no_index ) return 1; /* Determine regions where there is expected to be a peak */ - bgPkMask = make_BgMask(image, p, ir_out, cfs, css, ir_inn); + p_cfs = cfs - p->min_fs; + p_css = css - p->min_ss; /* Panel-relative coordinates */ + p_w = p->max_fs - p->min_fs + 1; + p_h = p->max_ss - p->min_ss + 1; + bgPkMask = make_BgMask(image, p, ir_out, p_cfs, p_css, ir_inn); if ( bgPkMask == NULL ) return 1; aduph = p->adu_per_eV * ph_lambda_to_eV(image->lambda); @@ -238,26 +250,27 @@ static int integrate_peak(struct image *image, int cfs, int css, out_lim_sq = pow(ir_out, 2.0); /* Estimate the background */ - for ( fs=-ir_out; fs<=+ir_out; fs++ ) { - for ( ss=-ir_out; ss<=+ir_out; ss++ ) { + for ( dfs=-ir_out; dfs<=+ir_out; dfs++ ) { + for ( dss=-ir_out; dss<=+ir_out; dss++ ) { double val; uint16_t flags; - struct panel *p2; int idx; /* Restrict to annulus */ - if ( fs*fs + ss*ss > out_lim_sq ) continue; - if ( fs*fs + ss*ss < mid_lim_sq ) continue; - - /* Check if there is a peak in the background region */ - if (bgPkMask[(int)(fs+ir_out+2*ir_out*(ss+ir_out))]) continue; + if ( dfs*dfs + dss*dss > out_lim_sq ) continue; + if ( dfs*dfs + dss*dss < mid_lim_sq ) continue; /* Strayed off one panel? */ - p2 = find_panel(image->det, fs+cfs, ss+css); - if ( p2 != p ) return 1; + if ( p_cfs+dfs >= p_w ) continue; + if ( p_css+dss >= p_h ) continue; + if ( p_cfs+dfs < 0 ) continue; + if ( p_css+dss < 0 ) continue; + + /* Check if there is a peak in the background region */ + if ( bgPkMask[(p_cfs+dfs) + p_w*(p_css+dss)] ) continue; - idx = fs+cfs+image->width*(ss+css); + idx = dfs+cfs+image->width*(dss+css); /* Veto this peak if we tried to integrate in a bad region */ if ( image->flags != NULL ) { @@ -293,22 +306,23 @@ static int integrate_peak(struct image *image, int cfs, int css, pk_total = 0.0; pk_counts = 0; fsct = 0.0; ssct = 0.0; - for ( fs=-ir_inn; fs<=+ir_inn; fs++ ) { - for ( ss=-ir_inn; ss<=+ir_inn; ss++ ) { + for ( dfs=-ir_inn; dfs<=+ir_inn; dfs++ ) { + for ( dss=-ir_inn; dss<=+ir_inn; dss++ ) { double val; uint16_t flags; - struct panel *p2; int idx; /* Inner mask radius */ - if ( fs*fs + ss*ss > lim_sq ) continue; + if ( dfs*dfs + dss*dss > lim_sq ) continue; /* Strayed off one panel? */ - p2 = find_panel(image->det, fs+cfs, ss+css); - if ( p2 != p ) return 1; + if ( p_cfs+dfs >= p_w ) continue; + if ( p_css+dss >= p_h ) continue; + if ( p_cfs+dfs < 0 ) continue; + if ( p_css+dss < 0 ) continue; - idx = fs+cfs+image->width*(ss+css); + idx = dfs+cfs+image->width*(dss+css); /* Veto this peak if we tried to integrate in a bad region */ if ( image->flags != NULL ) { @@ -332,8 +346,8 @@ static int integrate_peak(struct image *image, int cfs, int css, pk_counts++; pk_total += val; - fsct += val*(cfs+fs); - ssct += val*(css+ss); + fsct += val*(cfs+dfs); + ssct += val*(css+dss); } } @@ -726,7 +740,7 @@ void integrate_reflections(struct image *image, int use_closer, int bgsub, 1); /* I/sigma(I) cutoff */ - if ( !r && (intensity/sigma < min_snr) ) r = 1; + if ( intensity/sigma < min_snr ) r = 1; /* Record intensity and set redundancy to 1 on success */ if ( r == 0 ) { -- cgit v1.2.3 From 172293570b790ffcde14128404d9c0e294ca1f30 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Thu, 28 Jun 2012 17:23:46 +0200 Subject: Remove peak_sep from detector geometry file, use ir_inn instead --- libcrystfel/src/peaks.c | 26 +++++++++++--------------- 1 file changed, 11 insertions(+), 15 deletions(-) (limited to 'libcrystfel/src/peaks.c') diff --git a/libcrystfel/src/peaks.c b/libcrystfel/src/peaks.c index fe836c3f..121ecb6c 100644 --- a/libcrystfel/src/peaks.c +++ b/libcrystfel/src/peaks.c @@ -390,7 +390,6 @@ static void search_peaks_in_panel(struct image *image, float threshold, int nrej_snr = 0; int nacc = 0; int ncull; - const int pws = p->peak_sep/2; data = image->data; stride = image->width; @@ -433,16 +432,14 @@ static void search_peaks_in_panel(struct image *image, float threshold, max = data[mask_fs+stride*mask_ss]; did_something = 0; - for ( s_ss=biggest(mask_ss-pws/2, - p->min_ss); - s_ss<=smallest(mask_ss+pws/2, - p->max_ss); - s_ss++ ) { - for ( s_fs=biggest(mask_fs-pws/2, - p->min_fs); - s_fs<=smallest(mask_fs+pws/2, - p->max_fs); - s_fs++ ) { + for ( s_ss=biggest(mask_ss-ir_inn, p->min_ss); + s_ss<=smallest(mask_ss+ir_inn, p->max_ss); + s_ss++ ) + { + for ( s_fs=biggest(mask_fs-ir_inn, p->min_fs); + s_fs<=smallest(mask_fs+ir_inn, p->max_fs); + s_fs++ ) + { if ( data[s_fs+stride*s_ss] > max ) { max = data[s_fs+stride*s_ss]; @@ -455,8 +452,7 @@ static void search_peaks_in_panel(struct image *image, float threshold, } /* Abort if drifted too far from the foot point */ - if ( distance(mask_fs, mask_ss, fs, ss) > - p->peak_sep/2.0 ) + if ( distance(mask_fs, mask_ss, fs, ss) > ir_inn ) { break; } @@ -464,7 +460,7 @@ static void search_peaks_in_panel(struct image *image, float threshold, } while ( did_something ); /* Too far from foot point? */ - if ( distance(mask_fs, mask_ss, fs, ss) > p->peak_sep/2.0 ) { + if ( distance(mask_fs, mask_ss, fs, ss) > ir_inn ) { nrej_dis++; continue; } @@ -500,7 +496,7 @@ static void search_peaks_in_panel(struct image *image, float threshold, /* Check for a nearby feature */ image_feature_closest(image->features, f_fs, f_ss, &d, &idx); - if ( d < p->peak_sep/2.0 ) { + if ( d < 2.0*ir_inn ) { nrej_pro++; continue; } -- cgit v1.2.3 From 4a81d9f2d079651d8c7c0f555eba7903d54d7ca7 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Thu, 28 Jun 2012 17:24:15 +0200 Subject: Don't lie about the number of badrow culled peaks --- libcrystfel/src/peaks.c | 1 - 1 file changed, 1 deletion(-) (limited to 'libcrystfel/src/peaks.c') diff --git a/libcrystfel/src/peaks.c b/libcrystfel/src/peaks.c index 121ecb6c..cb98f4e9 100644 --- a/libcrystfel/src/peaks.c +++ b/libcrystfel/src/peaks.c @@ -99,7 +99,6 @@ static int cull_peaks_in_panel(struct image *image, struct panel *p) if ( ncol <= 3 ) continue; /* Yes? Delete them all... */ - nelim = 0; for ( j=0; jfeatures, j); -- cgit v1.2.3 From 18e3a7f47dfb9fd5ca732d4f21404e29e88873da Mon Sep 17 00:00:00 2001 From: Thomas White Date: Thu, 28 Jun 2012 17:24:44 +0200 Subject: Warn user if any peaks were badrow culled --- libcrystfel/src/peaks.c | 6 ++++++ 1 file changed, 6 insertions(+) (limited to 'libcrystfel/src/peaks.c') diff --git a/libcrystfel/src/peaks.c b/libcrystfel/src/peaks.c index cb98f4e9..ffe77290 100644 --- a/libcrystfel/src/peaks.c +++ b/libcrystfel/src/peaks.c @@ -521,6 +521,12 @@ static void search_peaks_in_panel(struct image *image, float threshold, // "%i in bad regions, %i with SNR < %g, %i badrow culled.\n", // nacc, nrej_dis, nrej_pro, nrej_fra, nrej_bad, // nrej_snr, min_snr, ncull); + + if ( ncull != 0 ) { + STATUS("WARNING: %i peaks were badrow culled. This feature" + " should not usually be used.\nConsider setting" + " badrow=- in the geometry file.\n", ncull); + } } -- cgit v1.2.3