Now tries the shortest vectors first.

1 files changed, 417 insertions, 398 deletions

diff --git a/libcrystfel/src/taketwo.c b/libcrystfel/src/taketwo.c
index b92743b4..a1e51c2a 100644
--- a/libcrystfel/src/taketwo.c
+++ b/libcrystfel/src/taketwo.c
@@ -73,7 +73,7 @@ struct taketwo_private
 
 
 /* Maximum distance between two rlp sizes to consider info for indexing */
-#define MAX_RECIP_DISTANCE (0.12*1e10)
+#define MAX_RECIP_DISTANCE (0.20*1e10)
 
 /* Tolerance for two lengths in reciprocal space to be considered the same */
 #define RECIP_TOLERANCE (0.001*1e10)
@@ -88,7 +88,7 @@ struct taketwo_private
 #define MAX_DEAD_ENDS (5)
 
 /* Tolerance for two angles to be considered the same */
-#define ANGLE_TOLERANCE (deg2rad(1.0))
+#define ANGLE_TOLERANCE (deg2rad(3.0))
 
 /* Tolerance for rot_mats_are_similar */
 #define TRACE_TOLERANCE (deg2rad(3.0))
@@ -129,8 +129,8 @@ static struct rvec new_rvec(double new_u, double new_v, double new_w)
 static struct rvec diff_vec(struct rvec from, struct rvec to)
 {
 	struct rvec diff = new_rvec(to.u - from.u,
-	                     to.v - from.v,
-	                     to.w - from.w);
+				    to.v - from.v,
+				    to.w - from.w);
 
 	return diff;
 }
@@ -151,10 +151,10 @@ static double rvec_length(struct rvec vec)
 
 static void normalise_rvec(struct rvec *vec)
 {
-        double length = rvec_length(*vec);
-        vec->u /= length;
-        vec->v /= length;
-        vec->w /= length;
+	double length = rvec_length(*vec);
+	vec->u /= length;
+	vec->v /= length;
+	vec->w /= length;
 }
 
 
@@ -193,8 +193,8 @@ static struct rvec rvec_cross(struct rvec a, struct rvec b)
 
 static void show_rvec(struct rvec r2)
 {
-        struct rvec r = r2;
-        normalise_rvec(&r);
+	struct rvec r = r2;
+	normalise_rvec(&r);
 	STATUS("[ %.3f %.3f %.3f ]\n", r.u, r.v, r.w);
 }
 
