aboutsummaryrefslogtreecommitdiff
path: root/libcrystfel
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2015-07-01 16:26:40 +0200
committerThomas White <taw@physics.org>2015-07-03 16:35:40 +0200
commita669675a497fffdb575788340a1ca52bb7473bfe (patch)
tree54a3a95420c130ce141c9224389e4ad93511e88d /libcrystfel
parentc6be0929ad9e4cb7289215a5ef9f48444314e0f0 (diff)
Delete trailing whitespace
Diffstat (limited to 'libcrystfel')
-rw-r--r--libcrystfel/src/asdf.c614
1 files changed, 307 insertions, 307 deletions
diff --git a/libcrystfel/src/asdf.c b/libcrystfel/src/asdf.c
index b5436f7f..df6d7500 100644
--- a/libcrystfel/src/asdf.c
+++ b/libcrystfel/src/asdf.c
@@ -68,21 +68,21 @@ struct asdf_private {
struct tvector {
gsl_vector *t;
int n; // number of fitting reflections
- int *fits;
+ int *fits;
};
struct fftw_vars fftw_vars_new()
-{
+{
int N = 1024;
-
+
struct fftw_vars fftw;
-
+
fftw.N = N;
fftw.in = fftw_malloc(N * sizeof(double));
fftw.out = fftw_malloc(N * sizeof(fftw_complex));
fftw.p = fftw_plan_dft_r2c_1d(N, fftw.in, fftw.out, FFTW_MEASURE);
-
+
return fftw;
}
@@ -99,60 +99,60 @@ static void fftw_vars_free(struct fftw_vars fftw)
struct asdf_cell {
gsl_vector *axes[3];
gsl_vector *reciprocal[3];
-
+
int n; // number of fitting reflections
- double volume;
-
+ 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
+ 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);
t.n = 0;
t.fits = malloc(sizeof(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 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]);
}
-
+
free(c->reflections);
for ( i = 0; i < c->N_refls; i++ ) {
free(c->indices[i]);
}
free(c->indices);
free(c);
-
+
return 1;
}
-static struct asdf_cell *asdf_cell_new(int n)
+static struct asdf_cell *asdf_cell_new(int n)
{
struct asdf_cell *c;
c = malloc(sizeof(struct asdf_cell));
@@ -162,48 +162,48 @@ static struct asdf_cell *asdf_cell_new(int n)
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->reflections = malloc(sizeof(int) * n);
if (c->reflections == NULL) return NULL;
-
+
c->indices = malloc(sizeof(double *) * n);
if (c->indices == NULL) return NULL;
-
+
for ( i = 0; i < n; i++ ) {
c->indices[i] = malloc(sizeof(double) * 3);
if (c->indices[i] == NULL) return NULL;
}
-
+
c->n = 0;
-
+
c->acl = 0;
c->n_max = 0;
-
+
return c;
}
-static int asdf_cell_memcpy(struct asdf_cell *dest, struct asdf_cell *src)
+static int asdf_cell_memcpy(struct asdf_cell *dest, struct asdf_cell *src)
{
int i;
for ( i = 0; i < 3; i++ ) {
gsl_vector_memcpy(dest->axes[i], src->axes[i]);
gsl_vector_memcpy(dest->reciprocal[i], src->reciprocal[i]);
}
-
+
dest->volume = src->volume;
-
+
int n = src->N_refls;
dest->N_refls = n;
-
+
dest->n = src->n;
memcpy(dest->reflections, src->reflections, sizeof(int) * n);
-
+
for (i = 0; i < n; i++ ) {
memcpy(dest->indices[i], src->indices[i], sizeof(double) * 3);
}
-
+
dest->acl = src->acl;
dest->n_max = src->n_max;
return 1;
@@ -211,8 +211,8 @@ 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)
+static int cross_product(gsl_vector *vec1, gsl_vector *vec2,
+ gsl_vector **result)
{
double c1[3], c2[3], p[3];
int i;
@@ -233,34 +233,34 @@ 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)
+/* 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);
-
+ 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)
+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;
}
@@ -311,7 +311,7 @@ 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);
}
@@ -321,41 +321,41 @@ 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)) -
+ 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)
+static int calc_normal(gsl_vector *c1, gsl_vector *c2, gsl_vector *c3,
+ gsl_vector *normal)
{
gsl_vector *c12 = gsl_vector_alloc(3);
gsl_vector *c23 = gsl_vector_alloc(3);
gsl_vector *c31 = gsl_vector_alloc(3);
gsl_vector *res = gsl_vector_alloc(3);
-
+
cross_product(c1, c2, &c12);
cross_product(c2, c3, &c23);
cross_product(c3, c1, &c31);
-
+
int i;
for ( i = 0; i < 3; i++ ) {
gsl_vector_set(res, i, gsl_vector_get(c12, i) +
gsl_vector_get(c23, i) +
gsl_vector_get(c31, i));
- }
-
+ }
+
gsl_vector_free(c12);
gsl_vector_free(c23);
gsl_vector_free(c31);
-
+
double norm = gsl_blas_dnrm2(res);
if ( norm < 0.0001 ) {
gsl_vector_free(res);
@@ -365,42 +365,42 @@ static int calc_normal(gsl_vector *c1, gsl_vector *c2, gsl_vector *c3,
gsl_vector_memcpy(normal, res);
gsl_vector_free(res);
}
-
+
return 1;
}
static float find_ds_fft(double *projections, int N_projections, double d_max,
- struct fftw_vars fftw)
-{
- int n = N_projections;
+ struct fftw_vars fftw)
+{
+ int n = N_projections;
double projections_sorted[n];
- memcpy(projections_sorted, projections, sizeof(double) * n);
+ memcpy(projections_sorted, projections, sizeof(double) * n);
qsort(projections_sorted, n, sizeof(double), compare_doubles);
-
+
int i, k;
-
+
int N = fftw.N; // number of points in fft calculation
double *in = fftw.in;
fftw_complex *out = fftw.out;
fftw_plan p = fftw.p;
-
+
for ( i=0; i<N; i++ ) {
in[i] = 0;
}
-
+
for ( i = 0; i < n; i++ ) {
- k = (int)((projections_sorted[i] - projections_sorted[0]) /
- (projections_sorted[n - 1] - projections_sorted[0]) *
+ 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);
-
- int i_max = (int)(d_max * (projections_sorted[n - 1] -
+
+ int i_max = (int)(d_max * (projections_sorted[n - 1] -
projections_sorted[0]));
-
+
int d = 1;
double max = 0;
double a;
@@ -411,55 +411,55 @@ static float find_ds_fft(double *projections, int N_projections, double d_max,
d = i;
}
}
-
+
double ds = (projections_sorted[n - 1] - projections_sorted[0]) / d;
-
+
return ds;
}
/* Returns number of reflections fitting ds.
- * A projected reflection fits a one-dimensional lattice with elementary
- * lattice vector d* if its absolute distance to the nearest lattice
+ * 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)
+static int check_refl_fitting_ds(double *projections, int N_projections,
+ double ds, double LevelFit)
{
if ( ds == 0 ) return 0;
-
+
int i;
int n = 0;
for ( i = 0; i < N_projections; i++ ) {
- if ( fabs(projections[i] -
- ds * round(projections[i]/ds)) < LevelFit )
+ 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* */
-static float refine_ds(double *projections, int N_projections, double ds,
- double LevelFit, int *fits)
+/* 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 )
+ 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);
@@ -469,49 +469,49 @@ static float refine_ds(double *projections, int N_projections, double ds,
fits[i] = 0;
}
}
-
-
- gsl_fit_mul(indices, 1, fit_refls, 1, N_fits, &c1, &cov11,
+
+
+ gsl_fit_mul(indices, 1, fit_refls, 1, N_fits, &c1, &cov11,
&sumsq);
- N_new = check_refl_fitting_ds(projections, N_projections, c1,
+ 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,
+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( 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) *
+ 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
+
+ /* 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] <
+
+ if ( dist[0]*dist[0] + dist[1]*dist[1] + dist[2]*dist[2] <
IndexFit * IndexFit ) {
c->reflections[i] = 1;
c->n++;
@@ -519,57 +519,57 @@ static int check_refl_fitting_cell(struct asdf_cell *c,
c->reflections[i] = 0;
}
}
-
+
return 1;
}
-/* Returns 0 when refinement doesn't converge (i.e. all fitting reflections
+/* 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,
+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 *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 ( 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_set(r[j], n,
gsl_vector_get(reflections[i], j));
}
n++;
}
-
- gsl_multifit_linear_workspace *work = gsl_multifit_linear_alloc(c->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_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);
@@ -577,33 +577,33 @@ static int refine_asdf_cell(struct asdf_cell *c, gsl_vector **reflections,
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]) )
+
+ 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)
+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 *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);
@@ -615,7 +615,7 @@ static int reduce_asdf_cell(struct asdf_cell *cl)
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 ) {
@@ -623,7 +623,7 @@ static int reduce_asdf_cell(struct asdf_cell *cl)
gamma = 180 - gamma;
alpha = 180 - alpha;
}
-
+
gsl_vector_add(vb, va);
bb = gsl_blas_dnrm2(vb);
if ( bb < b ) {
@@ -635,17 +635,17 @@ static int reduce_asdf_cell(struct asdf_cell *cl)
gsl_vector_memcpy(cl->axes[0], vb);
}
changed = 1;
- }
+ }
}
-
+
if ( changed == 0 ) {
-
+
if ( beta < 90 ) {
- gsl_vector_scale(vc, -1);
+ gsl_vector_scale(vc, -1);
beta = 180 - beta;
alpha = 180 - alpha;
}
-
+
gsl_vector_add(vc, va);
cc = gsl_blas_dnrm2(vc);
if ( cc < c ) {
@@ -659,18 +659,18 @@ static int reduce_asdf_cell(struct asdf_cell *cl)
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);
+ gsl_vector_scale(vc, -1);
beta = 180 - beta;
alpha = 180 - alpha;
}
-
+
gsl_vector_add(vc, vb);
cc = gsl_blas_dnrm2(vc);
if ( cc < c ) {
@@ -684,55 +684,55 @@ static int reduce_asdf_cell(struct asdf_cell *cl)
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)
+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 ||
+
+ 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)
+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;
@@ -743,121 +743,121 @@ static int find_acl(struct tvector t1, struct tvector t2, struct tvector t3)
}
-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,
+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);
+ 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],
+ 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
+
+ /* Refine cell until the number of fitting
* reflections stops increasing */
int n = 0;
int cell_correct = 1;
while ( c->n - n > 0 && cell_correct ) {
-
+
n = c->n;
- cell_correct = refine_asdf_cell(c, reflections, N_reflections,
- IndexFit);
- check_refl_fitting_cell(c, reflections, N_reflections,
+ 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)
-{
+static int find_cell(struct tvector *tvectors, int N_tvectors, double IndexFit,
+ double volume_min, double volume_max, int n_max,
+ gsl_vector **reflections, int N_reflections,
+ struct asdf_cell *result)
+{
int i, j, k;
-
- /* Only tvectors with the number of fitting reflections > acl are
+
+ /* Only tvectors with the number of fitting reflections > acl are
* considered */
int acl = N_reflections < 18 ? 6 : N_reflections/3;
-
+
struct asdf_cell *c = asdf_cell_new(N_reflections);
if (c == NULL) {
ERROR("Failed to allocate asdf_cell in find_cell!\n");
return 0;
- }
-
+ }
+
/* Traversing a 3d array in slices perpendicular to the main diagonal */
int sl;
for ( sl = 0; sl < 3 * N_tvectors - 1; sl++ ) {
-
+
int i_min = sl < 2 * N_tvectors ? 0 : sl - 2 * N_tvectors;
int i_max = sl < N_tvectors ? sl : N_tvectors;
-
+
for ( i = i_min; i < i_max; i++) if (tvectors[i].n > acl ) {
-
- int j_min = sl - N_tvectors - 2 * i - 1 < 0 ?
+
+ int j_min = sl - N_tvectors - 2 * i - 1 < 0 ?
i + 1 : sl - N_tvectors - i;
- int j_max = sl - N_tvectors - i < 0 ?
+ int j_max = sl - N_tvectors - i < 0 ?
sl - i : N_tvectors;
-
- for ( j = j_min; j < j_max; j++ )
+
+ 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) )
+ check_cell_angles(tvectors[i].t,
+ tvectors[j].t,
+ tvectors[k].t, 0.99) )
{
-
- if ( !create_cell(tvectors[i],
- tvectors[j],
+
+ if ( !create_cell(tvectors[i],
+ tvectors[j],
tvectors[k],
c, IndexFit,
- volume_min,
+ volume_min,
volume_max,
reflections,
- N_reflections) )
+ N_reflections) )
{
break;
- }
-
- acl = find_acl(tvectors[i],
- tvectors[j],
+ }
+
+ acl = find_acl(tvectors[i],
+ tvectors[j],
tvectors[k]);
c->acl = acl;
c->n_max = n_max;
-
+
reduce_asdf_cell(c);
-
- /* If the new cell has more fitting
+
+ /* If the new cell has more fitting
* reflections save it to result */
if ( result->n < c->n ) {
- asdf_cell_memcpy(result, c);
+ asdf_cell_memcpy(result, c);
}
- acl++;
-
+ acl++;
+
if ( acl > n_max ) break;
- if ( tvectors[j].n <= acl ||
+ if ( tvectors[j].n <= acl ||
tvectors[i].n <= acl ) break;
}
}
@@ -866,18 +866,18 @@ static int find_cell(struct tvector *tvectors, int N_tvectors, double IndexFit,
}
if ( acl > n_max ) break;
}
-
+
asdf_cell_free(c);
-
+
if ( result->n ) return 1;
return 0;
}
-static void shuffle_triplets(int **triplets, int n)
+static void shuffle_triplets(int **triplets, int n)
{
int i, j;
- int t[3];
+ 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));
@@ -887,24 +887,24 @@ static void shuffle_triplets(int **triplets, int n)
}
-static int index_refls(gsl_vector **reflections, int N_reflections,
+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 fftw_vars fftw)
+ double LevelFit, double IndexFit, int i_max,
+ struct asdf_cell *c, struct fftw_vars fftw)
{
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) *
+ int N_triplets = N_reflections * (N_reflections - 1) *
(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++ ) {
@@ -915,7 +915,7 @@ static int index_refls(gsl_vector **reflections, int N_reflections,
" in index_refls!\n");
return 0;
}
-
+
triplets[l][0] = i;
triplets[l][1] = j;
triplets[l][2] = k;
@@ -926,173 +926,173 @@ static int index_refls(gsl_vector **reflections, int N_reflections,
/* 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 (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
-
+
+ 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) )
+ normal) )
{
-
+
/* Calculate projections of reflections to normal */
for ( k = 0; k < N_reflections; k++ ) {
- gsl_blas_ddot(normal, 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,
+ ds = find_ds_fft(projections, N_reflections, d_max,
fftw);
-
- /* Refine ds, write 1 to fits[i] if reflections[i]
+
+ /* Refine ds, write 1 to fits[i] if reflections[i]
* fits ds */
- ds = refine_ds(projections, N_reflections, ds, LevelFit,
+ ds = refine_ds(projections, N_reflections, ds, LevelFit,
fits);
-
+
/* n - number of reflections fitting ds */
- n = check_refl_fitting_ds(projections, N_reflections,
+ 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,
+
+ gsl_vector_memcpy(tvectors[N_tvectors].t,
normal);
- memcpy(tvectors[N_tvectors].fits, fits,
+ 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),
+ 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
+
+ /* 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,
+
+ find_cell(tvectors, N_tvectors, IndexFit, volume_min,
+ volume_max, n_max, reflections,
N_reflections, c);
-
+
if ( c->n > 4 * n_max / 5 ) {
break;
}
- }
- }
+ }
+ }
free(fits);
-
+
for ( i = 0; i < N_tvectors; i++ ) {
tvector_free(tvectors[i]);
}
free(tvectors);
-
+
for ( i = 0; i < N_triplets; i++ ) {
free(triplets[i]);
}
free(triplets);
-
+
gsl_vector_free(normal);
-
+
if ( c->n ) return 1;
-
+
return 0;
}
-int run_asdf(struct image *image, IndexingPrivate *ipriv)
-{
+int run_asdf(struct image *image, IndexingPrivate *ipriv)
+{
int i, j;
-
+
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
-
+
+ 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);
- double vtol = (dp->ltl[0] + dp->ltl[1] + dp->ltl[2]) / 100 +
+ 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);
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);
+ gsl_vector_set(reflections[i], 2, f->rz/1e10);
}
-
+
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,
+
+ j = index_refls(reflections, N_reflections, d_max, volume_min,
volume_max, LevelFit, IndexFit, i_max, c, dp->fftw);
-
+
for ( i = 0; i < N_reflections; i++ ) {
gsl_vector_free(reflections[i]);
}
-
+
if ( j ) {
UnitCell *uc;
uc = cell_new();
-
- cell_set_cartesian(uc, gsl_vector_get(c->axes[0], 0) * 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,
@@ -1101,16 +1101,16 @@ int run_asdf(struct image *image, IndexingPrivate *ipriv)
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) ) {
+
+ if ( check_cell(dp, image, uc) ) {
asdf_cell_free(c);
cell_free(uc);
return 1;
}
-
+
cell_free(uc);
}
-
+
asdf_cell_free(c);
return 0;
}
@@ -1142,9 +1142,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;
}