diff options
-rw-r--r-- | libcrystfel/src/cell-utils.c | 121 | ||||
-rw-r--r-- | libcrystfel/src/cell-utils.h | 9 | ||||
-rw-r--r-- | src/whirligig.c | 111 |
3 files changed, 131 insertions, 110 deletions
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/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; |