diff options
Diffstat (limited to 'libcrystfel/src')
-rw-r--r-- | libcrystfel/src/cell-utils.c | 97 | ||||
-rw-r--r-- | libcrystfel/src/cell-utils.h | 7 | ||||
-rw-r--r-- | libcrystfel/src/crystfel-mille.c | 30 | ||||
-rw-r--r-- | libcrystfel/src/datatemplate.c | 61 | ||||
-rw-r--r-- | libcrystfel/src/detgeom.c | 1 | ||||
-rw-r--r-- | libcrystfel/src/geometry.c | 3 | ||||
-rw-r--r-- | libcrystfel/src/image-hdf5.c | 77 | ||||
-rw-r--r-- | libcrystfel/src/image.c | 3 | ||||
-rw-r--r-- | libcrystfel/src/index.c | 7 | ||||
-rw-r--r-- | libcrystfel/src/indexers/asdf.c | 2 | ||||
-rw-r--r-- | libcrystfel/src/indexers/pinkindexer.c | 1 | ||||
-rw-r--r-- | libcrystfel/src/indexers/taketwo.c | 82 | ||||
-rw-r--r-- | libcrystfel/src/integration.c | 7 | ||||
-rw-r--r-- | libcrystfel/src/peakfinder8.c | 1 | ||||
-rw-r--r-- | libcrystfel/src/predict-refine.c | 31 | ||||
-rw-r--r-- | libcrystfel/src/utils.c | 26 | ||||
-rw-r--r-- | libcrystfel/src/utils.h | 1 |
17 files changed, 301 insertions, 136 deletions
diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c index 80cdd955..6e5d65b7 100644 --- a/libcrystfel/src/cell-utils.c +++ b/libcrystfel/src/cell-utils.c @@ -3,12 +3,12 @@ * * Unit Cell utility functions * - * Copyright © 2012-2021 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2024 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * Copyright © 2012 Lorenzo Galli * * Authors: - * 2009-2021 Thomas White <taw@physics.org> + * 2009-2024 Thomas White <taw@physics.org> * 2012 Lorenzo Galli * * This file is part of CrystFEL. @@ -42,6 +42,7 @@ #include "cell-utils.h" #include "utils.h" #include "image.h" +#include "symmetry.h" /** @@ -2245,9 +2246,9 @@ static IntegerMatrix *check_permutations(UnitCell *cell_reduced, UnitCell *refer best_diff = diff; n_best = 0; - } + best_m[n_best++] = intmat_copy(m); - if ( diff < 1.001*best_diff ) { + } else if ( diff < 1.001*best_diff ) { /* If the new solution is the same as the * previous one, add it to the list */ @@ -2414,3 +2415,91 @@ UnitCell *compare_reindexed_cell_parameters(UnitCell *cell_in, return match; } + + +static char *add_unique_axis(const char *inp, char ua) +{ + char *pg = cfmalloc(64); + if ( pg == NULL ) return NULL; + snprintf(pg, 63, "%s_ua%c", inp, ua); + return pg; +} + + +static char *get_chiral_holohedry(UnitCell *cell) +{ + LatticeType lattice = cell_get_lattice_type(cell); + char *pg; + int add_ua = 1; + + switch (lattice) + { + case L_TRICLINIC: + pg = "1"; + add_ua = 0; + break; + + case L_MONOCLINIC: + pg = "2"; + break; + + case L_ORTHORHOMBIC: + pg = "222"; + add_ua = 0; + break; + + case L_TETRAGONAL: + pg = "422"; + break; + + case L_RHOMBOHEDRAL: + pg = "3_R"; + add_ua = 0; + break; + + case L_HEXAGONAL: + if ( cell_get_centering(cell) == 'H' ) { + pg = "3_H"; + add_ua = 0; + } else { + pg = "622"; + } + break; + + case L_CUBIC: + pg = "432"; + add_ua = 0; + break; + + default: + pg = "error"; + break; + } + + if ( add_ua ) { + return add_unique_axis(pg, cell_get_unique_axis(cell)); + } else { + return cfstrdup(pg); + } +} + + +/** + * \param cell: A UnitCell + * + * Determines the chiral point group of the lattice, based on the lattice type + * and unique axis of \p cell (the lattice parameters are not checked). + * + * This point group is often needed within indexing algorithms. + * + * \returns a newly allocated SymOpList. + * + */ +SymOpList *get_lattice_symmetry(UnitCell *cell) +{ + SymOpList *rawList; + char *pg = get_chiral_holohedry(cell); + rawList = get_pointgroup(pg); + cffree(pg); + return rawList; +} diff --git a/libcrystfel/src/cell-utils.h b/libcrystfel/src/cell-utils.h index b473a4c5..de200bcb 100644 --- a/libcrystfel/src/cell-utils.h +++ b/libcrystfel/src/cell-utils.h @@ -3,12 +3,12 @@ * * Unit Cell utility functions * - * Copyright © 2012-2021 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2024 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * Copyright © 2012 Lorenzo Galli * * Authors: - * 2009-2020 Thomas White <taw@physics.org> + * 2009-2024 Thomas White <taw@physics.org> * 2012 Lorenzo Galli * * This file is part of CrystFEL. @@ -34,6 +34,7 @@ #include <gsl/gsl_matrix.h> #include "cell.h" +#include "symmetry.h" #ifdef __cplusplus extern "C" { @@ -102,6 +103,8 @@ extern UnitCell *compare_reindexed_cell_parameters(UnitCell *cell_in, const double *tols, RationalMatrix **pmb); +extern SymOpList *get_lattice_symmetry(UnitCell *cell); + #ifdef __cplusplus } #endif diff --git a/libcrystfel/src/crystfel-mille.c b/libcrystfel/src/crystfel-mille.c index 72f89425..c2fe2a57 100644 --- a/libcrystfel/src/crystfel-mille.c +++ b/libcrystfel/src/crystfel-mille.c @@ -72,6 +72,10 @@ struct mille int *int_arr; int max_entries; int n; + + int *have_local; + int n_local; + FILE *fh; }; @@ -115,6 +119,17 @@ static void mille_add_measurement(Mille *m, m->max_entries = new_max_entries; } + if ( NLC > m->n_local ) { + int i; + int *new_have_local = cfrealloc(m->have_local, NLC*sizeof(int)); + if ( new_have_local == NULL ) return; + m->have_local = new_have_local; + for ( i=m->n_local; i<NLC; i++ ) { + m->have_local[i] = 0; + } + m->n_local = NLC; + } + /* The measurement */ m->float_arr[m->n] = rMeas; m->int_arr[m->n] = 0; @@ -126,6 +141,7 @@ static void mille_add_measurement(Mille *m, m->float_arr[m->n] = derLc[i]; m->int_arr[m->n] = i+1; m->n++; + m->have_local[i] = 1; } } @@ -292,6 +308,8 @@ Mille *crystfel_mille_new(const char *outFileName) m->n = 0; m->float_arr = NULL; m->int_arr = NULL; + m->have_local = NULL; + m->n_local = 0; m->fh = fopen(outFileName, "wb"); if ( m->fh == NULL ) { @@ -318,11 +336,15 @@ void crystfel_mille_free(Mille *m) void crystfel_mille_delete_last_record(Mille *m) { m->n = 0; + cffree(m->have_local); + m->have_local = NULL; + m->n_local = 0; } void crystfel_mille_write_record(Mille *m) { + int i; float nf = 0.0; int ni = 0; int nw = (m->n * 2)+2; @@ -330,6 +352,14 @@ void crystfel_mille_write_record(Mille *m) /* Don't write empty records */ if ( m->n == 0 ) return; + /* Don't write records with incomplete local gradients */ + for ( i=0; i<m->n_local; i++ ) { + if ( !m->have_local[i] ) { + crystfel_mille_delete_last_record(m); + return; + } + } + fwrite(&nw, sizeof(int), 1, m->fh); fwrite(&nf, sizeof(float), 1, m->fh); diff --git a/libcrystfel/src/datatemplate.c b/libcrystfel/src/datatemplate.c index a7b8950f..c97bd6d7 100644 --- a/libcrystfel/src/datatemplate.c +++ b/libcrystfel/src/datatemplate.c @@ -3,11 +3,11 @@ * * Data template structure * - * Copyright © 2019-2021 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2019-2024 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * * Authors: - * 2019-2021 Thomas White <taw@physics.org> + * 2019-2024 Thomas White <taw@physics.org> * * This file is part of CrystFEL. * @@ -830,12 +830,35 @@ static int parse_peak_layout(const char *val, return 1; } +#define MAX_FOR_LATER (1024) + +struct forlater +{ + char *keys[MAX_FOR_LATER]; + char *vals[MAX_FOR_LATER]; + int n_forlater; +}; + + +static void store_for_later(struct forlater *fl, const char *key, const char *val) +{ + if ( fl->n_forlater >= MAX_FOR_LATER ) { + ERROR("Too many lines stored.\n"); + return; + } + + fl->keys[fl->n_forlater] = strdup(key); + fl->vals[fl->n_forlater] = strdup(val); + fl->n_forlater++; +} + static int parse_toplevel(DataTemplate *dt, const char *key, const char *val, struct panel_template *defaults, - int *defaults_updated) + int *defaults_updated, + struct forlater *for_later) { if ( strcmp(key, "detector_shift_x") == 0 ) { dt->shift_x_from = cfstrdup(val); @@ -883,8 +906,12 @@ static int parse_toplevel(DataTemplate *dt, } else if ( strncmp(key, "group_", 6) == 0 ) { - if ( parse_group(key+6, dt, val) ) { - return 1; + if ( for_later != NULL ) { + store_for_later(for_later, key, val); + } else { + if ( parse_group(key+6, dt, val) ) { + return 1; + } } } else { @@ -1042,6 +1069,9 @@ DataTemplate *data_template_new_from_string(const char *string_in) size_t len; struct panel_template defaults; int have_unused_defaults = 0; + struct forlater lines_for_later; + + lines_for_later.n_forlater = 0; dt = cfcalloc(1, sizeof(DataTemplate)); if ( dt == NULL ) return NULL; @@ -1179,10 +1209,10 @@ DataTemplate *data_template_new_from_string(const char *string_in) /* Top-level option */ if ( parse_toplevel(dt, line, val, &defaults, - &have_unused_defaults) ) + &have_unused_defaults, + &lines_for_later) ) { - ERROR("Invalid top-level line '%s'\n", - line); + ERROR("Invalid top-level line '%s'\n", line); reject = 1; } cffree(line); @@ -1217,6 +1247,21 @@ DataTemplate *data_template_new_from_string(const char *string_in) } while ( !done ); + for ( i=0; i<lines_for_later.n_forlater; i++ ) { + char *key = lines_for_later.keys[i]; + char *val = lines_for_later.vals[i]; + if ( parse_toplevel(dt, key, val, + &defaults, + &have_unused_defaults, + NULL) ) + { + ERROR("Invalid top-level line '%s' = '%s'\n", key, val); + reject = 1; + } + cffree(key); + cffree(val); + } + if ( dt->n_panels == 0 ) { ERROR("No panel descriptions in geometry file.\n"); cffree(dt); diff --git a/libcrystfel/src/detgeom.c b/libcrystfel/src/detgeom.c index 39b00663..2b38664c 100644 --- a/libcrystfel/src/detgeom.c +++ b/libcrystfel/src/detgeom.c @@ -269,6 +269,7 @@ gsl_matrix **make_panel_minvs(struct detgeom *dg) gsl_matrix_set(M, 2, 2, p->pixel_pitch*p->ssz); Minvs[i] = matrix_invert(M); + gsl_matrix_free(M); if ( Minvs[i] == NULL ) { ERROR("Failed to calculate inverse panel matrix for %s\n", p->name); diff --git a/libcrystfel/src/geometry.c b/libcrystfel/src/geometry.c index 10a7fa19..7bcb8d8f 100644 --- a/libcrystfel/src/geometry.c +++ b/libcrystfel/src/geometry.c @@ -306,7 +306,7 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst, z -= g.kcen; /* Three because the general case fails in extreme cases */ - if ( w0 / w1 <= DBL_MIN ) { + if ( isinf(w0) || (w0 / w1 <= DBL_MIN) ) { /* 'Laue' corner case */ kpred = g.kcen; @@ -352,6 +352,7 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst, /* Revert the 'Lorentz' factor */ partiality *= sqrt( ( R*R + M2_k/sumw_k) / ( R*R ) ); + if ( isnan(partiality) ) partiality = 0.0; if ( (updateme == NULL) && ( partiality < min_partiality ) ) return NULL; diff --git a/libcrystfel/src/image-hdf5.c b/libcrystfel/src/image-hdf5.c index a148c307..93ea2d8d 100644 --- a/libcrystfel/src/image-hdf5.c +++ b/libcrystfel/src/image-hdf5.c @@ -414,6 +414,7 @@ static int load_hdf5_hyperslab(struct panel_template *p, if ( ndims < 0 ) { ERROR("Failed to get number of dimensions for panel %s\n", p->name); + H5Sclose(dataspace); H5Dclose(dh); return 1; } @@ -436,6 +437,7 @@ static int load_hdf5_hyperslab(struct panel_template *p, p->name, ndims, total_dt_dims, total_dt_dims - plh_dt_dims); H5Dclose(dh); + H5Sclose(dataspace); return 1; } } else { @@ -453,6 +455,7 @@ static int load_hdf5_hyperslab(struct panel_template *p, cffree(f_offset); cffree(f_count); H5Dclose(dh); + H5Sclose(dataspace); return 1; } @@ -502,6 +505,7 @@ static int load_hdf5_hyperslab(struct panel_template *p, cffree(f_offset); cffree(f_count); H5Dclose(dh); + H5Sclose(dataspace); return 1; } @@ -511,6 +515,8 @@ static int load_hdf5_hyperslab(struct panel_template *p, profile_start("H5Dread"); r = H5Dread(dh, el_type, memspace, dataspace, H5P_DEFAULT, data); + H5Sclose(memspace); + H5Sclose(dataspace); profile_end("H5Dread"); if ( r < 0 ) { ERROR("Couldn't read data for panel %s\n", @@ -538,6 +544,7 @@ static int load_hdf5_hyperslab(struct panel_template *p, static hid_t open_hdf5_file(const char *filename) { hid_t fh; + hid_t fapl; if ( access(filename, R_OK) == -1 ) { ERROR("File does not exist or cannot be read: %s\n", @@ -545,7 +552,10 @@ static hid_t open_hdf5_file(const char *filename) return -1; } - fh = H5Fopen(filename, H5F_ACC_RDONLY, H5P_DEFAULT); + fapl = H5Pcreate(H5P_FILE_ACCESS); + H5Pset_fclose_degree(fapl, H5F_CLOSE_STRONG); + fh = H5Fopen(filename, H5F_ACC_RDONLY, fapl); + H5Pclose(fapl); if ( fh < 0 ) { ERROR("Couldn't open HDF5 file: %s\n", filename); return -1; @@ -650,6 +660,7 @@ int image_hdf5_read_satmap(struct panel_template *p, sizeof(float), 1, map_location, NULL) ) { ERROR("Failed to load saturation map data\n"); + close_hdf5(fh); return 1; } @@ -687,6 +698,7 @@ int image_hdf5_read_mask(struct panel_template *p, sizeof(int), 1, mask_location, NULL) ) { ERROR("Failed to load mask data\n"); + close_hdf5(fh); cffree(mask); return 1; } @@ -718,6 +730,7 @@ static char *read_single_fixed_string(hid_t dh) sh = H5Dget_space(dh); if ( H5Sget_simple_extent_ndims(sh) ) { ERROR("Non-scalar string\n"); + H5Sclose(sh); return NULL; } @@ -727,6 +740,7 @@ static char *read_single_fixed_string(hid_t dh) tmp = cfmalloc(size+1); if ( tmp == NULL ) { H5Tclose(type); + H5Sclose(sh); return NULL; } r = H5Dread(dh, type, sh, H5S_ALL, H5P_DEFAULT, tmp); @@ -832,6 +846,7 @@ int image_hdf5_read_header_to_cache(struct image *image, const char *name) ndims = H5Sget_simple_extent_ndims(sh); if ( ndims > 64 ) { ERROR("Too many dimensions for numeric value\n"); + H5Sclose(sh); close_hdf5(fh); cffree(subst_name); return 1; @@ -856,11 +871,15 @@ int image_hdf5_read_header_to_cache(struct image *image, const char *name) ERROR("Couldn't read scalar value from %s.\n", subst_name); cffree(subst_name); + H5Sclose(sh); + H5Sclose(ms); close_hdf5(fh); return 1; } image_cache_header_float(image, name, val); cffree(subst_name); + H5Sclose(sh); + H5Sclose(ms); return 0; } else if ( class == H5T_INTEGER ) { @@ -868,6 +887,8 @@ int image_hdf5_read_header_to_cache(struct image *image, const char *name) int val; r = H5Dread(dh, H5T_NATIVE_INT, ms, sh, H5P_DEFAULT, &val); + H5Sclose(sh); + H5Sclose(ms); if ( r < 0 ) { ERROR("Couldn't read scalar value from %s.\n", subst_name); @@ -912,12 +933,16 @@ int image_hdf5_read_header_to_cache(struct image *image, const char *name) cffree(subst_name); close_hdf5(fh); + H5Sclose(sh); + H5Sclose(ms); return rv; } else { /* Should never be reached */ ERROR("Invalid HDF5 class %i\n", class); cffree(subst_name); + H5Sclose(sh); + H5Sclose(ms); return 1; } } @@ -927,6 +952,8 @@ int image_hdf5_read_header_to_cache(struct image *image, const char *name) ERROR("Couldn't parse event '%s'\n"); close_hdf5(fh); cffree(subst_name); + H5Sclose(sh); + H5Sclose(ms); return 1; } @@ -936,6 +963,8 @@ int image_hdf5_read_header_to_cache(struct image *image, const char *name) ERROR("Couldn't allocate dimension arrays\n"); close_hdf5(fh); cffree(subst_name); + H5Sclose(sh); + H5Sclose(ms); return 1; } @@ -953,6 +982,8 @@ int image_hdf5_read_header_to_cache(struct image *image, const char *name) subst_name, i, dim_vals[dim_val_pos], size[i]); close_hdf5(fh); + H5Sclose(sh); + H5Sclose(ms); cffree(subst_name); cffree(dim_vals); return 1; @@ -979,6 +1010,8 @@ int image_hdf5_read_header_to_cache(struct image *image, const char *name) cffree(f_offset); cffree(f_count); cffree(subst_name); + H5Sclose(sh); + H5Sclose(ms); close_hdf5(fh); return 1; } @@ -986,12 +1019,13 @@ int image_hdf5_read_header_to_cache(struct image *image, const char *name) cffree(f_offset); cffree(f_count); - ms = H5Screate_simple(1,msdims,NULL); check = H5Sselect_hyperslab(ms, H5S_SELECT_SET, m_offset, NULL, m_count, NULL); if ( check < 0 ) { ERROR("Error selecting memory dataspace for header value\n"); close_hdf5(fh); + H5Sclose(sh); + H5Sclose(ms); cffree(subst_name); return 1; } @@ -1000,6 +1034,8 @@ int image_hdf5_read_header_to_cache(struct image *image, const char *name) double val; r = H5Dread(dh, H5T_NATIVE_DOUBLE, ms, sh, H5P_DEFAULT, &val); + H5Sclose(sh); + H5Sclose(ms); if ( r < 0 ) { ERROR("Couldn't read value.\n"); close_hdf5(fh); @@ -1016,6 +1052,8 @@ int image_hdf5_read_header_to_cache(struct image *image, const char *name) int val; r = H5Dread(dh, H5T_NATIVE_INT, ms, sh, H5P_DEFAULT, &val); + H5Sclose(sh); + H5Sclose(ms); if ( r < 0 ) { ERROR("Couldn't read value.\n"); close_hdf5(fh); @@ -1041,6 +1079,8 @@ int image_hdf5_read_header_to_cache(struct image *image, const char *name) char *val; rv = H5Dread(dh, stype, ms, sh, H5P_DEFAULT, &val); + H5Sclose(sh); + H5Sclose(ms); if ( rv < 0 ) { ERROR("Can't read HDF5 vlen string from array - %s\n", subst_name); @@ -1069,10 +1109,14 @@ int image_hdf5_read_header_to_cache(struct image *image, const char *name) val = cfmalloc(ssize+1); if ( val == NULL ) { close_hdf5(fh); + H5Sclose(ms); + H5Sclose(sh); cffree(subst_name); return 1; } rv = H5Dread(dh, stype, ms, sh, H5P_DEFAULT, val); + H5Sclose(sh); + H5Sclose(ms); H5Tclose(stype); if ( rv < 0 ) { ERROR("Couldn't read HDF5 fixed string from array - %s\n", @@ -1097,6 +1141,8 @@ int image_hdf5_read_header_to_cache(struct image *image, const char *name) } else { /* Should never be reached */ ERROR("Invalid HDF5 class %i\n", class); + H5Sclose(sh); + H5Sclose(ms); close_hdf5(fh); cffree(subst_name); return 1; @@ -1303,7 +1349,7 @@ ImageFeatureList *image_hdf5_read_peaks_cxi(const DataTemplate *dtempl, int half_pixel_shift) { ImageFeatureList *features; - hid_t fh; + hid_t fh, fapl; char path_n[1024]; char path_x[1024]; char path_y[1024]; @@ -1355,7 +1401,10 @@ ImageFeatureList *image_hdf5_read_peaks_cxi(const DataTemplate *dtempl, snprintf(path_y, 1024, "%s/peakYPosRaw", subst_name); snprintf(path_i, 1024, "%s/peakTotalIntensity", subst_name); - fh = H5Fopen(filename, H5F_ACC_RDONLY, H5P_DEFAULT); + fapl = H5Pcreate(H5P_FILE_ACCESS); + H5Pset_fclose_degree(fapl, H5F_CLOSE_STRONG); + fh = H5Fopen(filename, H5F_ACC_RDONLY, fapl); + H5Pclose(fapl); if ( fh < 0 ) { ERROR("Couldn't open file (peaks/cxi): %s\n", filename); cffree(subst_name); @@ -1431,7 +1480,7 @@ ImageFeatureList *image_hdf5_read_peaks_hdf5(const DataTemplate *dtempl, const char *event, int half_pixel_shift) { - hid_t fh, dh, sh; + hid_t fh, dh, sh, fapl; hsize_t size[2]; hsize_t max_size[2]; int i; @@ -1453,7 +1502,10 @@ ImageFeatureList *image_hdf5_read_peaks_hdf5(const DataTemplate *dtempl, return NULL; } - fh = H5Fopen(filename, H5F_ACC_RDONLY, H5P_DEFAULT); + fapl = H5Pcreate(H5P_FILE_ACCESS); + H5Pset_fclose_degree(fapl, H5F_CLOSE_STRONG); + fh = H5Fopen(filename, H5F_ACC_RDONLY, fapl); + H5Pclose(fapl); if ( fh < 0 ) { ERROR("Couldn't open file (peaks/hdf5): %s\n", filename); return NULL; @@ -1486,6 +1538,7 @@ ImageFeatureList *image_hdf5_read_peaks_hdf5(const DataTemplate *dtempl, if ( H5Sget_simple_extent_ndims(sh) != 2 ) { ERROR("Peak list has the wrong dimensionality (%i).\n", H5Sget_simple_extent_ndims(sh)); + H5Sclose(sh); close_hdf5(fh); return NULL; } @@ -1860,7 +1913,7 @@ char **image_hdf5_expand_frames(const DataTemplate *dtempl, { char **path_evs; int n_path_evs; - hid_t fh; + hid_t fh, fapl; int i; int dims_expected; struct ev_list full_evs; @@ -1886,7 +1939,10 @@ char **image_hdf5_expand_frames(const DataTemplate *dtempl, return NULL; } - fh = H5Fopen(filename, H5F_ACC_RDONLY, H5P_DEFAULT); + fapl = H5Pcreate(H5P_FILE_ACCESS); + H5Pset_fclose_degree(fapl, H5F_CLOSE_STRONG); + fh = H5Fopen(filename, H5F_ACC_RDONLY, fapl); + H5Pclose(fapl); if ( fh < 0 ) { ERROR("Couldn't open file (hdf5_expand_frames): %s\n", filename); return NULL; @@ -1950,6 +2006,7 @@ char **image_hdf5_expand_frames(const DataTemplate *dtempl, "(%s has %i, expected %i)\n", path, dims, dims_expected); close_hdf5(fh); + H5Sclose(sh); return NULL; } @@ -1958,15 +2015,19 @@ char **image_hdf5_expand_frames(const DataTemplate *dtempl, if ( (size == NULL) || (placeholder_sizes == NULL) ) { ERROR("Failed to allocate dimensions\n"); close_hdf5(fh); + H5Sclose(sh); return NULL; } if ( H5Sget_simple_extent_dims(sh, size, NULL) < 0 ) { ERROR("Failed to get size\n"); close_hdf5(fh); + H5Sclose(sh); return NULL; } + H5Sclose(sh); + n_placeholder_dims = 0; for ( j=0; j<dims; j++ ) { if ( p->dims[j] == DIM_PLACEHOLDER ) { diff --git a/libcrystfel/src/image.c b/libcrystfel/src/image.c index e91c09d6..3b215c08 100644 --- a/libcrystfel/src/image.c +++ b/libcrystfel/src/image.c @@ -1431,6 +1431,9 @@ void image_free(struct image *image) for ( i=0; i<image->n_cached_headers; i++ ) { cffree(image->header_cache[i]->header_name); + if ( image->header_cache[i]->type == HEADER_STR ) { + cffree(image->header_cache[i]->val_str); + } cffree(image->header_cache[i]); } diff --git a/libcrystfel/src/index.c b/libcrystfel/src/index.c index 12cd9190..5b0f742f 100644 --- a/libcrystfel/src/index.c +++ b/libcrystfel/src/index.c @@ -50,6 +50,7 @@ #include "cell-utils.h" #include "predict-refine.h" #include "profile.h" +#include "crystfel-mille.h" #include "indexers/dirax.h" #include "indexers/asdf.h" #include "indexers/mosflm.h" @@ -774,6 +775,12 @@ static int try_indexer(struct image *image, IndexingMethod indm, } profile_end("cell-compare-to-others"); + if ( (mille != NULL) && !crystal_get_user_flag(cr) ) { + profile_start("mille-write"); + crystfel_mille_write_record(mille); + profile_end("mille-write"); + } + } n_bad = remove_flagged_crystals(image); diff --git a/libcrystfel/src/indexers/asdf.c b/libcrystfel/src/indexers/asdf.c index c136efc1..6e1e5781 100644 --- a/libcrystfel/src/indexers/asdf.c +++ b/libcrystfel/src/indexers/asdf.c @@ -1308,7 +1308,7 @@ int run_asdf(struct image *image, void *ipriv) } -void *asdf_prepare(IndexingMethod *indm, UnitCell *cell) +void *asdf_prepare(IndexingMethod *indm, UnitCell *cell, struct asdf_options *asdf_opts) { ERROR("This copy of CrystFEL was compiled without FFTW support.\n"); ERROR("To use asdf indexing, recompile with FFTW.\n"); diff --git a/libcrystfel/src/indexers/pinkindexer.c b/libcrystfel/src/indexers/pinkindexer.c index 180246c2..14e18548 100644 --- a/libcrystfel/src/indexers/pinkindexer.c +++ b/libcrystfel/src/indexers/pinkindexer.c @@ -33,7 +33,6 @@ #include <stdlib.h> -#include <sys/errno.h> #include <argp.h> #include "utils.h" diff --git a/libcrystfel/src/indexers/taketwo.c b/libcrystfel/src/indexers/taketwo.c index c87d6b3d..02fa9a23 100644 --- a/libcrystfel/src/indexers/taketwo.c +++ b/libcrystfel/src/indexers/taketwo.c @@ -522,84 +522,6 @@ static double matrix_trace(gsl_matrix *a) return tr; } -static char *add_unique_axis(const char *inp, char ua) -{ - char *pg = cfmalloc(64); - if ( pg == NULL ) return NULL; - snprintf(pg, 63, "%s_ua%c", inp, ua); - return pg; -} - - -static char *get_chiral_holohedry(UnitCell *cell) -{ - LatticeType lattice = cell_get_lattice_type(cell); - char *pg; - int add_ua = 1; - - switch (lattice) - { - case L_TRICLINIC: - pg = "1"; - add_ua = 0; - break; - - case L_MONOCLINIC: - pg = "2"; - break; - - case L_ORTHORHOMBIC: - pg = "222"; - add_ua = 0; - break; - - case L_TETRAGONAL: - pg = "422"; - break; - - case L_RHOMBOHEDRAL: - pg = "3_R"; - add_ua = 0; - break; - - case L_HEXAGONAL: - if ( cell_get_centering(cell) == 'H' ) { - pg = "3_H"; - add_ua = 0; - } else { - pg = "622"; - } - break; - - case L_CUBIC: - pg = "432"; - add_ua = 0; - break; - - default: - pg = "error"; - break; - } - - if ( add_ua ) { - return add_unique_axis(pg, cell_get_unique_axis(cell)); - } else { - return cfstrdup(pg); - } -} - - -static SymOpList *sym_ops_for_cell(UnitCell *cell) -{ - SymOpList *rawList; - - char *pg = get_chiral_holohedry(cell); - rawList = get_pointgroup(pg); - cffree(pg); - - return rawList; -} - static int rot_mats_are_similar(gsl_matrix *rot1, gsl_matrix *rot2, gsl_matrix *sub, gsl_matrix *mul, double *score, struct TakeTwoCell *cell) @@ -1619,7 +1541,7 @@ static void set_gsl_matrix(gsl_matrix *mat, static int generate_rotation_sym_ops(struct TakeTwoCell *ttCell) { - SymOpList *rawList = sym_ops_for_cell(ttCell->cell); + SymOpList *rawList = get_lattice_symmetry(ttCell->cell); /* Now we must convert these into rotation matrices rather than hkl * transformations (affects triclinic, monoclinic, rhombohedral and @@ -1858,7 +1780,7 @@ static int gen_theoretical_vecs(UnitCell *cell, struct TheoryVec **cell_vecs, &bsx, &bsy, &bsz, &csx, &csy, &csz); - SymOpList *rawList = sym_ops_for_cell(cell); + SymOpList *rawList = get_lattice_symmetry(cell); cell_get_parameters(cell, &a, &b, &c, &alpha, &beta, &gamma); diff --git a/libcrystfel/src/integration.c b/libcrystfel/src/integration.c index d2c06db7..9ad22960 100644 --- a/libcrystfel/src/integration.c +++ b/libcrystfel/src/integration.c @@ -1467,11 +1467,6 @@ int integrate_rings_once(Reflection *refl, * appropriate size */ if ( intensity < -aduph ) { sig2_poisson = -aduph*intensity; - } else if ( intensity < 0.0 ) { - /* If the intensity is negative (by less than one - * photon), assume the reflection is very weak and - * therefore has a Poisson error of one photon. */ - sig2_poisson = aduph; } sigma = sqrt(sig2_poisson + bx->m*sig2_bg); @@ -1816,6 +1811,8 @@ IntegrationMethod integration_method(const char *str, int *err) meth = INTEGRATION_DEFAULTS_PROF2D; } else if ( strcmp(methods[i], "none") == 0 ) { + cffree(methods[i]); + cffree(methods); return INTEGRATION_NONE; } else if ( strcmp(methods[i], "sat") == 0 ) { diff --git a/libcrystfel/src/peakfinder8.c b/libcrystfel/src/peakfinder8.c index 1def35d3..28bee8d7 100644 --- a/libcrystfel/src/peakfinder8.c +++ b/libcrystfel/src/peakfinder8.c @@ -141,7 +141,6 @@ static struct radial_stats_pixels *compute_rstats_pixels(struct radius_maps *rma cffree(rsp); return NULL; } - srand(0); int n_pixels_per_bin = 100; // Can make this a parameter diff --git a/libcrystfel/src/predict-refine.c b/libcrystfel/src/predict-refine.c index 6da0d3d1..083d86c5 100644 --- a/libcrystfel/src/predict-refine.c +++ b/libcrystfel/src/predict-refine.c @@ -142,15 +142,15 @@ double r_gradient(int param, Reflection *refl, UnitCell *cell, double wavelength /* Spot position gradients for diffraction physics (anything that changes the * diffracted ray direction) */ int fs_ss_gradient_physics(int param, Reflection *refl, UnitCell *cell, - struct detgeom_panel *p, gsl_matrix *Minv, - double fs, double ss, double mu, - float *fsg, float *ssg) + struct detgeom_panel *p, gsl_matrix *Minv, + double fs, double ss, double mu, + float *fsg, float *ssg) { signed int h, k, l; gsl_vector *dRdp; gsl_vector *v; - get_indices(refl, &h, &k, &l); + get_symmetric_indices(refl, &h, &k, &l); dRdp = gsl_vector_calloc(3); @@ -302,7 +302,7 @@ int fs_ss_gradient(int param, Reflection *refl, UnitCell *cell, double mu; double fs, ss; - get_indices(refl, &h, &k, &l); + get_symmetric_indices(refl, &h, &k, &l); kpred = get_kpred(refl); cell_get_reciprocal(cell, &asx, &asy, &asz, &bsx, &bsy, &bsz, @@ -811,7 +811,7 @@ int refine_prediction(struct image *image, Crystal *cr, reflist = reflist_new(); n = pair_peaks(image, cr, reflist, rps); - if ( n < 10 ) { + if ( n < 3 ) { cffree(rps); reflist_free(reflist); return 1; @@ -871,7 +871,8 @@ int refine_prediction(struct image *image, Crystal *cr, total_shifts[0], total_shifts[1], total_shifts[2]); crystal_add_notes(cr, tmp); - if ( mille != NULL ) { + if ( (mille != NULL) && (n>4) ) { + crystfel_mille_delete_last_record(mille); profile_start("mille-calc"); write_mille(mille, n, crystal_get_cell(cr), rps, image, max_mille_level, Minvs); @@ -887,18 +888,10 @@ int refine_prediction(struct image *image, Crystal *cr, n = pair_peaks(image, cr, NULL, rps); free_rps_noreflist(rps, n); - if ( n < 10 ) { - if ( mille != NULL ) { - crystfel_mille_delete_last_record(mille); - } - return 1; - } - if ( mille != NULL ) { - profile_start("mille-write"); - crystfel_mille_write_record(mille); - profile_end("mille-write"); + if ( n < 3 ) { + return 1; + } else { + return 0; } - - return 0; } diff --git a/libcrystfel/src/utils.c b/libcrystfel/src/utils.c index 7bfac487..a35ba272 100644 --- a/libcrystfel/src/utils.c +++ b/libcrystfel/src/utils.c @@ -135,39 +135,53 @@ gsl_matrix *matrix_invert(gsl_matrix *m) perm = gsl_permutation_alloc(m->size1); if ( perm == NULL ) { ERROR("Couldn't allocate permutation\n"); - gsl_matrix_free(m); return NULL; } inv = gsl_matrix_alloc(m->size1, m->size2); if ( inv == NULL ) { ERROR("Couldn't allocate inverse\n"); - gsl_matrix_free(m); gsl_permutation_free(perm); return NULL; } if ( gsl_linalg_LU_decomp(m, perm, &s) ) { ERROR("Couldn't decompose matrix\n"); - gsl_matrix_free(m); gsl_permutation_free(perm); return NULL; } if ( gsl_linalg_LU_invert(m, perm, inv) ) { - ERROR("Couldn't invert cell matrix:\n"); - gsl_matrix_free(m); + ERROR("Couldn't invert matrix:\n"); gsl_permutation_free(perm); return NULL; } gsl_permutation_free(perm); - gsl_matrix_free(m); return inv; } +gsl_vector *solve_inv(gsl_vector *v, gsl_matrix *M) +{ + gsl_matrix *Minv; + gsl_vector *s; + int n; + + n = v->size; + if ( v->size != M->size1 ) return NULL; + if ( v->size != M->size2 ) return NULL; + + Minv = matrix_invert(M); + if ( Minv == NULL ) return NULL; + + s = gsl_vector_calloc(n); + gsl_blas_dgemv(CblasNoTrans, 1.0, M, v, 0.0, s); + return s; +} + + static int check_eigen(gsl_vector *e_val, int verbose) { int i; diff --git a/libcrystfel/src/utils.h b/libcrystfel/src/utils.h index 30a7c4ef..08e059a7 100644 --- a/libcrystfel/src/utils.h +++ b/libcrystfel/src/utils.h @@ -75,6 +75,7 @@ extern "C" { extern void show_matrix_eqn(gsl_matrix *M, gsl_vector *v); extern void show_matrix(gsl_matrix *M); extern void show_vector(gsl_vector *M); +extern gsl_vector *solve_inv(gsl_vector *v, gsl_matrix *M); extern gsl_vector *solve_svd(gsl_vector *v, gsl_matrix *M, int *n_filt, int verbose); extern gsl_matrix *matrix_mult2(gsl_matrix *A, gsl_matrix *B); |