diff options
| -rw-r--r-- | src/ambigator.c | 29 |
1 files changed, 24 insertions, 5 deletions
diff --git a/src/ambigator.c b/src/ambigator.c index 96309a8d..9aed3b3f 100644 --- a/src/ambigator.c +++ b/src/ambigator.c @@ -40,6 +40,9 @@ #include <unistd.h> #include <getopt.h> #include <assert.h> +#include <gsl/gsl_rng.h> +#include <gsl/gsl_permutation.h> +#include <gsl/gsl_randist.h> #include <image.h> #include <utils.h> @@ -310,10 +313,11 @@ struct cc_list static struct cc_list *calc_ccs(struct flist **crystals, int n_crystals, - int ncorr, SymOpList *amb) + int ncorr, SymOpList *amb, gsl_rng *rng) { struct cc_list *ccs; int i; + gsl_permutation *p; assert(n_crystals >= ncorr); ncorr++; /* Extra value at end for sentinel */ @@ -321,9 +325,12 @@ static struct cc_list *calc_ccs(struct flist **crystals, int n_crystals, ccs = malloc(n_crystals*sizeof(struct cc_list)); if ( ccs == NULL ) return NULL; + p = gsl_permutation_alloc(n_crystals); + gsl_permutation_init(p); + for ( i=0; i<n_crystals; i++ ) { - int j, k; + int k, l; ccs[i].ind = malloc(ncorr*sizeof(int)); ccs[i].cc = malloc(ncorr*sizeof(float)); @@ -335,11 +342,17 @@ static struct cc_list *calc_ccs(struct flist **crystals, int n_crystals, } k = 0; - for ( j=(i+1)%n_crystals; j!=i; j=(j+1)%n_crystals ) { + gsl_ran_shuffle(rng, p->data, n_crystals, sizeof(size_t)); + + for ( l=0; l<n_crystals; l++ ) { int n; + int j; float cc; + j = gsl_permutation_get(p, l); + if ( i == j ) continue; + cc = corr(crystals[i], crystals[j], &n, 0); if ( n < 4 ) continue; @@ -356,11 +369,15 @@ static struct cc_list *calc_ccs(struct flist **crystals, int n_crystals, if ( amb != NULL ) { k = 0; - for ( j=(i+1)%n_crystals; j!=i; j=(j+1)%n_crystals ) { + for ( l=0; l<n_crystals; l++ ) { int n; + int j; float cc; + j = gsl_permutation_get(p, l); + if ( i == j ) continue; + cc = corr(crystals[i], crystals[j], &n, 1); if ( n < 4 ) continue; @@ -380,6 +397,8 @@ static struct cc_list *calc_ccs(struct flist **crystals, int n_crystals, } + gsl_permutation_free(p); + return ccs; } @@ -713,7 +732,7 @@ int main(int argc, char *argv[]) } } - ccs = calc_ccs(crystals, n_crystals, ncorr, amb); + ccs = calc_ccs(crystals, n_crystals, ncorr, amb, rng); if ( ccs == NULL ) { ERROR("Failed to allocate CCs\n"); return 1; |