Add partial matching (against 'a' and 'b' only)

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 */