aboutsummaryrefslogtreecommitdiff
path: root/libcrystfel
diff options
context:
space:
mode:
Diffstat (limited to 'libcrystfel')
-rw-r--r--libcrystfel/Makefile.am4
-rw-r--r--libcrystfel/src/asdf.c1175
-rw-r--r--libcrystfel/src/asdf.h28
-rw-r--r--libcrystfel/src/cell-utils.c28
-rw-r--r--libcrystfel/src/cell-utils.h2
-rw-r--r--libcrystfel/src/index.c29
-rw-r--r--libcrystfel/src/index.h4
7 files changed, 1261 insertions, 9 deletions
diff --git a/libcrystfel/Makefile.am b/libcrystfel/Makefile.am
index 9fbd4483..989ce434 100644
--- a/libcrystfel/Makefile.am
+++ b/libcrystfel/Makefile.am
@@ -10,7 +10,7 @@ libcrystfel_la_SOURCES = src/reflist.c src/utils.c src/cell.c src/detector.c \
src/render.c src/index.c src/dirax.c src/mosflm.c \
src/cell-utils.c src/integer_matrix.c src/crystal.c \
src/grainspotter.c src/xds.c src/integration.c \
- src/histogram.c src/events.c
+ src/histogram.c src/events.c src/asdf.c
if HAVE_FFTW
libcrystfel_la_SOURCES += src/reax.c
@@ -30,7 +30,7 @@ libcrystfel_la_include_HEADERS = ${top_srcdir}/version.h \
src/integer_matrix.h src/crystal.h \
src/grainspotter.h src/xds.h \
src/integration.h src/histogram.h \
- src/events.h
+ src/events.h src/asdf.h
AM_CPPFLAGS = -DDATADIR=\""$(datadir)"\" -I$(top_builddir)/lib -Wall
AM_CPPFLAGS += -I$(top_srcdir)/lib @LIBCRYSTFEL_CFLAGS@
diff --git a/libcrystfel/src/asdf.c b/libcrystfel/src/asdf.c
new file mode 100644
index 00000000..d28fcf6a
--- /dev/null
+++ b/libcrystfel/src/asdf.c
@@ -0,0 +1,1175 @@
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+#include <string.h>
+#include <gsl/gsl_matrix.h>
+#include <gsl/gsl_vector.h>
+#include <gsl/gsl_blas.h>
+#include <gsl/gsl_linalg.h>
+#include <gsl/gsl_multifit.h>
+#include <fftw3.h>
+
+#include "image.h"
+#include "dirax.h"
+#include "utils.h"
+#include "peaks.h"
+#include "cell-utils.h"
+#include "asdf.h"
+
+
+struct fftw_vars {
+ int N;
+ fftw_plan p;
+ double *in;
+ fftw_complex *out;
+};
+
+
+struct asdf_private {
+ IndexingMethod indm;
+ float *ltl;
+ UnitCell *template;
+ struct fftw_vars fftw;
+};
+
+
+/* Possible direct vector */
+struct tvector {
+ gsl_vector *t;
+ int n; // number of fitting reflections
+ int *fits;
+};
+
+
+struct fftw_vars fftw_vars_new()
+{
+ int N = 1024;
+
+ struct fftw_vars fftw;
+
+ fftw.N = N;
+ fftw.in = fftw_malloc(N * sizeof(double));
+ fftw.out = fftw_malloc(N * sizeof(fftw_complex));
+ fftw.p = fftw_plan_dft_r2c_1d(N, fftw.in, fftw.out, FFTW_MEASURE);
+
+ return fftw;
+}
+
+
+static int fftw_vars_free(struct fftw_vars fftw)
+{
+ fftw_free(fftw.in);
+ fftw_free(fftw.out);
+ fftw_destroy_plan(fftw.p);
+ fftw_cleanup();
+}
+
+
+struct asdf_cell {
+ gsl_vector *axes[3];
+ gsl_vector *reciprocal[3];
+
+ int n; // number of fitting reflections
+ double volume;
+
+ int N_refls; // total number of reflections
+ int *reflections; // reflections[i] = 1 if reflections fits
+ double **indices; // indices[i] = [h, k, l] for all reflections (not rounded)
+
+ int acl; // minimum number of reflections fitting to one of the axes[]
+ int n_max; // maximum number of reflections fitting to some t-vector
+};
+
+
+struct tvector tvector_new(int n)
+{
+ struct tvector t;
+
+ t.t = gsl_vector_alloc(3);
+ t.n = 0;
+ t.fits = malloc(sizeof(int) * n);
+
+ return t;
+}
+
+
+static int tvector_free(struct tvector t)
+{
+ gsl_vector_free(t.t);
+ free(t.fits);
+
+ return 1;
+}
+
+
+static int tvector_memcpy(struct tvector *dest, struct tvector *src, int n)
+{
+ gsl_vector_memcpy(dest->t, src->t);
+ dest->n = src->n;
+ memcpy(dest->fits, src->fits, sizeof(int) * n);
+
+ return 1;
+}
+
+
+static int asdf_cell_free(struct asdf_cell *c)
+{
+ int i;
+ for ( i = 0; i < 3; i++ ) {
+ gsl_vector_free(c->axes[i]);
+ gsl_vector_free(c->reciprocal[i]);
+ }
+
+ free(c->reflections);
+ for ( i = 0; i < c->N_refls; i++ ) {
+ free(c->indices[i]);
+ }
+ free(c->indices);
+ free(c);
+
+ return 1;
+}
+
+
+static struct asdf_cell *asdf_cell_new(int n)
+{
+ struct asdf_cell *c;
+ c = malloc(sizeof(struct asdf_cell));
+
+ int i;
+ for ( i = 0; i < 3; i++ ) {
+ c->axes[i] = gsl_vector_alloc(3);
+ c->reciprocal[i] = gsl_vector_alloc(3);
+ }
+
+ c->N_refls = n;
+ c->reflections = malloc(sizeof(int) * n);
+ if (c->reflections == NULL) return NULL;
+
+ c->indices = malloc(sizeof(double *) * n);
+ if (c->indices == NULL) return NULL;
+
+ for ( i = 0; i < n; i++ ) {
+ c->indices[i] = malloc(sizeof(double) * 3);
+ if (c->indices[i] == NULL) return NULL;
+ }
+
+ c->n = 0;
+
+ c->acl = 0;
+ c->n_max = 0;
+
+ return c;
+}
+
+
+static int asdf_cell_memcpy(struct asdf_cell *dest, struct asdf_cell *src)
+{
+ int i;
+ for ( i = 0; i < 3; i++ ) {
+ gsl_vector_memcpy(dest->axes[i], src->axes[i]);
+ gsl_vector_memcpy(dest->reciprocal[i], src->reciprocal[i]);
+ }
+
+ dest->volume = src->volume;
+
+ int n = src->N_refls;
+ dest->N_refls = n;
+
+ dest->n = src->n;
+ memcpy(dest->reflections, src->reflections, sizeof(int) * n);
+
+ for (i = 0; i < n; i++ ) {
+ memcpy(dest->indices[i], src->indices[i], sizeof(double) * 3);
+ }
+
+ dest->acl = src->acl;
+ dest->n_max = src->n_max;
+ return 1;
+}
+
+
+/* result = vec1 cross vec2 */
+static int cross_product(gsl_vector *vec1, gsl_vector *vec2,
+ gsl_vector **result)
+{
+ double c1[3], c2[3], p[3];
+ int i;
+ for ( i = 0; i < 3; i++ ) {
+ c1[i] = gsl_vector_get(vec1, i);
+ c2[i] = gsl_vector_get(vec2, i);
+ }
+
+ p[0] = c1[1] * c2[2] - c1[2] * c2[1];
+ p[1] = - c1[0] * c2[2] + c1[2] * c2[0];
+ p[2] = c1[0] * c2[1] - c1[1] * c2[0];
+
+ for ( i = 0; i < 3; i++ ) {
+ gsl_vector_set(*result, i, p[i]);
+ }
+
+ return 1;
+}
+
+
+/* Returns triple product of three gsl_vectors */
+static double calc_volume(gsl_vector *vec1, gsl_vector *vec2, gsl_vector *vec3)
+{
+ double volume;
+ gsl_vector *cross = gsl_vector_alloc(3);
+
+ cross_product(vec1, vec2, &cross);
+ gsl_blas_ddot(vec3, cross, &volume);
+
+ gsl_vector_free(cross);
+ return volume;
+}
+
+
+static int calc_reciprocal(gsl_vector **direct, gsl_vector **reciprocal)
+{
+ double volume;
+
+ cross_product(direct[1], direct[2], &reciprocal[0]);
+ gsl_blas_ddot(direct[0], reciprocal[0], &volume);
+ gsl_vector_scale(reciprocal[0], 1/volume);
+
+ cross_product(direct[2], direct[0], &reciprocal[1]);
+ gsl_vector_scale(reciprocal[1], 1/volume);
+
+ cross_product(direct[0], direct[1], &reciprocal[2]);
+ gsl_vector_scale(reciprocal[2], 1/volume);
+
+ return 1;
+}
+
+
+static int check_cell(struct asdf_private *dp, struct image *image,
+ UnitCell *cell)
+{
+ UnitCell *out;
+ Crystal *cr;
+
+ if ( dp->indm & INDEXING_CHECK_CELL_COMBINATIONS ) {
+
+ out = match_cell(cell, dp->template, 0, dp->ltl, 1);
+ if ( out == NULL ) return 0;
+
+ } else if ( dp->indm & INDEXING_CHECK_CELL_AXES ) {
+
+ out = match_cell(cell, dp->template, 0, dp->ltl, 0);
+ if ( out == NULL ) return 0;
+
+ } else {
+ out = cell_new_from_cell(cell);
+ }
+
+ cr = crystal_new();
+ if ( cr == NULL ) {
+ ERROR("Failed to allocate crystal.\n");
+ return 0;
+ }
+
+ crystal_set_cell(cr, out);
+
+ if ( dp->indm & INDEXING_CHECK_PEAKS ) {
+ if ( !peak_sanity_check(image, &cr, 1) ) {
+ crystal_free(cr); /* Frees the cell as well */
+ cell_free(out);
+ return 0;
+ }
+ }
+
+ image_add_crystal(image, cr);
+
+ return 1;
+}
+
+
+static int compare_doubles(const void *a, const void *b)
+{
+ const double *da = (const double *) a;
+ const double *db = (const double *) b;
+
+ return (*da > *db) - (*da < *db);
+}
+
+
+/* Compares tvectors by length */
+static int compare_tvectors(const void *a, const void *b)
+{
+ struct tvector *ta = (struct tvector *) a;
+ struct tvector *tb = (struct tvector *) b;
+
+ //~ if (ta->n == tb->n) {
+ return (gsl_blas_dnrm2(ta->t) > gsl_blas_dnrm2(tb->t)) -
+ (gsl_blas_dnrm2(ta->t) < gsl_blas_dnrm2(tb->t));
+ //~ }
+ //~
+ //~ return (ta->n > tb->n) - (ta->n < tb->n);
+}
+
+
+/* Calculates normal to a triplet c1, c2, c3. Returns 0 if reflections are on
+ * the same line */
+static int calc_normal(gsl_vector *c1, gsl_vector *c2, gsl_vector *c3,
+ gsl_vector *normal)
+{
+ gsl_vector *c12 = gsl_vector_alloc(3);
+ gsl_vector *c23 = gsl_vector_alloc(3);
+ gsl_vector *c31 = gsl_vector_alloc(3);
+ gsl_vector *res = gsl_vector_alloc(3);
+
+ cross_product(c1, c2, &c12);
+ cross_product(c2, c3, &c23);
+ cross_product(c3, c1, &c31);
+
+ int i;
+ for ( i = 0; i < 3; i++ ) {
+ gsl_vector_set(res, i, gsl_vector_get(c12, i) +
+ gsl_vector_get(c23, i) +
+ gsl_vector_get(c31, i));
+ }
+
+ gsl_vector_free(c12);
+ gsl_vector_free(c23);
+ gsl_vector_free(c31);
+
+ double norm = gsl_blas_dnrm2(res);
+ if ( norm < 0.0001 ) {
+ gsl_vector_free(res);
+ return 0;
+ } else {
+ gsl_vector_scale(res, 1/norm);
+ gsl_vector_memcpy(normal, res);
+ gsl_vector_free(res);
+ }
+
+ return 1;
+}
+
+
+static float find_ds_fft(double *projections, int N_projections, double d_max,
+ struct fftw_vars fftw)
+{
+ int n = N_projections;
+ double projections_sorted[n];
+ memcpy(projections_sorted, projections, sizeof(double) * n);
+ qsort(projections_sorted, n, sizeof(double), compare_doubles);
+
+ int i, k;
+
+ int N = fftw.N; // number of points in fft calculation
+ double *in = fftw.in;
+ fftw_complex *out = fftw.out;
+ fftw_plan p = fftw.p;
+
+ for ( i = 0; i < N; i++ ) {
+ in[i] = 0;
+ }
+
+ for ( i = 0; i < n; i++ ) {
+
+ k = (int)((projections_sorted[i] - projections_sorted[0]) /
+ (projections_sorted[n - 1] - projections_sorted[0]) *
+ (N - 1));
+ in[k] ++;
+ }
+
+ fftw_execute_dft_r2c(p, in, out);
+
+ int i_max = (int)(d_max * (projections_sorted[n - 1] -
+ projections_sorted[0]));
+
+ int d = 1;
+ double max = 0;
+ double a;
+ for ( i = 1; i <= i_max; i++ ) {
+ a = sqrt(out[i][0] * out[i][0] + out[i][1] * out[i][1]);
+ if (a > max) {
+ max = a;
+ d = i;
+ }
+ }
+
+ double ds = (projections_sorted[n - 1] - projections_sorted[0]) / d;
+
+ return ds;
+}
+
+
+/* Returns number of reflections fitting ds.
+ * A projected reflection fits a one-dimensional lattice with elementary
+ * lattice vector d* if its absolute distance to the nearest lattice
+ * point is less than LevelFit. */
+static int check_refl_fitting_ds(double *projections, int N_projections,
+ double ds, double LevelFit)
+{
+ if ( ds == 0 ) return 0;
+
+ int i;
+ int n = 0;
+ for ( i = 0; i < N_projections; i++ ) {
+ if ( fabs(projections[i] -
+ ds * round(projections[i]/ds)) < LevelFit )
+ {
+ n += 1;
+ }
+ }
+
+ return n;
+}
+
+
+/* Refines d*, writes 1 to fits[i] if the i-th projection fits d* */
+static float refine_ds(double *projections, int N_projections, double ds,
+ double LevelFit, int *fits)
+{
+ double fit_refls[N_projections];
+ double indices[N_projections];
+
+ int i;
+
+ int N_fits = 0;
+ int N_new = N_projections;
+
+ double c1, cov11, sumsq;
+ double ds_new = ds;
+ while ( N_fits < N_new ) {
+ N_fits = 0;
+ for ( i = 0; i < N_projections; i++ ) {
+ if ( fabs(projections[i] - ds_new *
+ round(projections[i] / ds_new)) < LevelFit )
+ {
+ fit_refls[N_fits] = projections[i];
+ indices[N_fits] = round(projections[i]/ds_new);
+ N_fits ++;
+ fits[i] = 1;
+ } else {
+ fits[i] = 0;
+ }
+ }
+
+
+ gsl_fit_mul(indices, 1, fit_refls, 1, N_fits, &c1, &cov11,
+ &sumsq);
+ N_new = check_refl_fitting_ds(projections, N_projections, c1,
+ LevelFit);
+ if ( N_new >= N_fits ) {
+ ds_new = c1;
+ }
+ }
+
+ return ds_new;
+}
+
+
+static int check_refl_fitting_cell(struct asdf_cell *c,
+ gsl_vector **reflections,
+ int N_reflections, double IndexFit)
+{
+ double dist[3];
+
+ calc_reciprocal(c->axes, c->reciprocal);
+ c->n = 0;
+ int i, j, k;
+ for( i = 0; i < N_reflections; i += 1 ) {
+
+ for ( j = 0; j < 3; j++ ) dist[j] = 0;
+
+ for ( j = 0; j < 3; j++ ) {
+ gsl_blas_ddot(reflections[i], c->axes[j],
+ &c->indices[i][j]);
+
+ for ( k = 0; k < 3; k++ ) {
+ dist[k] += gsl_vector_get(c->reciprocal[j], k) *
+ (c->indices[i][j] - round(c->indices[i][j]));
+ }
+ }
+
+ /* A reflection fits if the distance (in reciprocal space)
+ * between the observed and calculated reflection position
+ * is less than Indexfit */
+
+ if ( dist[0]*dist[0] + dist[1]*dist[1] + dist[2]*dist[2] <
+ IndexFit * IndexFit ) {
+ c->reflections[i] = 1;
+ c->n++;
+ } else {
+ c->reflections[i] = 0;
+ }
+ }
+
+ return 1;
+}
+
+
+static void print_asdf_cell(struct asdf_cell *cc)
+{
+ double a, b, c, alpha, beta, gamma, ab, bc, ca;
+
+ a = gsl_blas_dnrm2(cc->axes[0]);
+ b = gsl_blas_dnrm2(cc->axes[1]);
+ c = gsl_blas_dnrm2(cc->axes[2]);
+
+ gsl_blas_ddot(cc->axes[0], cc->axes[1], &ab);
+ gsl_blas_ddot(cc->axes[1], cc->axes[2], &bc);
+ gsl_blas_ddot(cc->axes[0], cc->axes[2], &ca);
+
+ alpha = acos(bc/b/c)/M_PI*180;
+ beta = acos(ca/a/c)/M_PI*180;
+ gamma = acos(ab/a/b)/M_PI*180;
+
+ //~ int i, j;
+ //~ for (i = 0; i < 3; i ++) {
+ //~ for (j = 0; j < 3; j ++) {
+ //~ printf("%f ", gsl_vector_get(cc.axes[i], j));
+ //~ }
+ //~ printf("\n");
+ //~ }
+
+ printf("%.2f %.2f %.2f %.2f %.2f %.2f %.0f %d \n", a, b, c,
+ alpha, beta, gamma,
+ cc->volume, cc->n);
+
+}
+
+
+/* Returns 0 when refinement doesn't converge (i.e. all fitting reflections
+ * lie in the same plane) */
+static int refine_asdf_cell(struct asdf_cell *c, gsl_vector **reflections,
+ int N_reflections, double IndexFit)
+{
+ gsl_matrix *X = gsl_matrix_alloc(c->n, 3);
+
+ gsl_vector *r[] = {gsl_vector_alloc(c->n),
+ gsl_vector_alloc(c->n),
+ gsl_vector_alloc(c->n)};
+
+ gsl_vector *res = gsl_vector_alloc(3);
+ gsl_matrix *cov = gsl_matrix_alloc (3, 3);
+ double chisq;
+
+ int i, j;
+ int n = 0;
+ for ( i = 0; i < N_reflections; i++ ) if ( c->reflections[i] == 1 )
+ {
+ for ( j = 0; j < 3; j++ ) {
+ gsl_matrix_set(X, n, j, round(c->indices[i][j]));
+ gsl_vector_set(r[j], n,
+ gsl_vector_get(reflections[i], j));
+ }
+ n++;
+ }
+
+ gsl_multifit_linear_workspace *work = gsl_multifit_linear_alloc(c->n,
+ 3);
+
+ for ( i = 0; i < 3; i++ ) {
+ gsl_multifit_linear (X, r[i], res, cov, &chisq, work);
+
+ for (j = 0; j < 3; j++ ) {
+ gsl_vector_set(c->reciprocal[j], i,
+ gsl_vector_get(res, j));
+ }
+ }
+
+ calc_reciprocal(c->reciprocal, c->axes);
+
+ double a[3];
+ for ( i = 0; i < 3; i++ ) {
+ a[i] = gsl_blas_dnrm2(c->axes[i]);
+ }
+
+ gsl_multifit_linear_free(work);
+ gsl_vector_free(res);
+ gsl_matrix_free(cov);
+ gsl_matrix_free(X);
+ for ( i = 0; i < 3; i++ ) {
+ gsl_vector_free(r[i]);
+ }
+
+ if ( fabs(a[0]) > 10000 || fabs(a[1]) > 10000 ||
+ fabs(a[2]) > 10000 || isnan(a[0]) )
+ {
+ return 0;
+ }
+
+ return 1;
+}
+
+
+static int reduce_asdf_cell(struct asdf_cell *cl)
+{
+ double a, b, c, alpha, beta, gamma, ab, bc, ca, bb, cc;
+
+ gsl_vector *va = gsl_vector_alloc(3);
+ gsl_vector *vb = gsl_vector_alloc(3);
+ gsl_vector *vc = gsl_vector_alloc(3);
+
+ int changed = 1;
+ while ( changed ) {
+ changed = 0;
+
+ gsl_vector_memcpy(va, cl->axes[0]);
+ gsl_vector_memcpy(vb, cl->axes[1]);
+ gsl_vector_memcpy(vc, cl->axes[2]);
+
+ a = gsl_blas_dnrm2(va);
+ b = gsl_blas_dnrm2(vb);
+ c = gsl_blas_dnrm2(vc);
+
+ gsl_blas_ddot(va, vb, &ab);
+ gsl_blas_ddot(vb, vc, &bc);
+ gsl_blas_ddot(vc, va, &ca);
+
+ alpha = acos(bc/b/c)/M_PI*180;
+ beta = acos(ca/a/c)/M_PI*180;
+ gamma = acos(ab/a/b)/M_PI*180;
+
+ if ( changed == 0 ) {
+
+ if ( gamma < 90 ) {
+ gsl_vector_scale(vb, -1);
+ gamma = 180 - gamma;
+ alpha = 180 - alpha;
+ }
+
+ gsl_vector_add(vb, va);
+ bb = gsl_blas_dnrm2(vb);
+ if ( bb < b ) {
+ b = bb;
+ if ( a < b ) {
+ gsl_vector_memcpy(cl->axes[1], vb);
+ } else {
+ gsl_vector_memcpy(cl->axes[1], va);
+ gsl_vector_memcpy(cl->axes[0], vb);
+ }
+ changed = 1;
+ }
+ }
+
+ if ( changed == 0 ) {
+
+ if ( beta < 90 ) {
+ gsl_vector_scale(vc, -1);
+ beta = 180 - beta;
+ alpha = 180 - alpha;
+ }
+
+ gsl_vector_add(vc, va);
+ cc = gsl_blas_dnrm2(vc);
+ if ( cc < c ) {
+ c = cc;
+ if ( b < c ) {
+ gsl_vector_memcpy(cl->axes[2], vc);
+ } else if ( a < c ) {
+ gsl_vector_memcpy(cl->axes[1], vc);
+ gsl_vector_memcpy(cl->axes[2], vb);
+ } else {
+ gsl_vector_memcpy(cl->axes[0], vc);
+ gsl_vector_memcpy(cl->axes[1], va);
+ gsl_vector_memcpy(cl->axes[2], vb);
+ }
+ changed = 1;
+ }
+ }
+
+ if ( changed == 0 ) {
+ if ( alpha < 90 ) {
+ gsl_vector_scale(vc, -1);
+ beta = 180 - beta;
+ alpha = 180 - alpha;
+ }
+
+ gsl_vector_add(vc, vb);
+ cc = gsl_blas_dnrm2(vc);
+ if ( cc < c ) {
+ c = cc;
+ if ( b < c ) {
+ gsl_vector_memcpy(cl->axes[2], vc);
+ } else if ( a < c ) {
+ gsl_vector_memcpy(cl->axes[1], vc);
+ gsl_vector_memcpy(cl->axes[2], vb);
+ } else {
+ gsl_vector_memcpy(cl->axes[0], vc);
+ gsl_vector_memcpy(cl->axes[1], va);
+ gsl_vector_memcpy(cl->axes[2], vb);
+ }
+ changed = 1;
+ }
+ }
+ }
+
+ cross_product(cl->axes[0], cl->axes[1], &vc);
+ gsl_blas_ddot(vc, cl->axes[2], &cl->volume);
+ if ( cl->volume < 0 ) {
+ gsl_vector_scale(cl->axes[2], -1);
+ cl->volume *= -1;
+ }
+
+ gsl_vector_free(va);
+ gsl_vector_free(vb);
+ gsl_vector_free(vc);
+
+ return 1;
+}
+
+
+static int check_cell_angles(gsl_vector *va, gsl_vector *vb, gsl_vector *vc,
+ double max_cos)
+{
+ double a, b, c, cosa, cosb, cosg, ab, bc, ca;
+
+ a = gsl_blas_dnrm2(va);
+ b = gsl_blas_dnrm2(vb);
+ c = gsl_blas_dnrm2(vc);
+
+ gsl_blas_ddot(va, vb, &ab);
+ gsl_blas_ddot(vb, vc, &bc);
+ gsl_blas_ddot(vc, va, &ca);
+
+ cosa = bc/b/c;
+ cosb = ca/a/c;
+ cosg = ab/a/b;
+
+ if ( fabs(cosa) > max_cos || fabs(cosb) > max_cos ||
+ fabs(cosg) > max_cos ) {
+ return 0;
+ }
+
+ return 1;
+}
+
+
+/* Returns min(t1.n, t2.n, t3.n) */
+static int find_acl(struct tvector t1, struct tvector t2, struct tvector t3)
+{
+ int i = t1.n, j = t2.n, k = t3.n;
+ if ( i <= j && i <= k ) return i;
+ if ( j <= i && j <= k ) return j;
+ if ( k <= i && k <= j ) return k;
+}
+
+
+static int create_cell(struct tvector tvec1, struct tvector tvec2,
+ struct tvector tvec3, struct asdf_cell *c,
+ double IndexFit, double volume_min, double volume_max,
+ gsl_vector **reflections, int N_reflections)
+{
+
+ double volume = calc_volume(tvec1.t, tvec2.t, tvec3.t);
+ if ( fabs(volume) < volume_min || fabs(volume) > volume_max ) return 0;
+
+ gsl_vector_memcpy(c->axes[0], tvec1.t);
+ gsl_vector_memcpy(c->axes[1], tvec2.t);
+ gsl_vector_memcpy(c->axes[2], tvec3.t);
+
+ c->volume = volume;
+ check_refl_fitting_cell(c, reflections, N_reflections, IndexFit);
+
+ if ( c->n < 6 ) return 0;
+
+ reduce_asdf_cell(c);
+
+ /* If one of the cell angles > 135 or < 45 return 0 */
+ if ( !check_cell_angles(c->axes[0], c->axes[1],
+ c->axes[2], 0.71) ) return 0;
+
+ /* Index reflections with new cell axes */
+ check_refl_fitting_cell(c, reflections, N_reflections, IndexFit);
+
+ /* Refine cell until the number of fitting
+ * reflections stops increasing */
+ int n = 0;
+ int cell_correct = 1;
+ while ( c->n - n > 0 && cell_correct ) {
+
+ n = c->n;
+ cell_correct = refine_asdf_cell(c, reflections, N_reflections,
+ IndexFit);
+ check_refl_fitting_cell(c, reflections, N_reflections,
+ IndexFit);
+ }
+
+ return cell_correct;
+}
+
+
+static int find_cell(struct tvector *tvectors, int N_tvectors, double IndexFit,
+ double volume_min, double volume_max, int n_max,
+ gsl_vector **reflections, int N_reflections,
+ struct asdf_cell *result)
+{
+ int i, j, k;
+
+ /* Only tvectors with the number of fitting reflections > acl are
+ * considered */
+ int acl = N_reflections < 18 ? 6 : N_reflections/3;
+
+ struct asdf_cell *c = asdf_cell_new(N_reflections);
+ if (c == NULL) {
+ ERROR("Failed to allocate asdf_cell in find_cell!\n");
+ return 0;
+ }
+
+ /* Traversing a 3d array in slices perpendicular to the main diagonal */
+ int sl;
+ for ( sl = 0; sl < 3 * N_tvectors - 1; sl++ ) {
+
+ int i_min = sl < 2 * N_tvectors ? 0 : sl - 2 * N_tvectors;
+ int i_max = sl < N_tvectors ? sl : N_tvectors;
+
+ for ( i = i_min; i < i_max; i++) if (tvectors[i].n > acl ) {
+
+ int j_min = sl - N_tvectors - 2 * i - 1 < 0 ?
+ i + 1 : sl - N_tvectors - i;
+ int j_max = sl - N_tvectors - i < 0 ?
+ sl - i : N_tvectors;
+
+ for ( j = j_min; j < j_max; j++ )
+ if ( tvectors[j].n > acl ) {
+
+ k = sl - i - j - 1;
+
+ if ( k > j && tvectors[k].n > acl &&
+ check_cell_angles(tvectors[i].t,
+ tvectors[j].t,
+ tvectors[k].t, 0.99) )
+ {
+
+ if ( !create_cell(tvectors[i],
+ tvectors[j],
+ tvectors[k],
+ c, IndexFit,
+ volume_min,
+ volume_max,
+ reflections,
+ N_reflections) )
+ {
+ break;
+ }
+
+ acl = find_acl(tvectors[i],
+ tvectors[j],
+ tvectors[k]);
+ c->acl = acl;
+ c->n_max = n_max;
+
+ reduce_asdf_cell(c);
+
+ /* If the new cell has more fitting
+ * reflections save it to result */
+ if ( result->n < c->n ) {
+ asdf_cell_memcpy(result, c);
+ }
+ acl++;
+
+ if ( acl > n_max ) break;
+ if ( tvectors[j].n <= acl ||
+ tvectors[i].n <= acl ) break;
+ }
+ }
+ if ( acl > n_max ) break;
+ if ( tvectors[i].n <= acl ) break;
+ }
+ if ( acl > n_max ) break;
+ }
+
+ asdf_cell_free(c);
+
+ if ( result->n ) return 1;
+ return 0;
+}
+
+
+static void shuffle_triplets(int **triplets, int n)
+{
+ int i, j;
+ int t[3];
+ for ( i = 0; i < n - 1; i++ ) {
+ j = i + rand() / (RAND_MAX / (n - i) + 1);
+ memcpy(t, triplets[j], 3 * sizeof(int));
+ memcpy(triplets[j], triplets[i], 3 * sizeof(int));
+ memcpy(triplets[i], t, 3 * sizeof(int));
+ }
+}
+
+
+static double angle_between_gsl(gsl_vector *a, gsl_vector *b)
+{
+ double ab;
+ gsl_blas_ddot(a, b, &ab);
+ return acos(ab/gsl_blas_dnrm2(a)/gsl_blas_dnrm2(b)) * 180 / M_PI;
+}
+
+
+static int index_refls(gsl_vector **reflections, int N_reflections,
+ double d_max, double volume_min, double volume_max,
+ double LevelFit, double IndexFit, int i_max,
+ struct asdf_cell *c, struct fftw_vars fftw)
+{
+
+ int i, j, k, l, n;
+
+ /* Number of triplets = c_n^3 if n - number of reflections */
+ int N_triplets = N_reflections * (N_reflections - 1) *
+ (N_reflections - 2) / 6;
+
+ int **triplets = malloc(N_triplets * sizeof(int *));
+ if (triplets == NULL) {
+ ERROR("Failed to allocate triplets in index_refls!\n");
+ return 0;
+ }
+
+ l = 0;
+ for ( i = 0; i < N_reflections; i++ ) {
+ for ( j = i + 1; j < N_reflections; j++ ) {
+ for ( k = j + 1; k < N_reflections; k++ ) {
+ triplets[l] = malloc(3 * sizeof(int));
+ if (triplets[l] == NULL) {
+ ERROR("Failed to allocate triplets "
+ " in index_refls!\n");
+ return 0;
+ }
+
+ triplets[l][0] = i;
+ triplets[l][1] = j;
+ triplets[l][2] = k;
+ l++;
+ }
+ }
+ }
+
+ /* Triplets are processed in a random sequence if N_triplets > 10000 */
+ if ( N_reflections > 40 ) shuffle_triplets(triplets, N_triplets);
+
+ gsl_vector *normal = gsl_vector_alloc(3);
+
+ double projections[N_reflections];
+ double ds;
+
+ int *fits = malloc(N_reflections * sizeof(int));
+ if (fits == NULL) {
+ ERROR("Failed to allocate fits in index_refls!\n");
+ return 0;
+ }
+
+ if ( i_max > N_triplets ) i_max = N_triplets;
+
+ struct tvector *tvectors = malloc(i_max * sizeof(struct tvector));
+ if (tvectors == NULL) {
+ ERROR("Failed to allocate tvectors in index_refls!\n");
+ return 0;
+ }
+
+ int N_tvectors = 0;
+
+ int n_max = 0; // maximum number of reflections fitting one of tvectors
+
+ for ( i = 0; i < i_max; i++ ) {
+ if ( calc_normal(reflections[triplets[i][0]],
+ reflections[triplets[i][1]],
+ reflections[triplets[i][2]],
+ normal) )
+ {
+
+ /* Calculate projections of reflections to normal */
+ for ( k = 0; k < N_reflections; k++ ) {
+ gsl_blas_ddot(normal, reflections[k],
+ &projections[k]);
+ }
+
+ /* Find ds - period in 1d lattice of projections */
+ ds = find_ds_fft(projections, N_reflections, d_max,
+ fftw);
+
+ /* Refine ds, write 1 to fits[i] if reflections[i]
+ * fits ds */
+ ds = refine_ds(projections, N_reflections, ds, LevelFit,
+ fits);
+
+ /* n - number of reflections fitting ds */
+ n = check_refl_fitting_ds(projections, N_reflections,
+ ds, LevelFit);
+
+ /* normal/ds - possible direct vector */
+ gsl_vector_scale(normal, 1/ds);
+
+ if ( n > N_reflections/3 && n > 6 ) {
+
+ tvectors[N_tvectors] = tvector_new(N_reflections);
+
+ gsl_vector_memcpy(tvectors[N_tvectors].t,
+ normal);
+ memcpy(tvectors[N_tvectors].fits, fits,
+ N_reflections * sizeof(int));
+
+ tvectors[N_tvectors].n = n;
+
+ N_tvectors++;
+
+ if (n > n_max) n_max = n;
+ }
+ }
+
+ if ( (i != 0 && i % 10000 == 0) || i == i_max - 1 ) {
+
+ /* Sort tvectors by length */
+ qsort(tvectors, N_tvectors, sizeof(struct tvector),
+ compare_tvectors);
+
+ /* Three shortest independent tvectors with t.n > acl
+ * determine the final cell. acl is selected for the
+ * solution with the maximum number of fitting
+ * reflections */
+
+ find_cell(tvectors, N_tvectors, IndexFit, volume_min,
+ volume_max, n_max, reflections,
+ N_reflections, c);
+
+ if ( c->n > 4 * n_max / 5 ) {
+ break;
+ }
+ }
+ }
+ free(fits);
+
+ for ( i = 0; i < N_tvectors; i++ ) {
+ tvector_free(tvectors[i]);
+ }
+ free(tvectors);
+
+ for ( i = 0; i < N_triplets; i++ ) {
+ free(triplets[i]);
+ }
+ free(triplets);
+
+ gsl_vector_free(normal);
+
+ if ( c->n ) return 1;
+
+ return 0;
+}
+
+
+int run_asdf(struct image *image, IndexingPrivate *ipriv)
+{
+ int i, j;
+
+ double LevelFit = 1./1000;
+ double IndexFit = 1./500;
+ double d_max = 220.; // thrice the maximum expected axis length
+ double volume_min = 100.;
+ double volume_max = 1000000.;
+
+ int i_max = 10000; // maximum number of triplets
+
+ struct asdf_private *dp = (struct asdf_private *)ipriv;
+
+ if ( dp->indm & INDEXING_CHECK_CELL_AXES ) {
+ double volume = cell_get_volume(dp->template);
+ double vtol = (dp->ltl[0] + dp->ltl[1] + dp->ltl[2]) / 100 +
+ dp->ltl[3] / 180 * M_PI;
+ volume_min = volume * (1 - vtol);
+ volume_max = volume * (1 + vtol);
+ }
+
+ int N_reflections = image_feature_count(image->features);
+ gsl_vector *reflections[N_reflections];
+
+ for ( i = 0; i < N_reflections; i++ ) {
+ struct imagefeature *f;
+
+ f = image_get_feature(image->features, i);
+ if ( f == NULL ) continue;
+
+ reflections[i] = gsl_vector_alloc(3);
+ gsl_vector_set(reflections[i], 0, f->rx/1e10);
+ gsl_vector_set(reflections[i], 1, f->ry/1e10);
+ gsl_vector_set(reflections[i], 2, f->rz/1e10);
+ }
+
+ struct asdf_cell *c = asdf_cell_new(N_reflections);
+ if (c == NULL) {
+ ERROR("Failed to allocate asdf_cell in run_asdf!\n");
+ return 0;
+ }
+
+ if ( N_reflections == 0 ) return 0;
+
+ j = index_refls(reflections, N_reflections, d_max, volume_min,
+ volume_max, LevelFit, IndexFit, i_max, c, dp->fftw);
+
+ for ( i = 0; i < N_reflections; i++ ) {
+ gsl_vector_free(reflections[i]);
+ }
+
+ if ( j ) {
+ UnitCell *uc;
+ uc = cell_new();
+
+ cell_set_cartesian(uc, gsl_vector_get(c->axes[0], 0) * 1e-10,
+ gsl_vector_get(c->axes[0], 1) * 1e-10,
+ gsl_vector_get(c->axes[0], 2) * 1e-10,
+ gsl_vector_get(c->axes[1], 0) * 1e-10,
+ gsl_vector_get(c->axes[1], 1) * 1e-10,
+ gsl_vector_get(c->axes[1], 2) * 1e-10,
+ gsl_vector_get(c->axes[2], 0) * 1e-10,
+ gsl_vector_get(c->axes[2], 1) * 1e-10,
+ gsl_vector_get(c->axes[2], 2) * 1e-10);
+
+ if ( check_cell(dp, image, uc) ) {
+ asdf_cell_free(c);
+ cell_free(uc);
+ return 1;
+ }
+
+ cell_free(uc);
+ }
+
+ asdf_cell_free(c);
+ return 0;
+}
+
+
+IndexingPrivate *asdf_prepare(IndexingMethod *indm, UnitCell *cell,
+ struct detector *det, float *ltl)
+{
+ struct asdf_private *dp;
+ int need_cell = 0;
+
+ if ( *indm & INDEXING_CHECK_CELL_COMBINATIONS ) need_cell = 1;
+ if ( *indm & INDEXING_CHECK_CELL_AXES ) need_cell = 1;
+
+ if ( need_cell && !cell_has_parameters(cell) ) {
+ ERROR("Altering your asdf flags because cell parameters were"
+ " not provided.\n");
+ *indm &= ~INDEXING_CHECK_CELL_COMBINATIONS;
+ *indm &= ~INDEXING_CHECK_CELL_AXES;
+ }
+
+ /* Flags that asdf knows about */
+ *indm &= INDEXING_METHOD_MASK | INDEXING_CHECK_CELL_COMBINATIONS
+ | INDEXING_CHECK_CELL_AXES | INDEXING_CHECK_PEAKS;
+
+ dp = malloc(sizeof(struct asdf_private));
+ if ( dp == NULL ) return NULL;
+
+ dp->ltl = ltl;
+ dp->template = cell;
+ dp->indm = *indm;
+
+ dp->fftw = fftw_vars_new();
+
+ return (IndexingPrivate *)dp;
+}
+
+
+void asdf_cleanup(IndexingPrivate *pp)
+{
+ struct asdf_private *p;
+ p = (struct asdf_private *)pp;
+ fftw_vars_free(p->fftw);
+ free(p);
+}
diff --git a/libcrystfel/src/asdf.h b/libcrystfel/src/asdf.h
new file mode 100644
index 00000000..402636d1
--- /dev/null
+++ b/libcrystfel/src/asdf.h
@@ -0,0 +1,28 @@
+
+
+#ifndef ASDF_H
+#define ASDF_H
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include "index.h"
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+extern int run_asdf(struct image *image, IndexingPrivate *ipriv);
+
+extern IndexingPrivate *asdf_prepare(IndexingMethod *indm,
+ UnitCell *cell, struct detector *det,
+ float *ltl);
+
+extern void asdf_cleanup(IndexingPrivate *pp);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif /* DIRAX_H */
diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c
index 9676aade..75419c99 100644
--- a/libcrystfel/src/cell-utils.c
+++ b/libcrystfel/src/cell-utils.c
@@ -1576,3 +1576,31 @@ int forbidden_reflection(UnitCell *cell,
return 0;
}
+
+
+/* Returns cell volume in A^3 */
+double cell_get_volume(UnitCell *cell)
+{
+ double asx, asy, asz;
+ double bsx, bsy, bsz;
+ double csx, csy, csz;
+ struct rvec aCb;
+ double rec_volume;
+
+ if ( cell_get_reciprocal(cell, &asx, &asy, &asz,
+ &bsx, &bsy, &bsz,
+ &csx, &csy, &csz) ) {
+ ERROR("Couldn't get reciprocal cell.\n");
+ return 0;
+ }
+
+ /* "a" cross "b" */
+ aCb.u = asy*bsz - asz*bsy;
+ aCb.v = - (asx*bsz - asz*bsx);
+ aCb.w = asx*bsy - asy*bsx;
+
+ /* "a cross b" dot "c" */
+ rec_volume = (aCb.u*csx + aCb.v*csy + aCb.w*csz)/1e30;
+
+ return 1/rec_volume;
+}
diff --git a/libcrystfel/src/cell-utils.h b/libcrystfel/src/cell-utils.h
index 0b900096..98bc667c 100644
--- a/libcrystfel/src/cell-utils.h
+++ b/libcrystfel/src/cell-utils.h
@@ -74,6 +74,8 @@ extern LatticeType lattice_from_str(const char *s);
extern int forbidden_reflection(UnitCell *cell,
signed int h, signed int k, signed int l);
+extern double cell_get_volume(UnitCell *cell);
+
#ifdef __cplusplus
}
#endif
diff --git a/libcrystfel/src/index.c b/libcrystfel/src/index.c
index 6f046ab6..6b36af3e 100644
--- a/libcrystfel/src/index.c
+++ b/libcrystfel/src/index.c
@@ -44,6 +44,7 @@
#include "utils.h"
#include "peaks.h"
#include "dirax.h"
+#include "asdf.h"
#include "mosflm.h"
#include "xds.h"
#include "detector.h"
@@ -92,12 +93,16 @@ IndexingPrivate **prepare_indexing(IndexingMethod *indm, UnitCell *cell,
iprivs[n] = dirax_prepare(&indm[n], cell, det, ltl);
break;
+ case INDEXING_ASDF :
+ iprivs[n] = asdf_prepare(&indm[n], cell, det, ltl);
+ break;
+
case INDEXING_MOSFLM :
iprivs[n] = mosflm_prepare(&indm[n], cell, det, ltl);
break;
case INDEXING_XDS :
- iprivs[n] = xds_prepare(&indm[n], cell, det, ltl);
+ iprivs[n] = xds_prepare(&indm[n], cell, det, ltl);
break;
case INDEXING_REAX :
@@ -167,6 +172,10 @@ void cleanup_indexing(IndexingMethod *indms, IndexingPrivate **privs)
dirax_cleanup(privs[n]);
break;
+ case INDEXING_ASDF :
+ asdf_cleanup(privs[n]);
+ break;
+
case INDEXING_MOSFLM :
mosflm_cleanup(privs[n]);
break;
@@ -236,6 +245,10 @@ static int try_indexer(struct image *image, IndexingMethod indm,
return run_dirax(image, ipriv);
break;
+ case INDEXING_ASDF :
+ return run_asdf(image, ipriv);
+ break;
+
case INDEXING_MOSFLM :
return run_mosflm(image, ipriv);
break;
@@ -272,11 +285,6 @@ void index_pattern(struct image *image,
if ( indms == NULL ) return;
- if ( image_feature_count(image->features) > 10000 ) {
- STATUS("WARNING: The number of peaks is very large for '%s'.\n",
- image->filename);
- }
-
map_all_peaks(image);
image->crystals = NULL;
image->n_crystals = 0;
@@ -373,6 +381,10 @@ char *indexer_str(IndexingMethod indm)
strcpy(str, "dirax");
break;
+ case INDEXING_ASDF :
+ strcpy(str, "asdf");
+ break;
+
case INDEXING_MOSFLM :
strcpy(str, "mosflm");
break;
@@ -451,13 +463,16 @@ IndexingMethod *build_indexer_list(const char *str)
if ( strcmp(methods[i], "dirax") == 0) {
list[++nmeth] = INDEXING_DEFAULTS_DIRAX;
+ } else if ( strcmp(methods[i], "asdf") == 0) {
+ list[++nmeth] = INDEXING_DEFAULTS_ASDF;
+
} else if ( strcmp(methods[i], "mosflm") == 0) {
list[++nmeth] = INDEXING_DEFAULTS_MOSFLM;
} else if ( strcmp(methods[i], "grainspotter") == 0) {
list[++nmeth] = INDEXING_DEFAULTS_GRAINSPOTTER;
- } else if ( strcmp(methods[i], "xds") == 0) {
+ } else if ( strcmp(methods[i], "xds") == 0) {
list[++nmeth] = INDEXING_DEFAULTS_XDS;
} else if ( strcmp(methods[i], "reax") == 0) {
diff --git a/libcrystfel/src/index.h b/libcrystfel/src/index.h
index dcb21858..ff1cd5df 100644
--- a/libcrystfel/src/index.h
+++ b/libcrystfel/src/index.h
@@ -41,6 +41,9 @@
#define INDEXING_DEFAULTS_DIRAX (INDEXING_DIRAX | INDEXING_CHECK_PEAKS \
| INDEXING_CHECK_CELL_COMBINATIONS)
+#define INDEXING_DEFAULTS_ASDF (INDEXING_ASDF | INDEXING_CHECK_PEAKS \
+ | INDEXING_CHECK_CELL_COMBINATIONS)
+
#define INDEXING_DEFAULTS_MOSFLM (INDEXING_MOSFLM | INDEXING_CHECK_PEAKS \
| INDEXING_CHECK_CELL_COMBINATIONS \
| INDEXING_USE_LATTICE_TYPE)
@@ -97,6 +100,7 @@ typedef enum {
INDEXING_XDS = 5,
INDEXING_SIMULATION = 6,
INDEXING_DEBUG = 7,
+ INDEXING_ASDF = 8,
/* Bits at the top of the IndexingMethod are flags which modify the
* behaviour of the indexer. */
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-05 21:44:08 +0000