From 95cce19a7aeb6be7568ab73a5263913fe476591e Mon Sep 17 00:00:00 2001 From: Thomas White Date: Mon, 5 Aug 2019 14:25:04 +0200 Subject: Remove match_cell and match_cell_ab --- libcrystfel/src/cell-utils.c | 398 ------------------------------------------- libcrystfel/src/cell-utils.h | 5 - tests/cellcompare_check.c | 10 +- 3 files changed, 5 insertions(+), 408 deletions(-) diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c index e59f1fa1..2af3f293 100644 --- a/libcrystfel/src/cell-utils.c +++ b/libcrystfel/src/cell-utils.c @@ -574,404 +574,6 @@ UnitCell *uncenter_cell(UnitCell *in, IntegerMatrix **pC, RationalMatrix **pCi) } -#define MAX_CAND (1024) - -static int right_handed_vec(struct rvec a, struct rvec b, struct rvec c) -{ - struct rvec aCb; - double aCb_dot_c; - - /* "a" cross "b" */ - aCb.u = a.v*b.w - a.w*b.v; - aCb.v = - (a.u*b.w - a.w*b.u); - aCb.w = a.u*b.v - a.v*b.u; - - /* "a cross b" dot "c" */ - aCb_dot_c = aCb.u*c.u + aCb.v*c.v + aCb.w*c.w; - - if ( aCb_dot_c > 0.0 ) return 1; - return 0; -} - - -struct cvec { - struct rvec vec; - float na; - float nb; - float nc; - float fom; -}; - - -static int same_vector(struct cvec a, struct cvec b) -{ - if ( a.na != b.na ) return 0; - if ( a.nb != b.nb ) return 0; - if ( a.nc != b.nc ) return 0; - return 1; -} - - -/* Attempt to make 'cell' fit into 'template' somehow */ -UnitCell *match_cell(UnitCell *cell_in, UnitCell *template_in, int verbose, - const float *tols, int reduce) -{ - signed int n1l, n2l, n3l; - double asx, asy, asz; - double bsx, bsy, bsz; - double csx, csy, csz; - int i, j; - double lengths[3]; - double angles[3]; - struct cvec *cand[3]; - UnitCell *new_cell = NULL; - float best_fom = +999999999.9; /* Large number.. */ - int ncand[3] = {0,0,0}; - signed int ilow, ihigh; - float angtol = deg2rad(tols[3]); - UnitCell *cell; - UnitCell *template; - IntegerMatrix *centering; - UnitCell *new_cell_trans; - - /* "Un-center" the template unit cell to make the comparison easier */ - template = uncenter_cell(template_in, ¢ering, NULL); - if ( template == NULL ) return NULL; - - /* The candidate cell is also uncentered, because it might be centered - * if it came from (e.g.) MOSFLM */ - cell = uncenter_cell(cell_in, NULL, NULL); - if ( cell == NULL ) return NULL; - - if ( cell_get_reciprocal(template, &asx, &asy, &asz, - &bsx, &bsy, &bsz, - &csx, &csy, &csz) ) { - ERROR("Couldn't get reciprocal cell for template.\n"); - cell_free(template); - cell_free(cell); - intmat_free(centering); - return NULL; - } - - lengths[0] = modulus(asx, asy, asz); - lengths[1] = modulus(bsx, bsy, bsz); - lengths[2] = modulus(csx, csy, csz); - - angles[0] = angle_between(bsx, bsy, bsz, csx, csy, csz); - angles[1] = angle_between(asx, asy, asz, csx, csy, csz); - angles[2] = angle_between(asx, asy, asz, bsx, bsy, bsz); - - cand[0] = malloc(MAX_CAND*sizeof(struct cvec)); - cand[1] = malloc(MAX_CAND*sizeof(struct cvec)); - cand[2] = malloc(MAX_CAND*sizeof(struct cvec)); - - if ( cell_get_reciprocal(cell, &asx, &asy, &asz, - &bsx, &bsy, &bsz, - &csx, &csy, &csz) ) { - ERROR("Couldn't get reciprocal cell.\n"); - cell_free(template); - cell_free(cell); - intmat_free(centering); - return NULL; - } - - if ( reduce ) { - ilow = -2; ihigh = 4; - } else { - ilow = 0; ihigh = 1; - } - - /* Negative values mean 1/n, positive means n, zero means zero */ - for ( n1l=ilow; n1l<=ihigh; n1l++ ) { - for ( n2l=ilow; n2l<=ihigh; n2l++ ) { - for ( n3l=ilow; n3l<=ihigh; n3l++ ) { - - float n1, n2, n3; - signed int b1, b2, b3; - - n1 = (n1l>=0) ? (n1l) : (1.0/n1l); - n2 = (n2l>=0) ? (n2l) : (1.0/n2l); - n3 = (n3l>=0) ? (n3l) : (1.0/n3l); - - if ( !reduce ) { - if ( n1l + n2l + n3l > 1 ) continue; - } - - /* 'bit' values can be +1 or -1 */ - for ( b1=-1; b1<=1; b1+=2 ) { - for ( b2=-1; b2<=1; b2+=2 ) { - for ( b3=-1; b3<=1; b3+=2 ) { - - double tx, ty, tz; - double tlen; - int i; - - n1 *= b1; n2 *= b2; n3 *= b3; - - tx = n1*asx + n2*bsx + n3*csx; - ty = n1*asy + n2*bsy + n3*csy; - tz = n1*asz + n2*bsz + n3*csz; - tlen = modulus(tx, ty, tz); - - /* Test modulus for agreement with moduli of template */ - for ( i=0; i<3; i++ ) { - - if ( !within_tolerance(lengths[i], tlen, - tols[i]) ) - { - continue; - } - - if ( ncand[i] == MAX_CAND ) { - ERROR("Too many cell candidates - "); - ERROR("consider tightening the unit "); - ERROR("cell tolerances.\n"); - } else { - - double fom; - - fom = fabs(lengths[i] - tlen); - - cand[i][ncand[i]].vec.u = tx; - cand[i][ncand[i]].vec.v = ty; - cand[i][ncand[i]].vec.w = tz; - cand[i][ncand[i]].na = n1; - cand[i][ncand[i]].nb = n2; - cand[i][ncand[i]].nc = n3; - cand[i][ncand[i]].fom = fom; - - ncand[i]++; - - } - } - - } - } - } - - } - } - } - - if ( verbose ) { - STATUS("Candidates: %i %i %i\n", ncand[0], ncand[1], ncand[2]); - } - - for ( i=0; i angtol ) continue; - - fom1 = fabs(ang - angles[2]); - - for ( k=0; k angtol ) continue; - - fom2 = fom1 + fabs(ang - angles[1]); - - /* Finally, the angle between the current candidate for - * axis 1 and the kth candidate for axis 2 */ - ang = angle_between(cand[1][j].vec.u, cand[1][j].vec.v, - cand[1][j].vec.w, cand[2][k].vec.u, - cand[2][k].vec.v, cand[2][k].vec.w); - - /* ... it should be angle 0 ... */ - if ( fabs(ang - angles[0]) > angtol ) continue; - - /* Unit cell must be right-handed */ - if ( !right_handed_vec(cand[0][i].vec, cand[1][j].vec, - cand[2][k].vec) ) continue; - - fom3 = fom2 + fabs(ang - angles[0]); - fom3 += LWEIGHT * (cand[0][i].fom + cand[1][j].fom - + cand[2][k].fom); - - if ( fom3 < best_fom ) { - if ( new_cell != NULL ) free(new_cell); - new_cell = cell_new_from_reciprocal_axes( - cand[0][i].vec, cand[1][j].vec, - cand[2][k].vec); - best_fom = fom3; - } - - } - - } - } - - free(cand[0]); - free(cand[1]); - free(cand[2]); - - cell_free(cell); - - /* Reverse the de-centering transformation */ - if ( new_cell != NULL ) { - - new_cell_trans = cell_transform_intmat(new_cell, centering); - cell_free(new_cell); - cell_set_lattice_type(new_cell_trans, - cell_get_lattice_type(template_in)); - cell_set_centering(new_cell_trans, - cell_get_centering(template_in)); - cell_set_unique_axis(new_cell_trans, - cell_get_unique_axis(template_in)); - - cell_free(template); - intmat_free(centering); - - return new_cell_trans; - - } else { - cell_free(template); - intmat_free(centering); - return NULL; - } -} - - -UnitCell *match_cell_ab(UnitCell *cell_in, UnitCell *template_in) -{ - double ax, ay, az; - double bx, by, bz; - double cx, cy, cz; - int i; - double lengths[3]; - int used[3]; - struct rvec real_a, real_b, real_c; - struct rvec params[3]; - double alen, blen; - float ltl = 5.0; /* percent */ - int have_real_a; - int have_real_b; - int have_real_c; - UnitCell *cell; - UnitCell *template; - IntegerMatrix *to_given_cell; - UnitCell *new_cell; - UnitCell *new_cell_trans; - - /* "Un-center" the template unit cell to make the comparison easier */ - template = uncenter_cell(template_in, &to_given_cell, NULL); - - /* The candidate cell is also uncentered, because it might be centered - * if it came from (e.g.) MOSFLM */ - cell = uncenter_cell(cell_in, NULL, NULL); - - /* Get the lengths to match */ - if ( cell_get_cartesian(template, &ax, &ay, &az, - &bx, &by, &bz, - &cx, &cy, &cz) ) - { - ERROR("Couldn't get cell for template.\n"); - return NULL; - } - alen = modulus(ax, ay, az); - blen = modulus(bx, by, bz); - - /* Get the lengths from the cell and turn them into anonymous vectors */ - if ( cell_get_cartesian(cell, &ax, &ay, &az, - &bx, &by, &bz, - &cx, &cy, &cz) ) - { - ERROR("Couldn't get cell.\n"); - return NULL; - } - lengths[0] = modulus(ax, ay, az); - lengths[1] = modulus(bx, by, bz); - lengths[2] = modulus(cx, cy, cz); - used[0] = 0; used[1] = 0; used[2] = 0; - params[0].u = ax; params[0].v = ay; params[0].w = az; - params[1].u = bx; params[1].v = by; params[1].w = bz; - params[2].u = cx; params[2].v = cy; params[2].w = cz; - - real_a.u = 0.0; real_a.v = 0.0; real_a.w = 0.0; - real_b.u = 0.0; real_b.v = 0.0; real_b.w = 0.0; - real_c.u = 0.0; real_c.v = 0.0; real_c.w = 0.0; - - /* Check each vector against a and b */ - have_real_a = 0; - have_real_b = 0; - for ( i=0; i<3; i++ ) { - if ( within_tolerance(lengths[i], alen, ltl) - && !used[i] && !have_real_a ) - { - used[i] = 1; - memcpy(&real_a, ¶ms[i], sizeof(struct rvec)); - have_real_a = 1; - } - if ( within_tolerance(lengths[i], blen, ltl) - && !used[i] && !have_real_b ) - { - used[i] = 1; - memcpy(&real_b, ¶ms[i], sizeof(struct rvec)); - have_real_b = 1; - } - } - - /* Have we matched both a and b? */ - if ( !(have_real_a && have_real_b) ) return NULL; - - /* "c" is "the other one" */ - have_real_c = 0; - for ( i=0; i<3; i++ ) { - if ( !used[i] ) { - memcpy(&real_c, ¶ms[i], sizeof(struct rvec)); - have_real_c = 1; - } - } - - if ( !have_real_c ) { - ERROR("Huh? Couldn't find the third vector.\n"); - ERROR("Matches: %i %i %i\n", used[0], used[1], used[2]); - return NULL; - } - - /* Flip c if not right-handed */ - if ( !right_handed_vec(real_a, real_b, real_c) ) { - real_c.u = -real_c.u; - real_c.v = -real_c.v; - real_c.w = -real_c.w; - } - - new_cell = cell_new_from_direct_axes(real_a, real_b, real_c); - - /* Reverse the de-centering transformation */ - new_cell_trans = cell_transform_intmat_inverse(new_cell, to_given_cell); - cell_free(new_cell); - cell_set_lattice_type(new_cell, cell_get_lattice_type(template_in)); - cell_set_centering(new_cell, cell_get_centering(template_in)); - cell_set_unique_axis(new_cell, cell_get_unique_axis(template_in)); - - return new_cell_trans; -} - - /* Return sin(theta)/lambda = 1/2d. Multiply by two if you want 1/d */ double resolution(UnitCell *cell, signed int h, signed int k, signed int l) { diff --git a/libcrystfel/src/cell-utils.h b/libcrystfel/src/cell-utils.h index bb1e5cd0..5b705247 100644 --- a/libcrystfel/src/cell-utils.h +++ b/libcrystfel/src/cell-utils.h @@ -59,11 +59,6 @@ extern UnitCell *rotate_cell(UnitCell *in, double omega, double phi, extern void cell_print(UnitCell *cell); extern void cell_print_full(UnitCell *cell); -extern UnitCell *match_cell(UnitCell *cell, UnitCell *tempcell, int verbose, - const float *ltl, int reduce); - -extern UnitCell *match_cell_ab(UnitCell *cell, UnitCell *tempcell); - extern UnitCell *load_cell_from_pdb(const char *filename); extern UnitCell *load_cell_from_file(const char *filename); extern void write_cell(UnitCell *cell, FILE *fh); diff --git a/tests/cellcompare_check.c b/tests/cellcompare_check.c index 107702c4..7aa59fe4 100644 --- a/tests/cellcompare_check.c +++ b/tests/cellcompare_check.c @@ -217,7 +217,7 @@ static void yaro_test() UnitCell *cell; UnitCell *reference; UnitCell *cmatch; - float tols[] = {5, 5, 5, 1.5}; + //float tols[] = {5, 5, 5, 1.5}; double dtols[] = {0.05, 0.05, 0.05, deg2rad(5.0), deg2rad(5.0), deg2rad(5.0)}; cell = cell_new_from_parameters(63.24e-10, 63.94e-10, 64.61e-10, @@ -236,10 +236,10 @@ static void yaro_test() cell_print(cell); STATUS("The reference:\n"); cell_print(reference); - cmatch = match_cell(cell, reference, 0, tols, 1); - STATUS("The match:\n"); - cell_print(cmatch); - cell_free(cmatch); + //cmatch = match_cell(cell, reference, 0, tols, 1); + //STATUS("The match:\n"); + //cell_print(cmatch); + //cell_free(cmatch); RationalMatrix *m = NULL; cmatch = compare_reindexed_cell_parameters(cell, reference, dtols, &m); -- cgit v1.2.3