indexamajig: Add other cell reduction possibilities, and remove --no-match

5 files changed, 79 insertions, 11 deletions

diff --git a/src/cell.c b/src/cell.c
index b6accd58..7002f9fa 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -759,6 +759,28 @@ UnitCell *match_cell(UnitCell *cell, UnitCell *template, int verbose)
 }
 
 
+int cells_similar(UnitCell *cell1, UnitCell *cell2)
+{
+	double a1, b1, c1, al1, be1, ga1;
+	double a2, b2, c2, al2, be2, ga2;
+	double ltl = 5.0;     /* percent */
+	double angtol = deg2rad(1.5);
+
+	cell_get_parameters(cell1, &a1, &b1, &c1, &al1, &be1, &ga1);
+	cell_get_parameters(cell2, &a2, &b2, &c2, &al2, &be2, &ga2);
+
+	if ( !within_tolerance(a1, a2, ltl) ) return 0;
+	if ( !within_tolerance(b1, b2, ltl) ) return 0;
+	if ( !within_tolerance(c1, c2, ltl) ) return 0;
+
+	if ( fabs(al1-al2) > angtol ) return 0;
+	if ( fabs(be1-be2) > angtol ) return 0;
+	if ( fabs(ga1-ga2) > angtol ) return 0;
+
+	return 1;
+}
+
+
 /* 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 a4d32bf8..d3edb503 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -76,6 +76,7 @@ extern double resolution(UnitCell *cell,
 extern void cell_print(UnitCell *cell);
 
 extern UnitCell *match_cell(UnitCell *cell, UnitCell *template, int verbose);
+extern int cells_similar(UnitCell *c1, UnitCell *c2);
 
 extern UnitCell *load_cell_from_pdb(const char *filename);
 
diff --git a/src/index.c b/src/index.c
index 0c656042..f285cccb 100644
--- a/src/index.c
+++ b/src/index.c
@@ -124,7 +124,7 @@ void map_all_peaks(struct image *image)
 
 
 void index_pattern(struct image *image, UnitCell *cell, IndexingMethod indm,
-                   int no_match, int verbose, IndexingPrivate *ipriv)
+                   int cellr, int verbose, IndexingPrivate *ipriv)
 {
 	int i;
 
@@ -150,7 +150,7 @@ void index_pattern(struct image *image, UnitCell *cell, IndexingMethod indm,
 		return;
 	}
 
-	if ( no_match || (indm == INDEXING_TEMPLATE) ) {
+	if ( (cellr == CELLR_NONE) || (indm == INDEXING_TEMPLATE) ) {
 		image->indexed_cell = image->candidate_cells[0];
 		if ( verbose ) {
 			STATUS("--------------------\n");
@@ -172,7 +172,22 @@ void index_pattern(struct image *image, UnitCell *cell, IndexingMethod indm,
 			STATUS("--------------------\n");
 		}
 
-		new_cell = match_cell(image->candidate_cells[i], cell, verbose);
+		/* Match or reduce the cell as appropriate */
+		switch ( cellr ) {
+		case CELLR_NONE :
+			/* Never happens */
+			break;
+		case CELLR_REDUCE :
+			new_cell = match_cell(image->candidate_cells[i],
+			                      cell, verbose);
+			break;
+		case CELLR_COMPARE :
+			if ( cells_similar(image->candidate_cells[i], cell) ) {
+				new_cell = image->candidate_cells[i];
+			}
+			break;
+		}
+
 		image->indexed_cell = new_cell;
 		if ( new_cell != NULL ) {
 			STATUS("Matched on attempt %i.\n", i);
diff --git a/src/index.h b/src/index.h
index 0b740995..81549e00 100644
--- a/src/index.h
+++ b/src/index.h
@@ -22,6 +22,8 @@
 #include "image.h"
 #include "detector.h"
 
+
+/* Indexing methods */
 typedef enum {
 	INDEXING_NONE,
 	INDEXING_DIRAX,
@@ -29,6 +31,14 @@ typedef enum {
 } IndexingMethod;
 
 
+/* Cell reduction methods */
+enum {
+	CELLR_NONE,
+	CELLR_REDUCE,
+	CELLR_COMPARE
+};
+
+
 typedef struct _indexingprivate IndexingPrivate;
 
 extern IndexingPrivate *prepare_indexing(IndexingMethod indm, UnitCell *cell,
@@ -39,7 +49,7 @@ extern IndexingPrivate *prepare_indexing(IndexingMethod indm, UnitCell *cell,
 extern void map_all_peaks(struct image *image);
 
 extern void index_pattern(struct image *image, UnitCell *cell,
-                          IndexingMethod indm, int no_match, int verbose,
+                          IndexingMethod indm, int cellr, int verbose,
                           IndexingPrivate *priv);
 
 extern void cleanup_indexing(IndexingPrivate *priv);
diff --git a/src/indexamajig.c b/src/indexamajig.c
index 324e090f..2b733f04 100644
--- a/src/indexamajig.c
+++ b/src/indexamajig.c
@@ -59,7 +59,6 @@ struct static_index_args
 	int config_nearbragg;
 	int config_gpu;
 	int config_simulate;
-	int config_nomatch;
 	int config_polar;
 	int config_sanity;
 	int config_satcorr;
@@ -72,6 +71,7 @@ struct static_index_args
 	const double *intensities;
 	struct gpu_context *gctx;
 	int peaks;
+	int cellr;
 	double nominal_photon_energy;
 
 	/* Output stream */
@@ -313,7 +313,6 @@ static void process_image(void *pp, int cookie)
 	int config_nearbragg = pargs->static_args.config_nearbragg;
 	int config_gpu = pargs->static_args.config_gpu;
 	int config_simulate = pargs->static_args.config_simulate;
-	int config_nomatch = pargs->static_args.config_nomatch;
 	int config_polar = pargs->static_args.config_polar;
 	IndexingMethod indm = pargs->static_args.indm;
 	const double *intensities = pargs->static_args.intensities;
@@ -397,7 +396,7 @@ static void process_image(void *pp, int cookie)
 
 	/* Calculate orientation matrix (by magic) */
 	if ( config_writedrx || (indm != INDEXING_NONE) ) {
-		index_pattern(&image, cell, indm, config_nomatch,
+		index_pattern(&image, cell, indm, pargs->static_args.cellr,
 		              config_verbose, pargs->static_args.ipriv);
 	}
 
@@ -522,7 +521,6 @@ int main(int argc, char *argv[])
 	int config_simulate = 0;
 	int config_cmfilter = 0;
 	int config_noisefilter = 0;
-	int config_nomatch = 0;
 	int config_gpu = 0;
 	int config_verbose = 0;
 	int config_alternate = 0;
@@ -543,6 +541,8 @@ int main(int argc, char *argv[])
 	char *pdb = NULL;
 	char *prefix = NULL;
 	char *speaks = NULL;
+	char *scellr = NULL;
+	int cellr;
 	int peaks;
 	int nthreads = 1;
 	int i;
@@ -564,6 +564,7 @@ int main(int argc, char *argv[])
 		{"no-index",           0, &config_noindex,     1},
 		{"dump-peaks",         0, &config_dumpfound,   1},
 		{"peaks",              1, NULL,                2},
+		{"cell-reduction",     1, NULL,                3},
 		{"near-bragg",         0, &config_nearbragg,   1},
 		{"write-drx",          0, &config_writedrx,    1},
 		{"indexing",           1, NULL,               'z'},
@@ -572,7 +573,6 @@ int main(int argc, char *argv[])
 		{"simulate",           0, &config_simulate,    1},
 		{"filter-cm",          0, &config_cmfilter,    1},
 		{"filter-noise",       0, &config_noisefilter, 1},
-		{"no-match",           0, &config_nomatch,     1},
 		{"verbose",            0, &config_verbose,     1},
 		{"alternate",          0, &config_alternate,   1},
 		{"intensities",        1, NULL,               'q'},
@@ -647,6 +647,10 @@ int main(int argc, char *argv[])
 			speaks = strdup(optarg);
 			break;
 
+		case 3 :
+			scellr = strdup(optarg);
+			break;
+
 		case 0 :
 			break;
 
@@ -743,6 +747,22 @@ int main(int argc, char *argv[])
 	}
 	free(indm_str);
 
+	if ( scellr == NULL ) {
+		STATUS("You didn't specify a cell reduction method, so I'm"
+		       " going to use 'reduce'.\n");
+		cellr = CELLR_REDUCE;
+	} else if ( strcmp(scellr, "none") == 0 ) {
+		cellr = CELLR_NONE;
+	} else if ( strcmp(scellr, "reduce") == 0) {
+		cellr = CELLR_REDUCE;
+	} else if ( strcmp(scellr, "compare") == 0) {
+		cellr = CELLR_COMPARE;
+	} else {
+		ERROR("Unrecognised cell reduction method '%s'\n", scellr);
+		return 1;
+	}
+	free(scellr);  /* free(NULL) is OK. */
+
 	if ( geometry == NULL ) {
 		ERROR("You need to specify a geometry file with --geometry\n");
 		return 1;
@@ -755,7 +775,7 @@ int main(int argc, char *argv[])
 	}
 	free(geometry);
 
-	if ( (!config_nomatch) || (indm == INDEXING_TEMPLATE) ) {
+	if ( (cellr != CELLR_NONE) || (indm == INDEXING_TEMPLATE) ) {
 		cell = load_cell_from_pdb(pdb);
 		if ( cell == NULL ) {
 			ERROR("Couldn't read unit cell (from %s)\n", pdb);
@@ -814,12 +834,12 @@ int main(int argc, char *argv[])
 	qargs.static_args.config_nearbragg = config_nearbragg;
 	qargs.static_args.config_gpu = config_gpu;
 	qargs.static_args.config_simulate = config_simulate;
-	qargs.static_args.config_nomatch = config_nomatch;
 	qargs.static_args.config_polar = config_polar;
 	qargs.static_args.config_sanity = config_sanity;
 	qargs.static_args.config_satcorr = config_satcorr;
 	qargs.static_args.config_sa = config_sa;
 	qargs.static_args.config_closer = config_closer;
+	qargs.static_args.cellr = cellr;
 	qargs.static_args.threshold = threshold;
 	qargs.static_args.det = det;
 	qargs.static_args.indm = indm;