diff options
| -rw-r--r-- | libcrystfel/src/taketwo.c | 74 |
1 files changed, 47 insertions, 27 deletions
diff --git a/libcrystfel/src/taketwo.c b/libcrystfel/src/taketwo.c index cc3d83f0..c9314395 100644 --- a/libcrystfel/src/taketwo.c +++ b/libcrystfel/src/taketwo.c @@ -108,7 +108,7 @@ struct TakeTwoCell #define MAX_NETWORK_MEMBERS (NETWORK_MEMBER_THRESHOLD + 3) /* Maximum dead ends for a single branch extension during indexing */ -#define MAX_DEAD_ENDS (0) +#define MAX_DEAD_ENDS (1) /* Tolerance for two angles to be considered the same */ #define ANGLE_TOLERANCE (deg2rad(0.5)) @@ -227,12 +227,11 @@ static void show_rvec(struct rvec r2) * functions called under the core functions, still specialised (Level 3) * ------------------------------------------------------------------------*/ -static gsl_matrix *rotation_around_axis(struct rvec c, double th) +static gsl_matrix *rotation_around_axis(struct rvec c, double th, + gsl_matrix *res) { double omc = 1.0 - cos(th); double s = sin(th); - gsl_matrix *res = gsl_matrix_alloc(3, 3); - gsl_matrix_set(res, 0, 0, cos(th) + c.u*c.u*omc); gsl_matrix_set(res, 0, 1, c.u*c.v*omc - c.w*s); gsl_matrix_set(res, 0, 2, c.u*c.w*omc + c.v*s); @@ -253,7 +252,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, gsl_matrix *twizzle) { /* Let's have unit vectors */ normalise_rvec(&vec1); @@ -303,7 +302,7 @@ static gsl_matrix *closest_rot_mat(struct rvec vec1, struct rvec vec2, /* Return an identity matrix which has been rotated by * theta around "axis" */ - return rotation_around_axis(axis, bestAngle); + rotation_around_axis(axis, bestAngle, twizzle); } @@ -461,12 +460,13 @@ static int symm_rot_mats_are_similar(gsl_matrix *rot1, gsl_matrix *rot2, return 0; } -static gsl_matrix *rotation_between_vectors(struct rvec a, struct rvec b) +static gsl_matrix *rotation_between_vectors(struct rvec a, struct rvec b, + gsl_matrix *twizzle) { double th = rvec_angle(a, b); struct rvec c = rvec_cross(a, b); normalise_rvec(&c); - return rotation_around_axis(c, th); + rotation_around_axis(c, th, twizzle); } @@ -495,10 +495,9 @@ 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 *twiz1, gsl_matrix *twiz2) { - gsl_matrix *rotateSpotDiffMatrix; - gsl_matrix *secondTwizzleMatrix; gsl_matrix *fullMat; gsl_vector *cell2v = rvec_to_gsl(cell2); gsl_vector *cell2vr = gsl_vector_calloc(3); @@ -510,30 +509,28 @@ static gsl_matrix *generate_rot_mat(struct rvec obs1, struct rvec obs2, /* Rotate reciprocal space so that the first simulated vector lines up * with the observed vector. */ - rotateSpotDiffMatrix = rotation_between_vectors(cell1, obs1); + rotation_between_vectors(cell1, obs1, twiz1); normalise_rvec(&obs1); /* Multiply cell2 by rotateSpotDiffMatrix --> cell2vr */ - gsl_blas_dgemv(CblasNoTrans, 1.0, rotateSpotDiffMatrix, cell2v, + gsl_blas_dgemv(CblasNoTrans, 1.0, twiz1, cell2v, 0.0, cell2vr); /* 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); + closest_rot_mat(gsl_to_rvec(cell2vr), obs2, obs1, twiz2); /* 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); + twiz1, twiz2, 0.0, fullMat); gsl_matrix_transpose(fullMat); gsl_vector_free(cell2v); gsl_vector_free(cell2vr); - gsl_matrix_free(secondTwizzleMatrix); - gsl_matrix_free(rotateSpotDiffMatrix); return fullMat; } @@ -711,11 +708,12 @@ static signed int finish_solution(gsl_matrix *rot, struct SpotVec *obs_vecs, int *obs_members, int *match_members, int member_num) { - gsl_matrix *sub; - gsl_matrix *mul; - sub = gsl_matrix_calloc(3, 3); - mul = gsl_matrix_calloc(3, 3); + gsl_matrix *sub = gsl_matrix_calloc(3, 3); + gsl_matrix *mul = gsl_matrix_calloc(3, 3); + gsl_matrix *twiz1 = gsl_matrix_calloc(3, 3); + gsl_matrix *twiz2 = gsl_matrix_calloc(3, 3); + gsl_matrix **rotations = malloc(sizeof(*rotations)* pow(member_num, 2) - member_num); int i, j, count; @@ -733,7 +731,8 @@ static signed int finish_solution(gsl_matrix *rot, struct SpotVec *obs_vecs, struct rvec j_cellvec = j_vec.matches[match_members[j]]; rotations[count] = generate_rot_mat(i_obsvec, j_obsvec, - i_cellvec, j_cellvec); + i_cellvec, j_cellvec, + twiz1, twiz2); count++; } @@ -768,6 +767,8 @@ static signed int finish_solution(gsl_matrix *rot, struct SpotVec *obs_vecs, free(rotations); gsl_matrix_free(sub); gsl_matrix_free(mul); + gsl_matrix_free(twiz1); + gsl_matrix_free(twiz2); return 1; } @@ -783,7 +784,9 @@ static signed int find_next_index(gsl_matrix *rot, struct SpotVec *obs_vecs, gsl_matrix *mul; sub = gsl_matrix_calloc(3, 3); mul = gsl_matrix_calloc(3, 3); - + gsl_matrix *twiz1 = gsl_matrix_calloc(3, 3); + gsl_matrix *twiz2 = gsl_matrix_calloc(3, 3); + for ( i=start; i<obs_vec_count && i < start + 1000; i++ ) { @@ -836,7 +839,8 @@ static signed int find_next_index(gsl_matrix *rot, struct SpotVec *obs_vecs, test_rot = generate_rot_mat(obs_vecs[idx_k].obsvec, obs_vecs[i].obsvec, member_match, - a_match); + a_match, + twiz1, twiz2); double trace = 0; int ok = rot_mats_are_similar(rot, test_rot, @@ -857,6 +861,8 @@ static signed int find_next_index(gsl_matrix *rot, struct SpotVec *obs_vecs, free(test_idx); gsl_matrix_free(sub); gsl_matrix_free(mul); + gsl_matrix_free(twiz1); + gsl_matrix_free(twiz2); return i; } @@ -869,7 +875,8 @@ static signed int find_next_index(gsl_matrix *rot, struct SpotVec *obs_vecs, gsl_matrix_free(sub); gsl_matrix_free(mul); - + gsl_matrix_free(twiz1); + gsl_matrix_free(twiz2); /* give up. */ @@ -977,12 +984,16 @@ 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) { + gsl_matrix *twiz1 = gsl_matrix_calloc(3, 3); + gsl_matrix *twiz2 = gsl_matrix_calloc(3, 3); + 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]); + obs_vecs[j].matches[j_match], + twiz1, twiz2); /* Try to expand this rotation matrix to a larger network */ @@ -991,6 +1002,8 @@ static int start_seed(struct SpotVec *obs_vecs, int obs_vec_count, int i, /* return this matrix and if it was immediately successful */ *rotation = rot_mat; + gsl_matrix_free(twiz1); + gsl_matrix_free(twiz2); return success; } @@ -1003,6 +1016,9 @@ static int weed_duplicate_matches(struct SpotVec *her_obs, int num_occupied = 0; gsl_matrix **old_mats = calloc(*match_count, sizeof(gsl_matrix *)); + gsl_matrix *twiz1 = gsl_matrix_calloc(3, 3); + gsl_matrix *twiz2 = gsl_matrix_calloc(3, 3); + if (old_mats == NULL) { apologise(); @@ -1022,7 +1038,8 @@ static int weed_duplicate_matches(struct SpotVec *her_obs, struct rvec j_cellvec = his_obs->matches[his_match]; gsl_matrix *mat = generate_rot_mat(i_obsvec, j_obsvec, - i_cellvec, j_cellvec); + i_cellvec, j_cellvec, + twiz1, twiz2); int found = 0; @@ -1056,6 +1073,9 @@ static int weed_duplicate_matches(struct SpotVec *her_obs, free(old_mats); + gsl_matrix_free(twiz1); + gsl_matrix_free(twiz2); + return 1; } |