From ef7157ef517cf66dfc5b2c2cfc6602e30a31d060 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 21 Sep 2022 16:22:47 +0200 Subject: Move create_detgeom to DataTemplate module It seems to make more sense here, because it's all about interpreting the contents of the DataTemplate structure. --- libcrystfel/src/datatemplate.c | 227 ++++++++++++++++++++++++++++++++++++ libcrystfel/src/datatemplate_priv.h | 2 + libcrystfel/src/image.c | 225 ----------------------------------- libcrystfel/src/image.h | 3 - libcrystfel/src/stream.c | 1 + 5 files changed, 230 insertions(+), 228 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/datatemplate.c b/libcrystfel/src/datatemplate.c index 1ffa2777..94a8de2a 100644 --- a/libcrystfel/src/datatemplate.c +++ b/libcrystfel/src/datatemplate.c @@ -36,6 +36,7 @@ #include "utils.h" #include "datatemplate.h" +#include "image.h" #include "datatemplate_priv.h" @@ -1734,3 +1735,229 @@ double data_template_get_wavelength_if_possible(const DataTemplate *dt) return NAN; } } + + +static int separate_value_and_units(const char *from, + char **pvalue, + char **punits) +{ + char *sp; + char *fromcpy; + char *unitscpy; + + if ( from == NULL ) return 1; + + fromcpy = strdup(from); + if ( fromcpy == NULL ) return 1; + + sp = strchr(fromcpy, ' '); + if ( sp == NULL ) { + unitscpy = NULL; + } else { + unitscpy = strdup(sp+1); + sp[0] = '\0'; + } + + *pvalue = fromcpy; + *punits = unitscpy; + return 0; +} + + +/* default_scale is a value to be used if both of the following + * conditions are met: + * + * 1. The value is a reference to image headers/metadata, + * rather than a literal number. + * 2. No units are specified in the number. + * + * This is totally horrible. Sorry. Blame history. + */ +static int im_get_length(struct image *image, const char *from, + double default_scale, double *pval) +{ + char *value_str; + char *units; + + if ( from == NULL ) return 1; + + if ( separate_value_and_units(from, &value_str, &units) ) return 1; + + if ( units == NULL ) { + + /* No units given */ + + if ( convert_float(value_str, pval) == 0 ) { + + /* Literal value with no units */ + free(value_str); + return 0; + + } else { + + int r; + r = image_read_header_float(image, value_str, pval); + free(value_str); + + if ( r == 0 ) { + /* Value read from headers with no units */ + *pval *= default_scale; + return 0; + } else { + /* Failed to read value from headers */ + return 1; + } + } + + } else { + + /* Units are specified */ + + double scale; + + if ( strcmp(units, "mm") == 0 ) { + scale = 1e-3; + } else if ( strcmp(units, "m") == 0 ) { + scale = 1.0; + } else { + ERROR("Invalid length unit '%s'\n", units); + free(value_str); + free(units); + return 1; + } + + if ( convert_float(value_str, pval) == 0 ) { + + /* Literal value, units specified */ + free(value_str); + free(units); + *pval *= scale; + return 0; + + } else { + + int r; + r = image_read_header_float(image, value_str, pval); + free(value_str); + + if ( r == 0 ) { + /* Value read from headers, units specified */ + *pval *= scale; + return 0; + } else { + /* Failed to read value from headers */ + return 1; + } + } + } +} + + + +int create_detgeom(struct image *image, const DataTemplate *dtempl) +{ + struct detgeom *detgeom; + int i; + + if ( dtempl == NULL ) { + ERROR("NULL data template!\n"); + return 1; + } + + detgeom = malloc(sizeof(struct detgeom)); + if ( detgeom == NULL ) return 1; + + detgeom->panels = malloc(dtempl->n_panels*sizeof(struct detgeom_panel)); + if ( detgeom->panels == NULL ) { + free(detgeom); + return 1; + } + + detgeom->n_panels = dtempl->n_panels; + + for ( i=0; in_panels; i++ ) { + + struct detgeom_panel *p = &detgeom->panels[i]; + struct panel_template *tmpl = &dtempl->panels[i]; + double shift_x, shift_y; + + p->name = safe_strdup(tmpl->name); + + p->pixel_pitch = tmpl->pixel_pitch; + + /* NB cnx,cny are in pixels, cnz is in m */ + p->cnx = tmpl->cnx; + p->cny = tmpl->cny; + if ( im_get_length(image, tmpl->cnz_from, + 1e-3, &p->cnz) ) + { + ERROR("Failed to read length from '%s'\n", tmpl->cnz_from); + return 1; + } + + /* Apply offset (in m) and then convert cnz from + * m to pixels */ + p->cnz += tmpl->cnz_offset; + p->cnz /= p->pixel_pitch; + + /* Apply overall shift (already in m) */ + if ( dtempl->shift_x_from != NULL ) { + if ( im_get_length(image, dtempl->shift_x_from, 1.0, &shift_x) ) { + ERROR("Failed to read length from '%s'\n", + dtempl->shift_x_from); + return 1; + } + if ( im_get_length(image, dtempl->shift_y_from, 1.0, &shift_y) ) { + ERROR("Failed to read length from '%s'\n", + dtempl->shift_y_from); + return 1; + } + } else { + shift_x = 0.0; + shift_y = 0.0; + } + + if ( !isnan(shift_x) ) { + p->cnx += shift_x / p->pixel_pitch; + } + if ( !isnan(shift_y) ) { + p->cny += shift_y / p->pixel_pitch; + } + + p->max_adu = tmpl->max_adu; + + switch ( tmpl->adu_scale_unit ) { + + case ADU_PER_PHOTON: + p->adu_per_photon = tmpl->adu_scale; + break; + + case ADU_PER_EV: + p->adu_per_photon = tmpl->adu_scale + * ph_lambda_to_eV(image->lambda); + break; + + default: + p->adu_per_photon = 1.0; + ERROR("Invalid ADU/ph scale unit (%i)\n", + tmpl->adu_scale_unit); + break; + + } + + p->w = tmpl->orig_max_fs - tmpl->orig_min_fs + 1; + p->h = tmpl->orig_max_ss - tmpl->orig_min_ss + 1; + + p->fsx = tmpl->fsx; + p->fsy = tmpl->fsy; + p->fsz = tmpl->fsz; + p->ssx = tmpl->ssx; + p->ssy = tmpl->ssy; + p->ssz = tmpl->ssz; + + } + + image->detgeom = detgeom; + + return 0; +} diff --git a/libcrystfel/src/datatemplate_priv.h b/libcrystfel/src/datatemplate_priv.h index ff383705..7ade9ff8 100644 --- a/libcrystfel/src/datatemplate_priv.h +++ b/libcrystfel/src/datatemplate_priv.h @@ -234,5 +234,7 @@ struct _datatemplate }; extern double convert_to_m(double val, int units); +extern int create_detgeom(struct image *image, + const DataTemplate *dtempl); #endif /* DATATEMPLATE_PRIV_H */ diff --git a/libcrystfel/src/image.c b/libcrystfel/src/image.c index 3448421e..3ce8914e 100644 --- a/libcrystfel/src/image.c +++ b/libcrystfel/src/image.c @@ -486,231 +486,6 @@ static DataSourceType file_type(const char *filename) } -static int separate_value_and_units(const char *from, - char **pvalue, - char **punits) -{ - char *sp; - char *fromcpy; - char *unitscpy; - - if ( from == NULL ) return 1; - - fromcpy = strdup(from); - if ( fromcpy == NULL ) return 1; - - sp = strchr(fromcpy, ' '); - if ( sp == NULL ) { - unitscpy = NULL; - } else { - unitscpy = strdup(sp+1); - sp[0] = '\0'; - } - - *pvalue = fromcpy; - *punits = unitscpy; - return 0; -} - - -/* default_scale is a value to be used if both of the following - * conditions are met: - * - * 1. The value is a reference to image headers/metadata, - * rather than a literal number. - * 2. No units are specified in the number. - * - * This is totally horrible. Sorry. Blame history. - */ -static int im_get_length(struct image *image, const char *from, - double default_scale, double *pval) -{ - char *value_str; - char *units; - - if ( from == NULL ) return 1; - - if ( separate_value_and_units(from, &value_str, &units) ) return 1; - - if ( units == NULL ) { - - /* No units given */ - - if ( convert_float(value_str, pval) == 0 ) { - - /* Literal value with no units */ - free(value_str); - return 0; - - } else { - - int r; - r = image_read_header_float(image, value_str, pval); - free(value_str); - - if ( r == 0 ) { - /* Value read from headers with no units */ - *pval *= default_scale; - return 0; - } else { - /* Failed to read value from headers */ - return 1; - } - } - - } else { - - /* Units are specified */ - - double scale; - - if ( strcmp(units, "mm") == 0 ) { - scale = 1e-3; - } else if ( strcmp(units, "m") == 0 ) { - scale = 1.0; - } else { - ERROR("Invalid length unit '%s'\n", units); - free(value_str); - free(units); - return 1; - } - - if ( convert_float(value_str, pval) == 0 ) { - - /* Literal value, units specified */ - free(value_str); - free(units); - *pval *= scale; - return 0; - - } else { - - int r; - r = image_read_header_float(image, value_str, pval); - free(value_str); - - if ( r == 0 ) { - /* Value read from headers, units specified */ - *pval *= scale; - return 0; - } else { - /* Failed to read value from headers */ - return 1; - } - } - } -} - - -int create_detgeom(struct image *image, const DataTemplate *dtempl) -{ - struct detgeom *detgeom; - int i; - - if ( dtempl == NULL ) { - ERROR("NULL data template!\n"); - return 1; - } - - detgeom = malloc(sizeof(struct detgeom)); - if ( detgeom == NULL ) return 1; - - detgeom->panels = malloc(dtempl->n_panels*sizeof(struct detgeom_panel)); - if ( detgeom->panels == NULL ) { - free(detgeom); - return 1; - } - - detgeom->n_panels = dtempl->n_panels; - - for ( i=0; in_panels; i++ ) { - - struct detgeom_panel *p = &detgeom->panels[i]; - struct panel_template *tmpl = &dtempl->panels[i]; - double shift_x, shift_y; - - p->name = safe_strdup(tmpl->name); - - p->pixel_pitch = tmpl->pixel_pitch; - - /* NB cnx,cny are in pixels, cnz is in m */ - p->cnx = tmpl->cnx; - p->cny = tmpl->cny; - if ( im_get_length(image, tmpl->cnz_from, - 1e-3, &p->cnz) ) - { - ERROR("Failed to read length from '%s'\n", tmpl->cnz_from); - return 1; - } - - /* Apply offset (in m) and then convert cnz from - * m to pixels */ - p->cnz += tmpl->cnz_offset; - p->cnz /= p->pixel_pitch; - - /* Apply overall shift (already in m) */ - if ( dtempl->shift_x_from != NULL ) { - if ( im_get_length(image, dtempl->shift_x_from, 1.0, &shift_x) ) { - ERROR("Failed to read length from '%s'\n", - dtempl->shift_x_from); - return 1; - } - if ( im_get_length(image, dtempl->shift_y_from, 1.0, &shift_y) ) { - ERROR("Failed to read length from '%s'\n", - dtempl->shift_y_from); - return 1; - } - } else { - shift_x = 0.0; - shift_y = 0.0; - } - - if ( !isnan(shift_x) ) { - p->cnx += shift_x / p->pixel_pitch; - } - if ( !isnan(shift_y) ) { - p->cny += shift_y / p->pixel_pitch; - } - - p->max_adu = tmpl->max_adu; - - switch ( tmpl->adu_scale_unit ) { - - case ADU_PER_PHOTON: - p->adu_per_photon = tmpl->adu_scale; - break; - - case ADU_PER_EV: - p->adu_per_photon = tmpl->adu_scale - * ph_lambda_to_eV(image->lambda); - break; - - default: - p->adu_per_photon = 1.0; - ERROR("Invalid ADU/ph scale unit (%i)\n", - tmpl->adu_scale_unit); - break; - - } - - p->w = tmpl->orig_max_fs - tmpl->orig_min_fs + 1; - p->h = tmpl->orig_max_ss - tmpl->orig_min_ss + 1; - - p->fsx = tmpl->fsx; - p->fsy = tmpl->fsy; - p->fsz = tmpl->fsz; - p->ssx = tmpl->ssx; - p->ssy = tmpl->ssy; - p->ssz = tmpl->ssz; - - } - - image->detgeom = detgeom; - - return 0; -} - - int image_set_zero_data(struct image *image, const DataTemplate *dtempl) { diff --git a/libcrystfel/src/image.h b/libcrystfel/src/image.h index 3f1ed219..10ae0ba7 100644 --- a/libcrystfel/src/image.h +++ b/libcrystfel/src/image.h @@ -270,9 +270,6 @@ extern int image_write(const struct image *image, const DataTemplate *dtempl, const char *filename); -/* Use within libcrystfel only */ -extern int create_detgeom(struct image *image, - const DataTemplate *dtempl); #ifdef __cplusplus } diff --git a/libcrystfel/src/stream.c b/libcrystfel/src/stream.c index 95345f40..960d0a02 100644 --- a/libcrystfel/src/stream.c +++ b/libcrystfel/src/stream.c @@ -49,6 +49,7 @@ #include "reflist.h" #include "reflist-utils.h" #include "datatemplate.h" +#include "datatemplate_priv.h" #include "detgeom.h" #include "libcrystfel-version.h" -- cgit v1.2.3 From 7c42f8e2b675e017ab1144ca38c9e74c24d68266 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 21 Sep 2022 16:32:20 +0200 Subject: create_detgeom: Return detgeom structure rather than altering image argument --- libcrystfel/src/datatemplate.c | 23 +++++++++++------------ libcrystfel/src/datatemplate_priv.h | 6 ++++-- libcrystfel/src/image.c | 7 ++++--- libcrystfel/src/stream.c | 3 ++- 4 files changed, 21 insertions(+), 18 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/datatemplate.c b/libcrystfel/src/datatemplate.c index 94a8de2a..8f909eb8 100644 --- a/libcrystfel/src/datatemplate.c +++ b/libcrystfel/src/datatemplate.c @@ -1853,24 +1853,24 @@ static int im_get_length(struct image *image, const char *from, } - -int create_detgeom(struct image *image, const DataTemplate *dtempl) +struct detgeom *create_detgeom(struct image *image, + const DataTemplate *dtempl) { struct detgeom *detgeom; int i; if ( dtempl == NULL ) { ERROR("NULL data template!\n"); - return 1; + return NULL; } detgeom = malloc(sizeof(struct detgeom)); - if ( detgeom == NULL ) return 1; + if ( detgeom == NULL ) return NULL; detgeom->panels = malloc(dtempl->n_panels*sizeof(struct detgeom_panel)); if ( detgeom->panels == NULL ) { free(detgeom); - return 1; + return NULL; } detgeom->n_panels = dtempl->n_panels; @@ -1888,11 +1888,10 @@ int create_detgeom(struct image *image, const DataTemplate *dtempl) /* NB cnx,cny are in pixels, cnz is in m */ p->cnx = tmpl->cnx; p->cny = tmpl->cny; - if ( im_get_length(image, tmpl->cnz_from, - 1e-3, &p->cnz) ) + if ( im_get_length(image, tmpl->cnz_from, 1e-3, &p->cnz) ) { ERROR("Failed to read length from '%s'\n", tmpl->cnz_from); - return 1; + return NULL; } /* Apply offset (in m) and then convert cnz from @@ -1905,12 +1904,12 @@ int create_detgeom(struct image *image, const DataTemplate *dtempl) if ( im_get_length(image, dtempl->shift_x_from, 1.0, &shift_x) ) { ERROR("Failed to read length from '%s'\n", dtempl->shift_x_from); - return 1; + return NULL; } if ( im_get_length(image, dtempl->shift_y_from, 1.0, &shift_y) ) { ERROR("Failed to read length from '%s'\n", dtempl->shift_y_from); - return 1; + return NULL; } } else { shift_x = 0.0; @@ -1957,7 +1956,7 @@ int create_detgeom(struct image *image, const DataTemplate *dtempl) } - image->detgeom = detgeom; + return detgeom; +} - return 0; } diff --git a/libcrystfel/src/datatemplate_priv.h b/libcrystfel/src/datatemplate_priv.h index 7ade9ff8..3e61f216 100644 --- a/libcrystfel/src/datatemplate_priv.h +++ b/libcrystfel/src/datatemplate_priv.h @@ -32,6 +32,8 @@ #ifndef DATATEMPLATE_PRIV_H #define DATATEMPLATE_PRIV_H +#include "detgeom.h" + /* Maximum number of dimensions expected in data files */ #define MAX_DIMS (16) @@ -234,7 +236,7 @@ struct _datatemplate }; extern double convert_to_m(double val, int units); -extern int create_detgeom(struct image *image, - const DataTemplate *dtempl); +extern struct detgeom *create_detgeom(struct image *image, + const DataTemplate *dtempl); #endif /* DATATEMPLATE_PRIV_H */ diff --git a/libcrystfel/src/image.c b/libcrystfel/src/image.c index 3ce8914e..45b61bec 100644 --- a/libcrystfel/src/image.c +++ b/libcrystfel/src/image.c @@ -1074,7 +1074,8 @@ struct image *image_create_for_simulation(const DataTemplate *dtempl) return NULL; } - if ( create_detgeom(image, dtempl) ) { + image->detgeom = create_detgeom(image, dtempl); + if ( image->detgeom == NULL ) { image_free(image); return NULL; } @@ -1125,9 +1126,9 @@ static int do_image_read(struct image *image, const DataTemplate *dtempl, } profile_start("create-detgeom"); - r = create_detgeom(image, dtempl); + image->detgeom = create_detgeom(image, dtempl); profile_end("create-detgeom"); - if ( r ) { + if ( image->detgeom == NULL ) { ERROR("Failed to read geometry information\n"); return 1; } diff --git a/libcrystfel/src/stream.c b/libcrystfel/src/stream.c index 960d0a02..3a4f5dae 100644 --- a/libcrystfel/src/stream.c +++ b/libcrystfel/src/stream.c @@ -897,7 +897,8 @@ struct image *stream_read_chunk(Stream *st, StreamFlags srf) if ( have_filename && have_ev ) { /* Success */ if ( srf & STREAM_DATA_DETGEOM ) { - if ( create_detgeom(image, st->dtempl_read) ) { + image->detgeom = create_detgeom(image, st->dtempl_read); + if ( image->detgeom == NULL ) { image_free(image); return NULL; } -- cgit v1.2.3 From e8053e5386146bf88b2b45d028774da9f4439069 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Thu, 22 Sep 2022 14:07:34 +0200 Subject: Implement data_template_get_2d_detgeom_if_possible It is horrible. But it's contained inside the DataTemplate module. --- libcrystfel/src/datatemplate.c | 134 +++++++++++++++++++++++++++++++++++- libcrystfel/src/datatemplate_priv.h | 3 +- libcrystfel/src/image.c | 11 ++- libcrystfel/src/stream.c | 2 +- 4 files changed, 138 insertions(+), 12 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/datatemplate.c b/libcrystfel/src/datatemplate.c index 8f909eb8..1075b349 100644 --- a/libcrystfel/src/datatemplate.c +++ b/libcrystfel/src/datatemplate.c @@ -1853,8 +1853,108 @@ static int im_get_length(struct image *image, const char *from, } +static int all_panels_reference_same_clen(const DataTemplate *dtempl) +{ + int i; + char *first_val = NULL; + char *first_units = NULL; + int fail = 0; + + for ( i=0; in_panels; i++ ) { + struct panel_template *p = &dtempl->panels[i]; + char *val; + char *units; + if ( separate_value_and_units(p->cnz_from, &val, &units) ) { + /* Parse error */ + return 0; + } + if ( i == 0 ) { + first_val = val; + first_units = units; + } else { + if ( strcmp(val, first_val) != 0 ) fail = 1; + if ( strcmp(units, first_units) != 0 ) fail = 1; + free(val); + free(units); + } + } + + free(first_val); + free(first_units); + return fail; +} + + +static int all_coffsets_small(const DataTemplate *dtempl) +{ + int i; + + for ( i=0; in_panels; i++ ) { + struct panel_template *p = &dtempl->panels[i]; + if ( p->cnz_offset > 10.0*p->pixel_pitch ) return 0; + } + + return 1; +} + + +static int all_panels_same_clen(const DataTemplate *dtempl) +{ + int i; + double *zvals; + double total = 0.0; + double mean; + + zvals = malloc(sizeof(double)*dtempl->n_panels); + if ( zvals == NULL ) return 0; + + for ( i=0; in_panels; i++ ) { + struct panel_template *p = &dtempl->panels[i]; + if ( im_get_length(NULL, p->cnz_from, 1e-3, &zvals[i]) ) { + /* Can't get length because it used a header reference */ + free(zvals); + return 0; + } + total += zvals[i]; + } + + mean = total/dtempl->n_panels; + for ( i=0; in_panels; i++ ) { + struct panel_template *p = &dtempl->panels[i]; + if ( fabs(zvals[i] - mean) > 10.0*p->pixel_pitch ) return 0; + } + + free(zvals); + + return 1; +} + + +static int all_panels_perpendicular_to_beam(const DataTemplate *dtempl) +{ + int i; + + for ( i=0; in_panels; i++ ) { + double z_diff; + struct panel_template *p = &dtempl->panels[i]; + z_diff = p->fsz*PANEL_WIDTH(p) + p->ssz*PANEL_HEIGHT(p); + if ( z_diff > 10.0*p->pixel_pitch ) return 0; + } + return 1; +} + + +static int detector_flat(const DataTemplate *dtempl) +{ + return all_panels_perpendicular_to_beam(dtempl) + && ( (all_panels_reference_same_clen(dtempl) && all_coffsets_small(dtempl)) + || all_panels_same_clen(dtempl) ); +} + + struct detgeom *create_detgeom(struct image *image, - const DataTemplate *dtempl) + const DataTemplate *dtempl, + int two_d_only) { struct detgeom *detgeom; int i; @@ -1875,6 +1975,12 @@ struct detgeom *create_detgeom(struct image *image, detgeom->n_panels = dtempl->n_panels; + if ( two_d_only ) { + if ( !detector_flat(dtempl) ) return NULL; + if ( dtempl->shift_x_from != NULL ) return NULL; + if ( dtempl->shift_y_from != NULL ) return NULL; + } + for ( i=0; in_panels; i++ ) { struct detgeom_panel *p = &detgeom->panels[i]; @@ -1888,10 +1994,15 @@ struct detgeom *create_detgeom(struct image *image, /* NB cnx,cny are in pixels, cnz is in m */ p->cnx = tmpl->cnx; p->cny = tmpl->cny; + if ( im_get_length(image, tmpl->cnz_from, 1e-3, &p->cnz) ) { - ERROR("Failed to read length from '%s'\n", tmpl->cnz_from); - return NULL; + if ( two_d_only ) { + p->cnz = NAN; + } else { + ERROR("Failed to read length from '%s'\n", tmpl->cnz_from); + return NULL; + } } /* Apply offset (in m) and then convert cnz from @@ -1959,4 +2070,21 @@ struct detgeom *create_detgeom(struct image *image, return detgeom; } + +/** + * Create a detgeom structure from the DataTemplate, if possible, and ignoring + * 3D information. + * + * This procedure will create a detgeom structure provided that the detector + * is close to lying in a single flat plane perpendicular to the beam + * direction. If certain things (e.g. panel z-positions) refer to headers, + * it might not be possible to determine that the detector is really flat + * until an image is loaded. Therefore you must gracefully handle a NULL + * return value from this routine. + * + * \returns the detgeom structure, or NULL if impossible. + */ +struct detgeom *data_template_get_2d_detgeom_if_possible(const DataTemplate *dt) +{ + return create_detgeom(NULL, dt, 1); } diff --git a/libcrystfel/src/datatemplate_priv.h b/libcrystfel/src/datatemplate_priv.h index 3e61f216..26b90d6b 100644 --- a/libcrystfel/src/datatemplate_priv.h +++ b/libcrystfel/src/datatemplate_priv.h @@ -237,6 +237,7 @@ struct _datatemplate extern double convert_to_m(double val, int units); extern struct detgeom *create_detgeom(struct image *image, - const DataTemplate *dtempl); + const DataTemplate *dtempl, + int two_d_only); #endif /* DATATEMPLATE_PRIV_H */ diff --git a/libcrystfel/src/image.c b/libcrystfel/src/image.c index 45b61bec..84088ead 100644 --- a/libcrystfel/src/image.c +++ b/libcrystfel/src/image.c @@ -412,11 +412,6 @@ static struct header_cache_entry *cached_header(struct image *image, const char { struct header_cache_entry *ce; - if ( image == NULL ) { - ERROR("Attempt to retrieve a header value without an image\n"); - return NULL; - } - ce = find_cache_entry(image, from); if ( ce != NULL ) return ce; @@ -434,6 +429,8 @@ int image_read_header_float(struct image *image, const char *from, double *val) { struct header_cache_entry *ce; + if ( image == NULL ) return 1; + ce = cached_header(image, from); if ( ce == NULL ) return 1; @@ -1074,7 +1071,7 @@ struct image *image_create_for_simulation(const DataTemplate *dtempl) return NULL; } - image->detgeom = create_detgeom(image, dtempl); + image->detgeom = create_detgeom(image, dtempl, 0); if ( image->detgeom == NULL ) { image_free(image); return NULL; @@ -1126,7 +1123,7 @@ static int do_image_read(struct image *image, const DataTemplate *dtempl, } profile_start("create-detgeom"); - image->detgeom = create_detgeom(image, dtempl); + image->detgeom = create_detgeom(image, dtempl, 0); profile_end("create-detgeom"); if ( image->detgeom == NULL ) { ERROR("Failed to read geometry information\n"); diff --git a/libcrystfel/src/stream.c b/libcrystfel/src/stream.c index 3a4f5dae..0195e0ff 100644 --- a/libcrystfel/src/stream.c +++ b/libcrystfel/src/stream.c @@ -897,7 +897,7 @@ struct image *stream_read_chunk(Stream *st, StreamFlags srf) if ( have_filename && have_ev ) { /* Success */ if ( srf & STREAM_DATA_DETGEOM ) { - image->detgeom = create_detgeom(image, st->dtempl_read); + image->detgeom = create_detgeom(image, st->dtempl_read, 0); if ( image->detgeom == NULL ) { image_free(image); return NULL; -- cgit v1.2.3 From 369a951b5bf6bb94fbffb2002281fa4d1a67506d Mon Sep 17 00:00:00 2001 From: Thomas White Date: Thu, 22 Sep 2022 14:32:15 +0200 Subject: Add missing prototype --- libcrystfel/src/datatemplate.h | 2 ++ 1 file changed, 2 insertions(+) (limited to 'libcrystfel') diff --git a/libcrystfel/src/datatemplate.h b/libcrystfel/src/datatemplate.h index 28cabcf5..ddc210e6 100644 --- a/libcrystfel/src/datatemplate.h +++ b/libcrystfel/src/datatemplate.h @@ -96,6 +96,8 @@ extern struct rg_collection *data_template_get_rigid_groups(const DataTemplate * extern double data_template_get_wavelength_if_possible(const DataTemplate *dt); +extern struct detgeom *data_template_get_2d_detgeom_if_possible(const DataTemplate *dt); + #ifdef __cplusplus } #endif -- cgit v1.2.3 From 025f9e9c9022bb1118783a089fa5854b956eb04a Mon Sep 17 00:00:00 2001 From: Thomas White Date: Thu, 22 Sep 2022 14:42:01 +0200 Subject: all_panels_reference_same_clen: Handle missing units --- libcrystfel/src/datatemplate.c | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/datatemplate.c b/libcrystfel/src/datatemplate.c index 1075b349..bbcf1aa5 100644 --- a/libcrystfel/src/datatemplate.c +++ b/libcrystfel/src/datatemplate.c @@ -1853,6 +1853,14 @@ static int im_get_length(struct image *image, const char *from, } +static int safe_strcmp(const char *a, const char *b) +{ + if ( (a==NULL) && (b==NULL) ) return 0; + if ( (a!=NULL) && (b!=NULL) ) return strcmp(a, b); + return 1; +} + + static int all_panels_reference_same_clen(const DataTemplate *dtempl) { int i; @@ -1872,8 +1880,8 @@ static int all_panels_reference_same_clen(const DataTemplate *dtempl) first_val = val; first_units = units; } else { - if ( strcmp(val, first_val) != 0 ) fail = 1; - if ( strcmp(units, first_units) != 0 ) fail = 1; + if ( safe_strcmp(val, first_val) != 0 ) fail = 1; + if ( safe_strcmp(units, first_units) != 0 ) fail = 1; free(val); free(units); } -- cgit v1.2.3 From 46d4aab96fd1703a0c4d7ea36b6c1cfd315857d9 Mon Sep 17 00:00:00 2001 From: Alexandra Tolstikova Date: Thu, 22 Sep 2022 17:02:37 +0200 Subject: Add fast mode for peakfinder8 --- libcrystfel/src/peakfinder8.c | 297 ++++++++++++++++++++++++++++++++++++++---- libcrystfel/src/peakfinder8.h | 14 +- libcrystfel/src/peaks.c | 12 +- libcrystfel/src/peaks.h | 3 +- 4 files changed, 294 insertions(+), 32 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/peakfinder8.c b/libcrystfel/src/peakfinder8.c index 537c2365..3420a7a7 100644 --- a/libcrystfel/src/peakfinder8.c +++ b/libcrystfel/src/peakfinder8.c @@ -34,6 +34,7 @@ #include #include #include +#include #include "peakfinder8.h" #include "detgeom.h" @@ -43,13 +44,24 @@ /** \file peakfinder8.h */ // CrystFEL-only block 1 + struct radius_maps { float **r_maps; + int *n_pixels; int n_rmaps; }; +struct radial_stats_pixels +{ + int n_panels; + int *n_pixels; // n_pixels[panel] + int **pidx; // pixel_index[panel][0..n_pixels] + int **radius; // pixel_radius[panel][0..n_pixels] +}; + + struct peakfinder_mask { char **masks; @@ -101,7 +113,138 @@ struct peakfinder_peak_data }; -// CrystFEL-only block 2 +static struct radial_stats_pixels *compute_rstats_pixels(struct radius_maps *rmaps) +{ + int p; + int i; + + struct radial_stats_pixels *rsp = NULL; + rsp = (struct radial_stats_pixels *)malloc(sizeof(struct radial_stats_pixels)); + if ( rsp == NULL ) { + return NULL; + } + rsp->n_pixels = (int *)malloc(rmaps->n_rmaps * sizeof(int)); + if ( rsp->n_pixels == NULL ) { + free(rsp); + return NULL; + } + rsp->pidx = (int **)malloc(rmaps->n_rmaps * sizeof(int *)); + if ( rsp->pidx == NULL ) { + free(rsp->n_pixels); + free(rsp); + return NULL; + } + rsp->radius = (int **)malloc(rmaps->n_rmaps * sizeof(int *)); + if ( rsp->radius == NULL ) { + free(rsp->n_pixels); + free(rsp->pidx); + free(rsp); + return NULL; + } + srand(0); + + int n_pixels_per_bin = 100; // Can make this a parameter + + // Assuming 5000 is the maximum possible radius + int n_pixels[5000] = {0}; // selected pixels per bin + int n_tot_pixels[5000] = {0}; // total pixels per bin + int panel[5000][n_pixels_per_bin]; // panel ID of selected pixels + int idx[5000][n_pixels_per_bin]; // index of selected pixels + int radius; + + for ( p = 0; p < rmaps->n_rmaps; p++ ) { + rsp->n_pixels[p] = 0; + for ( i = 0; i < rmaps->n_pixels[p]; i++ ) { + // Reservoir sampling: + radius = (int)rint(rmaps->r_maps[p][i]); + n_tot_pixels[radius] += 1; + + if ( n_pixels[radius] < n_pixels_per_bin ) { + panel[radius][n_pixels[radius]] = p; + idx[radius][n_pixels[radius]] = i; + + n_pixels[radius] += 1; + rsp->n_pixels[p] += 1; + } else { + int rand_i = rand() % n_tot_pixels[radius]; + if ( rand_i < n_pixels_per_bin ) { + rsp->n_pixels[panel[radius][rand_i]] -= 1; + rsp->n_pixels[p] += 1; + + panel[radius][rand_i] = p; + idx[radius][rand_i] = i; + } + } + } + } + + int *sidx = (int *)malloc(rmaps->n_rmaps * sizeof(int)); + if ( sidx == NULL ) { + free(rsp->n_pixels); + free(rsp->pidx); + free(rsp->radius); + free(rsp); + return NULL; + } + for ( p = 0; p < rmaps->n_rmaps; p++ ) { + rsp->pidx[p] = (int *)malloc(rsp->n_pixels[p] * sizeof(int)); + if ( rsp->pidx[p] == NULL ) { + for ( i = 0; i < p; i++ ) { + free(rsp->pidx[i]); + free(rsp->radius[i]); + } + free(rsp->pidx); + free(rsp->radius); + free(rsp->n_pixels); + free(rsp); + free(sidx); + return NULL; + } + rsp->radius[p] = (int *)malloc(rsp->n_pixels[p] * sizeof(int)); + if ( rsp->radius[p] == NULL ) { + for ( i = 0; i < p; i++ ) { + free(rsp->pidx[i]); + free(rsp->radius[i]); + } + free(rsp->pidx[p]); + free(rsp->pidx); + free(rsp->radius); + free(rsp->n_pixels); + free(rsp); + free(sidx); + return NULL; + } + sidx[p] = 0; + } + + for ( radius = 0; radius < 5000; radius++ ) { + for ( i = 0; i < n_pixels[radius]; i++ ) { + p = panel[radius][i]; + rsp->pidx[p][sidx[p]] = idx[radius][i]; + rsp->radius[p][sidx[p]] = radius; + sidx[p] += 1; + } + } + free(sidx); + + rsp->n_panels = rmaps->n_rmaps; + return rsp; +} + +static void free_rstats_pixels(struct radial_stats_pixels *rsp) +{ + int i; + for ( i = 0; i < rsp->n_panels; i++ ) { + free(rsp->pidx[i]); + free(rsp->radius[i]); + } + free(rsp->pidx); + free(rsp->radius); + free(rsp->n_pixels); + free(rsp); +} + + static struct radius_maps *compute_radius_maps(struct detgeom *det) { int i, u, iss, ifs; @@ -119,6 +262,12 @@ static struct radius_maps *compute_radius_maps(struct detgeom *det) return NULL; } + rm->n_pixels = (int *)malloc(det->n_panels*sizeof(int*)); + if ( rm->r_maps == NULL ) { + free(rm); + return NULL; + } + rm->n_rmaps = det->n_panels; for( i=0 ; in_panels ; i++ ) { @@ -133,7 +282,7 @@ static struct radius_maps *compute_radius_maps(struct detgeom *det) free(rm); return NULL; } - + rm->n_pixels[i] = p.h * p.w; for ( iss=0 ; issr_maps[i]); } free(r_maps->r_maps); + free(r_maps->n_pixels); free(r_maps); } +// CrystFEL-only block 2 +struct pf8_private_data *prepare_peakfinder8(struct detgeom *det, int fast_mode) +{ + struct pf8_private_data *data = NULL; + if ( det == NULL ) { + return NULL; + } + + data = (struct pf8_private_data *)malloc(sizeof(struct pf8_private_data)); + if ( data == NULL ) { + return NULL; + } + data->rmaps = compute_radius_maps(det); + if ( data->rmaps == NULL ) { + free(data); + return NULL; + } + if ( fast_mode ) { + data->rpixels = compute_rstats_pixels(data->rmaps); + if ( data->rpixels == NULL ) { + free_radius_maps(data->rmaps); + free(data); + return NULL; + } + } else { + data->rpixels = NULL; + } + data->fast_mode = fast_mode; + return data; +} + + +void free_pf8_private_data(struct pf8_private_data *data) +{ + free_radius_maps(data->rmaps); + if ( data->fast_mode ) { + free_rstats_pixels(data->rpixels); + } + free(data); +} + + static struct peakfinder_mask *create_peakfinder_mask(struct image *img, struct radius_maps *rmps, int min_res, @@ -385,6 +577,31 @@ static void fill_radial_bins(float *data, } } +static void fill_radial_bins_fast(float *data, int w, int h, int n_pixels, + int *pidx, int *radius, char *mask, + float *rthreshold, float *lthreshold, + float *roffset, float *rsigma, int *rcount) +{ + int i; + + int curr_r; + float value; + + for (i = 0; i < n_pixels; i++) + { + if (mask[pidx[i]] != 0) + { + curr_r = radius[i]; + value = data[pidx[i]]; + if (value < rthreshold[curr_r] && value > lthreshold[curr_r]) + { + roffset[curr_r] += value; + rsigma[curr_r] += (value * value); + rcount[curr_r] += 1; + } + } + } +} static void compute_radial_stats(float *rthreshold, float *lthreshold, @@ -1021,9 +1238,13 @@ int peakfinder8(struct image *img, int max_n_peaks, float threshold, float min_snr, int min_pix_count, int max_pix_count, int local_bg_radius, int min_res, - int max_res, int use_saturated) + int max_res, int use_saturated, + int fast_mode, struct pf8_private_data *private_data) { + struct pf8_private_data *geomdata; struct radius_maps *rmaps; + struct radial_stats_pixels *rspixels; + struct peakfinder_mask *pfmask; struct peakfinder_panel_data *pfdata; struct radial_stats *rstats; @@ -1040,18 +1261,28 @@ int peakfinder8(struct image *img, int max_n_peaks, if ( img->detgeom == NULL) return 1; - rmaps = compute_radius_maps(img->detgeom); - if ( rmaps == NULL ) return 1; + profile_start("pf8-rmaps"); + if ( private_data == NULL ) { + geomdata = prepare_peakfinder8(img->detgeom, fast_mode); + } else { + geomdata = private_data; + } + rmaps = geomdata->rmaps; + rspixels = geomdata->rpixels; + profile_end("pf8-rmaps"); + if (geomdata == NULL) return 1; + profile_start("pf8-mask"); pfmask = create_peakfinder_mask(img, rmaps, min_res, max_res); + profile_end("pf8-mask"); if ( pfmask == NULL ) { - free_radius_maps(rmaps); + if ( private_data == NULL ) free_pf8_private_data(geomdata); return 1; } pfdata = allocate_panel_data(img->detgeom->n_panels); if ( pfdata == NULL) { - free_radius_maps(rmaps); + if ( private_data == NULL ) free_pf8_private_data(geomdata); free_peakfinder_mask(pfmask); return 1; } @@ -1077,7 +1308,7 @@ int peakfinder8(struct image *img, int max_n_peaks, rstats = allocate_radial_stats(num_rad_bins); if ( rstats == NULL ) { - free_radius_maps(rmaps); + if ( private_data == NULL ) free_pf8_private_data(geomdata); free_peakfinder_mask(pfmask); free_panel_data(pfdata); return 1; @@ -1087,7 +1318,7 @@ int peakfinder8(struct image *img, int max_n_peaks, rstats->rthreshold[i] = 1e9; rstats->lthreshold[i] = -1e9; } - + profile_start("pf8-rstats"); for ( it_counter=0 ; it_counternum_panels ; pi++ ) { - - fill_radial_bins(pfdata->panel_data[pi], - pfdata->panel_w[pi], - pfdata->panel_h[pi], - rmaps->r_maps[pi], - pfmask->masks[pi], - rstats->rthreshold, - rstats->lthreshold, - rstats->roffset, - rstats->rsigma, - rstats->rcount); - + if ( fast_mode ) { + fill_radial_bins_fast(pfdata->panel_data[pi], + pfdata->panel_w[pi], + pfdata->panel_h[pi], + rspixels->n_pixels[pi], + rspixels->pidx[pi], + rspixels->radius[pi], + pfmask->masks[pi], + rstats->rthreshold, + rstats->lthreshold, + rstats->roffset, + rstats->rsigma, + rstats->rcount); + } else { + fill_radial_bins(pfdata->panel_data[pi], + pfdata->panel_w[pi], + pfdata->panel_h[pi], + rmaps->r_maps[pi], + pfmask->masks[pi], + rstats->rthreshold, + rstats->lthreshold, + rstats->roffset, + rstats->rsigma, + rstats->rcount); + } } compute_radial_stats(rstats->rthreshold, @@ -1121,10 +1365,11 @@ int peakfinder8(struct image *img, int max_n_peaks, threshold); } + profile_end("pf8-rstats"); pkdata = allocate_peak_data(max_n_peaks); if ( pkdata == NULL ) { - free_radius_maps(rmaps); + if ( private_data == NULL ) free_pf8_private_data(geomdata); free_peakfinder_mask(pfmask); free_panel_data(pfdata); free_radial_stats(rstats); @@ -1132,7 +1377,7 @@ int peakfinder8(struct image *img, int max_n_peaks, } remaining_max_num_peaks = max_n_peaks; - + profile_start("pf8-search"); for ( pi=0 ; pidetgeom->n_panels ; pi++) { int peaks_to_add; @@ -1165,10 +1410,11 @@ int peakfinder8(struct image *img, int max_n_peaks, NULL); if ( ret != 0 ) { - free_radius_maps(rmaps); + if ( private_data == NULL ) free_pf8_private_data(geomdata); free_peakfinder_mask(pfmask); free_panel_data(pfdata); free_radial_stats(rstats); + profile_end("pf8-search"); return 1; } @@ -1199,8 +1445,9 @@ int peakfinder8(struct image *img, int max_n_peaks, NULL); } } + profile_end("pf8-search"); - free_radius_maps(rmaps); + if ( private_data == NULL ) free_pf8_private_data(geomdata); free_peakfinder_mask(pfmask); free_panel_data(pfdata); free_radial_stats(rstats); diff --git a/libcrystfel/src/peakfinder8.h b/libcrystfel/src/peakfinder8.h index 28212e49..e5e86038 100644 --- a/libcrystfel/src/peakfinder8.h +++ b/libcrystfel/src/peakfinder8.h @@ -43,11 +43,23 @@ extern "C" { * The "peakfinder8" peak search algorithm, originally implemented in Cheetah */ +struct pf8_private_data +{ + int fast_mode; + struct radius_maps *rmaps; + struct radial_stats_pixels *rpixels; +}; + +struct pf8_private_data *prepare_peakfinder8(struct detgeom *det, int fast_mode); + +void free_pf8_private_data(struct pf8_private_data *data); + extern int peakfinder8(struct image *img, int max_n_peaks, float threshold, float min_snr, int mix_pix_count, int max_pix_count, int local_bg_radius, int min_res, - int max_res, int use_saturated); + int max_res, int use_saturated, + int fast_mode, struct pf8_private_data *private_data); #ifdef __cplusplus } diff --git a/libcrystfel/src/peaks.c b/libcrystfel/src/peaks.c index d260b66d..57d2c92e 100644 --- a/libcrystfel/src/peaks.c +++ b/libcrystfel/src/peaks.c @@ -454,10 +454,11 @@ void search_peaks(struct image *image, float threshold, float min_sq_gradient, * the actual \ref peakfinder8 function, found in \ref peakfinder8.h. */ int search_peaks_peakfinder8(struct image *image, int max_n_peaks, - float threshold, float min_snr, - int min_pix_count, int max_pix_count, - int local_bg_radius, int min_res, - int max_res, int use_saturated) + float threshold, float min_snr, + int min_pix_count, int max_pix_count, + int local_bg_radius, int min_res, + int max_res, int use_saturated, + int fast_mode, void *private_data) { if ( image->features != NULL ) { image_feature_list_free(image->features); @@ -467,7 +468,8 @@ int search_peaks_peakfinder8(struct image *image, int max_n_peaks, return peakfinder8(image, max_n_peaks, threshold, min_snr, min_pix_count, max_pix_count, local_bg_radius, min_res, - max_res, use_saturated); + max_res, use_saturated, + fast_mode, private_data); } diff --git a/libcrystfel/src/peaks.h b/libcrystfel/src/peaks.h index 2f100db2..d6d5c1f9 100644 --- a/libcrystfel/src/peaks.h +++ b/libcrystfel/src/peaks.h @@ -73,7 +73,8 @@ extern int search_peaks_peakfinder8(struct image *image, int max_n_peaks, float threshold, float min_snr, int mix_pix_count, int max_pix_count, int local_bg_radius, int min_res, - int max_res, int use_saturated); + int max_res, int use_saturated, + int fast_mode, void *private_data); extern int search_peaks_peakfinder9(struct image *image, float min_snr_biggest_pix, -- cgit v1.2.3