From f094e70fd60899579a8452c1c08121f25e0be885 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Sat, 5 Dec 2015 13:35:56 -0800 Subject: Add option for per-pixel saturation values --- libcrystfel/src/detector.c | 7 +++ libcrystfel/src/detector.h | 2 + libcrystfel/src/hdf5-file.c | 116 ++++++++++++++++++++++++++++++++++++++++-- libcrystfel/src/image.h | 1 + libcrystfel/src/integration.c | 7 ++- 5 files changed, 127 insertions(+), 6 deletions(-) diff --git a/libcrystfel/src/detector.c b/libcrystfel/src/detector.c index cb0bff95..7689f5f8 100644 --- a/libcrystfel/src/detector.c +++ b/libcrystfel/src/detector.c @@ -956,6 +956,11 @@ static int parse_field_for_panel(struct panel *panel, const char *key, } else if ( strcmp(key, "mask_file") == 0 ) { panel->mask_file = strdup(val); + } else if ( strcmp(key, "saturation_map") == 0 ) { + panel->satmap = strdup(val); + } else if ( strcmp(key, "saturation_map_file") == 0 ) { + panel->satmap_file = strdup(val); + } else if ( strcmp(key, "coffset") == 0) { panel->coffset = atof(val); } else if ( strcmp(key, "res") == 0 ) { @@ -1263,6 +1268,8 @@ struct detector *get_detector_geometry(const char *filename, det->defaults.max_adu = +INFINITY; det->defaults.mask = NULL; det->defaults.mask_file = NULL; + det->defaults.satmap = NULL; + det->defaults.satmap_file = NULL; det->defaults.data = NULL; det->defaults.dim_structure = NULL; strncpy(det->defaults.name, "", 1023); diff --git a/libcrystfel/src/detector.h b/libcrystfel/src/detector.h index 5e544b75..540db4b6 100644 --- a/libcrystfel/src/detector.h +++ b/libcrystfel/src/detector.h @@ -97,6 +97,8 @@ struct panel char *clen_from; char *mask; char *mask_file; + char *satmap; + char *satmap_file; 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 */ diff --git a/libcrystfel/src/hdf5-file.c b/libcrystfel/src/hdf5-file.c index b159327e..09a34c31 100644 --- a/libcrystfel/src/hdf5-file.c +++ b/libcrystfel/src/hdf5-file.c @@ -1126,13 +1126,16 @@ static void debodge_saturation(struct hdfile *f, struct image *image) static int unpack_panels(struct image *image, struct detector *det, - float *data, uint16_t *flags) + float *data, uint16_t *flags, float *sat) { int pi; image->dp = malloc(det->n_panels * sizeof(float *)); image->bad = malloc(det->n_panels * sizeof(int *)); - if ( (image->dp == NULL) || (image->bad == NULL) ) { + image->sat = malloc(det->n_panels * sizeof(float *)); + if ( (image->dp == NULL) || (image->bad == NULL) + || (image->sat == NULL) ) + { ERROR("Failed to allocate panels.\n"); return 1; } @@ -1145,7 +1148,10 @@ static int unpack_panels(struct image *image, struct detector *det, p = &det->panels[pi]; image->dp[pi] = malloc(p->w*p->h*sizeof(float)); image->bad[pi] = calloc(p->w*p->h, sizeof(int)); - if ( (image->dp[pi] == NULL) || (image->bad[pi] == NULL) ) { + image->sat[pi] = malloc(p->w*p->h*sizeof(float)); + if ( (image->dp[pi] == NULL) || (image->bad[pi] == NULL) + || (image->sat[pi] == NULL) ) + { ERROR("Failed to allocate panel\n"); return 1; } @@ -1163,6 +1169,10 @@ static int unpack_panels(struct image *image, struct detector *det, image->dp[pi][fs+p->w*ss] = data[idx]; + if ( sat != NULL ) { + image->sat[pi][fs+p->w*ss] = sat[idx]; + } + if ( p->no_index ) bad = 1; if ( in_bad_region(det, cfs, css) ) { @@ -1452,7 +1462,7 @@ int hdf5_read(struct hdfile *f, struct image *image, const char *element, } image->det = simple_geometry(image); - unpack_panels(image, image->det, buf, NULL); + unpack_panels(image, image->det, buf, NULL, NULL); if ( satcorr ) debodge_saturation(f, image); if ( image->beam != NULL ) { @@ -1551,15 +1561,95 @@ err: } +static int load_satmap(struct hdfile *f, struct event *ev, char *satmap, + const char *satmap_file, + const char *pname, struct image *image, + size_t p_w, size_t sum_p_h, float *smap, + hsize_t *f_offset, hsize_t *f_count, + hsize_t *m_offset, hsize_t *m_count) +{ + hid_t satmap_dataspace, satmap_dh; + int exists; + int check, r; + hid_t memspace; + hsize_t dimsm[2]; + hid_t fh; + + if ( satmap_file != NULL ) { + fh = H5Fopen(satmap_file, H5F_ACC_RDONLY, H5P_DEFAULT); + if ( fh < 0 ) { + ERROR("Couldn't open satmap file '%s'\n", satmap_file); + return 1; + } + } else { + fh = f->fh; + } + + if ( ev != NULL ) { + satmap = retrieve_full_path(ev, satmap); + } + + exists = check_path_existence(fh, satmap); + if ( !exists ) { + ERROR("Cannot find satmap for panel %s\n", pname); + goto err; + } + + satmap_dh = H5Dopen2(fh, satmap, H5P_DEFAULT); + if ( satmap_dh <= 0 ) { + ERROR("Couldn't open satmap for panel %s\n", pname); + goto err; + } + + satmap_dataspace = H5Dget_space(satmap_dh); + check = H5Sselect_hyperslab(satmap_dataspace, H5S_SELECT_SET, + f_offset, NULL, f_count, NULL); + if ( check < 0 ) { + ERROR("Error selecting satmap dataspace for panel %s\n", pname); + goto err; + } + + dimsm[0] = sum_p_h; + dimsm[1] = p_w; + memspace = H5Screate_simple(2, dimsm, NULL); + check = H5Sselect_hyperslab(memspace, H5S_SELECT_SET, + m_offset, NULL, m_count, NULL); + if ( check < 0 ) { + ERROR("Error selecting memory dataspace for panel %s\n", pname); + goto err; + } + + r = H5Dread(satmap_dh, H5T_NATIVE_FLOAT, memspace, + satmap_dataspace, H5P_DEFAULT, smap); + if ( r < 0 ) { + ERROR("Couldn't read satmap for panel %s\n", pname); + goto err; + } + + H5Sclose(satmap_dataspace); + H5Dclose(satmap_dh); + if ( ev != NULL ) free(satmap); + + return 0; + +err: + if ( satmap_file != NULL ) H5Fclose(fh); + if ( ev != NULL ) free(satmap); + return 1; +} + + int hdf5_read2(struct hdfile *f, struct image *image, struct event *ev, int satcorr) { herr_t r; float *buf; uint16_t *flags; + float *smap; int sum_p_h; int p_w; int pi; + int i; if ( image->det == NULL ) { ERROR("Geometry not available\n"); @@ -1594,6 +1684,12 @@ int hdf5_read2(struct hdfile *f, struct image *image, struct event *ev, return 1; } + smap = calloc(p_w*sum_p_h,sizeof(float)); + if ( smap == NULL ) { + ERROR("Failed to allocate memory for satmap\n"); + return 1; + } + for ( i=0; idet->n_panels; pi++ ) { @@ -1742,6 +1838,15 @@ int hdf5_read2(struct hdfile *f, struct image *image, struct event *ev, } } + if ( p->satmap != NULL ) { + if ( load_satmap(f, ev, p->satmap, p->satmap_file, + p->name, image, p_w, sum_p_h, smap, + f_offset, f_count, m_offset, m_count) ) + { + ERROR("Error loading saturation map!\n"); + } + } + free(f_offset); free(f_count); @@ -1749,7 +1854,7 @@ int hdf5_read2(struct hdfile *f, struct image *image, struct event *ev, fill_in_values(image->det, f, ev); - unpack_panels(image, image->det, buf, flags); + unpack_panels(image, image->det, buf, flags, smap); if ( satcorr ) debodge_saturation(f, image); if ( image->beam != NULL ) { @@ -1768,6 +1873,7 @@ int hdf5_read2(struct hdfile *f, struct image *image, struct event *ev, free(buf); free(flags); + free(smap); return 0; } diff --git a/libcrystfel/src/image.h b/libcrystfel/src/image.h index 52614d06..37904308 100644 --- a/libcrystfel/src/image.h +++ b/libcrystfel/src/image.h @@ -161,6 +161,7 @@ struct image { float **dp; /* Data in panel */ int **bad; /* Bad pixels by panel */ + float **sat; /* Per-pixel saturation values */ Crystal **crystals; int n_crystals; diff --git a/libcrystfel/src/integration.c b/libcrystfel/src/integration.c index 730a14ae..b3bee120 100644 --- a/libcrystfel/src/integration.c +++ b/libcrystfel/src/integration.c @@ -804,6 +804,7 @@ static int check_box(struct intcontext *ic, struct peak_box *bx, int *sat) double hd, kd, ld; signed int h, k, l; struct rvec dv; + float lsat; fs = bx->cfs + p; ss = bx->css + q; @@ -844,9 +845,13 @@ static int check_box(struct intcontext *ic, struct peak_box *bx, int *sat) } } + /* Per-pixel saturation value */ + lsat = ic->image->sat[bx->pn][fs + bx->p->w*ss]; if ( (bx->bm[p+ic->w*q] != BM_IG) && (bx->bm[p+ic->w*q] != BM_BH) - && (boxi(ic, bx, p, q) > bx->p->max_adu) ) { + && ((boxi(ic, bx, p, q) > bx->p->max_adu) + || (boxi(ic, bx, p, q) > lsat)) ) + { if ( sat != NULL ) *sat = 1; } -- cgit v1.2.3