aboutsummaryrefslogtreecommitdiff
path: root/libcrystfel
diff options
context:
space:
mode:
authorAlexandra Tolstikova <alexandra.tolstikova@desy.de>2015-02-11 13:44:59 +0100
committerAlexandra Tolstikova <alexandra.tolstikova@desy.de>2015-02-11 13:44:59 +0100
commitf8d230451e7e0f71d5bf90d5f1079b89d19c8272 (patch)
treed8ffc5a24298393ab4a6e43b2e61ae440989e449 /libcrystfel
parenta80c4c5d0fbc26062084dde2c380ae5370dcd600 (diff)
asdf.c: codestyle fixed
Diffstat (limited to 'libcrystfel')
-rw-r--r--libcrystfel/src/asdf.c910
1 files changed, 469 insertions, 441 deletions
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;