From a80c4c5d0fbc26062084dde2c380ae5370dcd600 Mon Sep 17 00:00:00 2001 From: Alexandra Tolstikova Date: Wed, 11 Feb 2015 11:52:53 +0100 Subject: 'asdf' indexing method added --- libcrystfel/Makefile.am | 4 +- libcrystfel/src/asdf.c | 1045 +++++++++++++++++++++++++++++++++++++++++++++++ libcrystfel/src/asdf.h | 28 ++ libcrystfel/src/index.c | 29 +- libcrystfel/src/index.h | 4 + 5 files changed, 1101 insertions(+), 9 deletions(-) create mode 100644 libcrystfel/src/asdf.c create mode 100644 libcrystfel/src/asdf.h (limited to 'libcrystfel') diff --git a/libcrystfel/Makefile.am b/libcrystfel/Makefile.am index 2cdfba10..79699e4a 100644 --- a/libcrystfel/Makefile.am +++ b/libcrystfel/Makefile.am @@ -10,7 +10,7 @@ libcrystfel_la_SOURCES = src/reflist.c src/utils.c src/cell.c src/detector.c \ src/render.c src/index.c src/dirax.c src/mosflm.c \ src/cell-utils.c src/integer_matrix.c src/crystal.c \ src/grainspotter.c src/xds.c src/integration.c \ - src/histogram.c src/events.c + src/histogram.c src/events.c src/asdf.c if HAVE_FFTW libcrystfel_la_SOURCES += src/reax.c @@ -30,7 +30,7 @@ libcrystfel_la_include_HEADERS = ${top_srcdir}/version.h \ src/integer_matrix.h src/crystal.h \ src/grainspotter.h src/xds.h \ src/integration.h src/histogram.h \ - src/events.h + src/events.h src/asdf.h AM_CPPFLAGS = -DDATADIR=\""$(datadir)"\" -I$(top_builddir)/lib -Wall AM_CPPFLAGS += -I$(top_srcdir)/lib @LIBCRYSTFEL_CFLAGS@ diff --git a/libcrystfel/src/asdf.c b/libcrystfel/src/asdf.c new file mode 100644 index 00000000..0735a14f --- /dev/null +++ b/libcrystfel/src/asdf.c @@ -0,0 +1,1045 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "image.h" +#include "dirax.h" +#include "utils.h" +#include "peaks.h" +#include "cell-utils.h" +#include "asdf.h" + + +struct asdf_private { + IndexingMethod indm; + float *ltl; + UnitCell *template; +}; + + +/* Possible direct vector */ +struct tvector { + gsl_vector *t; + int n; // number of fitting reflections + int *fits; +}; + +struct asdf_cell { + gsl_vector *axes[3]; + gsl_vector *reciprocal[3]; + + int n; // number of fitting reflections + double volume; + + int N_refls; // total number of reflections + int *reflections; // reflections[i] = 1 if reflections fits + double **indices; // indices[i] = [h, k, l] for all reflections (not rounded) + + int acl; // minimum number of reflections fitting to one of the axes[] + int n_max; // maximum number of reflections fitting to some t-vector +}; + +struct tvector tvector_new(int n) { + struct tvector t; + + t.t = gsl_vector_alloc(3); + t.n = 0; + t.fits = malloc(sizeof(int) * n); + + return t; +} + +static int tvector_free(struct tvector t) { + gsl_vector_free(t.t); + free(t.fits); + + return 1; +} + +static int tvector_memcpy(struct tvector *dest, struct tvector *src, int n) { + gsl_vector_memcpy(dest->t, src->t); + dest->n = src->n; + memcpy(dest->fits, src->fits, sizeof(int) * n); + + return 1; +} + +static int asdf_cell_free(struct asdf_cell c) { + int i; + for (i = 0; i < 3; i ++) { + gsl_vector_free(c.axes[i]); + gsl_vector_free(c.reciprocal[i]); + } + + free(c.reflections); + for (i = 0; i < c.N_refls; i ++) { + free(c.indices[i]); + } + free(c.indices); + + return 1; +} + +static struct asdf_cell asdf_cell_new(int n) { + + struct asdf_cell c; + + int i; + for (i = 0; i < 3; i ++) { + c.axes[i] = gsl_vector_alloc(3); + c.reciprocal[i] = gsl_vector_alloc(3); + } + + c.N_refls = n; + c.reflections = malloc(sizeof(int) * n); + + c.indices = malloc(sizeof(int *) * n); + for (i = 0; i < n; i ++) { + c.indices[i] = malloc(sizeof(int) * 3); + } + + c.n = 0; + + c.acl = 0; + c.n_max = 0; + + return c; +} + +static int asdf_cell_memcpy(struct asdf_cell *dest, struct asdf_cell *src) { + int i; + for (i = 0; i < 3; i ++) { + gsl_vector_memcpy(dest->axes[i], src->axes[i]); + gsl_vector_memcpy(dest->reciprocal[i], src->reciprocal[i]); + } + + dest->volume = src->volume; + + int n = src->N_refls; + dest->N_refls = n; + + dest->n = src->n; + memcpy(dest->reflections, src->reflections, sizeof(int) * n); + + memcpy(dest->indices, src->indices, sizeof(int *) * n); + + for (i = 0; i < n; i ++) { + memcpy(dest->indices[i], src->indices[i], sizeof(int) * 3); + } + + dest->acl = src->acl; + dest->n_max = src->n_max; + return 1; +} + +/* result = vec1 cross vec2 */ +static int cross_product(gsl_vector *vec1, gsl_vector *vec2, gsl_vector **result) { + double c1[3], c2[3], p[3]; + int i; + for (i = 0; i < 3; i += 1) { + c1[i] = gsl_vector_get(vec1, i); + c2[i] = gsl_vector_get(vec2, i); + } + + p[0] = c1[1] * c2[2] - c1[2] * c2[1]; + p[1] = - c1[0] * c2[2] + c1[2] * c2[0]; + p[2] = c1[0] * c2[1] - c1[1] * c2[0]; + + for (i = 0; i < 3; i += 1) { + gsl_vector_set(*result, i, p[i]); + } + + return 1; +} + +/* Returns triple product of three gsl_vectors */ +static double calc_volume(gsl_vector *vec1, gsl_vector *vec2, gsl_vector *vec3) { + double volume; + gsl_vector *cross = gsl_vector_alloc(3); + + cross_product(vec1, vec2, &cross); + gsl_blas_ddot(vec3, cross, &volume); + + gsl_vector_free(cross); + return volume; +} + +static int calc_reciprocal(gsl_vector **direct, gsl_vector **reciprocal) { + double volume; + + cross_product(direct[1], direct[2], &reciprocal[0]); + gsl_blas_ddot(direct[0], reciprocal[0], &volume); + gsl_vector_scale(reciprocal[0], 1/volume); + + cross_product(direct[2], direct[0], &reciprocal[1]); + gsl_vector_scale(reciprocal[1], 1/volume); + + cross_product(direct[0], direct[1], &reciprocal[2]); + gsl_vector_scale(reciprocal[2], 1/volume); + + return 1; +} + +static int check_cell(struct asdf_private *dp, struct image *image, + UnitCell *cell) +{ + UnitCell *out; + Crystal *cr; + + if ( dp->indm & INDEXING_CHECK_CELL_COMBINATIONS ) { + + out = match_cell(cell, dp->template, 0, dp->ltl, 1); + if ( out == NULL ) return 0; + + } else if ( dp->indm & INDEXING_CHECK_CELL_AXES ) { + + out = match_cell(cell, dp->template, 0, dp->ltl, 0); + if ( out == NULL ) return 0; + + } else { + out = cell_new_from_cell(cell); + } + + cr = crystal_new(); + if ( cr == NULL ) { + ERROR("Failed to allocate crystal.\n"); + return 0; + } + + crystal_set_cell(cr, out); + + if ( dp->indm & INDEXING_CHECK_PEAKS ) { + if ( !peak_sanity_check(image, &cr, 1) ) { + crystal_free(cr); /* Frees the cell as well */ + cell_free(out); + return 0; + } + } + + image_add_crystal(image, cr); + + return 1; +} + +static int compare_doubles (const void *a, const void *b) +{ + const double *da = (const double *) a; + const double *db = (const double *) b; + + return (*da > *db) - (*da < *db); +} + +/* Compares tvectors by length */ +static int compare_tvectors (const void *a, const void *b) +{ + struct tvector *ta = (struct tvector *) a; + struct tvector *tb = (struct tvector *) b; + + //~ if (ta->n == tb->n) { + return (gsl_blas_dnrm2(ta->t) > gsl_blas_dnrm2(tb->t)) - + (gsl_blas_dnrm2(ta->t) < gsl_blas_dnrm2(tb->t)); + //~ } + //~ + //~ return (ta->n > tb->n) - (ta->n < tb->n); +} + +/* Calculates normal to a triplet c1, c2, c3. Returns 0 if reflections are on + * the same line */ +static int calc_normal(gsl_vector *c1, gsl_vector *c2, gsl_vector *c3, + gsl_vector *normal) { + + gsl_vector *c12 = gsl_vector_alloc(3); + gsl_vector *c23 = gsl_vector_alloc(3); + gsl_vector *c31 = gsl_vector_alloc(3); + gsl_vector *res = gsl_vector_alloc(3); + + cross_product(c1, c2, &c12); + cross_product(c2, c3, &c23); + cross_product(c3, c1, &c31); + + int i; + for (i = 0; i < 3; i += 1) { + gsl_vector_set(res, i, gsl_vector_get(c12, i) + + gsl_vector_get(c23, i) + + gsl_vector_get(c31, i)); + } + + gsl_vector_free(c12); + gsl_vector_free(c23); + gsl_vector_free(c31); + + double norm = gsl_blas_dnrm2(res); + if (norm < 0.0001) { + gsl_vector_free(res); + return 0; + } else { + gsl_vector_scale(res, 1/norm); + gsl_vector_memcpy(normal, res); + gsl_vector_free(res); + } + + return 1; +} + +static float find_ds_fft(double *projections, int N_projections, double d_max) { + + int n = N_projections; + double projections_sorted[n]; + memcpy(projections_sorted, projections, sizeof(double) * n); + qsort(projections_sorted, n, sizeof(double), compare_doubles); + + int i; + + int N = 1024; // number of points in fft calculation + + double *in = malloc(sizeof(double) * N); + fftw_complex *out = fftw_malloc(sizeof (fftw_complex) * N); + + for (i = 0; i < N; i++) { + in[i] = 0; + } + + for (i = 0; i < n; i++) { + in[(int)((projections_sorted[i] - projections_sorted[0]) / + (projections_sorted[n - 1] - projections_sorted[0]) * N)] ++; + } + + fftw_plan p = fftw_plan_dft_r2c_1d(N, in, out, FFTW_MEASURE); + fftw_execute(p); + fftw_destroy_plan(p); + + int i_max = (int)(d_max * (projections_sorted[n - 1] - + projections_sorted[0])); + + int d = 1; + double max = 0; + double a; + for (i = 1; i <= i_max; i ++) { + a = sqrt(out[i][0] * out[i][0] + out[i][1] * out[i][1]); + if (a > max) { + max = a; + d = i; + } + } + + double ds = (projections_sorted[n - 1] - projections_sorted[0]) / d; + + free(in); + fftw_free(out); + return ds; +} + +/* Returns number of reflections fitting ds. + * A projected reflection fits a one-dimensional lattice with elementary + * lattice vector d* if its absolute distance to the nearest lattice + * point is less than LevelFit. */ +static int check_refl_fitting_ds(double *projections, int N_projections, + double ds, double LevelFit) { + + if (ds == 0) return 0; + + int i; + int n = 0; + for (i = 0; i < N_projections; i += 1) { + if (fabs(projections[i] - ds * round(projections[i]/ds)) < LevelFit) { + n += 1; + } + } + + return n; +} + +/* Refines d*, writes 1 to fits[i] if the i-th projection fits d* */ +static float refine_ds(double *projections, int N_projections, double ds, + double LevelFit, int *fits) { + + double fit_refls[N_projections]; + double indices[N_projections]; + + int i; + + int N_fits = 0; + int N_new = N_projections; + + double c1, cov11, sumsq; + double ds_new = ds; + while (N_fits < N_new) { + N_fits = 0; + for (i = 0; i < N_projections; i ++) { + if (fabs(projections[i] - ds_new * round(projections[i]/ds_new)) < + LevelFit) { + fit_refls[N_fits] = projections[i]; + indices[N_fits] = round(projections[i]/ds_new); + N_fits ++; + fits[i] = 1; + } else { + fits[i] = 0; + } + } + + + gsl_fit_mul(indices, 1, fit_refls, 1, N_fits, &c1, &cov11, &sumsq); + N_new = check_refl_fitting_ds(projections, N_projections, c1, LevelFit); + if (N_new >= N_fits) { + ds_new = c1; + } + } + + return ds_new; +} + +static int check_refl_fitting_cell(struct asdf_cell *c, + gsl_vector **reflections, + int N_reflections, double IndexFit) { + + double dist[3]; + + calc_reciprocal(c->axes, c->reciprocal); + c->n = 0; + int i, j, k; + for(i = 0; i < N_reflections; i += 1) { + + for (j = 0; j < 3; j += 1) dist[j] = 0; + + for (j = 0; j < 3; j += 1) { + gsl_blas_ddot(reflections[i], c->axes[j], &c->indices[i][j]); + + for (k = 0; k < 3; k += 1) { + dist[k] += gsl_vector_get(c->reciprocal[j], k) * + (c->indices[i][j] - round(c->indices[i][j])); + } + } + + /* A reflection fits if the distance (in reciprocal space) between the + * observed and calculated reflection position is less than Indexfit */ + + if (dist[0]*dist[0] + dist[1]*dist[1] + dist[2]*dist[2] < + IndexFit * IndexFit ) { + c->reflections[i] = 1; + c->n += 1; + } else { + c->reflections[i] = 0; + } + } + + return 1; +} + +static void print_asdf_cell(struct asdf_cell cc) { + double a, b, c, alpha, beta, gamma, ab, bc, ca; + + a = gsl_blas_dnrm2(cc.axes[0]); + b = gsl_blas_dnrm2(cc.axes[1]); + c = gsl_blas_dnrm2(cc.axes[2]); + + gsl_blas_ddot(cc.axes[0], cc.axes[1], &ab); + gsl_blas_ddot(cc.axes[1], cc.axes[2], &bc); + gsl_blas_ddot(cc.axes[0], cc.axes[2], &ca); + + alpha = acos(bc/b/c)/M_PI*180; + beta = acos(ca/a/c)/M_PI*180; + gamma = acos(ab/a/b)/M_PI*180; + + //~ int i, j; + //~ for (i = 0; i < 3; i ++) { + //~ for (j = 0; j < 3; j ++) { + //~ printf("%f ", gsl_vector_get(cc.axes[i], j)); + //~ } + //~ printf("\n"); + //~ } + + printf("%.2f %.2f %.2f %.2f %.2f %.2f %.0f %d \n", a, b, c, + alpha, beta, gamma, + cc.volume, cc.n); + +} + +static int refine_asdf_cell(struct asdf_cell c, gsl_vector **reflections, + int N_reflections, double IndexFit) { + + gsl_matrix *X = gsl_matrix_alloc(c.n, 3); + + gsl_vector *r[] = {gsl_vector_alloc(c.n), + gsl_vector_alloc(c.n), + gsl_vector_alloc(c.n)}; + + gsl_vector *res = gsl_vector_alloc(3); + gsl_matrix *cov = gsl_matrix_alloc (3, 3); + double chisq; + + int i, j; + int n = 0; + for (i = 0; i < N_reflections; i += 1) { + if (c.reflections[i] == 1) { + for (j = 0; j < 3; j += 1) { + gsl_matrix_set(X, n, j, round(c.indices[i][j])); + gsl_vector_set(r[j], n, gsl_vector_get(reflections[i], j)); + } + n += 1; + } + } + + gsl_multifit_linear_workspace *work = gsl_multifit_linear_alloc(c.n, 3); + + for (i = 0; i < 3; i += 1) { + gsl_multifit_linear (X, r[i], res, cov, &chisq, work); + + for (j = 0; j < 3; j += 1) { + gsl_vector_set(c.reciprocal[j], i, gsl_vector_get(res, j)); + } + } + + calc_reciprocal(c.reciprocal, c.axes); + + double a[3]; + for (i = 0; i < 3; i ++) { + a[i] = gsl_blas_dnrm2(c.axes[i]); + } + + gsl_multifit_linear_free(work); + gsl_vector_free(res); + gsl_matrix_free(cov); + gsl_matrix_free(X); + for (i = 0; i < 3; i += 1) { + gsl_vector_free(r[i]); + } + + if (fabs(a[0]) > 10000 || fabs(a[1]) > 10000 || fabs(a[2]) > 10000 + || isnan(a[0])) { + return 0; + } + + return 1; +} + +static int reduce_asdf_cell(struct asdf_cell *cl) { + double a, b, c, alpha, beta, gamma, ab, bc, ca, bb, cc; + + gsl_vector *va = gsl_vector_alloc(3); + + gsl_vector *vb = gsl_vector_alloc(3); + + gsl_vector *vc = gsl_vector_alloc(3); + + int changed = 1; + while (changed) { + changed = 0; + + gsl_vector_memcpy(va, cl->axes[0]); + gsl_vector_memcpy(vb, cl->axes[1]); + gsl_vector_memcpy(vc, cl->axes[2]); + + a = gsl_blas_dnrm2(va); + b = gsl_blas_dnrm2(vb); + c = gsl_blas_dnrm2(vc); + + gsl_blas_ddot(va, vb, &ab); + gsl_blas_ddot(vb, vc, &bc); + gsl_blas_ddot(vc, va, &ca); + + alpha = acos(bc/b/c)/M_PI*180; + beta = acos(ca/a/c)/M_PI*180; + gamma = acos(ab/a/b)/M_PI*180; + + if (changed == 0) { + + if (gamma < 90) { + gsl_vector_scale(vb, -1); + gamma = 180 - gamma; + alpha = 180 - alpha; + } + + gsl_vector_add(vb, va); + bb = gsl_blas_dnrm2(vb); + if (bb < b) { + b = bb; + if (a < b) { + gsl_vector_memcpy(cl->axes[1], vb); + } else { + gsl_vector_memcpy(cl->axes[1], va); + gsl_vector_memcpy(cl->axes[0], vb); + } + changed = 1; + } + } + + if (changed == 0) { + + if (beta < 90) { + gsl_vector_scale(vc, -1); + beta = 180 - beta; + alpha = 180 - alpha; + } + + gsl_vector_add(vc, va); + cc = gsl_blas_dnrm2(vc); + if (cc < c) { + c = cc; + if (b < c) { + gsl_vector_memcpy(cl->axes[2], vc); + } else if (a < c) { + gsl_vector_memcpy(cl->axes[1], vc); + gsl_vector_memcpy(cl->axes[2], vb); + } else { + gsl_vector_memcpy(cl->axes[0], vc); + gsl_vector_memcpy(cl->axes[1], va); + gsl_vector_memcpy(cl->axes[2], vb); + } + changed = 1; + } + } + + if (changed == 0) { + if (alpha < 90) { + gsl_vector_scale(vc, -1); + beta = 180 - beta; + alpha = 180 - alpha; + } + + gsl_vector_add(vc, vb); + cc = gsl_blas_dnrm2(vc); + if (cc < c) { + c = cc; + if (b < c) { + gsl_vector_memcpy(cl->axes[2], vc); + } else if (a < c) { + gsl_vector_memcpy(cl->axes[1], vc); + gsl_vector_memcpy(cl->axes[2], vb); + } else { + gsl_vector_memcpy(cl->axes[0], vc); + gsl_vector_memcpy(cl->axes[1], va); + gsl_vector_memcpy(cl->axes[2], vb); + } + changed = 1; + } + } + } + + cross_product(cl->axes[0], cl->axes[1], &vc); + gsl_blas_ddot(vc, cl->axes[2], &cl->volume); + if (cl->volume < 0) { + gsl_vector_scale(cl->axes[2], -1); + cl->volume *= -1; + } + + gsl_vector_free(va); + gsl_vector_free(vb); + gsl_vector_free(vc); + + return 1; +} + +static int check_cell_angles(gsl_vector *va, gsl_vector *vb, gsl_vector *vc, + double max_cos) { + double a, b, c, cosa, cosb, cosg, ab, bc, ca; + + a = gsl_blas_dnrm2(va); + b = gsl_blas_dnrm2(vb); + c = gsl_blas_dnrm2(vc); + + gsl_blas_ddot(va, vb, &ab); + gsl_blas_ddot(vb, vc, &bc); + gsl_blas_ddot(vc, va, &ca); + + cosa = bc/b/c; + cosb = ca/a/c; + cosg = ab/a/b; + + if (fabs(cosa) > max_cos || fabs(cosb) > max_cos || fabs(cosg) > max_cos) { + return 0; + } + + return 1; +} + +/* Returns min(t1.n, t2.n, t3.n) */ +static int find_acl(struct tvector t1, struct tvector t2, struct tvector t3) { + int i = t1.n, j = t2.n, k = t3.n; + if (i <= j && i <= k) return i; + if (j <= i && j <= k) return j; + if (k <= i && k <= j) return k; +} + +static int find_cell(struct tvector *tvectors, int N_tvectors, double IndexFit, + double volume_min, double volume_max, int n_max, + gsl_vector **reflections, int N_reflections, + struct asdf_cell *result) { + + int i, j, k, n, cell_correct; + + double volume; + + // only tvectors with number of fitting reflections > acl are considered + int acl = N_reflections < 18 ? 6 : N_reflections/3; + + struct asdf_cell c = asdf_cell_new(N_reflections); + + // traversing a 3d array in slices perpendicular to the main diagonal + int sl; + for (sl = 0; sl < 3 * N_tvectors - 1; sl++) { + + int i_min = sl < 2 * N_tvectors ? 0 : sl - 2 * N_tvectors; + int i_max = sl < N_tvectors ? sl : N_tvectors; + + for (i = i_min; i < i_max; i++) if (tvectors[i].n > acl) { + + int j_min = sl - N_tvectors - 2 * i - 1 < 0 ? i + 1 : sl - N_tvectors - i; + int j_max = sl - N_tvectors - i < 0 ? sl - i : N_tvectors; + + + for (j = j_min; j < j_max; j++) if (tvectors[j].n > acl) { + + k = sl - i - j - 1; + + if (k > j && tvectors[k].n > acl && + check_cell_angles(tvectors[i].t, + tvectors[j].t, + tvectors[k].t, 0.99) ) { + + volume = calc_volume(tvectors[i].t, + tvectors[j].t, + tvectors[k].t); + + if (fabs(volume) > volume_min && fabs(volume) < volume_max){ + + gsl_vector_memcpy(c.axes[0], tvectors[i].t); + gsl_vector_memcpy(c.axes[1], tvectors[j].t); + gsl_vector_memcpy(c.axes[2], tvectors[k].t); + + c.volume = volume; + check_refl_fitting_cell(&c, reflections, N_reflections, + IndexFit); + + if (c.n >= 6) { + reduce_asdf_cell(&c); + + // if one of the cell angles > 135 or < 45 do not continue + if (!check_cell_angles(c.axes[0], c.axes[1], + c.axes[2], 0.71)) break; + + // index reflections with new cell axes + check_refl_fitting_cell(&c, reflections, + N_reflections, IndexFit); + + acl = find_acl(tvectors[i], + tvectors[j], + tvectors[k]); + + c.acl = acl; + c.n_max = n_max; + + // refine cell until number of fitting reflections + // stops increasing + n = 0; + cell_correct = 1; + while (c.n - n && cell_correct) { + n = c.n; + cell_correct = refine_asdf_cell(c, reflections, + N_reflections, + IndexFit); + + check_refl_fitting_cell(&c, reflections, + N_reflections, IndexFit); + } + + if (cell_correct) { + reduce_asdf_cell(&c); + if (result->n < c.n) + asdf_cell_memcpy(result, &c); + + acl++; + + if (acl > n_max) break; + if (tvectors[j].n <= acl || + tvectors[i].n <= acl) break; + } + } + } + } + } + if (acl > n_max) break; + if (tvectors[i].n <= acl) break; + } + if (acl > n_max) break; + } + + asdf_cell_free(c); + + if (result->n) return 1; + return 0; +} + +static void shuffle_triplets(int **triplets, int n) { + int i, j; + int t[3]; + for (i = 0; i < n - 1; i++) { + j = i + rand() / (RAND_MAX / (n - i) + 1); + memcpy(t, triplets[j], 3 * sizeof(int)); + memcpy(triplets[j], triplets[i], 3 * sizeof(int)); + memcpy(triplets[i], t, 3 * sizeof(int)); + } +} + +static double angle_between_gsl(gsl_vector *a, gsl_vector *b) { + double ab; + gsl_blas_ddot(a, b, &ab); + return acos(ab/gsl_blas_dnrm2(a)/gsl_blas_dnrm2(b)) * 180 / M_PI; +} + +static int index_refls(gsl_vector **reflections, int N_reflections, double d_max, + double volume_min, double volume_max, double LevelFit, + double IndexFit, int i_max, struct asdf_cell *c) { + + int i, j, k, l, n; + + /* Number of triplets = c_n^3 if n - number of reflections */ + int N_triplets = N_reflections * (N_reflections - 1) * + (N_reflections - 2) / 6; + + int **triplets = malloc(N_triplets * sizeof(int *)); + l = 0; + for (i = 0; i < N_reflections; i++) { + for (j = i + 1; j < N_reflections; j++) { + for (k = j + 1; k < N_reflections; k++) { + triplets[l] = malloc(3 * sizeof(int)); + + triplets[l][0] = i; + triplets[l][1] = j; + triplets[l][2] = k; + l++; + } + } + } + + /* Triplets are processed in a random sequence if N_triplets > 10000 */ + if (N_reflections > 40) shuffle_triplets(triplets, N_triplets); + + gsl_vector *normal = gsl_vector_alloc(3); + + double projections[N_reflections]; + double ds; + + int *fits = malloc(N_reflections * sizeof(int)); + + if (i_max > N_triplets) i_max = N_triplets; + + struct tvector *tvectors = malloc(i_max * sizeof(struct tvector)); + int N_tvectors = 0; + + int n_max = 0; // maximum number of reflections fitting one of tvectors + + for (i = 0; i < i_max; i++) { + if (calc_normal(reflections[triplets[i][0]], + reflections[triplets[i][1]], + reflections[triplets[i][2]], + normal)) { + + // calculate projections of reflections to normal + for (k = 0; k < N_reflections; k++) { + gsl_blas_ddot(normal, reflections[k], &projections[k]); + } + + // find ds - period in 1d lattice of projections + ds = find_ds_fft(projections, N_reflections, d_max); + + // refine ds, write 1 to fits[i] if reflections[i] fits ds + ds = refine_ds(projections, N_reflections, ds, LevelFit, fits); + + // n - number of reflections fitting ds + n = check_refl_fitting_ds(projections, N_reflections, ds, LevelFit); + + // normal/ds - possible direct vector + gsl_vector_scale(normal, 1/ds); + + if (n > N_reflections/3 && n > 6) { + + tvectors[N_tvectors] = tvector_new(N_reflections); + + gsl_vector_memcpy(tvectors[N_tvectors].t, normal); + memcpy(tvectors[N_tvectors].fits, fits, + N_reflections * sizeof(int)); + + tvectors[N_tvectors].n = n; + + N_tvectors++; + + if (n > n_max) n_max = n; + } + } + + if ((i != 0 && i % 10000 == 0) || i == i_max - 1) { + + // sort tvectors by length + qsort(tvectors, N_tvectors, sizeof(struct tvector), + compare_tvectors); + + /* three shortest independent tvectors with t.n > acl determine + * the final cell. acl is selected for the solution with the + * maximum number of fitting reflections */ + + find_cell(tvectors, N_tvectors, IndexFit, volume_min, volume_max, + n_max, reflections, N_reflections, c); + + if (c->n > 4 * n_max / 5) { + break; + } + } + } + free(fits); + + for (i = 0; i < N_tvectors; i++) { + tvector_free(tvectors[i]); + } + free(tvectors); + + for (i = 0; i < N_triplets; i++) { + free(triplets[i]); + } + free(triplets); + + if (c->n) return 1; + + return 0; +} + +double cell_get_volume(UnitCell *cell) +{ + double asx, asy, asz; + double bsx, bsy, bsz; + double csx, csy, csz; + struct rvec aCb; + double volume; + + if ( cell_get_reciprocal(cell, &asx, &asy, &asz, + &bsx, &bsy, &bsz, + &csx, &csy, &csz) ) { + ERROR("Couldn't get cell cartesian.\n"); + return 0; + } + + /* "a" cross "b" */ + aCb.u = asy*bsz - asz*bsy; + aCb.v = - (asx*bsz - asz*bsx); + aCb.w = asx*bsy - asy*bsx; + + /* "a cross b" dot "c" */ + volume = (aCb.u*csx + aCb.v*csy + aCb.w*csz)/1e30; + + return 1/volume; +} + +int run_asdf(struct image *image, IndexingPrivate *ipriv) { + int i; + + double LevelFit = 1./1000; + double IndexFit = 1./500; + double d_max = 220.; // thrice the maximum expected axis length + double volume_min = 100.; + double volume_max = 1000000.; + + int i_max = 10000; // maximum number of triplets + + struct asdf_private *dp = (struct asdf_private *)ipriv; + + if ( dp->indm & INDEXING_CHECK_CELL_AXES ) { + double volume = cell_get_volume(dp->template); + volume_min = volume * 0.95; + volume_max = volume * 1.05; + } + + int N_reflections = image_feature_count(image->features); + gsl_vector *reflections[N_reflections]; + + for (i = 0; i < N_reflections; i++) { + struct imagefeature *f; + + f = image_get_feature(image->features, i); + if ( f == NULL ) continue; + + reflections[i] = gsl_vector_alloc(3); + gsl_vector_set(reflections[i], 0, f->rx/1e10); + gsl_vector_set(reflections[i], 1, f->ry/1e10); + gsl_vector_set(reflections[i], 2, f->rz/1e10); + } + + struct asdf_cell c = asdf_cell_new(N_reflections); + + if (N_reflections == 0) return 0; + + i = index_refls(reflections, N_reflections, d_max, volume_min, volume_max, + LevelFit, IndexFit, i_max, &c); + //~ + + for (i = 0; i < N_reflections; i ++) { + gsl_vector_free(reflections[i]); + } + + if (i) { + UnitCell *uc; + uc = cell_new(); + + cell_set_cartesian(uc, gsl_vector_get(c.axes[0], 0) * 1e-10, + gsl_vector_get(c.axes[0], 1) * 1e-10, + gsl_vector_get(c.axes[0], 2) * 1e-10, + gsl_vector_get(c.axes[1], 0) * 1e-10, + gsl_vector_get(c.axes[1], 1) * 1e-10, + gsl_vector_get(c.axes[1], 2) * 1e-10, + gsl_vector_get(c.axes[2], 0) * 1e-10, + gsl_vector_get(c.axes[2], 1) * 1e-10, + gsl_vector_get(c.axes[2], 2) * 1e-10); + + if ( check_cell(dp, image, uc) ) { + cell_free(uc); + return 1; + } + + cell_free(uc); + } + + return 0; +} + +IndexingPrivate *asdf_prepare(IndexingMethod *indm, UnitCell *cell, + struct detector *det, float *ltl) +{ + struct asdf_private *dp; + int need_cell = 0; + + if ( *indm & INDEXING_CHECK_CELL_COMBINATIONS ) need_cell = 1; + if ( *indm & INDEXING_CHECK_CELL_AXES ) need_cell = 1; + + if ( need_cell && !cell_has_parameters(cell) ) { + ERROR("Altering your asdf flags because cell parameters were" + " not provided.\n"); + *indm &= ~INDEXING_CHECK_CELL_COMBINATIONS; + *indm &= ~INDEXING_CHECK_CELL_AXES; + } + + /* Flags that asdf knows about */ + *indm &= INDEXING_METHOD_MASK | INDEXING_CHECK_CELL_COMBINATIONS + | INDEXING_CHECK_CELL_AXES | INDEXING_CHECK_PEAKS; + + dp = malloc(sizeof(struct asdf_private)); + if ( dp == NULL ) return NULL; + + dp->ltl = ltl; + dp->template = cell; + dp->indm = *indm; + + return (IndexingPrivate *)dp; +} + + +void asdf_cleanup(IndexingPrivate *pp) +{ + struct asdf_private *p; + p = (struct asdf_private *)pp; + free(p); +} diff --git a/libcrystfel/src/asdf.h b/libcrystfel/src/asdf.h new file mode 100644 index 00000000..402636d1 --- /dev/null +++ b/libcrystfel/src/asdf.h @@ -0,0 +1,28 @@ + + +#ifndef ASDF_H +#define ASDF_H + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include "index.h" + +#ifdef __cplusplus +extern "C" { +#endif + +extern int run_asdf(struct image *image, IndexingPrivate *ipriv); + +extern IndexingPrivate *asdf_prepare(IndexingMethod *indm, + UnitCell *cell, struct detector *det, + float *ltl); + +extern void asdf_cleanup(IndexingPrivate *pp); + +#ifdef __cplusplus +} +#endif + +#endif /* DIRAX_H */ diff --git a/libcrystfel/src/index.c b/libcrystfel/src/index.c index c3f7c35d..97187655 100644 --- a/libcrystfel/src/index.c +++ b/libcrystfel/src/index.c @@ -44,6 +44,7 @@ #include "utils.h" #include "peaks.h" #include "dirax.h" +#include "asdf.h" #include "mosflm.h" #include "xds.h" #include "detector.h" @@ -79,12 +80,16 @@ IndexingPrivate **prepare_indexing(IndexingMethod *indm, UnitCell *cell, iprivs[n] = dirax_prepare(&indm[n], cell, det, ltl); break; + case INDEXING_ASDF : + iprivs[n] = asdf_prepare(&indm[n], cell, det, ltl); + break; + case INDEXING_MOSFLM : iprivs[n] = mosflm_prepare(&indm[n], cell, det, ltl); break; case INDEXING_XDS : - iprivs[n] = xds_prepare(&indm[n], cell, det, ltl); + iprivs[n] = xds_prepare(&indm[n], cell, det, ltl); break; case INDEXING_REAX : @@ -150,6 +155,10 @@ void cleanup_indexing(IndexingMethod *indms, IndexingPrivate **privs) dirax_cleanup(privs[n]); break; + case INDEXING_ASDF : + asdf_cleanup(privs[n]); + break; + case INDEXING_MOSFLM : mosflm_cleanup(privs[n]); break; @@ -215,6 +224,10 @@ static int try_indexer(struct image *image, IndexingMethod indm, return run_dirax(image, ipriv); break; + case INDEXING_ASDF : + return run_asdf(image, ipriv); + break; + case INDEXING_MOSFLM : return run_mosflm(image, ipriv); break; @@ -248,11 +261,6 @@ void index_pattern(struct image *image, if ( indms == NULL ) return; - if ( image_feature_count(image->features) > 10000 ) { - STATUS("WARNING: The number of peaks is very large for '%s'.\n", - image->filename); - } - map_all_peaks(image); image->crystals = NULL; image->n_crystals = 0; @@ -349,6 +357,10 @@ char *indexer_str(IndexingMethod indm) strcpy(str, "dirax"); break; + case INDEXING_ASDF : + strcpy(str, "asdf"); + break; + case INDEXING_MOSFLM : strcpy(str, "mosflm"); break; @@ -423,13 +435,16 @@ IndexingMethod *build_indexer_list(const char *str) if ( strcmp(methods[i], "dirax") == 0) { list[++nmeth] = INDEXING_DEFAULTS_DIRAX; + } else if ( strcmp(methods[i], "asdf") == 0) { + list[++nmeth] = INDEXING_DEFAULTS_ASDF; + } else if ( strcmp(methods[i], "mosflm") == 0) { list[++nmeth] = INDEXING_DEFAULTS_MOSFLM; } else if ( strcmp(methods[i], "grainspotter") == 0) { list[++nmeth] = INDEXING_DEFAULTS_GRAINSPOTTER; - } else if ( strcmp(methods[i], "xds") == 0) { + } else if ( strcmp(methods[i], "xds") == 0) { list[++nmeth] = INDEXING_DEFAULTS_XDS; } else if ( strcmp(methods[i], "reax") == 0) { diff --git a/libcrystfel/src/index.h b/libcrystfel/src/index.h index 0e105aeb..709b7507 100644 --- a/libcrystfel/src/index.h +++ b/libcrystfel/src/index.h @@ -41,6 +41,9 @@ #define INDEXING_DEFAULTS_DIRAX (INDEXING_DIRAX | INDEXING_CHECK_PEAKS \ | INDEXING_CHECK_CELL_COMBINATIONS) +#define INDEXING_DEFAULTS_ASDF (INDEXING_ASDF | INDEXING_CHECK_PEAKS \ + | INDEXING_CHECK_CELL_COMBINATIONS) + #define INDEXING_DEFAULTS_MOSFLM (INDEXING_MOSFLM | INDEXING_CHECK_PEAKS \ | INDEXING_CHECK_CELL_COMBINATIONS \ | INDEXING_USE_LATTICE_TYPE) @@ -95,6 +98,7 @@ typedef enum { INDEXING_GRAINSPOTTER = 4, INDEXING_XDS = 5, INDEXING_SIMULATION = 6, + INDEXING_ASDF = 7, /* Bits at the top of the IndexingMethod are flags which modify the * behaviour of the indexer. */ -- cgit v1.2.3 From f8d230451e7e0f71d5bf90d5f1079b89d19c8272 Mon Sep 17 00:00:00 2001 From: Alexandra Tolstikova Date: Wed, 11 Feb 2015 13:44:59 +0100 Subject: asdf.c: codestyle fixed --- libcrystfel/src/asdf.c | 910 +++++++++++++++++++++++++------------------------ 1 file changed, 469 insertions(+), 441 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/asdf.c b/libcrystfel/src/asdf.c index 0735a14f..6ad29645 100644 --- a/libcrystfel/src/asdf.c +++ b/libcrystfel/src/asdf.c @@ -26,24 +26,24 @@ struct asdf_private { /* Possible direct vector */ struct tvector { - gsl_vector *t; - int n; // number of fitting reflections - int *fits; + gsl_vector *t; + int n; // number of fitting reflections + int *fits; }; struct asdf_cell { - gsl_vector *axes[3]; - gsl_vector *reciprocal[3]; - - int n; // number of fitting reflections - double volume; - - int N_refls; // total number of reflections - int *reflections; // reflections[i] = 1 if reflections fits - double **indices; // indices[i] = [h, k, l] for all reflections (not rounded) - - int acl; // minimum number of reflections fitting to one of the axes[] - int n_max; // maximum number of reflections fitting to some t-vector + 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) { @@ -73,13 +73,13 @@ static int tvector_memcpy(struct tvector *dest, struct tvector *src, int n) { static int asdf_cell_free(struct asdf_cell c) { int i; - for (i = 0; i < 3; 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 ++) { + for ( i = 0; i < c.N_refls; i++ ) { free(c.indices[i]); } free(c.indices); @@ -92,7 +92,7 @@ static struct asdf_cell asdf_cell_new(int n) { struct asdf_cell c; int i; - for (i = 0; i < 3; i ++) { + for ( i = 0; i < 3; i++ ) { c.axes[i] = gsl_vector_alloc(3); c.reciprocal[i] = gsl_vector_alloc(3); } @@ -101,7 +101,7 @@ static struct asdf_cell asdf_cell_new(int n) { c.reflections = malloc(sizeof(int) * n); c.indices = malloc(sizeof(int *) * n); - for (i = 0; i < n; i ++) { + for ( i = 0; i < n; i++ ) { c.indices[i] = malloc(sizeof(int) * 3); } @@ -115,7 +115,7 @@ static struct asdf_cell asdf_cell_new(int n) { static int asdf_cell_memcpy(struct asdf_cell *dest, struct asdf_cell *src) { int i; - for (i = 0; i < 3; 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]); } @@ -130,7 +130,7 @@ static int asdf_cell_memcpy(struct asdf_cell *dest, struct asdf_cell *src) { memcpy(dest->indices, src->indices, sizeof(int *) * n); - for (i = 0; i < n; i ++) { + for (i = 0; i < n; i++ ) { memcpy(dest->indices[i], src->indices[i], sizeof(int) * 3); } @@ -140,51 +140,52 @@ static int asdf_cell_memcpy(struct asdf_cell *dest, struct asdf_cell *src) { } /* result = vec1 cross vec2 */ -static int cross_product(gsl_vector *vec1, gsl_vector *vec2, gsl_vector **result) { - double c1[3], c2[3], p[3]; - int i; - for (i = 0; i < 3; i += 1) { - c1[i] = gsl_vector_get(vec1, i); - c2[i] = gsl_vector_get(vec2, i); - } - - p[0] = c1[1] * c2[2] - c1[2] * c2[1]; - p[1] = - c1[0] * c2[2] + c1[2] * c2[0]; - p[2] = c1[0] * c2[1] - c1[1] * c2[0]; - - for (i = 0; i < 3; i += 1) { - gsl_vector_set(*result, i, p[i]); - } - - return 1; +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; + 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; + double volume; + + cross_product(direct[1], direct[2], &reciprocal[0]); + gsl_blas_ddot(direct[0], reciprocal[0], &volume); + gsl_vector_scale(reciprocal[0], 1/volume); + + cross_product(direct[2], direct[0], &reciprocal[1]); + gsl_vector_scale(reciprocal[1], 1/volume); + + cross_product(direct[0], direct[1], &reciprocal[2]); + gsl_vector_scale(reciprocal[2], 1/volume); + + return 1; } static int check_cell(struct asdf_private *dp, struct image *image, @@ -230,21 +231,21 @@ static int check_cell(struct asdf_private *dp, struct image *image, 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); + const double *da = (const double *) a; + const double *db = (const double *) b; + + return (*da > *db) - (*da < *db); } /* Compares tvectors by length */ static int compare_tvectors (const void *a, const void *b) { - struct tvector *ta = (struct tvector *) a; - struct tvector *tb = (struct tvector *) b; + 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 (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); @@ -253,75 +254,75 @@ static int compare_tvectors (const void *a, const void *b) /* 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 *normal) { - gsl_vector *c12 = gsl_vector_alloc(3); - gsl_vector *c23 = gsl_vector_alloc(3); - gsl_vector *c31 = gsl_vector_alloc(3); - gsl_vector *res = gsl_vector_alloc(3); - - cross_product(c1, c2, &c12); - cross_product(c2, c3, &c23); - cross_product(c3, c1, &c31); - - int i; - for (i = 0; i < 3; i += 1) { - gsl_vector_set(res, i, gsl_vector_get(c12, i) + - gsl_vector_get(c23, i) + - gsl_vector_get(c31, i)); - } - - gsl_vector_free(c12); - gsl_vector_free(c23); - gsl_vector_free(c31); - - double norm = gsl_blas_dnrm2(res); - if (norm < 0.0001) { - gsl_vector_free(res); - return 0; - } else { - gsl_vector_scale(res, 1/norm); - gsl_vector_memcpy(normal, res); - gsl_vector_free(res); - } - - return 1; + 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) { - int n = N_projections; - double projections_sorted[n]; - memcpy(projections_sorted, projections, sizeof(double) * n); - qsort(projections_sorted, n, sizeof(double), compare_doubles); - - int i; - - int N = 1024; // number of points in fft calculation - - double *in = malloc(sizeof(double) * N); - fftw_complex *out = fftw_malloc(sizeof (fftw_complex) * N); - - for (i = 0; i < N; i++) { - in[i] = 0; - } - - for (i = 0; i < n; i++) { - in[(int)((projections_sorted[i] - projections_sorted[0]) / - (projections_sorted[n - 1] - projections_sorted[0]) * N)] ++; - } + int n = N_projections; + double projections_sorted[n]; + memcpy(projections_sorted, projections, sizeof(double) * n); + qsort(projections_sorted, n, sizeof(double), compare_doubles); + + int i; + + int N = 1024; // number of points in fft calculation + + double *in = malloc(sizeof(double) * N); + fftw_complex *out = fftw_malloc(sizeof (fftw_complex) * N); + + for ( i = 0; i < N; i++ ) { + in[i] = 0; + } + + for ( i = 0; i < n; i++ ) { + in[(int)((projections_sorted[i] - projections_sorted[0]) / + (projections_sorted[n - 1] - projections_sorted[0]) * N)] ++; + } - fftw_plan p = fftw_plan_dft_r2c_1d(N, in, out, FFTW_MEASURE); - fftw_execute(p); - fftw_destroy_plan(p); - - int i_max = (int)(d_max * (projections_sorted[n - 1] - - projections_sorted[0])); - - int d = 1; - double max = 0; - double a; - for (i = 1; i <= i_max; i ++) { + fftw_plan p = fftw_plan_dft_r2c_1d(N, in, out, FFTW_MEASURE); + fftw_execute(p); + fftw_destroy_plan(p); + + int i_max = (int)(d_max * (projections_sorted[n - 1] - + projections_sorted[0])); + + int d = 1; + double max = 0; + double a; + for ( i = 1; i <= i_max; i++ ) { a = sqrt(out[i][0] * out[i][0] + out[i][1] * out[i][1]); if (a > max) { max = a; @@ -343,17 +344,19 @@ static float find_ds_fft(double *projections, int N_projections, double d_max) { static int check_refl_fitting_ds(double *projections, int N_projections, double ds, double LevelFit) { - if (ds == 0) return 0; + if ( ds == 0 ) return 0; - int i; - int n = 0; - for (i = 0; i < N_projections; i += 1) { - if (fabs(projections[i] - ds * round(projections[i]/ds)) < LevelFit) { + 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; + return n; } /* Refines d*, writes 1 to fits[i] if the i-th projection fits d* */ @@ -370,11 +373,12 @@ static float refine_ds(double *projections, int N_projections, double ds, double c1, cov11, sumsq; double ds_new = ds; - while (N_fits < N_new) { + 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) { + 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 ++; @@ -385,9 +389,11 @@ static float refine_ds(double *projections, int N_projections, double ds, } - 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) { + 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; } } @@ -397,142 +403,146 @@ static float refine_ds(double *projections, int N_projections, double ds, static int check_refl_fitting_cell(struct asdf_cell *c, gsl_vector **reflections, - int N_reflections, double IndexFit) { - - double dist[3]; - - calc_reciprocal(c->axes, c->reciprocal); - c->n = 0; - int i, j, k; - for(i = 0; i < N_reflections; i += 1) { - - for (j = 0; j < 3; j += 1) dist[j] = 0; - - for (j = 0; j < 3; j += 1) { - gsl_blas_ddot(reflections[i], c->axes[j], &c->indices[i][j]); - - for (k = 0; k < 3; k += 1) { - dist[k] += gsl_vector_get(c->reciprocal[j], k) * - (c->indices[i][j] - round(c->indices[i][j])); - } - } - - /* A reflection fits if the distance (in reciprocal space) between the - * observed and calculated reflection position is less than Indexfit */ - - if (dist[0]*dist[0] + dist[1]*dist[1] + dist[2]*dist[2] < - IndexFit * IndexFit ) { - c->reflections[i] = 1; - c->n += 1; - } else { - c->reflections[i] = 0; - } - } - - return 1; + 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; } static void print_asdf_cell(struct asdf_cell cc) { - double a, b, c, alpha, beta, gamma, ab, bc, ca; - - a = gsl_blas_dnrm2(cc.axes[0]); - b = gsl_blas_dnrm2(cc.axes[1]); - c = gsl_blas_dnrm2(cc.axes[2]); - - gsl_blas_ddot(cc.axes[0], cc.axes[1], &ab); - gsl_blas_ddot(cc.axes[1], cc.axes[2], &bc); - gsl_blas_ddot(cc.axes[0], cc.axes[2], &ca); - - alpha = acos(bc/b/c)/M_PI*180; - beta = acos(ca/a/c)/M_PI*180; - gamma = acos(ab/a/b)/M_PI*180; - - //~ int i, j; - //~ for (i = 0; i < 3; i ++) { + double a, b, c, alpha, beta, gamma, ab, bc, ca; + + a = gsl_blas_dnrm2(cc.axes[0]); + b = gsl_blas_dnrm2(cc.axes[1]); + c = gsl_blas_dnrm2(cc.axes[2]); + + gsl_blas_ddot(cc.axes[0], cc.axes[1], &ab); + gsl_blas_ddot(cc.axes[1], cc.axes[2], &bc); + gsl_blas_ddot(cc.axes[0], cc.axes[2], &ca); + + alpha = acos(bc/b/c)/M_PI*180; + beta = acos(ca/a/c)/M_PI*180; + gamma = acos(ab/a/b)/M_PI*180; + + //~ int i, j; + //~ for (i = 0; i < 3; i ++) { //~ for (j = 0; j < 3; j ++) { //~ printf("%f ", gsl_vector_get(cc.axes[i], j)); //~ } //~ printf("\n"); //~ } - printf("%.2f %.2f %.2f %.2f %.2f %.2f %.0f %d \n", a, b, c, - alpha, beta, gamma, - cc.volume, cc.n); - + printf("%.2f %.2f %.2f %.2f %.2f %.2f %.0f %d \n", a, b, c, + alpha, beta, gamma, + cc.volume, cc.n); + } +/* 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 += 1) { - if (c.reflections[i] == 1) { - for (j = 0; j < 3; j += 1) { - gsl_matrix_set(X, n, j, round(c.indices[i][j])); - gsl_vector_set(r[j], n, gsl_vector_get(reflections[i], j)); - } - n += 1; - } - } - - gsl_multifit_linear_workspace *work = gsl_multifit_linear_alloc(c.n, 3); - - for (i = 0; i < 3; i += 1) { - gsl_multifit_linear (X, r[i], res, cov, &chisq, work); - - for (j = 0; j < 3; j += 1) { - gsl_vector_set(c.reciprocal[j], i, gsl_vector_get(res, j)); - } - } - - calc_reciprocal(c.reciprocal, c.axes); + 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 ++) { + for ( i = 0; i < 3; i++ ) { a[i] = gsl_blas_dnrm2(c.axes[i]); } - - gsl_multifit_linear_free(work); - gsl_vector_free(res); - gsl_matrix_free(cov); - gsl_matrix_free(X); - for (i = 0; i < 3; i += 1) { - gsl_vector_free(r[i]); - } - - if (fabs(a[0]) > 10000 || fabs(a[1]) > 10000 || fabs(a[2]) > 10000 - || isnan(a[0])) { + + 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; + 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); + gsl_vector *va = gsl_vector_alloc(3); + gsl_vector *vb = gsl_vector_alloc(3); + gsl_vector *vc = gsl_vector_alloc(3); int changed = 1; - while (changed) { + while ( changed ) { changed = 0; - gsl_vector_memcpy(va, cl->axes[0]); + gsl_vector_memcpy(va, cl->axes[0]); gsl_vector_memcpy(vb, cl->axes[1]); gsl_vector_memcpy(vc, cl->axes[2]); @@ -548,9 +558,9 @@ static int reduce_asdf_cell(struct asdf_cell *cl) { beta = acos(ca/a/c)/M_PI*180; gamma = acos(ab/a/b)/M_PI*180; - if (changed == 0) { + if ( changed == 0 ) { - if (gamma < 90) { + if ( gamma < 90 ) { gsl_vector_scale(vb, -1); gamma = 180 - gamma; alpha = 180 - alpha; @@ -558,9 +568,9 @@ static int reduce_asdf_cell(struct asdf_cell *cl) { gsl_vector_add(vb, va); bb = gsl_blas_dnrm2(vb); - if (bb < b) { + if ( bb < b ) { b = bb; - if (a < b) { + if ( a < b ) { gsl_vector_memcpy(cl->axes[1], vb); } else { gsl_vector_memcpy(cl->axes[1], va); @@ -570,9 +580,9 @@ static int reduce_asdf_cell(struct asdf_cell *cl) { } } - if (changed == 0) { + if ( changed == 0 ) { - if (beta < 90) { + if ( beta < 90 ) { gsl_vector_scale(vc, -1); beta = 180 - beta; alpha = 180 - alpha; @@ -580,11 +590,11 @@ static int reduce_asdf_cell(struct asdf_cell *cl) { gsl_vector_add(vc, va); cc = gsl_blas_dnrm2(vc); - if (cc < c) { + if ( cc < c ) { c = cc; - if (b < c) { + if ( b < c ) { gsl_vector_memcpy(cl->axes[2], vc); - } else if (a < c) { + } else if ( a < c ) { gsl_vector_memcpy(cl->axes[1], vc); gsl_vector_memcpy(cl->axes[2], vb); } else { @@ -596,8 +606,8 @@ static int reduce_asdf_cell(struct asdf_cell *cl) { } } - if (changed == 0) { - if (alpha < 90) { + if ( changed == 0 ) { + if ( alpha < 90 ) { gsl_vector_scale(vc, -1); beta = 180 - beta; alpha = 180 - alpha; @@ -605,11 +615,11 @@ static int reduce_asdf_cell(struct asdf_cell *cl) { gsl_vector_add(vc, vb); cc = gsl_blas_dnrm2(vc); - if (cc < c) { + if ( cc < c ) { c = cc; - if (b < c) { + if ( b < c ) { gsl_vector_memcpy(cl->axes[2], vc); - } else if (a < c) { + } else if ( a < c ) { gsl_vector_memcpy(cl->axes[1], vc); gsl_vector_memcpy(cl->axes[2], vb); } else { @@ -624,7 +634,7 @@ static int reduce_asdf_cell(struct asdf_cell *cl) { cross_product(cl->axes[0], cl->axes[1], &vc); gsl_blas_ddot(vc, cl->axes[2], &cl->volume); - if (cl->volume < 0) { + if ( cl->volume < 0 ) { gsl_vector_scale(cl->axes[2], -1); cl->volume *= -1; } @@ -637,34 +647,36 @@ static int reduce_asdf_cell(struct asdf_cell *cl) { } static int check_cell_angles(gsl_vector *va, gsl_vector *vb, gsl_vector *vc, - double max_cos) { + 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) { + 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; + + 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; + if ( i <= j && i <= k ) return i; + if ( j <= i && j <= k ) return j; + if ( k <= i && k <= j ) return k; } static int find_cell(struct tvector *tvectors, int N_tvectors, double IndexFit, @@ -676,110 +688,118 @@ static int find_cell(struct tvector *tvectors, int N_tvectors, double IndexFit, double volume; - // only tvectors with number of fitting reflections > acl are considered + /* 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); - // 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; + /* Traversing a 3d array in slices perpendicular to the main diagonal */ + int sl; + for ( sl = 0; sl < 3 * N_tvectors - 1; sl++ ) { - 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) ) { + int i_min = sl < 2 * N_tvectors ? 0 : sl - 2 * N_tvectors; + int i_max = sl < N_tvectors ? sl : N_tvectors; + + for ( i = i_min; i < i_max; i++) if (tvectors[i].n > acl ) { + + int j_min = sl - N_tvectors - 2 * i - 1 < 0 ? i + 1 : sl - N_tvectors - i; + int j_max = sl - N_tvectors - i < 0 ? sl - i : N_tvectors; + + for ( j = j_min; j < j_max; j++) if (tvectors[j].n > acl ) { + + k = sl - i - j - 1; + + if ( k > j && tvectors[k].n > acl && + check_cell_angles(tvectors[i].t, + tvectors[j].t, + tvectors[k].t, 0.99) ) + { + + volume = calc_volume(tvectors[i].t, + tvectors[j].t, + tvectors[k].t); + + if ( fabs(volume) > volume_min && + fabs(volume) < volume_max ) + { + + gsl_vector_memcpy(c.axes[0], tvectors[i].t); + gsl_vector_memcpy(c.axes[1], tvectors[j].t); + gsl_vector_memcpy(c.axes[2], tvectors[k].t); + + c.volume = volume; + check_refl_fitting_cell(&c, reflections, + N_reflections, + IndexFit); - volume = calc_volume(tvectors[i].t, - tvectors[j].t, - tvectors[k].t); + if ( c.n < 6 ) break; + + reduce_asdf_cell(&c); + + /* if one of the cell angles > 135 or < 45 + * do not continue */ + if ( !check_cell_angles(c.axes[0], c.axes[1], + c.axes[2], 0.71) ) break; + + /* index reflections with new cell axes */ + check_refl_fitting_cell(&c, reflections, + N_reflections, + IndexFit); + + acl = find_acl(tvectors[i], + tvectors[j], + tvectors[k]); + + c.acl = acl; + c.n_max = n_max; + + /* refine cell until the number of fitting + * reflections stops increasing */ + n = 0; + cell_correct = 1; + while ( c.n - n && cell_correct ) { + n = c.n; + cell_correct = refine_asdf_cell(c, + reflections, + N_reflections, + IndexFit); - if (fabs(volume) > volume_min && fabs(volume) < volume_max){ - - gsl_vector_memcpy(c.axes[0], tvectors[i].t); - gsl_vector_memcpy(c.axes[1], tvectors[j].t); - gsl_vector_memcpy(c.axes[2], tvectors[k].t); - - c.volume = volume; - check_refl_fitting_cell(&c, reflections, N_reflections, - IndexFit); - - if (c.n >= 6) { - reduce_asdf_cell(&c); - - // if one of the cell angles > 135 or < 45 do not continue - if (!check_cell_angles(c.axes[0], c.axes[1], - c.axes[2], 0.71)) break; - - // index reflections with new cell axes - check_refl_fitting_cell(&c, reflections, - N_reflections, IndexFit); - - acl = find_acl(tvectors[i], - tvectors[j], - tvectors[k]); - - c.acl = acl; - c.n_max = n_max; - - // refine cell until number of fitting reflections - // stops increasing - n = 0; - cell_correct = 1; - while (c.n - n && cell_correct) { - n = c.n; - cell_correct = refine_asdf_cell(c, reflections, - N_reflections, - IndexFit); - - check_refl_fitting_cell(&c, reflections, - N_reflections, IndexFit); - } - - if (cell_correct) { - reduce_asdf_cell(&c); - if (result->n < c.n) - asdf_cell_memcpy(result, &c); - - acl++; - - if (acl > n_max) break; - if (tvectors[j].n <= acl || - tvectors[i].n <= acl) break; - } - } + check_refl_fitting_cell(&c, reflections, + N_reflections, + IndexFit); + } + + if ( cell_correct ) { + reduce_asdf_cell(&c); + if ( result->n < c.n ) { + asdf_cell_memcpy(result, &c); } + acl++; + + if (acl > n_max) break; + if (tvectors[j].n <= acl || + tvectors[i].n <= acl) break; } } - if (acl > n_max) break; - if (tvectors[i].n <= acl) break; } - if (acl > n_max) break; + } + 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; + if ( result->n ) return 1; return 0; } static void shuffle_triplets(int **triplets, int n) { int i, j; int t[3]; - for (i = 0; i < n - 1; i++) { + for ( i = 0; i < n - 1; i++ ) { j = i + rand() / (RAND_MAX / (n - i) + 1); memcpy(t, triplets[j], 3 * sizeof(int)); memcpy(triplets[j], triplets[i], 3 * sizeof(int)); @@ -793,9 +813,10 @@ static double angle_between_gsl(gsl_vector *a, gsl_vector *b) { return acos(ab/gsl_blas_dnrm2(a)/gsl_blas_dnrm2(b)) * 180 / M_PI; } -static int index_refls(gsl_vector **reflections, int N_reflections, double d_max, - double volume_min, double volume_max, double LevelFit, - double IndexFit, int i_max, struct asdf_cell *c) { +static int index_refls(gsl_vector **reflections, int N_reflections, + double d_max, double volume_min, double volume_max, + double LevelFit, double IndexFit, int i_max, + struct asdf_cell *c) { int i, j, k, l, n; @@ -805,9 +826,9 @@ static int index_refls(gsl_vector **reflections, int N_reflections, double d_max int **triplets = malloc(N_triplets * sizeof(int *)); l = 0; - for (i = 0; i < N_reflections; i++) { - for (j = i + 1; j < N_reflections; j++) { - for (k = j + 1; k < N_reflections; k++) { + for ( i = 0; i < N_reflections; i++ ) { + for ( j = i + 1; j < N_reflections; j++ ) { + for ( k = j + 1; k < N_reflections; k++ ) { triplets[l] = malloc(3 * sizeof(int)); triplets[l][0] = i; @@ -819,7 +840,7 @@ static int index_refls(gsl_vector **reflections, int N_reflections, double d_max } /* Triplets are processed in a random sequence if N_triplets > 10000 */ - if (N_reflections > 40) shuffle_triplets(triplets, N_triplets); + if ( N_reflections > 40 ) shuffle_triplets(triplets, N_triplets); gsl_vector *normal = gsl_vector_alloc(3); @@ -828,41 +849,47 @@ static int index_refls(gsl_vector **reflections, int N_reflections, double d_max int *fits = malloc(N_reflections * sizeof(int)); - if (i_max > N_triplets) i_max = N_triplets; + if ( i_max > N_triplets ) i_max = N_triplets; - struct tvector *tvectors = malloc(i_max * sizeof(struct tvector)); + struct tvector *tvectors = malloc(i_max * sizeof(struct tvector)); int N_tvectors = 0; int n_max = 0; // maximum number of reflections fitting one of tvectors - for (i = 0; i < i_max; i++) { - if (calc_normal(reflections[triplets[i][0]], - reflections[triplets[i][1]], - reflections[triplets[i][2]], - normal)) { + for ( i = 0; i < i_max; i++ ) { + if ( calc_normal(reflections[triplets[i][0]], + reflections[triplets[i][1]], + reflections[triplets[i][2]], + normal) ) + { - // calculate projections of reflections to normal - for (k = 0; k < N_reflections; k++) { - gsl_blas_ddot(normal, reflections[k], &projections[k]); + /* 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 + /* Find ds - period in 1d lattice of projections */ ds = find_ds_fft(projections, N_reflections, d_max); - // refine ds, write 1 to fits[i] if reflections[i] fits ds - ds = refine_ds(projections, N_reflections, ds, LevelFit, fits); + /* 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); + /* n - number of reflections fitting ds */ + n = check_refl_fitting_ds(projections, N_reflections, + ds, LevelFit); - // normal/ds - possible direct vector + /* normal/ds - possible direct vector */ gsl_vector_scale(normal, 1/ds); - if (n > N_reflections/3 && n > 6) { + if ( n > N_reflections/3 && n > 6 ) { tvectors[N_tvectors] = tvector_new(N_reflections); - gsl_vector_memcpy(tvectors[N_tvectors].t, normal); + gsl_vector_memcpy(tvectors[N_tvectors].t, + normal); memcpy(tvectors[N_tvectors].fits, fits, N_reflections * sizeof(int)); @@ -874,37 +901,39 @@ static int index_refls(gsl_vector **reflections, int N_reflections, double d_max } } - if ((i != 0 && i % 10000 == 0) || i == i_max - 1) { + if ( (i != 0 && i % 10000 == 0) || i == i_max - 1 ) { - // sort tvectors by length + /* 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); + /* 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) { + if ( c->n > 4 * n_max / 5 ) { break; } } } free(fits); - for (i = 0; i < N_tvectors; i++) { + for ( i = 0; i < N_tvectors; i++ ) { tvector_free(tvectors[i]); } free(tvectors); - for (i = 0; i < N_triplets; i++) { + for ( i = 0; i < N_triplets; i++ ) { free(triplets[i]); } free(triplets); - if (c->n) return 1; + if ( c->n ) return 1; return 0; } @@ -936,30 +965,30 @@ double cell_get_volume(UnitCell *cell) } int run_asdf(struct image *image, IndexingPrivate *ipriv) { - int i; + int i; - double LevelFit = 1./1000; - double IndexFit = 1./500; - double d_max = 220.; // thrice the maximum expected axis length - double volume_min = 100.; - double volume_max = 1000000.; - - int i_max = 10000; // maximum number of triplets - - struct asdf_private *dp = (struct asdf_private *)ipriv; - - if ( dp->indm & INDEXING_CHECK_CELL_AXES ) { + double LevelFit = 1./1000; + double IndexFit = 1./500; + double d_max = 220.; // thrice the maximum expected axis length + double volume_min = 100.; + double volume_max = 1000000.; + + int i_max = 10000; // maximum number of triplets + + struct asdf_private *dp = (struct asdf_private *)ipriv; + + if ( dp->indm & INDEXING_CHECK_CELL_AXES ) { double volume = cell_get_volume(dp->template); volume_min = volume * 0.95; volume_max = volume * 1.05; } - int N_reflections = image_feature_count(image->features); - gsl_vector *reflections[N_reflections]; - - for (i = 0; i < N_reflections; i++) { + int N_reflections = image_feature_count(image->features); + gsl_vector *reflections[N_reflections]; + + for ( i = 0; i < N_reflections; i++ ) { struct imagefeature *f; - + f = image_get_feature(image->features, i); if ( f == NULL ) continue; @@ -969,38 +998,37 @@ int run_asdf(struct image *image, IndexingPrivate *ipriv) { gsl_vector_set(reflections[i], 2, f->rz/1e10); } - struct asdf_cell c = asdf_cell_new(N_reflections); - - if (N_reflections == 0) return 0; - - i = index_refls(reflections, N_reflections, d_max, volume_min, volume_max, - LevelFit, IndexFit, i_max, &c); - //~ - - for (i = 0; i < N_reflections; i ++) { + struct asdf_cell c = asdf_cell_new(N_reflections); + + if ( N_reflections == 0 ) return 0; + + i = index_refls(reflections, N_reflections, d_max, volume_min, volume_max, + LevelFit, IndexFit, i_max, &c); + + for ( i = 0; i < N_reflections; i++ ) { gsl_vector_free(reflections[i]); } - - if (i) { + + if ( i ) { UnitCell *uc; uc = cell_new(); cell_set_cartesian(uc, gsl_vector_get(c.axes[0], 0) * 1e-10, - gsl_vector_get(c.axes[0], 1) * 1e-10, - gsl_vector_get(c.axes[0], 2) * 1e-10, - gsl_vector_get(c.axes[1], 0) * 1e-10, - gsl_vector_get(c.axes[1], 1) * 1e-10, - gsl_vector_get(c.axes[1], 2) * 1e-10, - gsl_vector_get(c.axes[2], 0) * 1e-10, - gsl_vector_get(c.axes[2], 1) * 1e-10, - gsl_vector_get(c.axes[2], 2) * 1e-10); - + gsl_vector_get(c.axes[0], 1) * 1e-10, + gsl_vector_get(c.axes[0], 2) * 1e-10, + gsl_vector_get(c.axes[1], 0) * 1e-10, + gsl_vector_get(c.axes[1], 1) * 1e-10, + gsl_vector_get(c.axes[1], 2) * 1e-10, + gsl_vector_get(c.axes[2], 0) * 1e-10, + gsl_vector_get(c.axes[2], 1) * 1e-10, + gsl_vector_get(c.axes[2], 2) * 1e-10); + if ( check_cell(dp, image, uc) ) { cell_free(uc); return 1; } - cell_free(uc); + cell_free(uc); } return 0; -- cgit v1.2.3 From 3605e92c4461d4742df22e98e9d1a4bc8e7cbe33 Mon Sep 17 00:00:00 2001 From: Alexandra Tolstikova Date: Wed, 11 Feb 2015 14:00:58 +0100 Subject: Use one fftw_plan for all fft calculations --- libcrystfel/src/asdf.c | 27 +++++++++++++++++---------- 1 file changed, 17 insertions(+), 10 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/asdf.c b/libcrystfel/src/asdf.c index 6ad29645..41973c9a 100644 --- a/libcrystfel/src/asdf.c +++ b/libcrystfel/src/asdf.c @@ -161,7 +161,8 @@ static int cross_product(gsl_vector *vec1, gsl_vector *vec2, } /* Returns triple product of three gsl_vectors */ -static double calc_volume(gsl_vector *vec1, gsl_vector *vec2, gsl_vector *vec3) { +static double calc_volume(gsl_vector *vec1, gsl_vector *vec2, gsl_vector *vec3) +{ double volume; gsl_vector *cross = gsl_vector_alloc(3); @@ -289,7 +290,8 @@ static int calc_normal(gsl_vector *c1, gsl_vector *c2, gsl_vector *c3, return 1; } -static float find_ds_fft(double *projections, int N_projections, double d_max) { +static float find_ds_fft(double *projections, int N_projections, double d_max, + fftw_plan p) { int n = N_projections; double projections_sorted[n]; @@ -300,7 +302,7 @@ static float find_ds_fft(double *projections, int N_projections, double d_max) { int N = 1024; // number of points in fft calculation - double *in = malloc(sizeof(double) * N); + double in[N]; fftw_complex *out = fftw_malloc(sizeof (fftw_complex) * N); for ( i = 0; i < N; i++ ) { @@ -311,10 +313,12 @@ static float find_ds_fft(double *projections, int N_projections, double d_max) { in[(int)((projections_sorted[i] - projections_sorted[0]) / (projections_sorted[n - 1] - projections_sorted[0]) * N)] ++; } - - fftw_plan p = fftw_plan_dft_r2c_1d(N, in, out, FFTW_MEASURE); - fftw_execute(p); - fftw_destroy_plan(p); + + if ( !p ) { + p = fftw_plan_dft_r2c_1d(N, in, out, FFTW_MEASURE); + } + + fftw_execute_dft_r2c(p, in, out); int i_max = (int)(d_max * (projections_sorted[n - 1] - projections_sorted[0])); @@ -332,7 +336,6 @@ static float find_ds_fft(double *projections, int N_projections, double d_max) { double ds = (projections_sorted[n - 1] - projections_sorted[0]) / d; - free(in); fftw_free(out); return ds; } @@ -855,7 +858,9 @@ static int index_refls(gsl_vector **reflections, int N_reflections, int N_tvectors = 0; int n_max = 0; // maximum number of reflections fitting one of tvectors - + + fftw_plan p; + for ( i = 0; i < i_max; i++ ) { if ( calc_normal(reflections[triplets[i][0]], reflections[triplets[i][1]], @@ -870,7 +875,7 @@ static int index_refls(gsl_vector **reflections, int N_reflections, } /* Find ds - period in 1d lattice of projections */ - ds = find_ds_fft(projections, N_reflections, d_max); + ds = find_ds_fft(projections, N_reflections, d_max, p); /* Refine ds, write 1 to fits[i] if reflections[i] * fits ds */ @@ -933,6 +938,8 @@ static int index_refls(gsl_vector **reflections, int N_reflections, } free(triplets); + fftw_destroy_plan(p); + if ( c->n ) return 1; return 0; -- cgit v1.2.3 From a095923e5ee59201c43e3722903bae8745e40bf3 Mon Sep 17 00:00:00 2001 From: Alexandra Tolstikova Date: Wed, 11 Feb 2015 15:42:24 +0100 Subject: Move creation of a cell from tvectors into a separate function --- libcrystfel/src/asdf.c | 309 +++++++++++++++++++++++++++++-------------------- 1 file changed, 181 insertions(+), 128 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/asdf.c b/libcrystfel/src/asdf.c index 41973c9a..28afe82d 100644 --- a/libcrystfel/src/asdf.c +++ b/libcrystfel/src/asdf.c @@ -31,6 +31,7 @@ struct tvector { int *fits; }; + struct asdf_cell { gsl_vector *axes[3]; gsl_vector *reciprocal[3]; @@ -46,7 +47,9 @@ struct asdf_cell { int n_max; // maximum number of reflections fitting to some t-vector }; -struct tvector tvector_new(int n) { + +struct tvector tvector_new(int n) +{ struct tvector t; t.t = gsl_vector_alloc(3); @@ -56,14 +59,18 @@ struct tvector tvector_new(int n) { return t; } -static int tvector_free(struct tvector t) { + +static int tvector_free(struct tvector t) +{ gsl_vector_free(t.t); free(t.fits); return 1; } -static int tvector_memcpy(struct tvector *dest, struct tvector *src, int n) { + +static int tvector_memcpy(struct tvector *dest, struct tvector *src, int n) +{ gsl_vector_memcpy(dest->t, src->t); dest->n = src->n; memcpy(dest->fits, src->fits, sizeof(int) * n); @@ -71,7 +78,9 @@ static int tvector_memcpy(struct tvector *dest, struct tvector *src, int n) { return 1; } -static int asdf_cell_free(struct asdf_cell c) { + +static int asdf_cell_free(struct asdf_cell c) +{ int i; for ( i = 0; i < 3; i++ ) { gsl_vector_free(c.axes[i]); @@ -87,8 +96,9 @@ static int asdf_cell_free(struct asdf_cell c) { return 1; } -static struct asdf_cell asdf_cell_new(int n) { - + +static struct asdf_cell asdf_cell_new(int n) +{ struct asdf_cell c; int i; @@ -113,7 +123,9 @@ static struct asdf_cell asdf_cell_new(int n) { return c; } -static int asdf_cell_memcpy(struct asdf_cell *dest, struct asdf_cell *src) { + +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]); @@ -139,9 +151,11 @@ static int asdf_cell_memcpy(struct asdf_cell *dest, struct asdf_cell *src) { return 1; } + /* result = vec1 cross vec2 */ static int cross_product(gsl_vector *vec1, gsl_vector *vec2, - gsl_vector **result) { + gsl_vector **result) +{ double c1[3], c2[3], p[3]; int i; for ( i = 0; i < 3; i++ ) { @@ -160,6 +174,7 @@ static int cross_product(gsl_vector *vec1, gsl_vector *vec2, return 1; } + /* Returns triple product of three gsl_vectors */ static double calc_volume(gsl_vector *vec1, gsl_vector *vec2, gsl_vector *vec3) { @@ -173,7 +188,9 @@ static double calc_volume(gsl_vector *vec1, gsl_vector *vec2, gsl_vector *vec3) return volume; } -static int calc_reciprocal(gsl_vector **direct, gsl_vector **reciprocal) { + +static int calc_reciprocal(gsl_vector **direct, gsl_vector **reciprocal) +{ double volume; cross_product(direct[1], direct[2], &reciprocal[0]); @@ -189,6 +206,7 @@ static int calc_reciprocal(gsl_vector **direct, gsl_vector **reciprocal) { return 1; } + static int check_cell(struct asdf_private *dp, struct image *image, UnitCell *cell) { @@ -230,6 +248,7 @@ static int check_cell(struct asdf_private *dp, struct image *image, return 1; } + static int compare_doubles (const void *a, const void *b) { const double *da = (const double *) a; @@ -238,6 +257,7 @@ static int compare_doubles (const void *a, const void *b) return (*da > *db) - (*da < *db); } + /* Compares tvectors by length */ static int compare_tvectors (const void *a, const void *b) { @@ -252,11 +272,12 @@ static int compare_tvectors (const void *a, const void *b) //~ 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 *normal) +{ gsl_vector *c12 = gsl_vector_alloc(3); gsl_vector *c23 = gsl_vector_alloc(3); gsl_vector *c31 = gsl_vector_alloc(3); @@ -290,9 +311,10 @@ static int calc_normal(gsl_vector *c1, gsl_vector *c2, gsl_vector *c3, return 1; } + static float find_ds_fft(double *projections, int N_projections, double d_max, - fftw_plan p) { - + fftw_plan p) +{ int n = N_projections; double projections_sorted[n]; memcpy(projections_sorted, projections, sizeof(double) * n); @@ -340,13 +362,14 @@ static float find_ds_fft(double *projections, int N_projections, double d_max, 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) { - + double ds, double LevelFit) +{ if ( ds == 0 ) return 0; int i; @@ -362,10 +385,11 @@ static int check_refl_fitting_ds(double *projections, int N_projections, 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 LevelFit, int *fits) +{ double fit_refls[N_projections]; double indices[N_projections]; @@ -404,10 +428,11 @@ static float refine_ds(double *projections, int N_projections, double ds, return ds_new; } + static int check_refl_fitting_cell(struct asdf_cell *c, gsl_vector **reflections, - int N_reflections, double IndexFit) { - + int N_reflections, double IndexFit) +{ double dist[3]; calc_reciprocal(c->axes, c->reciprocal); @@ -443,7 +468,9 @@ static int check_refl_fitting_cell(struct asdf_cell *c, return 1; } -static void print_asdf_cell(struct asdf_cell cc) { + +static void print_asdf_cell(struct asdf_cell cc) +{ double a, b, c, alpha, beta, gamma, ab, bc, ca; a = gsl_blas_dnrm2(cc.axes[0]); @@ -472,16 +499,17 @@ static void print_asdf_cell(struct asdf_cell cc) { } + /* 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); +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 *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); @@ -489,32 +517,33 @@ static int refine_asdf_cell(struct asdf_cell c, gsl_vector **reflections, int i, j; int n = 0; - for ( i = 0; i < N_reflections; i++ ) if ( c.reflections[i] == 1 ) + 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_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); + 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_set(c->reciprocal[j], i, gsl_vector_get(res, j)); } } - calc_reciprocal(c.reciprocal, c.axes); + calc_reciprocal(c->reciprocal, c->axes); double a[3]; for ( i = 0; i < 3; i++ ) { - a[i] = gsl_blas_dnrm2(c.axes[i]); + a[i] = gsl_blas_dnrm2(c->axes[i]); } gsl_multifit_linear_free(work); @@ -534,7 +563,9 @@ static int refine_asdf_cell(struct asdf_cell c, gsl_vector **reflections, return 1; } -static int reduce_asdf_cell(struct asdf_cell *cl) { + +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); @@ -649,6 +680,7 @@ static int reduce_asdf_cell(struct asdf_cell *cl) { return 1; } + static int check_cell_angles(gsl_vector *va, gsl_vector *vb, gsl_vector *vc, double max_cos) { @@ -674,123 +706,134 @@ static int check_cell_angles(gsl_vector *va, gsl_vector *vb, gsl_vector *vc, return 1; } + /* Returns min(t1.n, t2.n, t3.n) */ -static int find_acl(struct tvector t1, struct tvector t2, struct tvector t3) { +static int find_acl(struct tvector t1, struct tvector t2, struct tvector t3) +{ int i = t1.n, j = t2.n, k = t3.n; if ( i <= j && i <= k ) return i; if ( j <= i && j <= k ) return j; if ( k <= i && k <= j ) return k; } + +static int 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 && 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, n, cell_correct; - - double volume; - + 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); - + /* 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; + int i_min = sl < 2 * N_tvectors ? 0 : sl - 2 * N_tvectors; + int i_max = sl < N_tvectors ? sl : 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) ) - { + for ( i = i_min; i < i_max; i++) if (tvectors[i].n > acl ) { - volume = calc_volume(tvectors[i].t, - tvectors[j].t, - tvectors[k].t); - - if ( fabs(volume) > volume_min && - fabs(volume) < volume_max ) - { - - gsl_vector_memcpy(c.axes[0], tvectors[i].t); - gsl_vector_memcpy(c.axes[1], tvectors[j].t); - gsl_vector_memcpy(c.axes[2], tvectors[k].t); - - c.volume = volume; - check_refl_fitting_cell(&c, reflections, - N_reflections, - IndexFit); - - if ( c.n < 6 ) break; - - reduce_asdf_cell(&c); + 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; - /* if one of the cell angles > 135 or < 45 - * do not continue */ - if ( !check_cell_angles(c.axes[0], c.axes[1], - c.axes[2], 0.71) ) break; - - /* index reflections with new cell axes */ - check_refl_fitting_cell(&c, reflections, - N_reflections, - IndexFit); - - acl = find_acl(tvectors[i], - tvectors[j], - tvectors[k]); - - c.acl = acl; - c.n_max = n_max; + for ( j = j_min; j < j_max; j++ ) + if ( tvectors[j].n > acl ) { + + k = sl - i - j - 1; - /* refine cell until the number of fitting - * reflections stops increasing */ - n = 0; - cell_correct = 1; - while ( c.n - n && cell_correct ) { - n = c.n; - cell_correct = refine_asdf_cell(c, - reflections, - N_reflections, - IndexFit); + 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; - check_refl_fitting_cell(&c, reflections, - N_reflections, - IndexFit); - } - - if ( cell_correct ) { 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; + 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; - if ( tvectors[i].n <= acl ) break; - } - if ( acl > n_max ) break; + if ( acl > n_max ) break; } asdf_cell_free(c); @@ -799,7 +842,9 @@ static int find_cell(struct tvector *tvectors, int N_tvectors, double IndexFit, return 0; } -static void shuffle_triplets(int **triplets, int n) { + +static void shuffle_triplets(int **triplets, int n) +{ int i, j; int t[3]; for ( i = 0; i < n - 1; i++ ) { @@ -810,16 +855,20 @@ static void shuffle_triplets(int **triplets, int n) { } } -static double angle_between_gsl(gsl_vector *a, gsl_vector *b) { + +static double angle_between_gsl(gsl_vector *a, gsl_vector *b) +{ double ab; gsl_blas_ddot(a, b, &ab); return acos(ab/gsl_blas_dnrm2(a)/gsl_blas_dnrm2(b)) * 180 / M_PI; } + static int index_refls(gsl_vector **reflections, int N_reflections, double d_max, double volume_min, double volume_max, double LevelFit, double IndexFit, int i_max, - struct asdf_cell *c) { + struct asdf_cell *c) +{ int i, j, k, l, n; @@ -945,6 +994,7 @@ static int index_refls(gsl_vector **reflections, int N_reflections, return 0; } + double cell_get_volume(UnitCell *cell) { double asx, asy, asz; @@ -971,8 +1021,10 @@ double cell_get_volume(UnitCell *cell) return 1/volume; } -int run_asdf(struct image *image, IndexingPrivate *ipriv) { - int i; + +int run_asdf(struct image *image, IndexingPrivate *ipriv) +{ + int i, j; double LevelFit = 1./1000; double IndexFit = 1./500; @@ -1009,14 +1061,14 @@ int run_asdf(struct image *image, IndexingPrivate *ipriv) { if ( N_reflections == 0 ) return 0; - i = index_refls(reflections, N_reflections, d_max, volume_min, volume_max, - LevelFit, IndexFit, i_max, &c); + j = index_refls(reflections, N_reflections, d_max, volume_min, + volume_max, LevelFit, IndexFit, i_max, &c); for ( i = 0; i < N_reflections; i++ ) { gsl_vector_free(reflections[i]); } - if ( i ) { + if ( j ) { UnitCell *uc; uc = cell_new(); @@ -1041,6 +1093,7 @@ int run_asdf(struct image *image, IndexingPrivate *ipriv) { return 0; } + IndexingPrivate *asdf_prepare(IndexingMethod *indm, UnitCell *cell, struct detector *det, float *ltl) { -- cgit v1.2.3 From 21b67e8dc422ee5f437bf2e31c5616982fc41c98 Mon Sep 17 00:00:00 2001 From: Alexandra Tolstikova Date: Wed, 11 Feb 2015 16:18:18 +0100 Subject: Move cell_get_volume from asdf.c to cell_utils.c --- libcrystfel/src/asdf.c | 27 --------------------------- libcrystfel/src/cell-utils.c | 28 ++++++++++++++++++++++++++++ libcrystfel/src/cell-utils.h | 2 ++ 3 files changed, 30 insertions(+), 27 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/asdf.c b/libcrystfel/src/asdf.c index 28afe82d..0b890ef8 100644 --- a/libcrystfel/src/asdf.c +++ b/libcrystfel/src/asdf.c @@ -995,33 +995,6 @@ static int index_refls(gsl_vector **reflections, int N_reflections, } -double cell_get_volume(UnitCell *cell) -{ - double asx, asy, asz; - double bsx, bsy, bsz; - double csx, csy, csz; - struct rvec aCb; - double volume; - - if ( cell_get_reciprocal(cell, &asx, &asy, &asz, - &bsx, &bsy, &bsz, - &csx, &csy, &csz) ) { - ERROR("Couldn't get cell cartesian.\n"); - return 0; - } - - /* "a" cross "b" */ - aCb.u = asy*bsz - asz*bsy; - aCb.v = - (asx*bsz - asz*bsx); - aCb.w = asx*bsy - asy*bsx; - - /* "a cross b" dot "c" */ - volume = (aCb.u*csx + aCb.v*csy + aCb.w*csz)/1e30; - - return 1/volume; -} - - int run_asdf(struct image *image, IndexingPrivate *ipriv) { int i, j; diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c index 95e7353e..5fc885e5 100644 --- a/libcrystfel/src/cell-utils.c +++ b/libcrystfel/src/cell-utils.c @@ -1567,3 +1567,31 @@ int forbidden_reflection(UnitCell *cell, return 0; } + + +/* Returns cell volume in A^3 */ +double cell_get_volume(UnitCell *cell) +{ + double asx, asy, asz; + double bsx, bsy, bsz; + double csx, csy, csz; + struct rvec aCb; + double rec_volume; + + if ( cell_get_reciprocal(cell, &asx, &asy, &asz, + &bsx, &bsy, &bsz, + &csx, &csy, &csz) ) { + ERROR("Couldn't get reciprocal cell.\n"); + return 0; + } + + /* "a" cross "b" */ + aCb.u = asy*bsz - asz*bsy; + aCb.v = - (asx*bsz - asz*bsx); + aCb.w = asx*bsy - asy*bsx; + + /* "a cross b" dot "c" */ + rec_volume = (aCb.u*csx + aCb.v*csy + aCb.w*csz)/1e30; + + return 1/rec_volume; +} diff --git a/libcrystfel/src/cell-utils.h b/libcrystfel/src/cell-utils.h index 0b900096..98bc667c 100644 --- a/libcrystfel/src/cell-utils.h +++ b/libcrystfel/src/cell-utils.h @@ -74,6 +74,8 @@ extern LatticeType lattice_from_str(const char *s); extern int forbidden_reflection(UnitCell *cell, signed int h, signed int k, signed int l); +extern double cell_get_volume(UnitCell *cell); + #ifdef __cplusplus } #endif -- cgit v1.2.3 From d6f855a933aa900f47be27649d42895393420761 Mon Sep 17 00:00:00 2001 From: Alexandra Tolstikova Date: Wed, 11 Feb 2015 17:11:21 +0100 Subject: Get cell volume range from tolerances --- libcrystfel/src/asdf.c | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/asdf.c b/libcrystfel/src/asdf.c index 0b890ef8..1d3d2c89 100644 --- a/libcrystfel/src/asdf.c +++ b/libcrystfel/src/asdf.c @@ -1011,8 +1011,10 @@ int run_asdf(struct image *image, IndexingPrivate *ipriv) if ( dp->indm & INDEXING_CHECK_CELL_AXES ) { double volume = cell_get_volume(dp->template); - volume_min = volume * 0.95; - volume_max = volume * 1.05; + double vtol = (dp->ltl[0] + dp->ltl[1] + dp->ltl[2]) / 100 + + dp->ltl[3] / 180 * M_PI; + volume_min = volume * (1 - vtol); + volume_max = volume * (1 + vtol); } int N_reflections = image_feature_count(image->features); @@ -1060,7 +1062,7 @@ int run_asdf(struct image *image, IndexingPrivate *ipriv) return 1; } - cell_free(uc); + cell_free(uc); } return 0; -- cgit v1.2.3 From d2155d4115786c3520af6d9b76a85d3921fa034d Mon Sep 17 00:00:00 2001 From: Alexandra Tolstikova Date: Wed, 25 Mar 2015 10:19:07 +0100 Subject: Fix malloc and memcpy of asdf_cell.indices (int -> double) --- libcrystfel/src/asdf.c | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/asdf.c b/libcrystfel/src/asdf.c index 1d3d2c89..53214a12 100644 --- a/libcrystfel/src/asdf.c +++ b/libcrystfel/src/asdf.c @@ -110,9 +110,9 @@ static struct asdf_cell asdf_cell_new(int n) c.N_refls = n; c.reflections = malloc(sizeof(int) * n); - c.indices = malloc(sizeof(int *) * n); + c.indices = malloc(sizeof(double *) * n); for ( i = 0; i < n; i++ ) { - c.indices[i] = malloc(sizeof(int) * 3); + c.indices[i] = malloc(sizeof(double) * 3); } c.n = 0; @@ -139,11 +139,9 @@ static int asdf_cell_memcpy(struct asdf_cell *dest, struct asdf_cell *src) dest->n = src->n; memcpy(dest->reflections, src->reflections, sizeof(int) * n); - - memcpy(dest->indices, src->indices, sizeof(int *) * n); - + for (i = 0; i < n; i++ ) { - memcpy(dest->indices[i], src->indices[i], sizeof(int) * 3); + memcpy(dest->indices[i], src->indices[i], sizeof(double) * 3); } dest->acl = src->acl; @@ -1012,7 +1010,7 @@ int run_asdf(struct image *image, IndexingPrivate *ipriv) if ( dp->indm & INDEXING_CHECK_CELL_AXES ) { double volume = cell_get_volume(dp->template); double vtol = (dp->ltl[0] + dp->ltl[1] + dp->ltl[2]) / 100 + - dp->ltl[3] / 180 * M_PI; + dp->ltl[3] / 180 * M_PI; volume_min = volume * (1 - vtol); volume_max = volume * (1 + vtol); } -- cgit v1.2.3 From 218a6373e90f7e1d5b29afa2c582a90d3b0214c6 Mon Sep 17 00:00:00 2001 From: Alexandra Tolstikova Date: Mon, 20 Apr 2015 21:03:47 +0200 Subject: Create fftw_plan and in and out arrays for fft in asdf_prepare(), which also partially solves memory leak problem --- libcrystfel/src/asdf.c | 72 +++++++++++++++++++++++++++++++++++--------------- 1 file changed, 51 insertions(+), 21 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/asdf.c b/libcrystfel/src/asdf.c index 53214a12..2f3ddd04 100644 --- a/libcrystfel/src/asdf.c +++ b/libcrystfel/src/asdf.c @@ -17,10 +17,19 @@ #include "asdf.h" +struct fftw_vars { + int N; + fftw_plan p; + double *in; + fftw_complex *out; +}; + + struct asdf_private { IndexingMethod indm; float *ltl; UnitCell *template; + struct fftw_vars fftw; }; @@ -32,6 +41,30 @@ struct tvector { }; +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 int 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]; @@ -311,7 +344,7 @@ static int calc_normal(gsl_vector *c1, gsl_vector *c2, gsl_vector *c3, static float find_ds_fft(double *projections, int N_projections, double d_max, - fftw_plan p) + struct fftw_vars fftw) { int n = N_projections; double projections_sorted[n]; @@ -320,10 +353,10 @@ static float find_ds_fft(double *projections, int N_projections, double d_max, int i; - int N = 1024; // number of points in fft calculation - - double in[N]; - fftw_complex *out = fftw_malloc(sizeof (fftw_complex) * N); + 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; @@ -334,10 +367,6 @@ static float find_ds_fft(double *projections, int N_projections, double d_max, (projections_sorted[n - 1] - projections_sorted[0]) * N)] ++; } - if ( !p ) { - p = fftw_plan_dft_r2c_1d(N, in, out, FFTW_MEASURE); - } - fftw_execute_dft_r2c(p, in, out); int i_max = (int)(d_max * (projections_sorted[n - 1] - @@ -356,7 +385,6 @@ static float find_ds_fft(double *projections, int N_projections, double d_max, double ds = (projections_sorted[n - 1] - projections_sorted[0]) / d; - fftw_free(out); return ds; } @@ -865,7 +893,7 @@ static double angle_between_gsl(gsl_vector *a, gsl_vector *b) static int index_refls(gsl_vector **reflections, int N_reflections, double d_max, double volume_min, double volume_max, double LevelFit, double IndexFit, int i_max, - struct asdf_cell *c) + struct asdf_cell *c, struct fftw_vars fftw) { int i, j, k, l, n; @@ -906,8 +934,6 @@ static int index_refls(gsl_vector **reflections, int N_reflections, int n_max = 0; // maximum number of reflections fitting one of tvectors - fftw_plan p; - for ( i = 0; i < i_max; i++ ) { if ( calc_normal(reflections[triplets[i][0]], reflections[triplets[i][1]], @@ -922,7 +948,8 @@ static int index_refls(gsl_vector **reflections, int N_reflections, } /* Find ds - period in 1d lattice of projections */ - ds = find_ds_fft(projections, N_reflections, d_max, p); + ds = find_ds_fft(projections, N_reflections, d_max, + fftw); /* Refine ds, write 1 to fits[i] if reflections[i] * fits ds */ @@ -952,7 +979,7 @@ static int index_refls(gsl_vector **reflections, int N_reflections, if (n > n_max) n_max = n; } } - + if ( (i != 0 && i % 10000 == 0) || i == i_max - 1 ) { /* Sort tvectors by length */ @@ -984,9 +1011,7 @@ static int index_refls(gsl_vector **reflections, int N_reflections, free(triplets[i]); } free(triplets); - - fftw_destroy_plan(p); - + if ( c->n ) return 1; return 0; @@ -1035,7 +1060,7 @@ int run_asdf(struct image *image, IndexingPrivate *ipriv) if ( N_reflections == 0 ) return 0; j = index_refls(reflections, N_reflections, d_max, volume_min, - volume_max, LevelFit, IndexFit, i_max, &c); + volume_max, LevelFit, IndexFit, i_max, &c, dp->fftw); for ( i = 0; i < N_reflections; i++ ) { gsl_vector_free(reflections[i]); @@ -1055,7 +1080,8 @@ int run_asdf(struct image *image, IndexingPrivate *ipriv) gsl_vector_get(c.axes[2], 1) * 1e-10, gsl_vector_get(c.axes[2], 2) * 1e-10); - if ( check_cell(dp, image, uc) ) { + if ( check_cell(dp, image, uc) ) { + asdf_cell_free(c); cell_free(uc); return 1; } @@ -1063,6 +1089,7 @@ int run_asdf(struct image *image, IndexingPrivate *ipriv) cell_free(uc); } + asdf_cell_free(c); return 0; } @@ -1093,7 +1120,9 @@ IndexingPrivate *asdf_prepare(IndexingMethod *indm, UnitCell *cell, dp->ltl = ltl; dp->template = cell; dp->indm = *indm; - + + dp->fftw = fftw_vars_new(); + return (IndexingPrivate *)dp; } @@ -1102,5 +1131,6 @@ void asdf_cleanup(IndexingPrivate *pp) { struct asdf_private *p; p = (struct asdf_private *)pp; + fftw_vars_free(p->fftw); free(p); } -- cgit v1.2.3 From a2367c209492e8b91b21f435ea1b3cb99fc69cc7 Mon Sep 17 00:00:00 2001 From: Alexandra Tolstikova Date: Wed, 22 Apr 2015 16:25:33 +0200 Subject: Free gsl_vector normal in index_refls() --- libcrystfel/src/asdf.c | 2 ++ 1 file changed, 2 insertions(+) (limited to 'libcrystfel') diff --git a/libcrystfel/src/asdf.c b/libcrystfel/src/asdf.c index 2f3ddd04..427b6e52 100644 --- a/libcrystfel/src/asdf.c +++ b/libcrystfel/src/asdf.c @@ -1012,6 +1012,8 @@ static int index_refls(gsl_vector **reflections, int N_reflections, } free(triplets); + gsl_vector_free(normal); + if ( c->n ) return 1; return 0; -- cgit v1.2.3 From 067842bd0480ad7f4e5680158e80c40d23d08041 Mon Sep 17 00:00:00 2001 From: Alexandra Tolstikova Date: Tue, 16 Jun 2015 15:16:42 +0200 Subject: Check malloc return values --- libcrystfel/src/asdf.c | 129 +++++++++++++++++++++++++++++++------------------ 1 file changed, 81 insertions(+), 48 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/asdf.c b/libcrystfel/src/asdf.c index 427b6e52..fb2f8279 100644 --- a/libcrystfel/src/asdf.c +++ b/libcrystfel/src/asdf.c @@ -112,46 +112,52 @@ static int tvector_memcpy(struct tvector *dest, struct tvector *src, int n) } -static int asdf_cell_free(struct asdf_cell c) +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]); + 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->reflections); + for ( i = 0; i < c->N_refls; i++ ) { + free(c->indices[i]); } - free(c.indices); + free(c->indices); + free(c); return 1; } -static struct asdf_cell asdf_cell_new(int n) +static struct asdf_cell *asdf_cell_new(int n) { - struct asdf_cell c; + 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->axes[i] = gsl_vector_alloc(3); + c->reciprocal[i] = gsl_vector_alloc(3); } - c.N_refls = n; - c.reflections = malloc(sizeof(int) * n); + 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; - c.indices = malloc(sizeof(double *) * n); for ( i = 0; i < n; i++ ) { - c.indices[i] = malloc(sizeof(double) * 3); + c->indices[i] = malloc(sizeof(double) * 3); + if (c->indices[i] == NULL) return NULL; } - c.n = 0; + c->n = 0; - c.acl = 0; - c.n_max = 0; + c->acl = 0; + c->n_max = 0; return c; } @@ -495,17 +501,17 @@ static int check_refl_fitting_cell(struct asdf_cell *c, } -static void print_asdf_cell(struct asdf_cell cc) +static void print_asdf_cell(struct asdf_cell *cc) { double a, b, c, alpha, beta, gamma, ab, bc, ca; - a = gsl_blas_dnrm2(cc.axes[0]); - b = gsl_blas_dnrm2(cc.axes[1]); - c = gsl_blas_dnrm2(cc.axes[2]); + a = gsl_blas_dnrm2(cc->axes[0]); + b = gsl_blas_dnrm2(cc->axes[1]); + c = gsl_blas_dnrm2(cc->axes[2]); - gsl_blas_ddot(cc.axes[0], cc.axes[1], &ab); - gsl_blas_ddot(cc.axes[1], cc.axes[2], &bc); - gsl_blas_ddot(cc.axes[0], cc.axes[2], &ca); + gsl_blas_ddot(cc->axes[0], cc->axes[1], &ab); + gsl_blas_ddot(cc->axes[1], cc->axes[2], &bc); + gsl_blas_ddot(cc->axes[0], cc->axes[2], &ca); alpha = acos(bc/b/c)/M_PI*180; beta = acos(ca/a/c)/M_PI*180; @@ -521,7 +527,7 @@ static void print_asdf_cell(struct asdf_cell cc) printf("%.2f %.2f %.2f %.2f %.2f %.2f %.0f %d \n", a, b, c, alpha, beta, gamma, - cc.volume, cc.n); + cc->volume, cc->n); } @@ -534,8 +540,8 @@ static int refine_asdf_cell(struct asdf_cell *c, gsl_vector **reflections, 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_alloc(c->n), + gsl_vector_alloc(c->n)}; gsl_vector *res = gsl_vector_alloc(3); gsl_matrix *cov = gsl_matrix_alloc (3, 3); @@ -725,7 +731,7 @@ static int check_cell_angles(gsl_vector *va, gsl_vector *vb, gsl_vector *vc, cosg = ab/a/b; if ( fabs(cosa) > max_cos || fabs(cosb) > max_cos || - fabs(cosg) > max_cos ) { + fabs(cosg) > max_cos ) { return 0; } @@ -797,7 +803,11 @@ static int find_cell(struct tvector *tvectors, int N_tvectors, double IndexFit, * considered */ int acl = N_reflections < 18 ? 6 : N_reflections/3; - struct asdf_cell c = asdf_cell_new(N_reflections); + 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; @@ -827,7 +837,7 @@ static int find_cell(struct tvector *tvectors, int N_tvectors, double IndexFit, if ( !create_cell(tvectors[i], tvectors[j], tvectors[k], - &c, IndexFit, + c, IndexFit, volume_min, volume_max, reflections, @@ -839,15 +849,15 @@ static int find_cell(struct tvector *tvectors, int N_tvectors, double IndexFit, acl = find_acl(tvectors[i], tvectors[j], tvectors[k]); - c.acl = acl; - c.n_max = n_max; + c->acl = acl; + c->n_max = n_max; - reduce_asdf_cell(&c); + 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); + if ( result->n < c->n ) { + asdf_cell_memcpy(result, c); } acl++; @@ -903,12 +913,22 @@ static int index_refls(gsl_vector **reflections, int N_reflections, (N_reflections - 2) / 6; int **triplets = malloc(N_triplets * sizeof(int *)); + if (triplets == NULL) { + ERROR("Failed to allocate triplets in index_refls!\n"); + return 0; + } + l = 0; for ( i = 0; i < N_reflections; i++ ) { for ( j = i + 1; j < N_reflections; j++ ) { for ( k = j + 1; k < N_reflections; k++ ) { triplets[l] = malloc(3 * sizeof(int)); - + if (triplets[l] == NULL) { + ERROR("Failed to allocate triplets " + " in index_refls!\n"); + return 0; + } + triplets[l][0] = i; triplets[l][1] = j; triplets[l][2] = k; @@ -926,10 +946,19 @@ static int index_refls(gsl_vector **reflections, int 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; + } if ( i_max > N_triplets ) i_max = N_triplets; struct tvector *tvectors = malloc(i_max * 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 @@ -1057,12 +1086,16 @@ int run_asdf(struct image *image, IndexingPrivate *ipriv) gsl_vector_set(reflections[i], 2, f->rz/1e10); } - struct asdf_cell c = asdf_cell_new(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, i_max, &c, dp->fftw); + volume_max, LevelFit, IndexFit, i_max, c, dp->fftw); for ( i = 0; i < N_reflections; i++ ) { gsl_vector_free(reflections[i]); @@ -1072,15 +1105,15 @@ int run_asdf(struct image *image, IndexingPrivate *ipriv) UnitCell *uc; uc = cell_new(); - cell_set_cartesian(uc, gsl_vector_get(c.axes[0], 0) * 1e-10, - gsl_vector_get(c.axes[0], 1) * 1e-10, - gsl_vector_get(c.axes[0], 2) * 1e-10, - gsl_vector_get(c.axes[1], 0) * 1e-10, - gsl_vector_get(c.axes[1], 1) * 1e-10, - gsl_vector_get(c.axes[1], 2) * 1e-10, - gsl_vector_get(c.axes[2], 0) * 1e-10, - gsl_vector_get(c.axes[2], 1) * 1e-10, - gsl_vector_get(c.axes[2], 2) * 1e-10); + cell_set_cartesian(uc, gsl_vector_get(c->axes[0], 0) * 1e-10, + gsl_vector_get(c->axes[0], 1) * 1e-10, + gsl_vector_get(c->axes[0], 2) * 1e-10, + gsl_vector_get(c->axes[1], 0) * 1e-10, + gsl_vector_get(c->axes[1], 1) * 1e-10, + gsl_vector_get(c->axes[1], 2) * 1e-10, + gsl_vector_get(c->axes[2], 0) * 1e-10, + gsl_vector_get(c->axes[2], 1) * 1e-10, + gsl_vector_get(c->axes[2], 2) * 1e-10); if ( check_cell(dp, image, uc) ) { asdf_cell_free(c); -- cgit v1.2.3 From 10e85f43e59949eec1c8a479b37b029d8d495c71 Mon Sep 17 00:00:00 2001 From: Alexandra Tolstikova Date: Wed, 17 Jun 2015 12:49:04 +0200 Subject: Fix fft in[N] array subscript calculation --- libcrystfel/src/asdf.c | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/asdf.c b/libcrystfel/src/asdf.c index fb2f8279..da12e3c1 100644 --- a/libcrystfel/src/asdf.c +++ b/libcrystfel/src/asdf.c @@ -286,7 +286,7 @@ static int check_cell(struct asdf_private *dp, struct image *image, } -static int compare_doubles (const void *a, const void *b) +static int compare_doubles(const void *a, const void *b) { const double *da = (const double *) a; const double *db = (const double *) b; @@ -296,7 +296,7 @@ static int compare_doubles (const void *a, const void *b) /* Compares tvectors by length */ -static int compare_tvectors (const void *a, const void *b) +static int compare_tvectors(const void *a, const void *b) { struct tvector *ta = (struct tvector *) a; struct tvector *tb = (struct tvector *) b; @@ -357,7 +357,7 @@ static float find_ds_fft(double *projections, int N_projections, double d_max, memcpy(projections_sorted, projections, sizeof(double) * n); qsort(projections_sorted, n, sizeof(double), compare_doubles); - int i; + int i, k; int N = fftw.N; // number of points in fft calculation double *in = fftw.in; @@ -369,8 +369,11 @@ static float find_ds_fft(double *projections, int N_projections, double d_max, } for ( i = 0; i < n; i++ ) { - in[(int)((projections_sorted[i] - projections_sorted[0]) / - (projections_sorted[n - 1] - projections_sorted[0]) * N)] ++; + + k = (int)((projections_sorted[i] - projections_sorted[0]) / + (projections_sorted[n - 1] - projections_sorted[0]) * + (N - 1)); + in[k] ++; } fftw_execute_dft_r2c(p, in, out); -- cgit v1.2.3 From 704d778287ea73543cf0c7c7a32de97a818b6a09 Mon Sep 17 00:00:00 2001 From: Alexandra Tolstikova Date: Wed, 17 Jun 2015 18:08:07 +0200 Subject: Check c->n after cell reduction --- libcrystfel/src/asdf.c | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/asdf.c b/libcrystfel/src/asdf.c index da12e3c1..49d87d92 100644 --- a/libcrystfel/src/asdf.c +++ b/libcrystfel/src/asdf.c @@ -545,7 +545,7 @@ static int refine_asdf_cell(struct asdf_cell *c, gsl_vector **reflections, 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; @@ -778,7 +778,8 @@ static int create_cell(struct tvector tvec1, struct tvector tvec2, /* Index reflections with new cell axes */ check_refl_fitting_cell(c, reflections, N_reflections, IndexFit); - + if ( c->n < 6 ) return 0; + /* Refine cell until the number of fitting * reflections stops increasing */ int n = 0; @@ -1053,7 +1054,7 @@ static int index_refls(gsl_vector **reflections, int N_reflections, int run_asdf(struct image *image, IndexingPrivate *ipriv) -{ +{ int i, j; double LevelFit = 1./1000; -- cgit v1.2.3 From 091e3378a95edab85a6c655629df4270fb041b04 Mon Sep 17 00:00:00 2001 From: Alexandra Tolstikova Date: Thu, 18 Jun 2015 14:19:20 +0200 Subject: Fix bug in cell refinement --- libcrystfel/src/asdf.c | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/asdf.c b/libcrystfel/src/asdf.c index 49d87d92..d28fcf6a 100644 --- a/libcrystfel/src/asdf.c +++ b/libcrystfel/src/asdf.c @@ -778,13 +778,13 @@ static int create_cell(struct tvector tvec1, struct tvector tvec2, /* Index reflections with new cell axes */ check_refl_fitting_cell(c, reflections, N_reflections, IndexFit); - if ( c->n < 6 ) return 0; /* Refine cell until the number of fitting * reflections stops increasing */ int n = 0; int cell_correct = 1; - while ( c->n - n && cell_correct ) { + while ( c->n - n > 0 && cell_correct ) { + n = c->n; cell_correct = refine_asdf_cell(c, reflections, N_reflections, IndexFit); -- cgit v1.2.3