diff options
author | Thomas White <taw@physics.org> | 2017-09-15 15:46:37 +0200 |
---|---|---|
committer | Thomas White <taw@physics.org> | 2017-09-15 15:46:37 +0200 |
commit | 5ca956f67456831501cf413259790eec6d6c6c1f (patch) | |
tree | 7ae69fb902f78af2828d4bbd4609a934c0a09dc7 | |
parent | f6de9f595620b5d4a1ecf3d8d6d4cde9b3c179e0 (diff) | |
parent | ca0d7a9e1982b55ae3891d8dbe7f05a79ed3e7f1 (diff) |
Merge branch 'tom/index'
-rw-r--r-- | libcrystfel/src/asdf.c | 42 | ||||
-rw-r--r-- | libcrystfel/src/cell-utils.c | 121 | ||||
-rw-r--r-- | libcrystfel/src/cell-utils.h | 9 | ||||
-rw-r--r-- | libcrystfel/src/dirax.c | 31 | ||||
-rw-r--r-- | libcrystfel/src/felix.c | 16 | ||||
-rw-r--r-- | libcrystfel/src/hdf5-file.c | 6 | ||||
-rw-r--r-- | libcrystfel/src/image.c | 5 | ||||
-rw-r--r-- | libcrystfel/src/image.h | 2 | ||||
-rw-r--r-- | libcrystfel/src/index.c | 276 | ||||
-rw-r--r-- | libcrystfel/src/index.h | 77 | ||||
-rw-r--r-- | libcrystfel/src/mosflm.c | 63 | ||||
-rw-r--r-- | libcrystfel/src/predict-refine.c | 17 | ||||
-rw-r--r-- | libcrystfel/src/taketwo.c | 19 | ||||
-rw-r--r-- | libcrystfel/src/xds.c | 52 | ||||
-rw-r--r-- | src/indexamajig.c | 38 | ||||
-rw-r--r-- | src/whirligig.c | 111 |
16 files changed, 396 insertions, 489 deletions
diff --git a/libcrystfel/src/asdf.c b/libcrystfel/src/asdf.c index 0c43fe7a..bbfd1a56 100644 --- a/libcrystfel/src/asdf.c +++ b/libcrystfel/src/asdf.c @@ -272,39 +272,14 @@ static int calc_reciprocal(gsl_vector **direct, gsl_vector **reciprocal) static int check_cell(struct asdf_private *dp, struct image *image, UnitCell *cell) { - UnitCell *out; Crystal *cr; - if ( dp->indm & INDEXING_CHECK_CELL_COMBINATIONS ) { - - out = match_cell(cell, dp->template, 0, dp->ltl, 1); - if ( out == NULL ) return 0; - - } else if ( dp->indm & INDEXING_CHECK_CELL_AXES ) { - - out = match_cell(cell, dp->template, 0, dp->ltl, 0); - if ( out == NULL ) return 0; - - } else { - out = cell_new_from_cell(cell); - } - cr = crystal_new(); if ( cr == NULL ) { ERROR("Failed to allocate crystal.\n"); return 0; } - - crystal_set_cell(cr, out); - - if ( dp->indm & INDEXING_CHECK_PEAKS ) { - if ( !peak_sanity_check(image, &cr, 1) ) { - crystal_free(cr); /* Frees the cell as well */ - cell_free(out); - return 0; - } - } - + crystal_set_cell(cr, cell); image_add_crystal(image, cr); return 1; @@ -1206,22 +1181,9 @@ void *asdf_prepare(IndexingMethod *indm, UnitCell *cell, struct detector *det, float *ltl) { struct asdf_private *dp; - int need_cell = 0; - - if ( *indm & INDEXING_CHECK_CELL_COMBINATIONS ) need_cell = 1; - if ( *indm & INDEXING_CHECK_CELL_AXES ) need_cell = 1; - - if ( need_cell && !cell_has_parameters(cell) ) { - ERROR("Altering your asdf flags because cell parameters were" - " not provided.\n"); - *indm &= ~INDEXING_CHECK_CELL_COMBINATIONS; - *indm &= ~INDEXING_CHECK_CELL_AXES; - } /* Flags that asdf knows about */ - *indm &= INDEXING_METHOD_MASK | INDEXING_CHECK_CELL_COMBINATIONS - | INDEXING_CHECK_CELL_AXES | INDEXING_CHECK_PEAKS - | INDEXING_CONTROL_FLAGS; + *indm &= INDEXING_METHOD_MASK; dp = malloc(sizeof(struct asdf_private)); if ( dp == NULL ) return NULL; diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c index 6610abc3..41436bb4 100644 --- a/libcrystfel/src/cell-utils.c +++ b/libcrystfel/src/cell-utils.c @@ -3,12 +3,12 @@ * * Unit Cell utility functions * - * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * Copyright © 2012 Lorenzo Galli * * Authors: - * 2009-2012,2014-2015 Thomas White <taw@physics.org> + * 2009-2012,2014-2017 Thomas White <taw@physics.org> * 2012 Lorenzo Galli * * This file is part of CrystFEL. @@ -1686,3 +1686,120 @@ double cell_get_volume(UnitCell *cell) return 1/rec_volume; } + + +static double moduli_check(double ax, double ay, double az, + double bx, double by, double bz) +{ + double ma = modulus(ax, ay, az); + double mb = modulus(bx, by, bz); + return fabs(ma-mb)/ma; +} + + +static int cells_are_similar(UnitCell *cell1, UnitCell *cell2, + const double ltl, const double atl) +{ + double asx1, asy1, asz1, bsx1, bsy1, bsz1, csx1, csy1, csz1; + double asx2, asy2, asz2, bsx2, bsy2, bsz2, csx2, csy2, csz2; + UnitCell *pcell1, *pcell2; + + /* Compare primitive cells, not centered */ + pcell1 = uncenter_cell(cell1, NULL); + pcell2 = uncenter_cell(cell2, NULL); + + cell_get_reciprocal(pcell1, &asx1, &asy1, &asz1, + &bsx1, &bsy1, &bsz1, + &csx1, &csy1, &csz1); + + cell_get_reciprocal(pcell2, &asx2, &asy2, &asz2, + &bsx2, &bsy2, &bsz2, + &csx2, &csy2, &csz2); + + + cell_free(pcell1); + cell_free(pcell2); + + if ( angle_between(asx1, asy1, asz1, asx2, asy2, asz2) > atl ) return 0; + if ( angle_between(bsx1, bsy1, bsz1, bsx2, bsy2, bsz2) > atl ) return 0; + if ( angle_between(csx1, csy1, csz1, csx2, csy2, csz2) > atl ) return 0; + + if ( moduli_check(asx1, asy1, asz1, asx2, asy2, asz2) > ltl ) return 0; + if ( moduli_check(bsx1, bsy1, bsz1, bsx2, bsy2, bsz2) > ltl ) return 0; + if ( moduli_check(csx1, csy1, csz1, csx2, csy2, csz2) > ltl ) return 0; + + return 1; +} + + + +/** + * compare_cells: + * @a: A %UnitCell + * @b: Another %UnitCell + * @ltl: Maximum allowable fractional difference in reciprocal axis lengths + * @atl: Maximum allowable difference in reciprocal angles (in radians) + * @pmb: Place to store pointer to matrix + * + * Compare the two units cells. If they agree to within @ltl and @atl, using + * any change of axes, returns non-zero and stores the transformation to map @b + * onto @a. + * + * Returns: non-zero if the cells match. + * + */ +int compare_cells(UnitCell *a, UnitCell *b, double ltl, double atl, + IntegerMatrix **pmb) +{ + IntegerMatrix *m; + int i[9]; + + m = intmat_new(3, 3); + + for ( i[0]=-1; i[0]<=+1; i[0]++ ) { + for ( i[1]=-1; i[1]<=+1; i[1]++ ) { + for ( i[2]=-1; i[2]<=+1; i[2]++ ) { + for ( i[3]=-1; i[3]<=+1; i[3]++ ) { + for ( i[4]=-1; i[4]<=+1; i[4]++ ) { + for ( i[5]=-1; i[5]<=+1; i[5]++ ) { + for ( i[6]=-1; i[6]<=+1; i[6]++ ) { + for ( i[7]=-1; i[7]<=+1; i[7]++ ) { + for ( i[8]=-1; i[8]<=+1; i[8]++ ) { + + UnitCellTransformation *tfn; + UnitCell *nc; + int j, k; + int l = 0; + + for ( j=0; j<3; j++ ) + for ( k=0; k<3; k++ ) + intmat_set(m, j, k, i[l++]); + + if ( intmat_det(m) != +1 ) continue; + + tfn = tfn_from_intmat(m); + nc = cell_transform(b, tfn); + + if ( cells_are_similar(a, nc, ltl, atl) ) { + if ( pmb != NULL ) *pmb = m; + tfn_free(tfn); + cell_free(nc); + return 1; + } + + tfn_free(tfn); + cell_free(nc); + + } + } + } + } + } + } + } + } + } + + intmat_free(m); + return 0; +} diff --git a/libcrystfel/src/cell-utils.h b/libcrystfel/src/cell-utils.h index efb4b25b..5e2b2825 100644 --- a/libcrystfel/src/cell-utils.h +++ b/libcrystfel/src/cell-utils.h @@ -3,13 +3,13 @@ * * Unit Cell utility functions * - * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * Copyright © 2012 Lorenzo Galli * * Authors: - * 2009-2013,2014 Thomas White <taw@physics.org> - * 2012 Lorenzo Galli + * 2009-2013,2014,2017 Thomas White <taw@physics.org> + * 2012 Lorenzo Galli * * This file is part of CrystFEL. * @@ -80,6 +80,9 @@ extern int forbidden_reflection(UnitCell *cell, extern double cell_get_volume(UnitCell *cell); +extern int compare_cells(UnitCell *a, UnitCell *b, double ltl, double atl, + IntegerMatrix **pmb); + #ifdef __cplusplus } #endif diff --git a/libcrystfel/src/dirax.c b/libcrystfel/src/dirax.c index b263eea0..29b735ad 100644 --- a/libcrystfel/src/dirax.c +++ b/libcrystfel/src/dirax.c @@ -108,39 +108,14 @@ struct dirax_data { static int check_cell(struct dirax_private *dp, struct image *image, UnitCell *cell) { - UnitCell *out; Crystal *cr; - if ( dp->indm & INDEXING_CHECK_CELL_COMBINATIONS ) { - - out = match_cell(cell, dp->template, 0, dp->ltl, 1); - if ( out == NULL ) return 0; - - } else if ( dp->indm & INDEXING_CHECK_CELL_AXES ) { - - out = match_cell(cell, dp->template, 0, dp->ltl, 0); - if ( out == NULL ) return 0; - - } else { - out = cell_new_from_cell(cell); - } - cr = crystal_new(); if ( cr == NULL ) { ERROR("Failed to allocate crystal.\n"); return 0; } - - crystal_set_cell(cr, out); - - if ( dp->indm & INDEXING_CHECK_PEAKS ) { - if ( !peak_sanity_check(image, &cr, 1) ) { - crystal_free(cr); /* Frees the cell as well */ - cell_free(out); - return 0; - } - } - + crystal_set_cell(cr, cell); image_add_crystal(image, cr); return 1; @@ -655,9 +630,7 @@ void *dirax_prepare(IndexingMethod *indm, UnitCell *cell, } /* Flags that DirAx knows about */ - *indm &= INDEXING_METHOD_MASK | INDEXING_CHECK_CELL_COMBINATIONS - | INDEXING_CHECK_CELL_AXES | INDEXING_CHECK_PEAKS - | INDEXING_CONTROL_FLAGS; + *indm &= INDEXING_METHOD_MASK; dp = malloc(sizeof(struct dirax_private)); if ( dp == NULL ) return NULL; diff --git a/libcrystfel/src/felix.c b/libcrystfel/src/felix.c index 6daf4408..4c12778b 100644 --- a/libcrystfel/src/felix.c +++ b/libcrystfel/src/felix.c @@ -222,16 +222,7 @@ static int read_felix(struct felix_private *gp, struct image *image, fclose(fh); - if ( gp->indm & INDEXING_CHECK_PEAKS ) { - if ( !peak_sanity_check(image, image->crystals, - image->n_crystals) ) - { - free_all_crystals(image); - return 0; - } - } - - return n_crystals; + return n_crystals; } @@ -672,9 +663,8 @@ void *felix_prepare(IndexingMethod *indm, UnitCell *cell, if ( gp == NULL ) return NULL; /* Flags that Felix knows about */ - *indm &= INDEXING_METHOD_MASK | INDEXING_CHECK_PEAKS - | INDEXING_USE_LATTICE_TYPE | INDEXING_USE_CELL_PARAMETERS - | INDEXING_CONTROL_FLAGS; + *indm &= INDEXING_METHOD_MASK + | INDEXING_USE_LATTICE_TYPE | INDEXING_USE_CELL_PARAMETERS; gp->cell = cell; gp->indm = *indm; diff --git a/libcrystfel/src/hdf5-file.c b/libcrystfel/src/hdf5-file.c index 13339506..6fc10357 100644 --- a/libcrystfel/src/hdf5-file.c +++ b/libcrystfel/src/hdf5-file.c @@ -409,6 +409,8 @@ int get_peaks_cxi_2(struct image *image, struct hdfile *f, const char *p, } image->features = image_feature_list_new(); + image->num_peaks = 0; + image->num_saturated_peaks = 0; for ( pk=0; pk<num_peaks; pk++ ) { float fs, ss, val; @@ -427,6 +429,7 @@ int get_peaks_cxi_2(struct image *image, struct hdfile *f, const char *p, ss = ss - p->orig_min_ss; image_add_feature(image->features, fs, ss, p, image, val, NULL); + image->num_peaks++; } @@ -554,6 +557,8 @@ int get_peaks_2(struct image *image, struct hdfile *f, const char *p, } image->features = image_feature_list_new(); + image->num_peaks = 0; + image->num_saturated_peaks = 0; for ( i=0; i<size[0]; i++ ) { float fs, ss, val; @@ -573,6 +578,7 @@ int get_peaks_2(struct image *image, struct hdfile *f, const char *p, image_add_feature(image->features, fs, ss, p, image, val, NULL); + image->num_peaks++; } diff --git a/libcrystfel/src/image.c b/libcrystfel/src/image.c index aad5017c..8de16bd5 100644 --- a/libcrystfel/src/image.c +++ b/libcrystfel/src/image.c @@ -288,9 +288,10 @@ void image_add_crystal(struct image *image, Crystal *cryst) } -void remove_flagged_crystals(struct image *image) +int remove_flagged_crystals(struct image *image) { int i; + int n_bad = 0; for ( i=0; i<image->n_crystals; i++ ) { if ( crystal_get_user_flag(image->crystals[i]) ) { @@ -302,9 +303,11 @@ void remove_flagged_crystals(struct image *image) image->crystals[j] = image->crystals[j+1]; } image->n_crystals--; + n_bad++; } } + return n_bad; } diff --git a/libcrystfel/src/image.h b/libcrystfel/src/image.h index d7884e4c..c900bd29 100644 --- a/libcrystfel/src/image.h +++ b/libcrystfel/src/image.h @@ -257,7 +257,7 @@ extern struct imagefeature *image_get_feature(ImageFeatureList *flist, int idx); extern ImageFeatureList *sort_peaks(ImageFeatureList *flist); extern void image_add_crystal(struct image *image, Crystal *cryst); -extern void remove_flagged_crystals(struct image *image); +extern int remove_flagged_crystals(struct image *image); extern void free_all_crystals(struct image *image); /* Image files */ diff --git a/libcrystfel/src/index.c b/libcrystfel/src/index.c index 644ea31f..8c1f52e5 100644 --- a/libcrystfel/src/index.c +++ b/libcrystfel/src/index.c @@ -61,6 +61,7 @@ struct _indexingprivate { + IndexingFlags flags; UnitCell *target_cell; float tolerance[4]; @@ -72,6 +73,38 @@ struct _indexingprivate }; +static const char *onoff(int a) +{ + if ( a ) return "on"; + return "off"; +} + + +static void show_indexing_flags(IndexingFlags flags) +{ + char check[64]; + + assert( !((flags & INDEXING_CHECK_CELL_COMBINATIONS) + && (flags & INDEXING_CHECK_CELL_AXES)) ); + strcpy(check, onoff(flags & (INDEXING_CHECK_CELL_COMBINATIONS | INDEXING_CHECK_CELL_AXES))); + if ( flags & INDEXING_CHECK_CELL_AXES ) { + strcat(check, " (axis permutations only)"); + } + if ( flags & INDEXING_CHECK_CELL_COMBINATIONS ) { + strcat(check, " (axis combinations)"); + } + STATUS(" Check unit cell parameters: %s\n", check); + STATUS(" Check peak alignment: %s\n", + onoff(flags & INDEXING_CHECK_PEAKS)); + STATUS(" Refine indexing solutions: %s\n", + onoff(flags & INDEXING_REFINE)); + STATUS(" Multi-lattice indexing (\"delete and retry\"): %s\n", + onoff(flags & INDEXING_MULTI)); + STATUS(" Retry indexing: %s\n", + onoff(flags & INDEXING_RETRY)); +} + + static int debug_index(struct image *image) { Crystal *cr = crystal_new(); @@ -159,7 +192,7 @@ static void *prepare_method(IndexingMethod *m, UnitCell *cell, IndexingPrivate *setup_indexing(const char *method_list, UnitCell *cell, struct detector *det, float *ltl, - int no_refine, const char *options, + IndexingFlags flags, const char *options, struct taketwo_options *ttopts) { int i, n; @@ -178,10 +211,30 @@ IndexingPrivate *setup_indexing(const char *method_list, UnitCell *cell, for ( i=0; i<n; i++ ) { methods[i] = get_indm_from_string(method_strings[i]); if ( methods[i] == INDEXING_ERROR ) { + ERROR("----- Notice -----\n"); + ERROR("The way indexing options are used has changed in this CrystFEL version.\n"); + ERROR("To disable prediction refinement ('norefine'), use --no-refine.\n"); + ERROR("To check cell axes only ('axes'), use --no-cell-combinations.\n"); + ERROR("To disable all unit cell checks ('raw'), use --no-check-cell.\n"); + ERROR("To disable indexing retry ('noretry'), use --no-retry.\n"); + ERROR("Multi-lattice indexing ('multi') is now the default: " + "use --no-multi to disable it.\n"); + ERROR("------------------\n"); free(methods); return NULL; } - if ( no_refine ) methods[i] &= ~INDEXING_REFINE; + } + + if ( cell == NULL ) { + if ( (flags & INDEXING_CHECK_CELL_COMBINATIONS) + || (flags & INDEXING_CHECK_CELL_COMBINATIONS) ) + { + ERROR("WARNING: Forcing --no-cell-combinations " + "and --no-check-cell because you didn't " + "provide a unit cell.\n"); + flags &= ~INDEXING_CHECK_CELL_COMBINATIONS; + flags &= ~INDEXING_CHECK_CELL_AXES; + } } ipriv = malloc(sizeof(struct _indexingprivate)); @@ -212,6 +265,7 @@ IndexingPrivate *setup_indexing(const char *method_list, UnitCell *cell, ipriv->methods = methods; ipriv->n_methods = n; + ipriv->flags = flags; if ( cell != NULL ) { ipriv->target_cell = cell_new_from_cell(cell); @@ -222,6 +276,8 @@ IndexingPrivate *setup_indexing(const char *method_list, UnitCell *cell, ipriv->ttopts = ttopts; + show_indexing_flags(flags); + return ipriv; } @@ -308,6 +364,36 @@ void map_all_peaks(struct image *image) } +static int check_cell(IndexingFlags flags, Crystal *cr, UnitCell *target, + float *tolerance) +{ + if ( (flags & INDEXING_CHECK_CELL_COMBINATIONS) + || (flags & INDEXING_CHECK_CELL_AXES) ) + { + UnitCell *out; + int reduce; + + if ( flags & INDEXING_CHECK_CELL_COMBINATIONS ) + { + reduce = 1; + } else { + reduce = 0; + } + + out = match_cell(crystal_get_cell(cr), + target, 0, tolerance, reduce); + + if ( out == NULL ) { + return 1; + } + + cell_free(crystal_get_cell(cr)); + crystal_set_cell(cr, out); + } + return 0; +} + + /* Return non-zero for "success" */ static int try_indexer(struct image *image, IndexingMethod indm, IndexingPrivate *ipriv, void *mpriv) @@ -354,53 +440,86 @@ static int try_indexer(struct image *image, IndexingMethod indm, } + /* For all the crystals found this time ... */ for ( i=0; i<r; i++ ) { - Crystal *cr = image->crystals[image->n_crystals-i-1]; + int j; + int this_crystal = image->n_crystals - i - 1; + + /* ... starting at the end of the (complete) list ... */ + Crystal *cr = image->crystals[this_crystal]; + crystal_set_image(cr, image); crystal_set_user_flag(cr, 0); crystal_set_profile_radius(cr, 0.02e9); crystal_set_mosaicity(cr, 0.0); - /* Prediction refinement if requested */ - if ( indm & INDEXING_REFINE ) { + assert( !((ipriv->flags & INDEXING_CHECK_CELL_COMBINATIONS) + && (ipriv->flags & INDEXING_CHECK_CELL_AXES)) ); - UnitCell *out; + /* Pre-refinement unit cell check if requested */ + if ( check_cell(ipriv->flags, cr, ipriv->target_cell, + ipriv->tolerance) ) + { + crystal_set_user_flag(cr, 1); + n_bad++; + continue; + } + /* Prediction refinement if requested */ + if ( ipriv->flags & INDEXING_REFINE ) + { if ( refine_prediction(image, cr) ) { crystal_set_user_flag(cr, 1); n_bad++; continue; } + } - if ( (indm & INDEXING_CHECK_CELL_COMBINATIONS) - || (indm & INDEXING_CHECK_CELL_AXES) ) - { + /* After refinement unit cell check if requested */ + if ( check_cell(ipriv->flags, cr, ipriv->target_cell, + ipriv->tolerance) ) + { + crystal_set_user_flag(cr, 1); + n_bad++; + continue; + } - /* Check that the cell parameters are still - * within the tolerance */ - out = match_cell(crystal_get_cell(cr), - ipriv->target_cell, 0, - ipriv->tolerance, 0); + /* Peak alignment check if requested */ + if ( ipriv->flags & INDEXING_CHECK_PEAKS ) + { + if ( !peak_sanity_check(image, &cr, 1) ) { + crystal_set_user_flag(cr, 1); + n_bad++; + continue; + } + } - if ( out == NULL ) { - crystal_set_user_flag(cr, 1); - n_bad++; - } + /* Don't do similarity check if this crystal is bad */ + if ( crystal_get_user_flag(cr) ) continue; - cell_free(out); + /* Check if cell is too similar to existing ones */ + for ( j=0; j<this_crystal; j++ ) { - } + Crystal *that_cr = image->crystals[j]; + /* Don't do similarity check against bad crystals */ + if ( crystal_get_user_flag(that_cr) ) continue; + + if ( compare_cells(crystal_get_cell(cr), + crystal_get_cell(that_cr), + 0.1, deg2rad(5.0), NULL) ) + { + crystal_set_user_flag(cr, 1); + } } } - remove_flagged_crystals(image); - - if ( n_bad == r ) return 0; + n_bad = remove_flagged_crystals(image); + assert(r >= n_bad); - return r; + return r - n_bad; } @@ -490,14 +609,15 @@ static int delete_explained_peaks(struct image *image, Crystal *cr) * * Returns false for "try again", true for "no, stop now" */ -static int finished_retry(IndexingMethod indm, int r, struct image *image) +static int finished_retry(IndexingMethod indm, IndexingFlags flags, + int r, struct image *image) { if ( r == 0 ) { /* Indexing failed on the previous attempt. Maybe try again * after poking the peak list a bit */ - if ( indm & INDEXING_RETRY ) { + if ( flags & INDEXING_RETRY ) { /* Retry with fewer peaks */ return delete_weakest_peaks(image->features); } else { @@ -510,13 +630,10 @@ static int finished_retry(IndexingMethod indm, int r, struct image *image) /* Indexing succeeded on previous attempt. Maybe try again * after deleting the explained peaks */ - if ( indm & INDEXING_MULTI ) { - /* Remove "used" spots and try for - * another lattice */ + if ( flags & INDEXING_MULTI ) { + /* Remove "used" spots and try for another lattice */ Crystal *cr; cr = image->crystals[image->n_crystals-1]; - STATUS("WARNING: Multi-lattice indexing does not work" - " well in this version.\n"); return delete_explained_peaks(image, cr); } else { return 1; @@ -525,6 +642,7 @@ static int finished_retry(IndexingMethod indm, int r, struct image *image) } } + void index_pattern(struct image *image, IndexingPrivate *ipriv) { index_pattern_2(image, ipriv, NULL); @@ -559,7 +677,8 @@ void index_pattern_2(struct image *image, IndexingPrivate *ipriv, int *ping) ipriv, ipriv->engine_private[n]); success += r; ntry++; - done = finished_retry(ipriv->methods[n], r, image); + done = finished_retry(ipriv->methods[n], ipriv->flags, + r, image); if ( ntry > 5 ) done = 1; if ( ping != NULL ) (*ping)++; @@ -583,39 +702,6 @@ void index_pattern_2(struct image *image, IndexingPrivate *ipriv, int *ping) } -/* Set the indexer flags for "raw mode" ("--cell-reduction=none") */ -static IndexingMethod set_raw(IndexingMethod a) -{ - /* Disable all unit cell checks */ - a &= ~(INDEXING_CHECK_CELL_COMBINATIONS | INDEXING_CHECK_CELL_AXES); - return a; -} - - -/* Set the indexer flags for "bad mode" ("--insane) */ -static IndexingMethod set_bad(IndexingMethod a) -{ - /* Disable the peak check */ - return a & ~INDEXING_CHECK_PEAKS; -} - - -/* Set the indexer flags for "axes mode" ("--cell-reduction=compare") */ -static IndexingMethod set_axes(IndexingMethod a) -{ - return (a & ~INDEXING_CHECK_CELL_COMBINATIONS) - | INDEXING_CHECK_CELL_AXES; -} - - -/* Set the indexer flags for "combination mode" ("--cell-reduction=reduce") */ -static IndexingMethod set_comb(IndexingMethod a) -{ - return (a & ~INDEXING_CHECK_CELL_AXES) - | INDEXING_CHECK_CELL_COMBINATIONS; -} - - /* Set the indexer flags for "use no lattice type information" */ static IndexingMethod set_nolattice(IndexingMethod a) { @@ -703,18 +789,6 @@ char *indexer_str(IndexingMethod indm) if ( (indm & INDEXING_METHOD_MASK) == INDEXING_SIMULATION ) return str; - if ( indm & INDEXING_CHECK_CELL_COMBINATIONS ) { - strcat(str, "-comb"); - } else if ( indm & INDEXING_CHECK_CELL_AXES ) { - strcat(str, "-axes"); - } else { - strcat(str, "-raw"); - } - - if ( !(indm & INDEXING_CHECK_PEAKS) ) { - strcat(str, "-bad"); - } - if ( indm & INDEXING_USE_LATTICE_TYPE ) { strcat(str, "-latt"); } else { @@ -727,24 +801,6 @@ char *indexer_str(IndexingMethod indm) strcat(str, "-nocell"); } - if ( indm & INDEXING_RETRY ) { - strcat(str, "-retry"); - } else { - strcat(str, "-noretry"); - } - - if ( indm & INDEXING_MULTI ) { - strcat(str, "-multi"); - } else { - strcat(str, "-nomulti"); - } - - if ( indm & INDEXING_REFINE ) { - strcat(str, "-refine"); - } else { - strcat(str, "-norefine"); - } - return str; } @@ -813,18 +869,6 @@ IndexingMethod get_indm_from_string(const char *str) method = INDEXING_DEBUG; return method; - } else if ( strcmp(bits[i], "raw") == 0) { - method = set_raw(method); - - } else if ( strcmp(bits[i], "bad") == 0) { - method = set_bad(method); - - } else if ( strcmp(bits[i], "comb") == 0) { - method = set_comb(method); /* Default */ - - } else if ( strcmp(bits[i], "axes") == 0) { - method = set_axes(method); - } else if ( strcmp(bits[i], "latt") == 0) { method = set_lattice(method); @@ -837,24 +881,6 @@ IndexingMethod get_indm_from_string(const char *str) } else if ( strcmp(bits[i], "nocell") == 0) { method = set_nocellparams(method); - } else if ( strcmp(bits[i], "retry") == 0) { - method |= INDEXING_RETRY; - - } else if ( strcmp(bits[i], "noretry") == 0) { - method &= ~INDEXING_RETRY; - - } else if ( strcmp(bits[i], "multi") == 0) { - method |= INDEXING_MULTI; - - } else if ( strcmp(bits[i], "nomulti") == 0) { - method &= ~INDEXING_MULTI; - - } else if ( strcmp(bits[i], "refine") == 0) { - method |= INDEXING_REFINE; - - } else if ( strcmp(bits[i], "norefine") == 0) { - method &= ~INDEXING_REFINE; - } else { ERROR("Bad list of indexing methods: '%s'\n", str); return INDEXING_ERROR; diff --git a/libcrystfel/src/index.h b/libcrystfel/src/index.h index 8b3aadeb..2aad377b 100644 --- a/libcrystfel/src/index.h +++ b/libcrystfel/src/index.h @@ -39,36 +39,22 @@ #endif -#define INDEXING_DEFAULTS_DIRAX (INDEXING_DIRAX | INDEXING_CHECK_PEAKS \ - | INDEXING_CHECK_CELL_COMBINATIONS \ - | INDEXING_RETRY | INDEXING_REFINE) - -#define INDEXING_DEFAULTS_ASDF (INDEXING_ASDF | INDEXING_CHECK_PEAKS \ - | INDEXING_CHECK_CELL_COMBINATIONS \ - | INDEXING_RETRY | INDEXING_REFINE) - -#define INDEXING_DEFAULTS_MOSFLM (INDEXING_MOSFLM | INDEXING_CHECK_PEAKS \ - | INDEXING_CHECK_CELL_COMBINATIONS \ - | INDEXING_USE_LATTICE_TYPE \ - | INDEXING_USE_CELL_PARAMETERS \ - | INDEXING_RETRY | INDEXING_REFINE) - -#define INDEXING_DEFAULTS_FELIX (INDEXING_FELIX \ - | INDEXING_USE_LATTICE_TYPE \ - | INDEXING_USE_CELL_PARAMETERS \ - | INDEXING_RETRY | INDEXING_REFINE) - -#define INDEXING_DEFAULTS_TAKETWO (INDEXING_TAKETWO | INDEXING_CHECK_PEAKS \ - | INDEXING_USE_CELL_PARAMETERS \ - | INDEXING_USE_LATTICE_TYPE \ - | INDEXING_RETRY | INDEXING_REFINE) - -/* Axis check is needed for XDS, because it likes to permute the axes */ -#define INDEXING_DEFAULTS_XDS (INDEXING_XDS | INDEXING_USE_LATTICE_TYPE \ - | INDEXING_USE_CELL_PARAMETERS \ - | INDEXING_CHECK_CELL_AXES \ - | INDEXING_CHECK_PEAKS \ - | INDEXING_RETRY | INDEXING_REFINE) +#define INDEXING_DEFAULTS_DIRAX (INDEXING_DIRAX) + +#define INDEXING_DEFAULTS_ASDF (INDEXING_ASDF) + +#define INDEXING_DEFAULTS_MOSFLM (INDEXING_MOSFLM | INDEXING_USE_LATTICE_TYPE \ + | INDEXING_USE_CELL_PARAMETERS) + +#define INDEXING_DEFAULTS_FELIX (INDEXING_FELIX | INDEXING_USE_LATTICE_TYPE \ + | INDEXING_USE_CELL_PARAMETERS) + +#define INDEXING_DEFAULTS_TAKETWO (INDEXING_TAKETWO \ + | INDEXING_USE_CELL_PARAMETERS \ + | INDEXING_USE_LATTICE_TYPE) + +#define INDEXING_DEFAULTS_XDS (INDEXING_XDS | INDEXING_USE_LATTICE_TYPE \ + | INDEXING_USE_CELL_PARAMETERS) /** * IndexingMethod: @@ -82,21 +68,10 @@ * @INDEXING_ASDF: Use in-built "asdf" indexer * @INDEXING_TAKETWO: Use in-built "taketwo" 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 - * cell, and permute them if necessary. - * @INDEXING_CHECK_PEAKS: Check that the peaks can be explained by the indexing - * result. * @INDEXING_USE_LATTICE_TYPE: Use lattice type and centering information to * guide the indexing process. * @INDEXING_USE_CELL_PARAMETERS: Use the unit cell parameters to guide the * indexing process. - * @INDEXING_RETRY: If the indexer doesn't succeed, delete the weakest peaks - * and try again. - * @INDEXING_MULTI: If the indexer succeeds, delete the peaks explained by the - * lattice and try again in the hope of finding another crystal. - * @INDEXING_REFINE: Perform "prediction refinement" after indexing. * * An enumeration of all the available indexing methods. The dummy value * @INDEXING_SIMULATION is used by partial_sim to indicate that no indexing was @@ -120,14 +95,8 @@ typedef enum { /* Bits at the top of the IndexingMethod are flags which modify the * behaviour of the indexer. */ - INDEXING_CHECK_CELL_COMBINATIONS = 256, - INDEXING_CHECK_CELL_AXES = 512, - INDEXING_CHECK_PEAKS = 1024, INDEXING_USE_LATTICE_TYPE = 2048, INDEXING_USE_CELL_PARAMETERS = 4096, - INDEXING_RETRY = 8192, - INDEXING_MULTI = 16384, - INDEXING_REFINE = 32768, } IndexingMethod; @@ -135,8 +104,16 @@ typedef enum { * core of the indexing method */ #define INDEXING_METHOD_MASK (0xff) -/* Indexing flags which the indexing method does not need to know about */ -#define INDEXING_CONTROL_FLAGS (INDEXING_RETRY | INDEXING_MULTI | INDEXING_REFINE) +typedef enum { + + INDEXING_RETRY = 1, + INDEXING_MULTI = 2, + INDEXING_REFINE = 4, + INDEXING_CHECK_CELL_COMBINATIONS = 8, + INDEXING_CHECK_CELL_AXES = 16, + INDEXING_CHECK_PEAKS = 32, + +} IndexingFlags; #ifdef __cplusplus extern "C" { @@ -162,7 +139,7 @@ extern IndexingMethod get_indm_from_string(const char *method); extern IndexingPrivate *setup_indexing(const char *methods, UnitCell *cell, struct detector *det, float *ltl, - int no_refine, const char *options, + IndexingFlags flags, const char *options, struct taketwo_options *ttopts); extern void index_pattern(struct image *image, IndexingPrivate *ipriv); diff --git a/libcrystfel/src/mosflm.c b/libcrystfel/src/mosflm.c index edf8524e..110243cb 100644 --- a/libcrystfel/src/mosflm.c +++ b/libcrystfel/src/mosflm.c @@ -128,7 +128,6 @@ struct mosflm_data { static int check_cell(struct mosflm_private *mp, struct image *image, UnitCell *cell) { - UnitCell *out; Crystal *cr; /* If we sent lattice information, make sure that we got back what we @@ -163,35 +162,13 @@ static int check_cell(struct mosflm_private *mp, struct image *image, } - if ( mp->indm & INDEXING_CHECK_CELL_COMBINATIONS ) { - - out = match_cell(cell, mp->template, 0, mp->ltl, 1); - if ( out == NULL ) return 0; - - } else if ( mp->indm & INDEXING_CHECK_CELL_AXES ) { - - out = match_cell(cell, mp->template, 0, mp->ltl, 0); - if ( out == NULL ) return 0; - - } else { - out = cell_new_from_cell(cell); - } - cr = crystal_new(); if ( cr == NULL ) { ERROR("Failed to allocate crystal.\n"); return 0; } - crystal_set_cell(cr, out); - - if ( mp->indm & INDEXING_CHECK_PEAKS ) { - if ( !peak_sanity_check(image, &cr, 1) ) { - cell_free(out); - crystal_free(cr); - return 0; - } - } + crystal_set_cell(cr, cell); image_add_crystal(image, cr); @@ -340,8 +317,6 @@ static int read_newmat(struct mosflm_data *mosflm, const char *filename, mosflm->done = 1; } - cell_free(cell); - return 0; } @@ -850,14 +825,10 @@ void *mosflm_prepare(IndexingMethod *indm, UnitCell *cell, int need_cell = 0; /* Check if cell parameters are needed/provided */ - if ( *indm & INDEXING_CHECK_CELL_COMBINATIONS ) need_cell = 1; - if ( *indm & INDEXING_CHECK_CELL_AXES ) need_cell = 1; if ( *indm & INDEXING_USE_CELL_PARAMETERS ) need_cell = 1; if ( need_cell && !cell_has_parameters(cell) ) { ERROR("Altering your MOSFLM flags because cell parameters were" " not provided.\n"); - *indm &= ~INDEXING_CHECK_CELL_COMBINATIONS; - *indm &= ~INDEXING_CHECK_CELL_AXES; *indm &= ~INDEXING_USE_CELL_PARAMETERS; } @@ -869,36 +840,8 @@ void *mosflm_prepare(IndexingMethod *indm, UnitCell *cell, } /* Flags that MOSFLM knows about */ - *indm &= INDEXING_METHOD_MASK | INDEXING_CHECK_CELL_COMBINATIONS - | INDEXING_CHECK_CELL_AXES | INDEXING_CHECK_PEAKS - | INDEXING_USE_LATTICE_TYPE| INDEXING_USE_CELL_PARAMETERS - | INDEXING_CONTROL_FLAGS; - - if ( (*indm & INDEXING_USE_LATTICE_TYPE) - && !((*indm & INDEXING_CHECK_CELL_COMBINATIONS) - || (*indm & INDEXING_CHECK_CELL_AXES)) - && (cell_has_parameters(cell) - || (cell_get_unique_axis(cell) == 'a') - || (cell_get_unique_axis(cell) == 'b') - || (cell_get_unique_axis(cell) == 'c')) ) - { - ERROR("WARNING: The unit cell from %s might have had " - "its axes permuted from the unit cell you gave.\n" - "If this is a problem, consider using " - "mosflm-axes-latt or mosflm-comb-latt instead of " - "mosflm-raw-latt.\n", indexer_str(*indm)); - /* It'll be OK if INDEXING_CHECK_CELL_PARAMETERS and - * MOSFLM 7.2.0 or later, but that's getting a bit too - * complicated to explain to the user. */ - } - - if ( (*indm & INDEXING_USE_CELL_PARAMETERS) - && !((*indm & INDEXING_CHECK_CELL_COMBINATIONS) - || (*indm & INDEXING_CHECK_CELL_AXES)) ) - { - ERROR("WARNING: Using 'mosflm-raw' but providing cell " - "parameters as prior information to mosflm.\n"); - } + *indm &= INDEXING_METHOD_MASK + | INDEXING_USE_LATTICE_TYPE | INDEXING_USE_CELL_PARAMETERS; mp = malloc(sizeof(struct mosflm_private)); if ( mp == NULL ) return NULL; diff --git a/libcrystfel/src/predict-refine.c b/libcrystfel/src/predict-refine.c index 61ce3629..e702cf33 100644 --- a/libcrystfel/src/predict-refine.c +++ b/libcrystfel/src/predict-refine.c @@ -138,10 +138,12 @@ static void UNUSED write_pairs(const char *filename, struct reflpeak *rps, double write_fs, write_ss; double fs, ss; struct panel *p; + signed int h, k, l; fs = rps[i].peak->fs; ss = rps[i].peak->ss; p = rps[i].panel; + get_indices(rps[i].refl, &h, &k, &l); write_fs = fs + p->orig_min_fs; write_ss = ss + p->orig_min_ss; @@ -152,6 +154,9 @@ static void UNUSED write_pairs(const char *filename, struct reflpeak *rps, x_dev(&rps[i], det), y_dev(&rps[i], det)); + //fprintf(fh, "%4i %4i %4i 0.0 - 0.0 1 %7.2f %7.2f %s\n", + // h, k, l, write_fs, write_ss, p->name); + } fclose(fh); @@ -483,6 +488,18 @@ static int iterate(struct reflpeak *rps, int n, UnitCell *cell, } + int k; + for ( k=0; k<num_params; k++ ) { + double M_curr; + M_curr = gsl_matrix_get(M, k, k); + if ( (rv[k] == GPARAM_DETX) || (rv[k] == GPARAM_DETY) ) { + M_curr += 10.0; + } else { + M_curr += 1e-18; + } + gsl_matrix_set(M, k, k, M_curr); + } + //show_matrix_eqn(M, v); shifts = solve_svd(v, M, NULL, 0); if ( shifts == NULL ) { diff --git a/libcrystfel/src/taketwo.c b/libcrystfel/src/taketwo.c index 87da44a1..f7a2fda9 100644 --- a/libcrystfel/src/taketwo.c +++ b/libcrystfel/src/taketwo.c @@ -1460,7 +1460,7 @@ static void cleanup_taketwo_cell(struct TakeTwoCell *ttCell) * successful. **/ static UnitCell *run_taketwo(UnitCell *cell, const struct taketwo_options *opts, - struct rvec *rlps, int rlp_count) + struct rvec *rlps, int rlp_count) { int cell_vec_count = 0; int asym_vec_count = 0; @@ -1581,18 +1581,6 @@ int taketwo_index(struct image *image, const struct taketwo_options *opts, crystal_set_cell(cr, cell); - if ( tp->indm & INDEXING_CHECK_PEAKS ) { - if ( !peak_sanity_check(image, &cr, 1) ) { - cell_free(cell); - crystal_free(cr); - // STATUS("Rubbish!!\n"); - - return 0; - } else { - // STATUS("That's good!\n"); - } - } - image_add_crystal(image, cr); return 1; @@ -1605,9 +1593,8 @@ void *taketwo_prepare(IndexingMethod *indm, UnitCell *cell, struct taketwo_private *tp; /* Flags that TakeTwo knows about */ - *indm &= INDEXING_METHOD_MASK | INDEXING_CHECK_PEAKS - | INDEXING_USE_LATTICE_TYPE | INDEXING_USE_CELL_PARAMETERS - | INDEXING_CONTROL_FLAGS; + *indm &= INDEXING_METHOD_MASK | INDEXING_USE_LATTICE_TYPE + | INDEXING_USE_CELL_PARAMETERS; if ( !( (*indm & INDEXING_USE_LATTICE_TYPE) && (*indm & INDEXING_USE_CELL_PARAMETERS)) ) diff --git a/libcrystfel/src/xds.c b/libcrystfel/src/xds.c index 093b9160..39578c60 100644 --- a/libcrystfel/src/xds.c +++ b/libcrystfel/src/xds.c @@ -187,39 +187,14 @@ static int xds_readable(struct image *image, struct xds_data *xds) static int check_cell(struct xds_private *xp, struct image *image, UnitCell *cell) { - UnitCell *out; Crystal *cr; - if ( xp->indm & INDEXING_CHECK_CELL_COMBINATIONS ) { - - out = match_cell(cell, xp->cell, 0, xp->ltl, 1); - if ( out == NULL ) return 0; - - } else if ( xp->indm & INDEXING_CHECK_CELL_AXES ) { - - out = match_cell(cell, xp->cell, 0, xp->ltl, 0); - if ( out == NULL ) return 0; - - } else { - out = cell_new_from_cell(cell); - } - cr = crystal_new(); if ( cr == NULL ) { ERROR("Failed to allocate crystal.\n"); return 0; } - - crystal_set_cell(cr, out); - - if ( xp->indm & INDEXING_CHECK_PEAKS ) { - if ( !peak_sanity_check(image, &cr, 1) ) { - cell_free(out); - crystal_free(cr); - return 0; - } - } - + crystal_set_cell(cr, cell); image_add_crystal(image, cr); return 1; @@ -626,14 +601,10 @@ void *xds_prepare(IndexingMethod *indm, UnitCell *cell, int need_cell = 0; /* Check if cell parameters are needed/provided */ - if ( *indm & INDEXING_CHECK_CELL_COMBINATIONS ) need_cell = 1; - if ( *indm & INDEXING_CHECK_CELL_AXES ) need_cell = 1; if ( *indm & INDEXING_USE_CELL_PARAMETERS ) need_cell = 1; if ( need_cell && !cell_has_parameters(cell) ) { ERROR("Altering your XDS flags because cell parameters were" " not provided.\n"); - *indm &= ~INDEXING_CHECK_CELL_COMBINATIONS; - *indm &= ~INDEXING_CHECK_CELL_AXES; *indm &= ~INDEXING_USE_CELL_PARAMETERS; } @@ -652,36 +623,35 @@ void *xds_prepare(IndexingMethod *indm, UnitCell *cell, } if ( (*indm & INDEXING_USE_LATTICE_TYPE) - && !(*indm & INDEXING_USE_CELL_PARAMETERS) ) { + && !(*indm & INDEXING_USE_CELL_PARAMETERS) ) + { ERROR("Invalid XDS options (-latt-nocell): " "try xds-nolatt-nocell.\n"); return NULL; } if ( (*indm & INDEXING_USE_CELL_PARAMETERS) - && !(*indm & INDEXING_USE_LATTICE_TYPE) ) { + && !(*indm & INDEXING_USE_LATTICE_TYPE) ) + { ERROR("Invalid XDS options (-cell-nolatt): " "try xds-nolatt-nocell.\n"); return NULL; } if ( ((*indm & INDEXING_USE_CELL_PARAMETERS) - || (*indm & INDEXING_USE_LATTICE_TYPE)) - && !(*indm & INDEXING_CHECK_CELL_AXES) - && !(*indm & INDEXING_CHECK_CELL_COMBINATIONS) ) { - ERROR("The cell from xds-raw-cell or xds-raw-latt may have had" + || (*indm & INDEXING_USE_LATTICE_TYPE))) + { + ERROR("The cell from xds-cell or xds-latt may have had" " its axes permuted from the cell you provided. If this" - " is a problem, consider using xds-axes-cell.\n"); + " is a problem, consider using xds-cell.\n"); } xp = calloc(1, sizeof(*xp)); if ( xp == NULL ) return NULL; /* Flags that XDS knows about */ - *indm &= INDEXING_METHOD_MASK | INDEXING_CHECK_CELL_COMBINATIONS - | INDEXING_CHECK_CELL_AXES | INDEXING_USE_LATTICE_TYPE - | INDEXING_CHECK_PEAKS | INDEXING_USE_CELL_PARAMETERS - | INDEXING_CONTROL_FLAGS; + *indm &= INDEXING_METHOD_MASK | INDEXING_USE_LATTICE_TYPE + | INDEXING_USE_CELL_PARAMETERS; xp->ltl = ltl; xp->cell = cell; diff --git a/src/indexamajig.c b/src/indexamajig.c index a79234d3..44b91338 100644 --- a/src/indexamajig.c +++ b/src/indexamajig.c @@ -227,8 +227,13 @@ int main(int argc, char *argv[]) struct beam_params beam; int have_push_res = 0; int len; - int no_refine = 0; char *command_line_peak_path = NULL; + int if_refine = 1; + int if_comb = 1; + int if_axes = 1; + int if_peaks = 0; + int if_multi = 1; + int if_retry = 1; /* Defaults */ iargs.cell = NULL; @@ -312,9 +317,14 @@ int main(int argc, char *argv[]) {"no-use-saturated", 0, &iargs.use_saturated, 0}, {"no-revalidate", 0, &iargs.no_revalidate, 1}, {"check-hdf5-snr", 0, &iargs.check_hdf5_snr, 1}, - {"no-refine", 0, &no_refine, 1}, {"profile", 0, &iargs.profile, 1}, {"no-half-pixel-shift",0, &iargs.half_pixel_shift, 0}, + {"no-refine", 0, &if_refine, 0}, + {"no-cell-combinations",0,&if_comb, 0}, + {"no-check-cell", 0, &if_axes, 0}, + {"check-peaks", 0, &if_peaks, 1}, + {"no-retry", 0, &if_retry, 0}, + {"no-multi", 0, &if_multi, 0}, /* Long-only options which don't actually do anything */ {"no-sat-corr", 0, &iargs.satcorr, 0}, @@ -856,8 +866,30 @@ int main(int argc, char *argv[]) } else { + IndexingFlags flags = 0; + + if ( if_axes ) { + flags |= INDEXING_CHECK_CELL_AXES; + } + if ( if_comb ) { + flags |= INDEXING_CHECK_CELL_COMBINATIONS; + flags &= ~INDEXING_CHECK_CELL_AXES; /* Not needed */ + } + if ( if_refine ) { + flags |= INDEXING_REFINE; + } + if ( if_peaks ) { + flags |= INDEXING_CHECK_PEAKS; + } + if ( if_multi ) { + flags |= INDEXING_MULTI; + } + if ( if_retry ) { + flags |= INDEXING_RETRY; + } + iargs.ipriv = setup_indexing(indm_str, iargs.cell, iargs.det, - iargs.tols, no_refine, + iargs.tols, flags, iargs.felix_options, &iargs.taketwo_opts); if ( iargs.ipriv == NULL ) { diff --git a/src/whirligig.c b/src/whirligig.c index c5e8084b..e8829880 100644 --- a/src/whirligig.c +++ b/src/whirligig.c @@ -279,107 +279,6 @@ static void find_and_process_series(struct window *win, int is_last_frame, } -static double moduli_check(double ax, double ay, double az, - double bx, double by, double bz) -{ - double ma = modulus(ax, ay, az); - double mb = modulus(bx, by, bz); - return fabs(ma-mb)/ma; -} - - -static int cells_are_similar(UnitCell *cell1, UnitCell *cell2) -{ - double asx1, asy1, asz1, bsx1, bsy1, bsz1, csx1, csy1, csz1; - double asx2, asy2, asz2, bsx2, bsy2, bsz2, csx2, csy2, csz2; - UnitCell *pcell1, *pcell2; - const double atl = deg2rad(5.0); - const double ltl = 0.1; - - /* Compare primitive cells, not centered */ - pcell1 = uncenter_cell(cell1, NULL); - pcell2 = uncenter_cell(cell2, NULL); - - cell_get_reciprocal(pcell1, &asx1, &asy1, &asz1, - &bsx1, &bsy1, &bsz1, - &csx1, &csy1, &csz1); - - cell_get_reciprocal(pcell2, &asx2, &asy2, &asz2, - &bsx2, &bsy2, &bsz2, - &csx2, &csy2, &csz2); - - - cell_free(pcell1); - cell_free(pcell2); - - if ( angle_between(asx1, asy1, asz1, asx2, asy2, asz2) > atl ) return 0; - if ( angle_between(bsx1, bsy1, bsz1, bsx2, bsy2, bsz2) > atl ) return 0; - if ( angle_between(csx1, csy1, csz1, csx2, csy2, csz2) > atl ) return 0; - - if ( moduli_check(asx1, asy1, asz1, asx2, asy2, asz2) > ltl ) return 0; - if ( moduli_check(bsx1, bsy1, bsz1, bsx2, bsy2, bsz2) > ltl ) return 0; - if ( moduli_check(csx1, csy1, csz1, csx2, csy2, csz2) > ltl ) return 0; - - return 1; -} - - -static int gatinator(UnitCell *a, UnitCell *b, IntegerMatrix **pmb) -{ - IntegerMatrix *m; - int i[9]; - - m = intmat_new(3, 3); - - for ( i[0]=-1; i[0]<=+1; i[0]++ ) { - for ( i[1]=-1; i[1]<=+1; i[1]++ ) { - for ( i[2]=-1; i[2]<=+1; i[2]++ ) { - for ( i[3]=-1; i[3]<=+1; i[3]++ ) { - for ( i[4]=-1; i[4]<=+1; i[4]++ ) { - for ( i[5]=-1; i[5]<=+1; i[5]++ ) { - for ( i[6]=-1; i[6]<=+1; i[6]++ ) { - for ( i[7]=-1; i[7]<=+1; i[7]++ ) { - for ( i[8]=-1; i[8]<=+1; i[8]++ ) { - - UnitCellTransformation *tfn; - UnitCell *nc; - int j, k; - int l = 0; - - for ( j=0; j<3; j++ ) - for ( k=0; k<3; k++ ) - intmat_set(m, j, k, i[l++]); - - if ( intmat_det(m) != +1 ) continue; - - tfn = tfn_from_intmat(m); - nc = cell_transform(b, tfn); - - if ( cells_are_similar(a, nc) ) { - *pmb = m; - tfn_free(tfn); - cell_free(nc); - return 1; - } - - tfn_free(tfn); - cell_free(nc); - - } - } - } - } - } - } - } - } - } - - intmat_free(m); - return 0; -} - - static int crystal_used(struct window *win, int pos, int cn) { int i; @@ -411,8 +310,9 @@ static IntegerMatrix *try_all(struct window *win, int n1, int n2, for ( i=0; i<i1->n_crystals; i++ ) { for ( j=0; j<i2->n_crystals; j++ ) { - if ( gatinator(crystal_get_cell(i1->crystals[i]), - crystal_get_cell(i2->crystals[j]), &m) ) + if ( compare_cells(crystal_get_cell(i1->crystals[i]), + crystal_get_cell(i2->crystals[j]), + 0.1, deg2rad(5.0), &m) ) { if ( !crystal_used(win, n1, i) && !crystal_used(win, n2, j) ) @@ -478,8 +378,9 @@ static int try_join(struct window *win, int sn) for ( j=0; j<win->img[win->join_ptr].n_crystals; j++ ) { Crystal *cr2; cr2 = win->img[win->join_ptr].crystals[j]; - if ( gatinator(ref, crystal_get_cell(cr2), - &win->mat[sn][win->join_ptr]) ) { + if ( compare_cells(ref, crystal_get_cell(cr2), + 0.1, deg2rad(5.0), + &win->mat[sn][win->join_ptr]) ) { win->ser[sn][win->join_ptr] = j; cell_free(ref); return 1; |