diff options
Diffstat (limited to 'libcrystfel/src')
-rw-r--r-- | libcrystfel/src/taketwo.c | 815 |
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); } + |