aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2011-05-16 17:46:09 +0200
committerThomas White <taw@physics.org>2012-02-22 15:27:27 +0100
commitb90ec67959ca72cab7a9339600eeacba37370fa8 (patch)
treeba048d10d38ccd81a745b5db8a2ab33a76871ed0
parentc8057fc5efc6275ac05351d1ea79c49a295efc05 (diff)
Add partial matching (against 'a' and 'b' only)
-rw-r--r--src/cell.c111
-rw-r--r--src/cell.h2
-rw-r--r--src/index.c13
-rw-r--r--src/index.h3
-rw-r--r--src/indexamajig.c45
5 files changed, 155 insertions, 19 deletions
diff --git a/src/cell.c b/src/cell.c
index 144d9140..b1579251 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -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)
{
diff --git a/src/cell.h b/src/cell.h
index fc283837..23b28da5 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -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 */