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 +++++++----- tests/integration_check.c | 10 +++++----- 2 files changed, 12 insertions(+), 10 deletions(-) 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 ) { diff --git a/tests/integration_check.c b/tests/integration_check.c index 2d18ac5a..21c10f79 100644 --- a/tests/integration_check.c +++ b/tests/integration_check.c @@ -64,7 +64,7 @@ static void third_integration_check(struct image *image, int n_trials, } r = integrate_peak(image, 64, 64, &fsp, &ssp, - &intensity, &sigma, 10.0, 15.0, 17.0); + &intensity, &sigma, 10.0, 15.0, 17.0, 0); if ( r == 0 ) { mean_intensity += intensity; @@ -125,7 +125,7 @@ static void fourth_integration_check(struct image *image, int n_trials, } r = integrate_peak(image, 64, 64, &fsp, &ssp, - &intensity, &sigma, 10.0, 15.0, 17.0); + &intensity, &sigma, 10.0, 15.0, 17.0, 0); if ( r == 0 ) { mean_intensity += intensity; @@ -205,7 +205,7 @@ int main(int argc, char *argv[]) /* 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); + 10.0, 15.0, 17.0, 0); STATUS(" First check: integrate_peak() returned %i", r); if ( r == 0 ) { @@ -231,7 +231,7 @@ int main(int argc, char *argv[]) } r = integrate_peak(&image, 64, 64, &fsp, &ssp, &intensity, &sigma, - 10.0, 15.0, 17.0); + 10.0, 15.0, 17.0, 0); if ( r ) { ERROR(" Second check: integrate_peak() returned %i (wrong).\n", r); @@ -273,7 +273,7 @@ int main(int argc, char *argv[]) } r = integrate_peak(&image, 64, 64, &fsp, &ssp, &intensity, &sigma, - 10.0, 15.0, 17.0); + 10.0, 15.0, 17.0, 0); if ( r ) { ERROR(" Fifth check: integrate_peak() returned %i (wrong).\n", r); -- 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 +++++++++++++++++++++++++++++++++++++++++++++++ tests/integration_check.c | 2 ++ 2 files changed, 59 insertions(+) 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; } diff --git a/tests/integration_check.c b/tests/integration_check.c index 2d18ac5a..7fe98fa0 100644 --- a/tests/integration_check.c +++ b/tests/integration_check.c @@ -203,6 +203,8 @@ int main(int argc, char *argv[]) image.height = 128; memset(image.data, 0, 128*128*sizeof(float)); + image.reflections=reflist_new(); + /* 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); -- 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(-) 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(+) 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 --- AUTHORS | 3 +++ README | 1 + libcrystfel/src/peaks.c | 1 + 3 files changed, 5 insertions(+) diff --git a/AUTHORS b/AUTHORS index 54a29e7b..5ce6a9c1 100644 --- a/AUTHORS +++ b/AUTHORS @@ -1,6 +1,9 @@ * Thomas White Lead author, architecture, "vision" and so on. +* Kenneth Beyerlein + Peak integration improvements + * Richard Kirian MOSFLM auto-indexing, bug fixing. diff --git a/README b/README index 700e646c..23237849 100644 --- a/README +++ b/README @@ -7,6 +7,7 @@ Copyright © 2012 Deutsches Elektronen-Synchrotron DESY, Authors: Thomas White Richard Kirian + Kenneth Beyerlein Andrew Aquila Andrew Martin Lorenzo Galli 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 +++++++++++++++++------------------- libcrystfel/src/reflist.c | 5 ++++- src/indexamajig.c | 1 + 3 files changed, 22 insertions(+), 20 deletions(-) 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; diff --git a/libcrystfel/src/reflist.c b/libcrystfel/src/reflist.c index cea25bf3..64154774 100644 --- a/libcrystfel/src/reflist.c +++ b/libcrystfel/src/reflist.c @@ -900,6 +900,7 @@ struct _reflistiterator { **/ Reflection *first_refl(RefList *list, RefListIterator **piter) { + Reflection *refl; RefListIterator *iter; iter = malloc(sizeof(struct _reflistiterator)); @@ -908,7 +909,9 @@ Reflection *first_refl(RefList *list, RefListIterator **piter) iter->stack_ptr = 0; *piter = iter; - Reflection *refl = list->head; + if ( list == NULL ) return NULL; + + refl = list->head; do { diff --git a/src/indexamajig.c b/src/indexamajig.c index 8c073e2c..c3e1bedf 100644 --- a/src/indexamajig.c +++ b/src/indexamajig.c @@ -264,6 +264,7 @@ static void process_image(void *pp, int cookie) image.data = NULL; image.flags = NULL; image.indexed_cell = NULL; + image.reflections = NULL; image.id = cookie; image.filename = filename; image.det = copy_geom(pargs->static_args.det); -- cgit v1.2.3 From 8e6a4919ca0aee13ee98ffefc9054ec5f8661740 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Thu, 21 Jun 2012 12:20:24 +0200 Subject: Fix skip_some_files() Without chomp(), the string comparison never works. --- libcrystfel/src/stream.c | 1 + 1 file changed, 1 insertion(+) diff --git a/libcrystfel/src/stream.c b/libcrystfel/src/stream.c index 81f4d722..ce0ee33e 100644 --- a/libcrystfel/src/stream.c +++ b/libcrystfel/src/stream.c @@ -474,6 +474,7 @@ int skip_some_files(FILE *fh, int n) rval = fgets(line, 1023, fh); if ( rval == NULL ) continue; + chomp(line); if ( strcmp(line, CHUNK_END_MARKER) == 0 ) n_patterns++; } while ( rval != NULL ); -- 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(-) 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 714ceebf6b568a94b1e2ddbfd562307697abbfd7 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Sat, 23 Jun 2012 21:47:37 +0200 Subject: Fix missing word in crystfel_geometry(5) --- doc/man/crystfel_geometry.5 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/man/crystfel_geometry.5 b/doc/man/crystfel_geometry.5 index 006af722..4fa22da1 100644 --- a/doc/man/crystfel_geometry.5 +++ b/doc/man/crystfel_geometry.5 @@ -47,7 +47,7 @@ of pixel values in the HDF5 file, defined in terms only of the "fast scan" and +x completes the right-handed coordinate system. .PP -Naively speaking, this means that CrystFEL at the images from the "into the +Naively speaking, this means that CrystFEL looks at the images from the "into the beam" perspective, but please avoid thinking of things in this way. It's much better to consider the precise way in which the coordinates are mapped. -- cgit v1.2.3 From b9d70f03c58ad8b2479b0e30a296b41d75a11a61 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Tue, 26 Jun 2012 15:14:28 +0200 Subject: Make units of corner_* explicity in crystfel_geometry(5) --- doc/man/crystfel_geometry.5 | 2 ++ 1 file changed, 2 insertions(+) diff --git a/doc/man/crystfel_geometry.5 b/doc/man/crystfel_geometry.5 index 4fa22da1..68abb341 100644 --- a/doc/man/crystfel_geometry.5 +++ b/doc/man/crystfel_geometry.5 @@ -116,6 +116,8 @@ panel0/ss = -x .br ; the HDF5 file, is now given a position in the lab coordinate system. .br +; Units are pixels. +.br ; Note that "first point in the panel" is a conceptual simplification. We refer .br ; to that corner, and to the very corner of the pixel - NOT, for example, to the -- cgit v1.2.3 From 51f5ed0461208f813a7c4ab45546516223dd4d5b Mon Sep 17 00:00:00 2001 From: Thomas White Date: Thu, 28 Jun 2012 17:23:14 +0200 Subject: Restore bandwidth and subsampling to pattern_sim --- data/diffraction.cl | 15 ++-- src/diffraction-gpu.c | 216 ++++++++++++++++++++++++++++++-------------------- src/pattern_sim.c | 2 + 3 files changed, 141 insertions(+), 92 deletions(-) diff --git a/data/diffraction.cl b/data/diffraction.cl index f6318378..d128446e 100644 --- a/data/diffraction.cl +++ b/data/diffraction.cl @@ -33,7 +33,7 @@ const sampler_t sampler_c = CLK_NORMALIZED_COORDS_TRUE float4 get_q(float fs, float ss, float res, float clen, float k, - float *ttp, float corner_x, float corner_y, + float corner_x, float corner_y, float fsx, float fsy, float ssx, float ssy) { float rx, ry, r; @@ -51,7 +51,6 @@ float4 get_q(float fs, float ss, float res, float clen, float k, r = sqrt(pow(rx, 2.0f) + pow(ry, 2.0f)); tt = atan2(r, clen); - *ttp = tt; az = atan2(ry, rx); @@ -211,7 +210,7 @@ float molecule_factor(global float *intensities, global float *flags, } -kernel void diffraction(global float *diff, global float *ttp, float k, +kernel void diffraction(global float *diff, float k, int w, float corner_x, float corner_y, float res, float clen, float16 cell, global float *intensities, @@ -219,7 +218,8 @@ kernel void diffraction(global float *diff, global float *ttp, float k, read_only image2d_t func_b, read_only image2d_t func_c, global float *flags, - float fsx, float fsy, float ssx, float ssy) + float fsx, float fsy, float ssx, float ssy, + float xo, float yo) { float tt; float fs, ss; @@ -230,11 +230,11 @@ kernel void diffraction(global float *diff, global float *ttp, float k, int idx; /* Calculate fractional coordinates in fs/ss */ - fs = convert_float(get_global_id(0)); - ss = convert_float(get_global_id(1)); + fs = convert_float(get_global_id(0)) + xo; + ss = convert_float(get_global_id(1)) + yo; /* Get the scattering vector */ - q = get_q(fs, ss, res, clen, k, &tt, + q = get_q(fs, ss, res, clen, k, corner_x, corner_y, fsx, fsy, ssx, ssy); /* Calculate the diffraction */ @@ -245,5 +245,4 @@ kernel void diffraction(global float *diff, global float *ttp, float k, /* Write the value to memory */ idx = convert_int_rtz(fs) + w*convert_int_rtz(ss); diff[idx] = I_molecule * I_lattice; - ttp[idx] = tt; } diff --git a/src/diffraction-gpu.c b/src/diffraction-gpu.c index 2b8e97e5..15ff39cc 100644 --- a/src/diffraction-gpu.c +++ b/src/diffraction-gpu.c @@ -167,57 +167,13 @@ static int set_arg_mem(struct gpu_context *gctx, int idx, cl_mem val) } -void get_diffraction_gpu(struct gpu_context *gctx, struct image *image, - int na, int nb, int nc, UnitCell *ucell) +static void do_panels(struct gpu_context *gctx, struct image *image, + int *n_inf, int *n_neg, int *n_nan, double k, double frac, + double xo, double yo) { - cl_int err; - double ax, ay, az; - double bx, by, bz; - double cx, cy, cz; int i; - cl_float16 cell; - cl_int4 ncells; - int n_inf = 0; - int n_neg = 0; - int n_nan = 0; - if ( gctx == NULL ) { - ERROR("GPU setup failed.\n"); - return; - } - - cell_get_cartesian(ucell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz); - cell.s[0] = ax; cell.s[1] = ay; cell.s[2] = az; - cell.s[3] = bx; cell.s[4] = by; cell.s[5] = bz; - cell.s[6] = cx; cell.s[7] = cy; cell.s[8] = cz; - - ncells.s[0] = na; - ncells.s[1] = nb; - ncells.s[2] = nc; - ncells.s[3] = 0; /* unused */ - - /* Ensure all required LUTs are available */ - check_sinc_lut(gctx, na); - check_sinc_lut(gctx, nb); - check_sinc_lut(gctx, nc); - - if ( set_arg_float(gctx, 2, 1.0/image->lambda) ) return; - if ( set_arg_mem(gctx, 9, gctx->intensities) ) return; - if ( set_arg_mem(gctx, 10, gctx->sinc_luts[na-1]) ) return; - if ( set_arg_mem(gctx, 11, gctx->sinc_luts[nb-1]) ) return; - if ( set_arg_mem(gctx, 12, gctx->sinc_luts[nc-1]) ) return; - if ( set_arg_mem(gctx, 13, gctx->flags) ) return; - - /* Unit cell */ - err = clSetKernelArg(gctx->kern, 8, sizeof(cl_float16), &cell); - if ( err != CL_SUCCESS ) { - ERROR("Couldn't set unit cell: %s\n", clError(err)); - return; - } - - /* Allocate memory for the result */ - image->data = calloc(image->width * image->height, sizeof(float)); - image->twotheta = calloc(image->width * image->height, sizeof(double)); + if ( set_arg_float(gctx, 1, k) ) return; /* Iterate over panels */ for ( i=0; idet->n_panels; i++ ) { @@ -225,14 +181,12 @@ void get_diffraction_gpu(struct gpu_context *gctx, struct image *image, size_t dims[2]; size_t ldims[2] = {1, 1}; struct panel *p; - cl_mem tt; - size_t tt_size; cl_mem diff; size_t diff_size; float *diff_ptr; - float *tt_ptr; int pan_width, pan_height; int fs, ss; + cl_int err; p = &image->det->panels[i]; @@ -247,25 +201,19 @@ void get_diffraction_gpu(struct gpu_context *gctx, struct image *image, ERROR("Couldn't allocate diffraction memory\n"); return; } - tt_size = pan_width * pan_height * sizeof(cl_float); - tt = clCreateBuffer(gctx->ctx, CL_MEM_WRITE_ONLY, tt_size, - NULL, &err); - if ( err != CL_SUCCESS ) { - ERROR("Couldn't allocate twotheta memory\n"); - return; - } if ( set_arg_mem(gctx, 0, diff) ) return; - if ( set_arg_mem(gctx, 1, tt) ) return; - if ( set_arg_int(gctx, 3, pan_width) ) return; - if ( set_arg_float(gctx, 4, p->cnx) ) return; - if ( set_arg_float(gctx, 5, p->cny) ) return; - if ( set_arg_float(gctx, 6, p->res) ) return; - if ( set_arg_float(gctx, 7, p->clen) ) return; - if ( set_arg_float(gctx, 14, p->fsx) ) return; - if ( set_arg_float(gctx, 15, p->fsy) ) return; - if ( set_arg_float(gctx, 16, p->ssx) ) return; - if ( set_arg_float(gctx, 17, p->ssy) ) return; + if ( set_arg_int(gctx, 2, pan_width) ) return; + if ( set_arg_float(gctx, 3, p->cnx) ) return; + if ( set_arg_float(gctx, 4, p->cny) ) return; + if ( set_arg_float(gctx, 5, p->res) ) return; + if ( set_arg_float(gctx, 6, p->clen) ) return; + if ( set_arg_float(gctx, 13, p->fsx) ) return; + if ( set_arg_float(gctx, 14, p->fsy) ) return; + if ( set_arg_float(gctx, 15, p->ssx) ) return; + if ( set_arg_float(gctx, 16, p->ssy) ) return; + if ( set_arg_float(gctx, 17, xo) ) return; + if ( set_arg_float(gctx, 18, yo) ) return; dims[0] = pan_width; dims[1] = pan_height; @@ -288,43 +236,124 @@ void get_diffraction_gpu(struct gpu_context *gctx, struct image *image, clError(err)); return; } - tt_ptr = clEnqueueMapBuffer(gctx->cq, tt, CL_TRUE, CL_MAP_READ, - 0, tt_size, 0, NULL, NULL, &err); - if ( err != CL_SUCCESS ) { - ERROR("Couldn't map tt buffer\n"); - return; - } for ( fs=0; fsmin_fs + fs; tss = p->min_ss + ss; - image->data[tfs + image->width*tss] = val; - image->twotheta[tfs + image->width*tss] = tt; + image->data[tfs + image->width*tss] += frac*val; } } clEnqueueUnmapMemObject(gctx->cq, diff, diff_ptr, 0, NULL, NULL); - clEnqueueUnmapMemObject(gctx->cq, tt, tt_ptr, - 0, NULL, NULL); clReleaseMemObject(diff); - clReleaseMemObject(tt); } +} + + +void get_diffraction_gpu(struct gpu_context *gctx, struct image *image, + int na, int nb, int nc, UnitCell *ucell) +{ + double ax, ay, az; + double bx, by, bz; + double cx, cy, cz; + cl_float16 cell; + cl_int4 ncells; + cl_int err; + int n_inf = 0; + int n_neg = 0; + int n_nan = 0; + int fs, ss; + const int nkstep = 10; + const int nxs = 4; + const int nys = 4; + int kstep; + double klow, khigh, kinc, w; + double frac; + int xs, ys; + int n, ntot; + + if ( gctx == NULL ) { + ERROR("GPU setup failed.\n"); + return; + } + cell_get_cartesian(ucell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz); + cell.s[0] = ax; cell.s[1] = ay; cell.s[2] = az; + cell.s[3] = bx; cell.s[4] = by; cell.s[5] = bz; + cell.s[6] = cx; cell.s[7] = cy; cell.s[8] = cz; + + ncells.s[0] = na; + ncells.s[1] = nb; + ncells.s[2] = nc; + ncells.s[3] = 0; /* unused */ + + /* Ensure all required LUTs are available */ + check_sinc_lut(gctx, na); + check_sinc_lut(gctx, nb); + check_sinc_lut(gctx, nc); + + if ( set_arg_mem(gctx, 8, gctx->intensities) ) return; + if ( set_arg_mem(gctx, 9, gctx->sinc_luts[na-1]) ) return; + if ( set_arg_mem(gctx, 10, gctx->sinc_luts[nb-1]) ) return; + if ( set_arg_mem(gctx, 11, gctx->sinc_luts[nc-1]) ) return; + if ( set_arg_mem(gctx, 12, gctx->flags) ) return; + + /* Unit cell */ + err = clSetKernelArg(gctx->kern, 7, sizeof(cl_float16), &cell); + if ( err != CL_SUCCESS ) { + ERROR("Couldn't set unit cell: %s\n", clError(err)); + return; + } + + /* Allocate memory for the result */ + image->data = calloc(image->width * image->height, sizeof(float)); + + w = image->lambda * image->bw; + klow = 1.0 / (image->lambda + (w/2.0)); /* Smallest k */ + khigh = 1.0 / (image->lambda - (w/2.0)); /* Largest k */ + kinc = (khigh-klow) / (nkstep+1); + + ntot = nkstep * nxs * nys; + frac = 1.0 / ntot; + + n = 0; + for ( xs=0; xstwotheta = calloc(image->width * image->height, sizeof(double)); + for ( fs=0; fswidth; fs++ ) { + for ( ss=0; ssheight; ss++ ) { + + struct rvec q; + double twotheta, k; + int idx; + + /* Calculate k this time round */ + k = 1.0/image->lambda; + + q = get_q(image, fs, ss, &twotheta, k); + + idx = fs + image->width*ss; + image->twotheta[idx] = twotheta; + + } + } } diff --git a/src/pattern_sim.c b/src/pattern_sim.c index 03b009b4..ea8a317f 100644 --- a/src/pattern_sim.c +++ b/src/pattern_sim.c @@ -497,6 +497,8 @@ int main(int argc, char *argv[]) image.width = image.det->max_fs + 1; image.height = image.det->max_ss + 1; image.lambda = ph_en_to_lambda(eV_to_J(image.beam->photon_energy)); + image.bw = image.beam->bandwidth; + image.div = image.beam->divergence; /* Load unit cell */ input_cell = load_cell_from_pdb(filename); -- 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/detector.c | 7 ------- libcrystfel/src/detector.h | 1 - libcrystfel/src/peaks.c | 26 +++++++++++--------------- 3 files changed, 11 insertions(+), 23 deletions(-) diff --git a/libcrystfel/src/detector.c b/libcrystfel/src/detector.c index b7d809df..f212821d 100644 --- a/libcrystfel/src/detector.c +++ b/libcrystfel/src/detector.c @@ -608,8 +608,6 @@ static int parse_field_for_panel(struct panel *panel, const char *key, panel->res = atof(val); } else if ( strcmp(key, "max_adu") == 0 ) { panel->max_adu = atof(val); - } else if ( strcmp(key, "peak_sep") == 0 ) { - panel->peak_sep = atof(val); } else if ( strcmp(key, "badrow_direction") == 0 ) { panel->badrow = val[0]; /* First character only */ if ( (panel->badrow != 'x') && (panel->badrow != 'y') @@ -687,8 +685,6 @@ static void parse_toplevel(struct detector *det, const char *key, det->mask_good = v; } - } else if ( strcmp(key, "peak_sep") == 0 ) { - det->defaults.peak_sep = atof(val); } else if ( strcmp(key, "coffset") == 0 ) { det->defaults.coffset = atof(val); } else if ( parse_field_for_panel(&det->defaults, key, val, det) ) { @@ -738,7 +734,6 @@ struct detector *get_detector_geometry(const char *filename) det->defaults.res = -1.0; det->defaults.badrow = '-'; det->defaults.no_index = 0; - det->defaults.peak_sep = 50.0; det->defaults.fsx = 1.0; det->defaults.fsy = 0.0; det->defaults.ssx = 0.0; @@ -887,7 +882,6 @@ struct detector *get_detector_geometry(const char *filename) } /* It's OK if the badrow direction is '0' */ /* It's not a problem if "no_index" is still zero */ - /* The default peak_sep is OK (maybe) */ /* The default transformation matrix is at least valid */ if ( det->panels[i].max_fs > max_fs ) { @@ -1250,7 +1244,6 @@ int write_detector_geometry(const char *filename, struct detector *det) fprintf(fh, "%s/max_ss = %d\n", p->name, p->max_ss); fprintf(fh, "%s/badrow_direction = %C\n", p->name, p->badrow); fprintf(fh, "%s/res = %g\n", p->name, p->res); - fprintf(fh, "%s/peak_sep = %g\n", p->name, p->peak_sep); fprintf(fh, "%s/clen = %s\n", p->name, p->clen_from); fprintf(fh, "%s/fs = %+fx %+fy\n", p->name, p->fsx, p->fsy); fprintf(fh, "%s/ss = %+fx %+fy\n", p->name, p->ssx, p->ssy); diff --git a/libcrystfel/src/detector.h b/libcrystfel/src/detector.h index 2b4ac5c6..b740965c 100644 --- a/libcrystfel/src/detector.h +++ b/libcrystfel/src/detector.h @@ -62,7 +62,6 @@ struct panel double res; /* Resolution in pixels per metre */ char badrow; /* 'x' or 'y' */ int no_index; /* Don't index peaks in this panel if non-zero */ - double peak_sep; /* Characteristic peak separation */ char *rigid_group; /* Rigid group, or -1 for none */ double adu_per_eV; /* Number of ADU per eV */ double max_adu; /* Treat pixel as unreliable if higher than this */ 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(-) 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(+) 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 From 9f86f285dd59583866ad69149fa0b2dd324b918a Mon Sep 17 00:00:00 2001 From: Thomas White Date: Thu, 28 Jun 2012 17:24:56 +0200 Subject: Formatting --- src/diffraction.c | 1 - 1 file changed, 1 deletion(-) diff --git a/src/diffraction.c b/src/diffraction.c index 870569b7..b76154c2 100644 --- a/src/diffraction.c +++ b/src/diffraction.c @@ -420,7 +420,6 @@ void get_diffraction(struct image *image, int na, int nb, int nc, image->data[idx] = I_lattice * I_molecule; image->twotheta[idx] = twotheta; - } progress_bar(fs, image->width-1, "Calculating diffraction"); } -- cgit v1.2.3 From d7cacb7871a7d49057453b768e420dfb0ab52077 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Thu, 28 Jun 2012 17:25:02 +0200 Subject: pattern_sim: Use hdf5_write() instead of hdf5_write_image() This means that the wavelength gets put into the HDF5 file. --- src/pattern_sim.c | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/pattern_sim.c b/src/pattern_sim.c index ea8a317f..e05dcdec 100644 --- a/src/pattern_sim.c +++ b/src/pattern_sim.c @@ -619,8 +619,7 @@ int main(int argc, char *argv[]) number++; /* Write the output file */ - hdf5_write(filename, image.data, - image.width, image.height, H5T_NATIVE_FLOAT); + hdf5_write_image(filename, &image); } -- cgit v1.2.3