Crashy for unknown reason, but have added the code for using custom params.

1 files changed, 41 insertions, 23 deletions

diff --git a/libcrystfel/src/taketwo.c b/libcrystfel/src/taketwo.c
index 307c6737..12b608f8 100644
--- a/libcrystfel/src/taketwo.c
+++ b/libcrystfel/src/taketwo.c
@@ -82,6 +82,7 @@ struct TakeTwoCell
 	unsigned int   numOps;
 	struct SpotVec **obs_vecs; // Pointer to an array
 	int obs_vec_count;
+	struct taketwo_options *options;
 	gsl_matrix *twiz1Tmp;
 	gsl_matrix *twiz2Tmp;
 	gsl_vector *vec1Tmp;
@@ -403,7 +404,7 @@ static SymOpList *sym_ops_for_cell(UnitCell *cell)
 
 static int rot_mats_are_similar(gsl_matrix *rot1, gsl_matrix *rot2,
   				gsl_matrix *sub, gsl_matrix *mul,
-				double *score)
+				double *score, struct TakeTwoCell *cell)
 {
 	double tr;
 
@@ -415,9 +416,7 @@ static int rot_mats_are_similar(gsl_matrix *rot1, gsl_matrix *rot2,
 	tr = matrix_trace(mul);
 	if (score != NULL) *score = tr;
 
-	double max = sqrt(4.0*(1.0-cos(TRACE_TOLERANCE)));
-
-	return (tr < max);
+	return (tr < cell->options->trace_tol);
 }
 
 static int symm_rot_mats_are_similar(gsl_matrix *rot1, gsl_matrix *rot2,
@@ -435,7 +434,7 @@ static int symm_rot_mats_are_similar(gsl_matrix *rot1, gsl_matrix *rot2,
 		gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, rot1, symOp,
 		0.0, testRot);
 		
-		if (rot_mats_are_similar(testRot, rot2, sub, mul, NULL)) {
+		if (rot_mats_are_similar(testRot, rot2, sub, mul, NULL, cell)) {
 			gsl_matrix_free(testRot);
 			gsl_matrix_free(sub);
 			gsl_matrix_free(mul);
@@ -559,7 +558,7 @@ static int obs_shares_spot_w_array(struct SpotVec *obs_vecs, int test_idx,
 static int obs_vecs_match_angles(struct SpotVec *her_obs,
                                  struct SpotVec *his_obs,
                                  int **her_match_idxs, int **his_match_idxs,
-                                 int *match_count)
+                                 int *match_count, struct TakeTwoCell *cell)
 {
         int i, j;
         *match_count = 0;
@@ -585,7 +584,7 @@ static int obs_vecs_match_angles(struct SpotVec *her_obs,
 
                 double angle_diff = fabs(theory_angle - obs_angle);
 
-                if ( angle_diff < ANGLE_TOLERANCE ) {
+                if ( angle_diff < cell->options->angle_tol ) {
                         // in the case of a brief check only            
                         if (!her_match_idxs || !his_match_idxs) {
                                 return 1;
@@ -694,7 +693,7 @@ static signed int finish_solution(gsl_matrix *rot, struct SpotVec *obs_vecs,
 			double addition;
 			rot_mats_are_similar(rotations[i], rotations[j],
 					     sub, mul,
-					     &addition);
+					     &addition, cell);
 
 			current_score += addition;
 		}
@@ -839,7 +838,7 @@ static signed int find_next_index(gsl_matrix *rot, int *obs_members,
 
 				double trace = 0;
 				int ok = rot_mats_are_similar(rot, test_rot,
-						       sub, mul, &trace);
+						       sub, mul, &trace, cell);
 				
 				gsl_matrix_free(test_rot);
 
@@ -954,13 +953,13 @@ static int grow_network(gsl_matrix *rot, int obs_idx1, int obs_idx2,
 		}
 
 		/* If member_num is high enough, we want to return a yes */
-		if ( member_num > NETWORK_MEMBER_THRESHOLD ) break;
+		if ( member_num > cell->options->member_thresh ) break;
 	}
 
 	finish_solution(rot, obs_vecs, obs_members,
 			  match_members, member_num, cell);
 	
-	return ( member_num > NETWORK_MEMBER_THRESHOLD );
+	return ( member_num > cell->options->member_thresh );
 }
 
 
@@ -1016,19 +1015,16 @@ static int find_seed(gsl_matrix **rotation, struct TakeTwoCell *cell)
 
 			if ( !shared ) continue;
 
-			/* cell vector index matches stored in i, j and total number
-			 * stored in int matches.
+			/* cell vector index matches stored in i, j and total
+			 * number stored in int matches.
 			 */
 			int *i_idx = 0;
 			int *j_idx = 0;
 			int matches;
 				
-		//	double dist1Pc = 100 * obs_vecs[i].distance / MAX_RECIP_DISTANCE;
-//			double dist2Pc = 100 * obs_vecs[j].distance / MAX_RECIP_DISTANCE;
-			
 			/* 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);
+					      &i_idx, &j_idx, &matches, cell);
 			if ( matches == 0 ) {
 				free(i_idx);
 				free(j_idx);
@@ -1192,7 +1188,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,
-				  int is_asymmetric)
+				  int is_asymmetric, struct TakeTwoCell *cell)
 {
 	int i, j;
 
@@ -1209,7 +1205,7 @@ static int match_obs_to_cell_vecs(struct rvec *cell_vecs, int cell_vec_count,
 			/* check if this matches the observed length */
 			double dist_diff = fabs(cell_length - obs_length);
 
-			if ( dist_diff > RECIP_TOLERANCE ) continue;
+			if ( dist_diff > cell->options->len_tol ) continue;
 
 			/* we have a match, add to array! */
 
@@ -1449,7 +1445,8 @@ static void cleanup_taketwo_cell(struct TakeTwoCell *ttCell)
  * @rot: pointer to be given an assignment (hopefully) if indexing is
  * successful.
  **/
-static UnitCell *run_taketwo(UnitCell *cell, struct rvec *rlps, int rlp_count)
+static UnitCell *run_taketwo(UnitCell *cell, struct taketwo_options *opts,
+			     struct rvec *rlps, int rlp_count)
 {
 	int cell_vec_count = 0;
 	int asym_vec_count = 0;
@@ -1463,6 +1460,7 @@ static UnitCell *run_taketwo(UnitCell *cell, struct rvec *rlps, int rlp_count)
 	
 	struct TakeTwoCell ttCell;
 	ttCell.cell = cell;
+	ttCell.options = opts;
 	ttCell.rotSymOps = NULL;
 	ttCell.twiz1Tmp = gsl_matrix_calloc(3, 3);
 	ttCell.twiz2Tmp = gsl_matrix_calloc(3, 3);
@@ -1483,10 +1481,10 @@ static UnitCell *run_taketwo(UnitCell *cell, struct rvec *rlps, int rlp_count)
 	ttCell.obs_vec_count = obs_vec_count;
 
 	success = match_obs_to_cell_vecs(asym_vecs, asym_vec_count,
-					 obs_vecs, obs_vec_count, 1);
+					 obs_vecs, obs_vec_count, 1, &ttCell);
 
 	success = match_obs_to_cell_vecs(cell_vecs, cell_vec_count,
-					 obs_vecs, obs_vec_count, 0);
+					 obs_vecs, obs_vec_count, 0, &ttCell);
 
 	free(cell_vecs);
 	free(asym_vecs);
@@ -1520,6 +1518,26 @@ int taketwo_index(struct image *image, struct taketwo_options *opts, void *priv)
 	int i;
 	struct taketwo_private *tp = (struct taketwo_private *)priv;
 
+	/* Replace negative values in options with defaults */
+
+	if (opts->member_thresh < 0) {
+		opts->member_thresh = NETWORK_MEMBER_THRESHOLD;
+	}
+
+	if (opts->len_tol < 0) {
+		opts->len_tol = RECIP_TOLERANCE;
+	}
+	
+	if (opts->angle_tol < 0) {
+		opts->angle_tol = ANGLE_TOLERANCE;
+	}
+
+	if (opts->trace_tol < 0) {
+		opts->trace_tol = TRACE_TOLERANCE;
+	}
+	
+	opts->trace_tol = sqrt(4.0*(1.0-cos(opts->trace_tol)));
+
 	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);
@@ -1533,7 +1551,7 @@ int taketwo_index(struct image *image, struct taketwo_options *opts, void *priv)
 	rlps[n_rlps].v = 0.0;
 	rlps[n_rlps++].w = 0.0;
 
-	cell = run_taketwo(tp->cell, rlps, n_rlps);
+	cell = run_taketwo(tp->cell, opts, rlps, n_rlps);
 	free(rlps);
 	if ( cell == NULL ) return 0;