diff options
Diffstat (limited to 'libcrystfel')
-rw-r--r-- | libcrystfel/src/cell-utils.c | 26 | ||||
-rw-r--r-- | libcrystfel/src/taketwo.c | 137 |
2 files changed, 132 insertions, 31 deletions
diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c index 74a6b905..841d88f9 100644 --- a/libcrystfel/src/cell-utils.c +++ b/libcrystfel/src/cell-utils.c @@ -1445,27 +1445,33 @@ UnitCell *transform_cell_gsl(UnitCell *in, gsl_matrix *m) c = gsl_matrix_alloc(3, 3); gsl_matrix_set(c, 0, 0, asx); - gsl_matrix_set(c, 0, 1, asy); - gsl_matrix_set(c, 0, 2, asz); - gsl_matrix_set(c, 1, 0, bsx); + gsl_matrix_set(c, 1, 0, asy); + gsl_matrix_set(c, 2, 0, asz); + gsl_matrix_set(c, 0, 1, bsx); gsl_matrix_set(c, 1, 1, bsy); - gsl_matrix_set(c, 1, 2, bsz); - gsl_matrix_set(c, 2, 0, csx); - gsl_matrix_set(c, 2, 1, csy); + gsl_matrix_set(c, 2, 1, bsz); + gsl_matrix_set(c, 0, 2, csx); + gsl_matrix_set(c, 1, 2, csy); gsl_matrix_set(c, 2, 2, csz); + STATUS("---\n"); + show_matrix(m); + STATUS("\n"); + show_matrix(c); + STATUS("\n"); res = gsl_matrix_calloc(3, 3); gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, m, c, 0.0, res); + show_matrix(res); out = cell_new_from_cell(in); cell_set_reciprocal(out, gsl_matrix_get(res, 0, 0), - gsl_matrix_get(res, 0, 1), - gsl_matrix_get(res, 0, 2), gsl_matrix_get(res, 1, 0), - gsl_matrix_get(res, 1, 1), - gsl_matrix_get(res, 1, 2), gsl_matrix_get(res, 2, 0), + gsl_matrix_get(res, 0, 1), + gsl_matrix_get(res, 1, 1), gsl_matrix_get(res, 2, 1), + gsl_matrix_get(res, 0, 2), + gsl_matrix_get(res, 1, 2), gsl_matrix_get(res, 2, 2)); gsl_matrix_free(res); diff --git a/libcrystfel/src/taketwo.c b/libcrystfel/src/taketwo.c index 036db5c6..9c173303 100644 --- a/libcrystfel/src/taketwo.c +++ b/libcrystfel/src/taketwo.c @@ -69,6 +69,9 @@ struct taketwo_private UnitCell *cell; }; +struct rvec *global_rlps; +int global_nrlps; + /* Maximum distance between two rlp sizes to consider info for indexing */ #define MAX_RECIP_DISTANCE (0.15*1e10) @@ -77,10 +80,10 @@ struct taketwo_private #define RECIP_TOLERANCE (0.001*1e10) /* Threshold for network members to consider a potential solution */ -#define NETWORK_MEMBER_THRESHOLD (20) +#define NETWORK_MEMBER_THRESHOLD (500) /* Maximum network members (obviously a solution so should stop) */ -#define MAX_NETWORK_MEMBERS (100) +#define MAX_NETWORK_MEMBERS (500) /* Maximum dead ends for a single branch extension during indexing */ #define MAX_DEAD_ENDS (5) @@ -189,6 +192,12 @@ static struct rvec rvec_cross(struct rvec a, struct rvec b) } +static void show_rvec(struct rvec r) +{ + STATUS("[ %+9.3e %+9.3e %+9.3e ]\n", r.u, r.v, r.w); +} + + /* ------------------------------------------------------------------------ * functions called under the core functions, still specialised (Level 3) * ------------------------------------------------------------------------*/ @@ -326,9 +335,18 @@ 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); + STATUS("--\n"); + show_rvec(obs1); + show_rvec(obs2); + show_rvec(cell1); + show_rvec(cell2); + STATUS(".....\n"); + /* Rotate reciprocal space so that the first simulated vector lines up * with the observed vector. */ rotateSpotDiffMatrix = rotation_between_vectors(cell1, obs1); + STATUS("rotateSpotDiffMatrix:\n"); + show_matrix(rotateSpotDiffMatrix); normalise_rvec(&obs1); @@ -338,14 +356,22 @@ static gsl_matrix *generate_rot_mat(struct rvec obs1, struct rvec obs2, /* 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); + STATUS("for secondTwizzleMatrix:\n"); + show_rvec(gsl_to_rvec(cell2vr)); + show_rvec(obs2); + show_rvec(obs1); + STATUS("..\n"); + secondTwizzleMatrix = closest_rot_mat(gsl_to_rvec(cell2vr), obs2, obs1); + STATUS("secondTwizzleMatrix:\n"); + show_matrix(secondTwizzleMatrix); /* 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(CblasNoTrans, CblasNoTrans, 1.0, - secondTwizzleMatrix, rotateSpotDiffMatrix, 0.0, fullMat); + gsl_blas_dgemm(CblasTrans, CblasTrans, 1.0, + rotateSpotDiffMatrix, secondTwizzleMatrix, 0.0, fullMat); + STATUS("final matrix:\n"); + show_matrix(fullMat); return fullMat; } @@ -448,6 +474,8 @@ static int obs_angles_match_array(struct SpotVec *obs_vecs, int test_idx, &idx1, &idx2); if (idx2 != match_members[i]) return 0; + //STATUS("match %i %i %i %i\n", idx1, idx2, test_idx, + // match_members[i]); } return 1; @@ -499,8 +527,7 @@ static signed int find_next_index(gsl_matrix *rot, struct SpotVec *obs_vecs, int j; /* if ok is set to 0, give up on this vector before - * checking the next value of j - */ + * thecking the next value of j */ int ok = 1; for ( j=0; j<2; j++ ) { @@ -510,6 +537,17 @@ static signed int find_next_index(gsl_matrix *rot, struct SpotVec *obs_vecs, int j_idx = obs_members[j]; member_match = obs_vecs[j_idx].matches[match_members[j]]; + STATUS("]]]]] %i %i %i %i\n", + j_idx, i, match_members[j], test_idx); + + //STATUS("test:\n"); + //struct rvec obs1 = {0, 1, 2}; + //struct rvec obs2 = {3, 4, 5}; + //struct rvec cell1 = {5, 7, 2}; + //struct rvec cell2 = {-1, 4, 3}; + //generate_rot_mat(obs1, obs2, cell1, cell2); + //STATUS("end of test\n"); + //exit(0); test_rot = generate_rot_mat(obs_vecs[j_idx].obsvec, obs_vecs[i].obsvec, @@ -517,6 +555,10 @@ static signed int find_next_index(gsl_matrix *rot, struct SpotVec *obs_vecs, test_match); ok = rot_mats_are_similar(rot, test_rot); + //exit(0); + //STATUS("--\n"); + //show_matrix(rot); + STATUS("ok = %i\n", ok); if ( !ok ) break; } @@ -555,9 +597,12 @@ 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; + STATUS("starting..\n"); while ( 1 ) { int match_found = -1; + STATUS("member_num = %i\n", member_num); + signed int next_index = find_next_index(rot, obs_vecs, obs_vec_count, obs_members, match_members, start, member_num, @@ -585,7 +630,7 @@ static int grow_network(gsl_matrix *rot, struct SpotVec *obs_vecs, obs_members[member_num] = next_index; match_members[member_num] = match_found; - + start = next_index + 1; member_num++; @@ -599,6 +644,17 @@ static int grow_network(gsl_matrix *rot, struct SpotVec *obs_vecs, } +static signed int spot_idx(struct rvec *rlp) +{ + int i; + for ( i=0; i<global_nrlps; i++ ) { + if ( rlp == &global_rlps[i] ) return i; + } + return -1; +} + + + static int find_seed_and_network(struct SpotVec *obs_vecs, int obs_vec_count, gsl_matrix **rotation) { @@ -615,6 +671,7 @@ static int find_seed_and_network(struct SpotVec *obs_vecs, int obs_vec_count, * by spot instead of by vector. */ int shared = obs_vecs_share_spot(&obs_vecs[i], &obs_vecs[j]); + //STATUS("checking %i %i\n", i, j); if ( !shared ) continue; @@ -627,6 +684,7 @@ static int find_seed_and_network(struct SpotVec *obs_vecs, int obs_vec_count, match = obs_vecs_match_angles(&obs_vecs[i], &obs_vecs[j], &i_idx, &j_idx); if ( !match ) continue; + STATUS("...ok\n"); /* We have a seed! Generate a matrix based on this solution */ gsl_matrix *rot_mat = gsl_matrix_calloc(3, 3); @@ -636,8 +694,14 @@ static int find_seed_and_network(struct SpotVec *obs_vecs, int obs_vec_count, obs_vecs[i].matches[i_idx], obs_vecs[j].matches[j_idx]); - /* try to expand this rotation matrix to a larger network */ - + /* Try to expand this rotation matrix to a larger network */ + STATUS("idx: %i %i %i %i spots %i %i and %i %i\n", + i_idx, j_idx, i, j, + spot_idx(obs_vecs[i].her_rlp), + spot_idx(obs_vecs[i].his_rlp), + spot_idx(obs_vecs[j].her_rlp), + spot_idx(obs_vecs[j].his_rlp)); + show_matrix(rot_mat); int success = grow_network(rot_mat, obs_vecs, obs_vec_count, i, j, i_idx, j_idx); @@ -680,6 +744,9 @@ static int match_obs_to_cell_vecs(struct rvec *cell_vecs, int cell_vec_count, int count = 0; struct sortme *for_sort = NULL; + if ( (i==0) || (i==1) || (i==2) ) { + STATUS("length of %i = %f\n", i, obs_vecs[i].distance/1e10); + } for ( j=0; j<cell_vec_count; j++ ) { @@ -702,12 +769,21 @@ 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++; - + + } + if ( (i==0) || (i==1) || (i==2) ) { + STATUS("%i matches for %i\n", count, i); } qsort(for_sort, count, sizeof(struct sortme), sort_func); obs_vecs[i].matches = malloc(count*sizeof(struct rvec)); for ( j=0; j<count; j++ ) { obs_vecs[i].matches[j] = for_sort[j].v; + if ( (i==0) || (i==1) || (i==2) ) { + STATUS("%i %14f %14f : ", j, + rvec_length(for_sort[j].v)/1e10, + for_sort[j].dist/1e10); + show_rvec(for_sort[j].v); + } } free(for_sort); @@ -729,6 +805,8 @@ 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++ ) { + //STATUS("%i %i %i %i\n", i, j, rlp_count, count); + /* calculate difference vector between rlps */ struct rvec diff = diff_vec(rlps[i], rlps[j]); @@ -871,11 +949,18 @@ static UnitCell *run_taketwo(UnitCell *cell, struct rvec *rlps, int rlp_count) int success = 0; gsl_matrix *solution = NULL; + STATUS("theoretical..\n"); success = gen_theoretical_vecs(cell, &cell_vecs, &cell_vec_count); if ( !success ) return NULL; + STATUS("cell_vec_count = %i\n", cell_vec_count); + +global_rlps = rlps; +global_nrlps = rlp_count; success = gen_observed_vecs(rlps, rlp_count, &obs_vecs, &obs_vec_count); if ( !success ) return NULL; + STATUS("obs_vec_count = %i\n", obs_vec_count); + STATUS("rlp_count = %i\n", rlp_count); success = match_obs_to_cell_vecs(cell_vecs, cell_vec_count, obs_vecs, obs_vec_count); @@ -903,21 +988,31 @@ int taketwo_index(struct image *image, IndexingPrivate *ipriv) UnitCell *cell; struct rvec *rlps; int n_rlps = 0; - 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++; + FILE *fh = fopen("../../spots.csv", "r"); + rlps = malloc(100*sizeof(struct rvec)); + while ( (fh != NULL) && !feof(fh) ) { + float x, y, z; + float dd; + int di; + char line[1024]; + fgets(line, 1024, fh); + + if ( sscanf(line, "%f,%f,%f,%f,%i,%i,%f,%f,%i", + &x, &y, &z, &dd, &di, &di, &dd, &dd, &di) == 9 ) { + if ( isnan(x) ) continue; + STATUS("got %f %f %f\n", x, y, z); + rlps[n_rlps].u = x*1e10; + rlps[n_rlps].v = y*1e10; + rlps[n_rlps].w = z*1e10; + n_rlps++; + } } rlps[n_rlps].u = 0.0; rlps[n_rlps].v = 0.0; rlps[n_rlps++].w = 0.0; + if ( fh != NULL ) fclose(fh); cell = run_taketwo(tp->cell, rlps, n_rlps); if ( cell == NULL ) return 0; |