From 89368cc3992e721c65d3783a76f854cab2643b9b Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 19 Jun 2015 11:50:28 +0200 Subject: Remove ReAx and GrainSpotter Neither of these ever really worked, and better algorithms exist now. --- libcrystfel/src/reax.c | 1212 ------------------------------------------------ 1 file changed, 1212 deletions(-) delete mode 100644 libcrystfel/src/reax.c (limited to 'libcrystfel/src/reax.c') diff --git a/libcrystfel/src/reax.c b/libcrystfel/src/reax.c deleted file mode 100644 index e2dd249a..00000000 --- a/libcrystfel/src/reax.c +++ /dev/null @@ -1,1212 +0,0 @@ -/* - * reax.c - * - * A new auto-indexer - * - * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY, - * a research centre of the Helmholtz Association. - * Copyright © 2014 Takanori Nakane - * - * Authors: - * 2011-2014 Thomas White - * 2014 Takanori Nakane - * - * 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 . - * - */ - -#ifdef HAVE_CONFIG_H -#include -#endif - - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include "image.h" -#include "utils.h" -#include "peaks.h" -#include "cell.h" -#include "cell-utils.h" -#include "index.h" -#include "reax.h" - - -/* Minimum number of standard deviations above the mean a peak must be in the - * 1D FT to qualify as a candidate vector */ -#define MIN_SIGMAS (3.0) - - -/* Maximum number of times the angular tolerance that vectors are permitted to - * be together before they get merged by squash_vectors() */ -#define INC_TOL_MULTIPLIER (3.0) - - -/* Maximum number of candidate vectors to find (we will take the best ones) */ -#define MAX_CANDIDATES (1024) - - -/* Choose the best solution from this many candidate cells */ -#define MAX_REAX_CELL_CANDIDATES (32) - -struct dvec -{ - double x; - double y; - double z; - double th; - double ph; -}; - - -struct reax_candidate -{ - struct dvec v; /* This is the vector for the candidate */ - double fom; -}; - - -struct reax_search_v -{ - unsigned int smin; - unsigned int smax; /* Search for vector in this range */ - - struct reax_candidate *cand; /* Candidate vectors go here */ - int n_cand; /* There are this many candidates */ - int max_warned; -}; - - -struct reax_search -{ - struct reax_search_v *search; /* Search for these vectors */ - int n_search; /* There are this many vectors to find */ - double pmax; /* The maximum feature resolution */ -}; - - -struct reax_private -{ - IndexingMethod indm; - struct dvec *directions; - int n_dir; - double angular_inc; - UnitCell *cell; - - double *fft_in; - fftw_complex *fft_out; - fftw_plan plan; - int nel; - - fftw_complex *r_fft_in; - fftw_complex *r_fft_out; - fftw_plan r_plan; - int ch; - int cw; -}; - - -static void fill_and_transform(struct dvec *dir, ImageFeatureList *flist, - int nel, double pmax, double *fft_in, - fftw_complex *fft_out, fftw_plan plan, - const struct rigid_group *rg, struct detector *det) -{ - int n, i; - - for ( i=0; ifs, f->ss); - assert(p != NULL); - - if ( !panel_is_in_rigid_group(rg, p) ) continue; - - } - - val = f->rx*dir->x + f->ry*dir->y + f->rz*dir->z; - - idx = nel/2 + nel*val/(2.0*pmax); - fft_in[idx]++; - - } - - fftw_execute_dft_r2c(plan, fft_in, fft_out); -} - - -static void add_candidate(struct reax_search_v *s, struct reax_candidate *c) -{ - int idx; - - if ( s->n_cand == MAX_CANDIDATES ) { - if ( !s->max_warned ) { - ERROR("WARNING: Too many candidates.\n"); - s->max_warned = 1; - } - return; - } - - idx = s->n_cand++; - - s->cand[idx].v = c->v; - s->cand[idx].fom = c->fom; -} - - -static double check_dir(struct dvec *dir, ImageFeatureList *flist, - int nel, double pmax, double *fft_in, - fftw_complex *fft_out, fftw_plan plan, - struct reax_search *s, - const struct rigid_group *rg, struct detector *det) -{ - int i; - double tot; - - fill_and_transform(dir, flist, nel, pmax, fft_in, fft_out, - plan, rg, det); - - tot = 0.0; - for ( i=0; in_search; i++ ) { - - double tot = 0.0; - double peak = 0.0; - double peak_mod = 0.0; - double mean; - double sd = 0.0; - int j; - int n = 0; - - for ( j=0; j= s->search[i].smin ) - && ( j <= s->search[i].smax ) ) { - if ( am > peak ) { - peak = am; - peak_mod = (double)j/(2.0*pmax); - } - } - - } - mean = tot/(double)n; - - for ( j=0; j mean+MIN_SIGMAS*sd ) { - - struct reax_candidate c; - - c.v.x = dir->x * peak_mod; - c.v.y = dir->y * peak_mod; - c.v.z = dir->z * peak_mod; - c.fom = peak; - - add_candidate(&s->search[i], &c); - - } - - } - - return tot; -} - - -/* Refine a direct space vector. From Clegg (1984) - * with added iteration because more reflections might get included as the - * refinement proceeds. */ -static double iterate_refine_vector(double *x, double *y, double *z, - ImageFeatureList *flist) -{ - int fi, n, err; - gsl_matrix *C; - gsl_vector *A; - gsl_vector *t; - gsl_matrix *s_vec; - gsl_vector *s_val; - double dtmax; - - A = gsl_vector_calloc(3); - C = gsl_matrix_calloc(3, 3); - - n = image_feature_count(flist); - fesetround(1); - for ( fi=0; firx*(*x) + f->ry*(*y) + f->rz*(*z); /* Sorry ... */ - kn = nearbyint(kno); - if ( fabs(kn - kno) > 0.3 ) continue; - - xv[0] = f->rx; xv[1] = f->ry; xv[2] = f->rz; - - for ( i=0; i<3; i++ ) { - - val = gsl_vector_get(A, i); - gsl_vector_set(A, i, val+xv[i]*kn); - - for ( j=0; j<3; j++ ) { - val = gsl_matrix_get(C, i, j); - gsl_matrix_set(C, i, j, val+xv[i]*xv[j]); - } - - } - - } - - s_val = gsl_vector_calloc(3); - s_vec = gsl_matrix_calloc(3, 3); - err = gsl_linalg_SV_decomp_jacobi(C, s_vec, s_val); - if ( err ) { - ERROR("SVD failed: %s\n", gsl_strerror(err)); - gsl_matrix_free(s_vec); - gsl_vector_free(s_val); - gsl_matrix_free(C); - gsl_vector_free(A); - return 0.0; - } - - t = gsl_vector_calloc(3); - err = gsl_linalg_SV_solve(C, s_vec, s_val, A, t); - if ( err ) { - ERROR("Matrix solution failed: %s\n", gsl_strerror(err)); - gsl_matrix_free(s_vec); - gsl_vector_free(s_val); - gsl_matrix_free(C); - gsl_vector_free(A); - gsl_vector_free(t); - return 0.0; - } - - gsl_matrix_free(s_vec); - gsl_vector_free(s_val); - - dtmax = fabs(*x - gsl_vector_get(t, 0)); - dtmax += fabs(*y - gsl_vector_get(t, 1)); - dtmax += fabs(*z - gsl_vector_get(t, 2)); - - *x = gsl_vector_get(t, 0); - *y = gsl_vector_get(t, 1); - *z = gsl_vector_get(t, 2); - - gsl_matrix_free(C); - gsl_vector_free(A); - gsl_vector_free(t); - - return dtmax; -} - - -static void refine_vector(ImageFeatureList *flist, struct dvec *dir) -{ - int i; - double sm; - - i = 0; - do { - - sm = iterate_refine_vector(&dir->x, &dir->y, &dir->z, flist); - i++; - - } while ( (sm > 0.001e-9) && (i<10) ); -} - - -static void squash_vectors(struct reax_search *s, double tol) -{ - int i; - - for ( i=0; in_search; i++ ) { - - struct reax_search_v *sv; - struct reax_candidate *new; - int j, k; - int n_invalid = 0; - int n_copied; - - sv = &s->search[i]; - - for ( j=0; jn_cand; j++ ) { - for ( k=0; kn_cand; k++ ) { - - struct reax_candidate *v1, *v2; - - if ( j == k ) continue; - - v1 = &sv->cand[j]; - v2 = &sv->cand[k]; - - if ( angle_between(v1->v.x, v1->v.y, v1->v.z, - v2->v.x, v2->v.y, v2->v.z) < tol ) - { - if ( !isnan(v1->fom) && !isnan(v2->fom ) ) { - if ( v1->fom > v2->fom ) { - v2->fom = NAN; - } else { - v1->fom = NAN; - } - n_invalid++; - } - } - - } - } - - new = calloc(sv->n_cand - n_invalid, - sizeof(struct reax_candidate)); - if ( new == NULL ) { - ERROR("Failed to allocate memory for squashed" - " candidate list.\n"); - return; - } - - n_copied = 0; - for ( j=0; jn_cand; j++ ) { - if ( !isnan(sv->cand[j].fom) ) { - - new[n_copied] = sv->cand[j]; - n_copied++; - - } - } - assert(sv->n_cand - n_invalid == n_copied); - - free(sv->cand); - //STATUS("Search vector %i:", i); - //STATUS(" squashed %i candidates down to %i\n", - // sv->n_cand, n_copied); - sv->n_cand = n_copied; - sv->cand = new; - - } -} - - -static UNUSED void show_vectors(struct reax_search *s, const char *pre) -{ - int i; - - /* For each direction being searched for */ - for ( i=0; in_search; i++ ) { - - int j; - - for ( j=0; jsearch[i].n_cand; j++ ) { - STATUS("%s: %i/%i: %+6.2f %+6.2f %+6.2f nm %.2f\n", pre, - i, j, s->search[i].cand[j].v.x*1e9, - s->search[i].cand[j].v.y*1e9, - s->search[i].cand[j].v.z*1e9, - s->search[i].cand[j].fom); - } - - } -} - - -static void find_candidates(struct reax_private *p, - ImageFeatureList *flist, double pmax, - double *fft_in, fftw_complex *fft_out, - struct reax_search *s, - const char *rg, struct detector *det) -{ - int i; - - for ( i=0; in_search; i++ ) { - s->search[i].cand = calloc(MAX_CANDIDATES, - sizeof(struct reax_candidate)); - s->search[i].n_cand = 0; - } - - for ( i=0; in_dir; i++ ) { - check_dir(&p->directions[i], flist, - p->nel, pmax, fft_in, fft_out, p->plan, - s, NULL, NULL); - } - - squash_vectors(s, INC_TOL_MULTIPLIER*p->angular_inc); - - //show_vectors(s, "BEFORE"); - - for ( i=0; in_search; i++ ) { - - struct reax_search_v *sv; - int j; - - sv = &s->search[i]; - - for ( j=0; jn_cand; j++ ) { - refine_vector(flist, &sv->cand[j].v); - } - - } - - //show_vectors(s, "FINAL"); -} - - -/* Set up search parameters to look for all three cell axes */ -static struct reax_search *search_all_axes(UnitCell *cell, double pmax) -{ - double ax, ay, az; - double bx, by, bz; - double cx, cy, cz; - double mod_a, mod_b, mod_c; - double amin, amax; - double bmin, bmax; - double cmin, cmax; - unsigned int smin, smax; - const double ltol = 10.0; /* Direct space axis length tolerance in % */ - struct reax_search *s; - - cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz); - mod_a = modulus(ax, ay, az); - amin = mod_a * (1.0-ltol/100.0); - amax = mod_a * (1.0+ltol/100.0); - - mod_b = modulus(bx, by, bz); - bmin = mod_b * (1.0-ltol/100.0); - bmax = mod_b * (1.0+ltol/100.0); - - mod_c = modulus(cx, cy, cz); - cmin = mod_c * (1.0-ltol/100.0); - cmax = mod_c * (1.0+ltol/100.0); - - s = malloc(3*sizeof(*s)); - s->pmax = pmax; - s->n_search = 3; - s->search = malloc(3*sizeof(struct reax_search_v)); - smin = 2.0*pmax * amin; smax = 2.0*pmax * amax; - s->search[0].smin = smin; s->search[0].smax = smax; - s->search[0].max_warned = 0; - smin = 2.0*pmax * bmin; smax = 2.0*pmax * bmax; - s->search[1].smin = smin; s->search[1].smax = smax; - s->search[1].max_warned = 0; - smin = 2.0*pmax * cmin; smax = 2.0*pmax * cmax; - s->search[2].smin = smin; s->search[2].smax = smax; - s->search[2].max_warned = 0; - - return s; -} - - -static double get_model_phase(double x, double y, double z, ImageFeatureList *f, - int nel, double pmax, double *fft_in, - fftw_complex *fft_out, fftw_plan plan, - int smin, int smax, const struct rigid_group *rg, - struct detector *det) -{ - struct dvec dir; - int s, i; - double max; - double re, im; - - dir.x = x; dir.y = y; dir.z = z; - - fill_and_transform(&dir, f, nel, pmax, fft_in, fft_out, plan, rg, det); - - s = -1; - max = 0.0; - for ( i=smin; i<=smax; i++ ) { - - double re, im, m; - - re = fft_out[i][0]; - im = fft_out[i][1]; - m = sqrt(re*re + im*im); - if ( m > max ) { - max = m; - s = i; - } - - } - - re = fft_out[s][0]; - im = fft_out[s][1]; - - return atan2(im, re); -} - - -static void refine_rigid_group(struct image *image, UnitCell *cell, - const struct rigid_group *rg, double pmax, - double *fft_in, fftw_complex *fft_out, - fftw_plan plan, int smin, int smax, - struct detector *det, struct reax_private *pr) -{ - double ax, ay, az, ma; - double bx, by, bz, mb; - double cx, cy, cz, mc; - double pha, phb, phc; - int i, j; - fftw_complex *r_fft_in; - fftw_complex *r_fft_out; - double m2m; - signed int aix, aiy; - signed int bix, biy; - signed int cix, ciy; - - cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz); - - ma = modulus(ax, ay, az); - mb = modulus(bx, by, bz); - mc = modulus(cx, cy, cz); - - pha = get_model_phase(ax/ma, ay/ma, az/ma, image->features, - pr->nel, pmax, fft_in, fft_out, plan, - smin, smax, rg, det); - phb = get_model_phase(bx/mb, by/mb, bz/mb, image->features, - pr->nel, pmax, fft_in, fft_out, plan, - smin, smax, rg, det); - phc = get_model_phase(cx/mc, cy/mc, cz/mc, image->features, - pr->nel, pmax, fft_in, fft_out, plan, - smin, smax, rg, det); - - r_fft_in = fftw_malloc(pr->cw*pr->ch*sizeof(fftw_complex)); - r_fft_out = fftw_malloc(pr->cw*pr->ch*sizeof(fftw_complex)); - for ( i=0; icw; i++ ) { - for ( j=0; jch; j++ ) { - r_fft_in[i+pr->cw*j][0] = 0.0; - r_fft_in[i+pr->cw*j][1] = 0.0; - } - } - - ma = modulus(ax, ay, 0.0); - mb = modulus(bx, by, 0.0); - mc = modulus(cx, cy, 0.0); - m2m = ma; - if ( mb > m2m ) m2m = mb; - if ( mc > m2m ) m2m = mc; - - aix = (pr->cw/2)*ax/m2m; aiy = (pr->ch/2)*ay/m2m; - bix = (pr->cw/2)*bx/m2m; biy = (pr->ch/2)*by/m2m; - cix = (pr->cw/2)*cx/m2m; ciy = (pr->ch/2)*cy/m2m; - - if ( aix < 0 ) aix += pr->cw/2; - if ( bix < 0 ) bix += pr->cw/2; - if ( cix < 0 ) cix += pr->cw/2; - - if ( aiy < 0 ) aiy += pr->ch/2; - if ( biy < 0 ) biy += pr->ch/2; - if ( ciy < 0 ) ciy += pr->ch/2; - - r_fft_in[aix + pr->cw*aiy][0] = cos(pha); - r_fft_in[aix + pr->cw*aiy][1] = sin(pha); - r_fft_in[pr->cw-aix + pr->cw*(pr->ch-aiy)][0] = cos(pha); - r_fft_in[pr->cw-aix + pr->cw*(pr->ch-aiy)][1] = -sin(pha); - - r_fft_in[bix + pr->cw*biy][0] = cos(phb); - r_fft_in[bix + pr->cw*biy][1] = sin(phb); - r_fft_in[pr->cw-bix + pr->cw*(pr->ch-biy)][0] = cos(phb); - r_fft_in[pr->cw-bix + pr->cw*(pr->ch-biy)][1] = -sin(phb); - - r_fft_in[cix + pr->cw*ciy][0] = cos(phc); - r_fft_in[cix + pr->cw*ciy][1] = sin(phc); - r_fft_in[pr->cw-cix + pr->cw*(pr->ch-ciy)][0] = cos(phc); - r_fft_in[pr->cw-cix + pr->cw*(pr->ch-ciy)][1] = -sin(phc); - - const int tidx = 1; - r_fft_in[tidx][0] = 1.0; - r_fft_in[tidx][1] = 0.0; - -// STATUS("%i %i\n", aix, aiy); -// STATUS("%i %i\n", bix, biy); -// STATUS("%i %i\n", cix, ciy); - - fftw_execute_dft(pr->r_plan, r_fft_in, r_fft_out); - -// max = 0.0; -// FILE *fh = fopen("centering.dat", "w"); -// for ( i=0; icw; i++ ) { -// for ( j=0; jch; j++ ) { -// -// double re, im, am, ph; -// -// re = r_fft_out[i + pr->cw*j][0]; -// im = r_fft_out[i + pr->cw*j][1]; -// am = sqrt(re*re + im*im); -// ph = atan2(im, re); -// -// if ( am > max ) { -// max = am; -// max_i = i; -// max_j = j; -// } -// -// fprintf(fh, "%f ", am); -// -// } -// fprintf(fh, "\n"); -// } -// STATUS("Max at %i, %i\n", max_i, max_j); -// fclose(fh); -// exit(1); - -// STATUS("Offsets for '%s': %.2f, %.2f pixels\n", rg, dx, dy); -} - - -static UNUSED void refine_all_rigid_groups(struct image *image, UnitCell *cell, - double pmax, - double *fft_in, - fftw_complex *fft_out, - fftw_plan plan, int smin, int smax, - struct detector *det, - struct reax_private *p) -{ - int i; - - for ( i=0; idet->n_rigid_groups; i++ ) { - refine_rigid_group(image, cell, image->det->rigid_groups[i], - pmax, fft_in, fft_out, plan, smin, smax, - det, p); - } -} - - -static double max_feature_resolution(ImageFeatureList *flist) -{ - double pmax; - int i, n; - - pmax = 0.0; - n = image_feature_count(flist); - for ( i=0; irx, f->ry, f->rz); - if ( val > pmax ) pmax = val; - - } - - return pmax; -} - - -static int right_handed_vec(struct rvec a, struct rvec b, struct rvec c) -{ - struct rvec aCb; - double aCb_dot_c; - - /* "a" cross "b" */ - aCb.u = a.v*b.w - a.w*b.v; - aCb.v = - (a.u*b.w - a.w*b.u); - aCb.w = a.u*b.v - a.v*b.u; - - /* "a cross b" dot "c" */ - aCb_dot_c = aCb.u*c.u + aCb.v*c.v + aCb.w*c.w; - - if ( aCb_dot_c > 0.0 ) return 1; - return 0; -} - - -struct cell_candidate -{ - UnitCell *cell; - double fom; -}; - - -struct cell_candidate_list -{ - struct cell_candidate *cand; - int n_cand; -}; - - -static int check_twinning(UnitCell *c1, UnitCell *c2, int verbose) -{ - int i; - int n_dup; - const int n_trials = 40; - - double asx, asy, asz; - double bsx, bsy, bsz; - double csx, csy, csz; - double ax, ay, az; - double bx, by, bz; - double cx, cy, cz; - - gsl_rng *rng; - - /* This is a rubbish RNG, but it serves for this purpose: nothing more - * than "I couldn't be bothered to think of n_trials sets of random - * indices". */ - rng = gsl_rng_alloc(gsl_rng_taus2); - - cell_get_reciprocal(c1, &asx, &asy, &asz, - &bsx, &bsy, &bsz, - &csx, &csy, &csz); - cell_get_cartesian(c2, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz); - - n_dup = 0; - for ( i=0; i %5.2f %5.2f %5.2f -> " - "%3i %3i %3i -> %5.2f\n", h, k, l, - h2, k2, l2, h2i, k2i, l2i, dev); - } - - if ( dev < 0.1 ) { - n_dup++; - } - - } - - if ( verbose ) { - STATUS("%i duplicates.\n", n_dup); - } - - gsl_rng_free(rng); - - if ( n_dup > 10 ) return 1; - return 0; -} - - -/* Return true if "cnew" accounts for more than 25% of the peaks predicted by - * any of the "ncells" cells in "cells". */ -static int twinned(UnitCell *cnew, struct cell_candidate_list *cl) -{ - int i; - - for ( i=0; in_cand; i++ ) { - if ( check_twinning(cnew, cl->cand[i].cell, 0) ) return 1; - } - - return 0; -} - - -static int check_vector_combination(struct dvec *vi, struct dvec *vj, - struct dvec *vk, UnitCell *cell) -{ - double ang; - double a, b, c, al, be, ga; - const double angtol = deg2rad(5.0); - - cell_get_parameters(cell, &a, &b, &c, &al, &be, &ga); - - ang = angle_between(vi->x, vi->y, vi->z, vj->x, vj->y, vj->z); - if ( fabs(ang-ga) > angtol ) return 0; - - ang = angle_between(vi->x, vi->y, vi->z, vk->x, vk->y, vk->z); - if ( fabs(ang-be) > angtol ) return 0; - - ang = angle_between(vj->x, vj->y, vj->z, vk->x, vk->y, vk->z); - if ( fabs(ang-al) > angtol ) return 0; - - return 1; -} - - -static void add_cell_candidate(struct cell_candidate_list *cl, UnitCell *cnew, - double fom) -{ - struct cell_candidate cshift; - int i, cpos; - - cpos = cl->n_cand; - for ( i=0; in_cand; i++ ) { - if ( fom > cl->cand[i].fom ) { - cpos = i; - break; - } - } - - cshift.cell = cnew; - cshift.fom = fom; - - for ( i=cpos; in_cand; i++ ) { - - struct cell_candidate cshift2; - cshift2 = cl->cand[i]; - cl->cand[i] = cshift; - cshift = cshift2; - - } - - if ( cl->n_cand >= MAX_REAX_CELL_CANDIDATES ) { - /* "cshift" just fell off the end of the list */ - } else { - cl->cand[cl->n_cand++] = cshift; - } -} - - -static int check_cell(struct reax_private *dp, struct image *image, - UnitCell *cell) -{ - UnitCell *out; - Crystal *cr; - - out = cell_new_from_cell(cell); - - cr = crystal_new(); - if ( cr == NULL ) { - ERROR("Failed to allocate crystal.\n"); - return 0; - } - - crystal_set_cell(cr, out); - - if ( dp->indm & INDEXING_CHECK_PEAKS ) { - if ( !peak_sanity_check(image, &cr, 1) ) { - crystal_free(cr); /* Frees the cell as well */ - return 0; - } - } - - image_add_crystal(image, cr); - - return 1; -} - - -static int assemble_cells_from_candidates(struct image *image, - struct reax_search *s, - UnitCell *cell, - struct reax_private *p) -{ - int i, j, k; - signed int ti, tj, tk; - int rval; - struct cell_candidate_list cl; - - cl.cand = calloc(MAX_REAX_CELL_CANDIDATES, - sizeof(struct cell_candidate)); - if ( cl.cand == NULL ) { - ERROR("Failed to allocate cell candidate list.\n"); - return 0; - } - cl.n_cand = 0; - - /* Find candidates for axes 0 and 1 which have the right angle */ - for ( i=0; isearch[0].n_cand; i++ ) { - for ( j=0; jsearch[1].n_cand; j++ ) { - for ( k=0; ksearch[2].n_cand; k++ ) { - for ( ti=-1; ti<=1; ti+=2 ) { - for ( tj=-1; tj<=1; tj+=2 ) { - for ( tk=-1; tk<=1; tk+=2 ) { - - struct dvec vi, vj, vk; - struct rvec ai, bi, ci; - UnitCell *cnew; - double fom; - - vi = s->search[0].cand[i].v; - vj = s->search[1].cand[j].v; - vk = s->search[2].cand[k].v; - - vi.x *= ti; vi.y *= ti; vi.z *= ti; - vj.x *= tj; vj.y *= tj; vj.z *= tj; - vk.x *= tk; vk.y *= tk; vk.z *= tk; - - if ( !check_vector_combination(&vi, &vj, &vk, cell) ) continue; - - ai.u = vi.x; ai.v = vi.y; ai.w = vi.z; - bi.u = vj.x; bi.v = vj.y; bi.w = vj.z; - ci.u = vk.x; ci.v = vk.y; ci.w = vk.z; - - if ( !right_handed_vec(ai, bi, ci) ) continue; - - /* We have three vectors with the right angles */ - cnew = cell_new_from_direct_axes(ai, bi, ci); - - if ( twinned(cnew, &cl) ) { - cell_free(cnew); - continue; - } - - fom = 1.0; - add_cell_candidate(&cl, cnew, fom); - - } - } - } - } - } - } - - for ( i=0; i aA/10.0 ) w = 1; - if ( fabs(b - bA) > bA/10.0 ) w = 1; - if ( fabs(c - cA) > cA/10.0 ) w = 1; - if ( fabs(al - alA) > deg2rad(5.0) ) w = 1; - if ( fabs(be - beA) > deg2rad(5.0) ) w = 1; - if ( fabs(ga - gaA) > deg2rad(5.0) ) w = 1; - if ( w ) { - STATUS("This cell is a long way from that sought:\n"); - cell_print(cl.cand[i].cell); - } - } - - rval = 0; - for ( i=0; inel*sizeof(double)); - fft_out = fftw_malloc((p->nel/2 + 1)*sizeof(fftw_complex)); - - pmax = max_feature_resolution(image->features); - - /* Sanity check */ - if ( pmax < 1e4 ) { - fftw_free(fft_in); - fftw_free(fft_out); - return 0; - } - - s = search_all_axes(p->cell, pmax); - find_candidates(p, image->features, pmax, fft_in, fft_out, s, - NULL, image->det); - -// refine_all_rigid_groups(image, image->candidate_cells[0], pmax, -// fft_in, fft_out, p->plan, smin, smax, -// image->det, p); - - rval = assemble_cells_from_candidates(image, s, p->cell, p); - - for ( i=0; in_search; i++ ) { - free(s->search[i].cand); - } - free(s->search); - free(s); - fftw_free(fft_in); - fftw_free(fft_out); - return rval; -} - - -IndexingPrivate *reax_prepare(IndexingMethod *indm, UnitCell *cell, - struct detector *det, float *ltl) -{ - struct reax_private *p; - int samp; - double th; - - if ( cell == NULL ) { - ERROR("ReAx needs a unit cell.\n"); - return NULL; - } - - p = calloc(1, sizeof(*p)); - if ( p == NULL ) return NULL; - - p->cell = cell; - - /* Flags that ReAx knows about */ - *indm &= INDEXING_METHOD_MASK | INDEXING_CHECK_PEAKS; - - /* Flags that ReAx requires */ - *indm |= INDEXING_USE_LATTICE_TYPE; - *indm |= INDEXING_USE_CELL_PARAMETERS; - - p->angular_inc = deg2rad(1.0); - - /* Reserve memory, over-estimating the number of directions */ - samp = 2.0*M_PI / p->angular_inc; - p->directions = malloc(samp*samp*sizeof(struct dvec)); - if ( p == NULL) { - free(p); - return NULL; - } - STATUS("Allocated space for %i directions\n", samp*samp); - - /* Generate vectors for 1D Fourier transforms */ - fesetround(1); /* Round to nearest */ - p->n_dir = 0; - for ( th=0.0; thangular_inc ) { - - double ph, phstep, n_phstep; - - n_phstep = 2.0*M_PI*sin(th)/p->angular_inc; - n_phstep = nearbyint(n_phstep); - phstep = 2.0*M_PI/n_phstep; - - for ( ph=0.0; ph<2.0*M_PI; ph+=phstep ) { - - struct dvec *dir; - - assert(p->n_dirdirections[p->n_dir++]; - - dir->x = cos(ph) * sin(th); - dir->y = sin(ph) * sin(th); - dir->z = cos(th); - dir->th = th; - dir->ph = ph; - - } - - } - STATUS("Generated %i directions (angular increment %.3f deg)\n", - p->n_dir, rad2deg(p->angular_inc)); - - p->nel = 1024; - - /* These arrays are not actually used */ - p->fft_in = fftw_malloc(p->nel*sizeof(double)); - p->fft_out = fftw_malloc((p->nel/2 + 1)*sizeof(fftw_complex)); - - p->plan = fftw_plan_dft_r2c_1d(p->nel, p->fft_in, p->fft_out, - FFTW_MEASURE); - - p->cw = 128; p->ch = 128; - - /* Also not used */ - p->r_fft_in = fftw_malloc(p->cw*p->ch*sizeof(fftw_complex)); - p->r_fft_out = fftw_malloc(p->cw*p->ch*sizeof(fftw_complex)); - - p->r_plan = fftw_plan_dft_2d(p->cw, p->ch, p->r_fft_in, p->r_fft_out, - 1, FFTW_MEASURE); - - p->indm = *indm; - - return (IndexingPrivate *)p; -} - - -void reax_cleanup(IndexingPrivate *pp) -{ - struct reax_private *p; - - p = (struct reax_private *)pp; - - free(p->directions); - - fftw_destroy_plan(p->plan); - fftw_free(p->fft_in); - fftw_free(p->fft_out); - - fftw_destroy_plan(p->r_plan); - fftw_free(p->r_fft_in); - fftw_free(p->r_fft_out); - - free(p); -} -- cgit v1.2.3