From 198bdddf99f1db5510dff63b7dc21d01cb98994a Mon Sep 17 00:00:00 2001 From: Thomas White Date: Tue, 25 Oct 2016 17:45:59 +0200 Subject: TEMPORARY horrible debugging stuff --- libcrystfel/src/cell-utils.c | 26 ++++---- 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; ifeatures)+1)*sizeof(struct rvec)); - for ( i=0; ifeatures); 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; -- cgit v1.2.3