@@ -229,7 +229,7 @@ static gsl_matrix *rotation_around_axis(struct rvec c, double th)
  * that @result has already been allocated. Will upload the maths to the
  * shared Google drive. */
 static gsl_matrix *closest_rot_mat(struct rvec vec1, struct rvec vec2,
-                                   struct rvec axis)
+				   struct rvec axis)
 {
 	/* Let's have unit vectors */
 	normalise_rvec(&vec1);
@@ -244,8 +244,8 @@ static gsl_matrix *closest_rot_mat(struct rvec vec1, struct rvec vec2,
 
 	/* Components in handwritten maths online when I upload it */
 	double A = a*(p*x*x - p + x*y*q + x*z*r) +
-	           b*(p*x*y + q*y*y - q + r*y*z) +
-	           c*(p*x*z + q*y*z + r*z*z - r);
+	b*(p*x*y + q*y*y - q + r*y*z) +
+	c*(p*x*z + q*y*z + r*z*z - r);
 
 	double B = a*(y*r - z*q) + b*(p*z - r*x) + c*(q*x - p*y);
 
@@ -256,26 +256,26 @@ static gsl_matrix *closest_rot_mat(struct rvec vec1, struct rvec vec2,
 	 * and we need to work out which one. This could potentially be
 	 * simplified - do we really need so many cos/sins? maybe check
 	 * the 2nd derivative instead? */
-        double cc = cos(theta);
-        double C = 1 - cc;
-        double s = sin(theta);
-        double occ = -cc;
-        double oC = 1 - occ;
-        double os = -s;
+	double cc = cos(theta);
+	double C = 1 - cc;
+	double s = sin(theta);
+	double occ = -cc;
+	double oC = 1 - occ;
+	double os = -s;
 
-        double pPrime = (x*x*C+cc)*p + (x*y*C-z*s)*q + (x*z*C+y*s)*r;
-        double qPrime = (y*x*C+z*s)*p + (y*y*C+cc)*q + (y*z*C-x*s)*r;
-        double rPrime = (z*x*C-y*s)*p + (z*y*C+x*s)*q + (z*z*C+cc)*r;
+	double pPrime = (x*x*C+cc)*p + (x*y*C-z*s)*q + (x*z*C+y*s)*r;
+	double qPrime = (y*x*C+z*s)*p + (y*y*C+cc)*q + (y*z*C-x*s)*r;
+	double rPrime = (z*x*C-y*s)*p + (z*y*C+x*s)*q + (z*z*C+cc)*r;
 
-        double pDbPrime = (x*x*oC+occ)*p + (x*y*oC-z*os)*q + (x*z*oC+y*os)*r;
-        double qDbPrime = (y*x*oC+z*os)*p + (y*y*oC+occ)*q + (y*z*oC-x*os)*r;
-        double rDbPrime = (z*x*oC-y*os)*p + (z*y*oC+x*os)*q + (z*z*oC+occ)*r;
+	double pDbPrime = (x*x*oC+occ)*p + (x*y*oC-z*os)*q + (x*z*oC+y*os)*r;
+	double qDbPrime = (y*x*oC+z*os)*p + (y*y*oC+occ)*q + (y*z*oC-x*os)*r;
+	double rDbPrime = (z*x*oC-y*os)*p + (z*y*oC+x*os)*q + (z*z*oC+occ)*r;
 
-        double cosAlpha = pPrime * a + qPrime * b + rPrime * c;
-        double cosAlphaOther = pDbPrime * a + qDbPrime * b + rDbPrime * c;
+	double cosAlpha = pPrime * a + qPrime * b + rPrime * c;
+	double cosAlphaOther = pDbPrime * a + qDbPrime * b + rDbPrime * c;
 
-        int addPi = (cosAlphaOther > cosAlpha);
-        double bestAngle = theta + addPi * M_PI;
+	int addPi = (cosAlphaOther > cosAlpha);
+	double bestAngle = theta + addPi * M_PI;
 
 	/* Return an identity matrix which has been rotated by
 	 * theta around "axis" */
@@ -297,7 +297,7 @@ static double matrix_trace(gsl_matrix *a)
 
 
 static int rot_mats_are_similar(gsl_matrix *rot1, gsl_matrix *rot2,
-                                double *score)
+				double *score)
 {
 	double tr;
 	gsl_matrix *sub;
@@ -312,12 +312,12 @@ static int rot_mats_are_similar(gsl_matrix *rot1, gsl_matrix *rot2,
 
 	tr = matrix_trace(mul);
 	if (score != NULL) *score = tr;
-	
+
 	gsl_matrix_free(mul);
 	gsl_matrix_free(sub);
-        double max = sqrt(4.0*(1.0-cos(TRACE_TOLERANCE)));
+	double max = sqrt(4.0*(1.0-cos(TRACE_TOLERANCE)));
 
-  	return (tr < max);
+	return (tr < max);
 }
 
 
@@ -355,7 +355,7 @@ struct rvec gsl_to_rvec(gsl_vector *a)
  * intensive on the number crunching side so simple angle checks are used
  * to 'pre-scan' vectors beforehand. */
 static gsl_matrix *generate_rot_mat(struct rvec obs1, struct rvec obs2,
-                                    struct rvec cell1, struct rvec cell2)
+				    struct rvec cell1, struct rvec cell2)
 {
 	gsl_matrix *rotateSpotDiffMatrix;
 	gsl_matrix *secondTwizzleMatrix;
@@ -363,33 +363,33 @@ static gsl_matrix *generate_rot_mat(struct rvec obs1, struct rvec obs2,
 	gsl_vector *cell2v = rvec_to_gsl(cell2);
 	gsl_vector *cell2vr = gsl_vector_calloc(3);
 
-        normalise_rvec(&obs1);
-        normalise_rvec(&obs2);
-        normalise_rvec(&cell1);
-        normalise_rvec(&cell2);
+	normalise_rvec(&obs1);
+	normalise_rvec(&obs2);
+	normalise_rvec(&cell1);
+	normalise_rvec(&cell2);
 
-        /* Rotate reciprocal space so that the first simulated vector lines up
-         * with the observed vector. */
-        rotateSpotDiffMatrix = rotation_between_vectors(cell1, obs1);
+	/* Rotate reciprocal space so that the first simulated vector lines up
+	 * with the observed vector. */
+	rotateSpotDiffMatrix = rotation_between_vectors(cell1, obs1);
 
-        normalise_rvec(&obs1);
+	normalise_rvec(&obs1);
 
-        /* Multiply cell2 by rotateSpotDiffMatrix --> cell2vr */
-        gsl_blas_dgemv(CblasNoTrans, 1.0, rotateSpotDiffMatrix, cell2v,
-                       0.0, cell2vr);
+	/* Multiply cell2 by rotateSpotDiffMatrix --> cell2vr */
+	gsl_blas_dgemv(CblasNoTrans, 1.0, rotateSpotDiffMatrix, cell2v,
+		       0.0, cell2vr);
 
-        /* Now we twirl around the firstAxisUnit until the rotated simulated
-         * vector matches the second observed vector as closely as possible. */
+	/* Now we twirl around the firstAxisUnit until the rotated simulated
+	 * vector matches the second observed vector as closely as possible. */
 
-        secondTwizzleMatrix = closest_rot_mat(gsl_to_rvec(cell2vr), obs2, obs1);
+	secondTwizzleMatrix = closest_rot_mat(gsl_to_rvec(cell2vr), obs2, obs1);
 
 	/* We want to apply the first matrix and then the second matrix,
 	 * so we multiply these. */
 	fullMat = gsl_matrix_calloc(3, 3);
 	gsl_blas_dgemm(CblasTrans, CblasTrans, 1.0,
-	               rotateSpotDiffMatrix, secondTwizzleMatrix, 0.0, fullMat);
+		       rotateSpotDiffMatrix, secondTwizzleMatrix, 0.0, fullMat);
 	gsl_matrix_transpose(fullMat);
-	
+
 	gsl_vector_free(cell2v);
 	gsl_vector_free(cell2vr);
 	gsl_matrix_free(secondTwizzleMatrix);
@@ -403,9 +403,9 @@ static int obs_vecs_share_spot(struct SpotVec *her_obs, struct SpotVec *his_obs)
 {
 	/* FIXME: Disgusting... can I tone this down a bit? */
 	if ( (her_obs->her_rlp == his_obs->her_rlp) ||
-	     (her_obs->her_rlp == his_obs->his_rlp) ||
-	     (her_obs->his_rlp == his_obs->her_rlp) ||
-	     (her_obs->his_rlp == his_obs->his_rlp) ) {
+	    (her_obs->her_rlp == his_obs->his_rlp) ||
+	    (her_obs->his_rlp == his_obs->her_rlp) ||
+	    (her_obs->his_rlp == his_obs->his_rlp) ) {
 		return 1;
 	}
 
@@ -414,7 +414,7 @@ static int obs_vecs_share_spot(struct SpotVec *her_obs, struct SpotVec *his_obs)
 
 
 static int obs_shares_spot_w_array(struct SpotVec *obs_vecs, int test_idx,
-                                   int *members, int num)
+				   int *members, int num)
 {
 	int i;
 	struct SpotVec *her_obs = &obs_vecs[test_idx];
@@ -437,15 +437,15 @@ static int obs_shares_spot_w_array(struct SpotVec *obs_vecs, int test_idx,
  * cosine graph are more sensitive than others, so may be a trade off... or not.
  */
 static int obs_vecs_match_angles(struct SpotVec *her_obs,
-                                 struct SpotVec *his_obs,
+				 struct SpotVec *his_obs,
 				 int **her_match_idxs, int **his_match_idxs,
 				 int *match_count)
 {
-    int i, j;
+	int i, j;
 	*match_count = 0;
 
-       // *her_match_idx = -1;
-       // *his_match_idx = -1;
+	// *her_match_idx = -1;
+	// *his_match_idx = -1;
 
 	/* calculate angle between observed vectors */
 	double obs_angle = rvec_angle(her_obs->obsvec, his_obs->obsvec);
@@ -453,42 +453,43 @@ static int obs_vecs_match_angles(struct SpotVec *her_obs,
 	/* calculate angle between all potential theoretical vectors */
 
 	for ( i=0; i<her_obs->match_num; i++ ) {
-	for ( j=0; j<his_obs->match_num; j++ ) {
+		for ( j=0; j<his_obs->match_num; j++ ) {
 
-		struct rvec *her_match = &her_obs->matches[i];
-		struct rvec *his_match = &his_obs->matches[j];
+			struct rvec *her_match = &her_obs->matches[i];
+			struct rvec *his_match = &his_obs->matches[j];
 
-		double theory_angle = rvec_angle(*her_match, *his_match);
+			double theory_angle = rvec_angle(*her_match,
+							 *his_match);
 
-		/* is this angle a match? */
+			/* is this angle a match? */
 
-		double angle_diff = fabs(theory_angle - obs_angle);
+			double angle_diff = fabs(theory_angle - obs_angle);
 
-                if ( angle_diff < ANGLE_TOLERANCE ) {
-                    size_t new_size = (*match_count + 1) *
-                                        sizeof(int);
-                    if (her_match_idxs && his_match_idxs)
-                    {
-                        /* Reallocate the array to fit in another match */ 
-                        int *temp_hers;
-                        temp_hers = realloc(*her_match_idxs, new_size);
-                        int *temp_his;
-                        temp_his = realloc(*his_match_idxs, new_size);
+			if ( angle_diff < ANGLE_TOLERANCE ) {
+				size_t new_size = (*match_count + 1) *
+				sizeof(int);
+				if (her_match_idxs && his_match_idxs)
+				{
+					/* Reallocate the array to fit in another match */
+					int *temp_hers;
+					temp_hers = realloc(*her_match_idxs, new_size);
+					int *temp_his;
+					temp_his = realloc(*his_match_idxs, new_size);
 
-                        if ( temp_hers == NULL || temp_his == NULL ) {
-                                apologise();
-                        }
+					if ( temp_hers == NULL || temp_his == NULL ) {
+						apologise();
+					}
 
-                        (*her_match_idxs) = temp_hers;
-                        (*his_match_idxs) = temp_his;
+					(*her_match_idxs) = temp_hers;
+					(*his_match_idxs) = temp_his;
 
-                        (*her_match_idxs)[*match_count] = i;
-                        (*his_match_idxs)[*match_count] = j;
-                    }
+					(*her_match_idxs)[*match_count] = i;
+					(*his_match_idxs)[*match_count] = j;
+				}
 
-                (*match_count)++;
-                }
-	}
+				(*match_count)++;
+			}
+		}
 	}
 
 	return (*match_count > 0);
@@ -496,7 +497,7 @@ static int obs_vecs_match_angles(struct SpotVec *her_obs,
 
 
 static int obs_angles_match_array(struct SpotVec *obs_vecs, int test_idx,
-                                  int *obs_members, int *match_members, int num)
+				  int *obs_members, int *match_members, int num)
 {
 	/* note: this is just a preliminary check to reduce unnecessary
 	 * computation later down the line, but is not entirely accurate.
@@ -511,13 +512,13 @@ static int obs_angles_match_array(struct SpotVec *obs_vecs, int test_idx,
 		struct SpotVec *his_obs = &obs_vecs[obs_members[i]];
 
 		/* placeholders, but results are ignored */
-                int match_count = 0;
+		int match_count = 0;
 
 		/* check our test vector matches existing network member */
 		int success = obs_vecs_match_angles(her_obs, his_obs,
-		                      NULL, NULL,
-		                      &match_count);
-		
+						    NULL, NULL,
+						    &match_count);
+
 		if (!success) return 0;
 	}
 
@@ -530,82 +531,82 @@ static int obs_angles_match_array(struct SpotVec *obs_vecs, int test_idx,
  * ------------------------------------------------------------------------*/
 
 static signed int finalise_solution(gsl_matrix *rot, struct SpotVec *obs_vecs,
-                                    int *obs_members, int *match_members,
-                                    int member_num)
+				    int *obs_members, int *match_members,
+				    int member_num)
 {
-        gsl_matrix **rotations = malloc(sizeof(*rotations)* pow(member_num, 2) - member_num);
-        
-        int i, j, count;
-        
-        for ( i=0; i<1; i++ ) {
-        for ( j=0; j<member_num; j++ ) {                
-                if (i == j) continue;
-                struct SpotVec i_vec = obs_vecs[obs_members[i]];
-                struct SpotVec j_vec = obs_vecs[obs_members[j]];
-                
-                struct rvec i_obsvec = i_vec.obsvec;
-                struct rvec j_obsvec = j_vec.obsvec;
-                struct rvec i_cellvec = i_vec.matches[match_members[i]];
-                struct rvec j_cellvec = j_vec.matches[match_members[j]];
-                
-                rotations[count] = generate_rot_mat(i_obsvec, j_obsvec,
-                                                    i_cellvec, j_cellvec);
-
-                count++;
-        }
-        }
-        
-        double min_score = FLT_MAX;
-        int min_rot_index = 0;
-        
-        for (i=0; i<count; i++) {
-                double current_score = 0;
-                for (j=0; j<count; j++) {
-                        double addition;
-                        rot_mats_are_similar(rotations[i], rotations[j],
-                                             &addition);
-                
-                        current_score += addition;
-                }
-                
-                if (current_score < min_score) {
-                        min_score = current_score;
-                        min_rot_index = i;
-                }
-        }
-
-        gsl_matrix_memcpy(rot, rotations[min_rot_index]);
-        
-        for (i=0; i<count; i++) {
-                gsl_matrix_free(rotations[i]);
-        }
-        
-        free(rotations);
-        
-        return 1;
+	gsl_matrix **rotations = malloc(sizeof(*rotations)* pow(member_num, 2) - member_num);
+
+	int i, j, count;
+
+	for ( i=0; i<1; i++ ) {
+		for ( j=0; j<member_num; j++ ) {
+			if (i == j) continue;
+			struct SpotVec i_vec = obs_vecs[obs_members[i]];
+			struct SpotVec j_vec = obs_vecs[obs_members[j]];
+
+			struct rvec i_obsvec = i_vec.obsvec;
+			struct rvec j_obsvec = j_vec.obsvec;
+			struct rvec i_cellvec = i_vec.matches[match_members[i]];
+			struct rvec j_cellvec = j_vec.matches[match_members[j]];
+
+			rotations[count] = generate_rot_mat(i_obsvec, j_obsvec,
+							    i_cellvec, j_cellvec);
+
+			count++;
+		}
+	}
+
+	double min_score = FLT_MAX;
+	int min_rot_index = 0;
+
+	for (i=0; i<count; i++) {
+		double current_score = 0;
+		for (j=0; j<count; j++) {
+			double addition;
+			rot_mats_are_similar(rotations[i], rotations[j],
+					     &addition);
+
+			current_score += addition;
+		}
+
+		if (current_score < min_score) {
+			min_score = current_score;
+			min_rot_index = i;
+		}
+	}
+
+	gsl_matrix_memcpy(rot, rotations[min_rot_index]);
+
+	for (i=0; i<count; i++) {
+		gsl_matrix_free(rotations[i]);
+	}
+
+	free(rotations);
+
+	return 1;
 }
 
 static signed int find_next_index(gsl_matrix *rot, struct SpotVec *obs_vecs,
-                                  int obs_vec_count, int *obs_members,
-                                  int *match_members, int start, int member_num,
-                                  int *match_found, int match_start)
+				  int obs_vec_count, int *obs_members,
+				  int *match_members, int start, int member_num,
+				  int *match_found, int match_start)
 {
 	int i;
- 
+
 	for ( i=start; i<obs_vec_count; i++ ) {
 
 		/* first we check for a shared spot - harshest condition */
 		int shared = obs_shares_spot_w_array(obs_vecs, i, obs_members,
-		                                     member_num);
+						     member_num);
 
 		if ( !shared ) continue;
 
 		/* now we check that angles between all vectors match */
-	/*	int matches = obs_angles_match_array(obs_vecs, i, obs_members,
-		                                     match_members, member_num);
+		/*	int matches = obs_angles_match_array(obs_vecs, i, obs_members,
+		 match_members, member_num);
 
-		if ( !matches ) continue;
-*/
+		 if ( !matches ) continue;
+		 */
 		/* final test: does the corresponding rotation matrix
 		 * match the others? NOTE: have not tested to see if
 		 * every combination of test/existing member has to be checked
@@ -617,55 +618,56 @@ static signed int find_next_index(gsl_matrix *rot, struct SpotVec *obs_vecs,
 		int *member_idx = 0;
 		int *test_idx = 0;
 		int match_num;
-		
+
 		obs_vecs_match_angles(&obs_vecs[obs_members[0]], &obs_vecs[i],
-                                      &member_idx, &test_idx, &match_num);
-		
+				      &member_idx, &test_idx, &match_num);
+
 		/* if ok is set to 0, give up on this vector before
 		 * checking the next value of j */
-                int j, k;
-
-                for ( j=match_start; j<match_num; j++ ) {
-                        int all_ok = 1;
-                        for ( k=0; k<member_num; k++ )
-                        {
-                                gsl_matrix *test_rot;
-                                struct rvec member_match;
-                        
-                                int idx_k = obs_members[k];
-
-                                /* First observed vector and matching theoretical */
-                                member_match = obs_vecs[idx_k].matches[match_members[k]];
-
-                                /* Potential match with another vector */
-                                struct rvec a_match = obs_vecs[i].matches[test_idx[j]];
-  
-                                test_rot = generate_rot_mat(obs_vecs[idx_k].obsvec,
-                                                            obs_vecs[i].obsvec,
-                                                            member_match,
-                                                            a_match);
-
-                                double trace = 0;
-                                int ok = rot_mats_are_similar(rot, test_rot, &trace);
-                                gsl_matrix_free(test_rot);
-                                
-                                if (!ok) {
-                                        all_ok = 0;
-                                        break;
-                                }
-                        
-                        }
-
-                        if (all_ok) {
-                                *match_found = test_idx[j];
-                                free(member_idx);
-                                free(test_idx);
-                                return i;
-                        }
-                }
-                
-                free(member_idx); member_idx = NULL;
-                free(test_idx); member_idx = NULL;
+		int j, k;
+
+		for ( j=match_start; j<match_num; j++ ) {
+			int all_ok = 1;
+			for ( k=0; k<member_num; k++ )
+			{
+				gsl_matrix *test_rot;
+				struct rvec member_match;
+
+				int idx_k = obs_members[k];
+
+				/* First observed vector and matching theoretical */
+				member_match = obs_vecs[idx_k].matches[match_members[k]];
+
+				/* Potential match with another vector */
+				struct rvec a_match = obs_vecs[i].matches[test_idx[j]];
+
+				test_rot = generate_rot_mat(obs_vecs[idx_k].obsvec,
+							    obs_vecs[i].obsvec,
+							    member_match,
+							    a_match);
+
+				double trace = 0;
+				int ok = rot_mats_are_similar(rot, test_rot,
+							      &trace);
+				gsl_matrix_free(test_rot);
+
+				if (!ok) {
+					all_ok = 0;
+					break;
+				}
+
+			}
+
+			if (all_ok) {
+				*match_found = test_idx[j];
+				free(member_idx);
+				free(test_idx);
+				return i;
+			}
+		}
+
+		free(member_idx); member_idx = NULL;
+		free(test_idx); member_idx = NULL;
 
 	}
 
@@ -676,12 +678,12 @@ static signed int find_next_index(gsl_matrix *rot, struct SpotVec *obs_vecs,
 
 
 static int grow_network(gsl_matrix *rot, struct SpotVec *obs_vecs,
-                        int obs_vec_count, int obs_idx1, int obs_idx2,
-                        int match_idx1, int match_idx2, int *max_members)
+			int obs_vec_count, int obs_idx1, int obs_idx2,
+			int match_idx1, int match_idx2, int *max_members)
 {
 	/* indices of members of the self-consistent network of vectors */
 	int obs_members[MAX_NETWORK_MEMBERS];
-        int match_members[MAX_NETWORK_MEMBERS];
+	int match_members[MAX_NETWORK_MEMBERS];
 
 	/* initialise the ones we know already */
 	obs_members[0] = obs_idx1;
@@ -689,7 +691,7 @@ static int grow_network(gsl_matrix *rot, struct SpotVec *obs_vecs,
 	match_members[0] = match_idx1;
 	match_members[1] = match_idx2;
 	int member_num = 2;
-        *max_members = 2;
+	*max_members = 2;
 
 	/* counter for dead ends which must not exceed MAX_DEAD_ENDS
 	 * before it is reset in an additional branch */
@@ -697,169 +699,172 @@ static int grow_network(gsl_matrix *rot, struct SpotVec *obs_vecs,
 
 	/* we can start from after the 2nd observed vector in the seed */
 	int start = obs_idx2 + 1;
-    
+
 	while ( 1 ) {
 
-                if (start > obs_vec_count) {
-                    return 0;
-                }
-
-                int match_found = -1;
-                
-                signed int next_index = find_next_index(rot, obs_vecs,
-                                                 obs_vec_count, obs_members,
-                                                 match_members,
-                                                 start, member_num,
-                                                 &match_found, 0);
-
-                if ( member_num < 2 ) {
-                    return 0;
-                }
-
-                if ( next_index < 0 ) {
-                    /* If there have been too many dead ends, give up
-                         * on indexing altogether.
-                         **/
-                        if ( dead_ends > MAX_DEAD_ENDS ) {
-                            break;
-                        }
-
-                        /* We have not had too many dead ends. Try removing
-                           the last member and continue. */
-                        start++;
-                        member_num--;
-                        dead_ends++;
-
-                        continue;
-                }
-
-                /* we have elongated membership - so reset dead_ends counter */
-        	dead_ends = 0;
-
-                obs_members[member_num] = next_index;
-                match_members[member_num] = match_found;
-
-                start = next_index + 1;
-                member_num++;
-        
-                if (member_num > *max_members) {
-                    *max_members = member_num;
-                }
-        
-                /* If member_num is high enough, we want to return a yes */
-                if ( member_num > NETWORK_MEMBER_THRESHOLD ) break;
+		if (start > obs_vec_count) {
+			return 0;
+		}
+
+		int match_found = -1;
+
+		signed int next_index = find_next_index(rot, obs_vecs,
+							obs_vec_count, obs_members,
+							match_members,
+							start, member_num,
+							&match_found, 0);
+
+		if ( member_num < 2 ) {
+			return 0;
+		}
+
+		if ( next_index < 0 ) {
+			/* If there have been too many dead ends, give up
+			 * on indexing altogether.
+			 **/
+			if ( dead_ends > MAX_DEAD_ENDS ) {
+				break;
+			}
+
+			/* We have not had too many dead ends. Try removing
+			 the last member and continue. */
+			start++;
+			member_num--;
+			dead_ends++;
+
+			continue;
+		}
+
+		/* we have elongated membership - so reset dead_ends counter */
+		dead_ends = 0;
+
+		obs_members[member_num] = next_index;
+		match_members[member_num] = match_found;
+
+		start = next_index + 1;
+		member_num++;
+
+		if (member_num > *max_members) {
+			*max_members = member_num;
+		}
+
+		/* If member_num is high enough, we want to return a yes */
+		if ( member_num > NETWORK_MEMBER_THRESHOLD ) break;
 	}
 
-        finalise_solution(rot, obs_vecs, obs_members,
-                          match_members, member_num);
+	finalise_solution(rot, obs_vecs, obs_members,
+			  match_members, member_num);
 
 	return ( member_num > NETWORK_MEMBER_THRESHOLD );
 }
 
 
 static int start_seed(struct SpotVec *obs_vecs, int obs_vec_count, int i,
-                      int j, int i_match, int j_match, gsl_matrix **rotation,
-                      int *max_members)
+		      int j, int i_match, int j_match, gsl_matrix **rotation,
+		      int *max_members)
 {
-        gsl_matrix *rot_mat;
-
-        rot_mat = generate_rot_mat(obs_vecs[i].obsvec,
-                                   obs_vecs[j].obsvec,
-                                   obs_vecs[i].matches[i_match],
-                                   obs_vecs[j].matches[j_match]);
-
-        /* Try to expand this rotation matrix to a larger network */
-        
-        int success = grow_network(rot_mat, obs_vecs, obs_vec_count,
-                                       i, j, i_match, j_match, max_members);
-
-        /* return this matrix and if it was immediately successful */
-        *rotation = rot_mat;
-        
-        return success;
+	gsl_matrix *rot_mat;
+
+	rot_mat = generate_rot_mat(obs_vecs[i].obsvec,
+				   obs_vecs[j].obsvec,
+				   obs_vecs[i].matches[i_match],
+				   obs_vecs[j].matches[j_match]);
+
+	/* Try to expand this rotation matrix to a larger network */
+
+	int success = grow_network(rot_mat, obs_vecs, obs_vec_count,
+				   i, j, i_match, j_match, max_members);
+
+	/* return this matrix and if it was immediately successful */
+	*rotation = rot_mat;
+
+	return success;
 }
 
 static int find_seed(struct SpotVec *obs_vecs, int obs_vec_count,
-                     gsl_matrix **rotation)
+		     gsl_matrix **rotation)
 {
-    /* META: Maximum achieved maximum network membership */
+	/* META: Maximum achieved maximum network membership */
 	int max_max_members = 0;
 	gsl_matrix *best_rotation = NULL;
-	
+
 	unsigned long start_time = time(NULL);
-	
+
 	/* loop round pairs of vectors to try and find a suitable
 	 * seed to start building a self-consistent network of vectors
 	 */
 	int i, j;
-	
-	for ( i=0; i<obs_vec_count-1; i++ ) {
-	for ( j=i+1; j<obs_vec_count; j++ ) {
-
-		/** Check to see if there is a shared spot - opportunity
-		  * for optimisation by generating a look-up table
-		  * by spot instead of by vector.
-		  */
-		int shared = obs_vecs_share_spot(&obs_vecs[i], &obs_vecs[j]);
-		
-		if ( !shared ) continue;
 
-		/* cell vector "matches" index for i, j respectively */
-		int *i_idx = 0;
-		int *j_idx = 0;
-		int matches;
-
-		/* Check to see if any angles match from the cell vectors */
-		obs_vecs_match_angles(&obs_vecs[i], &obs_vecs[j],
-		                              &i_idx, &j_idx, &matches);
-		if ( matches == 0 ) {
-		        free(i_idx);
-		        free(j_idx);
-		        continue;
-		}
-		
-		/* We have seeds! Pass each of them through the seed-starter  */
-		/* If a seed has the highest achieved membership, make note...*/
-		int k;
-		for ( k=0; k<matches && k < 10; k++ ) {
-		        int max_members = 0;
-
-		        int success = start_seed(obs_vecs, obs_vec_count, i, j,
-		                                 i_idx[k], j_idx[k],
-		                                 rotation, &max_members);
-		        
-		        if (success) {
-		                free(i_idx); free(j_idx);
-		                gsl_matrix_free(best_rotation);
-		                return success;
-		        } else {
-		            if (max_members > max_max_members) {
-		                max_max_members = max_members;
-		                gsl_matrix_free(best_rotation);
-		                best_rotation = *rotation;
-		                *rotation = NULL;
-		            } else {
-		                gsl_matrix_free(*rotation);
-		                *rotation = NULL;
-		            }
-		        }
-		        
-		        unsigned long now_time = time(NULL);
-		        unsigned int seconds = now_time - start_time;
-		        
-		        if (seconds > 30) {
-		                /* Heading towards CrystFEL cutoff so
-		                   return your best guess and run */
-		                free(i_idx); free(j_idx);
-		                
-		                *rotation = best_rotation;
-		                return (best_rotation != NULL);
-		        }
+	for ( i=0; i<obs_vec_count; i++ ) {
+		for ( j=0; j<i; j++ ) {
+
+			/** Check to see if there is a shared spot - opportunity
+			 * for optimisation by generating a look-up table
+			 * by spot instead of by vector.
+			 */
+			int shared = obs_vecs_share_spot(&obs_vecs[i], &obs_vecs[j]);
+
+			if ( !shared ) continue;
+
+			/* cell vector "matches" index for i, j respectively */
+			int *i_idx = 0;
+			int *j_idx = 0;
+			int matches;
+
+			/* Check to see if any angles match from the cell vectors */
+			obs_vecs_match_angles(&obs_vecs[i], &obs_vecs[j],
+					      &i_idx, &j_idx, &matches);
+			if ( matches == 0 ) {
+				free(i_idx);
+				free(j_idx);
+				continue;
+			}
+
+			/* We have seeds! Pass each of them through the seed-starter  */
+			/* If a seed has the highest achieved membership, make note...*/
+			int k;
+			for ( k=0; k<matches; k++ ) {
+				int max_members = 0;
+
+				int success = start_seed(obs_vecs, obs_vec_count, i, j,
+							 i_idx[k], j_idx[k],
+							 rotation, &max_members);
+
+				if (success) {
+					free(i_idx); free(j_idx);
+					gsl_matrix_free(best_rotation);
+					return success;
+				} else {
+					if (max_members > max_max_members) {
+						max_max_members = max_members;
+						gsl_matrix_free(best_rotation);
+						best_rotation = *rotation;
+						*rotation = NULL;
+					} else {
+						gsl_matrix_free(*rotation);
+						*rotation = NULL;
+					}
+				}
+
+				unsigned long now_time = time(NULL);
+				unsigned int seconds = now_time - start_time;
+
+				// Commented out for the time being.
+				/*
+				 if (seconds > 30) {
+				 // Heading towards CrystFEL cutoff so
+				 // return your best guess and run
+					free(i_idx); free(j_idx);
+
+				 *rotation = best_rotation;
+				 return (best_rotation != NULL);
+				 }
+				 */
+			}
+
+			free(i_idx);
+			free(j_idx);
 		}
-		
-                free(i_idx);
-                free(j_idx);
-	}
 	} /* yes this } is meant to be here */
 
 	*rotation = best_rotation;
@@ -882,7 +887,7 @@ static int sort_func(const void *av, const void *bv)
 
 
 static int match_obs_to_cell_vecs(struct rvec *cell_vecs, int cell_vec_count,
-                                  struct SpotVec *obs_vecs, int obs_vec_count)
+				  struct SpotVec *obs_vecs, int obs_vec_count)
 {
 	int i, j;
 
@@ -912,8 +917,10 @@ static int match_obs_to_cell_vecs(struct rvec *cell_vecs, int cell_vec_count,
 			for_sort[count].v = cell_vecs[j];
 			for_sort[count].dist = dist_diff;
 			count++;
-		
+
 		}
+
+		/* Sort in order to get most agreeable matches first */
 		qsort(for_sort, count, sizeof(struct sortme), sort_func);
 		obs_vecs[i].matches = malloc(count*sizeof(struct rvec));
 		obs_vecs[i].match_num = count;
@@ -927,9 +934,15 @@ static int match_obs_to_cell_vecs(struct rvec *cell_vecs, int cell_vec_count,
 	return 1;
 }
 
+static int compare_spot_vecs(const void *av, const void *bv)
+{
+	struct SpotVec *a = (struct SpotVec *)av;
+	struct SpotVec *b = (struct SpotVec *)bv;
+	return a->distance > b->distance;
+}
 
 static int gen_observed_vecs(struct rvec *rlps, int rlp_count,
-                             struct SpotVec **obs_vecs, int *obs_vec_count)
+			     struct SpotVec **obs_vecs, int *obs_vec_count)
 {
 	int i, j;
 	int count = 0;
@@ -940,7 +953,7 @@ static int gen_observed_vecs(struct rvec *rlps, int rlp_count,
 	for ( i=0; i<rlp_count-1; i++ ) {
 		for ( j=i+1; j<rlp_count; j++ ) {
 
-			
+
 			/* calculate difference vector between rlps */
 			struct rvec diff = diff_vec(rlps[i], rlps[j]);
 
@@ -948,12 +961,12 @@ static int gen_observed_vecs(struct rvec *rlps, int rlp_count,
 			double sqlength = sq_length(diff);
 
 			if ( sqlength > max_sq_length ) continue;
-            
+
 			count++;
 
 			struct SpotVec *temp_obs_vecs;
 			temp_obs_vecs = realloc(*obs_vecs,
-			                        count*sizeof(struct SpotVec));
+						count*sizeof(struct SpotVec));
 
 			if ( temp_obs_vecs == NULL ) {
 				return 0;
@@ -975,6 +988,11 @@ static int gen_observed_vecs(struct rvec *rlps, int rlp_count,
 		}
 	}
 
+	/* Sort such that the shortest and least error-prone distances
+	   are searched first.
+	 */
+	qsort(*obs_vecs, count, sizeof(struct SpotVec), compare_spot_vecs);
+
 	*obs_vec_count = count;
 
 	return 1;
@@ -982,7 +1000,7 @@ static int gen_observed_vecs(struct rvec *rlps, int rlp_count,
 
 
 static int gen_theoretical_vecs(UnitCell *cell, struct rvec **cell_vecs,
-                                int *vec_count)
+				int *vec_count)
 {
 	double a, b, c, alpha, beta, gamma;
 	int h_max, k_max, l_max;
@@ -1001,37 +1019,37 @@ static int gen_theoretical_vecs(UnitCell *cell, struct rvec **cell_vecs,
 	int count = 0;
 
 	cell_get_reciprocal(cell, &asx, &asy, &asz,
-	                          &bsx, &bsy, &bsz,
-	                          &csx, &csy, &csz);
+			    &bsx, &bsy, &bsz,
+			    &csx, &csy, &csz);
 
 	for ( h=-h_max; h<=+h_max; h++ ) {
-	for ( k=-k_max; k<=+k_max; k++ ) {
-	for ( l=-l_max; l<=+l_max; l++ ) {
+		for ( k=-k_max; k<=+k_max; k++ ) {
+			for ( l=-l_max; l<=+l_max; l++ ) {
 
-		struct rvec cell_vec;
+				struct rvec cell_vec;
 
-		/* Exclude systematic absences from centering concerns */
-		if ( forbidden_reflection(cell, h, k, l) ) continue;
+				/* Exclude systematic absences from centering concerns */
+				if ( forbidden_reflection(cell, h, k, l) ) continue;
 
-		cell_vec.u = h*asx + k*bsx + l*csx;
-		cell_vec.v = h*asy + k*bsy + l*csy;
-		cell_vec.w = h*asz + k*bsz + l*csz;
+				cell_vec.u = h*asx + k*bsx + l*csx;
+				cell_vec.v = h*asy + k*bsy + l*csy;
+				cell_vec.w = h*asz + k*bsz + l*csz;
 
-		/* add this to our array - which may require expanding */
-		count++;
+				/* add this to our array - which may require expanding */
+				count++;
 
-		struct rvec *temp_cell_vecs;
-		temp_cell_vecs = realloc(*cell_vecs, count*sizeof(struct rvec));
+				struct rvec *temp_cell_vecs;
+				temp_cell_vecs = realloc(*cell_vecs, count*sizeof(struct rvec));
 
-		if ( temp_cell_vecs == NULL ) {
-			return 0;
-		} else {
-			*cell_vecs = temp_cell_vecs;
-			(*cell_vecs)[count - 1] = cell_vec;
+				if ( temp_cell_vecs == NULL ) {
+					return 0;
+				} else {
+					*cell_vecs = temp_cell_vecs;
+					(*cell_vecs)[count - 1] = cell_vec;
+				}
+			}
 		}
 	}
-	}
-	}
 
 	*vec_count = count;
 
@@ -1050,7 +1068,7 @@ static void cleanup_taketwo_cell_vecs(struct rvec *cell_vecs)
 
 
 static void cleanup_taketwo_obs_vecs(struct SpotVec *obs_vecs,
-                                    int obs_vec_count)
+				     int obs_vec_count)
 {
 	int i;
 	for ( i=0; i<obs_vec_count; i++ ) {
@@ -1088,9 +1106,9 @@ static UnitCell *run_taketwo(UnitCell *cell, struct rvec *rlps, int rlp_count)
 
 	success = gen_observed_vecs(rlps, rlp_count, &obs_vecs, &obs_vec_count);
 	if ( !success ) return NULL;
-	
+
 	success = match_obs_to_cell_vecs(cell_vecs, cell_vec_count,
-	                                 obs_vecs, obs_vec_count);
+					 obs_vecs, obs_vec_count);
 
 	cleanup_taketwo_cell_vecs(cell_vecs);
 
@@ -1099,8 +1117,8 @@ static UnitCell *run_taketwo(UnitCell *cell, struct rvec *rlps, int rlp_count)
 	int threshold_reached = find_seed(obs_vecs, obs_vec_count, &solution);
 
 	if ( solution == NULL ) {
-                cleanup_taketwo_obs_vecs(obs_vecs, obs_vec_count);
-	        return NULL;
+		cleanup_taketwo_obs_vecs(obs_vecs, obs_vec_count);
+		return NULL;
 	}
 
 	result = transform_cell_gsl(cell, solution);
@@ -1122,21 +1140,21 @@ int taketwo_index(struct image *image, IndexingPrivate *ipriv)
 	int i;
 	struct taketwo_private *tp = (struct taketwo_private *)ipriv;
 
-        rlps = malloc((image_feature_count(image->features)+1)*sizeof(struct rvec));
-        for ( i=0; i<image_feature_count(image->features); i++ ) {
-            struct imagefeature *pk = image_get_feature(image->features, i);
-            if ( pk == NULL ) continue;
-            rlps[n_rlps].u = pk->rx;
-            rlps[n_rlps].v = pk->ry;
-            rlps[n_rlps].w = pk->rz;
-            n_rlps++;
-        }
-        rlps[n_rlps].u = 0.0;
-        rlps[n_rlps].v = 0.0;
-        rlps[n_rlps++].w = 0.0;
+	rlps = malloc((image_feature_count(image->features)+1)*sizeof(struct rvec));
+	for ( i=0; i<image_feature_count(image->features); i++ ) {
+		struct imagefeature *pk = image_get_feature(image->features, i);
+		if ( pk == NULL ) continue;
+		rlps[n_rlps].u = pk->rx;
+		rlps[n_rlps].v = pk->ry;
+		rlps[n_rlps].w = pk->rz;
+		n_rlps++;
+	}
+	rlps[n_rlps].u = 0.0;
+	rlps[n_rlps].v = 0.0;
+	rlps[n_rlps++].w = 0.0;
 
 	cell = run_taketwo(tp->cell, rlps, n_rlps);
-        free(rlps);
+	free(rlps);
 	if ( cell == NULL ) return 0;
 
 	cr = crystal_new();
@@ -1151,7 +1169,7 @@ int taketwo_index(struct image *image, IndexingPrivate *ipriv)
 		if ( !peak_sanity_check(image, &cr, 1) ) {
 			cell_free(cell);
 			crystal_free(cr);
-			
+
 			return 0;
 		}
 	}
@@ -1163,17 +1181,17 @@ int taketwo_index(struct image *image, IndexingPrivate *ipriv)
 
 
 IndexingPrivate *taketwo_prepare(IndexingMethod *indm, UnitCell *cell,
-                                 struct detector *det, float *ltl)
+				 struct detector *det, float *ltl)
 {
 	struct taketwo_private *tp;
 
 	/* Flags that TakeTwo knows about */
 	*indm &= INDEXING_METHOD_MASK | INDEXING_CHECK_PEAKS
-	       | INDEXING_USE_LATTICE_TYPE | INDEXING_USE_CELL_PARAMETERS
-	       | INDEXING_CONTROL_FLAGS;
+	| INDEXING_USE_LATTICE_TYPE | INDEXING_USE_CELL_PARAMETERS
+	| INDEXING_CONTROL_FLAGS;
 
 	if ( !( (*indm & INDEXING_USE_LATTICE_TYPE)
-	     && (*indm & INDEXING_USE_CELL_PARAMETERS)) )
+	       && (*indm & INDEXING_USE_CELL_PARAMETERS)) )
 	{
 		ERROR("TakeTwo indexing requires cell and lattice type "
 		      "information.\n");
@@ -1184,19 +1202,19 @@ IndexingPrivate *taketwo_prepare(IndexingMethod *indm, UnitCell *cell,
 		ERROR("TakeTwo indexing requires a unit cell.\n");
 		return NULL;
 	}
-	
-	STATUS("Welcome to TakeTwo\n");
-        STATUS("If you use these indexing results, please cite:\n");
-        STATUS("Ginn et al., Acta Cryst. (2016). D72, 956-965\n");
-
 
+	STATUS("Welcome to TakeTwo\n");
+	STATUS("If you use these indexing results, please cite:\n");
+	STATUS("Ginn et al., Acta Cryst. (2016). D72, 956-965\n");
+	
+	
 	tp = malloc(sizeof(struct taketwo_private));
 	if ( tp == NULL ) return NULL;
-
+	
 	tp->ltl = ltl;
 	tp->cell = cell;
 	tp->indm = *indm;
-
+	
 	return (IndexingPrivate *)tp;
 }
 
@@ -1206,3 +1224,4 @@ void taketwo_cleanup(IndexingPrivate *pp)
 	struct taketwo_private *tp = (struct taketwo_private *)pp;
 	free(tp);
 }
+