aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--libcrystfel/src/reax.c140
1 files changed, 134 insertions, 6 deletions
diff --git a/libcrystfel/src/reax.c b/libcrystfel/src/reax.c
index 4387d039..2217d6d9 100644
--- a/libcrystfel/src/reax.c
+++ b/libcrystfel/src/reax.c
@@ -242,6 +242,125 @@ static double check_dir(struct dvec *dir, ImageFeatureList *flist,
}
+/* Refine a direct space vector. From Clegg (1984) */
+static double iterate_refine_vector(double *x, double *y, double *z,
+ ImageFeatureList *flist)
+{
+ int fi, n, err;
+ gsl_matrix *C;
+ gsl_vector *A;
+ gsl_vector *t;
+ gsl_matrix *s_vec;
+ gsl_vector *s_val;
+ double dtmax;
+
+ A = gsl_vector_calloc(3);
+ C = gsl_matrix_calloc(3, 3);
+
+ n = image_feature_count(flist);
+ fesetround(1);
+ for ( fi=0; fi<n; fi++ ) {
+
+ struct imagefeature *f;
+ double val;
+ double kn, kno;
+ double xv[3];
+ int i, j;
+
+ f = image_get_feature(flist, fi);
+ if ( f == NULL ) continue;
+
+ kno = f->rx*(*x) + f->ry*(*y) + f->rz*(*z); /* Sorry ... */
+ kn = nearbyint(kno);
+ if ( kn - kno > 0.3 ) continue;
+
+ xv[0] = f->rx; xv[1] = f->ry; xv[2] = f->rz;
+
+ for ( i=0; i<3; i++ ) {
+
+ val = gsl_vector_get(A, i);
+ gsl_vector_set(A, i, val+xv[i]*kn);
+
+ for ( j=0; j<3; j++ ) {
+ val = gsl_matrix_get(C, i, j);
+ gsl_matrix_set(C, i, j, val+xv[i]*xv[j]);
+ }
+
+ }
+
+ }
+
+ s_val = gsl_vector_calloc(3);
+ s_vec = gsl_matrix_calloc(3, 3);
+ err = gsl_linalg_SV_decomp_jacobi(C, s_vec, s_val);
+ if ( err ) {
+ ERROR("SVD failed: %s\n", gsl_strerror(err));
+ gsl_matrix_free(s_vec);
+ gsl_vector_free(s_val);
+ gsl_matrix_free(C);
+ gsl_vector_free(A);
+ return 0.0;
+ }
+
+ t = gsl_vector_calloc(3);
+ err = gsl_linalg_SV_solve(C, s_vec, s_val, A, t);
+ if ( err ) {
+ ERROR("Matrix solution failed: %s\n", gsl_strerror(err));
+ gsl_matrix_free(s_vec);
+ gsl_vector_free(s_val);
+ gsl_matrix_free(C);
+ gsl_vector_free(A);
+ gsl_vector_free(t);
+ return 0.0;
+ }
+
+ gsl_matrix_free(s_vec);
+ gsl_vector_free(s_val);
+
+ dtmax = fabs(*x - gsl_vector_get(t, 0));
+ dtmax += fabs(*y - gsl_vector_get(t, 1));
+ dtmax += fabs(*z - gsl_vector_get(t, 2));
+
+ *x = gsl_vector_get(t, 0);
+ *y = gsl_vector_get(t, 1);
+ *z = gsl_vector_get(t, 2);
+
+ gsl_matrix_free(C);
+ gsl_vector_free(A);
+
+ return dtmax;
+}
+
+
+static void refine_cell(struct image *image, UnitCell *cell,
+ ImageFeatureList *flist)
+{
+ double ax, ay, az;
+ double bx, by, bz;
+ double cx, cy, cz;
+ int i;
+ double sm;
+
+ cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
+ i = 0;
+ do {
+
+ sm = iterate_refine_vector(&ax, &ay, &az, flist);
+ sm += iterate_refine_vector(&bx, &by, &bz, flist);
+ sm += iterate_refine_vector(&cx, &cy, &cz, flist);
+ i++;
+
+ } while ( (sm > 0.001e-9) && (i<10) );
+
+ cell_set_cartesian(cell, ax, ay, az, bx, by, bz, cx, cy, cz);
+
+ if ( i == 10 ) {
+ cell_free(image->indexed_cell);
+ image->indexed_cell = NULL;
+ }
+}
+
+
static void fine_search(struct reax_private *p, ImageFeatureList *flist,
double pmax, double *fft_in, fftw_complex *fft_out,
struct reax_search_v *sv, struct reax_candidate *c,
@@ -333,7 +452,7 @@ static void squash_vectors(struct reax_search *s, double tol)
}
}
- new = calloc(sv->n_cand-n_invalid,
+ new = calloc(sv->n_cand - n_invalid,
sizeof(struct reax_candidate));
if ( new == NULL ) {
ERROR("Failed to allocate memory for squashed"
@@ -359,9 +478,10 @@ static void squash_vectors(struct reax_search *s, double tol)
sv->cand = new;
for ( j=0; j<sv->n_cand; j++ ) {
- STATUS("%i: %+6.2f %+6.2f %+6.2f %.2f\n",
- j, sv->cand[j].v.x, sv->cand[j].v.y,
- sv->cand[j].v.z, sv->cand[j].fom);
+ STATUS("%i: %+6.2f %+6.2f %+6.2f nm %.2f\n",
+ j, sv->cand[j].v.x*1e9,
+ sv->cand[j].v.y*1e9,
+ sv->cand[j].v.z*1e9, sv->cand[j].fom);
}
}
@@ -412,6 +532,14 @@ static void find_candidates(struct reax_private *p,
for ( j=0; j<sv->n_cand; j++ ) {
fine_search(p, flist, pmax, fft_in, fft_out, sv,
&sv->cand[j], rg, det);
+ STATUS("%i: %+6.2f %+6.2f %+6.2f, mod %.2f nm, %.2f\n",
+ j, sv->cand[j].v.x*1e9,
+ sv->cand[j].v.y*1e9,
+ sv->cand[j].v.z*1e9,
+ modulus(sv->cand[j].v.x,
+ sv->cand[j].v.y,
+ sv->cand[j].v.z)*1e9,
+ sv->cand[j].fom);
}
}
}
@@ -668,6 +796,8 @@ static void assemble_cells_from_candidates(struct image *image,
struct reax_search *s,
UnitCell *cell)
{
+ /* FIXME: Implement this */
+ image->ncells = 0;
}
@@ -702,8 +832,6 @@ void reax_index(IndexingPrivate *pp, struct image *image, UnitCell *cell)
fftw_free(fft_in);
fftw_free(fft_out);
-
- image->ncells = 1;
}