diff options
Diffstat (limited to 'libcrystfel')
-rw-r--r-- | libcrystfel/src/asdf.c | 309 |
1 files changed, 181 insertions, 128 deletions
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) { |