diff options
Diffstat (limited to 'src/peaks.c')
| -rw-r--r-- | src/peaks.c | 58 |
1 files changed, 44 insertions, 14 deletions
diff --git a/src/peaks.c b/src/peaks.c index 69859c9a..e9a10ab7 100644 --- a/src/peaks.c +++ b/src/peaks.c @@ -32,6 +32,7 @@ #include "diffraction.h" +/* Maximum number of peaks which may be predicted by find_projected_peaks() */ #define MAX_HITS (1024) /* How close a peak must be to an indexed position to be considered "close" @@ -45,11 +46,16 @@ /* Degree of polarisation of X-ray beam */ #define POL (1.0) - +/* Window size for Zaefferer peak detection */ #define PEAK_WINDOW_SIZE (10) -#define MAX_PEAKS (2048) + +/* Integration radius for peaks */ #define INTEGRATION_RADIUS (10) +/* Integration radius for background estimation */ +#define OUT_INTEGRATION_RADIUS (15) + + static int in_streak(int x, int y) { if ( (y>512) && (y<600) && (abs(x-489)<15) ) return 1; @@ -179,16 +185,21 @@ 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 xp, int yp, float *xc, float *yc, float *intensity, + double *pbg, double *pmax, int do_polar, int do_sa) { signed int x, y; const int lim = INTEGRATION_RADIUS * INTEGRATION_RADIUS; - double total = 0; + const int out_lim = OUT_INTEGRATION_RADIUS * OUT_INTEGRATION_RADIUS; + double total = 0.0; int xct = 0; int yct = 0; + double noise = 0.0; + int noise_counts = 0; + double max = 0.0; - for ( x=-INTEGRATION_RADIUS; x<+INTEGRATION_RADIUS; x++ ) { - for ( y=-INTEGRATION_RADIUS; y<+INTEGRATION_RADIUS; y++ ) { + for ( x=-OUT_INTEGRATION_RADIUS; x<+OUT_INTEGRATION_RADIUS; x++ ) { + for ( y=-OUT_INTEGRATION_RADIUS; y<+OUT_INTEGRATION_RADIUS; y++ ) { struct panel *p = NULL; double val, sa, pix_area, Lsq, dsq, proj_area; @@ -196,8 +207,8 @@ int integrate_peak(struct image *image, int xp, int yp, double phi, pa, pb, pol; uint16_t flags; - /* Circular mask */ - if ( x*x + y*y > lim ) continue; + /* Outer mask radius */ + if ( x*x + y*y > out_lim ) continue; if ( ((x+xp)>=image->width) || ((x+xp)<0) ) continue; if ( ((y+yp)>=image->height) || ((y+yp)<0) ) continue; @@ -210,6 +221,16 @@ int integrate_peak(struct image *image, int xp, int yp, val = image->data[(x+xp)+image->width*(y+yp)]; + /* Inner mask */ + if ( x*x + y*y > lim ) { + /* Estimate noise from this region */ + noise += fabs(val); + noise_counts++; + continue; + } + + if ( val > max ) max = val; + p = find_panel(image->det, x+xp, y+yp); if ( p == NULL ) return 1; @@ -268,6 +289,13 @@ int integrate_peak(struct image *image, int xp, int yp, *intensity = 0; } + if ( pbg != NULL ) { + *pbg = (noise / noise_counts); + } + if ( pmax != NULL ) { + *pmax = max; + } + return 0; } @@ -383,7 +411,7 @@ void search_peaks(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_x, mask_y, - &fx, &fy, &intensity, 0, 0); + &fx, &fy, &intensity, NULL, NULL, 0, 0); if ( r ) { /* Bad region - don't detect peak */ nrej_bad++; @@ -701,6 +729,7 @@ void output_intensities(struct image *image, UnitCell *cell, float x, y, intensity; double d; int idx; + double bg, max; /* Wait.. is there a really close feature which was detected? */ if ( use_closer ) { @@ -724,7 +753,8 @@ void output_intensities(struct image *image, UnitCell *cell, * coordinates. This will produce further * revised coordinates. */ r = integrate_peak(image, f->x, f->y, &x, &y, - &intensity, polar, sa); + &intensity, &bg, &max, + polar, sa); if ( r ) { /* The original peak (which also went * through integrate_peak(), but with @@ -745,8 +775,8 @@ void output_intensities(struct image *image, UnitCell *cell, r = integrate_peak(image, image->hits[i].x, image->hits[i].y, - &x, &y, &intensity, polar, - sa); + &x, &y, &intensity, &bg, &max, + polar, sa); if ( r ) { /* Plain old ordinary peak veto */ n_veto++; @@ -766,7 +796,7 @@ void output_intensities(struct image *image, UnitCell *cell, r = integrate_peak(image, image->hits[i].x, image->hits[i].y, - &x, &y, &intensity, polar, + &x, &y, &intensity, &bg, &max, polar, sa); if ( r ) { /* Plain old ordinary peak veto */ @@ -777,9 +807,9 @@ void output_intensities(struct image *image, UnitCell *cell, } /* Write h,k,l, integrated intensity and centroid coordinates */ - fprintf(ofh, "%3i %3i %3i %6f (at %5.2f,%5.2f)\n", + fprintf(ofh, "%3i %3i %3i %6f (at %5.2f,%5.2f) max=%6f bg=%6f\n", image->hits[i].h, image->hits[i].k, image->hits[i].l, - intensity, x, y); + intensity, x, y, max, bg); image->hits[i].x = x; image->hits[i].y = y; |