From c1b825fd1040c89114d4a8c2c5d3fbba3f946d2a Mon Sep 17 00:00:00 2001 From: Thomas White Date: Tue, 7 Feb 2017 15:06:48 +0100 Subject: Fix mask path placeholder check The mask paths for all panels have to have the same number of placeholders, but the masks do not have to have the same number of placeholders as the panel data blocks. This also tidies up a few excess strdup() calls, and removes partial_event_substitution() because retrieve_full_path() can now handle the number of placeholders being too small. --- libcrystfel/src/detector.c | 23 ++++++++++++++----- libcrystfel/src/events.c | 54 ++++++++++++++++++++------------------------- libcrystfel/src/events.h | 1 - libcrystfel/src/hdf5-file.c | 4 ++-- 4 files changed, 43 insertions(+), 39 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/detector.c b/libcrystfel/src/detector.c index 42ffbe40..8c98f79e 100644 --- a/libcrystfel/src/detector.c +++ b/libcrystfel/src/detector.c @@ -1204,7 +1204,7 @@ struct detector *get_detector_geometry_2(const char *filename, int i; int rgi, rgci; int reject = 0; - int path_dim; + int path_dim, mask_path_dim; int dim_dim; int dim_reject = 0; int dim_dim_reject = 0; @@ -1389,6 +1389,7 @@ struct detector *get_detector_geometry_2(const char *filename, } + mask_path_dim = -1; for ( i=0; in_panels; i++ ) { int panel_mask_dim = 0; @@ -1398,8 +1399,7 @@ struct detector *get_detector_geometry_2(const char *filename, next_instance = det->panels[i].mask; - while(next_instance) - { + while ( next_instance ) { next_instance = strstr(next_instance, "%"); if ( next_instance != NULL ) { next_instance += 1*sizeof(char); @@ -1407,18 +1407,29 @@ struct detector *get_detector_geometry_2(const char *filename, } } - if ( panel_mask_dim != path_dim ) { - dim_reject = 1; + if ( mask_path_dim == -1 ) { + mask_path_dim = panel_mask_dim; + } else { + if ( panel_mask_dim != mask_path_dim ) { + dim_reject = 1; + } } + } } - if ( dim_reject == 1) { + if ( dim_reject == 1 ) { ERROR("All panels' data and mask entries must have the same " "number of placeholders\n"); reject = 1; } + if ( mask_path_dim > path_dim ) { + ERROR("Number of placeholders in mask cannot be larger than " + "for data\n"); + reject = 1; + } + det->path_dim = path_dim; dim_dim_reject = 0; diff --git a/libcrystfel/src/events.c b/libcrystfel/src/events.c index 25771a69..0f170bb5 100644 --- a/libcrystfel/src/events.c +++ b/libcrystfel/src/events.c @@ -3,10 +3,11 @@ * * Event properties * - * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * * Authors: + * 2017 Thomas White * 2014 Valerio Mariani * * This file is part of CrystFEL. @@ -30,6 +31,7 @@ #define _GNU_SOURCE #include "events.h" +#include "utils.h" #include #include @@ -585,46 +587,38 @@ char *retrieve_full_path(struct event *ev, const char *data) { int ei ; char *return_value; + char *pholder; return_value = strdup(data); + pholder = strstr(return_value,"%"); + ei = 0; - for ( ei=0; eipath_length; ei++ ) { + while ( pholder != NULL ) { char *tmp; - tmp = event_path_placeholder_subst(ev->path_entries[ei], - return_value); - - free(return_value); - return_value = tmp; - - } - - return return_value; - -} - - -char *partial_event_substitution(struct event *ev, const char *data) -{ - int ei ; - char *return_value; - char *pholder; - return_value = strdup(data); - pholder = strstr(return_value,"%"); - ei = 0; + /* Check we have enough things to put in the placeholders */ + if ( ei >= ev->path_length ) { + ERROR("Too many placeholders ('%%') in location.\n"); + return NULL; + } - while( pholder != NULL) { + /* Substitute one placeholder */ + tmp = event_path_placeholder_subst(ev->path_entries[ei++], + return_value); - char *tmp_subst_data; + if ( tmp == NULL ) { + ERROR("Couldn't substitute placeholder\n"); + return NULL; + } - tmp_subst_data = event_path_placeholder_subst(ev->path_entries[ei], - return_value); + /* Next time, we will substitute the next part of the path into + * the partially substituted string */ free(return_value); - return_value = strdup(tmp_subst_data); - free(tmp_subst_data); + return_value = tmp; + pholder = strstr(return_value, "%"); - ei += 1; + } return return_value; diff --git a/libcrystfel/src/events.h b/libcrystfel/src/events.h index 7f9c6731..4f717209 100644 --- a/libcrystfel/src/events.h +++ b/libcrystfel/src/events.h @@ -78,7 +78,6 @@ extern char *get_event_string(struct event *ev); extern struct event *get_event_from_event_string(const char *ev_string); extern char *event_path_placeholder_subst(const char *ev_name, const char *data); -extern char *partial_event_substitution(struct event *ev, const char *data); extern char *retrieve_full_path(struct event *ev, const char *data); diff --git a/libcrystfel/src/hdf5-file.c b/libcrystfel/src/hdf5-file.c index 6542b360..9bcdb1b7 100644 --- a/libcrystfel/src/hdf5-file.c +++ b/libcrystfel/src/hdf5-file.c @@ -1154,7 +1154,7 @@ static int get_ev_based_value(struct hdfile *f, const char *name, char *subst_name = NULL; if ( ev->path_length != 0 ) { - subst_name = partial_event_substitution(ev, name); + subst_name = retrieve_full_path(ev, name); } else { subst_name = strdup(name); } @@ -1966,7 +1966,7 @@ char *hdfile_get_string_value(struct hdfile *f, const char *name, char *tmp = NULL, *subst_name = NULL; if (ev != NULL && ev->path_length != 0 ) { - subst_name = partial_event_substitution(ev, name); + subst_name = retrieve_full_path(ev, name); } else { subst_name = strdup(name); } -- cgit v1.2.3 From e3c6d8e8aec84f708e1098371ecb4621ad17223f Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 3 Mar 2017 13:02:51 +0100 Subject: Return error code if cell can't be inverted --- libcrystfel/src/cell.c | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/cell.c b/libcrystfel/src/cell.c index ec591e24..3de61073 100644 --- a/libcrystfel/src/cell.c +++ b/libcrystfel/src/cell.c @@ -460,10 +460,12 @@ int cell_get_parameters(UnitCell *cell, double *a, double *b, double *c, case CELL_REP_RECIP: /* Convert reciprocal -> crystallographic. * Start by converting reciprocal -> cartesian */ - cell_invert(cell->axs, cell->ays, cell->azs, - cell->bxs, cell->bys, cell->bzs, - cell->cxs, cell->cys, cell->czs, - &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz); + if ( cell_invert(cell->axs, cell->ays, cell->azs, + cell->bxs, cell->bys, cell->bzs, + cell->cxs, cell->cys, cell->czs, + &ax, &ay, &az, + &bx, &by, &bz, + &cx, &cy, &cz) ) return 1; /* Now convert cartesian -> crystallographic */ *a = modulus(ax, ay, az); -- cgit v1.2.3 From 61565336125a999790fb4c36219e9c46c5eb30cc Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 3 Mar 2017 17:37:10 +0100 Subject: Allow indexing system to store its own data, independently of indexing methods Previously, the indexing system passed all the information on to the indexing engines and then forgot about it. That made it difficult to do things like check the indexing solution after prediction refinement, because the target unit cell was unavailable. Now, the indexing system itself can keep some information. Of course, that information includes the private pointers for the indexing engines themselves. I took the opportunity to streamline things a little bit. The caller can now set up the indexing system in one step, without having to separately parse the names of the indexing methods. The caller no longer has to keep track of a separate array of methods, instead just one structure which contains everything. --- libcrystfel/src/index.c | 358 ++++++++++++++++++++++++++++------------------- libcrystfel/src/index.h | 26 ++-- libcrystfel/src/stream.c | 8 +- 3 files changed, 232 insertions(+), 160 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/index.c b/libcrystfel/src/index.c index a2fdb861..89a98920 100644 --- a/libcrystfel/src/index.c +++ b/libcrystfel/src/index.c @@ -58,6 +58,17 @@ #include "predict-refine.h" +struct _indexingprivate +{ + UnitCell *target_cell; + float tolerance[4]; + + int n_methods; + IndexingMethod *methods; + void **engine_private; +}; + + static int debug_index(struct image *image) { Crystal *cr = crystal_new(); @@ -71,74 +82,105 @@ static int debug_index(struct image *image) } -IndexingPrivate **prepare_indexing(IndexingMethod *indm, UnitCell *cell, - struct detector *det, float *ltl, - const char *options) +static void *prepare_method(IndexingMethod *m, UnitCell *cell, + struct detector *det, float *ltl, + const char *options) { - int n; - int nm = 0; - IndexingPrivate **iprivs; + char *str; + IndexingMethod in = *m; + void *priv = NULL; - while ( indm[nm] != INDEXING_NONE ) nm++; - iprivs = malloc((nm+1) * sizeof(IndexingPrivate *)); + switch ( *m & INDEXING_METHOD_MASK ) { - for ( n=0; nengine_private = malloc((n+1) * sizeof(void *)); + + for ( i=0; iengine_private[i] = prepare_method(&methods[i], cell, + det, ltl, options); + for ( j=0; jmethods = methods; + ipriv->n_methods = n; + if ( cell != NULL ) { + ipriv->target_cell = cell_new_from_cell(cell); + } else { + ipriv->target_cell = NULL; } - iprivs[n] = NULL; + for ( i=0; i<4; i++ ) ipriv->tolerance[i] = ltl[i]; - return iprivs; + return ipriv; } -void cleanup_indexing(IndexingMethod *indms, IndexingPrivate **privs) +void cleanup_indexing(IndexingPrivate *ipriv) { - int n = 0; + int n; - if ( indms == NULL ) return; /* Nothing to do */ - if ( privs == NULL ) return; /* Nothing to do */ + if ( ipriv == NULL ) return; /* Nothing to do */ - while ( indms[n] != INDEXING_NONE ) { + for ( n=0; nn_methods; n++ ) { - switch ( indms[n] & INDEXING_METHOD_MASK ) { + switch ( ipriv->methods[n] & INDEXING_METHOD_MASK ) { case INDEXING_NONE : break; case INDEXING_DIRAX : - dirax_cleanup(privs[n]); + dirax_cleanup(ipriv->engine_private[n]); break; case INDEXING_ASDF : - asdf_cleanup(privs[n]); + asdf_cleanup(ipriv->engine_private[n]); break; case INDEXING_MOSFLM : - mosflm_cleanup(privs[n]); + mosflm_cleanup(ipriv->engine_private[n]); break; case INDEXING_XDS : - xds_cleanup(privs[n]); + xds_cleanup(ipriv->engine_private[n]); break; case INDEXING_FELIX : - felix_cleanup(privs[n]); + felix_cleanup(ipriv->engine_private[n]); break; case INDEXING_DEBUG : - free(privs[n]); + free(ipriv->engine_private[n]); break; default : ERROR("Don't know how to clean up indexing method %i\n", - indms[n]); + ipriv->methods[n]); break; } @@ -204,8 +253,10 @@ void cleanup_indexing(IndexingMethod *indms, IndexingPrivate **privs) } - free(indms); - free(privs); + free(ipriv->methods); + free(ipriv->engine_private); + cell_free(ipriv->target_cell); + free(ipriv); } @@ -234,7 +285,7 @@ void map_all_peaks(struct image *image) /* Return non-zero for "success" */ static int try_indexer(struct image *image, IndexingMethod indm, - IndexingPrivate *ipriv) + IndexingPrivate *ipriv, void *mpriv) { int i, r; int n_bad = 0; @@ -245,19 +296,19 @@ static int try_indexer(struct image *image, IndexingMethod indm, return 0; case INDEXING_DIRAX : - r = run_dirax(image, ipriv); + r = run_dirax(image, mpriv); break; case INDEXING_ASDF : - r = run_asdf(image, ipriv); + r = run_asdf(image, mpriv); break; case INDEXING_MOSFLM : - r = run_mosflm(image, ipriv); + r = run_mosflm(image, mpriv); break; case INDEXING_XDS : - r = run_xds(image, ipriv); + r = run_xds(image, mpriv); break; case INDEXING_DEBUG : @@ -265,7 +316,7 @@ static int try_indexer(struct image *image, IndexingMethod indm, break; case INDEXING_FELIX : - r = felix_index(image, ipriv); + r = felix_index(image, mpriv); break; default : @@ -421,20 +472,18 @@ static int finished_retry(IndexingMethod indm, int r, struct image *image) } } -void index_pattern(struct image *image, - IndexingMethod *indms, IndexingPrivate **iprivs) +void index_pattern(struct image *image, IndexingPrivate *ipriv) { - index_pattern_2(image, indms, iprivs, NULL); + index_pattern_2(image, ipriv, NULL); } -void index_pattern_2(struct image *image, IndexingMethod *indms, - IndexingPrivate **iprivs, int *ping) +void index_pattern_2(struct image *image, IndexingPrivate *ipriv, int *ping) { int n = 0; ImageFeatureList *orig; - if ( indms == NULL ) return; + if ( ipriv == NULL ) return; map_all_peaks(image); image->crystals = NULL; @@ -442,7 +491,7 @@ void index_pattern_2(struct image *image, IndexingMethod *indms, orig = image->features; - while ( indms[n] != INDEXING_NONE ) { + for ( n=0; nn_methods; n++ ) { int done = 0; int r; @@ -453,10 +502,11 @@ void index_pattern_2(struct image *image, IndexingMethod *indms, do { - r = try_indexer(image, indms[n], iprivs[n]); + r = try_indexer(image, ipriv->methods[n], + ipriv, ipriv->engine_private[n]); success += r; ntry++; - done = finished_retry(indms[n], r, image); + done = finished_retry(ipriv->methods[n], r, image); if ( ntry > 5 ) done = 1; if ( ping != NULL ) (*ping)++; @@ -468,11 +518,14 @@ void index_pattern_2(struct image *image, IndexingMethod *indms, * crystals with a different indexing method) */ if ( success ) break; - n++; + } + if ( n < ipriv->n_methods ) { + image->indexed_by = ipriv->methods[n]; + } else { + image->indexed_by = INDEXING_NONE; } - image->indexed_by = indms[n]; image->features = orig; } @@ -639,97 +692,118 @@ char *indexer_str(IndexingMethod indm) } -IndexingMethod *build_indexer_list(const char *str) +static IndexingMethod warn_method(const char *str) +{ + ERROR("Indexing method must contain exactly one engine name: '%s'\n", + str); + return INDEXING_ERROR; +} + + +IndexingMethod get_indm_from_string(const char *str) { int n, i; - char **methods; - IndexingMethod *list; - int nmeth = 0; + char **bits; + IndexingMethod method = INDEXING_NONE; + int have_method = 0; - n = assplode(str, ",-", &methods, ASSPLODE_NONE); - list = malloc((n+1)*sizeof(IndexingMethod)); + n = assplode(str, "-", &bits, ASSPLODE_NONE); - nmeth = -1; /* So that the first method is #0 */ for ( i=0; iindexed_by = list[0]; - free(list); } } -- cgit v1.2.3 From ad305e5a1a84a5e83edbdf30e3aba5265ca40665 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Mon, 6 Mar 2017 11:18:05 +0100 Subject: Check unit cell parameters after prediction refinement --- libcrystfel/src/index.c | 33 ++++++++++++++++++++++++++++----- 1 file changed, 28 insertions(+), 5 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/index.c b/libcrystfel/src/index.c index 89a98920..028ba5e9 100644 --- a/libcrystfel/src/index.c +++ b/libcrystfel/src/index.c @@ -334,11 +334,34 @@ static int try_indexer(struct image *image, IndexingMethod indm, crystal_set_mosaicity(cr, 0.0); /* Prediction refinement if requested */ - if ( (indm & INDEXING_REFINE) - && (refine_prediction(image, cr) != 0) ) - { - crystal_set_user_flag(cr, 1); - n_bad++; + if ( indm & INDEXING_REFINE ) { + + UnitCell *out; + + if ( refine_prediction(image, cr) ) { + crystal_set_user_flag(cr, 1); + n_bad++; + } + + if ( (indm & INDEXING_CHECK_CELL_COMBINATIONS) + || (indm & INDEXING_CHECK_CELL_AXES) ) + { + + /* Check that the cell parameters are still + * within the tolerance */ + out = match_cell(crystal_get_cell(cr), + ipriv->target_cell, 0, + ipriv->tolerance, 0); + + if ( out == NULL ) { + crystal_set_user_flag(cr, 1); + n_bad++; + } + + cell_free(out); + + } + } } -- cgit v1.2.3 From 772ed44b99c2429ac9bcaab327583e92464e2fd5 Mon Sep 17 00:00:00 2001 From: Valerio Mariani Date: Fri, 10 Mar 2017 14:29:55 +0100 Subject: Peakfinder8 in CrystFEL. Same results as Anton's Cheetah implementation --- libcrystfel/Makefile.am | 6 +- libcrystfel/src/peakfinder8.c | 1051 +++++++++++++++++++++++++++++++++++++++++ libcrystfel/src/peakfinder8.h | 42 ++ libcrystfel/src/peaks.c | 28 +- libcrystfel/src/peaks.h | 14 +- 5 files changed, 1132 insertions(+), 9 deletions(-) create mode 100644 libcrystfel/src/peakfinder8.c create mode 100644 libcrystfel/src/peakfinder8.h (limited to 'libcrystfel') diff --git a/libcrystfel/Makefile.am b/libcrystfel/Makefile.am index d7142f6a..8af82552 100644 --- a/libcrystfel/Makefile.am +++ b/libcrystfel/Makefile.am @@ -10,7 +10,8 @@ libcrystfel_la_SOURCES = src/reflist.c src/utils.c src/cell.c src/detector.c \ src/render.c src/index.c src/dirax.c src/mosflm.c \ src/cell-utils.c src/integer_matrix.c src/crystal.c \ src/xds.c src/integration.c src/predict-refine.c \ - src/histogram.c src/events.c src/felix.c + src/histogram.c src/events.c src/felix.c \ + src/peakfinder8.c if HAVE_FFTW libcrystfel_la_SOURCES += src/asdf.c @@ -30,7 +31,8 @@ libcrystfel_la_include_HEADERS = ${top_srcdir}/version.h \ src/integer_matrix.h src/crystal.h \ src/xds.h src/predict-refine.h \ src/integration.h src/histogram.h \ - src/events.h src/asdf.h src/felix.h + src/events.h src/asdf.h src/felix.h \ + src/peakfinder8.h AM_CPPFLAGS = -DDATADIR=\""$(datadir)"\" -I$(top_builddir)/lib -Wall AM_CPPFLAGS += -I$(top_srcdir)/lib @LIBCRYSTFEL_CFLAGS@ diff --git a/libcrystfel/src/peakfinder8.c b/libcrystfel/src/peakfinder8.c new file mode 100644 index 00000000..2af8608a --- /dev/null +++ b/libcrystfel/src/peakfinder8.c @@ -0,0 +1,1051 @@ +/* + * peakfinder8.c + * + * The processing pipeline for one image + * + * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. + * + * Authors: + * 2017 Valerio Mariani + * 2017 Anton Barty + * 2017 Oleksandr Yefanov + * + * This file is part of CrystFEL. + * + * CrystFEL is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * CrystFEL is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with CrystFEL. If not, see . + * + */ + + +#include "peakfinder8.h" + + +struct radius_maps +{ + float **r_maps; + int n_rmaps; +}; + + +struct peakfinder_mask +{ + char **masks; + int n_masks; +}; + + +struct radial_stats +{ + float *roffset; + float *rthreshold; + int n_rad_bins; +}; + + +struct peakfinder_panel_data +{ + float ** panel_data; + int *panel_h; + int *panel_w; + int num_panels; +}; + + +struct peakfinder_peak_data +{ + int num_found_peaks; + int *npix; + float *com_fs; + float *com_ss; + float *tot_i; + float *max_i; + float *sigma; + float *snr; +}; + + +static double compute_r_sigma(float *rsigma, int *rcount, float* roffset, + int i) +{ + return sqrt(rsigma[i] / rcount[i] - + ((roffset[i] / rcount[i]) * + (roffset[i] / rcount[i]))); +} + + +static void update_radial_stats (float *roffset, float *rsigma, int *rcount, + float data, int curr_radius) +{ + roffset[curr_radius] += data; + rsigma[curr_radius] += (data * data); + rcount[curr_radius] += 1; +} + + +static float get_radius(struct panel p, float fs, float ss) +{ + float x, y; + + /* Calculate 3D position of given position, in m */ + x = (p.cnx + fs * p.fsx + ss * p.ssx); + y = (p.cny + fs * p.fsy + ss * p.ssy); + + return sqrt(x * x + y * y); +} + + +static struct radius_maps* compute_radius_maps(struct detector *det) +{ + int i, u, iss, ifs; + struct panel p; + struct radius_maps *rm = NULL; + + rm = malloc(sizeof(struct radius_maps)); + if ( rm == NULL ) { + return NULL; + } + + rm->r_maps = malloc(det->n_panels*sizeof(float*)); + if ( rm->r_maps == NULL ) { + free(rm); + return NULL; + } + + rm->n_rmaps = det->n_panels; + + for( i=0 ; in_panels ; i++ ) { + + p = det->panels[i]; + rm->r_maps[i] = malloc(p.h*p.w*sizeof(float)); + + if ( rm->r_maps[i] == NULL ) { + for ( u = 0; ur_maps[u]); + } + free(rm); + return NULL; + } + + for ( iss=0 ; issr_maps[i][rmi] = get_radius(p, ifs, iss); + } + } + + } + + return rm; +} + + +static void free_radius_maps(struct radius_maps *r_maps) +{ + int i; + + for ( i=0 ; in_rmaps ; i++ ) { + free(r_maps->r_maps[i]); + } + free(r_maps->r_maps); + free(r_maps); +} + + +static struct peakfinder_mask *create_peakfinder_mask(struct image *img, + struct radius_maps *rmps, + int min_res, + int max_res) +{ + int i; + struct peakfinder_mask *msk; + + msk = malloc(sizeof(struct peakfinder_mask)); + msk->masks = malloc(img->det->n_panels*sizeof(char*)); + msk->n_masks = img->det->n_panels; + for ( i=0; idet->n_panels; i++) { + + struct panel p; + int iss, ifs; + + p = img->det->panels[i]; + + msk->masks[i] = calloc(p.w*p.h,sizeof(char)); + + for ( iss=0 ; issr_maps[i][idx] < max_res + && rmps->r_maps[i][idx] > min_res ) { + + if (! ( ( img->bad != NULL ) + && ( img->bad[i] != NULL ) + && ( img->bad[i][idx] != 0 ) ) ) { + msk->masks[i][idx] = 1; + } + + + } + + } + } + } + return msk; +} + + +static void free_peakfinder_mask(struct peakfinder_mask * pfmask) +{ + int i; + + for ( i=0 ; in_masks ; i++ ) { + free(pfmask->masks[i]); + } + free(pfmask->masks); + free(pfmask); +} + + +static int compute_num_radial_bins( int num_panels, int *w, int *h, + float **r_maps ) +{ + int m; + + float max_r; + int max_r_int; + + max_r = -1e9; + for ( m=0 ; m max_r ) { + max_r = r_maps[m][i]; + } + } + } + + max_r_int = (int)ceil(max_r) + 1; + + return max_r_int; +} + + +static struct peakfinder_panel_data* allocate_panel_data(struct image *img) +{ + + struct peakfinder_panel_data *pfdata; + + + pfdata = malloc(sizeof(struct peakfinder_panel_data)); + if ( pfdata == NULL ) { + return NULL; + } + + pfdata->panel_h = malloc(img->det->n_panels*sizeof(int)); + if ( pfdata->panel_h == NULL ) { + free(pfdata); + return NULL; + } + + pfdata->panel_w = malloc(img->det->n_panels*sizeof(int)); + if ( pfdata->panel_w == NULL ) { + free(pfdata); + free(pfdata->panel_w); + return NULL; + } + + pfdata->panel_data = malloc(img->det->n_panels*sizeof(float*)); + if ( pfdata->panel_data == NULL ) { + free(pfdata); + free(pfdata->panel_w); + free(pfdata->panel_h); + return NULL; + } + + pfdata->num_panels = img->det->n_panels; + + return pfdata; +} + + +static void free_panel_data(struct peakfinder_panel_data *pfdata) +{ + free(pfdata->panel_data); + free(pfdata->panel_w); + free(pfdata->panel_h); + free(pfdata); +} + + +static struct radial_stats* allocate_radial_stats(int num_rad_bins) +{ + struct radial_stats* rstats; + + rstats = malloc(sizeof(struct radial_stats)); + + rstats->roffset = malloc(num_rad_bins*sizeof(float)); + if ( rstats->roffset == NULL ) { + return NULL; + } + + rstats->rthreshold = malloc(num_rad_bins*sizeof(float)); + if ( rstats->rthreshold == NULL ) { + free(rstats->roffset); + return NULL; + } + + rstats->n_rad_bins = num_rad_bins; + + return rstats; +} + + +static void free_radial_stats(struct radial_stats *rstats) +{ + free(rstats->roffset); + free(rstats->rthreshold); + free(rstats); +} + + +static int compute_radial_stats(float ** panel_data, + int num_panels, + int *w, + int *h, + float **r_maps, + char **masks, + float *rthreshold, + float *roffset, + int num_rad_bins, + float min_snr, + float acd_threshold, + int iterations) +{ + + int i, ri, pi; + int it_counter; + int curr_r; + float data; + float this_offset, this_sigma; + float *rsigma; + int *rcount; + + + for ( i=0 ; inpix = malloc(max_num_peaks*sizeof(int)); + if ( pkdata->npix == NULL ) { + free(pkdata); + free(pkdata->npix); + return NULL; + } + + pkdata->com_fs = malloc(max_num_peaks*sizeof(float)); + if ( pkdata->com_fs == NULL ) { + free(pkdata->npix); + free(pkdata); + return NULL; + } + + pkdata->com_ss = malloc(max_num_peaks*sizeof(float)); + if ( pkdata->com_ss == NULL ) { + free(pkdata->npix); + free(pkdata->com_fs); + free(pkdata); + return NULL; + } + + pkdata->tot_i = malloc(max_num_peaks*sizeof(float)); + if ( pkdata->tot_i == NULL ) { + free(pkdata->npix); + free(pkdata->com_fs); + free(pkdata->com_ss); + free(pkdata); + return NULL; + } + + pkdata->max_i = malloc(max_num_peaks*sizeof(float)); + if ( pkdata->max_i == NULL ) { + free(pkdata->npix); + free(pkdata->com_fs); + free(pkdata->com_ss); + free(pkdata->tot_i); + free(pkdata); + return NULL; + } + + pkdata->sigma = malloc(max_num_peaks*sizeof(float)); + if ( pkdata->sigma == NULL ) { + free(pkdata->npix); + free(pkdata->com_fs); + free(pkdata->com_ss); + free(pkdata->tot_i); + free(pkdata->max_i); + free(pkdata); + return NULL; + } + + pkdata->snr = malloc(max_num_peaks*sizeof(float)); + if ( pkdata->snr == NULL ) { + free(pkdata->npix); + free(pkdata->com_fs); + free(pkdata->com_ss); + free(pkdata->tot_i); + free(pkdata->max_i); + free(pkdata->sigma); + free(pkdata); + return NULL; + } + + return pkdata; +} + + +static void free_peak_data(struct peakfinder_peak_data *pkdata) { + free(pkdata->npix); + free(pkdata->com_fs); + free(pkdata->com_ss); + free(pkdata->tot_i); + free(pkdata->max_i); + free(pkdata->sigma); + free(pkdata->snr); + free(pkdata); +} + + +static void peak_search(int p, int *inss, int *infs, int h, int w, + int multip_w, float *copy, char *mask, float *r_map, + float *rthreshold, float *roffset, + char *pix_in_peak_map, int *peak_pixels, + int *num_pix_in_peak, float *sum_com_fs, + float *sum_com_ss, float *sum_i, int max_pix_count) +{ + + int k, pi; + int curr_radius; + float curr_threshold; + int curr_fs; + int curr_ss; + float curr_i; + + int search_fs[9] = { 0, -1, 0, 1, -1, 1, -1, 0, 1 }; + int search_ss[9] = { 0, -1, -1, -1, 0, 0, 1, 1, 1 }; + int search_n = 9; + + for ( k=0; k= w ) + continue; + if ( (inss[p] + search_ss[k]) < 0 ) + continue; + if ( (inss[p] + search_ss[k]) >= h ) + continue; + + curr_fs = infs[p] + search_fs[k]; + curr_ss = inss[p] + search_ss[k]; + pi = curr_fs + curr_ss * multip_w; + + curr_radius = (int)rint(r_map[pi]); + curr_threshold = rthreshold[curr_radius]; + + if ( copy[pi] > curr_threshold + && pix_in_peak_map[pi] == 0 + && mask[pi] != 0 ) { + + curr_i = copy[pi] - roffset[curr_radius]; + + *sum_i += curr_i; + *sum_com_fs += curr_i * ((float)curr_fs); + *sum_com_ss += curr_i * ((float)curr_ss); + + inss[*num_pix_in_peak] = inss[p] + search_ss[k]; + infs[*num_pix_in_peak] = infs[p] + search_fs[k]; + + pix_in_peak_map[pi] = 1; + if ( *num_pix_in_peak < max_pix_count ) + peak_pixels[*num_pix_in_peak] = pi; + *num_pix_in_peak = *num_pix_in_peak+1; + } + } +} + + +static void search_in_ring(int ring_width, int com_fs_int, int com_ss_int, + int h, int w, int multip_w, float *copy, + float *r_map, float *rthreshold, + float *roffset, char *pix_in_peak_map, + char *mask, float *local_sigma, + float *local_offset, float *background_max_i, + int com_idx, + int local_bg_radius) +{ + int ssj, fsi; + float pix_radius; + int pi; + int curr_radius; + float curr_threshold; + float curr_i; + + int np_sigma; + int np_counted; + int local_radius; + + float sum_i; + float sum_i_squared; + + ring_width = 2 * local_bg_radius; + + sum_i = 0; + sum_i_squared = 0; + np_sigma = 0; + np_counted = 0; + local_radius = 0; + + for ( ssj = -ring_width ; ssj= w ) + continue; + if ( (com_ss_int + ssj) < 0 ) + continue; + if ( (com_ss_int + ssj) >= h ) + continue; + + pix_radius = sqrt(fsi * fsi + ssj * ssj); + if ( pix_radius>ring_width ) + continue; + + pi = (com_fs_int +fsi) + (com_ss_int +ssj) * multip_w; + + curr_radius = rint(r_map[pi]); + curr_threshold = rthreshold[curr_radius]; + curr_i = copy[pi]; + + if ( copy[pi] < curr_threshold + && pix_in_peak_map[pi] == 0 + && mask[pi] != 0 ) { + + np_sigma++; + sum_i += curr_i; + sum_i_squared += (curr_i * curr_i); + + if ( curr_i > *background_max_i ) { + *background_max_i = curr_i; + } + } + np_counted += 1; + } + } + + if ( np_sigma != 0 ) { + *local_offset = sum_i / np_sigma; + *local_sigma = sqrt(sum_i_squared / np_sigma - + ((sum_i / np_sigma) * (sum_i / np_sigma))); + } else { + local_radius = rint(r_map[lrint(com_idx)]); + *local_offset = roffset[local_radius]; + *local_sigma = 0.01; + } +} + + +static int peakfinder8_multipanel(float *roffset, float *rthreshold, + float *data, char *mask, float *r_map, int w, + int h, int max_n_peaks, int *num_found_peaks, + int *npix, float *com_fs, float *com_ss, + float *tot_i, float *max_i, float *sigma, + float *snr, int min_pix_count, + int max_pix_count, int local_bg_radius, + float min_snr, int multip_w) +{ + int i, p, iss, ifs; + int peak_count; + int peak_idx; + int num_pix_in_peak, lt_num_pix_in_pk; + int idx, com_idx; + int curr_radius; + int ring_width; + int curr_idx; + int curr_fs, curr_ss; + int *peak_pixels; + int *infs; + int *inss; + float peak_com_fs, peak_com_ss; + float sum_com_fs, sum_com_ss, sum_i; + float peak_com_fs_int, peak_com_ss_int; + float curr_threshold; + float peak_tot_i, pk_tot_i_raw; + float peak_max_i, pk_max_i_raw; + float local_sigma; + float local_offset; + float background_max_i; + float f_background_thresh; + float peak_snr; + float curr_i, curr_i_raw; + float *copy; + char *pix_in_peak_map; + + copy = (float*) malloc(w*h*sizeof(float)); + if ( copy == NULL ) { + return 1; + } + + memcpy(copy, data, w*h*sizeof(float)); + pix_in_peak_map = calloc(w*h, sizeof(char)); + if ( pix_in_peak_map == NULL ) { + free(copy); + return 1; + } + + infs = (int *) calloc(w*h, sizeof(int)); + if ( infs == NULL ) { + free(pix_in_peak_map); + free(copy); + return 1; + } + + inss = (int *) calloc(w*h, sizeof(int)); + if ( inss == NULL ) { + free(infs); + free(pix_in_peak_map); + free(copy); + return 1; + } + + // TODO Possible bug: should 0,0 be included? + peak_pixels = calloc(max_pix_count, sizeof(int)); + if ( peak_pixels == NULL ) { + free(inss); + free(infs); + free(pix_in_peak_map); + free(copy); + return 1; + } + + for ( i = 0; i < w*h; i++ ) { + copy[i] *= mask[i]; + } + + peak_count = 0; + num_pix_in_peak = 0; + + for ( iss=0 ; iss curr_threshold && + pix_in_peak_map[idx] == 0 ) { + + infs[0] = ifs; + inss[0] = iss; + + peak_pixels[0] = idx; + num_pix_in_peak = 1; + sum_i = 0; + sum_com_fs = 0; + sum_com_ss = 0; + + do { + lt_num_pix_in_pk = num_pix_in_peak; + + for ( p=0; p max_pix_count) { + continue; + } + + peak_com_fs = sum_com_fs / fabs(sum_i); + peak_com_ss = sum_com_ss / fabs(sum_i); + com_idx = rint(peak_com_fs) + + rint(peak_com_ss) * multip_w; + + peak_com_fs_int = (int)rint(peak_com_fs); + peak_com_ss_int = (int)rint(peak_com_ss); + + local_sigma = 0.0; + local_offset = 0.0; + background_max_i = 0.0; + + search_in_ring(ring_width, peak_com_fs_int, + peak_com_ss_int, h, w, multip_w, + copy, r_map, rthreshold, + roffset, pix_in_peak_map, + mask, &local_sigma, + &local_offset, + &background_max_i, + com_idx, local_bg_radius); + + peak_tot_i = 0; + pk_tot_i_raw = 0; + peak_max_i = 0; + pk_max_i_raw = 0; + sum_com_fs = 0; + sum_com_ss = 0; + + for ( peak_idx = 1 ; + peak_idx < num_pix_in_peak && + peak_idx <= max_pix_count ; + peak_idx++ ) { + + curr_idx = peak_pixels[peak_idx]; + + curr_i_raw = copy[curr_idx]; + curr_i = curr_i_raw - local_offset; + + peak_tot_i += curr_i; + pk_tot_i_raw += curr_i_raw; + + curr_fs = curr_idx % multip_w; + curr_ss = curr_idx / multip_w; + sum_com_fs += curr_i * ((float)curr_fs); + sum_com_ss += curr_i * ((float)curr_ss); + + if ( curr_i_raw > pk_max_i_raw ) + pk_max_i_raw = curr_i_raw; + if ( curr_i > peak_max_i ) + peak_max_i = curr_i; + } + + peak_com_fs = sum_com_fs / fabs(peak_tot_i); + peak_com_ss = sum_com_ss / fabs(peak_tot_i); + peak_snr = peak_tot_i / local_sigma; + + + if (peak_snr < min_snr) { + continue; + } + + f_background_thresh = 1; + f_background_thresh *= + background_max_i - local_offset; + if (peak_max_i < f_background_thresh) + continue; + + if ( num_pix_in_peak >= min_pix_count + && num_pix_in_peak <= max_pix_count ) { + + if ( peak_tot_i == 0 ) + continue; + + if ( peak_count < max_n_peaks ) { + + int pidx; + pidx = peak_count; + + npix[pidx] = num_pix_in_peak; + com_fs[pidx] = peak_com_fs; + com_ss[pidx] = peak_com_ss; + tot_i[pidx] = peak_tot_i; + max_i[pidx] = peak_max_i; + sigma[pidx] = local_sigma; + snr[pidx] = peak_snr; + peak_count++; + } + else { + peak_count++; + } + } + } + } + } + *num_found_peaks = peak_count; + + free(peak_pixels); + free(inss); + free(infs); + free(pix_in_peak_map); + free(copy); + + return 0; +} + +static int peakfinder8_singlepanel(float *roffset, float *rthreshold, + float *data, char *mask, float *r_map, int w, + int h, int max_n_peaks, int *num_found_peaks, + int *npix, float *com_fs, float *com_ss, + float *tot_i, float *max_i, float *sigma, + float *snr, int min_pix_count, + int max_pix_count, int local_bg_radius, + float min_snr) +{ + return peakfinder8_multipanel(roffset, rthreshold, data, mask, r_map, + w, h, max_n_peaks, num_found_peaks, npix, + com_fs, com_ss, tot_i, max_i, sigma, snr, + min_pix_count, max_pix_count, + local_bg_radius, min_snr, w); + +} + + +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) +{ + struct radius_maps *rmaps; + struct peakfinder_mask *pfmask; + struct peakfinder_panel_data *pfdata; + struct radial_stats *rstats; + struct peakfinder_peak_data *pkdata; + int num_rad_bins; + int pi; + int ret; + int num_found_peaks; + + if ( img-> det == NULL) { + return 1; + } + + rmaps = compute_radius_maps(img->det); + if ( rmaps == NULL ) { + return 1; + } + + pfmask = create_peakfinder_mask(img, rmaps, min_res, max_res); + if ( rmaps == NULL ) { + free_radius_maps(rmaps); + return 1; + } + + pfdata = allocate_panel_data(img); + if ( pfdata == NULL) { + free_radius_maps(rmaps); + free_peakfinder_mask(pfmask); + } + + for ( pi=0 ; pidet->n_panels ; pi++ ) { + pfdata->panel_h[pi] = img->det->panels[pi].h; + pfdata->panel_w[pi] = img->det->panels[pi].w; + pfdata->panel_data[pi] = img->dp[pi]; + } + + num_rad_bins = compute_num_radial_bins(img->det->n_panels, + pfdata->panel_w, + pfdata->panel_h, + rmaps->r_maps); + + rstats = allocate_radial_stats(num_rad_bins); + if ( rstats == NULL ) { + free_radius_maps(rmaps); + free_peakfinder_mask(pfmask); + free_panel_data(pfdata); + return 1; + } + + ret = compute_radial_stats(pfdata->panel_data, pfdata->num_panels, + pfdata->panel_w, pfdata->panel_h, + rmaps->r_maps, + pfmask->masks, rstats->rthreshold, + rstats->roffset, num_rad_bins, + min_snr, threshold, 5); + if ( ret != 0 ) { + free_radius_maps(rmaps); + free_peakfinder_mask(pfmask); + free_panel_data(pfdata); + free_radial_stats(rstats); + return 1; + } + + pkdata = allocate_peak_data(max_n_peaks); + if ( pkdata == NULL ) { + free_radius_maps(rmaps); + free_peakfinder_mask(pfmask); + free_panel_data(pfdata); + free_radial_stats(rstats); + return 1; + } + + for ( pi=0 ; pidet->n_panels ; pi++) { + + int peaks_to_add; + int pki; + + num_found_peaks = 0; + + ret = peakfinder8_singlepanel(rstats->roffset, + rstats->rthreshold, + pfdata->panel_data[pi], + pfmask->masks[pi], + rmaps->r_maps[pi], + pfdata->panel_w[pi], + pfdata->panel_h[pi], + max_n_peaks, + &num_found_peaks, + pkdata->npix, + pkdata->com_fs, + pkdata->com_ss, + pkdata->tot_i, + pkdata->max_i, + pkdata->sigma, + pkdata->snr, + min_pix_count, + max_pix_count, + local_bg_radius, + min_snr); + + peaks_to_add = num_found_peaks; + + if ( num_found_peaks > max_n_peaks ) { + peaks_to_add = max_n_peaks; + } + + for ( pki=0 ; pkidet->panels[pi]; + + image_add_feature(img->features, + pkdata->com_fs[pki], + pkdata->com_ss[pki], + p, + img, + pkdata->tot_i[pki], + NULL); + + } + } + + free_radius_maps(rmaps); + free_peakfinder_mask(pfmask); + free_panel_data(pfdata); + free_radial_stats(rstats); + free_peak_data(pkdata); + return 0; +} diff --git a/libcrystfel/src/peakfinder8.h b/libcrystfel/src/peakfinder8.h new file mode 100644 index 00000000..34cf5d62 --- /dev/null +++ b/libcrystfel/src/peakfinder8.h @@ -0,0 +1,42 @@ +/* + * peakfinder8.h + * + * The processing pipeline for one image + * + * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. + * + * Authors: + * 2017 Valerio Mariani + * 2017 Anton Barty + * 2017 Oleksandr Yefanov + * + * This file is part of CrystFEL. + * + * CrystFEL is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * CrystFEL is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with CrystFEL. If not, see . + * + */ + +#ifndef PEAKFINDER8_H +#define PEAKFINDER8_H + +#include "image.h" + +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); + +#endif // PEAKFINDER8_H diff --git a/libcrystfel/src/peaks.c b/libcrystfel/src/peaks.c index cb2a9f5a..ea2fe3e3 100644 --- a/libcrystfel/src/peaks.c +++ b/libcrystfel/src/peaks.c @@ -3,7 +3,7 @@ * * Peak search and other image analysis * - * Copyright © 2012-2016 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * Copyright © 2012 Richard Kirian * @@ -12,6 +12,7 @@ * 2012 Kenneth Beyerlein * 2011 Andrew Martin * 2011 Richard Kirian + * 2017 Valerio Mariani * * This file is part of CrystFEL. * @@ -52,6 +53,7 @@ #include "reflist-utils.h" #include "cell-utils.h" #include "geometry.h" +#include "peakfinder8.h" static int cull_peaks_in_panel(struct image *image, struct panel *p) @@ -517,8 +519,8 @@ static void search_peaks_in_panel(struct image *image, float threshold, void search_peaks(struct image *image, float threshold, float min_gradient, - float min_snr, double ir_inn, double ir_mid, double ir_out, - int use_saturated) + float min_snr, double ir_inn, double ir_mid, + double ir_out, int use_saturated) { int i; @@ -541,6 +543,26 @@ void search_peaks(struct image *image, float threshold, float min_gradient, } +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) +{ + if ( image->features != NULL ) { + image_feature_list_free(image->features); + } + image->features = image_feature_list_new(); + image->num_peaks = 0; + image->num_saturated_peaks = 0; + + return peakfinder8(image, max_n_peaks, threshold, min_snr, + min_pix_count, max_pix_count, + local_bg_radius, min_res, + max_res); +} + + int peak_sanity_check(struct image *image, Crystal **crystals, int n_cryst) { int n_feat = 0; diff --git a/libcrystfel/src/peaks.h b/libcrystfel/src/peaks.h index ba724419..dfc286cf 100644 --- a/libcrystfel/src/peaks.h +++ b/libcrystfel/src/peaks.h @@ -3,11 +3,12 @@ * * Peak search and other image analysis * - * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * * Authors: * 2010-2015 Thomas White + * 2017 Valerio Mariani * * This file is part of CrystFEL. * @@ -45,9 +46,14 @@ extern "C" { extern int *make_BgMask(struct image *image, struct panel *p, double ir_inn); extern void search_peaks(struct image *image, float threshold, - float min_gradient, float min_snr, - double ir_inn, double ir_mid, double ir_out, - int use_saturated); + float min_gradient, float min_snr, double ir_inn, + double ir_mid, double ir_out, int use_saturated); + +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); extern int peak_sanity_check(struct image *image, Crystal **crystals, int n_cryst); -- cgit v1.2.3 From c9e4b3cd8c1a9421787168ce5efb20c89e5286fd Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 10 Mar 2017 15:15:03 +0100 Subject: cell_print(): Show reciprocal angles --- libcrystfel/src/cell-utils.c | 6 ++++++ 1 file changed, 6 insertions(+) (limited to 'libcrystfel') diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c index 69d4174a..84a14b5c 100644 --- a/libcrystfel/src/cell-utils.c +++ b/libcrystfel/src/cell-utils.c @@ -273,6 +273,12 @@ void cell_print(UnitCell *cell) STATUS("c* = %10.3e %10.3e %10.3e m^-1 (modulus %10.3e m^-1)\n", csx, csy, csz, modulus(csx, csy, csz)); + STATUS("alpha* = %6.2f deg, beta* = %6.2f deg, " + "gamma* = %6.2f deg\n", + rad2deg(angle_between(bsx, bsy, bsz, csx, csy, csz)), + rad2deg(angle_between(asx, asy, asz, csx, csy, csz)), + rad2deg(angle_between(asx, asy, asz, bsx, bsy, bsz))); + STATUS("Cell representation is %s.\n", cell_rep(cell)); } else { -- cgit v1.2.3 From de2237fc3f193e4f4916e3bc9268e52378fad939 Mon Sep 17 00:00:00 2001 From: Valerio Mariani Date: Fri, 10 Mar 2017 20:04:25 +0100 Subject: Fixed a couple of bugs reported by Tom --- libcrystfel/src/peakfinder8.c | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/peakfinder8.c b/libcrystfel/src/peakfinder8.c index 2af8608a..c312e120 100644 --- a/libcrystfel/src/peakfinder8.c +++ b/libcrystfel/src/peakfinder8.c @@ -370,7 +370,7 @@ static int compute_radial_stats(float ** panel_data, for ( it_counter=0 ; it_counter curr_threshold && @@ -762,7 +763,8 @@ static int peakfinder8_multipanel(float *roffset, float *rthreshold, sum_com_fs = 0; sum_com_ss = 0; - do { + do { + lt_num_pix_in_pk = num_pix_in_peak; for ( p=0; p Date: Sat, 11 Mar 2017 01:56:56 +0100 Subject: Fixed some more bugs reported by Yaroslav --- libcrystfel/src/peakfinder8.c | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/peakfinder8.c b/libcrystfel/src/peakfinder8.c index c312e120..fecc61da 100644 --- a/libcrystfel/src/peakfinder8.c +++ b/libcrystfel/src/peakfinder8.c @@ -930,6 +930,9 @@ int peakfinder8(struct image *img, int max_n_peaks, int pi; int ret; int num_found_peaks; + int iterations; + + iterations = 5; if ( img-> det == NULL) { return 1; @@ -941,7 +944,7 @@ int peakfinder8(struct image *img, int max_n_peaks, } pfmask = create_peakfinder_mask(img, rmaps, min_res, max_res); - if ( rmaps == NULL ) { + if ( pfmask == NULL ) { free_radius_maps(rmaps); return 1; } @@ -950,6 +953,7 @@ int peakfinder8(struct image *img, int max_n_peaks, if ( pfdata == NULL) { free_radius_maps(rmaps); free_peakfinder_mask(pfmask); + return 1; } for ( pi=0 ; pidet->n_panels ; pi++ ) { @@ -976,7 +980,7 @@ int peakfinder8(struct image *img, int max_n_peaks, rmaps->r_maps, pfmask->masks, rstats->rthreshold, rstats->roffset, num_rad_bins, - min_snr, threshold, 5); + min_snr, threshold, iterations); if ( ret != 0 ) { free_radius_maps(rmaps); free_peakfinder_mask(pfmask); @@ -1001,6 +1005,10 @@ int peakfinder8(struct image *img, int max_n_peaks, num_found_peaks = 0; + if ( img->det->panels[pi].no_index ) { + continue; + } + ret = peakfinder8_singlepanel(rstats->roffset, rstats->rthreshold, pfdata->panel_data[pi], -- cgit v1.2.3 From 8986433c649932c651b11719a31ba842f283fa2e Mon Sep 17 00:00:00 2001 From: Thomas White Date: Mon, 13 Mar 2017 11:06:55 +0100 Subject: Ask for C99 in configure.ac, remove weird C99 declarations --- libcrystfel/src/asdf.c | 2 -- libcrystfel/src/detector.c | 2 -- libcrystfel/src/events.c | 3 --- libcrystfel/src/image.c | 1 - libcrystfel/src/reflist-utils.c | 2 -- libcrystfel/src/statistics.c | 1 - 6 files changed, 11 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/asdf.c b/libcrystfel/src/asdf.c index 49917ad9..033688ba 100644 --- a/libcrystfel/src/asdf.c +++ b/libcrystfel/src/asdf.c @@ -28,8 +28,6 @@ * */ -#define _ISOC99_SOURCE - #ifdef HAVE_CONFIG_H #include #endif diff --git a/libcrystfel/src/detector.c b/libcrystfel/src/detector.c index 8c98f79e..90dbe535 100644 --- a/libcrystfel/src/detector.c +++ b/libcrystfel/src/detector.c @@ -31,8 +31,6 @@ * */ -#define _ISOC99_SOURCE -#define _GNU_SOURCE #include #include #include diff --git a/libcrystfel/src/events.c b/libcrystfel/src/events.c index 0f170bb5..8e4eb861 100644 --- a/libcrystfel/src/events.c +++ b/libcrystfel/src/events.c @@ -27,9 +27,6 @@ * */ -#define _ISOC99_SOURCE -#define _GNU_SOURCE - #include "events.h" #include "utils.h" diff --git a/libcrystfel/src/image.c b/libcrystfel/src/image.c index 7bb4cede..09f4958d 100644 --- a/libcrystfel/src/image.c +++ b/libcrystfel/src/image.c @@ -27,7 +27,6 @@ * */ -#define _ISOC99_SOURCE #include #include #include diff --git a/libcrystfel/src/reflist-utils.c b/libcrystfel/src/reflist-utils.c index 4621f4f4..380a1b8a 100644 --- a/libcrystfel/src/reflist-utils.c +++ b/libcrystfel/src/reflist-utils.c @@ -27,8 +27,6 @@ * */ -#define _ISOC99_SOURCE -#define _GNU_SOURCE #include #include #include diff --git a/libcrystfel/src/statistics.c b/libcrystfel/src/statistics.c index 56273fdb..ccf35194 100644 --- a/libcrystfel/src/statistics.c +++ b/libcrystfel/src/statistics.c @@ -30,7 +30,6 @@ #include #endif -#define _ISOC99_SOURCE #include #include #include -- cgit v1.2.3 From b2de09452e8edf050a8679e726f5075abd37e961 Mon Sep 17 00:00:00 2001 From: Valerio Mariani Date: Mon, 13 Mar 2017 11:17:57 +0100 Subject: Added saturated peak management to peakfinder8 --- libcrystfel/src/peakfinder8.c | 10 +++++++++- libcrystfel/src/peakfinder8.h | 2 +- libcrystfel/src/peaks.c | 4 ++-- libcrystfel/src/peaks.h | 2 +- 4 files changed, 13 insertions(+), 5 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/peakfinder8.c b/libcrystfel/src/peakfinder8.c index fecc61da..82a7e316 100644 --- a/libcrystfel/src/peakfinder8.c +++ b/libcrystfel/src/peakfinder8.c @@ -919,7 +919,7 @@ 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 max_res, int use_saturated) { struct radius_maps *rmaps; struct peakfinder_mask *pfmask; @@ -1042,6 +1042,14 @@ int peakfinder8(struct image *img, int max_n_peaks, p = &img->det->panels[pi]; + img->num_peaks += 1; + if ( pkdata->max_i[pki] > p->max_adu ) { + img->num_saturated_peaks++; + if ( !use_saturated ) { + continue; + } + } + image_add_feature(img->features, pkdata->com_fs[pki], pkdata->com_ss[pki], diff --git a/libcrystfel/src/peakfinder8.h b/libcrystfel/src/peakfinder8.h index 34cf5d62..1103ef96 100644 --- a/libcrystfel/src/peakfinder8.h +++ b/libcrystfel/src/peakfinder8.h @@ -37,6 +37,6 @@ 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 max_res, int use_saturated); #endif // PEAKFINDER8_H diff --git a/libcrystfel/src/peaks.c b/libcrystfel/src/peaks.c index ea2fe3e3..28c79538 100644 --- a/libcrystfel/src/peaks.c +++ b/libcrystfel/src/peaks.c @@ -547,7 +547,7 @@ 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 max_res, int use_saturated) { if ( image->features != NULL ) { image_feature_list_free(image->features); @@ -559,7 +559,7 @@ 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); + max_res, use_saturated); } diff --git a/libcrystfel/src/peaks.h b/libcrystfel/src/peaks.h index dfc286cf..a5095127 100644 --- a/libcrystfel/src/peaks.h +++ b/libcrystfel/src/peaks.h @@ -53,7 +53,7 @@ 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 max_res, int use_saturated); extern int peak_sanity_check(struct image *image, Crystal **crystals, int n_cryst); -- cgit v1.2.3 From 984379c6c48a2325134b5f1edfff78934604e1ae Mon Sep 17 00:00:00 2001 From: Valerio Mariani Date: Mon, 13 Mar 2017 11:37:38 +0100 Subject: Added management of max_num_peaks per image instead of per panel --- libcrystfel/src/peakfinder8.c | 2 ++ 1 file changed, 2 insertions(+) (limited to 'libcrystfel') diff --git a/libcrystfel/src/peakfinder8.c b/libcrystfel/src/peakfinder8.c index 82a7e316..9ad0e4bd 100644 --- a/libcrystfel/src/peakfinder8.c +++ b/libcrystfel/src/peakfinder8.c @@ -1036,6 +1036,8 @@ int peakfinder8(struct image *img, int max_n_peaks, peaks_to_add = max_n_peaks; } + max_n_peaks -= peaks_to_add; + for ( pki=0 ; pki Date: Wed, 15 Mar 2017 14:42:58 +0100 Subject: largest_q(): Handle NULL detector --- libcrystfel/src/detector.c | 5 +++++ 1 file changed, 5 insertions(+) (limited to 'libcrystfel') diff --git a/libcrystfel/src/detector.c b/libcrystfel/src/detector.c index 90dbe535..52e4c330 100644 --- a/libcrystfel/src/detector.c +++ b/libcrystfel/src/detector.c @@ -2002,6 +2002,11 @@ double largest_q(struct image *image) struct rvec q; double tt; + if ( image->det == NULL ) { + ERROR("No detector geometry. assuming detector is infinite!\n"); + return INFINITY; + } + q = get_q_for_panel(image->det->furthest_out_panel, image->det->furthest_out_fs, image->det->furthest_out_ss, -- cgit v1.2.3 From 6aea5300ee1c70ac2e99a7f3e6de4d26fb7e243a Mon Sep 17 00:00:00 2001 From: Valerio Mariani Date: Wed, 22 Mar 2017 15:13:16 +0100 Subject: Completely revamped implementation of peakfinder8 --- libcrystfel/src/peakfinder8.c | 813 +++++++++++++++++++++++++----------------- 1 file changed, 484 insertions(+), 329 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/peakfinder8.c b/libcrystfel/src/peakfinder8.c index 9ad0e4bd..0e916a0c 100644 --- a/libcrystfel/src/peakfinder8.c +++ b/libcrystfel/src/peakfinder8.c @@ -1,7 +1,7 @@ /* * peakfinder8.c * - * The processing pipeline for one image + * The peakfinder8 algorithm * * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. @@ -50,6 +50,8 @@ struct radial_stats { float *roffset; float *rthreshold; + float *rsigma; + int *rcount; int n_rad_bins; }; @@ -63,12 +65,22 @@ struct peakfinder_panel_data }; +struct peakfinder_intern_data +{ + char *pix_in_peak_map; + int *infs; + int *inss; + int *peak_pixels; +}; + + struct peakfinder_peak_data { int num_found_peaks; int *npix; float *com_fs; float *com_ss; + int *com_index; float *tot_i; float *max_i; float *sigma; @@ -80,16 +92,16 @@ static double compute_r_sigma(float *rsigma, int *rcount, float* roffset, int i) { return sqrt(rsigma[i] / rcount[i] - - ((roffset[i] / rcount[i]) * - (roffset[i] / rcount[i]))); + ((roffset[i] / rcount[i]) * + (roffset[i] / rcount[i]))); } static void update_radial_stats (float *roffset, float *rsigma, int *rcount, - float data, int curr_radius) + float value, int curr_radius) { - roffset[curr_radius] += data; - rsigma[curr_radius] += (data * data); + roffset[curr_radius] += value; + rsigma[curr_radius] += (value * value); rcount[curr_radius] += 1; } @@ -112,12 +124,12 @@ static struct radius_maps* compute_radius_maps(struct detector *det) struct panel p; struct radius_maps *rm = NULL; - rm = malloc(sizeof(struct radius_maps)); + rm = (struct radius_maps *)malloc(sizeof(struct radius_maps)); if ( rm == NULL ) { return NULL; } - rm->r_maps = malloc(det->n_panels*sizeof(float*)); + rm->r_maps = (float **)malloc(det->n_panels*sizeof(float*)); if ( rm->r_maps == NULL ) { free(rm); return NULL; @@ -128,7 +140,7 @@ static struct radius_maps* compute_radius_maps(struct detector *det) for( i=0 ; in_panels ; i++ ) { p = det->panels[i]; - rm->r_maps[i] = malloc(p.h*p.w*sizeof(float)); + rm->r_maps[i] = (float *)malloc(p.h*p.w*sizeof(float)); if ( rm->r_maps[i] == NULL ) { for ( u = 0; umasks = malloc(img->det->n_panels*sizeof(char*)); + msk = (struct peakfinder_mask *)malloc(sizeof(struct peakfinder_mask)); + msk->masks =(char **) malloc(img->det->n_panels*sizeof(char*)); msk->n_masks = img->det->n_panels; for ( i=0; idet->n_panels; i++) { @@ -185,7 +197,7 @@ static struct peakfinder_mask *create_peakfinder_mask(struct image *img, p = img->det->panels[i]; - msk->masks[i] = calloc(p.w*p.h,sizeof(char)); + msk->masks[i] = (char *)calloc(p.w*p.h,sizeof(char)); for ( iss=0 ; issmasks[i][idx] = 1; } - } } @@ -225,22 +236,27 @@ static void free_peakfinder_mask(struct peakfinder_mask * pfmask) } -static int compute_num_radial_bins( int num_panels, int *w, int *h, - float **r_maps ) +static int compute_num_radial_bins(int num_panels, int *w, int *h, + float **r_maps ) { - int m; float max_r; int max_r_int; + int m; max_r = -1e9; + for ( m=0 ; m max_r ) { - max_r = r_maps[m][i]; + for ( iss=0 ; iss max_r ) { + max_r = r_maps[m][pidx]; + } } } } @@ -251,39 +267,40 @@ static int compute_num_radial_bins( int num_panels, int *w, int *h, } -static struct peakfinder_panel_data* allocate_panel_data(struct image *img) +static struct peakfinder_panel_data* allocate_panel_data(int num_panels) { struct peakfinder_panel_data *pfdata; - pfdata = malloc(sizeof(struct peakfinder_panel_data)); + pfdata = (struct peakfinder_panel_data *)malloc( + sizeof(struct peakfinder_panel_data)); if ( pfdata == NULL ) { return NULL; } - pfdata->panel_h = malloc(img->det->n_panels*sizeof(int)); + pfdata->panel_h = (int *)malloc(num_panels*sizeof(int)); if ( pfdata->panel_h == NULL ) { free(pfdata); return NULL; } - pfdata->panel_w = malloc(img->det->n_panels*sizeof(int)); + pfdata->panel_w = (int *)malloc(num_panels*sizeof(int)); if ( pfdata->panel_w == NULL ) { + free(pfdata->panel_h); free(pfdata); - free(pfdata->panel_w); return NULL; } - pfdata->panel_data = malloc(img->det->n_panels*sizeof(float*)); + pfdata->panel_data = (float **)malloc(num_panels*sizeof(float*)); if ( pfdata->panel_data == NULL ) { - free(pfdata); free(pfdata->panel_w); free(pfdata->panel_h); + free(pfdata); return NULL; } - pfdata->num_panels = img->det->n_panels; + pfdata->num_panels = num_panels; return pfdata; } @@ -302,16 +319,38 @@ static struct radial_stats* allocate_radial_stats(int num_rad_bins) { struct radial_stats* rstats; - rstats = malloc(sizeof(struct radial_stats)); + rstats = (struct radial_stats *)malloc(sizeof(struct radial_stats)); + if ( rstats == NULL ) { + return NULL; + } - rstats->roffset = malloc(num_rad_bins*sizeof(float)); + rstats->roffset = (float *)malloc(num_rad_bins*sizeof(float)); if ( rstats->roffset == NULL ) { + free(rstats); return NULL; } - rstats->rthreshold = malloc(num_rad_bins*sizeof(float)); + rstats->rthreshold = (float *)malloc(num_rad_bins*sizeof(float)); if ( rstats->rthreshold == NULL ) { + free(rstats); + free(rstats->roffset); + return NULL; + } + + rstats->rsigma = (float *)malloc(num_rad_bins*sizeof(float)); + if ( rstats->rsigma == NULL ) { + free(rstats); free(rstats->roffset); + free(rstats->rthreshold); + return NULL; + } + + rstats->rcount = (int *)malloc(num_rad_bins*sizeof(int)); + if ( rstats->rcount == NULL ) { + free(rstats); + free(rstats->roffset); + free(rstats->rthreshold); + free(rstats->rsigma); return NULL; } @@ -325,110 +364,92 @@ static void free_radial_stats(struct radial_stats *rstats) { free(rstats->roffset); free(rstats->rthreshold); + free(rstats->rsigma); + free(rstats->rcount); free(rstats); } -static int compute_radial_stats(float ** panel_data, - int num_panels, - int *w, - int *h, - float **r_maps, - char **masks, - float *rthreshold, - float *roffset, - int num_rad_bins, - float min_snr, - float acd_threshold, - int iterations) +static void fill_radial_bins(float *data, + int w, + int h, + float *r_map, + char *mask, + float *rthreshold, + float *roffset, + float *rsigma, + int *rcount) { + int iss, ifs; + int pidx; - int i, ri, pi; - int it_counter; int curr_r; - float data; - float this_offset, this_sigma; - float *rsigma; - int *rcount; + float value; + for ( iss = 0; issnpix = malloc(max_num_peaks*sizeof(int)); + pkdata->npix = (int *)malloc(max_num_peaks*sizeof(int)); if ( pkdata->npix == NULL ) { free(pkdata); free(pkdata->npix); return NULL; } - pkdata->com_fs = malloc(max_num_peaks*sizeof(float)); + pkdata->com_fs = (float *)malloc(max_num_peaks*sizeof(float)); if ( pkdata->com_fs == NULL ) { free(pkdata->npix); free(pkdata); return NULL; } - pkdata->com_ss = malloc(max_num_peaks*sizeof(float)); + pkdata->com_ss = (float *)malloc(max_num_peaks*sizeof(float)); if ( pkdata->com_ss == NULL ) { free(pkdata->npix); free(pkdata->com_fs); @@ -463,41 +485,55 @@ struct peakfinder_peak_data *allocate_peak_data(int max_num_peaks) return NULL; } - pkdata->tot_i = malloc(max_num_peaks*sizeof(float)); + pkdata->com_index = (int *)malloc(max_num_peaks*sizeof(int)); + if ( pkdata->com_ss == NULL ) { + free(pkdata->npix); + free(pkdata->com_fs); + free(pkdata->com_ss); + free(pkdata); + return NULL; + } + + + pkdata->tot_i = (float *)malloc(max_num_peaks*sizeof(float)); if ( pkdata->tot_i == NULL ) { free(pkdata->npix); free(pkdata->com_fs); free(pkdata->com_ss); + free(pkdata->com_index); free(pkdata); return NULL; } - pkdata->max_i = malloc(max_num_peaks*sizeof(float)); + pkdata->max_i = (float *)malloc(max_num_peaks*sizeof(float)); if ( pkdata->max_i == NULL ) { free(pkdata->npix); free(pkdata->com_fs); free(pkdata->com_ss); + free(pkdata->com_index); free(pkdata->tot_i); free(pkdata); return NULL; } - pkdata->sigma = malloc(max_num_peaks*sizeof(float)); + pkdata->sigma = (float *)malloc(max_num_peaks*sizeof(float)); if ( pkdata->sigma == NULL ) { free(pkdata->npix); free(pkdata->com_fs); free(pkdata->com_ss); + free(pkdata->com_index); free(pkdata->tot_i); free(pkdata->max_i); free(pkdata); return NULL; } - pkdata->snr = malloc(max_num_peaks*sizeof(float)); + pkdata->snr = (float *)malloc(max_num_peaks*sizeof(float)); if ( pkdata->snr == NULL ) { free(pkdata->npix); free(pkdata->com_fs); free(pkdata->com_ss); + free(pkdata->com_index); free(pkdata->tot_i); free(pkdata->max_i); free(pkdata->sigma); @@ -513,6 +549,7 @@ static void free_peak_data(struct peakfinder_peak_data *pkdata) { free(pkdata->npix); free(pkdata->com_fs); free(pkdata->com_ss); + free(pkdata->com_index); free(pkdata->tot_i); free(pkdata->max_i); free(pkdata->sigma); @@ -521,11 +558,70 @@ static void free_peak_data(struct peakfinder_peak_data *pkdata) { } -static void peak_search(int p, int *inss, int *infs, int h, int w, - int multip_w, float *copy, char *mask, float *r_map, +static struct peakfinder_intern_data *allocate_peakfinder_intern_data( + int data_size, int max_pix_count) +{ + + struct peakfinder_intern_data *intern_data; + + intern_data = (struct peakfinder_intern_data *)malloc( + sizeof(struct peakfinder_intern_data)); + if ( intern_data == NULL ) { + return NULL; + } + + intern_data->pix_in_peak_map =(char *)calloc(data_size, sizeof(char)); + if ( intern_data->pix_in_peak_map == NULL ) { + free(intern_data); + return NULL; + } + + intern_data->infs =(int *)calloc(data_size, sizeof(int)); + if ( intern_data->infs == NULL ) { + free(intern_data->pix_in_peak_map); + free(intern_data); + return NULL; + } + + intern_data->inss =(int *)calloc(data_size, sizeof(int)); + if ( intern_data->inss == NULL ) { + free(intern_data->pix_in_peak_map); + free(intern_data->infs); + free(intern_data); + return NULL; + } + + intern_data->peak_pixels =(int *)calloc(max_pix_count, sizeof(int)); + if ( intern_data->peak_pixels == NULL ) { + free(intern_data->pix_in_peak_map); + free(intern_data->infs); + free(intern_data->inss); + free(intern_data); + return NULL; + } + + return intern_data; +} + + +static void free_peakfinder_intern_data(struct peakfinder_intern_data *pfid) +{ + free(pfid->peak_pixels); + free(pfid->pix_in_peak_map); + free(pfid->infs); + free(pfid->inss); + free(pfid); +} + + + +static void peak_search(int p, + struct peakfinder_intern_data *pfinter, + float *copy, char *mask, float *r_map, float *rthreshold, float *roffset, - char *pix_in_peak_map, int *peak_pixels, - int *num_pix_in_peak, float *sum_com_fs, + int *num_pix_in_peak, int asic_size_fs, + int asic_size_ss, int aifs, int aiss, + int num_pix_fs, float *sum_com_fs, float *sum_com_ss, float *sum_i, int max_pix_count) { @@ -542,24 +638,24 @@ static void peak_search(int p, int *inss, int *infs, int h, int w, for ( k=0; kinfs[p] + search_fs[k]) < 0 ) continue; - if ( (infs[p] + search_fs[k]) >= w ) + if ( (pfinter->infs[p] + search_fs[k]) >= asic_size_fs ) continue; - if ( (inss[p] + search_ss[k]) < 0 ) + if ( (pfinter->inss[p] + search_ss[k]) < 0 ) continue; - if ( (inss[p] + search_ss[k]) >= h ) + if ( (pfinter->inss[p] + search_ss[k]) >= asic_size_ss ) continue; - curr_fs = infs[p] + search_fs[k]; - curr_ss = inss[p] + search_ss[k]; - pi = curr_fs + curr_ss * multip_w; + curr_fs = pfinter->infs[p] + search_fs[k] + aifs * asic_size_fs; + curr_ss = pfinter->inss[p] + search_ss[k] + aiss * asic_size_ss; + pi = curr_fs + curr_ss * num_pix_fs; curr_radius = (int)rint(r_map[pi]); curr_threshold = rthreshold[curr_radius]; if ( copy[pi] > curr_threshold - && pix_in_peak_map[pi] == 0 + && pfinter->pix_in_peak_map[pi] == 0 && mask[pi] != 0 ) { @@ -569,12 +665,15 @@ static void peak_search(int p, int *inss, int *infs, int h, int w, *sum_com_fs += curr_i * ((float)curr_fs); *sum_com_ss += curr_i * ((float)curr_ss); - inss[*num_pix_in_peak] = inss[p] + search_ss[k]; - infs[*num_pix_in_peak] = infs[p] + search_fs[k]; + pfinter->inss[*num_pix_in_peak] = pfinter->inss[p] + + search_ss[k]; + pfinter->infs[*num_pix_in_peak] = pfinter->infs[p] + + search_fs[k]; - pix_in_peak_map[pi] = 1; - if ( *num_pix_in_peak < max_pix_count ) - peak_pixels[*num_pix_in_peak] = pi; + pfinter->pix_in_peak_map[pi] = 1; + if ( *num_pix_in_peak < max_pix_count ) { + pfinter->peak_pixels[*num_pix_in_peak] = pi; + } *num_pix_in_peak = *num_pix_in_peak+1; } } @@ -582,16 +681,17 @@ static void peak_search(int p, int *inss, int *infs, int h, int w, static void search_in_ring(int ring_width, int com_fs_int, int com_ss_int, - int h, int w, int multip_w, float *copy, - float *r_map, float *rthreshold, - float *roffset, char *pix_in_peak_map, - char *mask, float *local_sigma, + float *copy, float *r_map, + float *rthreshold, float *roffset, + char *pix_in_peak_map, char *mask, int asic_size_fs, + int asic_size_ss, int aifs, int aiss, + int num_pix_fs,float *local_sigma, float *local_offset, float *background_max_i, - int com_idx, - int local_bg_radius) + int com_idx, int local_bg_radius) { int ssj, fsi; float pix_radius; + int curr_fs, curr_ss; int pi; int curr_radius; float curr_threshold; @@ -617,26 +717,28 @@ static void search_in_ring(int ring_width, int com_fs_int, int com_ss_int, if ( (com_fs_int + fsi) < 0 ) continue; - if ( (com_fs_int + fsi) >= w ) + if ( (com_fs_int + fsi) >= asic_size_fs ) continue; if ( (com_ss_int + ssj) < 0 ) continue; - if ( (com_ss_int + ssj) >= h ) + if ( (com_ss_int + ssj) >= asic_size_ss ) continue; pix_radius = sqrt(fsi * fsi + ssj * ssj); if ( pix_radius>ring_width ) continue; - pi = (com_fs_int +fsi) + (com_ss_int +ssj) * multip_w; + curr_fs = com_fs_int +fsi + aifs * asic_size_fs; + curr_ss = com_ss_int +ssj + aiss * asic_size_ss; + pi = curr_fs + curr_ss * num_pix_fs; curr_radius = rint(r_map[pi]); curr_threshold = rthreshold[curr_radius]; curr_i = copy[pi]; if ( copy[pi] < curr_threshold - && pix_in_peak_map[pi] == 0 - && mask[pi] != 0 ) { + && pix_in_peak_map[pi] == 0 + && mask[pi] != 0 ) { np_sigma++; sum_i += curr_i; @@ -653,7 +755,7 @@ static void search_in_ring(int ring_width, int com_fs_int, int com_ss_int, if ( np_sigma != 0 ) { *local_offset = sum_i / np_sigma; *local_sigma = sqrt(sum_i_squared / np_sigma - - ((sum_i / np_sigma) * (sum_i / np_sigma))); + ((sum_i / np_sigma) * (sum_i / np_sigma))); } else { local_radius = rint(r_map[(int)rint(com_idx)]); @@ -663,121 +765,75 @@ static void search_in_ring(int ring_width, int com_fs_int, int com_ss_int, } -static int peakfinder8_multipanel(float *roffset, float *rthreshold, - float *data, char *mask, float *r_map, int w, - int h, int max_n_peaks, int *num_found_peaks, - int *npix, float *com_fs, float *com_ss, - float *tot_i, float *max_i, float *sigma, - float *snr, int min_pix_count, - int max_pix_count, int local_bg_radius, - float min_snr, int multip_w) +static void process_panel(int asic_size_fs, int asic_size_ss, int num_pix_fs, + int aiss, int aifs, float *rthreshold, + float *roffset, int *peak_count, + float *copy, struct peakfinder_intern_data *pfinter, + float *r_map, char *mask, int *npix, float *com_fs, + float *com_ss, int *com_index, float *tot_i, + float *max_i, float *sigma, float *snr, + int min_pix_count, int max_pix_count, + int local_bg_radius, float min_snr, int max_n_peaks) { - int i, p, iss, ifs; - int peak_count; - int peak_idx; - int num_pix_in_peak, lt_num_pix_in_pk; - int idx, com_idx; - int curr_radius; - int ring_width; - int curr_idx; - int curr_fs, curr_ss; - int *peak_pixels; - int *infs; - int *inss; - float peak_com_fs, peak_com_ss; - float sum_com_fs, sum_com_ss, sum_i; - float peak_com_fs_int, peak_com_ss_int; - float curr_threshold; - float peak_tot_i, pk_tot_i_raw; - float peak_max_i, pk_max_i_raw; - float local_sigma; - float local_offset; - float background_max_i; - float f_background_thresh; - float peak_snr; - float curr_i, curr_i_raw; - float *copy; - char *pix_in_peak_map; - - copy = (float*) malloc(w*h*sizeof(float)); - if ( copy == NULL ) { - return 1; - } - - memcpy(copy, data, w*h*sizeof(float)); - pix_in_peak_map = calloc(w*h, sizeof(char)); - if ( pix_in_peak_map == NULL ) { - free(copy); - return 1; - } - - infs = (int *) calloc(w*h, sizeof(int)); - if ( infs == NULL ) { - free(pix_in_peak_map); - free(copy); - return 1; - } - - inss = (int *) calloc(w*h, sizeof(int)); - if ( inss == NULL ) { - free(infs); - free(pix_in_peak_map); - free(copy); - return 1; - } - - // TODO Possible bug: should 0,0 be included? - peak_pixels = calloc(max_pix_count, sizeof(int)); - if ( peak_pixels == NULL ) { - free(inss); - free(infs); - free(pix_in_peak_map); - free(copy); - return 1; - } - - for ( i = 0; i < w*h; i++ ) { - copy[i] *= mask[i]; - } - - peak_count = 0; - num_pix_in_peak = 0; - - for ( iss=0 ; iss curr_threshold && - pix_in_peak_map[idx] == 0 ) { - - infs[0] = ifs; - inss[0] = iss; - - peak_pixels[0] = idx; + int pxss, pxfs; + int num_pix_in_peak; + + for ( pxss=1 ; pxss curr_thresh + && pfinter->pix_in_peak_map[pxidx] == 0 ) { + + float sum_com_fs, sum_com_ss; + float sum_i; + float peak_com_fs, peak_com_ss; + float peak_com_fs_int, peak_com_ss_int; + float peak_tot_i, pk_tot_i_raw; + float peak_max_i, pk_max_i_raw; + float peak_snr; + float local_sigma, local_offset; + float background_max_i; + float f_background_thresh; + int lt_num_pix_in_pk; + int ring_width; + int peak_idx; + int com_idx; + int p; + + pfinter->infs[0] = pxfs; + pfinter->inss[0] = pxss; + pfinter->peak_pixels[0] = pxidx; num_pix_in_peak = 1; + sum_i = 0; sum_com_fs = 0; sum_com_ss = 0; - do { - + do { lt_num_pix_in_pk = num_pix_in_peak; for ( p=0; ppix_in_peak_map, + mask, asic_size_fs, + asic_size_ss, + aifs, aiss, + num_pix_fs, + &local_sigma, &local_offset, &background_max_i, com_idx, local_bg_radius); @@ -821,11 +886,22 @@ static int peakfinder8_multipanel(float *roffset, float *rthreshold, sum_com_ss = 0; for ( peak_idx = 1 ; - peak_idx < num_pix_in_peak && - peak_idx <= max_pix_count ; + peak_idx < num_pix_in_peak + && peak_idx <= max_pix_count ; peak_idx++ ) { - curr_idx = peak_pixels[peak_idx]; + int curr_idx; + float curr_i; + float curr_i_raw; + int curr_fs, curr_ss; + + // BUG HERE, I THINK. PEAK_PIXELS + // HAS SIZE MAX_PIX_COUNT, BUT + // IN THE FOLLOWING LINE PEAK_IDX + // CAN HAVE VALUE EXACTLY MAX_PEAK_COUNT + // (SEE THE FOR LOOP ABOVE) + curr_idx = + pfinter->peak_pixels[peak_idx]; curr_i_raw = copy[curr_idx]; curr_i = curr_i_raw - local_offset; @@ -833,8 +909,8 @@ static int peakfinder8_multipanel(float *roffset, float *rthreshold, peak_tot_i += curr_i; pk_tot_i_raw += curr_i_raw; - curr_fs = curr_idx % multip_w; - curr_ss = curr_idx / multip_w; + curr_fs = curr_idx % asic_size_fs; + curr_ss = curr_idx / asic_size_fs; sum_com_fs += curr_i * ((float)curr_fs); sum_com_ss += curr_i * ((float)curr_ss); @@ -844,10 +920,11 @@ static int peakfinder8_multipanel(float *roffset, float *rthreshold, peak_max_i = curr_i; } + peak_com_fs = sum_com_fs / fabs(peak_tot_i); peak_com_ss = sum_com_ss / fabs(peak_tot_i); - peak_snr = peak_tot_i / local_sigma; + peak_snr = peak_tot_i / local_sigma; if (peak_snr < min_snr) { continue; @@ -856,70 +933,113 @@ static int peakfinder8_multipanel(float *roffset, float *rthreshold, f_background_thresh = 1; f_background_thresh *= background_max_i - local_offset; - if (peak_max_i < f_background_thresh) + if (peak_max_i < f_background_thresh) { continue; + } if ( num_pix_in_peak >= min_pix_count && num_pix_in_peak <= max_pix_count ) { - if ( peak_tot_i == 0 ) + int peak_com_idx; + + if ( peak_tot_i == 0 ) { continue; + } + + peak_com_idx = rint(peak_com_fs) + + rint(peak_com_ss) * + asic_size_fs; - if ( peak_count < max_n_peaks ) { + if ( *peak_count < max_n_peaks ) { int pidx; - pidx = peak_count; + pidx = *peak_count; npix[pidx] = num_pix_in_peak; com_fs[pidx] = peak_com_fs; com_ss[pidx] = peak_com_ss; + com_index[pidx] = peak_com_idx; tot_i[pidx] = peak_tot_i; max_i[pidx] = peak_max_i; sigma[pidx] = local_sigma; snr[pidx] = peak_snr; - peak_count++; - } - else { - peak_count++; + *peak_count = *peak_count + 1; + } else { + *peak_count = *peak_count + 1; } } } } } +} + + +static int peakfinder8_base(float *roffset, float *rthreshold, + float *data, char *mask, float *r_map, + int asic_size_fs, int num_asics_fs, + int asic_size_ss, int num_asics_ss, + int max_n_peaks, int *num_found_peaks, + int *npix, float *com_fs, + float *com_ss, int *com_index, float *tot_i, + float *max_i, float *sigma, float *snr, + int min_pix_count, int max_pix_count, + int local_bg_radius, float min_snr) +{ + int num_pix_fs, num_pix_ss, num_pix_tot; + int i, aifs, aiss; + int peak_count; + float *copy; + struct peakfinder_intern_data *pfinter; + + num_pix_fs = asic_size_fs * num_asics_fs; + num_pix_ss = asic_size_ss * num_asics_ss; + num_pix_tot = num_pix_fs * num_pix_ss; + + copy = (float *)calloc(num_pix_tot, sizeof(float)); + if ( copy == NULL ) { + return 1; + } + + memcpy(copy, data, num_pix_tot*sizeof(float)); + + for (i = 0; i < num_pix_tot; i++) { + copy[i] *= mask[i]; + } + + pfinter = allocate_peakfinder_intern_data(num_pix_tot, max_pix_count); + if ( pfinter == NULL ) { + free(copy); + return 1; + } + + peak_count = 0; + + for ( aiss=0 ; aissdet->n_panels); if ( pfdata == NULL) { free_radius_maps(rmaps); free_peakfinder_mask(pfmask); @@ -960,6 +1081,7 @@ int peakfinder8(struct image *img, int max_n_peaks, pfdata->panel_h[pi] = img->det->panels[pi].h; pfdata->panel_w[pi] = img->det->panels[pi].w; pfdata->panel_data[pi] = img->dp[pi]; + pfdata->num_panels = img->det->n_panels; } num_rad_bins = compute_num_radial_bins(img->det->n_panels, @@ -975,18 +1097,40 @@ int peakfinder8(struct image *img, int max_n_peaks, return 1; } - ret = compute_radial_stats(pfdata->panel_data, pfdata->num_panels, - pfdata->panel_w, pfdata->panel_h, - rmaps->r_maps, - pfmask->masks, rstats->rthreshold, - rstats->roffset, num_rad_bins, - min_snr, threshold, iterations); - if ( ret != 0 ) { - free_radius_maps(rmaps); - free_peakfinder_mask(pfmask); - free_panel_data(pfdata); - free_radial_stats(rstats); - return 1; + for ( i=0 ; in_rad_bins ; i++) { + rstats->rthreshold[i] = 1e9; + } + + for ( it_counter=0 ; it_counterroffset[i] = 0; + rstats->rsigma[i] = 0; + rstats->rcount[i] = 0; + } + + for ( pi=0 ; pinum_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->roffset, + rstats->rsigma, + rstats->rcount); + + } + + compute_radial_stats(rstats->rthreshold, + rstats->roffset, + rstats->rsigma, + rstats->rcount, + num_rad_bins, + min_snr, + threshold); + } pkdata = allocate_peak_data(max_n_peaks); @@ -998,10 +1142,13 @@ int peakfinder8(struct image *img, int max_n_peaks, return 1; } + remaining_max_num_peaks = max_n_peaks; + for ( pi=0 ; pidet->n_panels ; pi++) { int peaks_to_add; int pki; + int ret; num_found_peaks = 0; @@ -1009,34 +1156,43 @@ int peakfinder8(struct image *img, int max_n_peaks, continue; } - ret = peakfinder8_singlepanel(rstats->roffset, - rstats->rthreshold, - pfdata->panel_data[pi], - pfmask->masks[pi], - rmaps->r_maps[pi], - pfdata->panel_w[pi], - pfdata->panel_h[pi], - max_n_peaks, - &num_found_peaks, - pkdata->npix, - pkdata->com_fs, - pkdata->com_ss, - pkdata->tot_i, - pkdata->max_i, - pkdata->sigma, - pkdata->snr, - min_pix_count, - max_pix_count, - local_bg_radius, - min_snr); + ret = peakfinder8_base(rstats->roffset, + rstats->rthreshold, + pfdata->panel_data[pi], + pfmask->masks[pi], + rmaps->r_maps[pi], + pfdata->panel_w[pi], 1, + pfdata->panel_h[pi], 1, + max_n_peaks, + &num_found_peaks, + pkdata->npix, + pkdata->com_fs, + pkdata->com_ss, + pkdata->com_index, + pkdata->tot_i, + pkdata->max_i, + pkdata->sigma, + pkdata->snr, + min_pix_count, + max_pix_count, + local_bg_radius, + min_snr); + + if ( ret != 0 ) { + free_radius_maps(rmaps); + free_peakfinder_mask(pfmask); + free_panel_data(pfdata); + free_radial_stats(rstats); + return 1; + } peaks_to_add = num_found_peaks; - if ( num_found_peaks > max_n_peaks ) { - peaks_to_add = max_n_peaks; + if ( num_found_peaks > remaining_max_num_peaks ) { + peaks_to_add = remaining_max_num_peaks; } - max_n_peaks -= peaks_to_add; + remaining_max_num_peaks -= peaks_to_add; for ( pki=0 ; pkitot_i[pki], NULL); - } } -- cgit v1.2.3 From 425619c37c3973cd354a667a61e84ccab7d81572 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Thu, 23 Mar 2017 14:59:13 +0100 Subject: Mostly fussiness --- libcrystfel/src/peakfinder8.c | 34 ++++++++++++++++++++-------------- libcrystfel/src/peakfinder8.h | 20 ++++++++++++++------ 2 files changed, 34 insertions(+), 20 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/peakfinder8.c b/libcrystfel/src/peakfinder8.c index 0e916a0c..f6e70223 100644 --- a/libcrystfel/src/peakfinder8.c +++ b/libcrystfel/src/peakfinder8.c @@ -28,6 +28,12 @@ * */ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include #include "peakfinder8.h" @@ -58,7 +64,7 @@ struct radial_stats struct peakfinder_panel_data { - float ** panel_data; + float **panel_data; int *panel_h; int *panel_w; int num_panels; @@ -88,7 +94,7 @@ struct peakfinder_peak_data }; -static double compute_r_sigma(float *rsigma, int *rcount, float* roffset, +static double compute_r_sigma(float *rsigma, int *rcount, float *roffset, int i) { return sqrt(rsigma[i] / rcount[i] - @@ -97,8 +103,8 @@ static double compute_r_sigma(float *rsigma, int *rcount, float* roffset, } -static void update_radial_stats (float *roffset, float *rsigma, int *rcount, - float value, int curr_radius) +static void update_radial_stats(float *roffset, float *rsigma, int *rcount, + float value, int curr_radius) { roffset[curr_radius] += value; rsigma[curr_radius] += (value * value); @@ -118,7 +124,7 @@ static float get_radius(struct panel p, float fs, float ss) } -static struct radius_maps* compute_radius_maps(struct detector *det) +static struct radius_maps *compute_radius_maps(struct detector *det) { int i, u, iss, ifs; struct panel p; @@ -237,7 +243,7 @@ static void free_peakfinder_mask(struct peakfinder_mask * pfmask) static int compute_num_radial_bins(int num_panels, int *w, int *h, - float **r_maps ) + float **r_maps) { float max_r; @@ -267,7 +273,7 @@ static int compute_num_radial_bins(int num_panels, int *w, int *h, } -static struct peakfinder_panel_data* allocate_panel_data(int num_panels) +static struct peakfinder_panel_data *allocate_panel_data(int num_panels) { struct peakfinder_panel_data *pfdata; @@ -315,9 +321,9 @@ static void free_panel_data(struct peakfinder_panel_data *pfdata) } -static struct radial_stats* allocate_radial_stats(int num_rad_bins) +static struct radial_stats *allocate_radial_stats(int num_rad_bins) { - struct radial_stats* rstats; + struct radial_stats *rstats; rstats = (struct radial_stats *)malloc(sizeof(struct radial_stats)); if ( rstats == NULL ) { @@ -545,7 +551,8 @@ struct peakfinder_peak_data *allocate_peak_data(int max_num_peaks) } -static void free_peak_data(struct peakfinder_peak_data *pkdata) { +static void free_peak_data(struct peakfinder_peak_data *pkdata) +{ free(pkdata->npix); free(pkdata->com_fs); free(pkdata->com_ss); @@ -614,9 +621,7 @@ static void free_peakfinder_intern_data(struct peakfinder_intern_data *pfid) } - -static void peak_search(int p, - struct peakfinder_intern_data *pfinter, +static void peak_search(int p, struct peakfinder_intern_data *pfinter, float *copy, char *mask, float *r_map, float *rthreshold, float *roffset, int *num_pix_in_peak, int asic_size_fs, @@ -656,7 +661,8 @@ static void peak_search(int p, if ( copy[pi] > curr_threshold && pfinter->pix_in_peak_map[pi] == 0 - && mask[pi] != 0 ) { + && mask[pi] != 0 ) + { curr_i = copy[pi] - roffset[curr_radius]; diff --git a/libcrystfel/src/peakfinder8.h b/libcrystfel/src/peakfinder8.h index 1103ef96..483eebdf 100644 --- a/libcrystfel/src/peakfinder8.h +++ b/libcrystfel/src/peakfinder8.h @@ -1,7 +1,7 @@ /* * peakfinder8.h * - * The processing pipeline for one image + * The peakfinder8 algorithm * * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. @@ -33,10 +33,18 @@ #include "image.h" -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); +#ifdef __cplusplus +extern "C" { +#endif + +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); + +#ifdef __cplusplus +} +#endif #endif // PEAKFINDER8_H -- cgit v1.2.3 From c5ab3a8e0305995ed97fc85a5c190c069794fdb0 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 24 Mar 2017 15:04:21 +0100 Subject: Propagate error if indexing method is unrecognised --- libcrystfel/src/index.c | 4 ++++ 1 file changed, 4 insertions(+) (limited to 'libcrystfel') diff --git a/libcrystfel/src/index.c b/libcrystfel/src/index.c index 028ba5e9..685d0b69 100644 --- a/libcrystfel/src/index.c +++ b/libcrystfel/src/index.c @@ -163,6 +163,10 @@ IndexingPrivate *setup_indexing(const char *method_list, UnitCell *cell, for ( i=0; i Date: Fri, 24 Mar 2017 15:04:45 +0100 Subject: Remove an old debugging message --- libcrystfel/src/detector.c | 1 - 1 file changed, 1 deletion(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/detector.c b/libcrystfel/src/detector.c index 52e4c330..e160d1e1 100644 --- a/libcrystfel/src/detector.c +++ b/libcrystfel/src/detector.c @@ -909,7 +909,6 @@ static int parse_field_for_panel(struct panel *panel, const char *key, panel->adu_per_eV = atof(val); } else if ( strcmp(key, "adu_per_photon") == 0 ) { panel->adu_per_photon = atof(val); - STATUS("got adu per photon: %s\n", val); } else if ( strcmp(key, "rigid_group") == 0 ) { add_to_rigid_group(find_or_add_rg(det, val), panel); } else if ( strcmp(key, "clen") == 0 ) { -- cgit v1.2.3 From 3746b4dbed15569764ec6de473806468951748b8 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 31 Mar 2017 16:27:13 +0200 Subject: Offset peak locations from HDF5 or CXI files by 0.5,0.5 CrystFEL considers all peak locations to be distances from the corner of the detector panel, in pixel units, consistent with its description of detector geometry. In contrast, Cheetah considers the peak locations to be pixel indices in the data array. Therefore, a half-pixel offset is needed when importing the peak lists. For users who need the old behaviour, this commit adds a new option indexamajig --no-half-pixel-shift to deactivate this offset. --- libcrystfel/src/hdf5-file.c | 123 +++++++++++++++++++++++++++++++++++++++++--- libcrystfel/src/hdf5-file.h | 11 +++- 2 files changed, 124 insertions(+), 10 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/hdf5-file.c b/libcrystfel/src/hdf5-file.c index 9bcdb1b7..76617ee4 100644 --- a/libcrystfel/src/hdf5-file.c +++ b/libcrystfel/src/hdf5-file.c @@ -44,6 +44,19 @@ #include "utils.h" +/** + * SECTION:hdf5-file + * @short_description: HDF5 file handling + * @title: HDF5 file handling + * @section_id: + * @see_also: + * @include: "hdf5-file.h" + * @Image: + * + * Routines for accessing HDF5 files. + **/ + + struct hdf5_write_location { const char *location; @@ -320,9 +333,34 @@ static float *read_hdf5_data(struct hdfile *f, char *path, int line) } -/* Get peaks from HDF5, in "CXI format" (as in "CXIDB") */ -int get_peaks_cxi(struct image *image, struct hdfile *f, const char *p, - struct filename_plus_event *fpe) +/** + * get_peaks_cxi_2: + * @image: An %image structure + * @f: An %hdfile structure + * @p: The HDF5 path to the peak data + * @fpe: A %filename_plus_event structure specifying the event + * @half_pixel_shift: Non-zero if 0.5 should be added to all peak coordinates + * + * Get peaks from HDF5, in "CXI format" (as in "CXIDB"). The data should be in + * a set of arrays under @p. The number of peaks should be in a 1D array at + * @p/nPeaks. The fast-scan and slow-scan coordinates should be in 2D arrays at + * @p/peakXPosRaw and @p/peakYPosRaw respectively (sorry about the naming). The + * first dimension of these arrays should be the event number (as given by + * @fpe). The intensities are expected to be at @p/peakTotalIntensity in a + * similar 2D array. + * + * CrystFEL considers all peak locations to be distances from the corner of the + * detector panel, in pixel units, consistent with its description of detector + * geometry (see 'man crystfel_geometry'). The software which generates the + * CXI files, including Cheetah, may instead consider the peak locations to be + * pixel indices in the data array. In this case, the peak coordinates should + * have 0.5 added to them. This will be done if @half_pixel_shift is non-zero. + * + * Returns: non-zero on error, zero otherwise. + * + */ +int get_peaks_cxi_2(struct image *image, struct hdfile *f, const char *p, + struct filename_plus_event *fpe, int half_pixel_shift) { char path_n[1024]; char path_x[1024]; @@ -338,6 +376,8 @@ int get_peaks_cxi(struct image *image, struct hdfile *f, const char *p, float *buf_y; float *buf_i; + double peak_offset = half_pixel_shift ? 0.5 : 0.0; + if ( (fpe != NULL) && (fpe->ev != NULL) && (fpe->ev->dim_entries != NULL) ) { @@ -375,8 +415,8 @@ int get_peaks_cxi(struct image *image, struct hdfile *f, const char *p, float fs, ss, val; struct panel *p; - fs = buf_x[pk]; - ss = buf_y[pk]; + fs = buf_x[pk] + peak_offset; + ss = buf_y[pk] + peak_offset; val = buf_i[pk]; p = find_orig_panel(image->det, fs, ss); @@ -395,7 +435,53 @@ int get_peaks_cxi(struct image *image, struct hdfile *f, const char *p, } -int get_peaks(struct image *image, struct hdfile *f, const char *p) +/** + * get_peaks_cxi: + * @image: An %image structure + * @f: An %hdfile structure + * @p: The HDF5 path to the peak data + * @fpe: A %filename_plus_event structure specifying the event + * + * This is a wrapper function to preserve API compatibility with older CrystFEL + * versions. Use get_peaks_cxi_2() instead. + * + * This function is equivalent to get_peaks_cxi_2(@image, @f, @p, @fpe, 1). + * + * Returns: non-zero on error, zero otherwise. + * + */ +int get_peaks_cxi(struct image *image, struct hdfile *f, const char *p, + struct filename_plus_event *fpe) +{ + return get_peaks_cxi_2(image, f, p, fpe, 1); +} + + +/** + * get_peaks_2: + * @image: An %image structure + * @f: An %hdfile structure + * @p: The HDF5 path to the peak data + * @half_pixel_shift: Non-zero if 0.5 should be added to all peak coordinates + * + * Get peaks from HDF5. The peak list should be located at @p in the HDF5 file, + * a 2D array where the first dimension equals the number of peaks and second + * dimension is three. The first two columns contain the fast scan and slow + * scan coordinates, respectively, of the peaks. The third column contains the + * intensities. + * + * CrystFEL considers all peak locations to be distances from the corner of the + * detector panel, in pixel units, consistent with its description of detector + * geometry (see 'man crystfel_geometry'). The software which generates the + * CXI files, including Cheetah, may instead consider the peak locations to be + * pixel indices in the data array. In this case, the peak coordinates should + * have 0.5 added to them. This will be done if @half_pixel_shift is non-zero. + * + * Returns: non-zero on error, zero otherwise. + * + */ +int get_peaks_2(struct image *image, struct hdfile *f, const char *p, + int half_pixel_shift) { hid_t dh, sh; hsize_t size[2]; @@ -405,6 +491,7 @@ int get_peaks(struct image *image, struct hdfile *f, const char *p) herr_t r; int tw; char *np; + double peak_offset = half_pixel_shift ? 0.5 : 0.0; if ( image->event != NULL ) { np = retrieve_full_path(image->event, p); @@ -473,8 +560,8 @@ int get_peaks(struct image *image, struct hdfile *f, const char *p) float fs, ss, val; struct panel *p; - fs = buf[tw*i+0]; - ss = buf[tw*i+1]; + fs = buf[tw*i+0] + peak_offset; + ss = buf[tw*i+1] + peak_offset; val = buf[tw*i+2]; p = find_orig_panel(image->det, fs, ss); @@ -499,6 +586,26 @@ int get_peaks(struct image *image, struct hdfile *f, const char *p) } +/** + * get_peaks: + * @image: An %image structure + * @f: An %hdfile structure + * @p: The HDF5 path to the peak data + * + * This is a wrapper function to preserve API compatibility with older CrystFEL + * versions. Use get_peaks_2() instead. + * + * This function is equivalent to get_peaks_2(@image, @f, @p, 1). + * + * Returns: non-zero on error, zero otherwise. + * + */ +int get_peaks(struct image *image, struct hdfile *f, const char *p) +{ + return get_peaks_2(image, f, p, 1); +} + + static void cleanup(hid_t fh) { int n_ids, i; diff --git a/libcrystfel/src/hdf5-file.h b/libcrystfel/src/hdf5-file.h index 8c89eb93..f22d4d54 100644 --- a/libcrystfel/src/hdf5-file.h +++ b/libcrystfel/src/hdf5-file.h @@ -3,11 +3,11 @@ * * Read/write HDF5 data files * - * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * * Authors: - * 2009-2015 Thomas White + * 2009-2017 Thomas White * 2014 Valerio Mariani * @@ -78,9 +78,16 @@ extern void hdfile_close(struct hdfile *f); extern int get_peaks(struct image *image, struct hdfile *f, const char *p); +extern int get_peaks_2(struct image *image, struct hdfile *f, const char *p, + int half_pixel_shift); + extern int get_peaks_cxi(struct image *image, struct hdfile *f, const char *p, struct filename_plus_event *fpe); +extern int get_peaks_cxi_2(struct image *image, struct hdfile *f, const char *p, + struct filename_plus_event *fpe, + int half_pixel_shift); + extern struct copy_hdf5_field *new_copy_hdf5_field_list(void); extern void free_copy_hdf5_field_list(struct copy_hdf5_field *f); -- cgit v1.2.3 From adaff63c98ce5a0e138999a5ee8410c9404f6c29 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Thu, 13 Apr 2017 13:10:20 +0200 Subject: hdf5_read2(): Free buffers when read fails --- libcrystfel/src/hdf5-file.c | 28 +++++++++++++++++++++++++++- 1 file changed, 27 insertions(+), 1 deletion(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/hdf5-file.c b/libcrystfel/src/hdf5-file.c index 76617ee4..c36a58b3 100644 --- a/libcrystfel/src/hdf5-file.c +++ b/libcrystfel/src/hdf5-file.c @@ -1741,6 +1741,9 @@ int hdf5_read2(struct hdfile *f, struct image *image, struct event *ev, if ( !exists ) { ERROR("Cannot find data for panel %s\n", p->name); + free(image->dp); + free(image->bad); + free(image->sat); return 1; } @@ -1761,6 +1764,9 @@ int hdf5_read2(struct hdfile *f, struct image *image, struct event *ev, if ( !exists ) { ERROR("Cannot find data for panel %s\n", p->name); + free(image->dp); + free(image->bad); + free(image->sat); return 1; } fail = hdfile_set_image(f, p->data, p); @@ -1771,6 +1777,9 @@ int hdf5_read2(struct hdfile *f, struct image *image, struct event *ev, if ( fail ) { ERROR("Couldn't select path for panel %s\n", p->name); + free(image->dp); + free(image->bad); + free(image->sat); return 1; } @@ -1780,6 +1789,9 @@ int hdf5_read2(struct hdfile *f, struct image *image, struct event *ev, f_count = malloc(hsd->num_dims*sizeof(hsize_t)); if ( (f_offset == NULL) || (f_count == NULL ) ) { ERROR("Failed to allocate offset or count.\n"); + free(image->dp); + free(image->bad); + free(image->sat); return 1; } for ( hsi=0; hsinum_dims; hsi++ ) { @@ -1807,6 +1819,9 @@ int hdf5_read2(struct hdfile *f, struct image *image, struct event *ev, if ( check < 0 ) { ERROR("Error selecting file dataspace for panel %s\n", p->name); + free(image->dp); + free(image->bad); + free(image->sat); return 1; } @@ -1820,6 +1835,9 @@ int hdf5_read2(struct hdfile *f, struct image *image, struct event *ev, ERROR("Failed to allocate panel %s\n", p->name); free(f_offset); free(f_count); + free(image->dp); + free(image->bad); + free(image->sat); return 1; } for ( i=0; iw*p->h; i++ ) image->sat[pi][i] = INFINITY; @@ -1831,6 +1849,14 @@ int hdf5_read2(struct hdfile *f, struct image *image, struct event *ev, p->name); free(f_offset); free(f_count); + for ( i=0; i<=pi; i++ ) { + free(image->dp[pi]); + free(image->sat[pi]); + free(image->bad[pi]); + } + free(image->dp); + free(image->bad); + free(image->sat); return 1; } @@ -1847,7 +1873,7 @@ int hdf5_read2(struct hdfile *f, struct image *image, struct event *ev, if ( load_satmap(f, ev, p, f_offset, f_count, hsd, image->sat[pi]) ) { - ERROR("Failed to laod sat map for panel %s\n", + ERROR("Failed to load sat map for panel %s\n", p->name); } } -- cgit v1.2.3 From 4227f0b190b08ecc50a49875e86dbb9b14714bdd Mon Sep 17 00:00:00 2001 From: Thomas White Date: Thu, 13 Apr 2017 16:09:58 +0200 Subject: Add half-pixel offset for peakfinder8 Like zaef, the pixel indices need to be converted to geometrical coordinates. This increases the indexing rate by about 5% in my test. --- libcrystfel/src/peakfinder8.c | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/peakfinder8.c b/libcrystfel/src/peakfinder8.c index f6e70223..76d7c13c 100644 --- a/libcrystfel/src/peakfinder8.c +++ b/libcrystfel/src/peakfinder8.c @@ -1215,8 +1215,8 @@ int peakfinder8(struct image *img, int max_n_peaks, } image_add_feature(img->features, - pkdata->com_fs[pki], - pkdata->com_ss[pki], + pkdata->com_fs[pki]+0.5, + pkdata->com_ss[pki]+0.5, p, img, pkdata->tot_i[pki], -- cgit v1.2.3 From 8fbfaf71b1efef4bfdb40ce85200e772e82e9773 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Tue, 14 Mar 2017 14:14:06 +0100 Subject: Initial CBF stuff --- libcrystfel/src/image.c | 4 ++-- libcrystfel/src/image.h | 17 +++++++++++++++-- libcrystfel/src/stream.c | 6 +++--- libcrystfel/src/stream.h | 9 +++++++-- 4 files changed, 27 insertions(+), 9 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/image.c b/libcrystfel/src/image.c index 09f4958d..2cc0d792 100644 --- a/libcrystfel/src/image.c +++ b/libcrystfel/src/image.c @@ -3,12 +3,12 @@ * * Handle images and image features * - * Copyright © 2012-2016 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * * Authors: * 2014 Kenneth Beyerlein - * 2011-2016 Thomas White + * 2011-2017 Thomas White * * This file is part of CrystFEL. * diff --git a/libcrystfel/src/image.h b/libcrystfel/src/image.h index 9fd9b495..2ed7140f 100644 --- a/libcrystfel/src/image.h +++ b/libcrystfel/src/image.h @@ -3,11 +3,11 @@ * * Handle images and image features * - * Copyright © 2012-2016 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * * Authors: - * 2009-2016 Thomas White + * 2009-2017 Thomas White * 2014 Valerio Mariani * * @@ -44,6 +44,8 @@ struct detector; struct imagefeature; struct sample; struct image; +struct imagefile; +struct imagefile_field_list; #include "utils.h" #include "cell.h" @@ -51,6 +53,7 @@ struct image; #include "reflist.h" #include "crystal.h" #include "index.h" +#include "events.h" /** * SpectrumType: @@ -233,6 +236,16 @@ extern void image_add_crystal(struct image *image, Crystal *cryst); extern void remove_flagged_crystals(struct image *image); extern void free_all_crystals(struct image *image); +/* Image files */ +extern struct imagefile *imagefile_open(const char *filename); +extern int imagefile_read(struct imagefile *imfile, struct image *image, + struct event *event); +extern struct hdfile *imagefile_get_hdfile(struct imagefile *imfile); +extern void imagefile_copy_fields(struct imagefile *imfile, + struct copy_hdf5_file *copyme, FILE *fh, + struct event *ev); +extern void imagefile_close(struct imagefile *imfile); + #ifdef __cplusplus } #endif diff --git a/libcrystfel/src/stream.c b/libcrystfel/src/stream.c index 17da74b2..fb4b0c70 100644 --- a/libcrystfel/src/stream.c +++ b/libcrystfel/src/stream.c @@ -806,8 +806,8 @@ static int write_crystal(Stream *st, Crystal *cr, int include_reflections) } -int write_chunk(Stream *st, struct image *i, struct hdfile *hdfile, - int include_peaks, int include_reflections, struct event* ev) +int write_chunk(Stream *st, struct image *i, struct imagefile *imfile, + int include_peaks, int include_reflections, struct event *ev) { int j; char *indexer; @@ -832,7 +832,7 @@ int write_chunk(Stream *st, struct image *i, struct hdfile *hdfile, fprintf(st->fh, "beam_divergence = %.2e rad\n", i->div); fprintf(st->fh, "beam_bandwidth = %.2e (fraction)\n", i->bw); - copy_hdf5_fields(hdfile, i->copyme, st->fh, ev); + imagefile_copy_fields(imfile, i->copyme, st->fh, ev); if ( i->det != NULL ) { diff --git a/libcrystfel/src/stream.h b/libcrystfel/src/stream.h index ff8628c0..a95c7df0 100644 --- a/libcrystfel/src/stream.h +++ b/libcrystfel/src/stream.h @@ -3,11 +3,11 @@ * * Stream tools * - * Copyright © 2013-2014 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2013-2017 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * * Authors: - * 2010-2014 Thomas White + * 2010-2017 Thomas White * 2014 Valerio Mariani * 2011 Andrew Aquila * @@ -110,6 +110,11 @@ extern int write_chunk(Stream *st, struct image *image, struct hdfile *hdfile, int include_peaks, int include_reflections, struct event *ev); +extern int write_chunk_2(Stream *st, struct image *image, + struct imagefile *imfile, + int include_peaks, int include_reflections, + struct event *ev); + extern void write_command(Stream *st, int argc, char *argv[]); extern void write_geometry_file(Stream *st, const char *geom_filename); extern int rewind_stream(Stream *st); -- cgit v1.2.3 From cb6389ae61e8f7e279ea16f8ab1a94969d6c0dc9 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 28 Apr 2017 17:15:13 +0200 Subject: Skeleton image file API --- libcrystfel/src/image.c | 129 +++++++++++++++++++++++++++++++++++++++++++++++ libcrystfel/src/image.h | 19 ++++--- libcrystfel/src/stream.h | 3 +- 3 files changed, 144 insertions(+), 7 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/image.c b/libcrystfel/src/image.c index 2cc0d792..f322cbdd 100644 --- a/libcrystfel/src/image.c +++ b/libcrystfel/src/image.c @@ -299,3 +299,132 @@ void free_all_crystals(struct image *image) free(image->crystals); image->n_crystals = 0; } + + +/**************************** Image field lists *******************************/ + +struct imagefile_field_list +{ + char **fields; + int n_fields; + int max_fields; +}; + + +struct imagefile_field_list *new_imagefile_field_list() +{ + struct imagefile_field_list *n; + + n = calloc(1, sizeof(struct imagefile_field_list)); + if ( n == NULL ) return NULL; + + n->max_fields = 32; + n->fields = malloc(n->max_fields*sizeof(char *)); + if ( n->fields == NULL ) { + free(n); + return NULL; + } + + return n; +} + + +void free_imagefile_field_list(struct imagefile_field_list *n) +{ + int i; + for ( i=0; in_fields; i++ ) { + free(n->fields[i]); + } + free(n->fields); + free(n); +} + + +void add_imagefile_field(struct imagefile_field_list *copyme, const char *name) +{ + int i; + + /* Already on the list? Don't re-add if so. */ + for ( i=0; in_fields; i++ ) { + if ( strcmp(copyme->fields[i], name) == 0 ) return; + } + + /* Need more space? */ + if ( copyme->n_fields == copyme->max_fields ) { + + char **nfields; + int nmax = copyme->max_fields + 32; + + nfields = realloc(copyme->fields, nmax*sizeof(char *)); + if ( nfields == NULL ) { + ERROR("Failed to allocate space for new HDF5 field.\n"); + return; + } + + copyme->max_fields = nmax; + copyme->fields = nfields; + + } + + copyme->fields[copyme->n_fields] = strdup(name); + if ( copyme->fields[copyme->n_fields] == NULL ) { + ERROR("Failed to add field for copying '%s'\n", name); + return; + } + + copyme->n_fields++; +} + + +/****************************** Image files ***********************************/ + +struct imagefile *imagefile_open(const char *filename) +{ +} + + +int imagefile_read(struct imagefile *imfile, struct image *image, + struct event *event) +{ +} + + +struct hdfile *imagefile_get_hdfile(struct imagefile *imfile) +{ +} + + +void imagefile_copy_fields(struct imagefile *f, + const struct imagefile_field_list *copyme, + FILE *fh, struct event *ev) +{ + int i; + + if ( copyme == NULL ) return; + + for ( i=0; in_fields; i++ ) { + +#warning FIXME Implement imagefile field copying +#if 0 + char *val; + char *field; + + field = copyme->fields[i]; + val = hdfile_get_string_value(f, field, ev); + + if ( field[0] == '/' ) { + fprintf(fh, "hdf5%s = %s\n", field, val); + } else { + fprintf(fh, "hdf5/%s = %s\n", field, val); + } + + free(val); +#endif + + } +} + + +void imagefile_close(struct imagefile *imfile) +{ +} diff --git a/libcrystfel/src/image.h b/libcrystfel/src/image.h index 2ed7140f..f310a09b 100644 --- a/libcrystfel/src/image.h +++ b/libcrystfel/src/image.h @@ -123,7 +123,7 @@ struct beam_params * struct detector *det; * struct beam_params *beam; * char *filename; - * const struct copy_hdf5_field *copyme; + * const struct imagefile_field_list *copyme; * * int id; * @@ -155,8 +155,8 @@ struct beam_params * returned by the low-level indexing system. n_crystals * is the number of crystals which were found in the image. * - * copyme represents a list of HDF5 fields to copy - * to the output stream. + * copyme represents a list of fields in the image + * file (e.g. HDF5 fields or CBF headers) to copy to the output stream. **/ struct image; @@ -174,7 +174,7 @@ struct image { struct beam_params *beam; /* The nominal beam parameters */ char *filename; struct event *event; - const struct copy_hdf5_field *copyme; + const struct imagefile_field_list *copyme; struct stuff_from_stream *stuff_from_stream; double avg_clen; /* Average camera length extracted @@ -242,10 +242,17 @@ extern int imagefile_read(struct imagefile *imfile, struct image *image, struct event *event); extern struct hdfile *imagefile_get_hdfile(struct imagefile *imfile); extern void imagefile_copy_fields(struct imagefile *imfile, - struct copy_hdf5_file *copyme, FILE *fh, - struct event *ev); + const struct imagefile_field_list *copyme, + FILE *fh, struct event *ev); extern void imagefile_close(struct imagefile *imfile); +/* Field lists */ +extern struct imagefile_field_list *new_imagefile_field_list(void); +extern void free_imagefile_field_list(struct imagefile_field_list *f); + +extern void add_imagefile_field(struct imagefile_field_list *copyme, + const char *name); + #ifdef __cplusplus } #endif diff --git a/libcrystfel/src/stream.h b/libcrystfel/src/stream.h index a95c7df0..764e3e36 100644 --- a/libcrystfel/src/stream.h +++ b/libcrystfel/src/stream.h @@ -39,6 +39,7 @@ struct image; struct hdfile; struct event; +struct imagefile; #include "cell.h" #define GEOM_START_MARKER "----- Begin geometry file -----" @@ -106,7 +107,7 @@ extern int read_chunk(Stream *st, struct image *image); extern int read_chunk_2(Stream *st, struct image *image, StreamReadFlags srf); -extern int write_chunk(Stream *st, struct image *image, struct hdfile *hdfile, +extern int write_chunk(Stream *st, struct image *image, struct imagefile *imfile, int include_peaks, int include_reflections, struct event *ev); -- cgit v1.2.3 From 9b168c16fd09407b2713af3dd72d983f16eec08e Mon Sep 17 00:00:00 2001 From: Thomas White Date: Tue, 2 May 2017 11:58:35 +0200 Subject: HDF5 reading under new API --- libcrystfel/src/image.c | 83 ++++++++++++++++++++++++++++++++++++++++++------- 1 file changed, 71 insertions(+), 12 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/image.c b/libcrystfel/src/image.c index f322cbdd..5aa2f9c3 100644 --- a/libcrystfel/src/image.c +++ b/libcrystfel/src/image.c @@ -31,9 +31,12 @@ #include #include #include +#include #include "image.h" #include "utils.h" +#include "events.h" +#include "hdf5-file.h" /** * SECTION:image @@ -378,19 +381,69 @@ void add_imagefile_field(struct imagefile_field_list *copyme, const char *name) /****************************** Image files ***********************************/ +enum imagefile_type +{ + IMAGEFILE_HDF5, + IMAGEFILE_CBF +}; + + +struct imagefile +{ + enum imagefile_type type; + struct hdfile *hdfile; +}; + + struct imagefile *imagefile_open(const char *filename) { + struct imagefile *f; + + f = malloc(sizeof(struct imagefile)); + if ( f == NULL ) return NULL; + + if ( H5Fis_hdf5(filename) > 0 ) { + + /* This is an HDF5, pass through to HDF5 layer */ + f->type = IMAGEFILE_HDF5; + f->hdfile = hdfile_open(filename); + + if ( f->hdfile == NULL ) { + free(f); + return NULL; + } + + } else { + + STATUS("Mock CBF check\n"); + return NULL; + + } + + return f; } -int imagefile_read(struct imagefile *imfile, struct image *image, - struct event *event) +int imagefile_read(struct imagefile *f, struct image *image, + struct event *event) { + if ( f->type == IMAGEFILE_HDF5 ) { + return hdf5_read2(f->hdfile, image, event, 0); + } else { + STATUS("Mock CBF read\n"); + return 0; + } } -struct hdfile *imagefile_get_hdfile(struct imagefile *imfile) +struct hdfile *imagefile_get_hdfile(struct imagefile *f) { + if ( f->type != IMAGEFILE_HDF5 ) { + ERROR("Not an HDF5 file!\n"); + return NULL; + } + + return f->hdfile; } @@ -404,27 +457,33 @@ void imagefile_copy_fields(struct imagefile *f, for ( i=0; in_fields; i++ ) { -#warning FIXME Implement imagefile field copying -#if 0 char *val; char *field; field = copyme->fields[i]; - val = hdfile_get_string_value(f, field, ev); - if ( field[0] == '/' ) { - fprintf(fh, "hdf5%s = %s\n", field, val); + if ( f->type == IMAGEFILE_HDF5 ) { + val = hdfile_get_string_value(f->hdfile, field, ev); + if ( field[0] == '/' ) { + fprintf(fh, "hdf5%s = %s\n", field, val); + } else { + fprintf(fh, "hdf5/%s = %s\n", field, val); + } + free(val); + } else { - fprintf(fh, "hdf5/%s = %s\n", field, val); + STATUS("Mock CBF variable\n"); + fprintf(fh, "cbf/%s = %s\n", field, "(FIXME)"); } - free(val); -#endif } } -void imagefile_close(struct imagefile *imfile) +void imagefile_close(struct imagefile *f) { + if ( f->type == IMAGEFILE_HDF5 ) { + hdfile_close(f->hdfile); + } } -- cgit v1.2.3 From 00f8d418534b30e71fd925487c0284560e90b2b2 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 3 May 2017 11:19:26 +0200 Subject: hdfile_set_image(): Remove panel argument Seems to have been added in 2014 and is not used at the moment --- libcrystfel/src/hdf5-file.c | 11 +++++------ libcrystfel/src/hdf5-file.h | 4 ++-- 2 files changed, 7 insertions(+), 8 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/hdf5-file.c b/libcrystfel/src/hdf5-file.c index c36a58b3..9a52fd4d 100644 --- a/libcrystfel/src/hdf5-file.c +++ b/libcrystfel/src/hdf5-file.c @@ -134,8 +134,7 @@ struct hdfile *hdfile_open(const char *filename) } -int hdfile_set_image(struct hdfile *f, const char *path, - struct panel *p) +int hdfile_set_image(struct hdfile *f, const char *path) { f->dh = H5Dopen2(f->fh, path, H5P_DEFAULT); if ( f->dh < 0 ) { @@ -1433,7 +1432,7 @@ int hdf5_read(struct hdfile *f, struct image *image, const char *element, if ( element == NULL ) { fail = hdfile_set_first_image(f, "/"); } else { - fail = hdfile_set_image(f, element, NULL); + fail = hdfile_set_image(f, element); } if ( fail ) { @@ -1747,7 +1746,7 @@ int hdf5_read2(struct hdfile *f, struct image *image, struct event *ev, return 1; } - fail = hdfile_set_image(f, panel_full_path, p); + fail = hdfile_set_image(f, panel_full_path); free(panel_full_path); @@ -1769,7 +1768,7 @@ int hdf5_read2(struct hdfile *f, struct image *image, struct event *ev, free(image->sat); return 1; } - fail = hdfile_set_image(f, p->data, p); + fail = hdfile_set_image(f, p->data); } @@ -2299,7 +2298,7 @@ int hdfile_set_first_image(struct hdfile *f, const char *group) for ( i=0; i Date: Wed, 3 May 2017 11:25:24 +0200 Subject: Detect CBF files, interface bits --- libcrystfel/src/image.c | 70 ++++++++++++++++++++++++++++++++++++++++++------- libcrystfel/src/image.h | 19 +++++++++++--- 2 files changed, 76 insertions(+), 13 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/image.c b/libcrystfel/src/image.c index 5aa2f9c3..29e5d5e1 100644 --- a/libcrystfel/src/image.c +++ b/libcrystfel/src/image.c @@ -32,6 +32,7 @@ #include #include #include +#include #include "image.h" #include "utils.h" @@ -379,14 +380,25 @@ void add_imagefile_field(struct imagefile_field_list *copyme, const char *name) } -/****************************** Image files ***********************************/ +/******************************* CBF files ************************************/ -enum imagefile_type +static int read_cbf(struct imagefile *f, struct image *image) { - IMAGEFILE_HDF5, - IMAGEFILE_CBF -}; + cbf_handle cbfh; + + if ( cbf_make_handle(&cbfh) ) { + ERROR("Failed to allocate CBF handle\n"); + return 1; + } + + ERROR("Mock CBF read\n"); + + cbf_free_handle(cbfh); + return 0; +} + +/****************************** Image files ***********************************/ struct imagefile { @@ -395,6 +407,25 @@ struct imagefile }; +static signed int is_cbf_file(const char *filename) +{ + FILE *fh; + char line[1024]; + + fh = fopen(filename, "r"); + if ( fh == NULL ) return -1; + + if ( fgets(line, 1024, fh) == NULL ) return -1; + fclose(fh); + + if ( strstr(line, "CBF") == NULL ) { + return 0; + } + + return 1; +} + + struct imagefile *imagefile_open(const char *filename) { struct imagefile *f; @@ -413,10 +444,9 @@ struct imagefile *imagefile_open(const char *filename) return NULL; } - } else { + } else if ( is_cbf_file(filename) > 0 ) { - STATUS("Mock CBF check\n"); - return NULL; + f->type = IMAGEFILE_CBF; } @@ -429,13 +459,35 @@ int imagefile_read(struct imagefile *f, struct image *image, { if ( f->type == IMAGEFILE_HDF5 ) { return hdf5_read2(f->hdfile, image, event, 0); + } else if ( f->type == IMAGEFILE_CBF ) { + return read_cbf(f, image); + } else { + ERROR("Unknown file type %i\n", f->type); + return 1; + } +} + + +/* Read a simple file, no multi-event, no prior geometry etc, and + * generate a geometry for it */ +int imagefile_read_simple(struct imagefile *f, struct image *image) +{ + if ( f->type == IMAGEFILE_HDF5 ) { + return hdf5_read(f->hdfile, image, NULL, 0); } else { - STATUS("Mock CBF read\n"); + STATUS("Mock CBF simple read\n"); return 0; } } +enum imagefile_type imagefile_get_type(struct imagefile *f) +{ + assert(f != NULL); + return f->type; +} + + struct hdfile *imagefile_get_hdfile(struct imagefile *f) { if ( f->type != IMAGEFILE_HDF5 ) { diff --git a/libcrystfel/src/image.h b/libcrystfel/src/image.h index f310a09b..ead5cd4d 100644 --- a/libcrystfel/src/image.h +++ b/libcrystfel/src/image.h @@ -90,6 +90,15 @@ struct imagefeature { const char *name; }; + +/* An enum representing the image file formats we can handle */ +enum imagefile_type +{ + IMAGEFILE_HDF5, + IMAGEFILE_CBF +}; + + /* An opaque type representing a list of image features */ typedef struct _imagefeaturelist ImageFeatureList; @@ -238,13 +247,15 @@ extern void free_all_crystals(struct image *image); /* Image files */ extern struct imagefile *imagefile_open(const char *filename); -extern int imagefile_read(struct imagefile *imfile, struct image *image, +extern int imagefile_read(struct imagefile *f, struct image *image, struct event *event); -extern struct hdfile *imagefile_get_hdfile(struct imagefile *imfile); -extern void imagefile_copy_fields(struct imagefile *imfile, +extern int imagefile_read_simple(struct imagefile *f, struct image *image); +extern struct hdfile *imagefile_get_hdfile(struct imagefile *f); +extern enum imagefile_type imagefile_get_type(struct imagefile *f); +extern void imagefile_copy_fields(struct imagefile *f, const struct imagefile_field_list *copyme, FILE *fh, struct event *ev); -extern void imagefile_close(struct imagefile *imfile); +extern void imagefile_close(struct imagefile *f); /* Field lists */ extern struct imagefile_field_list *new_imagefile_field_list(void); -- cgit v1.2.3 From 4c5ae5576d943350fd958db3674f4d30ae2d6eae Mon Sep 17 00:00:00 2001 From: Thomas White Date: Thu, 4 May 2017 17:17:27 +0200 Subject: Read CBF data --- libcrystfel/src/image.c | 235 ++++++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 228 insertions(+), 7 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/image.c b/libcrystfel/src/image.c index 29e5d5e1..a8c278fc 100644 --- a/libcrystfel/src/image.c +++ b/libcrystfel/src/image.c @@ -38,6 +38,7 @@ #include "utils.h" #include "events.h" #include "hdf5-file.h" +#include "detector.h" /** * SECTION:image @@ -54,6 +55,14 @@ */ +struct imagefile +{ + enum imagefile_type type; + char *filename; + struct hdfile *hdfile; +}; + + struct _imagefeaturelist { struct imagefeature *features; @@ -382,16 +391,231 @@ void add_imagefile_field(struct imagefile_field_list *copyme, const char *name) /******************************* CBF files ************************************/ + +static char *cbf_strerr(int e) +{ + char *err; + + err = malloc(1024); + if ( err == NULL ) return NULL; + + err[0] = '\0'; + + /* NB Sum of lengths of all strings must be less than 1024 */ + if ( e & CBF_FORMAT ) strcat(err, "Invalid format"); + if ( e & CBF_ALLOC ) strcat(err, "Memory allocation failed"); + if ( e & CBF_ARGUMENT ) strcat(err, "Invalid argument"); + if ( e & CBF_ASCII ) strcat(err, "Value is ASCII"); + if ( e & CBF_BINARY ) strcat(err, "Value is binary"); + if ( e & CBF_BITCOUNT ) strcat(err, "Wrong number of bits"); + if ( e & CBF_ENDOFDATA ) strcat(err, "End of data"); + if ( e & CBF_FILECLOSE ) strcat(err, "File close error"); + if ( e & CBF_FILEOPEN ) strcat(err, "File open error"); + if ( e & CBF_FILEREAD ) strcat(err, "File read error"); + if ( e & CBF_FILETELL ) strcat(err, "File tell error"); + if ( e & CBF_FILEWRITE ) strcat(err, "File write error"); + if ( e & CBF_IDENTICAL ) strcat(err, "Name already exists"); + if ( e & CBF_NOTFOUND ) strcat(err, "Not found"); + if ( e & CBF_OVERFLOW ) strcat(err, "Overflow"); + if ( e & CBF_UNDEFINED ) strcat(err, "Number undefined"); + if ( e & CBF_NOTIMPLEMENTED ) strcat(err, "Not implemented"); + + return err; +} + + +static int unpack_panels(struct image *image, signed int *data, int data_width) +{ + int pi; + + /* FIXME: Load these masks from an HDF5 file, if filenames are + * given in the geometry file */ + uint16_t *flags = NULL; + float *sat = NULL; + + image->dp = malloc(image->det->n_panels * sizeof(float *)); + image->bad = malloc(image->det->n_panels * sizeof(int *)); + image->sat = malloc(image->det->n_panels * sizeof(float *)); + if ( (image->dp == NULL) || (image->bad == NULL) + || (image->sat == NULL) ) + { + ERROR("Failed to allocate panels.\n"); + return 1; + } + + for ( pi=0; pidet->n_panels; pi++ ) { + + struct panel *p; + int fs, ss; + + p = &image->det->panels[pi]; + image->dp[pi] = malloc(p->w*p->h*sizeof(float)); + image->bad[pi] = calloc(p->w*p->h, sizeof(int)); + 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; + } + + if ( p->mask != NULL ) { + ERROR("WARNING: Bad pixel masks do not currently work " + "with CBF files\n"); + ERROR(" (bad pixel regions specified in the geometry " + "file will be used, however)\n"); + } + + if ( p->satmap != NULL ) { + ERROR("WARNING: Saturation maps do not currently work " + "with CBF files\n"); + } + + for ( ss=0; ssh; ss++ ) { + for ( fs=0; fsw; fs++ ) { + + int idx; + int cfs, css; + int bad = 0; + + cfs = fs+p->orig_min_fs; + css = ss+p->orig_min_ss; + idx = cfs + css*data_width; + + image->dp[pi][fs+p->w*ss] = data[idx]; + + if ( sat != NULL ) { + image->sat[pi][fs+p->w*ss] = sat[idx]; + } else { + image->sat[pi][fs+p->w*ss] = INFINITY; + } + + if ( p->no_index ) bad = 1; + + if ( in_bad_region(image->det, p, cfs, css) ) { + bad = 1; + } + + if ( flags != NULL ) { + + int f; + + f = flags[idx]; + + /* Bad if it's missing any of the "good" bits */ + if ( (f & image->det->mask_good) + != image->det->mask_good ) bad = 1; + + /* Bad if it has any of the "bad" bits. */ + if ( f & image->det->mask_bad ) bad = 1; + + } + image->bad[pi][fs+p->w*ss] = bad; + + } + } + + } + + return 0; +} + + static int read_cbf(struct imagefile *f, struct image *image) { cbf_handle cbfh; + FILE *fh; + int r; + unsigned int compression; + int binary_id, minelement, maxelement, elsigned, elunsigned; + size_t elsize, elements, elread, dimfast, dimmid, dimslow, padding; + const char *byteorder; + signed int *data; if ( cbf_make_handle(&cbfh) ) { ERROR("Failed to allocate CBF handle\n"); return 1; } - ERROR("Mock CBF read\n"); + fh = fopen(f->filename, "rb"); + if ( fh == NULL ) { + ERROR("Failed to open '%s'\n", f->filename); + return 1; + } + /* CBFlib calls fclose(fh) when it's ready */ + + if ( cbf_read_widefile(cbfh, fh, 0) ) { + ERROR("Failed to read CBF file\n"); + cbf_free_handle(cbfh); + return 1; + } + + /* Select row 0 in data column inside array_data */ + cbf_find_category(cbfh, "array_data"); + cbf_find_column(cbfh, "data"); + cbf_select_row(cbfh, 0); + + /* Get parameters for array read */ + r = cbf_get_integerarrayparameters_wdims(cbfh, &compression, &binary_id, + &elsize, &elsigned, &elunsigned, + &elements, + &minelement, &maxelement, + &byteorder, + &dimfast, &dimmid, &dimslow, + &padding); + if ( r ) { + char *err = cbf_strerr(r); + ERROR("Failed to read CBF array parameters: %s\n", err); + free(err); + cbf_free_handle(cbfh); + return 1; + } + + if ( dimslow != 0 ) { + ERROR("CBF data array is 3D - don't know what to do with it\n"); + cbf_free_handle(cbfh); + return 1; + } + + if ( dimfast*dimmid*elsize > 10e9 ) { + ERROR("CBF data is far too big (%i x %i x %i bytes).\n", + (int)dimfast, (int)dimmid, (int)elsize); + cbf_free_handle(cbfh); + return 1; + } + + if ( elsize != 4 ) { + STATUS("Don't know what to do with element size %i\n", + (int)elsize); + cbf_free_handle(cbfh); + return 1; + } + + if ( strcmp(byteorder, "little_endian") != 0 ) { + STATUS("Don't know what to do with non-little-endian datan\n"); + cbf_free_handle(cbfh); + return 1; + } + + data = malloc(elsize*dimfast*dimmid); + if ( data == NULL ) { + ERROR("Failed to allocate memory for CBF data\n"); + cbf_free_handle(cbfh); + return 1; + } + + r = cbf_get_integerarray(cbfh, &binary_id, data, elsize, 1, + elsize*dimfast*dimmid, &elread); + if ( r ) { + char *err = cbf_strerr(r); + ERROR("Failed to read CBF array: %s\n", err); + free(err); + cbf_free_handle(cbfh); + return 1; + } + + unpack_panels(image, data, dimfast); + free(data); cbf_free_handle(cbfh); return 0; @@ -400,12 +624,6 @@ static int read_cbf(struct imagefile *f, struct image *image) /****************************** Image files ***********************************/ -struct imagefile -{ - enum imagefile_type type; - struct hdfile *hdfile; -}; - static signed int is_cbf_file(const char *filename) { @@ -450,6 +668,7 @@ struct imagefile *imagefile_open(const char *filename) } + f->filename = strdup(filename); return f; } @@ -538,4 +757,6 @@ void imagefile_close(struct imagefile *f) if ( f->type == IMAGEFILE_HDF5 ) { hdfile_close(f->hdfile); } + free(f->filename); + free(f); } -- cgit v1.2.3 From bb552aba7ac926bfae331b5648abd8647ba1d9f5 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 5 May 2017 16:33:38 +0200 Subject: Allow location of photon energy (eg in HDF5 file) to start with any character --- libcrystfel/src/detector.c | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/detector.c b/libcrystfel/src/detector.c index e160d1e1..b2e3cc58 100644 --- a/libcrystfel/src/detector.c +++ b/libcrystfel/src/detector.c @@ -1076,12 +1076,17 @@ static void parse_toplevel(struct detector *det, struct beam_params *beam, } else if ( strcmp(key, "photon_energy") == 0 ) { if ( beam != NULL ) { - if ( strncmp(val, "/", 1) == 0 ) { + + double v; + char *end; + + v = strtod(val, &end); + if ( (val[0] != '\0') && (end[0] == '\0') ) { + beam->photon_energy = v; + beam->photon_energy_from = NULL; + } else { beam->photon_energy = 0.0; beam->photon_energy_from = strdup(val); - } else { - beam->photon_energy = atof(val); - beam->photon_energy_from = NULL; } } -- cgit v1.2.3 From 9104fb61c53f7d7d8fcb32627a7f87c2062ef496 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 5 May 2017 13:42:54 +0200 Subject: Fill in photon energy, clen and adu for CBFs This needed a bit of reorganisation and clarification of who is repsonsible for loading what. The low-level file loaders, e.g. hdf5_read and hdf5_read2 in hdf5-file.c or cbf_read in image.c, are responsible. There is a helper function adjust_centering_for_rail in detector.h which they can use. It seems like this could be done more cleanly at the imagefile layer. However, we need to make sure that the "old" hdfile API still works on its own, even when not accessed via the new imagefile API. --- libcrystfel/src/detector.c | 38 ++++++-------------------------- libcrystfel/src/detector.h | 3 +-- libcrystfel/src/hdf5-file.c | 46 +++++++++++++++++++++++++++++++++------ libcrystfel/src/image.c | 53 +++++++++++++++++++++++++++++++++++++++++++++ 4 files changed, 100 insertions(+), 40 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/detector.c b/libcrystfel/src/detector.c index b2e3cc58..e3b351ff 100644 --- a/libcrystfel/src/detector.c +++ b/libcrystfel/src/detector.c @@ -528,37 +528,15 @@ int panel_number(struct detector *det, struct panel *p) } -void fill_in_values(struct detector *det, struct hdfile *f, struct event* ev) +void adjust_centering_for_rail(struct panel *p) { - int i; - - for ( i=0; in_panels; i++ ) { - - double offs; - struct panel *p = &det->panels[i]; + double offs; - if ( p->clen_from != NULL ) { - - double val; - int r; - - r = hdfile_get_value(f, p->clen_from, ev, &val, - H5T_NATIVE_DOUBLE); - if ( r ) { - ERROR("Failed to read '%s'\n", p->clen_from); - } else { - p->clen = val * 1.0e-3; - } - - } - - /* Offset in +z direction from calibrated clen to actual */ - offs = p->clen - p->clen_for_centering; - p->cnx += p->rail_x * offs; - p->cny += p->rail_y * offs; - p->clen = p->clen_for_centering + p->coffset + p->rail_z * offs; - - } + /* Offset in +z direction from calibrated clen to actual */ + offs = p->clen - p->clen_for_centering; + p->cnx += p->rail_x * offs; + p->cny += p->rail_y * offs; + p->clen = p->clen_for_centering + p->coffset + p->rail_z * offs; } @@ -1076,10 +1054,8 @@ static void parse_toplevel(struct detector *det, struct beam_params *beam, } else if ( strcmp(key, "photon_energy") == 0 ) { if ( beam != NULL ) { - double v; char *end; - v = strtod(val, &end); if ( (val[0] != '\0') && (end[0] == '\0') ) { beam->photon_energy = v; diff --git a/libcrystfel/src/detector.h b/libcrystfel/src/detector.h index 04e3c1ec..e9dd1154 100644 --- a/libcrystfel/src/detector.h +++ b/libcrystfel/src/detector.h @@ -224,9 +224,8 @@ extern void get_pixel_extents(struct detector *det, double *min_x, double *min_y, double *max_x, double *max_y); -extern void fill_in_values(struct detector *det, struct hdfile *f, - struct event* ev); extern void fill_in_adu(struct image *image); +extern void adjust_centering_for_rail(struct panel *p); extern int panel_is_in_rigid_group(const struct rigid_group *rg, struct panel *p); diff --git a/libcrystfel/src/hdf5-file.c b/libcrystfel/src/hdf5-file.c index 9a52fd4d..dc86d641 100644 --- a/libcrystfel/src/hdf5-file.c +++ b/libcrystfel/src/hdf5-file.c @@ -1390,8 +1390,10 @@ int hdfile_get_value(struct hdfile *f, const char *name, } -void fill_in_beam_parameters(struct beam_params *beam, struct hdfile *f, - struct event *ev, struct image *image) +static void hdfile_fill_in_beam_parameters(struct beam_params *beam, + struct hdfile *f, + struct event *ev, + struct image *image) { double eV; @@ -1404,8 +1406,8 @@ void fill_in_beam_parameters(struct beam_params *beam, struct hdfile *f, int r; - r = hdfile_get_value(f, beam->photon_energy_from, ev, &eV, - H5T_NATIVE_DOUBLE); + r = hdfile_get_value(f, beam->photon_energy_from, + ev, &eV, H5T_NATIVE_DOUBLE); if ( r ) { ERROR("Failed to read '%s'\n", beam->photon_energy_from); @@ -1417,6 +1419,36 @@ void fill_in_beam_parameters(struct beam_params *beam, struct hdfile *f, } +static void hdfile_fill_in_clen(struct detector *det, struct hdfile *f, + struct event *ev) +{ + int i; + + for ( i=0; in_panels; i++ ) { + + struct panel *p = &det->panels[i]; + + if ( p->clen_from != NULL ) { + + double val; + int r; + + r = hdfile_get_value(f, p->clen_from, ev, &val, + H5T_NATIVE_DOUBLE); + if ( r ) { + ERROR("Failed to read '%s'\n", p->clen_from); + } else { + p->clen = val * 1.0e-3; + } + + } + + adjust_centering_for_rail(p); + + } +} + + int hdf5_read(struct hdfile *f, struct image *image, const char *element, int satcorr) { @@ -1481,7 +1513,7 @@ int hdf5_read(struct hdfile *f, struct image *image, const char *element, if ( image->beam != NULL ) { - fill_in_beam_parameters(image->beam, f, NULL, image); + hdfile_fill_in_beam_parameters(image->beam, f, NULL, image); if ( image->lambda > 1000 ) { /* Error message covers a silly value in the beam file @@ -1885,13 +1917,13 @@ int hdf5_read2(struct hdfile *f, struct image *image, struct event *ev, H5Dclose(f->dh); f->data_open = 0; - fill_in_values(image->det, f, ev); + hdfile_fill_in_clen(image->det, f, ev); if ( satcorr ) debodge_saturation(f, image); if ( image->beam != NULL ) { - fill_in_beam_parameters(image->beam, f, ev, image); + hdfile_fill_in_beam_parameters(image->beam, f, ev, image); if ( (image->lambda > 1.0) || (image->lambda < 1e-20) ) { diff --git a/libcrystfel/src/image.c b/libcrystfel/src/image.c index a8c278fc..849d95bf 100644 --- a/libcrystfel/src/image.c +++ b/libcrystfel/src/image.c @@ -521,6 +521,48 @@ static int unpack_panels(struct image *image, signed int *data, int data_width) } +static void cbf_fill_in_beam_parameters(struct beam_params *beam, + struct imagefile *f, + struct image *image) +{ + double eV; + + if ( beam->photon_energy_from == NULL ) { + + /* Explicit value given */ + eV = beam->photon_energy; + + } else { + + ERROR("Can't get photon energy from CBF yet.\n"); + eV = 0.0; + + } + + image->lambda = ph_en_to_lambda(eV_to_J(eV))*beam->photon_energy_scale; +} + + +static void cbf_fill_in_clen(struct detector *det, struct imagefile *f) +{ + int i; + + for ( i=0; in_panels; i++ ) { + + struct panel *p = &det->panels[i]; + + if ( p->clen_from != NULL ) { + + ERROR("Can't get clen from CBF yet.\n"); + + } + + adjust_centering_for_rail(p); + + } +} + + static int read_cbf(struct imagefile *f, struct image *image) { cbf_handle cbfh; @@ -617,6 +659,17 @@ static int read_cbf(struct imagefile *f, struct image *image) unpack_panels(image, data, dimfast); free(data); + if ( image->beam != NULL ) { + cbf_fill_in_beam_parameters(image->beam, f, image); + if ( image->lambda > 1000 ) { + ERROR("WARNING: Missing or nonsensical wavelength " + "(%e m) for %s.\n", + image->lambda, image->filename); + } + } + cbf_fill_in_clen(image->det, f); + fill_in_adu(image); + cbf_free_handle(cbfh); return 0; } -- cgit v1.2.3 From cdd75c7fadfcd2dadee101611a00dcf0d0e6ec5c Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 5 May 2017 17:06:12 +0200 Subject: Show filename of CBF file if it can't be read --- libcrystfel/src/image.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/image.c b/libcrystfel/src/image.c index 849d95bf..c03b60d8 100644 --- a/libcrystfel/src/image.c +++ b/libcrystfel/src/image.c @@ -587,7 +587,7 @@ static int read_cbf(struct imagefile *f, struct image *image) /* CBFlib calls fclose(fh) when it's ready */ if ( cbf_read_widefile(cbfh, fh, 0) ) { - ERROR("Failed to read CBF file\n"); + ERROR("Failed to read CBF file '%s'\n", f->filename); cbf_free_handle(cbfh); return 1; } -- cgit v1.2.3 From 66dc7ec0acff10d6f22cb9422bb430018d31722e Mon Sep 17 00:00:00 2001 From: Thomas White Date: Tue, 9 May 2017 17:03:16 +0200 Subject: Add missing checks --- libcrystfel/src/image.c | 18 ++++++++++++++++-- 1 file changed, 16 insertions(+), 2 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/image.c b/libcrystfel/src/image.c index c03b60d8..ef04cbf9 100644 --- a/libcrystfel/src/image.c +++ b/libcrystfel/src/image.c @@ -424,7 +424,8 @@ static char *cbf_strerr(int e) } -static int unpack_panels(struct image *image, signed int *data, int data_width) +static int unpack_panels(struct image *image, signed int *data, int data_width, + int data_height) { int pi; @@ -471,6 +472,14 @@ static int unpack_panels(struct image *image, signed int *data, int data_width) "with CBF files\n"); } + if ( (p->orig_min_fs + p->w > data_width) + || (p->orig_min_ss + p->h > data_height) ) + { + ERROR("Panel %s is outside range of data in CBF file\n", + p->name); + return 1; + } + for ( ss=0; ssh; ss++ ) { for ( fs=0; fsw; fs++ ) { @@ -574,6 +583,11 @@ static int read_cbf(struct imagefile *f, struct image *image) const char *byteorder; signed int *data; + if ( image->det == NULL ) { + ERROR("read_cbf() needs a geometry\n"); + return 1; + } + if ( cbf_make_handle(&cbfh) ) { ERROR("Failed to allocate CBF handle\n"); return 1; @@ -656,7 +670,7 @@ static int read_cbf(struct imagefile *f, struct image *image) return 1; } - unpack_panels(image, data, dimfast); + unpack_panels(image, data, dimfast, dimmid); free(data); if ( image->beam != NULL ) { -- cgit v1.2.3 From 5617c8273e26469689c66fb46d5f7eb34dbf2bf4 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Tue, 9 May 2017 17:03:36 +0200 Subject: Allow reading CBF without geometry (only for simple viewing) --- libcrystfel/src/image.c | 140 +++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 138 insertions(+), 2 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/image.c b/libcrystfel/src/image.c index ef04cbf9..4e1a8b86 100644 --- a/libcrystfel/src/image.c +++ b/libcrystfel/src/image.c @@ -689,6 +689,143 @@ static int read_cbf(struct imagefile *f, struct image *image) } +static float *convert_float(signed int *data, int w, int h) +{ + float *df; + long int i; + + df = malloc(sizeof(float)*w*h); + if ( df == NULL ) return NULL; + + for ( i=0; ifilename, "rb"); + if ( fh == NULL ) { + ERROR("Failed to open '%s'\n", f->filename); + return 1; + } + /* CBFlib calls fclose(fh) when it's ready */ + + if ( cbf_read_widefile(cbfh, fh, 0) ) { + ERROR("Failed to read CBF file '%s'\n", f->filename); + cbf_free_handle(cbfh); + return 1; + } + + /* Select row 0 in data column inside array_data */ + cbf_find_category(cbfh, "array_data"); + cbf_find_column(cbfh, "data"); + cbf_select_row(cbfh, 0); + + /* Get parameters for array read */ + r = cbf_get_integerarrayparameters_wdims(cbfh, &compression, &binary_id, + &elsize, &elsigned, &elunsigned, + &elements, + &minelement, &maxelement, + &byteorder, + &dimfast, &dimmid, &dimslow, + &padding); + if ( r ) { + char *err = cbf_strerr(r); + ERROR("Failed to read CBF array parameters: %s\n", err); + free(err); + cbf_free_handle(cbfh); + return 1; + } + + if ( dimslow != 0 ) { + ERROR("CBF data array is 3D - don't know what to do with it\n"); + cbf_free_handle(cbfh); + return 1; + } + + if ( dimfast*dimmid*elsize > 10e9 ) { + ERROR("CBF data is far too big (%i x %i x %i bytes).\n", + (int)dimfast, (int)dimmid, (int)elsize); + cbf_free_handle(cbfh); + return 1; + } + + if ( elsize != 4 ) { + STATUS("Don't know what to do with element size %i\n", + (int)elsize); + cbf_free_handle(cbfh); + return 1; + } + + if ( strcmp(byteorder, "little_endian") != 0 ) { + STATUS("Don't know what to do with non-little-endian datan\n"); + cbf_free_handle(cbfh); + return 1; + } + + data = malloc(elsize*dimfast*dimmid); + if ( data == NULL ) { + ERROR("Failed to allocate memory for CBF data\n"); + cbf_free_handle(cbfh); + return 1; + } + + r = cbf_get_integerarray(cbfh, &binary_id, data, elsize, 1, + elsize*dimfast*dimmid, &elread); + if ( r ) { + char *err = cbf_strerr(r); + ERROR("Failed to read CBF array: %s\n", err); + free(err); + cbf_free_handle(cbfh); + return 1; + } + + image->det = simple_geometry(image, dimfast, dimmid); + image->dp = malloc(sizeof(float *)); + if ( image->dp == NULL ) { + ERROR("Failed to allocate dp array\n"); + return 1; + } + image->dp[0] = convert_float(data, dimfast, dimmid); + if ( image->dp[0] == NULL ) { + ERROR("Failed to allocate dp array\n"); + return 1; + } + + if ( image->beam != NULL ) { + cbf_fill_in_beam_parameters(image->beam, f, image); + if ( image->lambda > 1000 ) { + ERROR("WARNING: Missing or nonsensical wavelength " + "(%e m) for %s.\n", + image->lambda, image->filename); + } + } + cbf_fill_in_clen(image->det, f); + fill_in_adu(image); + + cbf_free_handle(cbfh); + return 0; +} + + /****************************** Image files ***********************************/ @@ -761,8 +898,7 @@ int imagefile_read_simple(struct imagefile *f, struct image *image) if ( f->type == IMAGEFILE_HDF5 ) { return hdf5_read(f->hdfile, image, NULL, 0); } else { - STATUS("Mock CBF simple read\n"); - return 0; + return read_cbf_simple(f, image); } } -- cgit v1.2.3 From 054fefd62406a76711cb8ddaa5d2e5e9271570a0 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 26 May 2017 10:33:47 +0200 Subject: hdf5_read2(): Fix exit path Bad pixel map is not allocated at this point. --- libcrystfel/src/hdf5-file.c | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/hdf5-file.c b/libcrystfel/src/hdf5-file.c index dc86d641..77a69dbd 100644 --- a/libcrystfel/src/hdf5-file.c +++ b/libcrystfel/src/hdf5-file.c @@ -1881,9 +1881,8 @@ int hdf5_read2(struct hdfile *f, struct image *image, struct event *ev, free(f_offset); free(f_count); for ( i=0; i<=pi; i++ ) { - free(image->dp[pi]); - free(image->sat[pi]); - free(image->bad[pi]); + free(image->dp[i]); + free(image->sat[i]); } free(image->dp); free(image->bad); -- cgit v1.2.3 From af43cdaa8e9d53b71f7ca0dfc0f16040c96fef7c Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 26 May 2017 16:41:47 +0200 Subject: Close leftover HDF5 handles --- libcrystfel/src/hdf5-file.c | 3 +++ 1 file changed, 3 insertions(+) (limited to 'libcrystfel') diff --git a/libcrystfel/src/hdf5-file.c b/libcrystfel/src/hdf5-file.c index 77a69dbd..eafc624d 100644 --- a/libcrystfel/src/hdf5-file.c +++ b/libcrystfel/src/hdf5-file.c @@ -2647,6 +2647,9 @@ static int check_dims(struct hdfile *hdfile, struct panel *p, struct event *ev, return 1; } + H5Sclose(sh); + H5Dclose(dh); + return 0; } -- cgit v1.2.3 From 6e10e48b217ea7abfc0e9edc3ef6211f676e5c31 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 21 Jun 2017 14:34:19 +0200 Subject: Don't compile CBF stuff without CBFlib --- libcrystfel/src/image.c | 22 ++++++++++++++++++++++ 1 file changed, 22 insertions(+) (limited to 'libcrystfel') diff --git a/libcrystfel/src/image.c b/libcrystfel/src/image.c index 4e1a8b86..1baa33bc 100644 --- a/libcrystfel/src/image.c +++ b/libcrystfel/src/image.c @@ -32,7 +32,10 @@ #include #include #include + +#ifdef HAVE_CBFLIB #include +#endif #include "image.h" #include "utils.h" @@ -391,6 +394,7 @@ void add_imagefile_field(struct imagefile_field_list *copyme, const char *name) /******************************* CBF files ************************************/ +#ifdef HAVE_CBFLIB static char *cbf_strerr(int e) { @@ -825,6 +829,24 @@ static int read_cbf_simple(struct imagefile *f, struct image *image) return 0; } +#else /* HAVE_CBFLIB */ + +static int read_cbf_simple(struct imagefile *f, struct image *image) +{ + ERROR("This version of CrystFEL was compiled without CBF support.\n"); + return 1; +} + + +static int read_cbf(struct imagefile *f, struct image *image) +{ + ERROR("This version of CrystFEL was compiled without CBF support.\n"); + return 1; +} + + +#endif /* HAVE_CBFLIB */ + /****************************** Image files ***********************************/ -- cgit v1.2.3 From 3c80790a4ddd12f3d1d864d12a0126acf24ed47e Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 21 Jun 2017 15:04:01 +0200 Subject: Close more leftover HDF5 handles --- libcrystfel/src/hdf5-file.c | 3 +++ 1 file changed, 3 insertions(+) (limited to 'libcrystfel') diff --git a/libcrystfel/src/hdf5-file.c b/libcrystfel/src/hdf5-file.c index eafc624d..13339506 100644 --- a/libcrystfel/src/hdf5-file.c +++ b/libcrystfel/src/hdf5-file.c @@ -1233,6 +1233,9 @@ static int get_scalar_value(struct hdfile *f, const char *name, void *val, return 1; } + H5Tclose(type); + H5Dclose(dh); + return 0; } -- cgit v1.2.3 From 80c6061730c2a6a02b9eea77140026ff6b822856 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 23 Jun 2017 11:05:18 +0200 Subject: Accept either cbf/cbf.h or cbflib/cbf.h You get cbflib/cbf.h if you install CBFlib from source, and cbf/cbf.h if using the package manager in Fedora, probably other systems as well. We definitely don't want to include the cbf or cbflib folder in the include path, because CBFlib may bring its own HDF5 headers, which we don't want to touch with a bargepole. --- libcrystfel/src/image.c | 7 +++++++ 1 file changed, 7 insertions(+) (limited to 'libcrystfel') diff --git a/libcrystfel/src/image.c b/libcrystfel/src/image.c index 1baa33bc..aad5017c 100644 --- a/libcrystfel/src/image.c +++ b/libcrystfel/src/image.c @@ -27,6 +27,8 @@ * */ +#include + #include #include #include @@ -34,8 +36,13 @@ #include #ifdef HAVE_CBFLIB +#ifdef HAVE_CBFLIB_CBF_H #include #endif +#ifdef HAVE_CBF_CBF_H +#include +#endif +#endif #include "image.h" #include "utils.h" -- cgit v1.2.3 From 61d8d5cf74fb962a25742c7547b7450414ae5a5b Mon Sep 17 00:00:00 2001 From: Thomas White Date: Tue, 27 Jun 2017 10:28:45 +0200 Subject: If prediction refinement fails, don't carry on and check the cell Doing so results in one crystal being counted as bad twice, which messes up the logic which follows. --- libcrystfel/src/index.c | 1 + 1 file changed, 1 insertion(+) (limited to 'libcrystfel') diff --git a/libcrystfel/src/index.c b/libcrystfel/src/index.c index 685d0b69..9bcf5e15 100644 --- a/libcrystfel/src/index.c +++ b/libcrystfel/src/index.c @@ -345,6 +345,7 @@ static int try_indexer(struct image *image, IndexingMethod indm, if ( refine_prediction(image, cr) ) { crystal_set_user_flag(cr, 1); n_bad++; + continue; } if ( (indm & INDEXING_CHECK_CELL_COMBINATIONS) -- cgit v1.2.3 From 5bb6d2e69535498a2e4aaf41c51c7566cd15a170 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Tue, 27 Jun 2017 14:33:57 +0200 Subject: Dummy setup bits for INDEXING_NONE This doesn't do much except for preventing a possibly confusing warning message about "Failed to prepare indexing method none". You can actually do --indexing=none,asdf, but don't tell anyone... --- libcrystfel/src/index.c | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/index.c b/libcrystfel/src/index.c index 9bcf5e15..446f57a3 100644 --- a/libcrystfel/src/index.c +++ b/libcrystfel/src/index.c @@ -92,6 +92,10 @@ static void *prepare_method(IndexingMethod *m, UnitCell *cell, switch ( *m & INDEXING_METHOD_MASK ) { + case INDEXING_NONE : + priv = "none"; + break; + case INDEXING_DIRAX : priv = dirax_prepare(m, cell, det, ltl); break; @@ -130,7 +134,9 @@ static void *prepare_method(IndexingMethod *m, UnitCell *cell, return NULL; } - STATUS("Prepared indexing method %s\n", str); + if ( *m != INDEXING_NONE ) { + STATUS("Prepared indexing method %s\n", str); + } free(str); if ( in != *m ) { -- cgit v1.2.3 From 91fa7dc2b06a9abc19366514b27e3b514f78b74f Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 5 Jul 2017 10:20:28 +0200 Subject: Rearrange quaternion declarations so that gtk-doc can cope The utils file contains two gtk-doc sections, "utils" and "quaternion", with a change of section halfway through. It seems that gtk-doc doesn't like it if the declarations are all mixed together. --- libcrystfel/src/utils.c | 276 ++++++++++++++++++++++++------------------------ libcrystfel/src/utils.h | 86 +++++++-------- 2 files changed, 179 insertions(+), 183 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/utils.c b/libcrystfel/src/utils.c index 461da0cd..90d81a71 100644 --- a/libcrystfel/src/utils.c +++ b/libcrystfel/src/utils.c @@ -387,144 +387,6 @@ int poisson_noise(gsl_rng *rng, double expected) } -/** - * SECTION:quaternion - * @short_description: Simple quaternion handling - * @title: Quaternion - * @section_id: - * @see_also: - * @include: "utils.h" - * @Image: - * - * There is a simple quaternion structure in CrystFEL. At the moment, it is - * only used when simulating patterns, as an argument to cell_rotate() to - * orient the unit cell. - */ - -/** - * quaternion_modulus: - * @q: A %quaternion - * - * If a quaternion represents a pure rotation, its modulus should be unity. - * - * Returns: the modulus of the given quaternion. - **/ -double quaternion_modulus(struct quaternion q) -{ - return sqrt(q.w*q.w + q.x*q.x + q.y*q.y + q.z*q.z); -} - - -/** - * normalise_quaternion: - * @q: A %quaternion - * - * Rescales the quaternion such that its modulus is unity. - * - * Returns: the normalised version of @q - **/ -struct quaternion normalise_quaternion(struct quaternion q) -{ - double mod; - struct quaternion r; - - mod = quaternion_modulus(q); - - r.w = q.w / mod; - r.x = q.x / mod; - r.y = q.y / mod; - r.z = q.z / mod; - - return r; -} - - -/** - * random_quaternion: - * @rng: A GSL random number generator to use - * - * Returns: a randomly generated, normalised, quaternion. - **/ -struct quaternion random_quaternion(gsl_rng *rng) -{ - struct quaternion q; - - q.w = 2.0*gsl_rng_uniform(rng) - 1.0; - q.x = 2.0*gsl_rng_uniform(rng) - 1.0; - q.y = 2.0*gsl_rng_uniform(rng) - 1.0; - q.z = 2.0*gsl_rng_uniform(rng) - 1.0; - q = normalise_quaternion(q); - - return q; -} - - -/** - * quaternion_valid: - * @q: A %quaternion - * - * Checks if the given quaternion is normalised. - * - * This function performs a nasty floating point comparison of the form - * (modulus > 0.999) && (modulus < 1.001), and so should not be - * relied upon to spot anything other than the most obvious input error. - * - * Returns: 1 if the quaternion is normalised, 0 if not. - **/ -int quaternion_valid(struct quaternion q) -{ - double qmod; - - qmod = quaternion_modulus(q); - - /* Modulus = 1 to within some tolerance? - * Nasty allowance for floating-point accuracy follows... */ - if ( (qmod > 0.999) && (qmod < 1.001) ) return 1; - - return 0; -} - - -/** - * quat_rot - * @q: A vector (in the form of a "struct rvec") - * @z: A %quaternion - * - * Rotates a vector according to a quaternion. - * - * Returns: A rotated version of @p. - **/ -struct rvec quat_rot(struct rvec q, struct quaternion z) -{ - struct rvec res; - double t01, t02, t03, t11, t12, t13, t22, t23, t33; - - t01 = z.w*z.x; - t02 = z.w*z.y; - t03 = z.w*z.z; - t11 = z.x*z.x; - t12 = z.x*z.y; - t13 = z.x*z.z; - t22 = z.y*z.y; - t23 = z.y*z.z; - t33 = z.z*z.z; - - res.u = (1.0 - 2.0 * (t22 + t33)) * q.u - + (2.0 * (t12 + t03)) * q.v - + (2.0 * (t13 - t02)) * q.w; - - res.v = (2.0 * (t12 - t03)) * q.u - + (1.0 - 2.0 * (t11 + t33)) * q.v - + (2.0 * (t01 + t23)) * q.w; - - res.w = (2.0 * (t02 + t13)) * q.u - + (2.0 * (t23 - t01)) * q.v - + (1.0 - 2.0 * (t11 + t22)) * q.w; - - return res; -} - - /* Return non-zero if c is in delims */ static int assplode_isdelim(const char c, const char *delims) { @@ -691,3 +553,141 @@ void utils_fudge_gslcblas() { STATUS("%p\n", cblas_sgemm); } + + +/** + * SECTION:quaternion + * @short_description: Simple quaternion handling + * @title: Quaternion + * @section_id: + * @see_also: + * @include: "utils.h" + * @Image: + * + * There is a simple quaternion structure in CrystFEL. At the moment, it is + * only used when simulating patterns, as an argument to cell_rotate() to + * orient the unit cell. + */ + +/** + * quaternion_modulus: + * @q: A %quaternion + * + * If a quaternion represents a pure rotation, its modulus should be unity. + * + * Returns: the modulus of the given quaternion. + **/ +double quaternion_modulus(struct quaternion q) +{ + return sqrt(q.w*q.w + q.x*q.x + q.y*q.y + q.z*q.z); +} + + +/** + * normalise_quaternion: + * @q: A %quaternion + * + * Rescales the quaternion such that its modulus is unity. + * + * Returns: the normalised version of @q + **/ +struct quaternion normalise_quaternion(struct quaternion q) +{ + double mod; + struct quaternion r; + + mod = quaternion_modulus(q); + + r.w = q.w / mod; + r.x = q.x / mod; + r.y = q.y / mod; + r.z = q.z / mod; + + return r; +} + + +/** + * random_quaternion: + * @rng: A GSL random number generator to use + * + * Returns: a randomly generated, normalised, quaternion. + **/ +struct quaternion random_quaternion(gsl_rng *rng) +{ + struct quaternion q; + + q.w = 2.0*gsl_rng_uniform(rng) - 1.0; + q.x = 2.0*gsl_rng_uniform(rng) - 1.0; + q.y = 2.0*gsl_rng_uniform(rng) - 1.0; + q.z = 2.0*gsl_rng_uniform(rng) - 1.0; + q = normalise_quaternion(q); + + return q; +} + + +/** + * quaternion_valid: + * @q: A %quaternion + * + * Checks if the given quaternion is normalised. + * + * This function performs a nasty floating point comparison of the form + * (modulus > 0.999) && (modulus < 1.001), and so should not be + * relied upon to spot anything other than the most obvious input error. + * + * Returns: 1 if the quaternion is normalised, 0 if not. + **/ +int quaternion_valid(struct quaternion q) +{ + double qmod; + + qmod = quaternion_modulus(q); + + /* Modulus = 1 to within some tolerance? + * Nasty allowance for floating-point accuracy follows... */ + if ( (qmod > 0.999) && (qmod < 1.001) ) return 1; + + return 0; +} + + +/** + * quat_rot + * @q: A vector (in the form of a "struct rvec") + * @z: A %quaternion + * + * Rotates a vector according to a quaternion. + * + * Returns: A rotated version of @p. + **/ +struct rvec quat_rot(struct rvec q, struct quaternion z) +{ + struct rvec res; + double t01, t02, t03, t11, t12, t13, t22, t23, t33; + + t01 = z.w*z.x; + t02 = z.w*z.y; + t03 = z.w*z.z; + t11 = z.x*z.x; + t12 = z.x*z.y; + t13 = z.x*z.z; + t22 = z.y*z.y; + t23 = z.y*z.z; + t33 = z.z*z.z; + + res.u = (1.0 - 2.0 * (t22 + t33)) * q.u + + (2.0 * (t12 + t03)) * q.v + + (2.0 * (t13 - t02)) * q.w; + + res.v = (2.0 * (t12 - t03)) * q.u + + (1.0 - 2.0 * (t11 + t33)) * q.v + + (2.0 * (t01 + t23)) * q.w; + + res.w = (2.0 * (t02 + t13)) * q.u + + (2.0 * (t23 - t01)) * q.v + + (1.0 - 2.0 * (t11 + t22)) * q.w; + + return res; +} diff --git a/libcrystfel/src/utils.h b/libcrystfel/src/utils.h index 4955f875..a759ff15 100644 --- a/libcrystfel/src/utils.h +++ b/libcrystfel/src/utils.h @@ -61,44 +61,10 @@ #define THOMSON_LENGTH (2.81794e-15) -/* ------------------------------ Quaternions ------------------------------- */ - -/** - * quaternion: - * - * - * struct quaternion - * { - * double w - * double x - * double y - * double z - * }; - * - * - * A structure representing a quaternion. - * - **/ -struct quaternion; - -struct quaternion { - double w; - double x; - double y; - double z; -}; - #ifdef __cplusplus extern "C" { #endif -extern struct quaternion normalise_quaternion(struct quaternion q); -extern double quaternion_modulus(struct quaternion q); -extern struct quaternion random_quaternion(gsl_rng *rng); -extern int quaternion_valid(struct quaternion q); -extern struct rvec quat_rot(struct rvec q, struct quaternion z); - - /* --------------------------- Useful functions ----------------------------- */ extern void show_matrix_eqn(gsl_matrix *M, gsl_vector *v); @@ -126,17 +92,6 @@ extern double flat_noise(gsl_rng *rng, double expected, double width); extern double gaussian_noise(gsl_rng *rng, double expected, double stddev); extern int poisson_noise(gsl_rng *rng, double expected); -/* Keep these ones inline, to avoid function call overhead */ -static inline struct quaternion invalid_quaternion(void) -{ - struct quaternion quat; - quat.w = 0.0; - quat.x = 0.0; - quat.y = 0.0; - quat.z = 0.0; - return quat; -} - static inline double distance(double x1, double y1, double x2, double y2) { return sqrt((x2-x1)*(x2-x1) + (y2-y1)*(y2-y1)); @@ -275,6 +230,47 @@ extern char *safe_basename(const char *in); #define unlikely(x) (x) #endif + +/* ------------------------------ Quaternions ------------------------------- */ + +/** + * quaternion: + * @w: component + * @x: component + * @y: component + * @z: component + * + * A structure representing a quaternion. + * + **/ +struct quaternion; + +struct quaternion { + double w; + double x; + double y; + double z; +}; + +extern struct quaternion normalise_quaternion(struct quaternion q); +extern double quaternion_modulus(struct quaternion q); +extern struct quaternion random_quaternion(gsl_rng *rng); +extern int quaternion_valid(struct quaternion q); +extern struct rvec quat_rot(struct rvec q, struct quaternion z); + +/* Keep these ones inline, to avoid function call overhead */ +static inline struct quaternion invalid_quaternion(void) +{ + struct quaternion quat; + quat.w = 0.0; + quat.x = 0.0; + quat.y = 0.0; + quat.z = 0.0; + return quat; +} + + + #ifdef __cplusplus } #endif -- cgit v1.2.3 From bb42944e6c79d5a81dc30840ceebe70dc0d96658 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 5 Jul 2017 11:46:24 +0200 Subject: Update docs --- libcrystfel/src/detector.h | 46 ++++++++++++++++++++-- libcrystfel/src/image.h | 95 ++++++++++++++++++++++++++++------------------ libcrystfel/src/index.h | 1 + 3 files changed, 102 insertions(+), 40 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/detector.h b/libcrystfel/src/detector.h index e9dd1154..a2be2b47 100644 --- a/libcrystfel/src/detector.h +++ b/libcrystfel/src/detector.h @@ -3,12 +3,12 @@ * * Detector properties * - * Copyright © 2012-2016 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * Copyright © 2012 Richard Kirian * * Authors: - * 2009-2016 Thomas White + * 2009-2017 Thomas White * 2011-2012 Richard Kirian * 2014 Valerio Mariani * 2011 Andrew Aquila @@ -42,7 +42,6 @@ struct rg_collection; struct detector; struct panel; struct badregion; -struct detector; struct beam_params; struct hdfile; struct event; @@ -80,6 +79,47 @@ struct rg_collection }; +/** + * panel: + * @name: Text name for the panel (fixed length array) + * @cnx: Location of corner, in pixels, x coordinate + * @cny: Location of corner, in pixels, y coordinate + * @coffset: The offset to be applied from @clen (which may come from elsewhere) + * @clen: The distance from the interaction point to the corner of the first pixel + * @clen_from: Location to get @clen from, e.g. from HDF5 file + * @mask: Location of mask data + * @mask_file: Filename for mask data + * @satmap: Location of per-pixel saturation map + * @satmap_file: Filename for saturation map + * @res: Resolution of panel in pixels per metre + * @badrow: Readout direction (for filtering out clusters of peaks) + * @no_index: Non-zero if panel is entirely "bad" + * @adu_per_photon: Number of detector intensity units per photon + * @adu_per_eV: Number of detector intensity units per eV of photon energy + * @max_adu: Saturation value + * @dim_structure: Dimension structure + * @fsx: Real-space x-direction of data fast-scan direction + * @fsy: Real-space y-direction of data fast-scan direction + * @fsz: Real-space z-direction of data fast-scan direction + * @ssx: Real-space x-direction of data slow-scan direction + * @ssy: Real-space y-direction of data slow-scan direction + * @ssz: Real-space z-direction of data slow-scan direction + * @rail_x: x direction of camera length "rail" + * @rail_y: y direction of camera length "rail" + * @rail_z: z direction of camera length "rail" + * @clen_for_centering: Value of clen (without coffset) at which beam is centered + * @xfs: Data fast-scan direction of real-space x-direction + * @yfs: Data fast-scan direction of real-space y-direction + * @xss: Data slow-scan direction of real-space x-direction + * @yss: Data slow-scan direction of real-space y-direction + * @orig_min_fs: Minimum fs coordinate of data in file + * @orig_max_fs: Maximum fs coordinate of data in file + * @orig_min_ss: Minimum ss coordinate of data in file (inclusive) + * @orig_max_ss: Maximum ss coordinate of data in file (inclusive) + * @data: Location of data in file + * @w: Width of panel + * @h: Height of panel + */ struct panel { char name[1024]; /* Name for this panel */ diff --git a/libcrystfel/src/image.h b/libcrystfel/src/image.h index ead5cd4d..9719bb59 100644 --- a/libcrystfel/src/image.h +++ b/libcrystfel/src/image.h @@ -70,7 +70,26 @@ typedef enum { SPECTRUM_TWOCOLOUR } SpectrumType; -/* Structure describing a feature in an image */ + +/** + * imagefeature: + * @parent: Image this feature belongs to + * @fs: Fast scan coordinate + * @ss: Slow scan coordinate + * @p: Pointer to panel + * @intensity: Intensity of peak + * @rx: Reciprocal x coordinate in m^-1 + * @ry: Reciprocal y coordinate in m^-1 + * @rz: Reciprocal z coordinate in m^-1 + * @name: Text name for feature + * + * Represents a peak in an image. + * + * Note carefully that the @fs and @ss coordinates are the distances, measured + * in pixels, from the corner of the panel. They are NOT pixel indices. + * If the peak is in the middle of the first pixel, its coordinates would be + * 0.5,0.5. + */ struct imagefeature { struct image *parent; @@ -84,10 +103,10 @@ struct imagefeature { double ry; double rz; - /* Internal use only */ - int valid; - const char *name; + + /*< private >*/ + int valid; }; @@ -110,44 +129,46 @@ struct sample }; +/** + * beam_params: + * @photon_energy: eV per photon + * @photon_energy_from: HDF5 dataset name + * @photon_energy_scale: Scale factor for photon energy, if it comes from HDF5 + */ struct beam_params { - double photon_energy; /* eV per photon */ - char *photon_energy_from; /* HDF5 dataset name */ - double photon_energy_scale; /* Scale factor for photon energy, if the - * energy is to be from the HDF5 file */ + double photon_energy; + char *photon_energy_from; + double photon_energy_scale; }; /** * image: - * - * - * struct image - * { - * Crystal **crystals; - * int n_crystals; - * IndexingMethod indexed_by; - * - * struct detector *det; - * struct beam_params *beam; - * char *filename; - * const struct imagefile_field_list *copyme; - * - * int id; - * - * double lambda; - * double div; - * double bw; - * - * int width; - * int height; - * - * long long int num_peaks; - * long long int num_saturated_peaks; - * ImageFeatureList *features; - * }; - * + * @crystals: Array of crystals in the image + * @n_crystals: The number of crystals in the image + * @indexed_by: Indexing method which indexed this pattern + * @det: Detector structure + * @beam: Beam parameters structure + * @filename: Filename for the image file + * @copyme: Fields to copy from the image file to the stream + * @id: ID number of the thread handling this image + * @serial: Serial number for this image + * @lambda: Wavelength + * @div: Divergence + * @bw: Bandwidth + * @num_peaks: The number of peaks + * @num_saturated_peaks: The number of saturated peaks + * @features: The peaks found in the image + * @dp: The image data, by panel + * @bad: The bad pixel mask, array by panel + * @sat: The per-pixel saturation mask, array by panel + * @event: Event ID for the image + * @stuff_from_stream: Items read back from the stream + * @avg_clen: Mean of camera length values for all panels + * @spectrum: Spectrum information + * @nsamples: Number of spectrum samples + * @spectrum_size: SIze of spectrum array * * The field data contains the raw image data, if it * is currently available. The data might be available throughout the @@ -204,8 +225,8 @@ struct image { double bw; /* Bandwidth as a fraction */ /* Detected peaks */ - long long int num_peaks; - long long int num_saturated_peaks; + long long num_peaks; + long long num_saturated_peaks; ImageFeatureList *features; }; diff --git a/libcrystfel/src/index.h b/libcrystfel/src/index.h index 8ce553ae..8dd7de21 100644 --- a/libcrystfel/src/index.h +++ b/libcrystfel/src/index.h @@ -75,6 +75,7 @@ * @INDEXING_SIMULATION: Dummy value * @INDEXING_DEBUG: Results injector for debugging * @INDEXING_ASDF: Use in-built "asdf" indexer + * @INDEXING_ERROR: Special value for unrecognised indexing engine name * @INDEXING_CHECK_CELL_COMBINATIONS: Check linear combinations of unit cell * axes for agreement with given cell. * @INDEXING_CHECK_CELL_AXES: Check unit cell axes for agreement with given -- cgit v1.2.3 From 608cebf106489636a487d9a96092d073a70fc660 Mon Sep 17 00:00:00 2001 From: Valerio Mariani Date: Wed, 5 Jul 2017 17:14:58 +0200 Subject: Update to peakfinder8, with bug fixed and new functionality. Code synced with OnDA and Oleksandr's programs --- libcrystfel/src/peakfinder8.c | 436 +++++++++++++++++++----------------------- 1 file changed, 199 insertions(+), 237 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/peakfinder8.c b/libcrystfel/src/peakfinder8.c index 76d7c13c..417465df 100644 --- a/libcrystfel/src/peakfinder8.c +++ b/libcrystfel/src/peakfinder8.c @@ -38,6 +38,7 @@ #include "peakfinder8.h" +// CrystFEL-only block 1 struct radius_maps { float **r_maps; @@ -52,25 +53,27 @@ struct peakfinder_mask }; +struct peakfinder_panel_data +{ + float **panel_data; + int *panel_h; + int *panel_w; + int num_panels; +}; +// End of CrystFEL-only block 1 + + struct radial_stats { float *roffset; float *rthreshold; + float *lthreshold; float *rsigma; int *rcount; int n_rad_bins; }; -struct peakfinder_panel_data -{ - float **panel_data; - int *panel_h; - int *panel_w; - int num_panels; -}; - - struct peakfinder_intern_data { char *pix_in_peak_map; @@ -94,36 +97,7 @@ struct peakfinder_peak_data }; -static double compute_r_sigma(float *rsigma, int *rcount, float *roffset, - int i) -{ - return sqrt(rsigma[i] / rcount[i] - - ((roffset[i] / rcount[i]) * - (roffset[i] / rcount[i]))); -} - - -static void update_radial_stats(float *roffset, float *rsigma, int *rcount, - float value, int curr_radius) -{ - roffset[curr_radius] += value; - rsigma[curr_radius] += (value * value); - rcount[curr_radius] += 1; -} - - -static float get_radius(struct panel p, float fs, float ss) -{ - float x, y; - - /* Calculate 3D position of given position, in m */ - x = (p.cnx + fs * p.fsx + ss * p.ssx); - y = (p.cny + fs * p.fsy + ss * p.ssy); - - return sqrt(x * x + y * y); -} - - +// CrystFEL-only block 2 static struct radius_maps *compute_radius_maps(struct detector *det) { int i, u, iss, ifs; @@ -160,15 +134,17 @@ static struct radius_maps *compute_radius_maps(struct detector *det) for ( ifs=0; ifsr_maps[i][rmi] = get_radius(p, ifs, iss); + x = (p.cnx + ifs * p.fsx + iss * p.ssx); + y = (p.cny + ifs * p.fsy + iss * p.ssy); + + rm->r_maps[i][rmi] = sqrt(x * x + y * y); } } - } - return rm; } @@ -222,7 +198,6 @@ static struct peakfinder_mask *create_peakfinder_mask(struct image *img, } } - } } } @@ -242,45 +217,12 @@ static void free_peakfinder_mask(struct peakfinder_mask * pfmask) } -static int compute_num_radial_bins(int num_panels, int *w, int *h, - float **r_maps) -{ - - float max_r; - int max_r_int; - int m; - - max_r = -1e9; - - for ( m=0 ; m max_r ) { - max_r = r_maps[m][pidx]; - } - } - } - } - - max_r_int = (int)ceil(max_r) + 1; - - return max_r_int; -} - - static struct peakfinder_panel_data *allocate_panel_data(int num_panels) { struct peakfinder_panel_data *pfdata; - - pfdata = (struct peakfinder_panel_data *)malloc( - sizeof(struct peakfinder_panel_data)); + pfdata = (struct peakfinder_panel_data *)malloc(sizeof(struct peakfinder_panel_data)); if ( pfdata == NULL ) { return NULL; } @@ -321,9 +263,26 @@ static void free_panel_data(struct peakfinder_panel_data *pfdata) } -static struct radial_stats *allocate_radial_stats(int num_rad_bins) +static void compute_num_radial_bins(int w, int h, float *r_map, float *max_r) { - struct radial_stats *rstats; + int ifs, iss; + int pidx; + + for ( iss=0 ; iss *max_r ) { + *max_r = r_map[pidx]; + } + } + } +} +// End of CrystFEL-only block 2 + + +static struct radial_stats* allocate_radial_stats(int num_rad_bins) +{ + struct radial_stats* rstats; rstats = (struct radial_stats *)malloc(sizeof(struct radial_stats)); if ( rstats == NULL ) { @@ -343,11 +302,20 @@ static struct radial_stats *allocate_radial_stats(int num_rad_bins) return NULL; } + rstats->lthreshold = (float *)malloc(num_rad_bins*sizeof(float)); + if ( rstats->lthreshold == NULL ) { + free(rstats); + free(rstats->rthreshold); + free(rstats->roffset); + return NULL; + } + rstats->rsigma = (float *)malloc(num_rad_bins*sizeof(float)); if ( rstats->rsigma == NULL ) { free(rstats); free(rstats->roffset); free(rstats->rthreshold); + free(rstats->lthreshold); return NULL; } @@ -356,6 +324,7 @@ static struct radial_stats *allocate_radial_stats(int num_rad_bins) free(rstats); free(rstats->roffset); free(rstats->rthreshold); + free(rstats->lthreshold); free(rstats->rsigma); return NULL; } @@ -370,6 +339,7 @@ static void free_radial_stats(struct radial_stats *rstats) { free(rstats->roffset); free(rstats->rthreshold); + free(rstats->lthreshold); free(rstats->rsigma); free(rstats->rcount); free(rstats); @@ -382,6 +352,7 @@ static void fill_radial_bins(float *data, float *r_map, char *mask, float *rthreshold, + float *lthreshold, float *roffset, float *rsigma, int *rcount) @@ -393,24 +364,15 @@ static void fill_radial_bins(float *data, float value; for ( iss = 0; isslthreshold[curr_r]) { + roffset[curr_r] += value; + rsigma[curr_r] += (value * value); + rcount[curr_r] += 1; } } } @@ -419,6 +381,7 @@ static void fill_radial_bins(float *data, static void compute_radial_stats(float *rthreshold, + float *lthreshold, float *roffset, float *rsigma, int *rcount, @@ -426,7 +389,6 @@ static void compute_radial_stats(float *rthreshold, float min_snr, float acd_threshold) { - int ri; float this_offset, this_sigma; @@ -436,19 +398,18 @@ static void compute_radial_stats(float *rthreshold, roffset[ri] = 0; rsigma[ri] = 0; rthreshold[ri] = 1e9; + lthreshold[ri] = -1e9; } else { - this_offset = roffset[ri]/rcount[ri]; - - this_sigma = compute_r_sigma(rsigma, - rcount, - roffset, - ri); + this_offset = roffset[ri] / rcount[ri]; + this_sigma = rsigma[ri] / rcount[ri] - (this_offset * this_offset); + if ( this_sigma >= 0 ) { + this_sigma = sqrt(this_sigma); + } roffset[ri] = this_offset; rsigma[ri] = this_sigma; - - rthreshold[ri] = roffset[ri] + - min_snr*rsigma[ri]; + rthreshold[ri] = roffset[ri] + min_snr*rsigma[ri]; + lthreshold[ri] = roffset[ri] - min_snr*rsigma[ri]; if ( rthreshold[ri] < acd_threshold ) { rthreshold[ri] = acd_threshold; @@ -463,8 +424,7 @@ struct peakfinder_peak_data *allocate_peak_data(int max_num_peaks) { struct peakfinder_peak_data *pkdata; - pkdata = (struct peakfinder_peak_data*)malloc( - sizeof(struct peakfinder_peak_data)); + pkdata = (struct peakfinder_peak_data*)malloc(sizeof(struct peakfinder_peak_data)); if ( pkdata == NULL ) { return NULL; } @@ -500,7 +460,6 @@ struct peakfinder_peak_data *allocate_peak_data(int max_num_peaks) return NULL; } - pkdata->tot_i = (float *)malloc(max_num_peaks*sizeof(float)); if ( pkdata->tot_i == NULL ) { free(pkdata->npix); @@ -551,8 +510,7 @@ struct peakfinder_peak_data *allocate_peak_data(int max_num_peaks) } -static void free_peak_data(struct peakfinder_peak_data *pkdata) -{ +static void free_peak_data(struct peakfinder_peak_data *pkdata) { free(pkdata->npix); free(pkdata->com_fs); free(pkdata->com_ss); @@ -565,14 +523,13 @@ static void free_peak_data(struct peakfinder_peak_data *pkdata) } -static struct peakfinder_intern_data *allocate_peakfinder_intern_data( - int data_size, int max_pix_count) +static struct peakfinder_intern_data *allocate_peakfinder_intern_data(int data_size, + int max_pix_count) { struct peakfinder_intern_data *intern_data; - intern_data = (struct peakfinder_intern_data *)malloc( - sizeof(struct peakfinder_intern_data)); + intern_data = (struct peakfinder_intern_data *)malloc(sizeof(struct peakfinder_intern_data)); if ( intern_data == NULL ) { return NULL; } @@ -621,7 +578,9 @@ static void free_peakfinder_intern_data(struct peakfinder_intern_data *pfid) } -static void peak_search(int p, struct peakfinder_intern_data *pfinter, + +static void peak_search(int p, + struct peakfinder_intern_data *pfinter, float *copy, char *mask, float *r_map, float *rthreshold, float *roffset, int *num_pix_in_peak, int asic_size_fs, @@ -641,17 +600,15 @@ static void peak_search(int p, struct peakfinder_intern_data *pfinter, int search_ss[9] = { 0, -1, -1, -1, 0, 0, 1, 1, 1 }; int search_n = 9; + // Loop through search pattern for ( k=0; kinfs[p] + search_fs[k]) < 0 ) - continue; - if ( (pfinter->infs[p] + search_fs[k]) >= asic_size_fs ) - continue; - if ( (pfinter->inss[p] + search_ss[k]) < 0 ) - continue; - if ( (pfinter->inss[p] + search_ss[k]) >= asic_size_ss ) - continue; + if ( (pfinter->infs[p] + search_fs[k]) < 0 ) continue; + if ( (pfinter->infs[p] + search_fs[k]) >= asic_size_fs ) continue; + if ( (pfinter->inss[p] + search_ss[k]) < 0 ) continue; + if ( (pfinter->inss[p] + search_ss[k]) >= asic_size_ss ) continue; + // Neighbour point in big array curr_fs = pfinter->infs[p] + search_fs[k] + aifs * asic_size_fs; curr_ss = pfinter->inss[p] + search_ss[k] + aiss * asic_size_ss; pi = curr_fs + curr_ss * num_pix_fs; @@ -659,28 +616,23 @@ static void peak_search(int p, struct peakfinder_intern_data *pfinter, curr_radius = (int)rint(r_map[pi]); curr_threshold = rthreshold[curr_radius]; + // Above threshold? if ( copy[pi] > curr_threshold && pfinter->pix_in_peak_map[pi] == 0 - && mask[pi] != 0 ) - { - + && mask[pi] != 0 ) { curr_i = copy[pi] - roffset[curr_radius]; - *sum_i += curr_i; - *sum_com_fs += curr_i * ((float)curr_fs); - *sum_com_ss += curr_i * ((float)curr_ss); - - pfinter->inss[*num_pix_in_peak] = pfinter->inss[p] + - search_ss[k]; - pfinter->infs[*num_pix_in_peak] = pfinter->infs[p] + - search_fs[k]; + *sum_com_fs += curr_i * ((float)curr_fs); // for center of mass x + *sum_com_ss += curr_i * ((float)curr_ss); // for center of mass y + pfinter->inss[*num_pix_in_peak] = pfinter->inss[p] + search_ss[k]; + pfinter->infs[*num_pix_in_peak] = pfinter->infs[p] + search_fs[k]; pfinter->pix_in_peak_map[pi] = 1; if ( *num_pix_in_peak < max_pix_count ) { - pfinter->peak_pixels[*num_pix_in_peak] = pi; + pfinter->peak_pixels[*num_pix_in_peak] = pi; } - *num_pix_in_peak = *num_pix_in_peak+1; + *num_pix_in_peak = *num_pix_in_peak + 1; } } } @@ -691,9 +643,9 @@ static void search_in_ring(int ring_width, int com_fs_int, int com_ss_int, float *rthreshold, float *roffset, char *pix_in_peak_map, char *mask, int asic_size_fs, int asic_size_ss, int aifs, int aiss, - int num_pix_fs,float *local_sigma, - float *local_offset, float *background_max_i, - int com_idx, int local_bg_radius) + int num_pix_fs,float *local_sigma, float *local_offset, + float *background_max_i, int com_idx, + int local_bg_radius) { int ssj, fsi; float pix_radius; @@ -721,30 +673,30 @@ static void search_in_ring(int ring_width, int com_fs_int, int com_ss_int, for ( ssj = -ring_width ; ssj= asic_size_fs ) - continue; - if ( (com_ss_int + ssj) < 0 ) - continue; + // Within-ASIC check + if ( (com_fs_int + fsi) < 0 ) continue; + if ( (com_fs_int + fsi) >= asic_size_fs ) continue; + if ( (com_ss_int + ssj) < 0 ) continue; if ( (com_ss_int + ssj) >= asic_size_ss ) - continue; + continue; + // Within outer ring check pix_radius = sqrt(fsi * fsi + ssj * ssj); - if ( pix_radius>ring_width ) - continue; + if ( pix_radius>ring_width ) continue; - curr_fs = com_fs_int +fsi + aifs * asic_size_fs; - curr_ss = com_ss_int +ssj + aiss * asic_size_ss; + // Position of this point in data stream + curr_fs = com_fs_int + fsi + aifs * asic_size_fs; + curr_ss = com_ss_int + ssj + aiss * asic_size_ss; pi = curr_fs + curr_ss * num_pix_fs; - curr_radius = rint(r_map[pi]); + curr_radius = (int)rint(r_map[pi]); curr_threshold = rthreshold[curr_radius]; + + // Intensity above background ??? just intensity? curr_i = copy[pi]; - if ( copy[pi] < curr_threshold - && pix_in_peak_map[pi] == 0 - && mask[pi] != 0 ) { + // Keep track of value and value-squared for offset and sigma calculation + if ( curr_i < curr_threshold && pix_in_peak_map[pi] == 0 && mask[pi] != 0 ) { np_sigma++; sum_i += curr_i; @@ -758,13 +710,17 @@ static void search_in_ring(int ring_width, int com_fs_int, int com_ss_int, } } + // Calculate local background and standard deviation if ( np_sigma != 0 ) { *local_offset = sum_i / np_sigma; - *local_sigma = sqrt(sum_i_squared / np_sigma - - ((sum_i / np_sigma) * (sum_i / np_sigma))); + *local_sigma = sum_i_squared / np_sigma - (*local_offset * *local_offset); + if (*local_sigma >= 0) { + *local_sigma = sqrt(*local_sigma); + } else { + *local_sigma = 0.01; + } } else { - - local_radius = rint(r_map[(int)rint(com_idx)]); + local_radius = (int)rint(r_map[(int)rint(com_idx)]); *local_offset = roffset[local_radius]; *local_sigma = 0.01; } @@ -784,6 +740,7 @@ static void process_panel(int asic_size_fs, int asic_size_ss, int num_pix_fs, int pxss, pxfs; int num_pix_in_peak; + // Loop over pixels within a module for ( pxss=1 ; pxss curr_thresh - && pfinter->pix_in_peak_map[pxidx] == 0 ) { + && pfinter->pix_in_peak_map[pxidx] == 0 + && mask[pxidx] != 0 ) { //??? not sure if needed + // This might be the start of a new peak - start searching float sum_com_fs, sum_com_ss; float sum_i; float peak_com_fs, peak_com_ss; @@ -810,7 +768,6 @@ static void process_panel(int asic_size_fs, int asic_size_ss, int num_pix_fs, float peak_snr; float local_sigma, local_offset; float background_max_i; - float f_background_thresh; int lt_num_pix_in_pk; int ring_width; int peak_idx; @@ -820,16 +777,18 @@ static void process_panel(int asic_size_fs, int asic_size_ss, int num_pix_fs, pfinter->infs[0] = pxfs; pfinter->inss[0] = pxss; pfinter->peak_pixels[0] = pxidx; - num_pix_in_peak = 1; + num_pix_in_peak = 0; //y 1; sum_i = 0; sum_com_fs = 0; sum_com_ss = 0; + // Keep looping until the pixel count within this peak does not change do { lt_num_pix_in_pk = num_pix_in_peak; - for ( p=0; p max_pix_count) { - continue; - } + // Too many or too few pixels means ignore this 'peak'; move on now + if ( num_pix_in_peak < min_pix_count || num_pix_in_peak > max_pix_count ) continue; + + // If for some reason sum_i is 0 - it's better to skip + if ( fabs(sum_i) < 1e-10 ) continue; + // Calculate center of mass for this peak from initial peak search peak_com_fs = sum_com_fs / fabs(sum_i); peak_com_ss = sum_com_ss / fabs(sum_i); - com_idx = rint(peak_com_fs) + - rint(peak_com_ss) * num_pix_fs; + com_idx = (int)rint(peak_com_fs) + (int)rint(peak_com_ss) * num_pix_fs; - peak_com_fs_int = (int)rint(peak_com_fs) - - aifs * asic_size_fs; - peak_com_ss_int = (int)rint(peak_com_ss) - - aiss * asic_size_ss; + peak_com_fs_int = (int)rint(peak_com_fs) - aifs * asic_size_fs; + peak_com_ss_int = (int)rint(peak_com_ss) - aiss * asic_size_ss; + // Calculate the local signal-to-noise ratio and local background in an annulus around + // this peak (excluding pixels which look like they might be part of another peak) local_sigma = 0.0; local_offset = 0.0; background_max_i = 0.0; @@ -884,6 +844,7 @@ static void process_panel(int asic_size_fs, int asic_size_ss, int num_pix_fs, &background_max_i, com_idx, local_bg_radius); + // Re-integrate (and re-centroid) peak using local background estimates peak_tot_i = 0; pk_tot_i_raw = 0; peak_max_i = 0; @@ -891,71 +852,71 @@ static void process_panel(int asic_size_fs, int asic_size_ss, int num_pix_fs, sum_com_fs = 0; sum_com_ss = 0; - for ( peak_idx = 1 ; - peak_idx < num_pix_in_peak - && peak_idx <= max_pix_count ; - peak_idx++ ) { + for ( peak_idx = 0 ; + peak_idx < num_pix_in_peak && peak_idx < max_pix_count ; + peak_idx++ ) { int curr_idx; float curr_i; float curr_i_raw; int curr_fs, curr_ss; - // BUG HERE, I THINK. PEAK_PIXELS - // HAS SIZE MAX_PIX_COUNT, BUT - // IN THE FOLLOWING LINE PEAK_IDX - // CAN HAVE VALUE EXACTLY MAX_PEAK_COUNT - // (SEE THE FOR LOOP ABOVE) - curr_idx = - pfinter->peak_pixels[peak_idx]; - + curr_idx = pfinter->peak_pixels[peak_idx]; curr_i_raw = copy[curr_idx]; curr_i = curr_i_raw - local_offset; - peak_tot_i += curr_i; pk_tot_i_raw += curr_i_raw; - curr_fs = curr_idx % asic_size_fs; - curr_ss = curr_idx / asic_size_fs; + // Remember that curr_idx = curr_fs + curr_ss*num_pix_fs + curr_fs = curr_idx % num_pix_fs; + curr_ss = curr_idx / num_pix_fs; sum_com_fs += curr_i * ((float)curr_fs); sum_com_ss += curr_i * ((float)curr_ss); - if ( curr_i_raw > pk_max_i_raw ) - pk_max_i_raw = curr_i_raw; - if ( curr_i > peak_max_i ) - peak_max_i = curr_i; + if ( curr_i_raw > pk_max_i_raw ) pk_max_i_raw = curr_i_raw; + if ( curr_i > peak_max_i ) peak_max_i = curr_i; } + // This CAN happen! Better to skip... + if ( fabs(peak_tot_i) < 1e-10 ) continue; + peak_com_fs = sum_com_fs / fabs(peak_tot_i); peak_com_ss = sum_com_ss / fabs(peak_tot_i); - peak_snr = peak_tot_i / local_sigma; - - if (peak_snr < min_snr) { - continue; + // Calculate signal-to-noise and apply SNR criteria + if ( fabs(local_sigma) > 1e-10 ) { + peak_snr = peak_tot_i / local_sigma; + } else { + peak_snr = 0; } - f_background_thresh = 1; - f_background_thresh *= - background_max_i - local_offset; - if (peak_max_i < f_background_thresh) { - continue; - } + if (peak_snr < min_snr) continue; + + // Is the maximum intensity in the peak enough above intensity in background region to + // be a peak and not noise? The more pixels there are in the peak, the more relaxed we + // are about this criterion + //f_background_thresh = background_max_i - local_offset; //!!! Ofiget'! If I uncomment + // if (peak_max_i < f_background_thresh) { // these lines the result is + // different! + if (peak_max_i < background_max_i - local_offset) continue; + // This is a peak? If so, add info to peak list if ( num_pix_in_peak >= min_pix_count && num_pix_in_peak <= max_pix_count ) { - int peak_com_idx; - - if ( peak_tot_i == 0 ) { - continue; + // Bragg peaks in the mask + for ( peak_idx = 0 ; + peak_idx < num_pix_in_peak && + peak_idx < max_pix_count ; + peak_idx++ ) { + pfinter->pix_in_peak_map[pfinter->peak_pixels[peak_idx]] = 2; } - peak_com_idx = rint(peak_com_fs) + - rint(peak_com_ss) * - asic_size_fs; - + int peak_com_idx; + peak_com_idx = (int)rint(peak_com_fs) + (int)rint(peak_com_ss) * + num_pix_fs; + // Remember peak information if ( *peak_count < max_n_peaks ) { int pidx; @@ -969,10 +930,8 @@ static void process_panel(int asic_size_fs, int asic_size_ss, int num_pix_fs, max_i[pidx] = peak_max_i; sigma[pidx] = local_sigma; snr[pidx] = peak_snr; - *peak_count = *peak_count + 1; - } else { - *peak_count = *peak_count + 1; } + *peak_count += 1; } } } @@ -989,43 +948,32 @@ static int peakfinder8_base(float *roffset, float *rthreshold, float *com_ss, int *com_index, float *tot_i, float *max_i, float *sigma, float *snr, int min_pix_count, int max_pix_count, - int local_bg_radius, float min_snr) + int local_bg_radius, float min_snr, + char* outliersMask) { + int num_pix_fs, num_pix_ss, num_pix_tot; - int i, aifs, aiss; + int aifs, aiss; int peak_count; - float *copy; struct peakfinder_intern_data *pfinter; num_pix_fs = asic_size_fs * num_asics_fs; num_pix_ss = asic_size_ss * num_asics_ss; num_pix_tot = num_pix_fs * num_pix_ss; - copy = (float *)calloc(num_pix_tot, sizeof(float)); - if ( copy == NULL ) { - return 1; - } - - memcpy(copy, data, num_pix_tot*sizeof(float)); - - for (i = 0; i < num_pix_tot; i++) { - copy[i] *= mask[i]; - } - pfinter = allocate_peakfinder_intern_data(num_pix_tot, max_pix_count); if ( pfinter == NULL ) { - free(copy); return 1; } peak_count = 0; + // Loop over modules (nxn array) for ( aiss=0 ; aisspix_in_peak_map, num_pix_tot*sizeof(char)); + } + free_peakfinder_intern_data(pfinter); - free(copy); return 0; } @@ -1058,6 +1009,7 @@ int peakfinder8(struct image *img, int max_n_peaks, int num_found_peaks; int remaining_max_num_peaks; int iterations; + float max_r; iterations = 5; @@ -1090,10 +1042,17 @@ int peakfinder8(struct image *img, int max_n_peaks, pfdata->num_panels = img->det->n_panels; } - num_rad_bins = compute_num_radial_bins(img->det->n_panels, - pfdata->panel_w, - pfdata->panel_h, - rmaps->r_maps); + max_r = -1e9; + + for ( pi=0 ; pinum_panels ; pi++ ) { + + compute_num_radial_bins(pfdata->panel_w[pi], + pfdata->panel_h[pi], + rmaps->r_maps[pi], + &max_r); + } + + num_rad_bins = (int)ceil(max_r) + 1; rstats = allocate_radial_stats(num_rad_bins); if ( rstats == NULL ) { @@ -1123,6 +1082,7 @@ int peakfinder8(struct image *img, int max_n_peaks, rmaps->r_maps[pi], pfmask->masks[pi], rstats->rthreshold, + rstats->lthreshold, rstats->roffset, rstats->rsigma, rstats->rcount); @@ -1130,6 +1090,7 @@ int peakfinder8(struct image *img, int max_n_peaks, } compute_radial_stats(rstats->rthreshold, + rstats->lthreshold, rstats->roffset, rstats->rsigma, rstats->rcount, @@ -1182,7 +1143,8 @@ int peakfinder8(struct image *img, int max_n_peaks, min_pix_count, max_pix_count, local_bg_radius, - min_snr); + min_snr, + NULL); if ( ret != 0 ) { free_radius_maps(rmaps); -- cgit v1.2.3 From dc3395900fc3ce0d3961757628ff83ad6456be19 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Thu, 6 Jul 2017 17:04:40 +0200 Subject: Indexing engine private pointers should be void * --- libcrystfel/src/asdf.c | 14 +++++++------- libcrystfel/src/asdf.h | 24 ++++++++++++------------ libcrystfel/src/dirax.c | 12 ++++++------ libcrystfel/src/dirax.h | 13 ++++++------- libcrystfel/src/felix.c | 9 ++++----- libcrystfel/src/felix.h | 16 +++++++--------- libcrystfel/src/mosflm.c | 8 ++++---- libcrystfel/src/mosflm.h | 12 ++++++------ libcrystfel/src/xds.c | 14 +++++++------- libcrystfel/src/xds.h | 16 ++++++++-------- 10 files changed, 67 insertions(+), 71 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/asdf.c b/libcrystfel/src/asdf.c index 033688ba..0c43fe7a 100644 --- a/libcrystfel/src/asdf.c +++ b/libcrystfel/src/asdf.c @@ -4,12 +4,12 @@ * Alexandra's Superior Direction Finder, or * Algorithm Similar to DirAx, FFT-based * - * Copyright © 2014-2015 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2014-2017 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * * Authors: * 2014-2015 Alexandra Tolstikova - * 2015 Thomas White + * 2015,2017 Thomas White * * This file is part of CrystFEL. * @@ -1101,7 +1101,7 @@ static int index_refls(gsl_vector **reflections, int N_reflections, } -int run_asdf(struct image *image, IndexingPrivate *ipriv) +int run_asdf(struct image *image, void *ipriv) { int i, j; @@ -1202,8 +1202,8 @@ int run_asdf(struct image *image, IndexingPrivate *ipriv) } -IndexingPrivate *asdf_prepare(IndexingMethod *indm, UnitCell *cell, - struct detector *det, float *ltl) +void *asdf_prepare(IndexingMethod *indm, UnitCell *cell, + struct detector *det, float *ltl) { struct asdf_private *dp; int need_cell = 0; @@ -1232,11 +1232,11 @@ IndexingPrivate *asdf_prepare(IndexingMethod *indm, UnitCell *cell, dp->fftw = fftw_vars_new(); - return (IndexingPrivate *)dp; + return (void *)dp; } -void asdf_cleanup(IndexingPrivate *pp) +void asdf_cleanup(void *pp) { struct asdf_private *p; p = (struct asdf_private *)pp; diff --git a/libcrystfel/src/asdf.h b/libcrystfel/src/asdf.h index 1e492a6f..f130d63d 100644 --- a/libcrystfel/src/asdf.h +++ b/libcrystfel/src/asdf.h @@ -4,12 +4,12 @@ * Alexandra's Superior Direction Finder, or * Algorithm Similar to DirAx, FFT-based * - * Copyright © 2014-2015 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2014-2017 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * * Authors: * 2014-2015 Alexandra Tolstikova - * 2015 Thomas White + * 2015,2017 Thomas White * * This file is part of CrystFEL. * @@ -43,33 +43,33 @@ extern "C" { #ifdef HAVE_FFTW -extern int run_asdf(struct image *image, IndexingPrivate *ipriv); +extern int run_asdf(struct image *image, void *ipriv); -extern IndexingPrivate *asdf_prepare(IndexingMethod *indm, - UnitCell *cell, struct detector *det, - float *ltl); +extern void *asdf_prepare(IndexingMethod *indm, + UnitCell *cell, struct detector *det, + float *ltl); -extern void asdf_cleanup(IndexingPrivate *pp); +extern void asdf_cleanup(void *pp); #else /* HAVE_FFTW */ -int run_asdf(struct image *image, IndexingPrivate *ipriv) +int run_asdf(struct image *image, void *ipriv) { ERROR("This copy of CrystFEL was compiled without FFTW support.\n"); return 0; } -IndexingPrivate *asdf_prepare(IndexingMethod *indm, - UnitCell *cell, struct detector *det, - float *ltl) +void *asdf_prepare(IndexingMethod *indm, + UnitCell *cell, struct detector *det, + float *ltl) { ERROR("This copy of CrystFEL was compiled without FFTW support.\n"); ERROR("To use asdf indexing, recompile with FFTW.\n"); return NULL; } -void asdf_cleanup(IndexingPrivate *pp) +void asdf_cleanup(void *pp) { } diff --git a/libcrystfel/src/dirax.c b/libcrystfel/src/dirax.c index 19f35696..e9466a24 100644 --- a/libcrystfel/src/dirax.c +++ b/libcrystfel/src/dirax.c @@ -3,11 +3,11 @@ * * Invoke the DirAx auto-indexing program * - * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * * Authors: - * 2010-2014 Thomas White + * 2010-2014,2017 Thomas White * * This file is part of CrystFEL. * @@ -532,7 +532,7 @@ static void write_drx(struct image *image) } -int run_dirax(struct image *image, IndexingPrivate *ipriv) +int run_dirax(struct image *image, void *ipriv) { unsigned int opts; int status; @@ -638,8 +638,8 @@ int run_dirax(struct image *image, IndexingPrivate *ipriv) } -IndexingPrivate *dirax_prepare(IndexingMethod *indm, UnitCell *cell, - struct detector *det, float *ltl) +void *dirax_prepare(IndexingMethod *indm, UnitCell *cell, + struct detector *det, float *ltl) { struct dirax_private *dp; int need_cell = 0; @@ -670,7 +670,7 @@ IndexingPrivate *dirax_prepare(IndexingMethod *indm, UnitCell *cell, } -void dirax_cleanup(IndexingPrivate *pp) +void dirax_cleanup(void *pp) { struct dirax_private *p; p = (struct dirax_private *)pp; diff --git a/libcrystfel/src/dirax.h b/libcrystfel/src/dirax.h index 9776f3f0..96ba6dbc 100644 --- a/libcrystfel/src/dirax.h +++ b/libcrystfel/src/dirax.h @@ -3,11 +3,11 @@ * * Invoke the DirAx auto-indexing program * - * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * * Authors: - * 2010,2012-2014 Thomas White + * 2010,2012-2014,2017 Thomas White * * This file is part of CrystFEL. * @@ -39,13 +39,12 @@ extern "C" { #endif -extern int run_dirax(struct image *image, IndexingPrivate *ipriv); +extern int run_dirax(struct image *image, void *ipriv); -extern IndexingPrivate *dirax_prepare(IndexingMethod *indm, - UnitCell *cell, struct detector *det, - float *ltl); +extern void *dirax_prepare(IndexingMethod *indm, + UnitCell *cell, struct detector *det, float *ltl); -extern void dirax_cleanup(IndexingPrivate *pp); +extern void dirax_cleanup(void *pp); #ifdef __cplusplus } diff --git a/libcrystfel/src/felix.c b/libcrystfel/src/felix.c index 02d523c6..b531e3af 100644 --- a/libcrystfel/src/felix.c +++ b/libcrystfel/src/felix.c @@ -3,8 +3,8 @@ * * Invoke Felix for multi-crystal autoindexing. * - * Copyright © 2015 Deutsches Elektronen-Synchrotron DESY, - * a research centre of the Helmholtz Association. + * Copyright © 2015-2017 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. * * Authors: * 2015 Thomas White @@ -658,9 +658,8 @@ static void parse_options(const char *options, struct felix_private *gp) free(option); } -IndexingPrivate *felix_prepare(IndexingMethod *indm, UnitCell *cell, - struct detector *det, float *ltl, - const char *options) +void *felix_prepare(IndexingMethod *indm, UnitCell *cell, + struct detector *det, float *ltl, const char *options) { struct felix_private *gp; diff --git a/libcrystfel/src/felix.h b/libcrystfel/src/felix.h index 40771933..c06e963a 100644 --- a/libcrystfel/src/felix.h +++ b/libcrystfel/src/felix.h @@ -3,12 +3,12 @@ * * Invoke Felix for multi-crystal autoindexing * - * Copyright © 2013 Deutsches Elektronen-Synchrotron DESY, - * a research centre of the Helmholtz Association. + * Copyright © 2013-2017 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. * * Authors: - * 2010-2013 Thomas White - * 2013-2014 Kenneth Beyerlein + * 2010-2013,2017 Thomas White + * 2013-2014 Kenneth Beyerlein * * This file is part of CrystFEL. * @@ -36,11 +36,9 @@ #include "cell.h" -extern IndexingPrivate *felix_prepare(IndexingMethod *indm, - UnitCell *cell, - struct detector *det, - float *ltl, - const char *options); +extern void *felix_prepare(IndexingMethod *indm, UnitCell *cell, + struct detector *det, float *ltl, + const char *options); extern void felix_cleanup(IndexingPrivate *pp); diff --git a/libcrystfel/src/mosflm.c b/libcrystfel/src/mosflm.c index 05f853eb..a8e6e119 100644 --- a/libcrystfel/src/mosflm.c +++ b/libcrystfel/src/mosflm.c @@ -725,7 +725,7 @@ static int mosflm_readable(struct image *image, struct mosflm_data *mosflm) } -int run_mosflm(struct image *image, IndexingPrivate *ipriv) +int run_mosflm(struct image *image, void *ipriv) { struct mosflm_data *mosflm; unsigned int opts; @@ -843,8 +843,8 @@ int run_mosflm(struct image *image, IndexingPrivate *ipriv) } -IndexingPrivate *mosflm_prepare(IndexingMethod *indm, UnitCell *cell, - struct detector *det, float *ltl) +void *mosflm_prepare(IndexingMethod *indm, UnitCell *cell, + struct detector *det, float *ltl) { struct mosflm_private *mp; int need_cell = 0; @@ -911,7 +911,7 @@ IndexingPrivate *mosflm_prepare(IndexingMethod *indm, UnitCell *cell, } -void mosflm_cleanup(IndexingPrivate *pp) +void mosflm_cleanup(void *pp) { struct mosflm_private *p; p = (struct mosflm_private *)pp; diff --git a/libcrystfel/src/mosflm.h b/libcrystfel/src/mosflm.h index 572741dd..39ec5390 100644 --- a/libcrystfel/src/mosflm.h +++ b/libcrystfel/src/mosflm.h @@ -3,13 +3,13 @@ * * Invoke the DPS auto-indexing algorithm through MOSFLM * - * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * Copyright © 2012 Richard Kirian * * Authors: * 2010 Richard Kirian - * 2012-2014 Thomas White + * 2012-2014,2017 Thomas White * * This file is part of CrystFEL. * @@ -41,12 +41,12 @@ extern "C" { #endif -extern int run_mosflm(struct image *image, IndexingPrivate *ipriv); +extern int run_mosflm(struct image *image, void *ipriv); -extern IndexingPrivate *mosflm_prepare(IndexingMethod *indm, UnitCell *cell, - struct detector *det, float *ltl); +extern void *mosflm_prepare(IndexingMethod *indm, UnitCell *cell, + struct detector *det, float *ltl); -extern void mosflm_cleanup(IndexingPrivate *pp); +extern void mosflm_cleanup(void *pp); #ifdef __cplusplus } diff --git a/libcrystfel/src/xds.c b/libcrystfel/src/xds.c index 15826760..8a7daab7 100644 --- a/libcrystfel/src/xds.c +++ b/libcrystfel/src/xds.c @@ -3,12 +3,12 @@ * * Invoke xds for crystal autoindexing * - * Copyright © 2013-2016 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2013-2017 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * Copyright © 2013 Cornelius Gati * * Authors: - * 2010-2016 Thomas White + * 2010-2017 Thomas White * 2013 Cornelius Gati * * This file is part of CrystFEL. @@ -498,7 +498,7 @@ static int write_inp(struct image *image, struct xds_private *xp) } -int run_xds(struct image *image, IndexingPrivate *priv) +int run_xds(struct image *image, void *priv) { unsigned int opts; int status; @@ -619,8 +619,8 @@ int run_xds(struct image *image, IndexingPrivate *priv) } -IndexingPrivate *xds_prepare(IndexingMethod *indm, UnitCell *cell, - struct detector *det, float *ltl) +void *xds_prepare(IndexingMethod *indm, UnitCell *cell, + struct detector *det, float *ltl) { struct xds_private *xp; int need_cell = 0; @@ -687,11 +687,11 @@ IndexingPrivate *xds_prepare(IndexingMethod *indm, UnitCell *cell, xp->cell = cell; xp->indm = *indm; - return (IndexingPrivate *)xp; + return xp; } -void xds_cleanup(IndexingPrivate *pp) +void xds_cleanup(void *pp) { struct xds_private *xp; diff --git a/libcrystfel/src/xds.h b/libcrystfel/src/xds.h index a0db2054..094d6d2f 100644 --- a/libcrystfel/src/xds.h +++ b/libcrystfel/src/xds.h @@ -3,13 +3,13 @@ * * Invoke xds for crystal autoindexing * - * Copyright © 2013 Deutsches Elektronen-Synchrotron DESY, - * a research centre of the Helmholtz Association. + * Copyright © 2013-2017 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. * Copyright © 2013 Cornelius Gati * * Authors: - * 2010-2013 Thomas White - * 2013 Cornelius Gati + * 2010-2013,2017 Thomas White + * 2013 Cornelius Gati * * This file is part of CrystFEL. * @@ -42,12 +42,12 @@ extern "C" { #endif -extern int run_xds(struct image *image, IndexingPrivate *ipriv); +extern int run_xds(struct image *image, void *ipriv); -extern IndexingPrivate *xds_prepare(IndexingMethod *indm, UnitCell *cell, - struct detector *det, float *ltl); +extern void *xds_prepare(IndexingMethod *indm, UnitCell *cell, + struct detector *det, float *ltl); -extern void xds_cleanup(IndexingPrivate *pp); +extern void xds_cleanup(void *pp); #ifdef __cplusplus } -- cgit v1.2.3