aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--libcrystfel/Makefile.am4
-rw-r--r--libcrystfel/src/asdf.c1045
-rw-r--r--libcrystfel/src/asdf.h28
-rw-r--r--libcrystfel/src/index.c29
-rw-r--r--libcrystfel/src/index.h4
5 files changed, 1101 insertions, 9 deletions
diff --git a/libcrystfel/Makefile.am b/libcrystfel/Makefile.am
index 2cdfba10..79699e4a 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..0735a14f
--- /dev/null
+++ b/libcrystfel/src/asdf.c
@@ -0,0 +1,1045 @@
+#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 asdf_private {
+ IndexingMethod indm;
+ float *ltl;
+ UnitCell *template;
+};
+
+
+/* Possible direct vector */
+struct tvector {
+ gsl_vector *t;
+ int n; // number of fitting reflections
+ int *fits;
+};
+
+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);
+
+ return 1;
+}
+
+static struct asdf_cell asdf_cell_new(int n) {
+
+ struct asdf_cell c;
+
+ 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);
+
+ c.indices = malloc(sizeof(int *) * n);
+ for (i = 0; i < n; i ++) {
+ c.indices[i] = malloc(sizeof(int) * 3);
+ }
+
+ 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);
+
+ memcpy(dest->indices, src->indices, sizeof(int *) * n);
+
+ for (i = 0; i < n; i ++) {
+ memcpy(dest->indices[i], src->indices[i], sizeof(int) * 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 += 1) {
+ 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 += 1) {
+ 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 += 1) {
+ 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) {
+
+ 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;
+
+ int N = 1024; // number of points in fft calculation
+
+ double *in = malloc(sizeof(double) * N);
+ fftw_complex *out = fftw_malloc(sizeof (fftw_complex) * N);
+
+ for (i = 0; i < N; i++) {
+ in[i] = 0;
+ }
+
+ for (i = 0; i < n; i++) {
+ in[(int)((projections_sorted[i] - projections_sorted[0]) /
+ (projections_sorted[n - 1] - projections_sorted[0]) * N)] ++;
+ }
+
+ fftw_plan p = fftw_plan_dft_r2c_1d(N, in, out, FFTW_MEASURE);
+ fftw_execute(p);
+ fftw_destroy_plan(p);
+
+ 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;
+
+ free(in);
+ fftw_free(out);
+ 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 += 1) {
+ 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 += 1) dist[j] = 0;
+
+ for (j = 0; j < 3; j += 1) {
+ gsl_blas_ddot(reflections[i], c->axes[j], &c->indices[i][j]);
+
+ for (k = 0; k < 3; k += 1) {
+ 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 += 1;
+ } 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);
+
+}
+
+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 += 1) {
+ if (c.reflections[i] == 1) {
+ for (j = 0; j < 3; j += 1) {
+ gsl_matrix_set(X, n, j, round(c.indices[i][j]));
+ gsl_vector_set(r[j], n, gsl_vector_get(reflections[i], j));
+ }
+ n += 1;
+ }
+ }
+
+ gsl_multifit_linear_workspace *work = gsl_multifit_linear_alloc(c.n, 3);
+
+ for (i = 0; i < 3; i += 1) {
+ gsl_multifit_linear (X, r[i], res, cov, &chisq, work);
+
+ for (j = 0; j < 3; j += 1) {
+ 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 += 1) {
+ 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 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, n, cell_correct;
+
+ double volume;
+
+ // only tvectors with 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);
+
+ // 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) ) {
+
+ volume = calc_volume(tvectors[i].t,
+ tvectors[j].t,
+ tvectors[k].t);
+
+ if (fabs(volume) > volume_min && fabs(volume) < volume_max){
+
+ gsl_vector_memcpy(c.axes[0], tvectors[i].t);
+ gsl_vector_memcpy(c.axes[1], tvectors[j].t);
+ gsl_vector_memcpy(c.axes[2], tvectors[k].t);
+
+ c.volume = volume;
+ check_refl_fitting_cell(&c, reflections, N_reflections,
+ IndexFit);
+
+ if (c.n >= 6) {
+ reduce_asdf_cell(&c);
+
+ // if one of the cell angles > 135 or < 45 do not continue
+ if (!check_cell_angles(c.axes[0], c.axes[1],
+ c.axes[2], 0.71)) break;
+
+ // index reflections with new cell axes
+ check_refl_fitting_cell(&c, reflections,
+ N_reflections, IndexFit);
+
+ acl = find_acl(tvectors[i],
+ tvectors[j],
+ tvectors[k]);
+
+ c.acl = acl;
+ c.n_max = n_max;
+
+ // refine cell until number of fitting reflections
+ // stops increasing
+ n = 0;
+ cell_correct = 1;
+ while (c.n - n && cell_correct) {
+ n = c.n;
+ cell_correct = refine_asdf_cell(c, reflections,
+ N_reflections,
+ IndexFit);
+
+ check_refl_fitting_cell(&c, reflections,
+ N_reflections, IndexFit);
+ }
+
+ if (cell_correct) {
+ reduce_asdf_cell(&c);
+ 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) {
+
+ 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 *));
+ 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));
+
+ 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 (i_max > N_triplets) i_max = N_triplets;
+
+ struct tvector *tvectors = malloc(i_max * sizeof(struct tvector));
+ 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);
+
+ // 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);
+
+ if (c->n) return 1;
+
+ return 0;
+}
+
+double cell_get_volume(UnitCell *cell)
+{
+ double asx, asy, asz;
+ double bsx, bsy, bsz;
+ double csx, csy, csz;
+ struct rvec aCb;
+ double volume;
+
+ if ( cell_get_reciprocal(cell, &asx, &asy, &asz,
+ &bsx, &bsy, &bsz,
+ &csx, &csy, &csz) ) {
+ ERROR("Couldn't get cell cartesian.\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" */
+ volume = (aCb.u*csx + aCb.v*csy + aCb.w*csz)/1e30;
+
+ return 1/volume;
+}
+
+int run_asdf(struct image *image, IndexingPrivate *ipriv) {
+ int i;
+
+ 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);
+ volume_min = volume * 0.95;
+ volume_max = volume * 1.05;
+ }
+
+ 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 (N_reflections == 0) return 0;
+
+ i = index_refls(reflections, N_reflections, d_max, volume_min, volume_max,
+ LevelFit, IndexFit, i_max, &c);
+ //~
+
+ for (i = 0; i < N_reflections; i ++) {
+ gsl_vector_free(reflections[i]);
+ }
+
+ if (i) {
+ 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) ) {
+ cell_free(uc);
+ return 1;
+ }
+
+ cell_free(uc);
+ }
+
+ 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;
+
+ return (IndexingPrivate *)dp;
+}
+
+
+void asdf_cleanup(IndexingPrivate *pp)
+{
+ struct asdf_private *p;
+ p = (struct asdf_private *)pp;
+ 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/index.c b/libcrystfel/src/index.c
index c3f7c35d..97187655 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"
@@ -79,12 +80,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 :
@@ -150,6 +155,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;
@@ -215,6 +224,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;
@@ -248,11 +261,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;
@@ -349,6 +357,10 @@ char *indexer_str(IndexingMethod indm)
strcpy(str, "dirax");
break;
+ case INDEXING_ASDF :
+ strcpy(str, "asdf");
+ break;
+
case INDEXING_MOSFLM :
strcpy(str, "mosflm");
break;
@@ -423,13 +435,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 0e105aeb..709b7507 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)
@@ -95,6 +98,7 @@ typedef enum {
INDEXING_GRAINSPOTTER = 4,
INDEXING_XDS = 5,
INDEXING_SIMULATION = 6,
+ INDEXING_ASDF = 7,
/* 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 17:30:09 +0000