aboutsummaryrefslogtreecommitdiff
path: root/libcrystfel/src
diff options
context:
space:
mode:
Diffstat (limited to 'libcrystfel/src')
-rw-r--r--libcrystfel/src/detector.c4
-rw-r--r--libcrystfel/src/detector.h2
-rw-r--r--libcrystfel/src/reax.c15
-rw-r--r--libcrystfel/src/utils.c40
-rw-r--r--libcrystfel/src/utils.h17
5 files changed, 44 insertions, 34 deletions
diff --git a/libcrystfel/src/detector.c b/libcrystfel/src/detector.c
index 5cd0ce0f..174abefe 100644
--- a/libcrystfel/src/detector.c
+++ b/libcrystfel/src/detector.c
@@ -311,7 +311,7 @@ double get_tt(struct image *image, double fs, double ss, int *err)
}
-void record_image(struct image *image, int do_poisson)
+void record_image(struct image *image, int do_poisson, gsl_rng *rng)
{
int x, y;
double total_energy, energy_density;
@@ -371,7 +371,7 @@ void record_image(struct image *image, int do_poisson)
sa = proj_area / (dsq + Lsq);
if ( do_poisson ) {
- counts = poisson_noise(intensity * ph_per_e * sa);
+ counts = poisson_noise(rng, intensity * ph_per_e * sa);
} else {
cf = intensity * ph_per_e * sa;
counts = cf;
diff --git a/libcrystfel/src/detector.h b/libcrystfel/src/detector.h
index cfe0808a..eaa5618e 100644
--- a/libcrystfel/src/detector.h
+++ b/libcrystfel/src/detector.h
@@ -153,7 +153,7 @@ extern double get_tt(struct image *image, double xs, double ys, int *err);
extern int in_bad_region(struct detector *det, double fs, double ss);
-extern void record_image(struct image *image, int do_poisson);
+extern void record_image(struct image *image, int do_poisson, gsl_rng *rng);
extern struct panel *find_panel(struct detector *det, double fs, double ss);
diff --git a/libcrystfel/src/reax.c b/libcrystfel/src/reax.c
index 2ab4397a..110d8701 100644
--- a/libcrystfel/src/reax.c
+++ b/libcrystfel/src/reax.c
@@ -793,6 +793,13 @@ static int check_twinning(UnitCell *c1, UnitCell *c2, int verbose)
double bx, by, bz;
double cx, cy, cz;
+ gsl_rng *rng;
+
+ /* This is a rubbish RNG, but it serves for this purpose: nothing more
+ * than "I couldn't be bothered to think of n_trials sets of random
+ * indices". */
+ rng = gsl_rng_alloc(gsl_rng_taus2);
+
cell_get_reciprocal(c1, &asx, &asy, &asz,
&bsx, &bsy, &bsz,
&csx, &csy, &csz);
@@ -807,9 +814,9 @@ static int check_twinning(UnitCell *c1, UnitCell *c2, int verbose)
double dev;
signed int h2i, k2i, l2i;
- h = flat_noise(0, 10);
- k = flat_noise(0, 10);
- l = flat_noise(0, 10);
+ h = gsl_rng_uniform_int(rng, 10);
+ k = gsl_rng_uniform_int(rng, 10);
+ l = gsl_rng_uniform_int(rng, 10);
/* Position of this (randomly selected)
* reciprocal lattice point */
@@ -844,6 +851,8 @@ static int check_twinning(UnitCell *c1, UnitCell *c2, int verbose)
STATUS("%i duplicates.\n", n_dup);
}
+ gsl_rng_free(rng);
+
if ( n_dup > 10 ) return 1;
return 0;
}
diff --git a/libcrystfel/src/utils.c b/libcrystfel/src/utils.c
index 2f5107e3..4c85bf9d 100644
--- a/libcrystfel/src/utils.c
+++ b/libcrystfel/src/utils.c
@@ -3,11 +3,11 @@
*
* Utility stuff
*
- * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY,
- * a research centre of the Helmholtz Association.
+ * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY,
+ * a research centre of the Helmholtz Association.
*
* Authors:
- * 2009-2012 Thomas White <taw@physics.org>
+ * 2009-2014 Thomas White <taw@physics.org>
*
* This file is part of CrystFEL.
*
@@ -165,26 +165,26 @@ void progress_bar(int val, int total, const char *text)
}
-double random_flat(double max)
+double random_flat(gsl_rng *rng, double max)
{
- return max * (double)random()/RAND_MAX;
+ return max * gsl_rng_uniform(rng);
}
-double flat_noise(double expected, double width)
+double flat_noise(gsl_rng *rng, double expected, double width)
{
- double noise = random_flat(2.0*width);
+ double noise = random_flat(rng, 2.0*width);
return expected+noise-width;
}
-double gaussian_noise(double expected, double stddev)
+double gaussian_noise(gsl_rng *rng, double expected, double stddev)
{
double x1, x2, noise;
/* A uniformly distributed random number between 0 and 1 */
- x1 = ((double)random()/RAND_MAX);
- x2 = ((double)random()/RAND_MAX);
+ x1 = gsl_rng_uniform(rng);
+ x2 = gsl_rng_uniform(rng);
noise = sqrt(-2.0*log(x1)) * cos(2.0*M_PI*x2);
@@ -192,14 +192,14 @@ double gaussian_noise(double expected, double stddev)
}
-static int fake_poisson_noise(double expected)
+static int fake_poisson_noise(gsl_rng *rng, double expected)
{
- double rf = gaussian_noise(expected, sqrt(expected));
+ double rf = gaussian_noise(rng, expected, sqrt(expected));
return (int)rf;
}
-int poisson_noise(double expected)
+int poisson_noise(gsl_rng *rng, double expected)
{
double L;
int k = 0;
@@ -207,7 +207,7 @@ int poisson_noise(double expected)
/* For large values of the mean, we get big problems with arithmetic.
* In such cases, fall back on a Gaussian with the right variance. */
- if ( expected > 100.0 ) return fake_poisson_noise(expected);
+ if ( expected > 100.0 ) return fake_poisson_noise(rng, expected);
L = exp(-expected);
@@ -216,7 +216,7 @@ int poisson_noise(double expected)
double r;
k++;
- r = (double)random()/(double)RAND_MAX;
+ r = gsl_rng_uniform(rng);
p *= r;
} while ( p > L );
@@ -282,14 +282,14 @@ struct quaternion normalise_quaternion(struct quaternion q)
*
* Returns: a randomly generated, normalised, quaternion.
**/
-struct quaternion random_quaternion()
+struct quaternion random_quaternion(gsl_rng *rng)
{
struct quaternion q;
- q.w = 2.0*(double)random()/RAND_MAX - 1.0;
- q.x = 2.0*(double)random()/RAND_MAX - 1.0;
- q.y = 2.0*(double)random()/RAND_MAX - 1.0;
- q.z = 2.0*(double)random()/RAND_MAX - 1.0;
+ q.w = 2.0*gsl_rng_uniform(rng) - 1.0;
+ q.x = 2.0*gsl_rng_uniform(rng) - 1.0;
+ q.y = 2.0*gsl_rng_uniform(rng) - 1.0;
+ q.z = 2.0*gsl_rng_uniform(rng) - 1.0;
q = normalise_quaternion(q);
return q;
diff --git a/libcrystfel/src/utils.h b/libcrystfel/src/utils.h
index f4477717..81270916 100644
--- a/libcrystfel/src/utils.h
+++ b/libcrystfel/src/utils.h
@@ -3,11 +3,11 @@
*
* Utility stuff
*
- * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY,
- * a research centre of the Helmholtz Association.
+ * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY,
+ * a research centre of the Helmholtz Association.
*
* Authors:
- * 2009-2012 Thomas White <taw@physics.org>
+ * 2009-2014 Thomas White <taw@physics.org>
*
* This file is part of CrystFEL.
*
@@ -40,6 +40,7 @@
#include <pthread.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
+#include <gsl/gsl_rng.h>
#include "thread-pool.h"
@@ -89,7 +90,7 @@ struct quaternion {
extern struct quaternion normalise_quaternion(struct quaternion q);
extern double quaternion_modulus(struct quaternion q);
-extern struct quaternion random_quaternion(void);
+extern struct quaternion random_quaternion(gsl_rng *rng);
extern int quaternion_valid(struct quaternion q);
extern struct rvec quat_rot(struct rvec q, struct quaternion z);
@@ -114,10 +115,10 @@ extern int assplode(const char *a, const char *delims, char ***pbits,
AssplodeFlag flags);
extern void progress_bar(int val, int total, const char *text);
-extern double random_flat(double max);
-extern double flat_noise(double expected, double width);
-extern double gaussian_noise(double expected, double stddev);
-extern int poisson_noise(double expected);
+extern double random_flat(gsl_rng *rng, double max);
+extern double flat_noise(gsl_rng *rng, double expected, double width);
+extern double gaussian_noise(gsl_rng *rng, double expected, double stddev);
+extern int poisson_noise(gsl_rng *rng, double expected);
/* Keep these ones inline, to avoid function call overhead */
static inline struct quaternion invalid_quaternion(void)
generated by cgit v1.2.3 (git 2.25.1) at 2024-07-07 10:16:25 +0000