Move indexers out of API

1 files changed, 0 insertions, 1239 deletions

diff --git a/libcrystfel/src/asdf.c b/libcrystfel/src/asdf.c
deleted file mode 100644
index 1536a109..00000000
--- a/libcrystfel/src/asdf.c
+++ /dev/null
@@ -1,1239 +0,0 @@
-/*
- * 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 */