diff options
| -rw-r--r-- | src/cell.c | 111 | ||||
| -rw-r--r-- | src/cell.h | 2 | ||||
| -rw-r--r-- | src/index.c | 13 | ||||
| -rw-r--r-- | src/index.h | 3 | ||||
| -rw-r--r-- | src/indexamajig.c | 45 |
5 files changed, 155 insertions, 19 deletions
@@ -844,6 +844,117 @@ UnitCell *match_cell(UnitCell *cell, UnitCell *template, int verbose, } +/* FIXME: Unify with proper match_cell(), or work out if it's even possible. + * FIXME: Negative vectors are allowable, but keep a right-handed unit cell. + * FIXME: Make sure unit cell is right handed. */ +UnitCell *match_cell_ab(UnitCell *cell, UnitCell *template) +{ + double ax, ay, az; + double bx, by, bz; + double cx, cy, cz; + int i; + double lengths[3]; + int used[3]; + double real_ax, real_ay, real_az; + double real_bx, real_by, real_bz; + double real_cx, real_cy, real_cz; + double params[3][3]; + double alen, blen, clen; + float ltl = 5.0; /* percent */ + int have_real_a; + int have_real_b; + int have_real_c; + UnitCell *new_cell; + + /* 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); + clen = modulus(cx, cy, cz); + + /* 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][0] = ax; params[0][1] = ay; params[0][2] = az; + params[1][0] = bx; params[1][1] = by; params[1][2] = bz; + params[2][0] = cx; params[2][1] = cy; params[2][2] = cz; + + real_ax = 0.0; real_ay = 0.0; real_az = 0.0; + real_bx = 0.0; real_by = 0.0; real_bz = 0.0; + real_cx = 0.0; real_cy = 0.0; real_cz = 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; + real_ax = params[i][0]; + real_ay = params[i][1]; + real_az = params[i][2]; + have_real_a = 1; + } + if ( within_tolerance(lengths[i], blen, ltl) && !used[i] + && !have_real_b ) + { + used[i] = 1; + real_bx = params[i][0]; + real_by = params[i][1]; + real_bz = params[i][2]; + 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] ) { + real_cx = params[i][0]; + real_cy = params[i][1]; + real_cz = params[i][2]; + 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; + } + + new_cell = cell_new(); + cell_set_cartesian(new_cell, real_ax, real_ay, real_az, + real_bx, real_by, real_bz, + real_cx, real_cy, real_cz); + + return new_cell; +} + + /* 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) { @@ -98,6 +98,8 @@ extern void cell_print(UnitCell *cell); extern UnitCell *match_cell(UnitCell *cell, UnitCell *template, int verbose, int reduce); +extern UnitCell *match_cell_ab(UnitCell *cell, UnitCell *template); + extern UnitCell *load_cell_from_pdb(const char *filename); diff --git a/src/index.c b/src/index.c index 47866a18..2fa7f58d 100644 --- a/src/index.c +++ b/src/index.c @@ -154,6 +154,7 @@ void index_pattern(struct image *image, UnitCell *cell, IndexingMethod *indm, for ( i=0; i<image->ncells; i++ ) { UnitCell *new_cell = NULL; + UnitCell *cand = image->candidate_cells[i]; if ( verbose ) { STATUS("--------------------\n"); @@ -166,16 +167,16 @@ void index_pattern(struct image *image, UnitCell *cell, IndexingMethod *indm, /* Match or reduce the cell as appropriate */ switch ( cellr ) { case CELLR_NONE : - new_cell = cell_new_from_cell(image - ->candidate_cells[i]); + new_cell = cell_new_from_cell(cand); break; case CELLR_REDUCE : - new_cell = match_cell(image->candidate_cells[i], - cell, verbose, 1); + new_cell = match_cell(cand, cell, verbose, 1); break; case CELLR_COMPARE : - new_cell = match_cell(image->candidate_cells[i], - cell, verbose, 0); + new_cell = match_cell(cand, cell, verbose, 0); + break; + case CELLR_COMPARE_AB : + new_cell = match_cell_ab(cand, cell); break; } diff --git a/src/index.h b/src/index.h index 005f0d93..03e246f4 100644 --- a/src/index.h +++ b/src/index.h @@ -35,7 +35,8 @@ typedef enum { enum { CELLR_NONE, CELLR_REDUCE, - CELLR_COMPARE + CELLR_COMPARE, + CELLR_COMPARE_AB, }; diff --git a/src/indexamajig.c b/src/indexamajig.c index 569a41bd..ed6b25c2 100644 --- a/src/indexamajig.c +++ b/src/indexamajig.c @@ -168,6 +168,7 @@ static void show_help(const char *s) " reduce : full cell reduction.\n" " compare : match by at most changing the order of\n" " the indices.\n" +" compare_ab : compare 'a' and 'b' lengths only.\n" " --filter-cm Perform common-mode noise subtraction on images\n" " before proceeding. Intensities will be extracted\n" " from the image as it is after this processing.\n" @@ -464,6 +465,29 @@ static void finalise_image(void *qp, void *pp) } +static int parse_cell_reduction(const char *scellr, int *err, + int *reduction_needs_cell) +{ + if ( strcmp(scellr, "none") == 0 ) { + *reduction_needs_cell = 0; + return CELLR_NONE; + } else if ( strcmp(scellr, "reduce") == 0) { + *reduction_needs_cell = 1; + return CELLR_REDUCE; + } else if ( strcmp(scellr, "compare") == 0) { + *reduction_needs_cell = 1; + return CELLR_COMPARE; + } else if ( strcmp(scellr, "compare_ab") == 0) { + *reduction_needs_cell = 1; + return CELLR_COMPARE_AB; + } else { + *err = 1; + *reduction_needs_cell = 0; + return CELLR_NONE; + } +} + + int main(int argc, char *argv[]) { int c; @@ -759,20 +783,17 @@ int main(int argc, char *argv[]) " I'm going to use 'reduce'.\n"); cellr = CELLR_REDUCE; reduction_needs_cell = 1; - } else if ( strcmp(scellr, "none") == 0 ) { - cellr = CELLR_NONE; - } else if ( strcmp(scellr, "reduce") == 0) { - cellr = CELLR_REDUCE; - reduction_needs_cell = 1; - } else if ( strcmp(scellr, "compare") == 0) { - cellr = CELLR_COMPARE; - reduction_needs_cell = 1; } else { - ERROR("Unrecognised cell reduction method '%s'\n", - scellr); - return 1; + int err; + cellr = parse_cell_reduction(scellr, &err, + &reduction_needs_cell); + if ( err ) { + ERROR("Unrecognised cell reduction '%s'\n", + scellr); + return 1; + } + free(scellr); } - free(scellr); /* free(NULL) is OK. */ } /* No indexing -> no reduction */ |