Remove match_cell and match_cell_ab

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, &centering, 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<ncand[0]; i++ ) {
-	for ( j=0; j<ncand[1]; j++ ) {
-
-		double ang;
-		int k;
-		float fom1;
-
-		if ( same_vector(cand[0][i], cand[1][j]) ) continue;
-
-		/* Measure the angle between the ith candidate for axis 0
-		 * and the jth candidate for axis 1 */
-		ang = angle_between(cand[0][i].vec.u, cand[0][i].vec.v,
-		                    cand[0][i].vec.w, cand[1][j].vec.u,
-		                    cand[1][j].vec.v, cand[1][j].vec.w);
-
-		/* Angle between axes 0 and 1 should be angle 2 */
-		if ( fabs(ang - angles[2]) > angtol ) continue;
-
-		fom1 = fabs(ang - angles[2]);
-
-		for ( k=0; k<ncand[2]; k++ ) {
-
-			float fom2, fom3;
-
-			if ( same_vector(cand[1][j], cand[2][k]) ) continue;
-
-			/* Measure the angle between the current candidate for
-			 * axis 0 and the kth candidate for axis 2 */
-			ang = angle_between(cand[0][i].vec.u, cand[0][i].vec.v,
-			                    cand[0][i].vec.w, cand[2][k].vec.u,
-			                    cand[2][k].vec.v, cand[2][k].vec.w);
-
-			/* ... it should be angle 1 ... */
-			if ( fabs(ang - angles[1]) > 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, &params[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, &params[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, &params[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);