diff options
Diffstat (limited to 'libcrystfel/src/indexers/asdf.c')
| -rw-r--r-- | libcrystfel/src/indexers/asdf.c | 1239 |
1 files changed, 1239 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 */ |