diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/hdf5-file.c | 23 | ||||
| -rw-r--r-- | src/image.h | 1 | ||||
| -rw-r--r-- | src/peaks.c | 63 |
3 files changed, 74 insertions, 13 deletions
diff --git a/src/hdf5-file.c b/src/hdf5-file.c index 5573267b..d635923d 100644 --- a/src/hdf5-file.c +++ b/src/hdf5-file.c @@ -318,6 +318,11 @@ int hdf5_read(struct hdfile *f, struct image *image) { herr_t r; float *buf; + uint8_t *flags; + hid_t mask_dh; + + image->height = f->nx; + image->width = f->ny; buf = malloc(sizeof(float)*f->nx*f->ny); @@ -327,10 +332,22 @@ int hdf5_read(struct hdfile *f, struct image *image) ERROR("Couldn't read data\n"); return 1; } - image->data = buf; - image->height = f->nx; - image->width = f->ny; + + mask_dh = H5Dopen(f->fh, "/processing/hitfinder/masks", H5P_DEFAULT); + if ( mask_dh <= 0 ) { + ERROR("Couldn't open flags\n"); + } else { + flags = malloc(sizeof(uint8_t)*f->nx*f->ny); + r = H5Dread(mask_dh, H5T_NATIVE_B8, H5S_ALL, H5S_ALL, + H5P_DEFAULT, flags); + if ( r < 0 ) { + ERROR("Couldn't read flags\n"); + image->flags = NULL; + } else { + image->flags = flags; + } + } /* Read wavelength from file */ image->lambda = get_wavelength(f); diff --git a/src/image.h b/src/image.h index db1b491e..e3a520ca 100644 --- a/src/image.h +++ b/src/image.h @@ -81,6 +81,7 @@ struct reflhit { struct image { float *data; + uint8_t *flags; double *twotheta; UnitCell *indexed_cell; UnitCell *candidate_cells[MAX_CELL_CANDIDATES]; diff --git a/src/peaks.c b/src/peaks.c index 22e12949..0d4292d3 100644 --- a/src/peaks.c +++ b/src/peaks.c @@ -133,7 +133,9 @@ static void cull_peaks(struct image *image) } -static void integrate_peak(struct image *image, int xp, int yp, +/* Returns non-zero if peak has been vetoed. + * i.e. don't use result if return value is not zero. */ +static int integrate_peak(struct image *image, int xp, int yp, float *xc, float *yc, float *intensity, int do_polar) { @@ -150,6 +152,7 @@ static void integrate_peak(struct image *image, int xp, int yp, double val, sa, pix_area, Lsq, dsq, proj_area; float tt; double phi, pa, pb, pol; + uint8_t flags; /* Circular mask */ if ( x*x + y*y > lim ) continue; @@ -157,6 +160,12 @@ static void integrate_peak(struct image *image, int xp, int yp, if ( ((x+xp)>=image->width) || ((x+xp)<0) ) continue; if ( ((y+yp)>=image->height) || ((y+yp)<0) ) continue; + /* Veto this peak if we tried to integrate in a bad region */ + if ( image->flags != NULL ) { + flags = image->flags[(x+xp)+image->width*(y+yp)]; + if ( flags & 0x03 ) return 1; + } + p = find_panel(&image->det, x+xp, y+yp); /* Area of one pixel */ @@ -203,6 +212,8 @@ static void integrate_peak(struct image *image, int xp, int yp, *yc = (float)yp; *intensity = 0; } + + return 0; } @@ -219,6 +230,7 @@ void search_peaks(struct image *image) int nrej_hot = 0; int nrej_pro = 0; int nrej_fra = 0; + int nrej_bad = 0; int nacc = 0; data = image->data; @@ -240,6 +252,7 @@ void search_peaks(struct image *image) int sx, sy; double max; unsigned int did_something; + int r; /* Overall threshold */ if ( data[x+width*y] < 800 ) continue; @@ -313,7 +326,13 @@ void search_peaks(struct image *image) /* Centroid peak and get better coordinates. * Don't bother doing polarisation correction, because the * intensity of this peak is only an estimate at this stage. */ - integrate_peak(image, mask_x, mask_y, &fx, &fy, &intensity, 0); + r = integrate_peak(image, mask_x, mask_y, + &fx, &fy, &intensity, 0); + if ( r ) { + /* Bad region - don't detect peak */ + nrej_bad++; + continue; + } /* It is possible for the centroid to fall outside the image */ if ( (fx < 0.0) || (fx > image->width) @@ -337,8 +356,9 @@ void search_peaks(struct image *image) } } - STATUS("%i accepted, %i box, %i hot, %i proximity, %i outside frame\n", - nacc, nrej_dis, nrej_hot, nrej_pro, nrej_fra); + STATUS("%i accepted, %i box, %i hot, %i proximity, %i outside frame," + "%i in bad regions.\n", + nacc, nrej_dis, nrej_hot, nrej_pro, nrej_fra, nrej_bad); cull_peaks(image); } @@ -492,6 +512,8 @@ void output_intensities(struct image *image, UnitCell *cell, int n_found; int n_indclose = 0; int n_foundclose = 0; + int n_veto = 0; + int n_veto_second = 0; double asx, asy, asz; double bsx, bsy, bsz; double csx, csy, csz; @@ -543,17 +565,35 @@ void output_intensities(struct image *image, UnitCell *cell, &d, &idx); if ( (f != NULL) && (d < PEAK_REALLY_CLOSE) ) { + int r; + /* f->intensity was measured on the filtered pattern, * so instead re-integrate using old coordinates. * This will produce further revised coordinates. */ - integrate_peak(image, f->x, f->y, &x, &y, &intensity, - !unpolar); + r = integrate_peak(image, f->x, f->y, &x, &y, + &intensity, !unpolar); + if ( r ) { + /* The original peak (which also went through + * integrate_peak(), but with the mangled + * image data) would have been rejected if it + * was in a bad region. Integration of the same + * peak included a bad region this time. */ + n_veto_second++; + continue; + } intensity = f->intensity; } else { - integrate_peak(image, hits[i].x, hits[i].y, - &x, &y, &intensity, !unpolar); + int r; + + r = integrate_peak(image, hits[i].x, hits[i].y, + &x, &y, &intensity, !unpolar); + if ( r ) { + /* Plain old ordinary peak veto */ + n_veto++; + continue; + } } @@ -592,9 +632,12 @@ void output_intensities(struct image *image, UnitCell *cell, printf("Peak statistics: %i peaks found by the peak search out of " "%i were close to indexed positions. " - "%i indexed positions out of %i were close to detected peaks\n", + "%i indexed positions out of %i were close to detected peaks.\n", n_foundclose, n_found, n_indclose, n_hits); - + printf("%i integrations using indexed locations were aborted because " + "they hit one or more bad pixels.\n", n_veto); + printf("%i integrations using peak search locations were aborted " + "because they hit one or more bad pixels.\n", n_veto_second); /* Blank line at end */ printf("\n"); |