aboutsummaryrefslogtreecommitdiff
path: root/libcrystfel/src/indexers
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2020-08-07 17:43:18 +0200
committerThomas White <taw@physics.org>2020-08-07 18:07:07 +0200
commitc0b01532441407dc97eaa9d44b540f1dd0223990 (patch)
treea1526e8bb84d38e8e6dfab05d2d95e1a5b5a3d10 /libcrystfel/src/indexers
parent08327436744a05e68daf1676f0fa4a82fb74408f (diff)
Move indexers out of API
Diffstat (limited to 'libcrystfel/src/indexers')
-rw-r--r--libcrystfel/src/indexers/asdf.c1239
-rw-r--r--libcrystfel/src/indexers/asdf.h56
-rw-r--r--libcrystfel/src/indexers/dirax.c683
-rw-r--r--libcrystfel/src/indexers/dirax.h52
-rw-r--r--libcrystfel/src/indexers/felix.c963
-rw-r--r--libcrystfel/src/indexers/felix.h52
-rw-r--r--libcrystfel/src/indexers/mosflm.c945
-rw-r--r--libcrystfel/src/indexers/mosflm.h56
-rw-r--r--libcrystfel/src/indexers/pinkindexer.c646
-rw-r--r--libcrystfel/src/indexers/pinkindexer.h47
-rw-r--r--libcrystfel/src/indexers/taketwo.c2370
-rw-r--r--libcrystfel/src/indexers/taketwo.h47
-rw-r--r--libcrystfel/src/indexers/xds.c541
-rw-r--r--libcrystfel/src/indexers/xds.h58
-rw-r--r--libcrystfel/src/indexers/xgandalf.c497
-rw-r--r--libcrystfel/src/indexers/xgandalf.h51
16 files changed, 8303 insertions, 0 deletions
diff --git a/libcrystfel/src/indexers/asdf.c b/libcrystfel/src/indexers/asdf.c
new file mode 100644
index 00000000..1536a109
--- /dev/null
+++ b/libcrystfel/src/indexers/asdf.c
@@ -0,0 +1,1239 @@
+/*
+ * asdf.c
+ *
+ * Alexandra's Superior Direction Finder, or
+ * Algorithm Similar to DirAx, FFT-based
+ *
+ * Copyright © 2014-2020 Deutsches Elektronen-Synchrotron DESY,
+ * a research centre of the Helmholtz Association.
+ *
+ * Authors:
+ * 2014-2015 Alexandra Tolstikova <alexandra.tolstikova@desy.de>
+ * 2015,2017 Thomas White <taw@physics.org>
+ *
+ * This file is part of CrystFEL.
+ *
+ * CrystFEL is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * CrystFEL is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with CrystFEL. If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#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 <gsl/gsl_fit.h>
+
+#include "image.h"
+#include "utils.h"
+#include "peaks.h"
+#include "cell-utils.h"
+#include "asdf.h"
+#include "cell.h"
+
+/**
+ * \file asdf.h
+ */
+
+#ifdef HAVE_FFTW
+
+#define FFTW_NO_Complex /* Please use "double[2]", not C99 "complex",
+ * despite complex.h possibly already being
+ * included. For more information, refer to:
+ * http://www.fftw.org/doc/Complex-numbers.html */
+
+#include <fftw3.h>
+
+struct fftw_vars {
+ int N;
+ fftw_plan p;
+ double *in;
+ fftw_complex *out;
+};
+
+
+struct asdf_private {
+ IndexingMethod indm;
+ 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 void 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 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 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);
+}
+
+
+static double max(double a, double b, double c)
+{
+ double m = a;
+ if ( m < b ) m = b;
+ if ( m < c ) m = c;
+ return m;
+}
+
+
+/* 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));
+ if ( (k>=N) || (k<0) ) {
+ ERROR("Bad k value in find_ds_fft() (k=%i, N=%i)\n",
+ k, N);
+ return -1.0;
+ }
+ 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;
+}
+
+
+/* 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;
+
+ int n = 0;
+ while ( changed ) {
+ n += 1;
+ 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;
+ }
+ }
+
+ if (n > 30) changed = 0;
+ }
+
+ 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;
+ ERROR("This point never reached!\n");
+ abort();
+}
+
+
+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;
+}
+
+
+void swap(int *a, int *b) {
+ int c = *a;
+ *a = *b;
+ *b = c;
+}
+
+
+/* Generate triplets of integers < N_reflections in random sequence */
+static int **generate_triplets(int N_reflections, int N_triplets_max, int *N)
+{
+ 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;
+
+ if ( N_triplets > N_triplets_max || N_reflections > 1000 ) {
+ N_triplets = N_triplets_max;
+ }
+ *N = N_triplets;
+
+ int **triplets = malloc(N_triplets * sizeof(int *));
+
+ if (triplets == NULL) {
+ ERROR("Failed to allocate triplets in generate_triplets!\n");
+ return 0;
+ }
+
+ int is_in_triplets;
+ n = 0;
+
+ while ( n < N_triplets ) {
+ /* Generate three different integer numbers < N_reflections */
+ i = rand() % N_reflections;
+ j = i;
+ k = i;
+ while ( j == i ) {
+ j = rand() % N_reflections;
+ }
+ while ( k == i || k == j ) {
+ k = rand() % N_reflections;
+ }
+
+ /* Sort (i, j, k) */
+ if ( i > k ) swap(&i, &k);
+ if ( i > j ) swap(&i, &j);
+ if ( j > k ) swap(&j, &k);
+
+ /* Check if it's already in triplets[] */
+ is_in_triplets = 0;
+ for ( l = 0; l < n; l++ ) {
+ if ( triplets[l][0] == i &&
+ triplets[l][1] == j &&
+ triplets[l][2] == k ) {
+ is_in_triplets = 1;
+ break;
+ }
+ }
+
+ if ( is_in_triplets == 0 ) {
+ triplets[n] = malloc(3 * sizeof(int));
+ if (triplets[n] == NULL) {
+ ERROR("Failed to allocate triplets "
+ " in generate_triplets!\n");
+ return 0;
+ }
+ triplets[n][0] = i;
+ triplets[n][1] = j;
+ triplets[n][2] = k;
+
+ n++;
+ }
+ }
+
+ return triplets;
+}
+
+
+static int index_refls(gsl_vector **reflections, int N_reflections,
+ double d_max, double volume_min, double volume_max,
+ double LevelFit, double IndexFit, int N_triplets_max,
+ struct asdf_cell *c, struct fftw_vars fftw)
+{
+
+ int i, k, n;
+
+ int N_triplets;
+ int **triplets = generate_triplets(N_reflections, N_triplets_max,
+ &N_triplets);
+ if ( N_triplets == 0 ) {
+ return 0;
+ }
+
+ 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;
+ }
+
+ struct tvector *tvectors = malloc(N_triplets * 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 < N_triplets; 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);
+ if ( ds < 0.0 ) {
+ ERROR("find_ds_fft() failed.\n");
+ return 0;
+ }
+
+ /* 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 == N_triplets - 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, void *ipriv)
+{
+ int i, j;
+
+ double LevelFit = 1./1000;
+ double IndexFit = 1./500;
+ double d_max = 1000.; // thrice the maximum expected axis length
+ double volume_min = 100.;
+ double volume_max = 100000000.;
+
+ int N_triplets_max = 10000; // maximum number of triplets
+
+ struct asdf_private *dp = (struct asdf_private *)ipriv;
+
+ if ( dp->indm & INDEXING_USE_CELL_PARAMETERS ) {
+
+ double a, b, c, gamma, beta, alpha;
+ cell_get_parameters(dp->template, &a, &b, &c,
+ &alpha, &beta, &gamma);
+
+ d_max = max(a, b, c) * 3 * 1e10;
+
+ double volume = cell_get_volume(dp->template) / 1e30;
+
+ /* Divide volume constraints by number of lattice points per
+ * unit cell since asdf always finds primitive cell */
+ int latt_points_per_uc = 1;
+ char centering = cell_get_centering(dp->template);
+ if ( centering == 'A' ||
+ centering == 'B' ||
+ centering == 'C' ||
+ centering == 'I' ) latt_points_per_uc = 2;
+ else if ( centering == 'F' ) latt_points_per_uc = 4;
+
+ volume_min = volume * 0.9/latt_points_per_uc;
+ volume_max = volume * 1.1/latt_points_per_uc;
+ }
+
+ int n = image_feature_count(image->features);
+ int N_reflections = 0;
+ gsl_vector *reflections[n];
+
+ for ( i=0; i<n; i++ ) {
+ struct imagefeature *f;
+
+ f = image_get_feature(image->features, i);
+ if ( f == NULL ) continue;
+
+ reflections[N_reflections] = gsl_vector_alloc(3);
+ gsl_vector_set(reflections[N_reflections], 0, f->rx/1e10);
+ gsl_vector_set(reflections[N_reflections], 1, f->ry/1e10);
+ gsl_vector_set(reflections[N_reflections], 2, f->rz/1e10);
+ N_reflections++;
+ }
+
+ 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, N_triplets_max, c,
+ dp->fftw);
+
+ for ( i = 0; i < N_reflections; i++ ) {
+ gsl_vector_free(reflections[i]);
+ }
+
+ if ( j ) {
+
+ UnitCell *uc;
+ Crystal *cr;
+ 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);
+
+ cr = crystal_new();
+ if ( cr == NULL ) {
+ ERROR("Failed to allocate crystal.\n");
+ return 0;
+ }
+ crystal_set_cell(cr, uc);
+ image_add_crystal(image, cr);
+ asdf_cell_free(c);
+ return 1;
+
+ }
+
+ asdf_cell_free(c);
+ return 0;
+}
+
+
+/**
+ * Prepare the ASDF indexing algorithm
+ */
+void *asdf_prepare(IndexingMethod *indm, UnitCell *cell)
+{
+ struct asdf_private *dp;
+
+ /* Flags that asdf knows about */
+ *indm &= INDEXING_METHOD_MASK | INDEXING_USE_CELL_PARAMETERS;
+
+ dp = malloc(sizeof(struct asdf_private));
+ if ( dp == NULL ) return NULL;
+
+ dp->template = cell;
+ dp->indm = *indm;
+
+ dp->fftw = fftw_vars_new();
+
+ return (void *)dp;
+}
+
+
+void asdf_cleanup(void *pp)
+{
+ struct asdf_private *p;
+ p = (struct asdf_private *)pp;
+ fftw_vars_free(p->fftw);
+ free(p);
+}
+
+
+const char *asdf_probe(UnitCell *cell)
+{
+ return "asdf";
+}
+
+#else /* HAVE_FFTW */
+
+int run_asdf(struct image *image, void *ipriv)
+{
+ ERROR("This copy of CrystFEL was compiled without FFTW support.\n");
+ return 0;
+}
+
+
+void *asdf_prepare(IndexingMethod *indm, UnitCell *cell)
+{
+ ERROR("This copy of CrystFEL was compiled without FFTW support.\n");
+ ERROR("To use asdf indexing, recompile with FFTW.\n");
+ return NULL;
+}
+
+
+const char *asdf_probe(UnitCell *cell)
+{
+ return NULL;
+}
+
+
+void asdf_cleanup(void *pp)
+{
+}
+
+#endif /* HAVE_FFTW */
diff --git a/libcrystfel/src/indexers/asdf.h b/libcrystfel/src/indexers/asdf.h
new file mode 100644
index 00000000..b3ccffcc
--- /dev/null
+++ b/libcrystfel/src/indexers/asdf.h
@@ -0,0 +1,56 @@
+/*
+ * asdf.h
+ *
+ * Alexandra's Superior Direction Finder, or
+ * Algorithm Similar to DirAx, FFT-based
+ *
+ * Copyright © 2014-2020 Deutsches Elektronen-Synchrotron DESY,
+ * a research centre of the Helmholtz Association.
+ *
+ * Authors:
+ * 2014-2015 Alexandra Tolstikova <alexandra.tolstikova@desy.de>
+ * 2015,2017 Thomas White <taw@physics.org>
+ *
+ * This file is part of CrystFEL.
+ *
+ * CrystFEL is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * CrystFEL is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with CrystFEL. If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+#ifndef ASDF_H
+#define ASDF_H
+
+#include "index.h"
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+/**
+ * \file asdf.h
+ * The ASDF indexing algorithm.
+ */
+
+extern int run_asdf(struct image *image, void *ipriv);
+
+extern void *asdf_prepare(IndexingMethod *indm, UnitCell *cell);
+extern const char *asdf_probe(UnitCell *cell);
+
+extern void asdf_cleanup(void *pp);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif /* ASDF_H */
diff --git a/libcrystfel/src/indexers/dirax.c b/libcrystfel/src/indexers/dirax.c
new file mode 100644
index 00000000..24be87ba
--- /dev/null
+++ b/libcrystfel/src/indexers/dirax.c
@@ -0,0 +1,683 @@
+/*
+ * dirax.c
+ *
+ * Invoke the DirAx auto-indexing program
+ *
+ * Copyright © 2012-2020 Deutsches Elektronen-Synchrotron DESY,
+ * a research centre of the Helmholtz Association.
+ *
+ * Authors:
+ * 2010-2014,2017 Thomas White <taw@physics.org>
+ *
+ * This file is part of CrystFEL.
+ *
+ * CrystFEL is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * CrystFEL is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with CrystFEL. If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+#include <string.h>
+#include <unistd.h>
+#include <sys/wait.h>
+#include <fcntl.h>
+#include <assert.h>
+#include <sys/ioctl.h>
+#include <errno.h>
+
+#ifdef HAVE_FORKPTY_PTY_H
+#include <pty.h>
+#endif
+#ifdef HAVE_FORKPTY_UTIL_H
+#include <util.h>
+#endif
+
+#include "image.h"
+#include "dirax.h"
+#include "utils.h"
+#include "peaks.h"
+#include "cell-utils.h"
+
+
+/** \file dirax.h */
+
+#define DIRAX_VERBOSE 0
+
+#define MAX_DIRAX_CELL_CANDIDATES (5)
+
+
+typedef enum {
+ DIRAX_INPUT_NONE,
+ DIRAX_INPUT_LINE,
+ DIRAX_INPUT_PROMPT,
+ DIRAX_INPUT_ACL
+} DirAxInputType;
+
+
+struct dirax_private {
+ IndexingMethod indm;
+ UnitCell *template;
+};
+
+
+struct dirax_data {
+
+ /* DirAx auto-indexing low-level stuff */
+ int pty;
+ pid_t pid;
+ char *rbuffer;
+ int rbufpos;
+ int rbuflen;
+
+ /* DirAx auto-indexing high-level stuff */
+ int step;
+ int finished_ok;
+ int read_cell;
+ int best_acl;
+ int best_acl_nh;
+ int acls_tried[MAX_DIRAX_CELL_CANDIDATES];
+ int n_acls_tried;
+ int done;
+ int success;
+
+ float ax;
+ float ay;
+ float az;
+ float bx;
+ float by;
+ float bz;
+ float cx;
+ float cy;
+ float cz;
+
+ struct dirax_private *dp;
+
+};
+
+
+static void dirax_parseline(const char *line, struct image *image,
+ struct dirax_data *dirax)
+{
+ int rf, i, di, acl, acl_nh;
+ float d;
+
+ #if DIRAX_VERBOSE
+ char *copy;
+
+ copy = strdup(line);
+ for ( i=0; i<strlen(copy); i++ ) {
+ if ( copy[i] == '\r' ) copy[i]='r';
+ if ( copy[i] == '\n' ) copy[i]='\0';
+ }
+ STATUS("DirAx: %s\n", copy);
+ free(copy);
+ #endif
+
+ if ( strstr(line, "reflections from file") ) {
+ ERROR("DirAx can't understand this data.\n");
+ return;
+ }
+
+ /* Is this the first line of a unit cell specification? */
+ rf = 0; i = 0;
+ while ( (i<strlen(line)) && ((line[i] == 'R')
+ || (line[i] == 'D') || (line[i] == ' ')) ) {
+ if ( line[i] == 'R' ) rf = 1;
+ if ( (line[i] == 'D') && rf ) {
+ dirax->read_cell = 1;
+ return;
+ }
+ i++;
+ }
+
+ /* Parse unit cell vectors as appropriate */
+ if ( dirax->read_cell == 1 ) {
+
+ /* First row of unit cell values */
+ int r;
+ r = sscanf(line, "%f %f %f %f %f %f",
+ &d, &d, &d, &dirax->ax, &dirax->ay, &dirax->az);
+ if ( r != 6 ) {
+ ERROR("Couldn't understand cell line:\n");
+ ERROR("'%s'\n", line);
+ dirax->read_cell = 0;
+ return;
+ }
+ dirax->ax *= 1e-10;
+ dirax->ay *= 1e-10;
+ dirax->az *= 1e-10;
+ dirax->read_cell++;
+ return;
+
+ } else if ( dirax->read_cell == 2 ) {
+
+ /* Second row of unit cell values */
+ int r;
+ r = sscanf(line, "%f %f %f %f %f %f",
+ &d, &d, &d, &dirax->bx, &dirax->by, &dirax->bz);
+ if ( r != 6 ) {
+ ERROR("Couldn't understand cell line:\n");
+ ERROR("'%s'\n", line);
+ dirax->read_cell = 0;
+ return;
+ }
+ dirax->bx *= 1e-10;
+ dirax->by *= 1e-10;
+ dirax->bz *= 1e-10;
+ dirax->read_cell++;
+ return;
+
+ } else if ( dirax->read_cell == 3 ) {
+
+ /* Third row of unit cell values */
+ int r;
+ UnitCell *cell;
+ Crystal *cr;
+
+ r = sscanf(line, "%f %f %f %f %f %f",
+ &d, &d, &d, &dirax->cx, &dirax->cy, &dirax->cz);
+ if ( r != 6 ) {
+ ERROR("Couldn't understand cell line:\n");
+ ERROR("'%s'\n", line);
+ dirax->read_cell = 0;
+ return;
+ }
+ dirax->cx *= 1e-10;
+ dirax->cy *= 1e-10;
+ dirax->cz *= 1e-10;
+ dirax->read_cell = 0;
+
+ cell = cell_new();
+ cell_set_cartesian(cell, dirax->ax, dirax->ay, dirax->az,
+ dirax->bx, dirax->by, dirax->bz,
+ dirax->cx, dirax->cy, dirax->cz);
+
+ /* Finished reading a cell. */
+
+ cr = crystal_new();
+ if ( cr == NULL ) {
+ ERROR("Failed to allocate crystal.\n");
+ return;
+ }
+ crystal_set_cell(cr, cell);
+ image_add_crystal(image, cr);
+ dirax->done = 1;
+ dirax->success = 1;
+
+ return;
+
+ }
+
+ dirax->read_cell = 0;
+
+ if ( sscanf(line, "%i %i %f %f %f %f %f %f %i", &acl, &acl_nh,
+ &d, &d, &d, &d, &d, &d, &di) == 9 ) {
+ if ( acl_nh > dirax->best_acl_nh ) {
+
+ int i, found = 0;
+
+ for ( i=0; i<dirax->n_acls_tried; i++ ) {
+ if ( dirax->acls_tried[i] == acl ) found = 1;
+ }
+
+ if ( !found ) {
+ dirax->best_acl_nh = acl_nh;
+ dirax->best_acl = acl;
+ }
+
+ }
+ }
+}
+
+
+static void dirax_sendline(const char *line, struct dirax_data *dirax)
+{
+ #if DIRAX_VERBOSE
+ char *copy;
+ int i;
+
+ copy = strdup(line);
+ for ( i=0; i<strlen(copy); i++ ) {
+ if ( copy[i] == '\r' ) copy[i]='\0';
+ if ( copy[i] == '\n' ) copy[i]='\0';
+ }
+ STATUS("To DirAx: '%s'\n", copy);
+ free(copy);
+ #endif
+
+ if ( write(dirax->pty, line, strlen(line)) == -1 ) {
+ ERROR("write() To dirax failed: %s\n", strerror(errno));
+ }
+}
+
+
+static void dirax_send_next(struct image *image, struct dirax_data *dirax)
+{
+ char tmp[32];
+
+ switch ( dirax->step ) {
+
+ case 1 :
+ dirax_sendline("\\echo off\n", dirax);
+ break;
+
+ case 2 :
+ snprintf(tmp, 31, "read xfel.drx\n");
+ dirax_sendline(tmp, dirax);
+ break;
+
+ case 3 :
+ dirax_sendline("dmax 1000\n", dirax);
+ break;
+
+ case 4 :
+ dirax_sendline("indexfit 2\n", dirax);
+ break;
+
+ case 5 :
+ dirax_sendline("levelfit 1000\n", dirax);
+ break;
+
+ case 6 :
+ dirax_sendline("go\n", dirax);
+ dirax->finished_ok = 1;
+ break;
+
+ case 7 :
+ dirax_sendline("acl\n", dirax);
+ break;
+
+ case 8 :
+ if ( dirax->best_acl_nh == 0 ) {
+ /* At this point, DirAx is presenting its ACL prompt
+ * and waiting for a single number. Use an extra
+ * newline to choose automatic ACL selection before
+ * exiting. */
+ dirax_sendline("\nexit\n", dirax);
+ break;
+ }
+ snprintf(tmp, 31, "%i\n", dirax->best_acl);
+ dirax->acls_tried[dirax->n_acls_tried++] = dirax->best_acl;
+ dirax_sendline(tmp, dirax);
+ break;
+
+ case 9 :
+ dirax_sendline("cell\n", dirax);
+ break;
+
+ case 10 :
+ if ( dirax->n_acls_tried == MAX_DIRAX_CELL_CANDIDATES ) {
+ dirax_sendline("exit\n", dirax);
+ } else {
+ /* Go back round for another cell */
+ dirax->best_acl_nh = 0;
+ dirax->step = 7;
+ dirax_sendline("acl\n", dirax);
+ }
+ break;
+
+ default:
+ dirax_sendline("exit\n", dirax);
+ return;
+
+ }
+
+ dirax->step++;
+}
+
+
+static int dirax_readable(struct image *image, struct dirax_data *dirax)
+{
+ int rval;
+ int no_string = 0;
+
+ rval = read(dirax->pty, dirax->rbuffer+dirax->rbufpos,
+ dirax->rbuflen-dirax->rbufpos);
+
+ if ( (rval == -1) || (rval == 0) ) return 1;
+
+ dirax->rbufpos += rval;
+ assert(dirax->rbufpos <= dirax->rbuflen);
+
+ while ( (!no_string) && (dirax->rbufpos > 0) ) {
+
+ int i;
+ int block_ready = 0;
+ DirAxInputType type = DIRAX_INPUT_NONE;
+
+ /* See if there's a full line in the buffer yet */
+ for ( i=0; i<dirax->rbufpos-1; i++ ) {
+ /* Means the last value looked at is rbufpos-2 */
+
+ /* Is there a prompt in the buffer? */
+ if ( (i+7 <= dirax->rbufpos)
+ && (!strncmp(dirax->rbuffer+i, "Dirax> ", 7)) ) {
+ block_ready = 1;
+ type = DIRAX_INPUT_PROMPT;
+ break;
+ }
+
+ /* How about an ACL prompt? */
+ if ( (i+10 <= dirax->rbufpos)
+ && (!strncmp(dirax->rbuffer+i, "acl/auto [", 10)) ) {
+ block_ready = 1;
+ type = DIRAX_INPUT_ACL;
+ break;
+ }
+
+ if ( (dirax->rbuffer[i] == '\r')
+ && (dirax->rbuffer[i+1] == '\n') ) {
+ block_ready = 1;
+ type = DIRAX_INPUT_LINE;
+ break;
+ }
+
+ }
+
+ if ( block_ready ) {
+
+ unsigned int new_rbuflen;
+ unsigned int endbit_length;
+ char *block_buffer = NULL;
+
+ switch ( type ) {
+
+ case DIRAX_INPUT_LINE :
+ /* Make buffer a bit too big to keep Valgrind
+ * quiet about alignment errors */
+ block_buffer = malloc(i+4);
+ memcpy(block_buffer, dirax->rbuffer, i);
+ block_buffer[i] = '\0';
+
+ if ( block_buffer[0] == '\r' ) {
+ memmove(block_buffer, block_buffer+1, i);
+ }
+
+ dirax_parseline(block_buffer, image, dirax);
+ free(block_buffer);
+ endbit_length = i+2;
+ break;
+
+ case DIRAX_INPUT_PROMPT :
+ dirax_send_next(image, dirax);
+ endbit_length = i+7;
+ break;
+
+ case DIRAX_INPUT_ACL :
+ dirax_send_next(image,dirax);
+ endbit_length = i+10;
+ break;
+
+ default :
+ /* Obviously, this never happens :) */
+ ERROR("Unrecognised DirAx input mode! "
+ "I don't know how to understand DirAx\n");
+ return 1;
+
+ }
+
+ /* Now the block's been parsed, it should be
+ * forgotten about */
+ memmove(dirax->rbuffer,
+ dirax->rbuffer + endbit_length,
+ dirax->rbuflen - endbit_length);
+
+ /* Subtract the number of bytes removed */
+ dirax->rbufpos = dirax->rbufpos
+ - endbit_length;
+ new_rbuflen = dirax->rbuflen - endbit_length;
+ if ( new_rbuflen == 0 ) new_rbuflen = 256;
+ dirax->rbuffer = realloc(dirax->rbuffer,
+ new_rbuflen);
+ dirax->rbuflen = new_rbuflen;
+
+ } else {
+
+ if ( dirax->rbufpos==dirax->rbuflen ) {
+
+ /* More buffer space is needed */
+ dirax->rbuffer = realloc(
+ dirax->rbuffer,
+ dirax->rbuflen + 256);
+ dirax->rbuflen = dirax->rbuflen
+ + 256;
+ /* The new space gets used at the next
+ * read, shortly... */
+
+ }
+ no_string = 1;
+
+ }
+
+ }
+
+ return 0;
+}
+
+
+static void write_drx(struct image *image)
+{
+ FILE *fh;
+ int i;
+ char filename[1024];
+
+ snprintf(filename, 1023, "xfel.drx");
+
+ fh = fopen(filename, "w");
+ if ( !fh ) {
+ ERROR("Couldn't open temporary file '%s'\n", filename);
+ return;
+ }
+ fprintf(fh, "%f\n", 0.5); /* Lie about the wavelength. */
+
+ for ( i=0; i<image_feature_count(image->features); i++ ) {
+
+ struct imagefeature *f;
+
+ f = image_get_feature(image->features, i);
+ if ( f == NULL ) continue;
+
+ fprintf(fh, "%10f %10f %10f %8f\n",
+ f->rx/1e10, f->ry/1e10, f->rz/1e10, 1.0);
+
+ }
+ fclose(fh);
+}
+
+
+int run_dirax(struct image *image, void *ipriv)
+{
+ unsigned int opts;
+ int status;
+ int rval;
+ struct dirax_data *dirax;
+
+ write_drx(image);
+
+ dirax = malloc(sizeof(struct dirax_data));
+ if ( dirax == NULL ) {
+ ERROR("Couldn't allocate memory for DirAx data.\n");
+ return 0;
+ }
+
+ dirax->pid = forkpty(&dirax->pty, NULL, NULL, NULL);
+ if ( dirax->pid == -1 ) {
+ ERROR("Failed to fork for DirAx: %s\n", strerror(errno));
+ return 0;
+ }
+ if ( dirax->pid == 0 ) {
+
+ /* Child process: invoke DirAx */
+ struct termios t;
+
+ /* Turn echo off */
+ tcgetattr(STDIN_FILENO, &t);
+ t.c_lflag &= ~(ECHO | ECHOE | ECHOK | ECHONL);
+ tcsetattr(STDIN_FILENO, TCSANOW, &t);
+
+ execlp("dirax", "dirax", (char *)NULL);
+ ERROR("Failed to invoke DirAx.\n");
+ _exit(0);
+
+ }
+
+ dirax->rbuffer = malloc(256);
+ dirax->rbuflen = 256;
+ dirax->rbufpos = 0;
+
+ /* Set non-blocking */
+ opts = fcntl(dirax->pty, F_GETFL);
+ fcntl(dirax->pty, F_SETFL, opts | O_NONBLOCK);
+
+ dirax->step = 1; /* This starts the "initialisation" procedure */
+ dirax->finished_ok = 0;
+ dirax->read_cell = 0;
+ dirax->n_acls_tried = 0;
+ dirax->best_acl_nh = 0;
+ dirax->done = 0;
+ dirax->success = 0;
+ dirax->dp = (struct dirax_private *)ipriv;
+
+ do {
+
+ fd_set fds;
+ struct timeval tv;
+ int sval;
+
+ FD_ZERO(&fds);
+ FD_SET(dirax->pty, &fds);
+
+ tv.tv_sec = 30;
+ tv.tv_usec = 0;
+
+ sval = select(dirax->pty+1, &fds, NULL, NULL, &tv);
+
+ if ( sval == -1 ) {
+
+ const int err = errno;
+
+ switch ( err ) {
+
+ case EINTR:
+ STATUS("Restarting select()\n");
+ break;
+
+ default:
+ ERROR("select() failed: %s\n", strerror(err));
+ rval = 1;
+
+ }
+
+ } else if ( sval != 0 ) {
+ rval = dirax_readable(image, dirax);
+ } else {
+ ERROR("No response from DirAx..\n");
+ rval = 1;
+ }
+
+ } while ( !rval && !dirax->success );
+
+ close(dirax->pty);
+ free(dirax->rbuffer);
+ waitpid(dirax->pid, &status, 0);
+
+ if ( dirax->finished_ok == 0 ) {
+ ERROR("DirAx doesn't seem to be working properly.\n");
+ }
+
+ rval = dirax->success;
+ free(dirax);
+ return rval;
+}
+
+
+void *dirax_prepare(IndexingMethod *indm, UnitCell *cell)
+{
+ struct dirax_private *dp;
+
+ if ( dirax_probe(cell) == NULL ) {
+ ERROR("DirAx does not appear to run properly.\n");
+ ERROR("Please check your DirAx installation.\n");
+ return NULL;
+ }
+
+ /* Flags that DirAx knows about */
+ *indm &= INDEXING_METHOD_MASK;
+
+ dp = malloc(sizeof(struct dirax_private));
+ if ( dp == NULL ) return NULL;
+
+ dp->template = cell;
+ dp->indm = *indm;
+
+ return (IndexingPrivate *)dp;
+}
+
+
+void dirax_cleanup(void *pp)
+{
+ struct dirax_private *p;
+ p = (struct dirax_private *)pp;
+ free(p);
+}
+
+
+const char *dirax_probe(UnitCell *cell)
+{
+ pid_t pid;
+ int pty;
+ int status;
+ FILE *fh;
+ char line[1024];
+ int ok = 0;
+
+ pid = forkpty(&pty, NULL, NULL, NULL);
+ if ( pid == -1 ) {
+ return NULL;
+ }
+ if ( pid == 0 ) {
+
+ /* Child process: invoke DirAx */
+ struct termios t;
+
+ /* Turn echo off */
+ tcgetattr(STDIN_FILENO, &t);
+ t.c_lflag &= ~(ECHO | ECHOE | ECHOK | ECHONL);
+ tcsetattr(STDIN_FILENO, TCSANOW, &t);
+
+ execlp("dirax", "dirax", (char *)NULL);
+ _exit(1);
+
+ }
+
+ fh = fdopen(pty, "r");
+ if ( fgets(line, 1024, fh) != NULL ) {
+ if ( strncmp(line, "dirax", 5) == 0 ) {
+ ok = 1;
+ }
+ }
+
+ fclose(fh);
+ close(pty);
+ waitpid(pid, &status, 0);
+
+ if ( ok ) return "dirax";
+ return NULL;
+}
diff --git a/libcrystfel/src/indexers/dirax.h b/libcrystfel/src/indexers/dirax.h
new file mode 100644
index 00000000..da4ae3d5
--- /dev/null
+++ b/libcrystfel/src/indexers/dirax.h
@@ -0,0 +1,52 @@
+/*
+ * dirax.h
+ *
+ * Invoke the DirAx auto-indexing program
+ *
+ * Copyright © 2012-2020 Deutsches Elektronen-Synchrotron DESY,
+ * a research centre of the Helmholtz Association.
+ *
+ * Authors:
+ * 2010,2012-2014,2017 Thomas White <taw@physics.org>
+ *
+ * This file is part of CrystFEL.
+ *
+ * CrystFEL is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * CrystFEL is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with CrystFEL. If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+#ifndef DIRAX_H
+#define DIRAX_H
+
+#include "index.h"
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+/** \file dirax.h
+ * DirAx indexer interface
+ */
+extern int run_dirax(struct image *image, void *ipriv);
+
+extern void *dirax_prepare(IndexingMethod *indm, UnitCell *cell);
+extern const char *dirax_probe(UnitCell *cell);
+
+extern void dirax_cleanup(void *pp);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif /* DIRAX_H */
diff --git a/libcrystfel/src/indexers/felix.c b/libcrystfel/src/indexers/felix.c
new file mode 100644
index 00000000..524e68a7
--- /dev/null
+++ b/libcrystfel/src/indexers/felix.c
@@ -0,0 +1,963 @@
+/*
+ * felix.c
+ *
+ * Invoke Felix for multi-crystal autoindexing.
+ *
+ * Copyright © 2015-2020 Deutsches Elektronen-Synchrotron DESY,
+ * a research centre of the Helmholtz Association.
+ *
+ * Authors:
+ * 2015-2018 Thomas White <taw@physics.org>
+ * 2015 Kenneth Beyerlein <kenneth.beyerlein@desy.de>
+ *
+ * This file is part of CrystFEL.
+ *
+ * CrystFEL is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * CrystFEL is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with CrystFEL. If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+
+/** \file felix.h */
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+#include <string.h>
+#include <sys/types.h>
+#include <sys/wait.h>
+#include <sys/stat.h>
+#include <unistd.h>
+#include <assert.h>
+#include <fcntl.h>
+#include <errno.h>
+
+#ifdef HAVE_CLOCK_GETTIME
+#include <time.h>
+#else
+#include <sys/time.h>
+#endif
+
+#ifdef HAVE_FORKPTY_PTY_H
+#include <pty.h>
+#endif
+#ifdef HAVE_FORKPTY_UTIL_H
+#include <util.h>
+#endif
+
+#include "image.h"
+#include "utils.h"
+#include "peaks.h"
+#include "cell.h"
+#include "cell-utils.h"
+#include "felix.h"
+
+
+#define FELIX_VERBOSE 0
+
+
+struct felix_options
+{
+ double ttmin; /* radians */
+ double ttmax; /* radians */
+ int min_visits;
+ double min_completeness;
+ double max_uniqueness;
+ int n_voxels;
+ double fraction_max_visits;
+ double sigma;
+ double domega;
+ double max_internal_angle;
+};
+
+
+/* Global private data, prepared once */
+struct felix_private
+{
+ IndexingMethod indm;
+ UnitCell *cell;
+
+ /* Options specific to Felix */
+ int spacegroup;
+ float tthrange_min;
+ float tthrange_max;
+ float etarange_min;
+ float etarange_max;
+ float domega;
+ float omegarange_min;
+ float omegarange_max;
+ int min_visits; /* related to "cuts" */
+ float min_completeness; /* related to "cuts" */
+ float max_uniqueness; /* related to "cuts" */
+ int n_voxels; /* related to "frustsumsize" */
+ float fraction_max_visits; /* related to "frustsumsize" */
+ float sigma_tth; /* related to "uncertainties" */
+ float sigma_eta; /* related to "uncertainties" */
+ float sigma_omega; /* related to "uncertainties" */
+ int n_sigmas;
+ int force4frustums;
+ float max_internal_angle;
+
+ /*Felix v0.3 options*/
+ int orispace_frustum;
+ int orispace_octa;
+ char *readhkl_file;
+ float maxtime;
+
+};
+
+
+/* Data needed to call Felix */
+struct felix_data {
+
+ struct felix_private *gp;
+
+ /* Low-level stuff */
+ int pty;
+ pid_t pid;
+ char *rbuffer;
+ int rbufpos;
+ int rbuflen;
+
+};
+
+
+static int read_felix(struct felix_private *gp, struct image *image,
+ char *filename)
+{
+ FILE *fh;
+ int d1;
+ float d2;
+ float ubi11, ubi12, ubi13;
+ float ubi21, ubi22, ubi23;
+ float ubi31, ubi32, ubi33;
+ float mean_ia;
+ int ngv;
+ char line[1024];
+ int r;
+ int n_crystals = 0;
+
+ fh = fopen(filename, "r");
+ if ( fh == NULL ) {
+ ERROR("Can't open '%s'\n", filename);
+ return 0;
+ }
+
+ /* Read and discard first line */
+ if ( fgets( line, 1024, fh ) == NULL ) {
+ ERROR("Failed to read *.felix file.\n");
+ return 0;
+ }
+
+ do {
+
+ Crystal *cr;
+ UnitCell *cell;
+
+ /* One line per grain */
+ if ( fgets( line, 1024, fh ) == NULL ) {
+ break;
+ }
+ /* File format of the .felix files
+ * version 0.1 - 0.2
+ *
+ * r = sscanf(line, "%i %f %i %i %f %f %f %f %f %f %f %f %f"
+ * "%f %f %f %f %f %f %f %f %f %f %f %f",
+ * &d1, &mean_ia, &ngv, &ngv_unique, &d2, &d2, &d2,
+ * &d2, &d2, &d2, &d2, &d2, &d2, &d2, &d2, &d2,
+ * &ubi11, &ubi12, &ubi13,
+ * &ubi21, &ubi22, &ubi23,
+ * &ubi31, &ubi32, &ubi33);
+ *
+ * if ( r != 25 ) {
+ * ERROR("Only %i parameters in .felix file\n", r);
+ * return 1;
+ * }
+ */
+
+ /* version 0.3 - present */
+ r = sscanf(line, "%i %f %i %f %f %f %f %f %f"
+ "%f %f %f %f %f %f %f %f %f %f %f %f",
+ &d1, &mean_ia, &ngv, &d2, &d2,
+ &d2, &d2, &d2, &d2, &d2, &d2, &d2,
+ &ubi11, &ubi12, &ubi13,
+ &ubi21, &ubi22, &ubi23,
+ &ubi31, &ubi32, &ubi33);
+
+ if ( r != 21 ) {
+ ERROR("Only %i parameters in .felix file, "
+ "check version and format.\n", r);
+ return -1;
+ }
+
+ cell = cell_new();
+
+ cell_set_cartesian(cell, ubi12/1e10, ubi13/1e10, ubi11/1e10,
+ ubi22/1e10, ubi23/1e10, ubi21/1e10,
+ ubi32/1e10, ubi33/1e10, ubi31/1e10);
+ cell_set_lattice_type(cell, cell_get_lattice_type(gp->cell));
+ cell_set_centering(cell, cell_get_centering(gp->cell));
+ cell_set_unique_axis(cell, cell_get_unique_axis(gp->cell));
+
+ cr = crystal_new();
+ if ( cr == NULL ) {
+ ERROR( "Failed to allocate crystal.\n" );
+ return 0;
+ }
+
+ crystal_set_cell(cr, cell);
+
+ /* Poor indexing criterion for Felix v0.1
+ *
+ * if (mean_ia > MAX_MEAN_IA || ngv < MIN_NGV ||
+ * ngv_unique < MIN_NGV_UNIQUE ){
+ * crystal_set_user_flag(cr, 1);
+ */
+
+ /* Poor indexing criterion for Felix v0.3 */
+
+ if ( mean_ia > gp->max_internal_angle ){
+ crystal_set_user_flag(cr, 1);
+ }
+
+ /* All crystals are saved to the image,
+ * but only good ones will be later written to the stream file.
+ */
+
+ image_add_crystal(image, cr);
+ n_crystals++;
+
+ } while ( !feof(fh) );
+
+ fclose(fh);
+
+ return n_crystals;
+}
+
+
+static void gs_parseline(char *line, struct image *image,
+ struct felix_data *gs)
+{
+ #if FELIX_VERBOSE
+ STATUS("%s\n", line);
+ #endif
+}
+
+
+static int felix_readable(struct image *image, struct felix_data *gs)
+{
+ int rval;
+ int no_string = 0;
+
+ rval = read(gs->pty, gs->rbuffer+gs->rbufpos, gs->rbuflen-gs->rbufpos);
+
+ if ( (rval == -1) || (rval == 0) ) return 1;
+
+ gs->rbufpos += rval;
+ assert(gs->rbufpos <= gs->rbuflen);
+
+ while ( (!no_string) && (gs->rbufpos > 0) ) {
+
+ int i;
+ int block_ready = 0;
+
+ /* See if there's a full line in the buffer yet */
+ for ( i=0; i<gs->rbufpos-1; i++ ) {
+ /* Means the last value looked at is rbufpos-2 */
+
+ if ( (gs->rbuffer[i] == '\r')
+ && ( gs->rbuffer[i+1] == '\n' ) ) {
+ block_ready = 1;
+ break;
+ }
+
+ }
+
+ if ( block_ready ) {
+
+ unsigned int new_rbuflen;
+ unsigned int endbit_length;
+ char *block_buffer = NULL;
+
+ block_buffer = malloc(i+1);
+ memcpy(block_buffer, gs->rbuffer, i);
+ block_buffer[i] = '\0';
+
+ if ( block_buffer[0] == '\r' ) {
+ memmove(block_buffer, block_buffer+1, i);
+ }
+
+ gs_parseline(block_buffer, image, gs);
+ free(block_buffer);
+ endbit_length = i+2;
+
+ /* Now the block's been parsed, it should be
+ * forgotten about */
+ memmove(gs->rbuffer,
+ gs->rbuffer + endbit_length,
+ gs->rbuflen - endbit_length);
+
+ /* Subtract the number of bytes removed */
+ gs->rbufpos = gs->rbufpos - endbit_length;
+ new_rbuflen = gs->rbuflen - endbit_length;
+ if ( new_rbuflen == 0 ) new_rbuflen = 256;
+ gs->rbuffer = realloc(gs->rbuffer, new_rbuflen);
+ gs->rbuflen = new_rbuflen;
+
+ } else {
+
+ if ( gs->rbufpos == gs->rbuflen ) {
+
+ /* More buffer space is needed */
+ gs->rbuffer = realloc(gs->rbuffer,
+ gs->rbuflen + 256);
+ gs->rbuflen = gs->rbuflen + 256;
+ /* The new space gets used at the next
+ * read, shortly... */
+
+ }
+ no_string = 1;
+
+ }
+
+ }
+
+ return 0;
+}
+
+
+static void write_gve(struct image *image, struct felix_private *gp)
+{
+ FILE *fh;
+ int i;
+ char filename[1024];
+ double a, b, c, al, be, ga;
+ snprintf(filename, 1023, "xfel.gve");
+ fh = fopen(filename, "w");
+ if ( !fh ) {
+ ERROR("Couldn't open temporary file '%s'\n", filename);
+ return;
+ }
+
+ cell_get_parameters(gp->cell, &a, &b, &c, &al, &be, &ga);
+ fprintf(fh, "%.6f %.6f %.6f %.6f %.6f %.6f P\n", a*1e10, b*1e10, c*1e10,
+ rad2deg(al), rad2deg(be), rad2deg(ga));
+ fprintf(fh, "# wavelength = %.6f\n", image->lambda*1e10);
+ fprintf(fh, "# wedge = 0.000000\n");
+ fprintf(fh, "# ds h k l\n");
+ fprintf(fh, "# xr yr zr xc yc ds eta omega\n");
+
+ for ( i=0; i<image_feature_count(image->features); i++ ) {
+
+ struct imagefeature *f;
+
+ f = image_get_feature(image->features, i);
+ if ( f == NULL ) continue;
+
+ fprintf(fh, "%.6f %.6f %.6f 0 0 %.6f %.6f %.6f 0\n",
+ f->rz/1e10, f->rx/1e10, f->ry/1e10,
+ modulus(f->rx, f->ry, f->rz)/1e10, /* dstar */
+ rad2deg(atan2(f->ry, f->rx)), 0.0); /* eta, omega */
+
+ }
+ fclose(fh);
+}
+
+
+static char *write_ini(struct image *image, struct felix_private *gp)
+{
+ FILE *fh;
+ char *filename;
+ char gveFilename[1024];
+ char logFilename[1024];
+
+ filename = malloc(1024);
+ if ( filename == NULL ) return NULL;
+
+ snprintf(filename, 1023, "xfel.ini");
+ snprintf(gveFilename, 1023, "xfel.gve");
+ snprintf(logFilename, 1023, "xfel.log");
+
+ fh = fopen(filename, "w");
+ if ( !fh ) {
+ ERROR("Couldn't open temporary file '%s'\n", filename);
+ free(filename);
+ return NULL;
+ }
+
+ fprintf(fh, "spacegroup %i\n", gp->spacegroup);
+ fprintf(fh, "tthrange %f %f\n", rad2deg(gp->tthrange_min),
+ rad2deg(gp->tthrange_max));
+ fprintf(fh, "etarange %f %f\n", gp->etarange_min, gp->etarange_max);
+ fprintf(fh, "domega %f\n", gp->domega);
+ fprintf(fh, "omegarange %f %f\n", gp->omegarange_min, gp->omegarange_max);
+ fprintf(fh, "filespecs %s %s\n", gveFilename, logFilename);
+ fprintf(fh, "cuts %i %f %f\n", gp->min_visits, gp->min_completeness,
+ gp->max_uniqueness);
+ fprintf(fh, "frustumsize %i %f\n", gp->n_voxels,
+ gp->fraction_max_visits);
+ fprintf(fh, "uncertainties %f %f %f\n", gp->sigma_tth,
+ gp->sigma_eta, gp->sigma_omega);
+ fprintf(fh, "nsigmas %i\n", gp->n_sigmas);
+
+ if ( gp->force4frustums == 1 ){
+ fprintf(fh, "force4frustums\n");
+ }
+
+ if ( gp->orispace_frustum == 1 ){
+ fprintf(fh, "orispace frustum\n");
+ } else if ( gp->orispace_octa ==1 ){
+ fprintf(fh, "orispace octa\n");
+ } else{
+ ERROR("No felix supported orispace specified.\n");
+ free(filename);
+ filename = NULL;
+ }
+
+ /* If an hkl file is not specified, generate the peak list. */
+ if ( gp->readhkl_file != NULL ){
+ fprintf(fh, "readhkl %s\n", gp->readhkl_file);
+ } else{
+ fprintf(fh, "genhkl\n");
+ }
+
+ fprintf(fh, "maxtime %f\n", gp->maxtime);
+
+ fclose(fh);
+
+ return filename;
+}
+
+
+int felix_index(struct image *image, IndexingPrivate *ipriv)
+{
+ unsigned int opts;
+ int status;
+ int rval;
+ struct felix_data *felix;
+ struct felix_private *gp = (struct felix_private *) ipriv;
+ char *ini_filename;
+ char gff_filename[1024];
+
+ write_gve (image, gp);
+ ini_filename = write_ini (image, gp);
+
+ if ( ini_filename == NULL ) {
+ ERROR("Failed to write ini file for Felix.\n");
+ return 0;
+ }
+
+ felix = malloc(sizeof(struct felix_data));
+ if ( felix == NULL ) {
+ ERROR("Couldn't allocate memory for Felix data.\n");
+ return 0;
+ }
+
+ felix->gp = gp;
+
+ snprintf(gff_filename, 1023, "xfel.felix");
+ remove(gff_filename);
+
+ felix->pid = forkpty(&felix->pty, NULL, NULL, NULL);
+ if ( felix->pid == -1 ) {
+ ERROR("Failed to fork for Felix: %s\n", strerror(errno));
+ return 0;
+ }
+ if ( felix->pid == 0 ) {
+
+ /* Child process: invoke Felix */
+ struct termios t;
+
+ /* Turn echo off */
+ tcgetattr(STDIN_FILENO, &t);
+ t.c_lflag &= ~(ECHO | ECHOE | ECHOK | ECHONL);
+ tcsetattr(STDIN_FILENO, TCSANOW, &t);
+
+ STATUS("Running Felix '%s'\n", ini_filename);
+ execlp("Felix", "Felix", ini_filename, (char *)NULL);
+ ERROR("Failed to invoke Felix.\n");
+ _exit(0);
+
+ }
+
+ free(ini_filename);
+
+ felix->rbuffer = malloc(256);
+ felix->rbuflen = 256;
+ felix->rbufpos = 0;
+
+ /* Set non-blocking */
+ opts = fcntl(felix->pty, F_GETFL);
+ fcntl(felix->pty, F_SETFL, opts | O_NONBLOCK);
+
+ do {
+
+ fd_set fds;
+ struct timeval tv;
+ int sval;
+
+ FD_ZERO(&fds);
+ FD_SET(felix->pty, &fds);
+
+ tv.tv_sec = gp->maxtime + 1.0;
+ tv.tv_usec = 0;
+
+ sval = select(felix->pty+1, &fds, NULL, NULL, &tv);
+
+ if ( sval == -1 ) {
+
+ const int err = errno;
+
+ switch ( err ) {
+
+ case EINTR:
+ STATUS("Restarting select()\n");
+ rval = 0;
+ break;
+
+ default:
+ ERROR("select() failed: %s\n", strerror(err));
+ rval = 1;
+
+ }
+
+ } else if ( sval != 0 ) {
+ rval = felix_readable(image, felix);
+ } else {
+ ERROR("No response from Felix..\n");
+ rval = 1;
+ }
+
+ } while ( !rval );
+
+ close(felix->pty);
+ free(felix->rbuffer);
+ waitpid(felix->pid, &status, 0);
+
+ if ( status != 0 ) {
+ ERROR("Felix either timed out, or is not working properly.\n");
+ free(felix);
+ return 0;
+ }
+
+ rval = read_felix(gp, image, gff_filename);
+
+ free(felix);
+ return rval;
+
+}
+
+
+static int sg_number_for_cell(UnitCell *cell)
+{
+ LatticeType lattice = cell_get_lattice_type(cell);
+ char cen = cell_get_centering(cell);
+
+ switch (lattice)
+ {
+ case L_TRICLINIC:
+ return 1; /* P1 */
+
+ case L_MONOCLINIC:
+ switch ( cen ) {
+ case 'P' : return 3; /* P2 */
+ case 'C' : return 5; /* C2 */
+ default : return 0;
+ }
+
+ case L_ORTHORHOMBIC:
+ switch ( cen ) {
+ case 'P' : return 16; /* P222 */
+ case 'C' : return 21; /* C222 */
+ case 'F' : return 22; /* F222 */
+ case 'I' : return 23; /* I222 */
+ case 'A' : return 38; /* Amm2 */
+ default : return 0;
+ }
+
+ case L_TETRAGONAL:
+ switch ( cen ) {
+ case 'P' : return 89; /* P422 */
+ case 'I' : return 97; /* I422 */
+ default : return 0;
+ }
+
+ case L_RHOMBOHEDRAL:
+ return 155; /* R32 */
+
+ case L_HEXAGONAL:
+ switch ( cen ) {
+ case 'P' : return 177; /* P622 */
+ case 'H' : return 143; /* P3 */
+ default : return 0;
+ }
+
+ case L_CUBIC:
+ switch ( cen ) {
+ case 'P' : return 207; /* P432 */
+ case 'F' : return 209; /* F432 */
+ case 'I' : return 211; /* I432 */
+ default : return 0;
+ }
+
+ default:
+ return 0;
+ }
+}
+
+
+void *felix_prepare(IndexingMethod *indm, UnitCell *cell,
+ struct felix_options *opts)
+{
+ struct felix_private *gp;
+
+ if ( !cell_has_parameters(cell) ) {
+ ERROR("Felix needs a unit cell.\n");
+ return NULL;
+ }
+
+ if ( felix_probe(cell) == NULL ) {
+ ERROR("Felix does not appear to run properly.\n");
+ ERROR("Please check your Felix installation.\n");
+ return NULL;
+ }
+
+ gp = calloc(1, sizeof(*gp));
+ if ( gp == NULL ) return NULL;
+
+ /* Flags that Felix knows about */
+ *indm &= INDEXING_METHOD_MASK
+ | INDEXING_USE_LATTICE_TYPE | INDEXING_USE_CELL_PARAMETERS;
+
+ gp->cell = cell;
+ gp->indm = *indm;
+
+ /* Default values of felix options */
+ gp->spacegroup = sg_number_for_cell(cell);
+ if ( gp->spacegroup == 0 ) {
+ ERROR("Couldn't determine representative space group for your cell.\n");
+ ERROR("Try again with a more conventional cell.\n");
+ return NULL;
+ }
+
+ /* Default parameters */
+ gp->n_voxels = 100;
+ gp->etarange_min = 0;
+ gp->etarange_max = 360;
+ gp->domega = 2;
+ gp->omegarange_min = -1.0;
+ gp->omegarange_max = 1.0;
+ gp->min_visits = 15;
+ gp->min_completeness = 0.001;
+ gp->max_uniqueness = 0.5;
+ gp->fraction_max_visits = 0.75;
+ gp->sigma_tth = 0.2;
+ gp->sigma_eta = 0.2;
+ gp->sigma_omega = 0.2;
+ gp->n_sigmas = 1;
+ gp->force4frustums = 0;
+ gp->orispace_frustum = 1;
+ gp->orispace_octa = 0;
+ gp->readhkl_file = NULL;
+ gp->maxtime = 120.0;
+ gp->tthrange_min = deg2rad(0.0);
+ gp->tthrange_max = deg2rad(30.0);
+ gp->max_internal_angle = 0.25;
+
+ if ( opts->ttmin > 0.0 ) {
+ gp->tthrange_min = opts->ttmin;
+ }
+ if ( opts->ttmax > 0.0 ) {
+ gp->tthrange_max = opts->ttmax;
+ }
+ if ( opts->min_visits > 0 ) {
+ gp->min_visits = opts->min_visits;
+ }
+ if ( opts->min_completeness > 0.0 ) {
+ gp->min_completeness = opts->min_completeness;
+ }
+ if ( opts->max_uniqueness > 0.0 ) {
+ gp->max_uniqueness = opts->max_uniqueness;
+ }
+ if ( opts->n_voxels > 0 ) {
+ gp->n_voxels = opts->n_voxels;
+ }
+ if ( opts->fraction_max_visits > 0.0 ) {
+ gp->fraction_max_visits = opts->fraction_max_visits;
+ }
+ if ( opts->sigma > 0.0 ) {
+ gp->sigma_tth = opts->sigma;
+ gp->sigma_eta = opts->sigma;
+ gp->sigma_omega = opts->sigma;
+ }
+ if ( opts->domega > 0.0 ) {
+ gp->domega = opts->domega;
+ }
+ if ( opts->max_internal_angle > 0 ) {
+ gp->max_internal_angle = opts->max_internal_angle;
+ }
+
+ return (IndexingPrivate *)gp;
+}
+
+
+void felix_cleanup(IndexingPrivate *pp)
+{
+ struct felix_private *p;
+
+ p = (struct felix_private *) pp;
+ free(p->readhkl_file);
+ free(p);
+}
+
+
+static int file_exists(const char *filename)
+{
+ struct stat s;
+
+ if ( stat(filename, &s) != 0 ) {
+ if ( errno == ENOENT ) return 0;
+ ERROR("Failed to check for %s.\n", filename);
+ exit(1);
+ }
+
+ return 1;
+}
+
+
+const char *felix_probe(UnitCell *cell)
+{
+ pid_t pid;
+ int pty;
+ int status;
+ FILE *fh;
+ char line[1024];
+ int ok = 0;
+
+ if ( !cell_has_parameters(cell) ) {
+ return NULL;
+ }
+
+ /* Felix will write gmon.out when we test it, which we are
+ * are going to delete afterwards. Better check the file doesn't exist
+ * first, in case it was important. */
+ if ( file_exists("gmon.out") ) {
+ ERROR("Please move or delete gmon.out from the working "
+ "directory first.\n");
+ exit(1);
+ }
+
+ pid = forkpty(&pty, NULL, NULL, NULL);
+ if ( pid == -1 ) {
+ return NULL;
+ }
+ if ( pid == 0 ) {
+
+ /* Child process: invoke DirAx */
+ struct termios t;
+
+ /* Turn echo off */
+ tcgetattr(STDIN_FILENO, &t);
+ t.c_lflag &= ~(ECHO | ECHOE | ECHOK | ECHONL);
+ tcsetattr(STDIN_FILENO, TCSANOW, &t);
+
+ execlp("Felix", "Felix", (char *)NULL);
+ _exit(1);
+
+ }
+
+ fh = fdopen(pty, "r");
+ if ( fgets(line, 1024, fh) == NULL ) {
+ ok = 0;
+ } else {
+ if ( strncmp(line, "Felix", 5) == 0 ) {
+ ok = 1;
+ }
+ }
+
+ fclose(fh);
+ close(pty);
+ waitpid(pid, &status, 0);
+
+ unlink("gmon.out");
+
+ if ( ok ) return "felix";
+ return NULL;
+}
+
+
+static void felix_show_help()
+{
+ printf("Parameters for the Felix indexing algorithm:\n"
+" --felix-domega Degree range of omega (moscaicity) to consider.\n"
+" Default: 2\n"
+" --felix-fraction-max-visits\n"
+" Cutoff for minimum fraction of the max visits.\n"
+" Default: 0.75\n"
+" --felix-max-internal-angle\n"
+" Cutoff for maximum internal angle between observed\n"
+" spots and predicted spots. Default: 0.25\n"
+" --felix-max-uniqueness\n"
+" Cutoff for maximum fraction of found spots which\n"
+" can belong to other crystallites. Default: 0.5\n"
+" --felix-min-completeness\n"
+" Cutoff for minimum fraction of projected spots\n"
+" found in the pattern. Default: 0.001\n"
+" --felix-min-visits\n"
+" Cutoff for minimum number of voxel visits.\n"
+" Default: 15\n"
+" --felix-num-voxels Number of voxels for Rodrigues space search\n"
+" Default: 100\n"
+" --felix-sigma The sigma of the 2theta, eta and omega angles.\n"
+" Default: 0.2\n"
+" --felix-tthrange-max Maximum 2theta to consider for indexing (degrees)\n"
+" Default: 30\n"
+" --felix-tthrange-min Minimum 2theta to consider for indexing (degrees)\n"
+" Default: 0\n"
+);
+}
+
+
+static error_t felix_parse_arg(int key, char *arg,
+ struct argp_state *state)
+{
+ struct felix_options **opts_ptr = state->input;
+ float tmp;
+
+ switch ( key ) {
+
+ case ARGP_KEY_INIT :
+ *opts_ptr = malloc(sizeof(struct felix_options));
+ if ( *opts_ptr == NULL ) return ENOMEM;
+ (*opts_ptr)->ttmin = -1.0;
+ (*opts_ptr)->ttmax = -1.0;
+ (*opts_ptr)->min_visits = 0;
+ (*opts_ptr)->min_completeness = -1.0;
+ (*opts_ptr)->max_uniqueness = -1.0;
+ (*opts_ptr)->n_voxels = 0;
+ (*opts_ptr)->fraction_max_visits = -1.0;
+ (*opts_ptr)->sigma = -1.0;
+ (*opts_ptr)->domega = -1.0;
+ (*opts_ptr)->max_internal_angle = -1.0;
+ break;
+
+ case 1 :
+ felix_show_help();
+ return EINVAL;
+
+ case 2 :
+ if ( sscanf(arg, "%f", &tmp) != 1 ) {
+ ERROR("Invalid value for --felix-tthrange-min\n");
+ return EINVAL;
+ }
+ (*opts_ptr)->ttmin = deg2rad(tmp);
+ break;
+
+ case 3 :
+ if ( sscanf(arg, "%f", &tmp) != 1 ) {
+ ERROR("Invalid value for --felix-tthrange-max\n");
+ return EINVAL;
+ }
+ (*opts_ptr)->ttmax = deg2rad(tmp);
+ break;
+
+ case 4 :
+ if ( sscanf(arg, "%d", &(*opts_ptr)->min_visits) != 1 ) {
+ ERROR("Invalid value for --felix-min-visits\n");
+ return EINVAL;
+ }
+ break;
+
+ case 5 :
+ if ( sscanf(arg, "%lf", &(*opts_ptr)->min_completeness) != 1 ) {
+ ERROR("Invalid value for --felix-min-completeness\n");
+ return EINVAL;
+ }
+ break;
+
+ case 6 :
+ if ( sscanf(arg, "%lf", &(*opts_ptr)->max_uniqueness) != 1 ) {
+ ERROR("Invalid value for --felix-max-uniqueness\n");
+ return EINVAL;
+ }
+ break;
+
+ case 7 :
+ if ( sscanf(arg, "%d", &(*opts_ptr)->n_voxels) != 1 ) {
+ ERROR("Invalid value for --felix-num-voxels\n");
+ return EINVAL;
+ }
+ break;
+
+ case 8 :
+ if ( sscanf(arg, "%lf", &(*opts_ptr)->fraction_max_visits) != 1 ) {
+ ERROR("Invalid value for --felix-fraction-max-visits\n");
+ return EINVAL;
+ }
+ break;
+
+ case 9 :
+ if ( sscanf(arg, "%lf", &(*opts_ptr)->sigma) != 1 ) {
+ ERROR("Invalid value for --felix-sigma\n");
+ return EINVAL;
+ }
+ break;
+
+ case 10 :
+ if ( sscanf(arg, "%lf", &(*opts_ptr)->domega) != 1 ) {
+ ERROR("Invalid value for --felix-domega\n");
+ return EINVAL;
+ }
+ break;
+
+ case 11 :
+ if ( sscanf(arg, "%lf", &(*opts_ptr)->max_internal_angle) != 1 ) {
+ ERROR("Invalid value for --felix-max-internal-angle\n");
+ return EINVAL;
+ }
+ break;
+
+ default :
+ return ARGP_ERR_UNKNOWN;
+
+ }
+
+ return 0;
+}
+
+
+static struct argp_option felix_options[] = {
+
+ {"help-felix", 1, NULL, OPTION_NO_USAGE,
+ "Show options for Felix indexing algorithm", 99},
+ {"felix-tthrange-min", 2, "2theta", OPTION_HIDDEN, NULL},
+ {"felix-tthrange-max", 3, "2theta", OPTION_HIDDEN, NULL},
+ {"felix-min-visits", 4, "n", OPTION_HIDDEN, NULL},
+ {"felix-min-completeness", 5, "frac", OPTION_HIDDEN, NULL},
+ {"felix-max-uniqueness", 6, "n", OPTION_HIDDEN, NULL},
+ {"felix-num-voxels", 7, "n", OPTION_HIDDEN, NULL},
+ {"felix-fraction-max-visits", 8, "n", OPTION_HIDDEN, NULL},
+ {"felix-sigma", 9, "n", OPTION_HIDDEN, NULL},
+ {"felix-domega", 10, "n", OPTION_HIDDEN, NULL},
+ {"felix-max-internal-angle", 11, "ang", OPTION_HIDDEN, NULL},
+
+ {0}
+};
+
+
+struct argp felix_argp = { felix_options, felix_parse_arg,
+ NULL, NULL, NULL, NULL, NULL };
diff --git a/libcrystfel/src/indexers/felix.h b/libcrystfel/src/indexers/felix.h
new file mode 100644
index 00000000..3c9d4a94
--- /dev/null
+++ b/libcrystfel/src/indexers/felix.h
@@ -0,0 +1,52 @@
+/*
+ * felix.h
+ *
+ * Invoke Felix for multi-crystal autoindexing
+ *
+ * Copyright © 2013-2020 Deutsches Elektronen-Synchrotron DESY,
+ * a research centre of the Helmholtz Association.
+ *
+ * Authors:
+ * 2010-2013,2017 Thomas White <taw@physics.org>
+ * 2013-2014 Kenneth Beyerlein <kenneth.beyerlein@desy.de>
+ *
+ * This file is part of CrystFEL.
+ *
+ * CrystFEL is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * CrystFEL is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with CrystFEL. If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+#ifndef FELIX_H
+#define FELIX_H
+
+#include <argp.h>
+
+#include "cell.h"
+
+/**
+ * \file felix.h
+ * Felix indexer interface
+ */
+
+extern void *felix_prepare(IndexingMethod *indm, UnitCell *cell,
+ struct felix_options *opts);
+
+extern const char *felix_probe(UnitCell *cell);
+
+extern void felix_cleanup(IndexingPrivate *pp);
+
+extern int felix_index(struct image *image, IndexingPrivate *p);
+
+
+#endif /* FELIX_H */
diff --git a/libcrystfel/src/indexers/mosflm.c b/libcrystfel/src/indexers/mosflm.c
new file mode 100644
index 00000000..bacd345f
--- /dev/null
+++ b/libcrystfel/src/indexers/mosflm.c
@@ -0,0 +1,945 @@
+/*
+ * mosflm.c
+ *
+ * Invoke the DPS auto-indexing algorithm through MOSFLM
+ *
+ * Copyright © 2012-2020 Deutsches Elektronen-Synchrotron DESY,
+ * a research centre of the Helmholtz Association.
+ * Copyright © 2012 Richard Kirian
+ *
+ * Authors:
+ * 2010 Richard Kirian <rkirian@asu.edu>
+ * 2010-2018 Thomas White <taw@physics.org>
+ * 2014 Takanori Nakane <nakane.t@gmail.com>
+ *
+ * This file is part of CrystFEL.
+ *
+ * CrystFEL is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * CrystFEL is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with CrystFEL. If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+/* TODO:
+ *
+ * Properly read the newmat file (don't use fscanf-- spaces between numers
+ * are not guaranteed)
+ *
+ * "Success" is indicated by existence of NEWMAT file written by mosflm.
+ * Better to interact with mosflm directly in order to somehow verify success.
+ *
+ * Investigate how these keywords affect mosflms behavior:
+ *
+ * MOSAICITY
+ * DISPERSION
+ * DIVERGENCE
+ * POLARISATION
+ * POSTREF BEAM
+ * POSTREF USEBEAM OFF
+ * PREREFINE ON
+ * EXTRA ON
+ * POSTREF ON
+ *
+ * These did not seem to affect the results by my (Rick's) experience, probably
+ * because they are only used conjunction with image intensity data, but it's
+ * worth another look at the documentation.
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+#include <string.h>
+#include <sys/types.h>
+#include <sys/wait.h>
+#include <sys/stat.h>
+#include <unistd.h>
+#include <assert.h>
+#include <fcntl.h>
+#include <errno.h>
+
+#ifdef HAVE_FORKPTY_PTY_H
+#include <pty.h>
+#endif
+#ifdef HAVE_FORKPTY_UTIL_H
+#include <util.h>
+#endif
+
+#ifdef HAVE_CLOCK_GETTIME
+#include <time.h>
+#else
+#include <sys/time.h>
+#endif
+
+#include "image.h"
+#include "mosflm.h"
+#include "utils.h"
+#include "peaks.h"
+#include "cell-utils.h"
+
+/** \file mosflm.h */
+
+#define MOSFLM_VERBOSE 0
+#define FAKE_CLEN (0.1)
+
+
+typedef enum {
+ MOSFLM_INPUT_NONE,
+ MOSFLM_INPUT_LINE,
+ MOSFLM_INPUT_PROMPT
+} MOSFLMInputType;
+
+
+
+struct mosflm_private {
+ IndexingMethod indm;
+ UnitCell *template;
+};
+
+
+struct mosflm_data {
+
+ /* MOSFLM auto-indexing low-level stuff */
+ int pty;
+ pid_t pid;
+ char *rbuffer;
+ int rbufpos;
+ int rbuflen;
+
+ /* MOSFLM high-level stuff */
+ char newmatfile[128];
+ char imagefile[128];
+ char sptfile[128];
+ int step;
+ int finished_ok;
+ int done;
+ int success;
+
+ struct mosflm_private *mp;
+
+};
+
+static int check_mosflm_cell(struct mosflm_private *mp, struct image *image,
+ UnitCell *cell)
+{
+ Crystal *cr;
+
+ /* If we sent lattice information, make sure that we got back what we
+ * asked for, not (e.g.) some "H" version of a rhombohedral R cell */
+ if ( mp->indm & INDEXING_USE_LATTICE_TYPE ) {
+
+ LatticeType latt_m, latt_r;
+ char cen_m, cen_r;
+
+ /* What we asked for */
+ latt_r = cell_get_lattice_type(mp->template);
+ cen_r = cell_get_centering(mp->template);
+
+ /* What we got back */
+ latt_m = cell_get_lattice_type(cell);
+ cen_m = cell_get_centering(cell);
+
+ if ( latt_r != latt_m ) {
+ ERROR("Lattice type produced by MOSFLM (%i) does not "
+ "match what was requested (%i). "
+ "Please report this.\n", latt_m, latt_r);
+ return 0;
+ }
+
+ if ( (latt_m != L_MONOCLINIC) && (cen_r != cen_m) ) {
+ ERROR("Centering produced by MOSFLM (%c) does not "
+ "match what was requested (%c). "
+ "Please report this.\n", cen_m, cen_r);
+ return 0;
+ }
+ /* If it's monoclinic, see the warning in mosflm_prepare() */
+
+ }
+
+ cr = crystal_new();
+ if ( cr == NULL ) {
+ ERROR("Failed to allocate crystal.\n");
+ return 0;
+ }
+
+ crystal_set_cell(cr, cell);
+
+ image_add_crystal(image, cr);
+
+ return 1;
+}
+
+
+static void mosflm_parseline(const char *line, struct image *image,
+ struct mosflm_data *dirax)
+{
+ if ( MOSFLM_VERBOSE || (strncmp(line, "Invocation:", 11) == 0) ) {
+ char *copy;
+ int i;
+
+ copy = strdup(line);
+ for ( i=0; i<strlen(copy); i++ ) {
+ if ( copy[i] == '\r' ) copy[i]='r';
+ if ( copy[i] == '\n' ) copy[i]='\0';
+ }
+ STATUS("MOSFLM: %s\n", copy);
+ free(copy);
+ }
+}
+
+
+/* This is the opposite of mosflm_spacegroup_for_lattice() below.
+ * Note that this is not general, just a set of rules for interpreting MOSFLM's
+ * output. */
+static LatticeType mosflm_spacegroup_to_lattice(const char *sg,
+ char *ua, char *cen)
+{
+ LatticeType latt;
+
+ *cen = sg[0];
+
+ if ( sg[1] == '1' ) {
+ latt = L_TRICLINIC;
+ *ua = '*';
+ } else if ( strncmp(sg+1, "23", 2) == 0 ) {
+ latt = L_CUBIC;
+ *ua = '*';
+ } else if ( strncmp(sg+1, "222", 3) == 0 ) {
+ latt = L_ORTHORHOMBIC;
+ *ua = '*';
+ } else if ( sg[1] == '2' ) {
+ latt = L_MONOCLINIC;
+ *ua = 'b';
+ } else if ( sg[1] == '4' ) {
+ latt = L_TETRAGONAL;
+ *ua = 'c';
+ } else if ( sg[1] == '6' ) {
+ latt = L_HEXAGONAL;
+ *ua = 'c';
+ } else if ( sg[1] == '3' ) {
+ if ( (sg[0] == 'H') || (sg[0] == 'P') ) {
+ latt = L_HEXAGONAL;
+ *ua = 'c';
+ } else {
+ latt = L_RHOMBOHEDRAL;
+ *ua = '*';
+ }
+ } else {
+ ERROR("Couldn't understand '%s'\n", sg);
+ latt = L_TRICLINIC;
+ }
+
+ return latt;
+}
+
+
+static int read_newmat(struct mosflm_data *mosflm, const char *filename,
+ struct image *image)
+{
+ FILE *fh;
+ float asx, asy, asz;
+ float bsx, bsy, bsz;
+ float csx, csy, csz;
+ int n;
+ double c;
+ UnitCell *cell;
+ char symm[32];
+ char *rval;
+ int i;
+ char cen;
+ LatticeType latt;
+ char ua = '?';
+
+ fh = fopen(filename, "r");
+ if ( fh == NULL ) {
+ return 1;
+ }
+ n = fscanf(fh, "%f %f %f\n", &asx, &bsx, &csx);
+ n += fscanf(fh, "%f %f %f\n", &asy, &bsy, &csy);
+ n += fscanf(fh, "%f %f %f\n", &asz, &bsz, &csz);
+ if ( n != 9 ) {
+ STATUS("Fewer than 9 parameters found in NEWMAT file.\n");
+ return 1;
+ }
+
+ /* Skip the next six lines */
+ for ( i=0; i<6; i++ ) {
+ char tmp[1024];
+ rval = fgets(tmp, 1024, fh);
+ if ( rval == NULL ) {
+ ERROR("Failed to read newmat file.\n");
+ return 1;
+ }
+ }
+
+ rval = fgets(symm, 32, fh);
+ if ( rval == NULL ) {
+ ERROR("Failed to read newmat file.\n");
+ return 1;
+ }
+
+ fclose(fh);
+
+ chomp(symm);
+ if ( strncmp(symm, "SYMM ", 5) != 0 ) {
+ ERROR("Bad 'SYMM' line from MOSFLM.\n");
+ return 1;
+ }
+ //STATUS("MOSFLM says '%s'\n", symm);
+ latt = mosflm_spacegroup_to_lattice(symm+5, &ua, &cen);
+
+ /* MOSFLM "A" matrix is multiplied by lambda, so fix this */
+ c = 1.0/image->lambda;
+
+ cell = cell_new();
+
+ /* The relationship between the coordinates in the spot file and the
+ * resulting matrix is diabolically complicated. This transformation
+ * seems to work, but is not derived by working through all the
+ * transformations. */
+ cell_set_reciprocal(cell,
+ -asy*c, -asz*c, asx*c,
+ -bsy*c, -bsz*c, bsx*c,
+ -csy*c, -csz*c, csx*c);
+ cell_set_centering(cell, cen);
+ cell_set_lattice_type(cell, latt);
+ cell_set_unique_axis(cell, ua);
+ //STATUS("My cell:\n");
+ //cell_print(cell);
+
+ if ( check_mosflm_cell(mosflm->mp, image, cell) ) {
+ mosflm->success = 1;
+ mosflm->done = 1;
+ }
+
+ return 0;
+}
+
+
+/* Write .spt file for mosflm */
+static void write_spt(struct image *image, const char *filename)
+{
+ FILE *fh;
+ int i;
+ int n;
+
+ fh = fopen(filename, "w");
+ if ( !fh ) {
+ ERROR("Couldn't open temporary file '%s'\n", filename);
+ return;
+ }
+
+ /* Number of pixels in x, number of pixels in y, pixel size (mm),
+ * YSCALE, OMEGA */
+ fputs("1 1 0.0 1.0 0.0\n", fh);
+
+ /* INVERTX, ISWUNG */
+ fputs("0 1\n", fh);
+
+ /* XBEAM, YBEAM */
+ fputs("0.0 0.0\n", fh);
+
+ n = image_feature_count(image->features);
+ for ( i=0; i<n; i++ ) {
+
+ struct imagefeature *f;
+ double ttx, tty, x, y;
+
+ f = image_get_feature(image->features, i);
+ if ( f == NULL ) continue;
+
+ ttx = angle_between_2d(0.0, 1.0,
+ f->rx, 1.0/image->lambda + f->rz);
+ tty = angle_between_2d(0.0, 1.0,
+ f->ry, 1.0/image->lambda + f->rz);
+ if ( f->rx < 0.0 ) ttx *= -1.0;
+ if ( f->ry < 0.0 ) tty *= -1.0;
+ x = -tan(ttx)*FAKE_CLEN;
+ y = tan(tty)*FAKE_CLEN;
+
+ fprintf(fh, "%10.2f %10.2f 0.0 0.0 1000.0 10.0\n",
+ x*1e3, y*1e3);
+
+ }
+
+ fputs("-999.0 -999.0 -999.0 -999.0 -999.0 -999.0\n", fh);
+
+ fclose(fh);
+}
+
+
+/* Write a dummy 1x1 pixel image file for MOSFLM. Without post refinement,
+ * MOSFLM will ignore this, but it must be present. */
+static void write_img(struct image *image, const char *filename)
+{
+ FILE *fh;
+ unsigned short int *intimage;
+
+ intimage = malloc(sizeof(unsigned short int));
+ intimage[0] = 1;
+
+ fh = fopen(filename, "w");
+ if ( !fh ) {
+ ERROR("Couldn't open temporary file '%s'\n", filename);
+ return;
+ }
+
+ /* Write header */
+ fprintf(fh, "{\nHEADER_BYTES=512;\n");
+ fprintf(fh, "BYTE_ORDER=little_endian;\n");
+ fprintf(fh, "TYPE=unsigned_short;\n");
+ fprintf(fh, "DIM=2;\n");
+ fprintf(fh, "SIZE1=1;\n");
+ fprintf(fh, "SIZE2=1;\n");
+ fprintf(fh, "}\n");
+
+ /* Header padding */
+ while ( ftell(fh) < 512 ) fprintf(fh," ");
+
+ fwrite(intimage, sizeof(unsigned short int), 1, fh);
+ free(intimage);
+ fclose(fh);
+}
+
+
+static void mosflm_sendline(const char *line, struct mosflm_data *mosflm)
+{
+ #if MOSFLM_VERBOSE
+ char *copy;
+ int i;
+
+ copy = strdup(line);
+ for ( i=0; i<strlen(copy); i++ ) {
+ if ( copy[i] == '\r' ) copy[i]='\0';
+ if ( copy[i] == '\n' ) copy[i]='\0';
+ }
+ STATUS("To MOSFLM: '%s'\n", copy);
+ free(copy);
+ #endif
+
+ if ( write(mosflm->pty, line, strlen(line)) == -1 ) {
+ ERROR("write() to MOSFLM failed: %s\n", strerror(errno));
+ }
+}
+
+
+/* Turn what we know about the unit cell into something which we can give to
+ * MOSFLM to make it give us only indexing results compatible with the cell. */
+static char *mosflm_spacegroup_for_lattice(UnitCell *cell)
+{
+ LatticeType latt;
+ char centering;
+ char *g = NULL;
+ char *result;
+ char ua;
+
+ latt = cell_get_lattice_type(cell);
+ centering = cell_get_centering(cell);
+ ua = cell_get_unique_axis(cell);
+
+ switch ( latt )
+ {
+ case L_TRICLINIC :
+ g = "1";
+ break;
+
+ case L_MONOCLINIC :
+ switch ( ua ) {
+ case 'a' : g = "211"; break;
+ case 'b' : g = "121"; break;
+ case 'c' : g = "112"; break;
+ }
+ break;
+
+ case L_ORTHORHOMBIC :
+ g = "222";
+ break;
+
+ case L_TETRAGONAL :
+ g = "4";
+ break;
+
+ case L_RHOMBOHEDRAL :
+ g = "3";
+ break;
+
+ case L_HEXAGONAL :
+ if ( centering != 'H' ) {
+ g = "6";
+ } else {
+ g = "3";
+ }
+ break;
+
+ case L_CUBIC :
+ g = "23";
+ break;
+ }
+ assert(g != NULL);
+
+ result = malloc(32);
+ if ( result == NULL ) return NULL;
+
+ snprintf(result, 31, "%c%s", centering, g);
+
+ return result;
+}
+
+
+static void mosflm_send_next(struct image *image, struct mosflm_data *mosflm)
+{
+ char tmp[256];
+ char cen;
+ double wavelength;
+ double a = 0, b = 0, c = 0, alpha = 0, beta = 0, gamma = 0;
+
+ switch ( mosflm->step )
+ {
+ case 1 :
+ /* Backwards-compatible workaround for different Mosflm behaviour
+ * in version 7.2.2 */
+ mosflm_sendline("DETECTOR ADSC\n", mosflm);
+ break;
+
+ case 2 :
+ mosflm_sendline("DETECTOR ROTATION HORIZONTAL"
+ " ANTICLOCKWISE ORIGIN LL FAST HORIZONTAL"
+ " RECTANGULAR\n", mosflm);
+ break;
+
+ case 3 :
+ if ( (mosflm->mp->indm & INDEXING_USE_LATTICE_TYPE)
+ && (mosflm->mp->template != NULL) )
+ {
+ char *symm;
+
+ if ( cell_get_lattice_type(mosflm->mp->template)
+ == L_RHOMBOHEDRAL ) {
+ mosflm_sendline("CRYSTAL R\n", mosflm);
+ }
+
+ symm = mosflm_spacegroup_for_lattice(mosflm->mp->template);
+ snprintf(tmp, 255, "SYMM %s\n", symm);
+ //STATUS("Asking MOSFLM for '%s'\n", symm);
+ free(symm);
+ mosflm_sendline(tmp, mosflm);
+
+ } else {
+ mosflm_sendline("\n", mosflm);
+ }
+ break;
+
+ case 4 :
+ snprintf(tmp, 255, "DISTANCE %f\n", FAKE_CLEN*1e3);
+ mosflm_sendline(tmp, mosflm);
+ break;
+
+ case 5 :
+ mosflm_sendline("BEAM 0.0 0.0\n", mosflm);
+ break;
+
+ case 6 :
+ wavelength = image->lambda*1e10;
+ snprintf(tmp, 255, "WAVELENGTH %10.5f\n", wavelength);
+ mosflm_sendline(tmp, mosflm);
+ break;
+
+ case 7 :
+ snprintf(tmp, 255, "NEWMAT %s\n", mosflm->newmatfile);
+ mosflm_sendline(tmp, mosflm);
+ break;
+
+ case 8 :
+ snprintf(tmp, 255, "IMAGE %s phi 0 0\n", mosflm->imagefile);
+ mosflm_sendline(tmp, mosflm);
+ break;
+
+ case 9 :
+ if ( mosflm->mp->indm & INDEXING_USE_CELL_PARAMETERS ) {
+ cell_get_parameters(mosflm->mp->template,
+ &a, &b, &c, &alpha, &beta, &gamma);
+ cen = cell_get_centering(mosflm->mp->template);
+ /* Old MOSFLM simply ignores CELL and CENTERING subkeywords.
+ * So this doesn't do any harm.
+ * TODO: but still better to show WARNING if MOSFLM is old. */
+ snprintf(tmp, 255, "AUTOINDEX DPS FILE %s IMAGE 1 "
+ "MAXCELL 1000 REFINE "
+ "CELL %.1f %.1f %.1f %.1f %.1f %.1f "
+ "CENTERING %c\n",
+ mosflm->sptfile, a*1e10, b*1e10, c*1e10,
+ rad2deg(alpha), rad2deg(beta), rad2deg(gamma),
+ cen);
+ } else {
+ snprintf(tmp, 255, "AUTOINDEX DPS FILE %s IMAGE 1 "
+ "MAXCELL 1000 REFINE\n", mosflm->sptfile);
+ }
+ mosflm_sendline(tmp, mosflm);
+ break;
+
+ case 10 :
+ mosflm_sendline("GO\n", mosflm);
+ mosflm->finished_ok = 1;
+ break;
+
+ default:
+ mosflm_sendline("exit\n", mosflm);
+ return;
+ }
+
+ mosflm->step++;
+}
+
+
+static int mosflm_readable(struct image *image, struct mosflm_data *mosflm)
+{
+ int rval;
+ int no_string = 0;
+
+ rval = read(mosflm->pty, mosflm->rbuffer+mosflm->rbufpos,
+ mosflm->rbuflen-mosflm->rbufpos);
+ if ( (rval == -1) || (rval == 0) ) return 1;
+
+ mosflm->rbufpos += rval;
+ assert(mosflm->rbufpos <= mosflm->rbuflen);
+
+ while ( (!no_string) && (mosflm->rbufpos > 0) ) {
+
+ int i;
+ int block_ready = 0;
+ MOSFLMInputType type = MOSFLM_INPUT_NONE;
+
+ /* See if there's a full line in the buffer yet */
+ for ( i=0; i<mosflm->rbufpos-1; i++ ) {
+ /* Means the last value looked at is rbufpos-2 */
+
+ /* Is there a prompt in the buffer? */
+ if ( (i+10 <= mosflm->rbufpos)
+ && (!strncmp(mosflm->rbuffer+i, "MOSFLM => ", 10)) ) {
+ block_ready = 1;
+ type = MOSFLM_INPUT_PROMPT;
+ break;
+ }
+
+ if ( (mosflm->rbuffer[i] == '\r')
+ && (mosflm->rbuffer[i+1] == '\n') ) {
+ block_ready = 1;
+ type = MOSFLM_INPUT_LINE;
+ break;
+ }
+
+ }
+
+ if ( block_ready ) {
+
+ unsigned int new_rbuflen;
+ unsigned int endbit_length;
+ char *block_buffer = NULL;
+
+ switch ( type ) {
+
+ case MOSFLM_INPUT_LINE :
+ block_buffer = malloc(i+1);
+ memcpy(block_buffer, mosflm->rbuffer, i);
+ block_buffer[i] = '\0';
+
+ if ( block_buffer[0] == '\r' ) {
+ memmove(block_buffer, block_buffer+1, i);
+ }
+
+ mosflm_parseline(block_buffer, image, mosflm);
+ free(block_buffer);
+ endbit_length = i+2;
+ break;
+
+ case MOSFLM_INPUT_PROMPT :
+ mosflm_send_next(image, mosflm);
+ endbit_length = i+7;
+ break;
+
+ default :
+
+ /* Obviously, this never happens :) */
+ ERROR("Unrecognised MOSFLM input mode! "
+ "I don't know how to understand MOSFLM\n");
+ return 1;
+
+ }
+
+ /* Now the block's been parsed, it should be
+ * forgotten about */
+ memmove(mosflm->rbuffer,
+ mosflm->rbuffer + endbit_length,
+ mosflm->rbuflen - endbit_length);
+
+ /* Subtract the number of bytes removed */
+ mosflm->rbufpos = mosflm->rbufpos
+ - endbit_length;
+ new_rbuflen = mosflm->rbuflen - endbit_length;
+ if ( new_rbuflen == 0 ) new_rbuflen = 256;
+ mosflm->rbuffer = realloc(mosflm->rbuffer,
+ new_rbuflen);
+ mosflm->rbuflen = new_rbuflen;
+
+ } else {
+
+ if ( mosflm->rbufpos==mosflm->rbuflen ) {
+
+ /* More buffer space is needed */
+ mosflm->rbuffer = realloc(
+ mosflm->rbuffer,
+ mosflm->rbuflen + 256);
+ mosflm->rbuflen = mosflm->rbuflen + 256;
+ /* The new space gets used at the next
+ * read, shortly... */
+
+ }
+ no_string = 1;
+
+ }
+
+ }
+
+ return 0;
+}
+
+
+int run_mosflm(struct image *image, void *ipriv)
+{
+ struct mosflm_data *mosflm;
+ unsigned int opts;
+ int status;
+ int rval;
+
+ mosflm = malloc(sizeof(struct mosflm_data));
+ if ( mosflm == NULL ) {
+ ERROR("Couldn't allocate memory for MOSFLM data.\n");
+ return 0;
+ }
+
+ snprintf(mosflm->imagefile, 127, "xfel_001.img");
+ write_img(image, mosflm->imagefile); /* Dummy image */
+
+ snprintf(mosflm->sptfile, 127, "xfel_001.spt");
+ write_spt(image, mosflm->sptfile);
+
+ snprintf(mosflm->newmatfile, 127, "xfel.newmat");
+ remove(mosflm->newmatfile);
+
+ mosflm->pid = forkpty(&mosflm->pty, NULL, NULL, NULL);
+
+ if ( mosflm->pid == -1 ) {
+ ERROR("Failed to fork for MOSFLM: %s\n", strerror(errno));
+ free(mosflm);
+ return 0;
+ }
+ if ( mosflm->pid == 0 ) {
+
+ /* Child process: invoke MOSFLM */
+ struct termios t;
+
+ /* Turn echo off */
+ tcgetattr(STDIN_FILENO, &t);
+ t.c_lflag &= ~(ECHO | ECHOE | ECHOK | ECHONL);
+ tcsetattr(STDIN_FILENO, TCSANOW, &t);
+
+ execlp("mosflm", "mosflm", (char *)NULL);
+ execlp("ipmosflm", "ipmosflm", (char *)NULL);
+ ERROR("Invocation: Failed to invoke MOSFLM: %s\n",
+ strerror(errno));
+ _exit(0);
+
+ }
+
+ mosflm->rbuffer = malloc(256);
+ mosflm->rbuflen = 256;
+ mosflm->rbufpos = 0;
+
+ /* Set non-blocking */
+ opts = fcntl(mosflm->pty, F_GETFL);
+ fcntl(mosflm->pty, F_SETFL, opts | O_NONBLOCK);
+
+ mosflm->step = 1; /* This starts the "initialisation" procedure */
+ mosflm->finished_ok = 0;
+ mosflm->mp = (struct mosflm_private *)ipriv;
+ mosflm->done = 0;
+ mosflm->success = 0;
+
+ rval = 0;
+ do {
+
+ fd_set fds;
+ struct timeval tv;
+ int sval;
+
+ FD_ZERO(&fds);
+ FD_SET(mosflm->pty, &fds);
+
+ tv.tv_sec = 30;
+ tv.tv_usec = 0;
+
+ sval = select(mosflm->pty+1, &fds, NULL, NULL, &tv);
+
+ if ( sval == -1 ) {
+
+ const int err = errno;
+
+ switch ( err ) {
+
+ case EINTR:
+ STATUS("Restarting select()\n");
+ break;
+
+ default:
+ ERROR("select() failed: %s\n", strerror(err));
+ rval = 1;
+
+ }
+
+ } else if ( sval != 0 ) {
+ rval = mosflm_readable(image, mosflm);
+ } else {
+ ERROR("No response from MOSFLM..\n");
+ rval = 1;
+ }
+
+ } while ( !rval );
+
+ close(mosflm->pty);
+ free(mosflm->rbuffer);
+ waitpid(mosflm->pid, &status, 0);
+
+ if ( mosflm->finished_ok == 0 ) {
+ ERROR("MOSFLM doesn't seem to be working properly.\n");
+ } else {
+ /* Read the mosflm NEWMAT file and get cell if found */
+ read_newmat(mosflm, mosflm->newmatfile, image);
+ }
+
+ rval = mosflm->success;
+ free(mosflm);
+ return rval;
+}
+
+
+void *mosflm_prepare(IndexingMethod *indm, UnitCell *cell)
+{
+ struct mosflm_private *mp;
+
+ if ( mosflm_probe(cell) == NULL ) {
+ ERROR("Mosflm does not appear to run properly.\n");
+ ERROR("Please check your Mosflm installation.\n");
+ return NULL;
+ }
+
+ /* Flags that MOSFLM knows about */
+ *indm &= INDEXING_METHOD_MASK
+ | INDEXING_USE_LATTICE_TYPE | INDEXING_USE_CELL_PARAMETERS;
+
+ if ( (cell != NULL) && (cell_get_centering(cell) == 'I')
+ && (cell_get_lattice_type(cell) == L_MONOCLINIC) )
+ {
+ ERROR("WARNING: Mosflm always gives the monoclinic C cell in "
+ "preference to the monoclinic I cell choice.\n");
+ ERROR("To get a higher indexing rate, convert your cell to the "
+ "monoclinic C cell choice.\n");
+ }
+
+ mp = malloc(sizeof(struct mosflm_private));
+ if ( mp == NULL ) return NULL;
+
+ mp->template = cell;
+ mp->indm = *indm;
+
+ return (IndexingPrivate *)mp;
+}
+
+
+void mosflm_cleanup(void *pp)
+{
+ struct mosflm_private *p;
+ p = (struct mosflm_private *)pp;
+ free(p);
+}
+
+
+static void chop_word(char *s)
+{
+ int i;
+ size_t l = strlen(s);
+ for ( i=0; i<l; i++ ) {
+ if ( s[i] == ' ' ) {
+ s[i] = '\0';
+ return;
+ }
+ }
+}
+
+
+const char *mosflm_probe(UnitCell *cell)
+{
+ pid_t pid;
+ int pty;
+ int status;
+ FILE *fh;
+ char line[1024];
+ int ok = 0;
+ int l;
+
+ pid = forkpty(&pty, NULL, NULL, NULL);
+ if ( pid == -1 ) {
+ return NULL;
+ }
+ if ( pid == 0 ) {
+
+ /* Child process */
+ struct termios t;
+
+ /* Turn echo off */
+ tcgetattr(STDIN_FILENO, &t);
+ t.c_lflag &= ~(ECHO | ECHOE | ECHOK | ECHONL);
+ tcsetattr(STDIN_FILENO, TCSANOW, &t);
+
+ execlp("mosflm", "mosflm", (char *)NULL);
+ execlp("ipmosflm", "ipmosflm", (char *)NULL);
+ _exit(1);
+
+ }
+
+ fh = fdopen(pty, "r");
+
+ for ( l=0; l<10; l++ ) {
+ char *pos;
+ if ( fgets(line, 1024, fh) != NULL ) {
+ pos = strstr(line, "Mosflm version ");
+ if ( pos != NULL ) {
+ char *vers = pos+15;
+ ok = 1;
+ chop_word(vers);
+ /* FIXME: Set capabilities based on version */
+ }
+ }
+ }
+
+ fclose(fh);
+ close(pty);
+ waitpid(pid, &status, 0);
+
+ if ( !ok ) return NULL;
+
+ if ( cell_has_parameters(cell) ) return "mosflm-cell-nolatt,mosflm-latt-nocell";
+ if ( cell != NULL ) return "mosflm-latt-nocell";
+ return "mosflm-nolatt-nocell";
+}
diff --git a/libcrystfel/src/indexers/mosflm.h b/libcrystfel/src/indexers/mosflm.h
new file mode 100644
index 00000000..9e956066
--- /dev/null
+++ b/libcrystfel/src/indexers/mosflm.h
@@ -0,0 +1,56 @@
+/*
+ * mosflm.h
+ *
+ * Invoke the DPS auto-indexing algorithm through MOSFLM
+ *
+ * Copyright © 2012-2020 Deutsches Elektronen-Synchrotron DESY,
+ * a research centre of the Helmholtz Association.
+ * Copyright © 2012 Richard Kirian
+ *
+ * Authors:
+ * 2010 Richard Kirian <rkirian@asu.edu>
+ * 2012-2014,2017 Thomas White <taw@physics.org>
+ *
+ * This file is part of CrystFEL.
+ *
+ * CrystFEL is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * CrystFEL is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with CrystFEL. If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+#ifndef MOSFLM_H
+#define MOSFLM_H
+
+#include "index.h"
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+/**
+ * \file mosflm.h
+ * MOSFLM indexer interface
+ */
+
+extern int run_mosflm(struct image *image, void *ipriv);
+
+extern void *mosflm_prepare(IndexingMethod *indm, UnitCell *cell);
+extern const char *mosflm_probe(UnitCell *cell);
+
+extern void mosflm_cleanup(void *pp);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif /* MOSFLM_H */
diff --git a/libcrystfel/src/indexers/pinkindexer.c b/libcrystfel/src/indexers/pinkindexer.c
new file mode 100644
index 00000000..67bd48f5
--- /dev/null
+++ b/libcrystfel/src/indexers/pinkindexer.c
@@ -0,0 +1,646 @@
+/*
+ * pinkindexer.c
+ *
+ * Interface to PinkIndexer
+ *
+ * Copyright © 2017-2020 Deutsches Elektronen-Synchrotron DESY,
+ * a research centre of the Helmholtz Association.
+ *
+ * Authors:
+ * 2017-2019 Yaroslav Gevorkov <yaroslav.gevorkov@desy.de>
+ *
+ * This file is part of CrystFEL.
+ *
+ * CrystFEL is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * CrystFEL is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with CrystFEL. If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include "pinkindexer.h"
+
+
+#include <stdlib.h>
+#include <sys/errno.h>
+#include <argp.h>
+
+#include "utils.h"
+#include "cell-utils.h"
+#include "peaks.h"
+
+struct pinkIndexer_options {
+ unsigned int considered_peaks_count;
+ unsigned int angle_resolution;
+ unsigned int refinement_type;
+ float maxResolutionForIndexing_1_per_A;
+ float tolerance;
+ int multi;
+ int thread_count;
+ int min_peaks;
+ int no_check_indexed;
+ float reflectionRadius; /* In m^-1 */
+ float customPhotonEnergy;
+ float customBandwidth;
+ float maxRefinementDisbalance;
+};
+
+#ifdef HAVE_PINKINDEXER
+
+#include <pinkIndexer/adaptions/crystfel/Lattice.h>
+#include <pinkIndexer/adaptions/crystfel/ExperimentSettings.h>
+#include <pinkIndexer/adaptions/crystfel/PinkIndexer.h>
+
+#define MAX_MULTI_LATTICE_COUNT 8
+
+struct pinkIndexer_private_data {
+ PinkIndexer *pinkIndexer;
+ reciprocalPeaks_1_per_A_t reciprocalPeaks_1_per_A;
+ float *intensities;
+
+ IndexingMethod indm;
+ UnitCell *cellTemplate;
+ int threadCount;
+ int multi;
+ int min_peaks;
+
+ int no_check_indexed;
+
+ float maxRefinementDisbalance;
+
+ IntegerMatrix *centeringTransformation;
+ LatticeTransform_t latticeReductionTransform;
+};
+
+//static void reduceCell(UnitCell* cell, LatticeTransform_t* appliedReductionTransform);
+//static void restoreCell(UnitCell *cell, LatticeTransform_t* appliedReductionTransform);
+static void reduceReciprocalCell(UnitCell* cell, LatticeTransform_t* appliedReductionTransform);
+static void restoreReciprocalCell(UnitCell *cell, LatticeTransform_t* appliedReductionTransform);
+static void makeRightHanded(UnitCell* cell);
+static void update_detector(struct detector *det, double xoffs, double yoffs);
+
+int run_pinkIndexer(struct image *image, void *ipriv)
+{
+ struct pinkIndexer_private_data* pinkIndexer_private_data = (struct pinkIndexer_private_data*) ipriv;
+ reciprocalPeaks_1_per_A_t* reciprocalPeaks_1_per_A = &(pinkIndexer_private_data->reciprocalPeaks_1_per_A);
+ float *intensities = pinkIndexer_private_data->intensities;
+
+ int peakCountMax = image_feature_count(image->features);
+ if (peakCountMax < 5) {
+ int goodLatticesCount = 0;
+ return goodLatticesCount;
+ }
+ reciprocalPeaks_1_per_A->peakCount = 0;
+ for (int i = 0; i < peakCountMax && i < MAX_PEAK_COUNT_FOR_INDEXER; i++) {
+ struct imagefeature *f;
+ f = image_get_feature(image->features, i);
+ if (f == NULL) {
+ continue;
+ }
+
+ reciprocalPeaks_1_per_A->coordinates_x[reciprocalPeaks_1_per_A->peakCount] = f->rz * 1e-10;
+ reciprocalPeaks_1_per_A->coordinates_y[reciprocalPeaks_1_per_A->peakCount] = f->rx * 1e-10;
+ reciprocalPeaks_1_per_A->coordinates_z[reciprocalPeaks_1_per_A->peakCount] = f->ry * 1e-10;
+ intensities[reciprocalPeaks_1_per_A->peakCount] = (float) (f->intensity);
+ reciprocalPeaks_1_per_A->peakCount++;
+ }
+ int indexed = 0;
+ Lattice_t indexedLattice[MAX_MULTI_LATTICE_COUNT];
+ float center_shift[MAX_MULTI_LATTICE_COUNT][2];
+
+
+
+ do {
+ int peakCount = reciprocalPeaks_1_per_A->peakCount;
+ int matchedPeaksCount = PinkIndexer_indexPattern(pinkIndexer_private_data->pinkIndexer,
+ &(indexedLattice[indexed]), center_shift[indexed], reciprocalPeaks_1_per_A, intensities,
+ pinkIndexer_private_data->maxRefinementDisbalance,
+ pinkIndexer_private_data->threadCount);
+
+ if(matchedPeaksCount == -1){
+ STATUS("WARNING: Indexing solution was rejected due to too large disbalance of the refinement."
+ "If you see this message often, check the documentation for the parameter "
+ "--pinkIndexer-max-refinement-disbalance\n");
+
+ matchedPeaksCount = 0;
+ }
+
+ printf("matchedPeaksCount %d from %d\n",matchedPeaksCount,peakCount);
+ if ((matchedPeaksCount >= 25 && matchedPeaksCount >= peakCount * 0.30)
+ || matchedPeaksCount >= peakCount * 0.4
+ || matchedPeaksCount >= 70
+ || pinkIndexer_private_data->no_check_indexed == 1)
+ {
+ UnitCell *uc;
+ uc = cell_new();
+
+ Lattice_t *l = &(indexedLattice[indexed]);
+
+ cell_set_reciprocal(uc, l->ay * 1e10, l->az * 1e10, l->ax * 1e10,
+ l->by * 1e10, l->bz * 1e10, l->bx * 1e10,
+ l->cy * 1e10, l->cz * 1e10, l->cx * 1e10);
+
+ restoreReciprocalCell(uc, &pinkIndexer_private_data->latticeReductionTransform);
+
+ UnitCell *new_cell_trans = cell_transform_intmat(uc, pinkIndexer_private_data->centeringTransformation);
+ cell_free(uc);
+ uc = new_cell_trans;
+
+ cell_set_lattice_type(new_cell_trans, cell_get_lattice_type(pinkIndexer_private_data->cellTemplate));
+ cell_set_centering(new_cell_trans, cell_get_centering(pinkIndexer_private_data->cellTemplate));
+ cell_set_unique_axis(new_cell_trans, cell_get_unique_axis(pinkIndexer_private_data->cellTemplate));
+
+ if (validate_cell(uc)) {
+ ERROR("pinkIndexer: problem with returned cell!\n");
+ }
+
+ Crystal * cr = crystal_new();
+ if (cr == NULL) {
+ ERROR("Failed to allocate crystal.\n");
+ return 0;
+ }
+ crystal_set_cell(cr, uc);
+ crystal_set_det_shift(cr, center_shift[indexed][0], center_shift[indexed][1]);
+ update_detector(image->det, center_shift[indexed][0], center_shift[indexed][1]);
+ image_add_crystal(image, cr);
+ indexed++;
+
+ } else {
+ break;
+ }
+ } while (pinkIndexer_private_data->multi
+ && indexed <= MAX_MULTI_LATTICE_COUNT
+ && reciprocalPeaks_1_per_A->peakCount >= pinkIndexer_private_data->min_peaks);
+
+ return indexed;
+}
+
+void *pinkIndexer_prepare(IndexingMethod *indm, UnitCell *cell,
+ struct pinkIndexer_options *pinkIndexer_opts,
+ const DataTemplate *dtempl)
+{
+ if ( beam->photon_energy_from != NULL && pinkIndexer_opts->customPhotonEnergy <= 0) {
+ ERROR("For pinkIndexer, the photon_energy must be defined as a "
+ "constant in the geometry file or as a parameter (see --pinkIndexer-override-photon-energy)\n");
+ return NULL;
+ }
+ if ( (det->panels[0].clen_from != NULL) && pinkIndexer_opts->refinement_type ==
+ REFINEMENT_TYPE_firstFixedThenVariableLatticeParametersCenterAdjustmentMultiSeed) {
+ ERROR("Using center refinement makes it necessary to have the detector distance fixed in the geometry file!");
+ return NULL;
+ }
+ if(cell == NULL){
+ ERROR("Pink indexer needs a unit cell file to be specified!");
+ return NULL;
+ }
+
+ struct pinkIndexer_private_data* pinkIndexer_private_data = malloc(sizeof(struct pinkIndexer_private_data));
+ allocReciprocalPeaks(&(pinkIndexer_private_data->reciprocalPeaks_1_per_A));
+ pinkIndexer_private_data->intensities = malloc(MAX_PEAK_COUNT_FOR_INDEXER * sizeof(float));
+ pinkIndexer_private_data->indm = *indm;
+ pinkIndexer_private_data->cellTemplate = cell;
+ pinkIndexer_private_data->threadCount = pinkIndexer_opts->thread_count;
+ pinkIndexer_private_data->multi = pinkIndexer_opts->multi;
+ pinkIndexer_private_data->min_peaks = pinkIndexer_opts->min_peaks;
+ pinkIndexer_private_data->no_check_indexed = pinkIndexer_opts->no_check_indexed;
+ pinkIndexer_private_data->maxRefinementDisbalance = pinkIndexer_opts->maxRefinementDisbalance;
+
+ UnitCell* primitiveCell = uncenter_cell(cell, &pinkIndexer_private_data->centeringTransformation, NULL);
+
+ //reduceCell(primitiveCell, &pinkIndexer_private_data->latticeReductionTransform);
+ reduceReciprocalCell(primitiveCell, &pinkIndexer_private_data->latticeReductionTransform);
+
+ double asx, asy, asz, bsx, bsy, bsz, csx, csy, csz;
+ int ret = cell_get_reciprocal(primitiveCell, &asx, &asy, &asz, &bsx, &bsy, &bsz, &csx, &csy, &csz);
+ if (ret != 0) {
+ ERROR("cell_get_reciprocal did not finish properly!");
+ }
+
+ Lattice_t lattice = { .ax = asz * 1e-10, .ay = asx * 1e-10, .az = asy * 1e-10,
+ .bx = bsz * 1e-10, .by = bsx * 1e-10, .bz = bsy * 1e-10,
+ .cx = csz * 1e-10, .cy = csx * 1e-10, .cz = csy * 1e-10 };
+
+ float detectorDistance_m;
+ if ( det->panels[0].clen_from != NULL ) {
+ detectorDistance_m = 0.25; /* fake value */
+ } else {
+ detectorDistance_m = det->panels[0].clen + det->panels[0].coffset;
+ }
+
+ float beamEenergy_eV = beam->photon_energy;
+ float nonMonochromaticity = beam->bandwidth*5;
+ if(pinkIndexer_opts->customPhotonEnergy > 0){
+ beamEenergy_eV = pinkIndexer_opts->customPhotonEnergy;
+ }
+ if(pinkIndexer_opts->customBandwidth >= 0){
+ nonMonochromaticity = pinkIndexer_opts->customBandwidth;
+ }
+
+ float reflectionRadius_1_per_A;
+ if (pinkIndexer_opts->reflectionRadius < 0) {
+ reflectionRadius_1_per_A = 0.02
+ * sqrt(lattice.ax * lattice.ax + lattice.ay * lattice.ay + lattice.az * lattice.az);
+ }
+ else {
+ reflectionRadius_1_per_A = pinkIndexer_opts->reflectionRadius * 1e10; /* m^-1 to A^-1*/
+ }
+
+ if(beamEenergy_eV > 75000 && nonMonochromaticity < 0.02 && reflectionRadius_1_per_A < 0.0005){
+ STATUS("Trying to index electron diffraction? It might be helpful to set a higher reflection radius (see documentation for --pinkIndexer-reflection-radius)")
+ }
+
+ float divergenceAngle_deg = 0.01; //fake
+
+ float tolerance = pinkIndexer_opts->tolerance;
+ Lattice_t sampleReciprocalLattice_1_per_A = lattice;
+ float detectorRadius_m = 0.03; //fake, only for prediction
+ ExperimentSettings* experimentSettings = ExperimentSettings_new(beamEenergy_eV, detectorDistance_m,
+ detectorRadius_m, divergenceAngle_deg, nonMonochromaticity, sampleReciprocalLattice_1_per_A, tolerance,
+ reflectionRadius_1_per_A);
+
+ consideredPeaksCount_t consideredPeaksCount = pinkIndexer_opts->considered_peaks_count;
+ angleResolution_t angleResolution = pinkIndexer_opts->angle_resolution;
+ refinementType_t refinementType = pinkIndexer_opts->refinement_type;
+ float maxResolutionForIndexing_1_per_A = pinkIndexer_opts->maxResolutionForIndexing_1_per_A;
+ pinkIndexer_private_data->pinkIndexer = PinkIndexer_new(experimentSettings, consideredPeaksCount, angleResolution,
+ refinementType,
+ maxResolutionForIndexing_1_per_A);
+
+ ExperimentSettings_delete(experimentSettings);
+ cell_free(primitiveCell);
+
+ /* Flags that pinkIndexer knows about */
+ *indm &= INDEXING_METHOD_MASK
+ | INDEXING_USE_CELL_PARAMETERS;
+
+ return pinkIndexer_private_data;
+}
+
+//static void reduceCell(UnitCell *cell, LatticeTransform_t* appliedReductionTransform)
+//{
+// double ax, ay, az, bx, by, bz, cx, cy, cz;
+// cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
+//
+// Lattice_t l = { ax, ay, az, bx, by, bz, cx, cy, cz };
+//
+// reduceLattice(&l, appliedReductionTransform);
+//
+// cell_set_cartesian(cell, l.ax, l.ay, l.az,
+// l.bx, l.by, l.bz,
+// l.cx, l.cy, l.cz);
+//
+// makeRightHanded(cell);
+//}
+//
+//static void restoreCell(UnitCell *cell, LatticeTransform_t* appliedReductionTransform)
+//{
+//
+// double ax, ay, az, bx, by, bz, cx, cy, cz;
+// cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
+//
+// Lattice_t l = { ax, ay, az, bx, by, bz, cx, cy, cz };
+//
+// restoreLattice(&l, appliedReductionTransform);
+//
+// cell_set_cartesian(cell, l.ax, l.ay, l.az,
+// l.bx, l.by, l.bz,
+// l.cx, l.cy, l.cz);
+//
+// makeRightHanded(cell);
+//}
+
+static void reduceReciprocalCell(UnitCell *cell, LatticeTransform_t* appliedReductionTransform)
+{
+ double ax, ay, az, bx, by, bz, cx, cy, cz;
+ cell_get_reciprocal(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
+
+ Lattice_t l = { ax, ay, az, bx, by, bz, cx, cy, cz };
+
+ reduceLattice(&l, appliedReductionTransform);
+
+ cell_set_reciprocal(cell, l.ax, l.ay, l.az,
+ l.bx, l.by, l.bz,
+ l.cx, l.cy, l.cz);
+
+ makeRightHanded(cell);
+}
+
+static void restoreReciprocalCell(UnitCell *cell, LatticeTransform_t* appliedReductionTransform)
+{
+
+ double ax, ay, az, bx, by, bz, cx, cy, cz;
+ cell_get_reciprocal(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
+
+ Lattice_t l = { ax, ay, az, bx, by, bz, cx, cy, cz };
+
+ restoreLattice(&l, appliedReductionTransform);
+
+ cell_set_reciprocal(cell, l.ax, l.ay, l.az,
+ l.bx, l.by, l.bz,
+ l.cx, l.cy, l.cz);
+
+ makeRightHanded(cell);
+}
+
+static void makeRightHanded(UnitCell *cell)
+{
+ double ax, ay, az, bx, by, bz, cx, cy, cz;
+ cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
+
+ if (!right_handed(cell)) {
+ cell_set_cartesian(cell, -ax, -ay, -az, -bx, -by, -bz, -cx, -cy, -cz);
+ }
+}
+
+//hack for electron crystallography while crystal_set_det_shift is not working approprietly
+static void update_detector(struct detector *det, double xoffs, double yoffs)
+{
+ int i;
+
+ for (i = 0; i < det->n_panels; i++) {
+ struct panel *p = &det->panels[i];
+ p->cnx += xoffs * p->res;
+ p->cny += yoffs * p->res;
+ }
+}
+
+void pinkIndexer_cleanup(void *pp)
+{
+ struct pinkIndexer_private_data* pinkIndexer_private_data = (struct pinkIndexer_private_data*) pp;
+
+ freeReciprocalPeaks(pinkIndexer_private_data->reciprocalPeaks_1_per_A);
+ free(pinkIndexer_private_data->intensities);
+ intmat_free(pinkIndexer_private_data->centeringTransformation);
+ PinkIndexer_delete(pinkIndexer_private_data->pinkIndexer);
+}
+
+const char *pinkIndexer_probe(UnitCell *cell)
+{
+ return "pinkIndexer";
+}
+
+#else /* HAVE_PINKINDEXER */
+
+int run_pinkIndexer(struct image *image, void *ipriv)
+{
+ ERROR("This copy of CrystFEL was compiled without PINKINDEXER support.\n");
+ return 0;
+}
+
+extern void *pinkIndexer_prepare(IndexingMethod *indm, UnitCell *cell,
+ struct pinkIndexer_options *pinkIndexer_opts,
+ const DataTemplate *dtempl)
+{
+ ERROR("This copy of CrystFEL was compiled without PINKINDEXER support.\n");
+ ERROR("To use PINKINDEXER indexing, recompile with PINKINDEXER.\n");
+ return NULL;
+}
+
+void pinkIndexer_cleanup(void *pp)
+{
+}
+
+const char *pinkIndexer_probe(UnitCell *cell)
+{
+ return NULL;
+}
+
+#endif /* HAVE_PINKINDEXER */
+
+static void pinkIndexer_show_help()
+{
+ printf(
+"Parameters for the PinkIndexer indexing algorithm:\n"
+" --pinkIndexer-considered-peaks-count=n\n"
+" Considered peaks count, 0 (fewest) to 4 (most)\n"
+" Default: 4\n"
+" --pinkIndexer-angle-resolution=n\n"
+" Angle resolution, 0 (loosest) to 4 (most dense)\n"
+" Default: 2\n"
+" --pinkIndexer-refinement-type=n\n"
+" Refinement type, 0 (none) to 5 (most accurate)\n"
+" Default: 1\n"
+" --pinkIndexer-tolerance=n\n"
+" Relative tolerance of the lattice vectors.\n"
+" Default 0.06\n"
+" --pinkIndexer-reflection-radius=n\n"
+" Radius of the reflections in reciprocal space.\n"
+" Specified in 1/A. Default is 2%% of a*.\n"
+" --pinkIndexer-max-resolution-for-indexing=n\n"
+" Measured in 1/A\n"
+" --pinkIndexer-multi Use pinkIndexers own multi indexing.\n"
+" --pinkIndexer-thread-count=n\n"
+" Thread count for internal parallelization \n"
+" Default: 1\n"
+" --pinkIndexer-no-check-indexed\n"
+" Disable internal check for correct indexing\n"
+" solutions\n"
+" --pinkIndexer-max-refinement-disbalance=n\n"
+" Maximum disbalance after refinement:\n"
+" 0 (no disbalance) to 2 (extreme disbalance), default 0.4\n"
+" --pinkIndexer-override-photon-energy=ev\n"
+" Mean energy in eV to use for indexing.\n"
+" --pinkIndexer-override-bandwidth=n\n"
+" Bandwidth in (delta energy)/(mean energy) to use for indexing.\n"
+" --pinkIndexer-override-visible-energy-range=min-max\n"
+" Overrides photon energy and bandwidth according to a range of \n"
+" energies that have high enough intensity to produce \"visible\" \n"
+" Bragg spots on the detector.\n"
+" Min and max range borders are separated by a minus sign (no whitespace).\n"
+ );
+}
+
+
+static error_t pinkindexer_parse_arg(int key, char *arg,
+ struct argp_state *state)
+{
+ float tmp, tmp2;
+ struct pinkIndexer_options **opts_ptr = state->input;
+
+ switch ( key ) {
+
+ case ARGP_KEY_INIT :
+ *opts_ptr = malloc(sizeof(struct pinkIndexer_options));
+ if ( *opts_ptr == NULL ) return ENOMEM;
+ (*opts_ptr)->considered_peaks_count = 4;
+ (*opts_ptr)->angle_resolution = 2;
+ (*opts_ptr)->refinement_type = 1;
+ (*opts_ptr)->tolerance = 0.06;
+ (*opts_ptr)->maxResolutionForIndexing_1_per_A = +INFINITY;
+ (*opts_ptr)->thread_count = 1;
+ (*opts_ptr)->multi = 0;
+ (*opts_ptr)->no_check_indexed = 0;
+ (*opts_ptr)->min_peaks = 2;
+ (*opts_ptr)->reflectionRadius = -1;
+ (*opts_ptr)->customPhotonEnergy = -1;
+ (*opts_ptr)->customBandwidth = -1;
+ (*opts_ptr)->maxRefinementDisbalance = 0.4;
+ break;
+
+ case 1 :
+ pinkIndexer_show_help();
+ return EINVAL;
+
+ case 2 :
+ if (sscanf(arg, "%u", &(*opts_ptr)->considered_peaks_count) != 1)
+ {
+ ERROR("Invalid value for "
+ "--pinkIndexer-considered-peaks-count\n");
+ return EINVAL;
+ }
+ break;
+
+ case 3 :
+ if (sscanf(arg, "%u", &(*opts_ptr)->angle_resolution) != 1)
+ {
+ ERROR("Invalid value for "
+ "--pinkIndexer-angle_resolution\n");
+ return EINVAL;
+ }
+ break;
+
+ case 4 :
+ if (sscanf(arg, "%u", &(*opts_ptr)->refinement_type) != 1)
+ {
+ ERROR("Invalid value for "
+ "--pinkIndexer-refinement-type\n");
+ return EINVAL;
+ }
+ break;
+
+ case 5 :
+ if (sscanf(arg, "%d", &(*opts_ptr)->thread_count) != 1)
+ {
+ ERROR("Invalid value for --pinkIndexer-thread-count\n");
+ return EINVAL;
+ }
+ break;
+
+ case 6 :
+ if (sscanf(arg, "%f", &(*opts_ptr)->maxResolutionForIndexing_1_per_A) != 1)
+ {
+ ERROR("Invalid value for "
+ "--pinkIndexer-max-resolution-for-indexing\n");
+ return EINVAL;
+ }
+ break;
+
+ case 7 :
+ if (sscanf(arg, "%f", &(*opts_ptr)->tolerance) != 1)
+ {
+ ERROR("Invalid value for --pinkIndexer-tolerance\n");
+ return EINVAL;
+ }
+ break;
+
+ case 8 :
+ (*opts_ptr)->multi = 1;
+ break;
+
+ case 9 :
+ (*opts_ptr)->no_check_indexed = 1;
+ break;
+
+ case 10 :
+ if (sscanf(arg, "%f", &tmp) != 1) {
+ ERROR("Invalid value for --pinkIndexer-reflection-radius\n");
+ return EINVAL;
+ }
+ (*opts_ptr)->reflectionRadius = tmp / 1e10; /* A^-1 to m^-1 */
+ break;
+
+ case 11 :
+ if (sscanf(arg, "%f", &(*opts_ptr)->customPhotonEnergy) != 1)
+ {
+ ERROR("Invalid value for --pinkIndexer-override-photon-energy\n");
+ return EINVAL;
+ }
+ break;
+
+ case 12 :
+ if (sscanf(arg, "%f", &(*opts_ptr)->customBandwidth) != 1)
+ {
+ ERROR("Invalid value for --pinkIndexer-override-bandwidth\n");
+ return EINVAL;
+ }
+ break;
+ case 13 :
+ if (sscanf(arg, "%f-%f", &tmp, &tmp2) != 2)
+ {
+ ERROR("Invalid value for --pinkIndexer-override-visible-energy-range\n");
+ return EINVAL;
+ }
+ (*opts_ptr)->customPhotonEnergy = (tmp + tmp2)/2;
+ (*opts_ptr)->customBandwidth = (tmp2 - tmp)/(*opts_ptr)->customPhotonEnergy;
+ if((*opts_ptr)->customBandwidth < 0){
+ (*opts_ptr)->customBandwidth *= -1;
+ }
+ break;
+ case 14 :
+ if (sscanf(arg, "%f", &(*opts_ptr)->maxRefinementDisbalance) != 1)
+ {
+ ERROR("Invalid value for --pinkIndexer-max-refinement-disbalance\n");
+ return EINVAL;
+ }
+ }
+
+ return 0;
+}
+
+
+static struct argp_option pinkindexer_options[] = {
+
+ {"help-pinkindexer", 1, NULL, OPTION_NO_USAGE,
+ "Show options for PinkIndexer indexing algorithm", 99},
+
+ {"pinkIndexer-considered-peaks-count", 2, "n", OPTION_HIDDEN, NULL},
+ {"pinkIndexer-cpc", 2, "n", OPTION_HIDDEN, NULL},
+
+ {"pinkIndexer-angle-resolution", 3, "ang", OPTION_HIDDEN, NULL},
+ {"pinkIndexer-ar", 3, "ang", OPTION_HIDDEN, NULL},
+
+ {"pinkIndexer-refinement-type", 4, "t", OPTION_HIDDEN, NULL},
+ {"pinkIndexer-rt", 4, "t", OPTION_HIDDEN, NULL},
+
+ {"pinkIndexer-thread-count", 5, "n", OPTION_HIDDEN, NULL},
+ {"pinkIndexer-tc", 5, "n", OPTION_HIDDEN, NULL},
+
+ {"pinkIndexer-max-resolution-for-indexing", 6, "res", OPTION_HIDDEN, NULL},
+ {"pinkIndexer-mrfi", 6, "res", OPTION_HIDDEN, NULL},
+
+ {"pinkIndexer-tolerance", 7, "tol", OPTION_HIDDEN, NULL},
+ {"pinkIndexer-tol", 7, "tol", OPTION_HIDDEN, NULL},
+
+ {"pinkIndexer-multi", 8, NULL, OPTION_HIDDEN, NULL},
+
+ {"pinkIndexer-no-check-indexed", 9, NULL, OPTION_HIDDEN, NULL},
+
+ {"pinkIndexer-reflection-radius", 10, "r", OPTION_HIDDEN, NULL},
+
+ {"pinkIndexer-override-photon-energy", 11, "ev", OPTION_HIDDEN, NULL},
+
+ {"pinkIndexer-override-bandwidth", 12, "bw", OPTION_HIDDEN, NULL},
+
+ {"pinkIndexer-override-visible-energy-range", 13, "overridenVisibleEnergyRange", OPTION_HIDDEN, NULL},
+
+ {"pinkIndexer-max-refinement-disbalance", 14, "maxDisbalance", OPTION_HIDDEN, NULL},
+
+ {0}
+};
+
+
+struct argp pinkIndexer_argp = { pinkindexer_options,
+ pinkindexer_parse_arg,
+ NULL, NULL, NULL, NULL, NULL };
diff --git a/libcrystfel/src/indexers/pinkindexer.h b/libcrystfel/src/indexers/pinkindexer.h
new file mode 100644
index 00000000..cab88a75
--- /dev/null
+++ b/libcrystfel/src/indexers/pinkindexer.h
@@ -0,0 +1,47 @@
+/*
+ * pinkindexer.h
+ *
+ * Interface to PinkIndexer
+ *
+ * Copyright © 2017-2020 Deutsches Elektronen-Synchrotron DESY,
+ * a research centre of the Helmholtz Association.
+ *
+ * Authors:
+ * 2017-2019 Yaroslav Gevorkov <yaroslav.gevorkov@desy.de>
+ *
+ * This file is part of CrystFEL.
+ *
+ * CrystFEL is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * CrystFEL is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with CrystFEL. If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+#ifndef LIBCRYSTFEL_SRC_PINKINDEXER_H_
+#define LIBCRYSTFEL_SRC_PINKINDEXER_H_
+
+#include <stddef.h>
+
+#include "index.h"
+#include "datatemplate.h"
+
+extern int run_pinkIndexer(struct image *image, void *ipriv);
+
+extern void *pinkIndexer_prepare(IndexingMethod *indm, UnitCell *cell,
+ struct pinkIndexer_options *pinkIndexer_opts,
+ const DataTemplate *dtempl);
+
+extern void pinkIndexer_cleanup(void *pp);
+
+extern const char *pinkIndexer_probe(UnitCell *cell);
+
+#endif /* LIBCRYSTFEL_SRC_PINKINDEXER_H_ */
diff --git a/libcrystfel/src/indexers/taketwo.c b/libcrystfel/src/indexers/taketwo.c
new file mode 100644
index 00000000..65265b0e
--- /dev/null
+++ b/libcrystfel/src/indexers/taketwo.c
@@ -0,0 +1,2370 @@
+/*
+ * taketwo.c
+ *
+ * Rewrite of TakeTwo algorithm (Acta D72 (8) 956-965) for CrystFEL
+ *
+ * Copyright © 2016-2017 Helen Ginn
+ * Copyright © 2016-2020 Deutsches Elektronen-Synchrotron DESY,
+ * a research centre of the Helmholtz Association.
+ *
+ * Authors:
+ * 2016-2017 Helen Ginn <helen@strubi.ox.ac.uk>
+ * 2016-2017 Thomas White <taw@physics.org>
+ *
+ * This file is part of CrystFEL.
+ *
+ * CrystFEL is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * CrystFEL is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with CrystFEL. If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+/**
+ * \file taketwo.h
+
+ * ## Code outline
+
+ * ### Get ready for calculation
+ * * Pre-calculate symmetry operations (generate_rotation_symops())
+ * * Pre-calculate theoretical vectors from unit cell dimensions
+ * (gen_theoretical_vecs())
+ * * Generate observed vectors from data (gen_observed_vecs())
+ * * Match observed vectors to theoretical vectors (match_obs_to_cell_vecs())
+ *
+ * ### Business bit
+ *
+ * ... n.b. rearranging to find all seeds in advance.
+ *
+ * * Find starting seeds (find_seeds()):
+ * - Loop through pairs of observed vectors
+ * - If they share a spot, find matching pairs of theoretical vectors
+ * - Remove all duplicate matches due to symmetry operations
+ * - For the remainder, loop through the matches and extend the seeds
+ * (start_seed()).
+ * - If it returns a membership greater than the highest member threshold,
+ * return the matrix to CrystFEL.
+ *
+ * * Extending a seed (start_seed()):
+ * - Generate a rotation matrix which matches the chosen start seed.
+ * - Loop through all observed vectors starting from 0.
+ * - Find another vector to add to the network, if available
+ * (find_next_index()).
+ * - If the index is not available, then give up if there were too many dead
+ * ends. Otherwise, remove the last member and pretend like that didn't
+ * happen, so it won't happen again.
+ * - Add the vector to increment the membership list.
+ * - If the membership number exceeds the maximum, tidy up the solution and
+ * return a success.
+ * - If the membership does not, then resume the loop and search for the
+ * next vector.
+ *
+ * * Finding the next member (find_next_index()):
+ * - Go through the observed vectors, starting from the last index + 1 to
+ * explore only the "new" vectors.
+ * - If the vector does not share a spot with the current array of vectors,
+ * then skip it.
+ * - We must loop through all the current vectors in the network, and see if
+ * they match the newcomer for a given matching theoretical vector.
+ * - We only accept a match if the rotation matrix matches the seed matrix
+ * for a single matching theoretical vector.
+ * - If it does match, we can return a success.
+ *
+ * * Tidying the solution (finish_solution()):
+ * - This chooses the most common rotation matrix of the bunch to choose to
+ * send to CrystFEL. But this should probably take the average instead,
+ * which is very possible.
+ *
+ * * Clean up the mess (cleanup_taketwo_obs_vecs())
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <gsl/gsl_matrix.h>
+#include <gsl/gsl_blas.h>
+#include <float.h>
+#include <math.h>
+#include <assert.h>
+#include <time.h>
+
+#include "cell.h"
+#include "cell-utils.h"
+#include "index.h"
+#include "taketwo.h"
+#include "peaks.h"
+#include "symmetry.h"
+
+struct taketwo_options
+{
+ int member_thresh;
+ double len_tol;
+ double angle_tol;
+ double trace_tol;
+};
+
+/**
+ * \param obsvec an observed vector between two spots
+ * \param matches array of matching theoretical vectors from unit cell
+ * \param match_num number of matches
+ * \param distance length of obsvec (do I need this?)
+ * \param her_rlp pointer to first rlp position for difference vec
+ * \param his_rlp pointer to second rlp position for difference vec
+ *
+ * Structure representing 3D vector between two potential Bragg peaks
+ * in reciprocal space, and an array of potential matching theoretical
+ * vectors from unit cell/centering considerations.
+ **/
+struct SpotVec
+{
+ struct rvec obsvec;
+ struct TheoryVec *matches;
+ int match_num;
+ int assignment;
+ int in_network;
+ double distance;
+ struct rvec *her_rlp;
+ struct rvec *his_rlp;
+};
+
+/**
+ * theoryvec
+ */
+
+struct TheoryVec
+{
+ struct rvec vec;
+ int asym;
+};
+
+
+/**
+ * seed
+ */
+
+struct Seed
+{
+ int obs1;
+ int obs2;
+ int idx1;
+ int idx2;
+ double score;
+};
+
+struct taketwo_private
+{
+ IndexingMethod indm;
+ struct taketwo_options *opts;
+ UnitCell *cell;
+ int serial_num; /**< Serial of last image, -1 if unassigned */
+ unsigned int xtal_num; /**< last number of crystals recorded */
+
+ struct TheoryVec *theory_vecs; /**< Theoretical vectors for given unit cell */
+ unsigned int vec_count; /**< Number of theoretical vectors */
+
+ gsl_matrix **prevSols; /**< Previous solutions to be ignored */
+ unsigned int numPrevs; /**< Previous solution count */
+ double *prevScores; /**< previous solution scores */
+ unsigned int *membership; /**< previous solution was success or failure */
+
+};
+
+/**
+ * Internal structure which gets passed the various functions looking after
+ * the indexing bits and bobs. */
+struct TakeTwoCell
+{
+ UnitCell *cell; /**< Contains unit cell dimensions */
+ gsl_matrix **rotSymOps;
+ unsigned int numOps;
+
+ struct SpotVec *obs_vecs;
+ struct Seed *seeds;
+ int seed_count;
+ int obs_vec_count;
+
+ /* Options */
+ int member_thresh;
+ double len_tol; /**< In reciprocal metres */
+ double angle_tol; /**< In radians */
+ double trace_tol; /**< Contains sqrt(4*(1-cos(angle))) */
+
+ /** A given solution to refine */
+ gsl_matrix *solution;
+
+ double x_ang; /**< Rotations in radians to apply to x axis of solution */
+ double y_ang; /**< Rotations in radians to apply to y axis of solution */
+ double z_ang; /**< Rotations in radians to apply to z axis of solution */
+
+ /**< Temporary memory always allocated for calculations */
+ gsl_matrix *twiz1Tmp;
+ /**< Temporary memory always allocated for calculations */
+ gsl_matrix *twiz2Tmp;
+ /**< Temporary memory always allocated for calculations */
+ gsl_vector *vec1Tmp;
+ /**< Temporary memory always allocated for calculations */
+ gsl_vector *vec2Tmp;
+};
+
+
+/* Maximum distance between two rlp sizes to consider info for indexing */
+#define MAX_RECIP_DISTANCE (0.15*1e10)
+
+/* Tolerance for two lengths in reciprocal space to be considered the same */
+#define RECIP_TOLERANCE (0.0010*1e10)
+
+/* Threshold for network members to consider a potential solution */
+#define NETWORK_MEMBER_THRESHOLD (20)
+
+/* Minimum for network members to consider a potential solution */
+#define MINIMUM_MEMBER_THRESHOLD (5)
+
+/* Maximum dead ends for a single branch extension during indexing */
+#define MAX_DEAD_ENDS (10)
+
+/* Maximum observed vectors before TakeTwo gives up and deals with
+ * what is already there. */
+#define MAX_OBS_VECTORS 100000
+
+/* Tolerance for two angles to be considered the same */
+#define ANGLE_TOLERANCE (deg2rad(0.6))
+
+/* Tolerance for rot_mats_are_similar */
+#define TRACE_TOLERANCE (deg2rad(3.0))
+
+/* Initial step size for refinement of solutions */
+#define ANGLE_STEP_SIZE (deg2rad(0.5))
+
+/* Final required step size for refinement of solutions */
+#define ANGLE_CONVERGE_SIZE (deg2rad(0.01))
+
+/* TODO: Multiple lattices */
+
+
+/* ------------------------------------------------------------------------
+ * apologetic function
+ * ------------------------------------------------------------------------*/
+
+void apologise()
+{
+ printf("Error - could not allocate memory for indexing.\n");
+}
+
+/* ------------------------------------------------------------------------
+ * functions concerning aspects of rvec which are very likely to be
+ * incorporated somewhere else in CrystFEL and therefore may need to be
+ * deleted and references connected to a pre-existing function. (Lowest level)
+ * ------------------------------------------------------------------------*/
+
+static struct rvec new_rvec(double new_u, double new_v, double new_w)
+{
+ struct rvec new_rvector;
+ new_rvector.u = new_u;
+ new_rvector.v = new_v;
+ new_rvector.w = new_w;
+
+ return new_rvector;
+}
+
+static struct rvec rvec_add_rvec(struct rvec first, struct rvec second)
+{
+ struct rvec diff = new_rvec(second.u + first.u,
+ second.v + first.v,
+ second.w + first.w);
+
+ return diff;
+}
+
+
+static struct rvec diff_vec(struct rvec from, struct rvec to)
+{
+ struct rvec diff = new_rvec(to.u - from.u,
+ to.v - from.v,
+ to.w - from.w);
+
+ return diff;
+}
+
+static double sq_length(struct rvec vec)
+{
+ double sqlength = (vec.u * vec.u + vec.v * vec.v + vec.w * vec.w);
+
+ return sqlength;
+}
+
+
+static double rvec_length(struct rvec vec)
+{
+ return sqrt(sq_length(vec));
+}
+
+
+static void normalise_rvec(struct rvec *vec)
+{
+ double length = rvec_length(*vec);
+ vec->u /= length;
+ vec->v /= length;
+ vec->w /= length;
+}
+
+
+static double rvec_cosine(struct rvec v1, struct rvec v2)
+{
+ double dot_prod = v1.u * v2.u + v1.v * v2.v + v1.w * v2.w;
+ double v1_length = rvec_length(v1);
+ double v2_length = rvec_length(v2);
+
+ double cos_theta = dot_prod / (v1_length * v2_length);
+
+ return cos_theta;
+}
+
+
+static double rvec_angle(struct rvec v1, struct rvec v2)
+{
+ double cos_theta = rvec_cosine(v1, v2);
+ double angle = acos(cos_theta);
+
+ return angle;
+}
+
+
+static struct rvec rvec_cross(struct rvec a, struct rvec b)
+{
+ struct rvec c;
+
+ c.u = a.v*b.w - a.w*b.v;
+ c.v = -(a.u*b.w - a.w*b.u);
+ c.w = a.u*b.v - a.v*b.u;
+
+ return c;
+}
+
+/*
+static void show_rvec(struct rvec r2)
+{
+ struct rvec r = r2;
+ normalise_rvec(&r);
+ STATUS("[ %.3f %.3f %.3f ]\n", r.u, r.v, r.w);
+}
+*/
+
+
+/* ------------------------------------------------------------------------
+ * functions called under the core functions, still specialised (Level 3)
+ * ------------------------------------------------------------------------*/
+
+/* cell_transform_gsl_direct() doesn't do quite what we want here.
+ * The matrix m should be post-multiplied by a matrix of real or reciprocal
+ * basis vectors (it doesn't matter which because it's just a rotation).
+ * M contains the basis vectors written in columns: M' = mM */
+static UnitCell *cell_post_smiley_face(UnitCell *in, gsl_matrix *m)
+{
+ gsl_matrix *c;
+ double asx, asy, asz;
+ double bsx, bsy, bsz;
+ double csx, csy, csz;
+ gsl_matrix *res;
+ UnitCell *out;
+
+ cell_get_cartesian(in, &asx, &asy, &asz,
+ &bsx, &bsy, &bsz,
+ &csx, &csy, &csz);
+
+ c = gsl_matrix_alloc(3, 3);
+ gsl_matrix_set(c, 0, 0, asx);
+ 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, 2, 1, bsz);
+ gsl_matrix_set(c, 0, 2, csx);
+ gsl_matrix_set(c, 1, 2, csy);
+ gsl_matrix_set(c, 2, 2, csz);
+
+ res = gsl_matrix_calloc(3, 3);
+ gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, m, c, 0.0, res);
+
+ out = cell_new_from_cell(in);
+ cell_set_cartesian(out, gsl_matrix_get(res, 0, 0),
+ gsl_matrix_get(res, 1, 0),
+ 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);
+ gsl_matrix_free(c);
+ return out;
+}
+
+
+static void rotation_around_axis(struct rvec c, double th,
+ gsl_matrix *res)
+{
+ double omc = 1.0 - cos(th);
+ double s = sin(th);
+ gsl_matrix_set(res, 0, 0, cos(th) + c.u*c.u*omc);
+ gsl_matrix_set(res, 0, 1, c.u*c.v*omc - c.w*s);
+ gsl_matrix_set(res, 0, 2, c.u*c.w*omc + c.v*s);
+ gsl_matrix_set(res, 1, 0, c.u*c.v*omc + c.w*s);
+ gsl_matrix_set(res, 1, 1, cos(th) + c.v*c.v*omc);
+ gsl_matrix_set(res, 1, 2, c.v*c.w*omc - c.u*s);
+ gsl_matrix_set(res, 2, 0, c.w*c.u*omc - c.v*s);
+ gsl_matrix_set(res, 2, 1, c.w*c.v*omc + c.u*s);
+ gsl_matrix_set(res, 2, 2, cos(th) + c.w*c.w*omc);
+}
+
+/** Rotate GSL matrix by three angles along x, y and z axes */
+static void rotate_gsl_by_angles(gsl_matrix *sol, double x, double y,
+ double z, gsl_matrix *result)
+{
+ gsl_matrix *x_rot = gsl_matrix_alloc(3, 3);
+ gsl_matrix *y_rot = gsl_matrix_alloc(3, 3);
+ gsl_matrix *z_rot = gsl_matrix_alloc(3, 3);
+ gsl_matrix *xy_rot = gsl_matrix_alloc(3, 3);
+ gsl_matrix *xyz_rot = gsl_matrix_alloc(3, 3);
+
+ struct rvec x_axis = new_rvec(1, 0, 0);
+ struct rvec y_axis = new_rvec(1, 0, 0);
+ struct rvec z_axis = new_rvec(1, 0, 0);
+
+ rotation_around_axis(x_axis, x, x_rot);
+ rotation_around_axis(y_axis, y, y_rot);
+ rotation_around_axis(z_axis, z, z_rot);
+
+ /* Collapse the rotations in x and y onto z */
+ gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, x_rot,
+ y_rot, 0.0, xy_rot);
+ gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, xy_rot,
+ z_rot, 0.0, xyz_rot);
+
+ /* Apply the whole rotation offset to the solution */
+ gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, xyz_rot,
+ sol, 0.0, result);
+
+ gsl_matrix_free(x_rot);
+ gsl_matrix_free(y_rot);
+ gsl_matrix_free(z_rot);
+ gsl_matrix_free(xy_rot);
+ gsl_matrix_free(xyz_rot);
+}
+
+
+/* Rotate vector (vec1) around axis (axis) by angle theta. Find value of
+ * theta for which the angle between (vec1) and (vec2) is minimised. */
+static void closest_rot_mat(struct rvec vec1, struct rvec vec2,
+ struct rvec axis, gsl_matrix *twizzle)
+{
+ /* Let's have unit vectors */
+ normalise_rvec(&vec1);
+ normalise_rvec(&vec2);
+ normalise_rvec(&axis);
+
+ /* Redeclaring these to try and maintain readability and
+ * check-ability against the maths I wrote down */
+ double a = vec2.u; double b = vec2.v; double c = vec2.w;
+ double p = vec1.u; double q = vec1.v; double r = vec1.w;
+ double x = axis.u; double y = axis.v; double z = axis.w;
+
+ /* Components in handwritten maths online when I upload it */
+ double A = a*(p*x*x - p + x*y*q + x*z*r) +
+ b*(p*x*y + q*y*y - q + r*y*z) +
+ c*(p*x*z + q*y*z + r*z*z - r);
+
+ double B = a*(y*r - z*q) + b*(p*z - r*x) + c*(q*x - p*y);
+
+ double tan_theta = - B / A;
+ double theta = atan(tan_theta);
+
+ /* Now we have two possible solutions, theta or theta+pi
+ * and we need to work out which one. This could potentially be
+ * simplified - do we really need so many cos/sins? maybe check
+ * the 2nd derivative instead? */
+ double cc = cos(theta);
+ double C = 1 - cc;
+ double s = sin(theta);
+ double occ = -cc;
+ double oC = 1 - occ;
+ double os = -s;
+
+ double pPrime = (x*x*C+cc)*p + (x*y*C-z*s)*q + (x*z*C+y*s)*r;
+ double qPrime = (y*x*C+z*s)*p + (y*y*C+cc)*q + (y*z*C-x*s)*r;
+ double rPrime = (z*x*C-y*s)*p + (z*y*C+x*s)*q + (z*z*C+cc)*r;
+
+ double pDbPrime = (x*x*oC+occ)*p + (x*y*oC-z*os)*q + (x*z*oC+y*os)*r;
+ double qDbPrime = (y*x*oC+z*os)*p + (y*y*oC+occ)*q + (y*z*oC-x*os)*r;
+ double rDbPrime = (z*x*oC-y*os)*p + (z*y*oC+x*os)*q + (z*z*oC+occ)*r;
+
+ double cosAlpha = pPrime * a + qPrime * b + rPrime * c;
+ double cosAlphaOther = pDbPrime * a + qDbPrime * b + rDbPrime * c;
+
+ int addPi = (cosAlphaOther > cosAlpha);
+ double bestAngle = theta + addPi * M_PI;
+
+ /* Don't return an identity matrix which has been rotated by
+ * theta around "axis", but do assign it to twizzle. */
+ rotation_around_axis(axis, bestAngle, twizzle);
+}
+
+static double matrix_trace(gsl_matrix *a)
+{
+ int i;
+ double tr = 0.0;
+
+ assert(a->size1 == a->size2);
+ for ( i=0; i<a->size1; i++ ) {
+ tr += gsl_matrix_get(a, i, i);
+ }
+ return tr;
+}
+
+static char *add_ua(const char *inp, char ua)
+{
+ char *pg = malloc(64);
+ if ( pg == NULL ) return NULL;
+ snprintf(pg, 63, "%s_ua%c", inp, ua);
+ return pg;
+}
+
+
+static char *get_chiral_holohedry(UnitCell *cell)
+{
+ LatticeType lattice = cell_get_lattice_type(cell);
+ char *pg;
+ char *pgout = 0;
+
+ switch (lattice)
+ {
+ case L_TRICLINIC:
+ pg = "1";
+ break;
+
+ case L_MONOCLINIC:
+ pg = "2";
+ break;
+
+ case L_ORTHORHOMBIC:
+ pg = "222";
+ break;
+
+ case L_TETRAGONAL:
+ pg = "422";
+ break;
+
+ case L_RHOMBOHEDRAL:
+ pg = "3_R";
+ break;
+
+ case L_HEXAGONAL:
+ if ( cell_get_centering(cell) == 'H' ) {
+ pg = "3_H";
+ } else {
+ pg = "622";
+ }
+ break;
+
+ case L_CUBIC:
+ pg = "432";
+ break;
+
+ default:
+ pg = "error";
+ break;
+ }
+
+ switch (lattice)
+ {
+ case L_TRICLINIC:
+ case L_ORTHORHOMBIC:
+ case L_RHOMBOHEDRAL:
+ case L_CUBIC:
+ pgout = strdup(pg);
+ break;
+
+ case L_MONOCLINIC:
+ case L_TETRAGONAL:
+ case L_HEXAGONAL:
+ pgout = add_ua(pg, cell_get_unique_axis(cell));
+ break;
+
+ default:
+ break;
+ }
+
+ return pgout;
+}
+
+
+static SymOpList *sym_ops_for_cell(UnitCell *cell)
+{
+ SymOpList *rawList;
+
+ char *pg = get_chiral_holohedry(cell);
+ rawList = get_pointgroup(pg);
+ free(pg);
+
+ return rawList;
+}
+
+static int rot_mats_are_similar(gsl_matrix *rot1, gsl_matrix *rot2,
+ gsl_matrix *sub, gsl_matrix *mul,
+ double *score, struct TakeTwoCell *cell)
+{
+ double tr;
+
+ gsl_matrix_memcpy(sub, rot1);
+ gsl_matrix_sub(sub, rot2); /* sub = rot1 - rot2 */
+
+ gsl_blas_dgemm(CblasNoTrans, CblasTrans, 1.0, sub, sub, 0.0, mul);
+
+ tr = matrix_trace(mul);
+ if (score != NULL) *score = tr;
+
+ return (tr < cell->trace_tol);
+}
+
+static int symm_rot_mats_are_similar(gsl_matrix *rot1, gsl_matrix *rot2,
+ struct TakeTwoCell *cell)
+{
+ int i;
+
+ gsl_matrix *sub = gsl_matrix_calloc(3, 3);
+ gsl_matrix *mul = gsl_matrix_calloc(3, 3);
+
+ for (i = 0; i < cell->numOps; i++) {
+ gsl_matrix *testRot = gsl_matrix_alloc(3, 3);
+ gsl_matrix *symOp = cell->rotSymOps[i];
+
+ gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, rot1, symOp,
+ 0.0, testRot);
+
+ if (rot_mats_are_similar(testRot, rot2, sub, mul, NULL, cell)) {
+ gsl_matrix_free(testRot);
+ gsl_matrix_free(sub);
+ gsl_matrix_free(mul);
+ return 1;
+ }
+
+ gsl_matrix_free(testRot);
+ }
+
+ gsl_matrix_free(sub);
+ gsl_matrix_free(mul);
+
+ return 0;
+}
+
+static void rotation_between_vectors(struct rvec a, struct rvec b,
+ gsl_matrix *twizzle)
+{
+ double th = rvec_angle(a, b);
+ struct rvec c = rvec_cross(a, b);
+ normalise_rvec(&c);
+ rotation_around_axis(c, th, twizzle);
+}
+
+
+static void rvec_to_gsl(gsl_vector *newVec, struct rvec v)
+{
+ gsl_vector_set(newVec, 0, v.u);
+ gsl_vector_set(newVec, 1, v.v);
+ gsl_vector_set(newVec, 2, v.w);
+}
+
+
+struct rvec gsl_to_rvec(gsl_vector *a)
+{
+ struct rvec v;
+ v.u = gsl_vector_get(a, 0);
+ v.v = gsl_vector_get(a, 1);
+ v.w = gsl_vector_get(a, 2);
+ return v;
+}
+
+
+/* Code me in gsl_matrix language to copy the contents of the function
+ * in cppxfel (IndexingSolution::createSolution). This function is quite
+ * intensive on the number crunching side so simple angle checks are used
+ * to 'pre-scan' vectors beforehand. */
+static gsl_matrix *generate_rot_mat(struct rvec obs1, struct rvec obs2,
+ struct rvec cell1, struct rvec cell2,
+ struct TakeTwoCell *cell)
+{
+ gsl_matrix *fullMat;
+ rvec_to_gsl(cell->vec1Tmp, cell2);
+
+ normalise_rvec(&obs1);
+ normalise_rvec(&obs2);
+ normalise_rvec(&cell1);
+ normalise_rvec(&cell2);
+
+ /* Rotate reciprocal space so that the first simulated vector lines up
+ * with the observed vector. */
+ rotation_between_vectors(cell1, obs1, cell->twiz1Tmp);
+
+ normalise_rvec(&obs1);
+
+ /* Multiply cell2 by rotateSpotDiffMatrix --> cell2vr */
+ gsl_blas_dgemv(CblasNoTrans, 1.0, cell->twiz1Tmp, cell->vec1Tmp,
+ 0.0, cell->vec2Tmp);
+
+ /* Now we twirl around the firstAxisUnit until the rotated simulated
+ * vector matches the second observed vector as closely as possible. */
+
+ closest_rot_mat(gsl_to_rvec(cell->vec2Tmp), obs2, obs1, cell->twiz2Tmp);
+
+ /* 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(CblasTrans, CblasTrans, 1.0,
+ cell->twiz1Tmp, cell->twiz2Tmp, 0.0, fullMat);
+ gsl_matrix_transpose(fullMat);
+
+ return fullMat;
+}
+
+
+static int obs_vecs_share_spot(struct SpotVec *her_obs, struct SpotVec *his_obs)
+{
+ if ( (her_obs->her_rlp == his_obs->her_rlp) ||
+ (her_obs->her_rlp == his_obs->his_rlp) ||
+ (her_obs->his_rlp == his_obs->her_rlp) ||
+ (her_obs->his_rlp == his_obs->his_rlp) ) {
+ return 1;
+ }
+
+ return 0;
+}
+
+
+static int obs_shares_spot_w_array(struct SpotVec *obs_vecs, int test_idx,
+ int *members, int num)
+{
+ int i;
+
+ struct SpotVec *her_obs = &obs_vecs[test_idx];
+
+ for ( i=0; i<num; i++ ) {
+ struct SpotVec *his_obs = &obs_vecs[members[i]];
+
+ int shares = obs_vecs_share_spot(her_obs, his_obs);
+
+ if ( shares ) return 1;
+ }
+
+ return 0;
+}
+
+
+static int obs_vecs_match_angles(int her, int his,
+ struct Seed **seeds, int *match_count,
+ struct TakeTwoCell *cell)
+{
+ struct SpotVec *obs_vecs = cell->obs_vecs;
+ struct SpotVec *her_obs = &obs_vecs[her];
+ struct SpotVec *his_obs = &obs_vecs[his];
+
+ *match_count = 0;
+
+ double min_angle = deg2rad(2.5);
+ double max_angle = deg2rad(187.5);
+
+ /* calculate angle between observed vectors */
+ double obs_angle = rvec_angle(her_obs->obsvec, his_obs->obsvec);
+
+ /* calculate angle between all potential theoretical vectors */
+
+ int i, j;
+ for ( i=0; i<her_obs->match_num; i++ ) {
+ for ( j=0; j<his_obs->match_num; j++ ) {
+ double score = 0;
+
+ struct rvec *her_match = &her_obs->matches[i].vec;
+ struct rvec *his_match = &his_obs->matches[j].vec;
+
+ double her_dist = rvec_length(*her_match);
+ double his_dist = rvec_length(*his_match);
+
+ double theory_angle = rvec_angle(*her_match,
+ *his_match);
+
+ /* is this angle a match? */
+
+ double angle_diff = fabs(theory_angle - obs_angle);
+
+ /* within basic tolerance? */
+ if ( angle_diff != angle_diff ||
+ angle_diff > cell->angle_tol ) {
+ continue;
+ }
+
+ double add = angle_diff;
+ if (add == add) {
+ score += add * her_dist * his_dist;
+ }
+
+ /* If the angles are too close to 0
+ or 180, one axis ill-determined */
+ if (theory_angle < min_angle ||
+ theory_angle > max_angle) {
+ continue;
+ }
+
+ /* check that third vector adequately completes
+ * triangle */
+
+ struct rvec theory_diff = diff_vec(*his_match, *her_match);
+ struct rvec obs_diff = diff_vec(his_obs->obsvec,
+ her_obs->obsvec);
+
+ theory_angle = rvec_angle(*her_match,
+ theory_diff);
+ obs_angle = rvec_angle(her_obs->obsvec, obs_diff);
+ angle_diff = fabs(obs_angle - theory_angle);
+
+ double diff_dist = rvec_length(obs_diff);
+
+ if (angle_diff > ANGLE_TOLERANCE) {
+ continue;
+ }
+
+ add = angle_diff;
+ if (add == add) {
+ score += add * her_dist * diff_dist;
+ }
+
+ theory_angle = rvec_angle(*his_match,
+ theory_diff);
+ obs_angle = rvec_angle(his_obs->obsvec, obs_diff);
+
+ if (fabs(obs_angle - theory_angle) > ANGLE_TOLERANCE) {
+ continue;
+ }
+
+ add = angle_diff;
+ if (add == add) {
+ score += add * his_dist * diff_dist;
+ }
+
+ /* we add a new seed to the array */
+ size_t new_size = (*match_count + 1);
+ new_size *= sizeof(struct Seed);
+
+ /* Reallocate the array to fit in another match */
+ struct Seed *tmp_seeds = realloc(*seeds, new_size);
+
+ if ( tmp_seeds == NULL ) {
+ apologise();
+ }
+
+ (*seeds) = tmp_seeds;
+
+ (*seeds)[*match_count].obs1 = her;
+ (*seeds)[*match_count].obs2 = his;
+ (*seeds)[*match_count].idx1 = i;
+ (*seeds)[*match_count].idx2 = j;
+ (*seeds)[*match_count].score = score * 1000;
+
+ (*match_count)++;
+ }
+ }
+
+ return (*match_count > 0);
+}
+
+/* ------------------------------------------------------------------------
+ * core functions regarding the meat of the TakeTwo algorithm (Level 2)
+ * ------------------------------------------------------------------------*/
+
+static signed int finish_solution(gsl_matrix *rot, struct SpotVec *obs_vecs,
+ int *obs_members, int *match_members,
+ int member_num, struct TakeTwoCell *cell)
+{
+ gsl_matrix *sub = gsl_matrix_calloc(3, 3);
+ gsl_matrix *mul = gsl_matrix_calloc(3, 3);
+
+ gsl_matrix **rotations = malloc(sizeof(*rotations)* pow(member_num, 2)
+ - member_num);
+
+ int i, j, count;
+
+ count = 0;
+ for ( i=0; i<1; i++ ) {
+ for ( j=0; j<member_num; j++ ) {
+ if (i == j) continue;
+ struct SpotVec i_vec = obs_vecs[obs_members[i]];
+ struct SpotVec j_vec = obs_vecs[obs_members[j]];
+
+ struct rvec i_obsvec = i_vec.obsvec;
+ struct rvec j_obsvec = j_vec.obsvec;
+ struct rvec i_cellvec = i_vec.matches[match_members[i]].vec;
+ struct rvec j_cellvec = j_vec.matches[match_members[j]].vec;
+
+ rotations[count] = generate_rot_mat(i_obsvec, j_obsvec,
+ i_cellvec, j_cellvec,
+ cell);
+
+ count++;
+ }
+ }
+
+ double min_score = FLT_MAX;
+ int min_rot_index = 0;
+
+ for (i=0; i<count; i++) {
+ double current_score = 0;
+ for (j=0; j<count; j++) {
+ double addition;
+ rot_mats_are_similar(rotations[i], rotations[j],
+ sub, mul,
+ &addition, cell);
+
+ current_score += addition;
+ }
+
+ if (current_score < min_score) {
+ min_score = current_score;
+ min_rot_index = i;
+ }
+ }
+
+ gsl_matrix_memcpy(rot, rotations[min_rot_index]);
+
+ for (i=0; i<count; i++) {
+ gsl_matrix_free(rotations[i]);
+ }
+
+ free(rotations);
+ gsl_matrix_free(sub);
+ gsl_matrix_free(mul);
+
+ return 1;
+}
+
+gsl_matrix *rot_mat_from_indices(int her, int his,
+ int her_match, int his_match,
+ struct TakeTwoCell *cell)
+{
+ struct SpotVec *obs_vecs = cell->obs_vecs;
+ struct SpotVec *her_obs = &obs_vecs[her];
+ struct SpotVec *his_obs = &obs_vecs[his];
+ struct rvec i_obsvec = her_obs->obsvec;
+ struct rvec j_obsvec = his_obs->obsvec;
+ struct rvec i_cellvec = her_obs->matches[her_match].vec;
+ struct rvec j_cellvec = his_obs->matches[his_match].vec;
+
+ gsl_matrix *mat = generate_rot_mat(i_obsvec, j_obsvec,
+ i_cellvec, j_cellvec, cell);
+
+ return mat;
+}
+
+static int weed_duplicate_matches(struct Seed **seeds,
+ int *match_count, struct TakeTwoCell *cell)
+{
+ int num_occupied = 0;
+ gsl_matrix **old_mats = calloc(*match_count, sizeof(gsl_matrix *));
+
+ if (old_mats == NULL)
+ {
+ apologise();
+ return 0;
+ }
+
+ signed int i, j;
+
+ int duplicates = 0;
+
+ /* Now we weed out the self-duplicates from the remaining batch */
+
+ for (i = *match_count - 1; i >= 0; i--) {
+ int her_match = (*seeds)[i].idx1;
+ int his_match = (*seeds)[i].idx2;
+
+ gsl_matrix *mat;
+ mat = rot_mat_from_indices((*seeds)[i].obs1, (*seeds)[i].obs2,
+ her_match, his_match, cell);
+
+ int found = 0;
+
+ for (j = 0; j < num_occupied; j++) {
+ if (old_mats[j] && mat &&
+ symm_rot_mats_are_similar(old_mats[j], mat, cell))
+ {
+ // we have found a duplicate, so flag as bad.
+ (*seeds)[i].idx1 = -1;
+ (*seeds)[i].idx2 = -1;
+ found = 1;
+
+ duplicates++;
+
+ gsl_matrix_free(mat);
+ break;
+ }
+ }
+
+ if (!found) {
+ // we have not found a duplicate, add to list.
+ old_mats[num_occupied] = mat;
+ num_occupied++;
+ }
+ }
+
+ for (i = 0; i < num_occupied; i++) {
+ if (old_mats[i]) {
+ gsl_matrix_free(old_mats[i]);
+ }
+ }
+
+ free(old_mats);
+
+ return 1;
+}
+
+static signed int find_next_index(gsl_matrix *rot, int *obs_members,
+ int *match_members, int start, int member_num,
+ int *match_found, struct TakeTwoCell *cell)
+{
+ struct SpotVec *obs_vecs = cell->obs_vecs;
+ int obs_vec_count = cell->obs_vec_count;
+ gsl_matrix *sub = gsl_matrix_calloc(3, 3);
+ gsl_matrix *mul = gsl_matrix_calloc(3, 3);
+
+ int i, j, k;
+
+ for ( i=start; i<obs_vec_count; i++ ) {
+
+ /* If we've considered this vector before, ignore it */
+ if (obs_vecs[i].in_network == 1)
+ {
+ continue;
+ }
+
+ /* first we check for a shared spot - harshest condition */
+ int shared = obs_shares_spot_w_array(obs_vecs, i, obs_members,
+ member_num);
+
+ if ( !shared ) continue;
+
+ int all_ok = 1;
+ int matched = -1;
+
+ /* Check all existing members are happy to let in the newcomer */
+ for ( j=0; j<member_num && all_ok; j++ ) {
+
+ struct SpotVec *me = &obs_vecs[i];
+ struct SpotVec *you = &obs_vecs[obs_members[j]];
+ struct rvec you_cell;
+ you_cell = you->matches[match_members[j]].vec;
+
+ struct rvec me_obs = me->obsvec;
+ struct rvec you_obs = you->obsvec;
+
+ int one_is_okay = 0;
+
+ /* Loop through all possible theoretical vector
+ * matches for the newcomer.. */
+
+ for ( k=0; k<me->match_num; k++ ) {
+
+ gsl_matrix *test_rot;
+
+ struct rvec me_cell = me->matches[k].vec;
+
+ test_rot = generate_rot_mat(me_obs,
+ you_obs, me_cell, you_cell,
+ cell);
+
+ double trace = 0;
+ int ok = rot_mats_are_similar(rot, test_rot,
+ sub, mul, &trace, cell);
+
+ gsl_matrix_free(test_rot);
+
+ if (ok) {
+ one_is_okay = 1;
+
+ /* We are only happy if the vector
+ * matches for only one kind of
+ * theoretical vector. We don't want to
+ * accept mixtures of theoretical vector
+ * matches. */
+ if (matched >= 0 && k == matched) {
+ *match_found = k;
+ } else if (matched < 0) {
+ matched = k;
+ } else {
+ one_is_okay = 0;
+ break;
+ }
+ }
+ }
+
+ if (!one_is_okay) {
+ all_ok = 0;
+ break;
+ }
+ }
+
+
+ if (all_ok) {
+ gsl_matrix_free(sub);
+ gsl_matrix_free(mul);
+ return i;
+ }
+ }
+
+ /* give up. */
+ gsl_matrix_free(sub);
+ gsl_matrix_free(mul);
+ return -1;
+}
+
+/**
+ * Reward target function for refining solution to all vectors in a
+ * given image. Sum of exponentials of the negative distances, which
+ * means that the reward decays as the distance from the nearest
+ * theoretical vector reduces. */
+static double obs_to_sol_score(struct TakeTwoCell *ttCell)
+{
+ double total = 0;
+ int count = 0;
+ int i;
+ gsl_matrix *solution = ttCell->solution;
+ gsl_matrix *full_rot = gsl_matrix_alloc(3, 3);
+ rotate_gsl_by_angles(solution, ttCell->x_ang, ttCell->y_ang,
+ ttCell->z_ang, full_rot);
+
+ for (i = 0; i < ttCell->obs_vec_count; i++)
+ {
+ struct rvec *obs = &ttCell->obs_vecs[i].obsvec;
+ rvec_to_gsl(ttCell->vec1Tmp, *obs);
+
+ /* Rotate all the observed vectors by the modified soln */
+ /* ttCell->vec2Tmp = 1.0 * full_rot * ttCell->vec1Tmp */
+ gsl_blas_dgemv(CblasTrans, 1.0, full_rot, ttCell->vec1Tmp,
+ 0.0, ttCell->vec2Tmp);
+ struct rvec rotated = gsl_to_rvec(ttCell->vec2Tmp);
+
+ int j = ttCell->obs_vecs[i].assignment;
+
+ if (j < 0) continue;
+
+ struct rvec *match = &ttCell->obs_vecs[i].matches[j].vec;
+ struct rvec diff = diff_vec(rotated, *match);
+
+ double length = rvec_length(diff);
+
+ double addition = exp(-(1 / RECIP_TOLERANCE) * length);
+ total += addition;
+ count++;
+ }
+
+ total /= (double)-count;
+
+ gsl_matrix_free(full_rot);
+
+ return total;
+}
+
+/**
+ * Matches every observed vector in the image to its closest theoretical
+ * neighbour after applying the rotation matrix, in preparation for
+ * refining the rotation matrix to match these. */
+static void match_all_obs_to_sol(struct TakeTwoCell *ttCell)
+{
+ int i, j;
+ double total = 0;
+ int count = 0;
+ gsl_matrix *solution = ttCell->solution;
+
+ for (i = 0; i < ttCell->obs_vec_count; i++)
+ {
+ struct rvec *obs = &ttCell->obs_vecs[i].obsvec;
+ rvec_to_gsl(ttCell->vec1Tmp, *obs);
+
+ /* ttCell->vec2Tmp = 1.0 * solution * ttCell->vec1Tmp */
+ gsl_blas_dgemv(CblasTrans, 1.0, solution, ttCell->vec1Tmp,
+ 0.0, ttCell->vec2Tmp);
+ struct rvec rotated = gsl_to_rvec(ttCell->vec2Tmp);
+
+ double smallest = FLT_MAX;
+ int assigned = -1;
+
+ for (j = 0; j < ttCell->obs_vecs[i].match_num; j++)
+ {
+ struct rvec *match = &ttCell->obs_vecs[i].matches[j].vec;
+ struct rvec diff = diff_vec(rotated, *match);
+
+ double length = rvec_length(diff);
+ if (length < smallest)
+ {
+ smallest = length;
+ assigned = j;
+ }
+ }
+
+ ttCell->obs_vecs[i].assignment = assigned;
+
+ if (smallest != FLT_MAX)
+ {
+ double addition = exp(-(1 / RECIP_TOLERANCE) * smallest);
+ total += addition;
+ count++;
+
+ }
+ }
+
+ total /= (double)count;
+}
+
+/**
+ * Refines a matrix against all of the observed vectors against their
+ * closest theoretical neighbour, by perturbing the matrix along the principle
+ * axes until it maximises a reward function consisting of the sum of
+ * the distances of individual observed vectors to their closest
+ * theoretical neighbour. Reward function means that noise and alternative
+ * lattices do not dominate the equation!
+ **/
+static void refine_solution(struct TakeTwoCell *ttCell)
+{
+ match_all_obs_to_sol(ttCell);
+
+ int i, j, k;
+ const int total = 3 * 3 * 3;
+ const int middle = (total - 1) / 2;
+
+ struct rvec steps[total];
+ double start = obs_to_sol_score(ttCell);
+ const int max_tries = 100;
+
+ int count = 0;
+ double size = ANGLE_STEP_SIZE;
+
+ /* First we create our combinations of steps */
+ for (i = -1; i <= 1; i++) {
+ for (j = -1; j <= 1; j++) {
+ for (k = -1; k <= 1; k++) {
+ struct rvec vec = new_rvec(i, j, k);
+ steps[count] = vec;
+ count++;
+ }
+ }
+ }
+
+ struct rvec current = new_rvec(ttCell->x_ang, ttCell->y_ang,
+ ttCell->z_ang);
+
+ double best = start;
+ count = 0;
+
+ while (size > ANGLE_CONVERGE_SIZE && count < max_tries)
+ {
+ struct rvec sized[total];
+
+ int best_num = middle;
+ for (i = 0; i < total; i++)
+ {
+ struct rvec sized_step = steps[i];
+ sized_step.u *= size;
+ sized_step.v *= size;
+ sized_step.w *= size;
+
+ struct rvec new_angles = rvec_add_rvec(current,
+ sized_step);
+
+ sized[i] = new_angles;
+
+ ttCell->x_ang = sized[i].u;
+ ttCell->y_ang = sized[i].v;
+ ttCell->z_ang = sized[i].w;
+
+ double score = obs_to_sol_score(ttCell);
+
+ if (score < best)
+ {
+ best = score;
+ best_num = i;
+ }
+ }
+
+ if (best_num == middle)
+ {
+ size /= 2;
+ }
+
+ current = sized[best_num];
+ count++;
+ }
+
+ ttCell->x_ang = 0;
+ ttCell->y_ang = 0;
+ ttCell->z_ang = 0;
+
+ gsl_matrix *tmp = gsl_matrix_alloc(3, 3);
+ rotate_gsl_by_angles(ttCell->solution, current.u,
+ current.v, current.w, tmp);
+ gsl_matrix_free(ttCell->solution);
+ ttCell->solution = tmp;
+}
+
+
+static unsigned int grow_network(gsl_matrix *rot, int obs_idx1, int obs_idx2,
+ int match_idx1, int match_idx2,
+ struct TakeTwoCell *cell)
+{
+
+ struct SpotVec *obs_vecs = cell->obs_vecs;
+ int obs_vec_count = cell->obs_vec_count;
+ int *obs_members;
+ int *match_members;
+
+ /* Clear the in_network status of all vectors to start */
+ int i;
+ for (i = 0; i < obs_vec_count; i++)
+ {
+ obs_vecs[i].in_network = 0;
+ }
+
+ /* indices of members of the self-consistent network of vectors */
+ obs_members = malloc((cell->member_thresh+3)*sizeof(int));
+ match_members = malloc((cell->member_thresh+3)*sizeof(int));
+ if ( (obs_members == NULL) || (match_members == NULL) ) {
+ apologise();
+ return 0;
+ }
+
+ /* initialise the ones we know already */
+ obs_members[0] = obs_idx1;
+ obs_members[1] = obs_idx2;
+ match_members[0] = match_idx1;
+ match_members[1] = match_idx2;
+ int member_num = 2;
+
+ /* counter for dead ends which must not exceed MAX_DEAD_ENDS
+ * before it is reset in an additional branch */
+ int dead_ends = 0;
+
+ /* we start from 0 */
+ int start = 0;
+
+ while ( 1 ) {
+
+ if (start > obs_vec_count) {
+ free(obs_members);
+ free(match_members);
+ return 0;
+ }
+
+ int match_found = -1;
+
+ signed int next_index = find_next_index(rot, obs_members,
+ match_members,
+ 0, member_num,
+ &match_found, cell);
+
+ if ( member_num < 2 ) {
+ free(obs_members);
+ free(match_members);
+ return 0;
+ }
+
+ if ( next_index < 0 ) {
+ /* If there have been too many dead ends, give up
+ * on indexing altogether.
+ **/
+ if ( dead_ends > MAX_DEAD_ENDS ) {
+ break;
+ }
+
+ /* We have not had too many dead ends. Try removing
+ the last member and continue. */
+ member_num--;
+ dead_ends++;
+
+ continue;
+ }
+
+ /* Elongation of the network was successful */
+ obs_vecs[next_index].in_network = 1;
+ obs_members[member_num] = next_index;
+ match_members[member_num] = match_found;
+
+ member_num++;
+
+ /* If member_num is high enough, we want to return a yes */
+ if ( member_num > cell->member_thresh ) break;
+
+ }
+
+ finish_solution(rot, obs_vecs, obs_members,
+ match_members, member_num, cell);
+
+ free(obs_members);
+ free(match_members);
+
+ return ( member_num );
+}
+
+
+static unsigned int start_seed(int i, int j, int i_match, int j_match,
+ gsl_matrix **rotation, struct TakeTwoCell *cell)
+{
+ struct SpotVec *obs_vecs = cell->obs_vecs;
+
+ gsl_matrix *rot_mat;
+
+ rot_mat = generate_rot_mat(obs_vecs[i].obsvec,
+ obs_vecs[j].obsvec,
+ obs_vecs[i].matches[i_match].vec,
+ obs_vecs[j].matches[j_match].vec,
+ cell);
+
+ /* Try to expand this rotation matrix to a larger network */
+
+ int member_num = grow_network(rot_mat, i, j, i_match, j_match,
+ cell);
+
+ /* return this matrix and if it was immediately successful */
+ *rotation = rot_mat;
+
+ return member_num;
+}
+
+static int sort_seed_by_score(const void *av, const void *bv)
+{
+ struct Seed *a = (struct Seed *)av;
+ struct Seed *b = (struct Seed *)bv;
+ return a->score > b->score;
+}
+
+static void remove_old_solutions(struct TakeTwoCell *cell,
+ struct taketwo_private *tp)
+{
+ int duplicates = 0;
+ struct Seed *seeds = cell->seeds;
+ unsigned int total = cell->seed_count;
+
+ /* First we remove duplicates with previous solutions */
+
+ int i, j;
+ for (i = total - 1; i >= 0; i--) {
+ int her_match = seeds[i].idx1;
+ int his_match = seeds[i].idx2;
+
+ gsl_matrix *mat;
+ mat = rot_mat_from_indices(seeds[i].obs1, seeds[i].obs2,
+ her_match, his_match, cell);
+
+ if (mat == NULL)
+ {
+ continue;
+ }
+
+ for (j = 0; j < tp->numPrevs; j++)
+ {
+ int sim = symm_rot_mats_are_similar(tp->prevSols[j],
+ mat, cell);
+
+ /* Found a duplicate with a previous solution */
+ if (sim)
+ {
+ seeds[i].idx1 = -1;
+ seeds[i].idx2 = -1;
+ duplicates++;
+ break;
+ }
+ }
+
+ gsl_matrix_free(mat);
+ }
+
+// STATUS("Removing %i duplicates due to prev solutions.\n", duplicates);
+}
+
+static int find_seeds(struct TakeTwoCell *cell, struct taketwo_private *tp)
+{
+ struct SpotVec *obs_vecs = cell->obs_vecs;
+ int obs_vec_count = cell->obs_vec_count;
+
+ /* loop round pairs of vectors to try and find a suitable
+ * seed to start building a self-consistent network of vectors
+ */
+ int i, j;
+
+ for ( i=1; i<obs_vec_count; i++ ) {
+
+ for ( j=0; j<i; j++ ) {
+
+ /** Only check distances which are accumulatively less
+ * than the limit if we can easily generate seeds */
+ if (obs_vecs[j].distance + obs_vecs[i].distance >
+ MAX_RECIP_DISTANCE && cell->seed_count > 100) {
+ continue;
+ }
+
+ /** Check to see if there is a shared spot - opportunity
+ * for optimisation by generating a look-up table
+ * by spot instead of by vector.
+ */
+ int shared = obs_vecs_share_spot(&obs_vecs[i],
+ &obs_vecs[j]);
+ if ( !shared ) continue;
+
+ /* cell vector index matches stored in i, j and total
+ * number stored in int matches.
+ */
+ int seed_num = 0;
+ struct Seed *seeds = NULL;
+
+ /* Check to see if any angles match from the cell
+ * vectors */
+ obs_vecs_match_angles(i, j, &seeds, &seed_num, cell);
+
+ if (seed_num == 0)
+ {
+ /* Nothing to clean up here */
+ continue;
+ }
+
+ /* Weed out the duplicate seeds (from symmetric
+ * reflection pairs) */
+ weed_duplicate_matches(&seeds, &seed_num, cell);
+
+ /* Add all the new seeds to the full list */
+
+ size_t new_size = cell->seed_count + seed_num;
+ new_size *= sizeof(struct Seed);
+
+ struct Seed *tmp = realloc(cell->seeds, new_size);
+
+ if (tmp == NULL) {
+ apologise();
+ }
+
+ cell->seeds = tmp;
+
+ for (int i = 0; i < seed_num; i++)
+ {
+ if (seeds[i].idx1 < 0 || seeds[i].idx2 < 0)
+ {
+ continue;
+ }
+
+ cell->seeds[cell->seed_count] = seeds[i];
+ cell->seed_count++;
+ }
+
+ free(seeds);
+ }
+ }
+
+ qsort(cell->seeds, cell->seed_count, sizeof(struct Seed),
+ sort_seed_by_score);
+
+
+ return 1;
+}
+
+static unsigned int start_seeds(gsl_matrix **rotation, struct TakeTwoCell *cell)
+{
+ struct Seed *seeds = cell->seeds;
+ int seed_num = cell->seed_count;
+ int member_num = 0;
+ int max_members = 0;
+ gsl_matrix *rot = NULL;
+
+ /* We have seeds! Pass each of them through the seed-starter */
+ /* If a seed has the highest achieved membership, make note...*/
+ int k;
+ for ( k=0; k<seed_num; k++ ) {
+ int seed_idx1 = seeds[k].idx1;
+ int seed_idx2 = seeds[k].idx2;
+
+ if (seed_idx1 < 0 || seed_idx2 < 0) {
+ continue;
+ }
+
+ int seed_obs1 = seeds[k].obs1;
+ int seed_obs2 = seeds[k].obs2;
+
+ member_num = start_seed(seed_obs1, seed_obs2, seed_idx1,
+ seed_idx2, &rot, cell);
+
+ if (member_num > max_members)
+ {
+ if ( *rotation != NULL ) {
+ /* Free previous best */
+ gsl_matrix_free(*rotation);
+ }
+ *rotation = rot;
+ max_members = member_num;
+ } else {
+ gsl_matrix_free(rot);
+ }
+
+ if (member_num >= NETWORK_MEMBER_THRESHOLD) {
+ free(seeds);
+ return max_members;
+ }
+ }
+
+ free(seeds);
+
+ return max_members;
+}
+
+
+static void set_gsl_matrix(gsl_matrix *mat,
+ double asx, double asy, double asz,
+ double bsx, double bsy, double bsz,
+ double csx, double csy, double csz)
+{
+ gsl_matrix_set(mat, 0, 0, asx);
+ gsl_matrix_set(mat, 0, 1, asy);
+ gsl_matrix_set(mat, 0, 2, asz);
+ gsl_matrix_set(mat, 1, 0, bsx);
+ gsl_matrix_set(mat, 1, 1, bsy);
+ gsl_matrix_set(mat, 1, 2, bsz);
+ gsl_matrix_set(mat, 2, 0, csx);
+ gsl_matrix_set(mat, 2, 1, csy);
+ gsl_matrix_set(mat, 2, 2, csz);
+}
+
+static int generate_rotation_sym_ops(struct TakeTwoCell *ttCell)
+{
+ SymOpList *rawList = sym_ops_for_cell(ttCell->cell);
+
+ /* Now we must convert these into rotation matrices rather than hkl
+ * transformations (affects triclinic, monoclinic, rhombohedral and
+ * hexagonal space groups only) */
+
+ double asx, asy, asz;
+ double bsx, bsy, bsz;
+ double csx, csy, csz;
+
+ gsl_matrix *recip = gsl_matrix_alloc(3, 3);
+ gsl_matrix *cart = gsl_matrix_alloc(3, 3);
+
+ cell_get_reciprocal(ttCell->cell, &asx, &asy, &asz,
+ &bsx, &bsy, &bsz,
+ &csx, &csy, &csz);
+
+ set_gsl_matrix(recip, asx, asy, asz,
+ asx, bsy, bsz,
+ csx, csy, csz);
+
+ cell_get_cartesian(ttCell->cell, &asx, &asy, &asz,
+ &bsx, &bsy, &bsz,
+ &csx, &csy, &csz);
+
+ set_gsl_matrix(cart, asx, bsx, csx,
+ asy, bsy, csy,
+ asz, bsz, csz);
+
+
+ int i, j, k;
+ int numOps = num_equivs(rawList, NULL);
+
+ ttCell->rotSymOps = malloc(numOps * sizeof(gsl_matrix *));
+ ttCell->numOps = numOps;
+
+ if (ttCell->rotSymOps == NULL) {
+ apologise();
+ return 0;
+ }
+
+ for (i = 0; i < numOps; i++)
+ {
+ gsl_matrix *symOp = gsl_matrix_alloc(3, 3);
+ IntegerMatrix *op = get_symop(rawList, NULL, i);
+
+ for (j = 0; j < 3; j++) {
+ for (k = 0; k < 3; k++) {
+ gsl_matrix_set(symOp, j, k, intmat_get(op, j, k));
+ }
+ }
+
+ gsl_matrix *first = gsl_matrix_calloc(3, 3);
+ gsl_matrix *second = gsl_matrix_calloc(3, 3);
+
+ /* Each equivalence is of the form:
+ cartesian * symOp * reciprocal.
+ First multiplication: symOp * reciprocal */
+
+ gsl_blas_dgemm(CblasNoTrans, CblasNoTrans,
+ 1.0, symOp, recip,
+ 0.0, first);
+
+ /* Second multiplication: cartesian * first */
+
+ gsl_blas_dgemm(CblasNoTrans, CblasNoTrans,
+ 1.0, cart, first,
+ 0.0, second);
+
+ ttCell->rotSymOps[i] = second;
+
+ gsl_matrix_free(symOp);
+ gsl_matrix_free(first);
+ }
+
+ gsl_matrix_free(cart);
+ gsl_matrix_free(recip);
+
+ free_symoplist(rawList);
+
+ return 1;
+}
+
+struct sortme
+{
+ struct TheoryVec v;
+ double dist;
+};
+
+static int sort_theory_distances(const void *av, const void *bv)
+{
+ struct sortme *a = (struct sortme *)av;
+ struct sortme *b = (struct sortme *)bv;
+ return a->dist > b->dist;
+}
+
+static int match_obs_to_cell_vecs(struct TheoryVec *cell_vecs, int cell_vec_count,
+ struct TakeTwoCell *cell)
+{
+ struct SpotVec *obs_vecs = cell->obs_vecs;
+ int obs_vec_count = cell->obs_vec_count;
+ int i, j;
+
+ for ( i=0; i<obs_vec_count; i++ ) {
+
+ int count = 0;
+ struct sortme *for_sort = NULL;
+
+ for ( j=0; j<cell_vec_count; j++ ) {
+ /* get distance for unit cell vector */
+ double cell_length = rvec_length(cell_vecs[j].vec);
+ double obs_length = obs_vecs[i].distance;
+
+ /* check if this matches the observed length */
+ double dist_diff = fabs(cell_length - obs_length);
+
+ if ( dist_diff > cell->len_tol ) continue;
+
+ /* we have a match, add to array! */
+
+ size_t new_size = (count+1)*sizeof(struct sortme);
+ for_sort = realloc(for_sort, new_size);
+
+ if ( for_sort == NULL ) return 0;
+
+ for_sort[count].v = cell_vecs[j];
+ for_sort[count].dist = dist_diff;
+ count++;
+
+ }
+
+ /* Pointers to relevant things */
+
+ struct TheoryVec **match_array;
+ int *match_count;
+
+ match_array = &(obs_vecs[i].matches);
+ match_count = &(obs_vecs[i].match_num);
+
+ /* Sort in order to get most agreeable matches first */
+ qsort(for_sort, count, sizeof(struct sortme), sort_theory_distances);
+ *match_array = malloc(count*sizeof(struct TheoryVec));
+ *match_count = count;
+ for ( j=0; j<count; j++ ) {
+ (*match_array)[j] = for_sort[j].v;
+
+ }
+ free(for_sort);
+ }
+
+ return 1;
+}
+
+static int compare_spot_vecs(const void *av, const void *bv)
+{
+ struct SpotVec *a = (struct SpotVec *)av;
+ struct SpotVec *b = (struct SpotVec *)bv;
+ return a->distance > b->distance;
+}
+
+static int gen_observed_vecs(struct rvec *rlps, int rlp_count,
+ struct TakeTwoCell *cell)
+{
+ int i, j;
+ int count = 0;
+
+ /* maximum distance squared for comparisons */
+ double max_sq_length = pow(MAX_RECIP_DISTANCE, 2);
+
+ for ( i=0; i<rlp_count-1 && count < MAX_OBS_VECTORS; i++ ) {
+ for ( j=i+1; j<rlp_count; j++ ) {
+
+ /* calculate difference vector between rlps */
+ struct rvec diff = diff_vec(rlps[i], rlps[j]);
+
+ /* are these two far from each other? */
+ double sqlength = sq_length(diff);
+
+ if ( sqlength > max_sq_length ) continue;
+
+ count++;
+
+ struct SpotVec *temp_obs_vecs;
+ temp_obs_vecs = realloc(cell->obs_vecs,
+ count*sizeof(struct SpotVec));
+
+ if ( temp_obs_vecs == NULL ) {
+ return 0;
+ } else {
+ cell->obs_vecs = temp_obs_vecs;
+
+ /* initialise all SpotVec struct members */
+
+ struct SpotVec spot_vec;
+ spot_vec.obsvec = diff;
+ spot_vec.distance = sqrt(sqlength);
+ spot_vec.matches = NULL;
+ spot_vec.assignment = -1;
+ spot_vec.match_num = 0;
+ spot_vec.her_rlp = &rlps[i];
+ spot_vec.his_rlp = &rlps[j];
+
+ cell->obs_vecs[count - 1] = spot_vec;
+ }
+ }
+ }
+
+ /* Sort such that the shortest distances are searched first. */
+ qsort(cell->obs_vecs, count, sizeof(struct SpotVec), compare_spot_vecs);
+
+ cell->obs_vec_count = count;
+
+ return 1;
+}
+
+
+static int gen_theoretical_vecs(UnitCell *cell, struct TheoryVec **cell_vecs,
+ unsigned int *vec_count)
+{
+ double a, b, c, alpha, beta, gamma;
+ int h_max, k_max, l_max;
+ double asx, asy, asz;
+ double bsx, bsy, bsz;
+ double csx, csy, csz;
+
+ cell_get_reciprocal(cell, &asx, &asy, &asz,
+ &bsx, &bsy, &bsz,
+ &csx, &csy, &csz);
+
+ SymOpList *rawList = sym_ops_for_cell(cell);
+
+ cell_get_parameters(cell, &a, &b, &c, &alpha, &beta, &gamma);
+
+ /* find maximum Miller (h, k, l) indices for a given resolution */
+ h_max = MAX_RECIP_DISTANCE * a;
+ k_max = MAX_RECIP_DISTANCE * b;
+ l_max = MAX_RECIP_DISTANCE * c;
+
+ int h, k, l;
+ int _h, _k, _l;
+ int count = 0;
+
+ for ( h=-h_max; h<=+h_max; h++ ) {
+ for ( k=-k_max; k<=+k_max; k++ ) {
+ for ( l=-l_max; l<=+l_max; l++ ) {
+
+ struct rvec cell_vec;
+
+ /* Exclude systematic absences from centering concerns */
+ if ( forbidden_reflection(cell, h, k, l) ) continue;
+
+ int asymmetric = 0;
+ get_asymm(rawList, h, k, l, &_h, &_k, &_l);
+
+ if (h == _h && k == _k && l == _l) {
+ asymmetric = 1;
+ }
+
+ cell_vec.u = h*asx + k*bsx + l*csx;
+ cell_vec.v = h*asy + k*bsy + l*csy;
+ cell_vec.w = h*asz + k*bsz + l*csz;
+
+ struct TheoryVec theory;
+ theory.vec = cell_vec;
+ theory.asym = asymmetric;
+
+ /* add this to our array - which may require expanding */
+ count++;
+
+ struct TheoryVec *temp_cell_vecs;
+ temp_cell_vecs = realloc(*cell_vecs,
+ count*sizeof(struct TheoryVec));
+
+ if ( temp_cell_vecs == NULL ) {
+ return 0;
+ } else {
+ *cell_vecs = temp_cell_vecs;
+ (*cell_vecs)[count - 1] = theory;
+ }
+ }
+ }
+ }
+
+ *vec_count = count;
+
+ free_symoplist(rawList);
+
+ return 1;
+}
+
+
+/* ------------------------------------------------------------------------
+ * cleanup functions - called from run_taketwo().
+ * ------------------------------------------------------------------------*/
+
+static void cleanup_taketwo_obs_vecs(struct SpotVec *obs_vecs,
+ int obs_vec_count)
+{
+ int i;
+ for ( i=0; i<obs_vec_count; i++ ) {
+ free(obs_vecs[i].matches);
+ }
+
+ free(obs_vecs);
+}
+
+static void cleanup_taketwo_cell(struct TakeTwoCell *ttCell)
+{
+ /* n.b. solutions in ttCell are taken care of in the
+ * partial taketwo cleanup. */
+ int i;
+ for ( i=0; i<ttCell->numOps; i++ ) {
+ gsl_matrix_free(ttCell->rotSymOps[i]);
+ }
+ free(ttCell->rotSymOps);
+
+ cleanup_taketwo_obs_vecs(ttCell->obs_vecs,
+ ttCell->obs_vec_count);
+
+ gsl_vector_free(ttCell->vec1Tmp);
+ gsl_vector_free(ttCell->vec2Tmp);
+ gsl_matrix_free(ttCell->twiz1Tmp);
+ gsl_matrix_free(ttCell->twiz2Tmp);
+}
+
+
+/* ------------------------------------------------------------------------
+ * external functions - top level functions (Level 1)
+ * ------------------------------------------------------------------------*/
+
+/**
+ * @cell: target unit cell
+ * @rlps: spot positions on detector back-projected into recripocal
+ * space depending on detector geometry etc.
+ * @rlp_count: number of rlps in rlps array.
+ * @rot: pointer to be given an assignment (hopefully) if indexing is
+ * successful.
+ **/
+static UnitCell *run_taketwo(UnitCell *cell, const struct taketwo_options *opts,
+ struct rvec *rlps, int rlp_count,
+ struct taketwo_private *tp)
+{
+ UnitCell *result;
+ int success = 0;
+ gsl_matrix *solution = NULL;
+
+ /* Initialise TakeTwoCell */
+ struct TakeTwoCell ttCell;
+ ttCell.cell = cell;
+ ttCell.seeds = NULL;
+ ttCell.seed_count = 0;
+ ttCell.rotSymOps = NULL;
+ ttCell.obs_vecs = NULL;
+ ttCell.twiz1Tmp = gsl_matrix_calloc(3, 3);
+ ttCell.twiz2Tmp = gsl_matrix_calloc(3, 3);
+ ttCell.vec1Tmp = gsl_vector_calloc(3);
+ ttCell.vec2Tmp = gsl_vector_calloc(3);
+ ttCell.numOps = 0;
+ ttCell.obs_vec_count = 0;
+ ttCell.solution = NULL;
+ ttCell.x_ang = 0;
+ ttCell.y_ang = 0;
+ ttCell.z_ang = 0;
+
+ success = generate_rotation_sym_ops(&ttCell);
+ if ( !success ) {
+ cleanup_taketwo_cell(&ttCell);
+ return NULL;
+ }
+
+ success = gen_observed_vecs(rlps, rlp_count, &ttCell);
+ if ( !success ) {
+ cleanup_taketwo_cell(&ttCell);
+ return NULL;
+ }
+
+ if ( opts->member_thresh < 0 ) {
+ ttCell.member_thresh = NETWORK_MEMBER_THRESHOLD;
+ } else {
+ ttCell.member_thresh = opts->member_thresh;
+ }
+
+ if ( opts->len_tol < 0.0 ) {
+ ttCell.len_tol = RECIP_TOLERANCE;
+ } else {
+ ttCell.len_tol = opts->len_tol; /* Already in m^-1 */
+ }
+
+ if ( opts->angle_tol < 0.0 ) {
+ ttCell.angle_tol = ANGLE_TOLERANCE;
+ } else {
+ ttCell.angle_tol = opts->angle_tol; /* Already in radians */
+ }
+
+ if ( opts->trace_tol < 0.0 ) {
+ ttCell.trace_tol = sqrt(4.0*(1.0-cos(TRACE_TOLERANCE)));
+ } else {
+ ttCell.trace_tol = sqrt(4.0*(1.0-cos(opts->trace_tol)));
+ }
+
+ success = match_obs_to_cell_vecs(tp->theory_vecs, tp->vec_count,
+ &ttCell);
+
+ if ( !success ) {
+ cleanup_taketwo_cell(&ttCell);
+ return NULL;
+ }
+
+ /* Find all the seeds, then take each one and extend them, returning
+ * a solution if it exceeds the NETWORK_MEMBER_THRESHOLD. */
+ find_seeds(&ttCell, tp);
+ remove_old_solutions(&ttCell, tp);
+ unsigned int members = start_seeds(&solution, &ttCell);
+
+ if ( solution == NULL ) {
+ cleanup_taketwo_cell(&ttCell);
+ return NULL;
+ }
+
+ /* If we have a solution, refine against vectors in the entire image */
+ ttCell.solution = solution;
+ refine_solution(&ttCell);
+ solution = ttCell.solution;
+ double score = obs_to_sol_score(&ttCell);
+
+ /* Add the current solution to the previous solutions list */
+ int new_size = (tp->numPrevs + 1) * sizeof(gsl_matrix *);
+ gsl_matrix **tmp = realloc(tp->prevSols, new_size);
+ double *tmpScores = realloc(tp->prevScores,
+ (tp->numPrevs + 1) * sizeof(double));
+ unsigned int *tmpSuccesses;
+ tmpSuccesses = realloc(tp->membership,
+ (tp->numPrevs + 1) * sizeof(unsigned int));
+
+ if (!tmp) {
+ apologise();
+ }
+
+ tp->prevSols = tmp;
+ tp->prevScores = tmpScores;
+ tp->membership = tmpSuccesses;
+
+ tp->prevSols[tp->numPrevs] = solution;
+ tp->prevScores[tp->numPrevs] = score;
+ tp->membership[tp->numPrevs] = members;
+ tp->numPrevs++;
+
+ /* Prepare the solution for CrystFEL friendliness */
+ result = cell_post_smiley_face(cell, solution);
+ cleanup_taketwo_cell(&ttCell);
+
+ return result;
+}
+
+/* Cleans up the per-image information for taketwo */
+
+static void partial_taketwo_cleanup(struct taketwo_private *tp)
+{
+ if (tp->prevSols != NULL)
+ {
+ for (int i = 0; i < tp->numPrevs; i++)
+ {
+ gsl_matrix_free(tp->prevSols[i]);
+ }
+
+ free(tp->prevSols);
+ }
+
+ free(tp->prevScores);
+ free(tp->membership);
+ tp->prevScores = NULL;
+ tp->membership = NULL;
+ tp->xtal_num = 0;
+ tp->numPrevs = 0;
+ tp->prevSols = NULL;
+
+}
+
+/* CrystFEL interface hooks */
+
+int taketwo_index(struct image *image, void *priv)
+{
+ Crystal *cr;
+ UnitCell *cell;
+ struct rvec *rlps;
+ int n_rlps = 0;
+ int i;
+ struct taketwo_private *tp = (struct taketwo_private *)priv;
+
+ /* Check serial number against previous for solution tracking */
+ int this_serial = image->serial;
+
+ if (tp->serial_num == this_serial)
+ {
+ tp->xtal_num = image->n_crystals;
+ }
+ else
+ {
+ /*
+ for (i = 0; i < tp->numPrevs; i++)
+ {
+ STATUS("score, %i, %.5f, %i\n",
+ this_serial, tp->prevScores[i],
+ tp->membership[i]);
+ }
+ */
+
+ partial_taketwo_cleanup(tp);
+ tp->serial_num = this_serial;
+ tp->xtal_num = image->n_crystals;
+ }
+
+ /*
+ STATUS("Indexing %i with %i attempts, %i crystals\n", this_serial, tp->attempts,
+ image->n_crystals);
+ */
+
+ 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++;
+ }
+ rlps[n_rlps].u = 0.0;
+ rlps[n_rlps].v = 0.0;
+ rlps[n_rlps++].w = 0.0;
+
+ cell = run_taketwo(tp->cell, tp->opts, rlps, n_rlps, tp);
+ free(rlps);
+ if ( cell == NULL ) return 0;
+
+ cr = crystal_new();
+ if ( cr == NULL ) {
+ ERROR("Failed to allocate crystal.\n");
+ return 0;
+ }
+
+ crystal_set_cell(cr, cell);
+
+ image_add_crystal(image, cr);
+
+ return 1;
+}
+
+
+void *taketwo_prepare(IndexingMethod *indm, struct taketwo_options *opts,
+ UnitCell *cell)
+{
+ struct taketwo_private *tp;
+
+ /* Flags that TakeTwo knows about */
+ *indm &= INDEXING_METHOD_MASK | INDEXING_USE_LATTICE_TYPE
+ | INDEXING_USE_CELL_PARAMETERS;
+
+ if ( !( (*indm & INDEXING_USE_LATTICE_TYPE)
+ && (*indm & INDEXING_USE_CELL_PARAMETERS)) )
+ {
+ ERROR("TakeTwo indexing requires cell and lattice type "
+ "information.\n");
+ return NULL;
+ }
+
+ if ( cell == NULL ) {
+ ERROR("TakeTwo indexing requires a unit cell.\n");
+ return NULL;
+ }
+
+ STATUS("*******************************************************************\n");
+ STATUS("***** Welcome to TakeTwo *****\n");
+ STATUS("*******************************************************************\n");
+ STATUS(" If you use these indexing results, please keep a roof\n");
+ STATUS(" over the author's head by citing this paper.\n\n");
+
+ STATUS("o o o o o o o o o o o o\n");
+ STATUS(" o o o o o o o o o o o \n");
+ STATUS("o o\n");
+ STATUS(" o The citation is: o \n");
+ STATUS("o Ginn et al., Acta Cryst. (2016). D72, 956-965 o\n");
+ STATUS(" o Thank you! o \n");
+ STATUS("o o\n");
+ STATUS(" o o o o o o o o o o o \n");
+ STATUS("o o o o o o o o o o o o\n");
+
+
+ STATUS("\n");
+
+ tp = malloc(sizeof(struct taketwo_private));
+ if ( tp == NULL ) return NULL;
+
+ tp->cell = cell;
+ tp->indm = *indm;
+ tp->serial_num = -1;
+ tp->xtal_num = 0;
+ tp->prevSols = NULL;
+ tp->numPrevs = 0;
+ tp->prevScores = NULL;
+ tp->membership = NULL;
+ tp->vec_count = 0;
+ tp->theory_vecs = NULL;
+
+ gen_theoretical_vecs(cell, &tp->theory_vecs, &tp->vec_count);
+
+ return tp;
+}
+
+
+void taketwo_cleanup(IndexingPrivate *pp)
+{
+ struct taketwo_private *tp = (struct taketwo_private *)pp;
+
+ partial_taketwo_cleanup(tp);
+ free(tp->theory_vecs);
+
+ free(tp);
+}
+
+
+const char *taketwo_probe(UnitCell *cell)
+{
+ if ( cell_has_parameters(cell) ) return "taketwo";
+ return NULL;
+}
+
+
+static void taketwo_show_help()
+{
+ printf("Parameters for the TakeTwo indexing algorithm:\n"
+" --taketwo-member-threshold\n"
+" Minimum number of members in network\n"
+" --taketwo-len-tolerance\n"
+" Reciprocal space length tolerance (1/A)\n"
+" --taketwo-angle-tolerance\n"
+" Reciprocal space angle tolerance (in degrees)\n"
+" --taketwo-trace-tolerance\n"
+" Rotation matrix equivalence tolerance (in degrees)\n"
+);
+}
+
+
+static error_t taketwo_parse_arg(int key, char *arg,
+ struct argp_state *state)
+{
+ struct taketwo_options **opts_ptr = state->input;
+ float tmp;
+
+ switch ( key ) {
+
+ case ARGP_KEY_INIT :
+ *opts_ptr = malloc(sizeof(struct taketwo_options));
+ if ( *opts_ptr == NULL ) return ENOMEM;
+ (*opts_ptr)->member_thresh = -1.0;
+ (*opts_ptr)->len_tol = -1.0;
+ (*opts_ptr)->angle_tol = -1.0;
+ (*opts_ptr)->trace_tol = -1.0;
+ break;
+
+ case 1 :
+ taketwo_show_help();
+ return EINVAL;
+
+ case 2 :
+ if ( sscanf(arg, "%i", &(*opts_ptr)->member_thresh) != 1 )
+ {
+ ERROR("Invalid value for --taketwo-member-threshold\n");
+ return EINVAL;
+ }
+ break;
+
+ case 3 :
+ if ( sscanf(arg, "%f", &tmp) != 1 )
+ {
+ ERROR("Invalid value for --taketwo-len-tol\n");
+ return EINVAL;
+ }
+ (*opts_ptr)->len_tol = tmp * 1e10; /* Convert to m^-1 */
+ break;
+
+ case 4 :
+ if ( sscanf(arg, "%f", &tmp) != 1 )
+ {
+ ERROR("Invalid value for --taketwo-angle-tol\n");
+ return EINVAL;
+ }
+ (*opts_ptr)->angle_tol = deg2rad(tmp);
+ break;
+
+ case 5 :
+ if ( sscanf(arg, "%f", &tmp) != 1 )
+ {
+ ERROR("Invalid value for --taketwo-trace-tol\n");
+ return EINVAL;
+ }
+ (*opts_ptr)->trace_tol = deg2rad(tmp);
+ break;
+
+ default :
+ return ARGP_ERR_UNKNOWN;
+
+ }
+
+ return 0;
+}
+
+
+static struct argp_option taketwo_options[] = {
+
+ {"help-taketwo", 1, NULL, OPTION_NO_USAGE,
+ "Show options for TakeTwo indexing algorithm", 99},
+
+ {"taketwo-member-threshold", 2, "n", OPTION_HIDDEN, NULL},
+ {"taketwo-len-tolerance", 3, "one_over_A", OPTION_HIDDEN, NULL},
+ {"taketwo-angle-tolerance", 4, "deg", OPTION_HIDDEN, NULL},
+ {"taketwo-trace-tolerance", 5, "deg", OPTION_HIDDEN, NULL},
+ {0}
+};
+
+
+struct argp taketwo_argp = { taketwo_options, taketwo_parse_arg,
+ NULL, NULL, NULL, NULL, NULL };
diff --git a/libcrystfel/src/indexers/taketwo.h b/libcrystfel/src/indexers/taketwo.h
new file mode 100644
index 00000000..50fa221d
--- /dev/null
+++ b/libcrystfel/src/indexers/taketwo.h
@@ -0,0 +1,47 @@
+/*
+ * taketwo.h
+ *
+ * Rewrite of TakeTwo algorithm (Acta D72 (8) 956-965) for CrystFEL
+ *
+ * Copyright © 2016-2017 Helen Ginn
+ * Copyright © 2016-2020 Deutsches Elektronen-Synchrotron DESY,
+ * a research centre of the Helmholtz Association.
+ *
+ * Authors:
+ * 2016 Helen Ginn <helen@strubi.ox.ac.uk>
+ * 2016-2017 Thomas White <taw@physics.org>
+ *
+ * This file is part of CrystFEL.
+ *
+ * CrystFEL is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * CrystFEL is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with CrystFEL. If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+#ifndef TAKETWO_H
+#define TAKETWO_H
+
+#include <argp.h>
+
+#include "cell.h"
+#include "index.h"
+
+/** \file taketwo.h */
+
+extern void *taketwo_prepare(IndexingMethod *indm, struct taketwo_options *opts,
+ UnitCell *cell);
+extern const char *taketwo_probe(UnitCell *cell);
+extern int taketwo_index(struct image *image, void *priv);
+extern void taketwo_cleanup(IndexingPrivate *pp);
+
+#endif /* TAKETWO_H */
diff --git a/libcrystfel/src/indexers/xds.c b/libcrystfel/src/indexers/xds.c
new file mode 100644
index 00000000..02610d5e
--- /dev/null
+++ b/libcrystfel/src/indexers/xds.c
@@ -0,0 +1,541 @@
+/*
+ * xds.c
+ *
+ * Invoke xds for crystal autoindexing
+ *
+ * Copyright © 2013-2020 Deutsches Elektronen-Synchrotron DESY,
+ * a research centre of the Helmholtz Association.
+ * Copyright © 2013 Cornelius Gati
+ *
+ * Authors:
+ * 2010-2018 Thomas White <taw@physics.org>
+ * 2013 Cornelius Gati <cornelius.gati@cfel.de>
+ *
+ * This file is part of CrystFEL.
+ *
+ * CrystFEL is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * CrystFEL is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with CrystFEL. If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+#include <string.h>
+#include <unistd.h>
+#include <sys/wait.h>
+#include <fcntl.h>
+#include <assert.h>
+#include <sys/ioctl.h>
+#include <errno.h>
+
+#ifdef HAVE_FORKPTY_PTY_H
+#include <pty.h>
+#endif
+#ifdef HAVE_FORKPTY_UTIL_H
+#include <util.h>
+#endif
+
+#include "xds.h"
+#include "cell.h"
+#include "image.h"
+#include "utils.h"
+#include "peaks.h"
+#include "cell-utils.h"
+
+/** \file xds.h */
+
+/* Fake pixel size and camera length, both in metres */
+#define FAKE_PIXEL_SIZE (70e-6)
+#define FAKE_CLEN (0.1)
+
+
+/* Global private data, prepared once */
+struct xds_private
+{
+ IndexingMethod indm;
+ UnitCell *cell;
+};
+
+
+/* Essentially the reverse of xds_spacegroup_for_lattice(), below */
+static int convert_spacegroup_number(int spg, LatticeType *lt, char *cen,
+ char *ua)
+{
+ switch ( spg ) {
+
+ case 1: *lt = L_TRICLINIC; *cen = 'P'; *ua = '*'; return 0;
+ case 3: *lt = L_MONOCLINIC; *cen = 'P'; *ua = 'b'; return 0;
+ case 5: *lt = L_MONOCLINIC; *cen = 'C'; *ua = 'b'; return 0;
+ case 16: *lt = L_ORTHORHOMBIC; *cen = 'P'; *ua = '*'; return 0;
+ case 21: *lt = L_ORTHORHOMBIC; *cen = 'C'; *ua = '*'; return 0;
+ case 22: *lt = L_ORTHORHOMBIC; *cen = 'F'; *ua = '*'; return 0;
+ case 23: *lt = L_ORTHORHOMBIC; *cen = 'I'; *ua = '*'; return 0;
+ case 75: *lt = L_TETRAGONAL; *cen = 'P'; *ua = 'c'; return 0;
+ case 79: *lt = L_TETRAGONAL; *cen = 'I'; *ua = 'c'; return 0;
+ case 143: *lt = L_HEXAGONAL; *cen = 'P'; *ua = 'c'; return 0;
+ case 146: *lt = L_HEXAGONAL; *cen = 'H'; *ua = 'c'; return 0;
+ case 195: *lt = L_CUBIC; *cen = 'P'; *ua = '*'; return 0;
+ case 196: *lt = L_CUBIC; *cen = 'F'; *ua = '*'; return 0;
+ case 197: *lt = L_CUBIC; *cen = 'I'; *ua = '*'; return 0;
+ default: return 1;
+
+ }
+}
+
+
+static int read_cell(struct image *image)
+{
+ FILE * fh;
+ float ax, ay, az;
+ float bx, by, bz;
+ float cx, cy, cz;
+ int spg;
+ char *rval, line[1024];
+ UnitCell *cell;
+ LatticeType latticetype;
+ char centering, ua;
+ Crystal *cr;
+
+ fh = fopen("IDXREF.LP", "r");
+ if ( fh == NULL ) return 0; /* Not indexable */
+
+ do {
+ rval = fgets(line, 1023, fh);
+ if ( rval == NULL ) {
+ fclose(fh);
+ return 0;
+ }
+
+ } while ( strcmp(line, " ***** DIFFRACTION PARAMETERS USED AT START OF "
+ "INTEGRATION *****\n") != 0 );
+
+ /* Find and read space group number */
+ do {
+ rval = fgets(line, 1023, fh);
+ if ( rval == NULL ) {
+ fclose(fh);
+ return 0;
+ }
+ } while ( strncmp(line, " SPACE GROUP NUMBER ", 20) != 0 );
+ sscanf(line+20, "%i\n", &spg);
+
+ /* Find and read a */
+ do {
+ rval = fgets(line, 1023, fh);
+ if ( rval == NULL ) {
+ fclose(fh);
+ return 0;
+ }
+ } while ( strncmp(line, " COORDINATES OF UNIT CELL A-AXIS ", 33) != 0 );
+ if ( sscanf(line+33, "%f %f %f\n", &ax, &ay, &az) < 3 ) {
+ fclose(fh);
+ return 0;
+ }
+
+ /* Read b */
+ rval = fgets(line, 1023, fh);
+ if ( rval == NULL ) {
+ fclose(fh);
+ return 0;
+ }
+ if ( sscanf(line+33, "%f %f %f\n", &bx, &by, &bz) < 3 ) {
+ fclose(fh);
+ return 0;
+ }
+
+ /* Read c */
+ rval = fgets(line, 1023, fh);
+ if ( rval == NULL ) {
+ fclose(fh);
+ return 0;
+ }
+ if ( sscanf(line+33, "%f %f %f\n", &cx, &cy, &cz) < 3 ) {
+ fclose(fh);
+ return 0;
+ }
+
+ cell = cell_new();
+ cell_set_cartesian(cell,
+ -ax*1e-10, -ay*1e-10, -az*1e-10,
+ -bx*1e-10, -by*1e-10, -bz*1e-10,
+ cx*1e-10, cy*1e-10, cz*1e-10);
+ if ( convert_spacegroup_number(spg, &latticetype, &centering, &ua) ) {
+ ERROR("Failed to convert XDS space group number (%i)\n", spg);
+ return 0;
+ }
+ cell_set_lattice_type(cell, latticetype);
+ cell_set_centering(cell, centering);
+ cell_set_unique_axis(cell, ua);
+
+ cr = crystal_new();
+ if ( cr == NULL ) {
+ ERROR("Failed to allocate crystal.\n");
+ return 0;
+ }
+ crystal_set_cell(cr, cell);
+ image_add_crystal(image, cr);
+
+ return 1;
+}
+
+
+static void write_spot(struct image *image)
+{
+ FILE *fh;
+ int i;
+ int n;
+
+ fh = fopen("SPOT.XDS", "w");
+ if ( !fh ) {
+ ERROR("Couldn't open temporary file '%s'\n", "SPOT.XDS");
+ return;
+ }
+
+ n = image_feature_count(image->features);
+ for ( i=0; i<n; i++ )
+
+ {
+ struct imagefeature *f;
+ double ttx, tty, x, y;
+
+ f = image_get_feature(image->features, i);
+ if ( f == NULL ) continue;
+ if ( f->intensity <= 0 ) continue;
+
+ ttx = angle_between_2d(0.0, 1.0,
+ f->rx, 1.0/image->lambda + f->rz);
+ tty = angle_between_2d(0.0, 1.0,
+ f->ry, 1.0/image->lambda + f->rz);
+ if ( f->rx < 0.0 ) ttx *= -1.0;
+ if ( f->ry < 0.0 ) tty *= -1.0;
+ x = tan(ttx)*FAKE_CLEN;
+ y = tan(tty)*FAKE_CLEN;
+
+ x = (x / FAKE_PIXEL_SIZE) + 1500;
+ y = (y / FAKE_PIXEL_SIZE) + 1500;
+
+ fprintf(fh, "%10.2f %10.2f %10.2f %10.0f.\n",
+ x, y, 0.5, f->intensity);
+
+ }
+ fclose(fh);
+}
+
+
+/* Turn what we know about the unit cell into something which we can give to
+ * XDS to make it give us only indexing results compatible with the cell. */
+static const char *xds_spacegroup_for_lattice(UnitCell *cell)
+{
+ LatticeType latt;
+ char centering;
+ char *g = NULL;
+
+ latt = cell_get_lattice_type(cell);
+ centering = cell_get_centering(cell);
+
+ switch ( latt )
+ {
+ case L_TRICLINIC :
+ if ( centering == 'P' ) {
+ g = "1";
+ }
+ break;
+
+ case L_MONOCLINIC :
+ if ( centering == 'P' ) {
+ g = "3";
+ } else if ( centering == 'C' ) {
+ g = "5";
+ }
+ break;
+
+ case L_ORTHORHOMBIC :
+ if ( centering == 'P' ) {
+ g = "16";
+ } else if ( centering == 'C' ) {
+ g = "21";
+ } else if ( centering == 'F' ) {
+ g = "22";
+ } else if ( centering == 'I' ) {
+ g = "23";
+ }
+ break;
+
+ case L_TETRAGONAL :
+ if ( centering == 'P' ) {
+ g = "75";
+ } else if ( centering == 'I' ) {
+ g = "79";
+ }
+ break;
+
+ /* Unfortunately, XDS only does "hexagonal H" */
+ case L_RHOMBOHEDRAL :
+ return NULL;
+
+ case L_HEXAGONAL :
+ if ( centering == 'P' ) {
+ g = "143";
+ }
+ if ( centering == 'H' ) {
+ g = "146";
+ }
+ break;
+
+ case L_CUBIC :
+ if ( centering == 'P' ) {
+ g = "195";
+ } else if ( centering == 'F' ) {
+ g = "196";
+ } else if ( centering == 'I' ) {
+ g = "197";
+ }
+ break;
+ }
+
+ return g;
+}
+
+
+static int write_inp(struct image *image, struct xds_private *xp)
+{
+ FILE *fh;
+
+ fh = fopen("XDS.INP", "w");
+ if ( !fh ) {
+ ERROR("Couldn't open XDS.INP\n");
+ return 1;
+ }
+
+ fprintf(fh, "JOB= IDXREF\n");
+ fprintf(fh, "ORGX= 1500\n");
+ fprintf(fh, "ORGY= 1500\n");
+ fprintf(fh, "DETECTOR_DISTANCE= %f\n", FAKE_CLEN*1e3);
+ fprintf(fh, "OSCILLATION_RANGE= 0.300\n");
+ fprintf(fh, "X-RAY_WAVELENGTH= %.6f\n", image->lambda*1e10);
+ fprintf(fh, "NAME_TEMPLATE_OF_DATA_FRAMES=???.img \n");
+ fprintf(fh, "DATA_RANGE=1 1\n");
+ fprintf(fh, "SPOT_RANGE=1 1\n");
+
+ if ( xp->indm & INDEXING_USE_LATTICE_TYPE ) {
+ fprintf(fh, "SPACE_GROUP_NUMBER= %s\n",
+ xds_spacegroup_for_lattice(xp->cell));
+ } else {
+ fprintf(fh, "SPACE_GROUP_NUMBER= 0\n");
+ }
+
+ if ( xp->indm & INDEXING_USE_CELL_PARAMETERS ) {
+
+ double a, b, c, al, be, ga;
+
+ cell_get_parameters(xp->cell, &a, &b, &c, &al, &be, &ga);
+
+ fprintf(fh, "UNIT_CELL_CONSTANTS= "
+ "%.2f %.2f %.2f %.2f %.2f %.2f\n",
+ a*1e10, b*1e10, c*1e10,
+ rad2deg(al), rad2deg(be), rad2deg(ga));
+
+ } else {
+ fprintf(fh, "UNIT_CELL_CONSTANTS= 0 0 0 0 0 0\n");
+ }
+
+ fprintf(fh, "NX= 3000\n");
+ fprintf(fh, "NY= 3000\n");
+ fprintf(fh, "QX= %f\n", FAKE_PIXEL_SIZE*1e3); /* Pixel size in mm */
+ fprintf(fh, "QY= %f\n", FAKE_PIXEL_SIZE*1e3); /* Pixel size in mm */
+ fprintf(fh, "INDEX_ORIGIN=0 0 0\n");
+ fprintf(fh, "DIRECTION_OF_DETECTOR_X-AXIS=1 0 0\n");
+ fprintf(fh, "DIRECTION_OF_DETECTOR_Y-AXIS=0 1 0\n");
+ fprintf(fh, "INCIDENT_BEAM_DIRECTION=0 0 1\n");
+ fprintf(fh, "ROTATION_AXIS=0 1 0\n");
+ fprintf(fh, "DETECTOR= CSPAD\n");
+ fprintf(fh, "MINIMUM_VALID_PIXEL_VALUE= 1\n");
+ fprintf(fh, "OVERLOAD= 200000000\n");
+ fprintf(fh, "REFINE(IDXREF)= CELL ORIENTATION\n");
+ fprintf(fh, "INDEX_ERROR= 0.3\n");
+ fprintf(fh, "INDEX_MAGNITUDE= 15\n");
+ fprintf(fh, "INDEX_QUALITY= 0.2\n");
+ fprintf(fh, "MAXIMUM_ERROR_OF_SPOT_POSITION= 5.0\n");
+
+ fclose(fh);
+
+ return 0;
+}
+
+
+int run_xds(struct image *image, void *priv)
+{
+ int status;
+ int rval;
+ int n;
+ pid_t pid;
+ int pty;
+ struct xds_private *xp = (struct xds_private *)priv;
+
+ if ( write_inp(image, xp) ) {
+ ERROR("Failed to write XDS.INP file for XDS.\n");
+ return 0;
+ }
+
+ n = image_feature_count(image->features);
+ if ( n < 25 ) return 0;
+
+ write_spot(image);
+
+ /* Delete any old indexing result which may exist */
+ remove("IDXREF.LP");
+
+ pid = forkpty(&pty, NULL, NULL, NULL);
+
+ if ( pid == -1 ) {
+ ERROR("Failed to fork for XDS\n");
+ return 0;
+ }
+ if ( pid == 0 ) {
+
+ /* Child process: invoke XDS */
+ struct termios t;
+
+ /* Turn echo off */
+ tcgetattr(STDIN_FILENO, &t);
+ t.c_lflag &= ~(ECHO | ECHOE | ECHOK | ECHONL);
+ tcsetattr(STDIN_FILENO, TCSANOW, &t);
+
+ execlp("xds", "xds", (char *)NULL);
+ ERROR("Failed to invoke XDS.\n");
+ _exit(0);
+
+ }
+ waitpid(pid, &status, 0);
+
+ close(pty);
+ rval = read_cell(image);
+
+ return rval;
+}
+
+
+void *xds_prepare(IndexingMethod *indm, UnitCell *cell)
+{
+ struct xds_private *xp;
+
+ if ( xds_probe(cell) == NULL ) {
+ ERROR("XDS does not appear to run properly.\n");
+ ERROR("Please check your XDS installation.\n");
+ return NULL;
+ }
+
+ /* Either cell,latt and cell provided, or nocell-nolatt and no cell
+ * - complain about anything else. Could figure this out automatically,
+ * but we'd have to decide whether the user just forgot the cell, or
+ * forgot "-nolatt", or whatever. */
+ if ( (*indm & INDEXING_USE_LATTICE_TYPE)
+ && !(*indm & INDEXING_USE_CELL_PARAMETERS) )
+ {
+ ERROR("Invalid XDS options (-latt-nocell): "
+ "try xds-nolatt-nocell.\n");
+ return NULL;
+ }
+
+ if ( (*indm & INDEXING_USE_CELL_PARAMETERS)
+ && !(*indm & INDEXING_USE_LATTICE_TYPE) )
+ {
+ ERROR("Invalid XDS options (-cell-nolatt): "
+ "try xds-nolatt-nocell.\n");
+ return NULL;
+ }
+
+ if ( (cell != NULL) && (xds_spacegroup_for_lattice(cell) == NULL) ) {
+ ERROR("Don't know how to ask XDS for your cell.\n");
+ return NULL;
+ }
+
+ xp = calloc(1, sizeof(*xp));
+ if ( xp == NULL ) return NULL;
+
+ /* Flags that XDS knows about */
+ *indm &= INDEXING_METHOD_MASK | INDEXING_USE_LATTICE_TYPE
+ | INDEXING_USE_CELL_PARAMETERS;
+
+ xp->cell = cell;
+ xp->indm = *indm;
+
+ return xp;
+}
+
+
+void xds_cleanup(void *pp)
+{
+ struct xds_private *xp;
+
+ xp = (struct xds_private *)pp;
+ free(xp);
+}
+
+
+const char *xds_probe(UnitCell *cell)
+{
+ pid_t pid;
+ int pty;
+ int status;
+ FILE *fh;
+ char line[1024];
+ int ok = 0;
+ int l;
+
+ pid = forkpty(&pty, NULL, NULL, NULL);
+ if ( pid == -1 ) {
+ return NULL;
+ }
+ if ( pid == 0 ) {
+
+ /* Child process */
+ struct termios t;
+
+ /* Turn echo off */
+ tcgetattr(STDIN_FILENO, &t);
+ t.c_lflag &= ~(ECHO | ECHOE | ECHOK | ECHONL);
+ tcsetattr(STDIN_FILENO, TCSANOW, &t);
+
+ execlp("xds", "xds", (char *)NULL);
+ _exit(1);
+
+ }
+
+ fh = fdopen(pty, "r");
+
+ for ( l=0; l<10; l++ ) {
+ char *pos;
+ if ( fgets(line, 1024, fh) != NULL ) {
+ pos = strstr(line, "** XDS **");
+ if ( pos != NULL ) {
+ ok = 1;
+ }
+ }
+ }
+
+ fclose(fh);
+ close(pty);
+ waitpid(pid, &status, 0);
+
+ if ( !ok ) return NULL;
+
+ if ( cell_has_parameters(cell) ) return "xds-cell-latt";
+ return "xds-nocell-nolatt";
+}
diff --git a/libcrystfel/src/indexers/xds.h b/libcrystfel/src/indexers/xds.h
new file mode 100644
index 00000000..8c4dc6d0
--- /dev/null
+++ b/libcrystfel/src/indexers/xds.h
@@ -0,0 +1,58 @@
+/*
+ * xds.h
+ *
+ * Invoke xds for crystal autoindexing
+ *
+ * Copyright © 2013-2020 Deutsches Elektronen-Synchrotron DESY,
+ * a research centre of the Helmholtz Association.
+ * Copyright © 2013 Cornelius Gati
+ *
+ * Authors:
+ * 2010-2013,2017 Thomas White <taw@physics.org>
+ * 2013 Cornelius Gati <cornelius.gati@cfel.de>
+ *
+ * This file is part of CrystFEL.
+ *
+ * CrystFEL is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * CrystFEL is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with CrystFEL. If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+#ifndef XDS_H
+#define XDS_H
+
+#include "cell.h"
+#include "index.h"
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+/**
+ * \file xds.h
+ * XDS indexer interface
+ */
+
+extern int run_xds(struct image *image, void *ipriv);
+
+extern void *xds_prepare(IndexingMethod *indm, UnitCell *cell);
+
+extern const char *xds_probe(UnitCell *cell);
+
+extern void xds_cleanup(void *pp);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif /* XDS_H */
diff --git a/libcrystfel/src/indexers/xgandalf.c b/libcrystfel/src/indexers/xgandalf.c
new file mode 100644
index 00000000..2d2dca48
--- /dev/null
+++ b/libcrystfel/src/indexers/xgandalf.c
@@ -0,0 +1,497 @@
+/*
+ * xgandalf.c
+ *
+ * Interface to XGANDALF indexer
+ *
+ * Copyright © 2017-2020 Deutsches Elektronen-Synchrotron DESY,
+ * a research centre of the Helmholtz Association.
+ *
+ * Authors:
+ * 2017-2018 Yaroslav Gevorkov <yaroslav.gevorkov@desy.de>
+ *
+ * This file is part of CrystFEL.
+ *
+ * CrystFEL is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * CrystFEL is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with CrystFEL. If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include "xgandalf.h"
+
+#include <stdlib.h>
+
+#include "utils.h"
+#include "cell-utils.h"
+#include "peaks.h"
+
+#ifdef HAVE_XGANDALF
+#include "xgandalf/adaptions/crystfel/Lattice.h"
+#include "xgandalf/adaptions/crystfel/ExperimentSettings.h"
+#include "xgandalf/adaptions/crystfel/IndexerPlain.h"
+#endif
+
+/** \file xgandalf.h */
+
+struct xgandalf_options {
+ unsigned int sampling_pitch;
+ unsigned int grad_desc_iterations;
+ float tolerance;
+ unsigned int no_deviation_from_provided_cell;
+ float minLatticeVectorLength_A;
+ float maxLatticeVectorLength_A;
+ int maxPeaksForIndexing;
+};
+
+#ifdef HAVE_XGANDALF
+
+struct xgandalf_private_data {
+ IndexerPlain *indexer;
+ reciprocalPeaks_1_per_A_t reciprocalPeaks_1_per_A;
+
+ IndexingMethod indm;
+ UnitCell *cellTemplate;
+ Lattice_t sampleRealLattice_A; //same as cellTemplate
+ IntegerMatrix *centeringTransformation;
+ LatticeTransform_t latticeReductionTransform;
+};
+
+#define FAKE_DETECTOR_DISTANCE (0.1)
+#define FAKE_DETECTOR_RADIUS (0.1)
+#define FAKE_BEAM_ENERGY (1)
+#define FAKE_DIVERGENCE_ANGLE_DEG (0.05)
+#define FAKE_NON_MONOCHROMATICITY (0.005)
+#define FAKE_REFLECTION_RADIUS (0.0001)
+
+#define MAX_ASSEMBLED_LATTICES_COUNT (10)
+
+static void reduceCell(UnitCell* cell, LatticeTransform_t* appliedReductionTransform);
+static void restoreCell(UnitCell *cell, LatticeTransform_t* appliedReductionTransform);
+static void makeRightHanded(UnitCell* cell);
+
+int run_xgandalf(struct image *image, void *ipriv)
+{
+ struct xgandalf_private_data *xgandalf_private_data = (struct xgandalf_private_data*) ipriv;
+ reciprocalPeaks_1_per_A_t *reciprocalPeaks_1_per_A = &(xgandalf_private_data->reciprocalPeaks_1_per_A);
+
+ int peakCountMax = image_feature_count(image->features);
+ reciprocalPeaks_1_per_A->peakCount = 0;
+ for (int i = 0; i < peakCountMax && i < MAX_PEAK_COUNT_FOR_INDEXER; i++) {
+ struct imagefeature *f;
+ f = image_get_feature(image->features, i);
+ if (f == NULL) {
+ continue;
+ }
+
+ reciprocalPeaks_1_per_A->coordinates_x[reciprocalPeaks_1_per_A->peakCount] = f->rx * 1e-10;
+ reciprocalPeaks_1_per_A->coordinates_y[reciprocalPeaks_1_per_A->peakCount] = f->ry * 1e-10;
+ reciprocalPeaks_1_per_A->coordinates_z[reciprocalPeaks_1_per_A->peakCount] = f->rz * 1e-10;
+ reciprocalPeaks_1_per_A->peakCount++;
+ }
+
+ Lattice_t assembledLattices[MAX_ASSEMBLED_LATTICES_COUNT];
+ int assembledLatticesCount;
+ IndexerPlain_index(xgandalf_private_data->indexer,
+ assembledLattices,
+ &assembledLatticesCount,
+ MAX_ASSEMBLED_LATTICES_COUNT,
+ *reciprocalPeaks_1_per_A,
+ NULL);
+
+ if (assembledLatticesCount > 0) { //no multi-lattice at the moment
+ assembledLatticesCount = 1;
+ }
+
+ int goodLatticesCount = assembledLatticesCount;
+ for (int i = 0; i < assembledLatticesCount && i < 1; i++) {
+ reorderLattice(&(xgandalf_private_data->sampleRealLattice_A),
+ &assembledLattices[i]);
+
+ UnitCell *uc;
+ uc = cell_new();
+
+ Lattice_t *l = &assembledLattices[i];
+
+ cell_set_cartesian(uc, l->ax * 1e-10, l->ay * 1e-10, l->az * 1e-10,
+ l->bx * 1e-10, l->by * 1e-10, l->bz * 1e-10,
+ l->cx * 1e-10, l->cy * 1e-10, l->cz * 1e-10);
+ makeRightHanded(uc);
+
+ if(xgandalf_private_data->cellTemplate != NULL){
+ restoreCell(uc, &xgandalf_private_data->latticeReductionTransform);
+
+ UnitCell *new_cell_trans = cell_transform_intmat(uc, xgandalf_private_data->centeringTransformation);
+ cell_free(uc);
+ uc = new_cell_trans;
+
+ cell_set_lattice_type(new_cell_trans, cell_get_lattice_type(xgandalf_private_data->cellTemplate));
+ cell_set_centering(new_cell_trans, cell_get_centering(xgandalf_private_data->cellTemplate));
+ cell_set_unique_axis(new_cell_trans, cell_get_unique_axis(xgandalf_private_data->cellTemplate));
+ }
+
+ if (validate_cell(uc)) {
+ STATUS("Problem with returned cell!\n");
+ }
+
+ Crystal *cr = crystal_new();
+ if (cr == NULL) {
+ ERROR("Failed to allocate crystal.\n");
+ return 0;
+ }
+ crystal_set_cell(cr, uc);
+ image_add_crystal(image, cr);
+
+ }
+
+ return goodLatticesCount;
+}
+
+void *xgandalf_prepare(IndexingMethod *indm, UnitCell *cell,
+ struct xgandalf_options *xgandalf_opts)
+{
+ struct xgandalf_private_data *xgandalf_private_data = malloc(sizeof(struct xgandalf_private_data));
+ allocReciprocalPeaks(&(xgandalf_private_data->reciprocalPeaks_1_per_A));
+ xgandalf_private_data->indm = *indm;
+ xgandalf_private_data->cellTemplate = NULL;
+ xgandalf_private_data->centeringTransformation = NULL;
+
+ float tolerance = xgandalf_opts->tolerance;
+ samplingPitch_t samplingPitch = xgandalf_opts->sampling_pitch;
+ gradientDescentIterationsCount_t gradientDescentIterationsCount = xgandalf_opts->grad_desc_iterations;
+
+ if (*indm & INDEXING_USE_CELL_PARAMETERS) {
+
+ xgandalf_private_data->cellTemplate = cell;
+
+ UnitCell* primitiveCell = uncenter_cell(cell, &xgandalf_private_data->centeringTransformation, NULL);
+
+ reduceCell(primitiveCell, &xgandalf_private_data->latticeReductionTransform);
+
+ double asx, asy, asz, bsx, bsy, bsz, csx, csy, csz;
+ int ret = cell_get_reciprocal(primitiveCell, &asx, &asy, &asz,
+ &bsx, &bsy, &bsz,
+ &csx, &csy, &csz);
+ if (ret != 0) {
+ ERROR("cell_get_reciprocal did not finish properly!");
+ }
+
+ Lattice_t sampleReciprocalLattice_1_per_A = {
+ .ax = asx * 1e-10, .ay = asy * 1e-10, .az = asz * 1e-10,
+ .bx = bsx * 1e-10, .by = bsy * 1e-10, .bz = bsz * 1e-10,
+ .cx = csx * 1e-10, .cy = csy * 1e-10, .cz = csz * 1e-10 };
+
+ double ax, ay, az, bx, by, bz, cx, cy, cz;
+ ret = cell_get_cartesian(primitiveCell, &ax, &ay, &az,
+ &bx, &by, &bz,
+ &cx, &cy, &cz);
+ if (ret != 0) {
+ ERROR("cell_get_cartesian did not finish properly!");
+ }
+ Lattice_t sampleRealLattice_A = {
+ .ax = ax * 1e10, .ay = ay * 1e10, .az = az * 1e10,
+ .bx = bx * 1e10, .by = by * 1e10, .bz = bz * 1e10,
+ .cx = cx * 1e10, .cy = cy * 1e10, .cz = cz * 1e10 };
+ xgandalf_private_data->sampleRealLattice_A = sampleRealLattice_A;
+
+ ExperimentSettings *experimentSettings =
+ ExperimentSettings_new(FAKE_BEAM_ENERGY,
+ FAKE_DETECTOR_DISTANCE,
+ FAKE_DETECTOR_RADIUS,
+ FAKE_DIVERGENCE_ANGLE_DEG,
+ FAKE_NON_MONOCHROMATICITY,
+ sampleReciprocalLattice_1_per_A,
+ tolerance,
+ FAKE_REFLECTION_RADIUS);
+
+ xgandalf_private_data->indexer = IndexerPlain_new(experimentSettings);
+
+ if (xgandalf_opts->no_deviation_from_provided_cell) {
+ IndexerPlain_setRefineWithExactLattice(xgandalf_private_data->indexer, 1);
+ }
+
+ ExperimentSettings_delete(experimentSettings);
+ cell_free(primitiveCell);
+
+ } else {
+
+ Lattice_t sampleRealLattice_A = { 0, 0, 0, 0, 0, 0, 0, 0, 0 };
+ xgandalf_private_data->sampleRealLattice_A = sampleRealLattice_A;
+
+ ExperimentSettings *experimentSettings =
+ ExperimentSettings_new_nolatt(FAKE_BEAM_ENERGY,
+ FAKE_DETECTOR_DISTANCE,
+ FAKE_DETECTOR_RADIUS,
+ FAKE_DIVERGENCE_ANGLE_DEG,
+ FAKE_NON_MONOCHROMATICITY,
+ xgandalf_opts->minLatticeVectorLength_A,
+ xgandalf_opts->maxLatticeVectorLength_A,
+ FAKE_REFLECTION_RADIUS);
+
+ xgandalf_private_data->indexer = IndexerPlain_new(experimentSettings);
+
+ ExperimentSettings_delete(experimentSettings);
+ }
+
+ IndexerPlain_setSamplingPitch(xgandalf_private_data->indexer,
+ samplingPitch);
+ IndexerPlain_setGradientDescentIterationsCount(xgandalf_private_data->indexer,
+ gradientDescentIterationsCount);
+ IndexerPlain_setMaxPeaksToUseForIndexing(xgandalf_private_data->indexer,
+ xgandalf_opts->maxPeaksForIndexing);
+
+ /* Flags that XGANDALF knows about */
+ *indm &= INDEXING_METHOD_MASK | INDEXING_USE_CELL_PARAMETERS;
+
+ return xgandalf_private_data;
+}
+
+
+void xgandalf_cleanup(void *pp)
+{
+ struct xgandalf_private_data *xgandalf_private_data = pp;
+
+ freeReciprocalPeaks(xgandalf_private_data->reciprocalPeaks_1_per_A);
+ IndexerPlain_delete(xgandalf_private_data->indexer);
+ if(xgandalf_private_data->centeringTransformation != NULL){
+ intmat_free(xgandalf_private_data->centeringTransformation);
+ }
+ free(xgandalf_private_data);
+}
+
+static void reduceCell(UnitCell *cell, LatticeTransform_t* appliedReductionTransform)
+{
+ double ax, ay, az, bx, by, bz, cx, cy, cz;
+ cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
+
+ Lattice_t l = { ax, ay, az, bx, by, bz, cx, cy, cz };
+
+ reduceLattice(&l, appliedReductionTransform);
+
+ cell_set_cartesian(cell, l.ax, l.ay, l.az,
+ l.bx, l.by, l.bz,
+ l.cx, l.cy, l.cz);
+
+ makeRightHanded(cell);
+}
+
+static void restoreCell(UnitCell *cell, LatticeTransform_t* appliedReductionTransform)
+{
+
+ double ax, ay, az, bx, by, bz, cx, cy, cz;
+ cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
+
+ Lattice_t l = { ax, ay, az, bx, by, bz, cx, cy, cz };
+
+ restoreLattice(&l, appliedReductionTransform);
+
+ cell_set_cartesian(cell, l.ax, l.ay, l.az,
+ l.bx, l.by, l.bz,
+ l.cx, l.cy, l.cz);
+
+ makeRightHanded(cell);
+}
+
+static void makeRightHanded(UnitCell *cell)
+{
+ double ax, ay, az, bx, by, bz, cx, cy, cz;
+ cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
+
+ if ( !right_handed(cell) ) {
+ cell_set_cartesian(cell, -ax, -ay, -az, -bx, -by, -bz, -cx, -cy, -cz);
+ }
+}
+
+
+const char *xgandalf_probe(UnitCell *cell)
+{
+ return "xgandalf";
+}
+
+#else
+
+int run_xgandalf(struct image *image, void *ipriv)
+{
+ ERROR("This copy of CrystFEL was compiled without XGANDALF support.\n");
+ return 0;
+}
+
+
+void *xgandalf_prepare(IndexingMethod *indm, UnitCell *cell,
+ struct xgandalf_options *xgandalf_opts)
+{
+ ERROR("This copy of CrystFEL was compiled without XGANDALF support.\n");
+ ERROR("To use XGANDALF indexing, recompile with XGANDALF.\n");
+ return NULL;
+}
+
+
+void xgandalf_cleanup(void *pp)
+{
+}
+
+
+const char *xgandalf_probe(UnitCell *cell)
+{
+ return NULL;
+}
+
+#endif // HAVE_XGANDALF
+
+static void xgandalf_show_help()
+{
+ printf("Parameters for the TakeTwo indexing algorithm:\n"
+" --xgandalf-sampling-pitch\n"
+" Sampling pitch: 0 (loosest) to 4 (most dense)\n"
+" or with secondary Miller indices: 5 (loosest) to\n"
+" 7 (most dense). Default: 6\n"
+" --xgandalf-grad-desc-iterations\n"
+" Gradient descent iterations: 0 (few) to 5 (many)\n"
+" Default: 4\n"
+" --xgandalf-fast-execution Shortcut to set\n"
+" --xgandalf-sampling-pitch=2\n"
+" --xgandalf-grad-desc-iterations=3\n"
+" --xgandalf-tolerance Relative tolerance of the lattice vectors.\n"
+" Default is 0.02\n"
+" --xgandalf-no-deviation-from-provided-cell\n"
+" Force the fitted cell to have the same lattice\n"
+" parameters as the provided one\n"
+" --xgandalf-min-lattice-vector-length\n"
+" Minimum possible lattice vector length in A.\n"
+" Default: 30 A\n"
+" --xgandalf-max-lattice-vector-length\n"
+" Maximum possible lattice vector length in A.\n"
+" Default: 250 A\n"
+" --xgandalf-max-peaks\n"
+" Maximum number of peaks used for indexing.\n"
+" All peaks are used for refinement.\n"
+" Default: 250\n"
+);
+}
+
+
+static error_t xgandalf_parse_arg(int key, char *arg,
+ struct argp_state *state)
+{
+ struct xgandalf_options **opts_ptr = state->input;
+
+ switch ( key ) {
+
+ case ARGP_KEY_INIT :
+ *opts_ptr = malloc(sizeof(struct xgandalf_options));
+ if ( *opts_ptr == NULL ) return ENOMEM;
+ (*opts_ptr)->sampling_pitch = 6;
+ (*opts_ptr)->grad_desc_iterations = 4;
+ (*opts_ptr)->tolerance = 0.02;
+ (*opts_ptr)->no_deviation_from_provided_cell = 0;
+ (*opts_ptr)->minLatticeVectorLength_A = 30;
+ (*opts_ptr)->maxLatticeVectorLength_A = 250;
+ (*opts_ptr)->maxPeaksForIndexing = 250;
+ break;
+
+ case 1 :
+ xgandalf_show_help();
+ return EINVAL;
+
+ case 2 :
+ if (sscanf(arg, "%u", &(*opts_ptr)->sampling_pitch) != 1) {
+ ERROR("Invalid value for --xgandalf-sampling-pitch\n");
+ return EINVAL;
+ }
+ break;
+
+ case 3 :
+ if (sscanf(arg, "%u", &(*opts_ptr)->grad_desc_iterations) != 1) {
+ ERROR("Invalid value for --xgandalf-grad-desc-iterations\n");
+ return EINVAL;
+ }
+ break;
+
+ case 4 :
+ if (sscanf(arg, "%f", &(*opts_ptr)->tolerance) != 1) {
+ ERROR("Invalid value for --xgandalf-tolerance\n");
+ return EINVAL;
+ }
+ break;
+
+ case 5 :
+ (*opts_ptr)->no_deviation_from_provided_cell = 1;
+ break;
+
+ case 6 :
+ if (sscanf(arg, "%f", &(*opts_ptr)->minLatticeVectorLength_A) != 1) {
+ ERROR("Invalid value for --xgandalf-min-lattice-vector-length\n");
+ return EINVAL;
+ }
+ break;
+
+ case 7 :
+ if (sscanf(arg, "%f", &(*opts_ptr)->maxLatticeVectorLength_A) != 1) {
+ ERROR("Invalid value for --xgandalf-max-lattice-vector-length\n");
+ return EINVAL;
+ }
+ break;
+
+ case 8 :
+ (*opts_ptr)->sampling_pitch = 2;
+ (*opts_ptr)->grad_desc_iterations = 3;
+ break;
+
+ case 9 :
+ if (sscanf(arg, "%i", &(*opts_ptr)->maxPeaksForIndexing) != 1) {
+ ERROR("Invalid value for --xgandalf-max-peaks\n");
+ return EINVAL;
+ }
+ break;
+
+ }
+
+ return 0;
+}
+
+
+static struct argp_option xgandalf_options[] = {
+
+ {"help-xgandalf", 1, NULL, OPTION_NO_USAGE,
+ "Show options for XGANDALF indexing algorithm", 99},
+
+ {"xgandalf-sampling-pitch", 2, "pitch", OPTION_HIDDEN, NULL},
+ {"xgandalf-sps", 2, "pitch", OPTION_HIDDEN, NULL},
+
+ {"xgandalf-grad-desc-iterations", 3, "n", OPTION_HIDDEN, NULL},
+ {"xgandalf-gdis", 3, "n", OPTION_HIDDEN, NULL},
+
+ {"xgandalf-tolerance", 4, "t", OPTION_HIDDEN, NULL},
+ {"xgandalf-tol", 4, "t", OPTION_HIDDEN, NULL},
+
+ {"xgandalf-no-deviation-from-provided-cell", 5, NULL, OPTION_HIDDEN, NULL},
+ {"xgandalf-ndfpc", 5, NULL, OPTION_HIDDEN, NULL},
+
+ {"xgandalf-min-lattice-vector-length", 6, "len", OPTION_HIDDEN, NULL},
+ {"xgandalf-min-lvl", 6, "len", OPTION_HIDDEN, NULL},
+
+ {"xgandalf-max-lattice-vector-length", 7, "len", OPTION_HIDDEN, NULL},
+ {"xgandalf-max-lvl", 7, "len", OPTION_HIDDEN, NULL},
+
+ {"xgandalf-fast-execution", 8, NULL, OPTION_HIDDEN, NULL},
+
+ {"xgandalf-max-peaks", 9, "n", OPTION_HIDDEN, NULL},
+
+ {0}
+};
+
+
+struct argp xgandalf_argp = { xgandalf_options, xgandalf_parse_arg,
+ NULL, NULL, NULL, NULL, NULL };
diff --git a/libcrystfel/src/indexers/xgandalf.h b/libcrystfel/src/indexers/xgandalf.h
new file mode 100644
index 00000000..288d141a
--- /dev/null
+++ b/libcrystfel/src/indexers/xgandalf.h
@@ -0,0 +1,51 @@
+/*
+ * xgandalf.h
+ *
+ * Interface to XGANDALF indexer
+ *
+ * Copyright © 2017-2020 Deutsches Elektronen-Synchrotron DESY,
+ * a research centre of the Helmholtz Association.
+ *
+ * Authors:
+ * 2017-2018 Yaroslav Gevorkov <yaroslav.gevorkov@desy.de>
+ *
+ * This file is part of CrystFEL.
+ *
+ * CrystFEL is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * CrystFEL is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with CrystFEL. If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+#ifndef LIBCRYSTFEL_SRC_XGANDALF_H
+#define LIBCRYSTFEL_SRC_XGANDALF_H
+
+#include <stddef.h>
+#include <argp.h>
+
+/**
+ * \file xgandalf.h
+ * XGANDALF indexer interface
+ */
+
+#include "index.h"
+
+extern int run_xgandalf(struct image *image, void *ipriv);
+
+extern void *xgandalf_prepare(IndexingMethod *indm, UnitCell *cell,
+ struct xgandalf_options *xgandalf_opts);
+
+extern void xgandalf_cleanup(void *pp);
+extern const char *xgandalf_probe(UnitCell *cell);
+
+
+#endif /* LIBCRYSTFEL_SRC_XGANDALF_H */
generated by cgit v1.2.3 (git 2.25.1) at 2024-06-23 02:53:06 +0000