diff options
Diffstat (limited to 'src/utils.c')
| -rw-r--r-- | src/utils.c | 488 |
1 files changed, 0 insertions, 488 deletions
diff --git a/src/utils.c b/src/utils.c deleted file mode 100644 index 6b29627f..00000000 --- a/src/utils.c +++ /dev/null @@ -1,488 +0,0 @@ -/* - * utils.c - * - * Utility stuff - * - * (c) 2006-2010 Thomas White <taw@physics.org> - * - * Part of CrystFEL - crystallography with a FEL - * - */ - -#include <libgen.h> -#include <math.h> -#include <string.h> -#include <stdio.h> -#include <unistd.h> -#include <sys/types.h> -#include <sys/stat.h> -#include <gsl/gsl_matrix.h> -#include <gsl/gsl_vector.h> -#include <gsl/gsl_blas.h> - -#include "utils.h" -#include "image.h" - - -/** - * SECTION:utils - * @short_description: Miscellaneous utilities - * @title: Utilities - * @section_id: - * @see_also: - * @include: "utils.h" - * @Image: - * - * Wibble - */ - -/** - * show_matrix_eqn: - * @M: A matrix - * @v: A vector - * @r: The number of elements in @v and the side length of @M. - * - * Displays a matrix equation of the form @M.a = @v. - * - * @M must be square. - **/ -void show_matrix_eqn(gsl_matrix *M, gsl_vector *v, int r) -{ - int i, j; - - for ( i=0; i<r; i++ ) { - STATUS("[ "); - for ( j=0; j<r; j++ ) { - STATUS("%+9.3e ", gsl_matrix_get(M, i, j)); - } - STATUS("][ a%2i ] = [ %+9.3e ]\n", i, gsl_vector_get(v, i)); - } -} - - -size_t notrail(char *s) -{ - size_t i; - size_t munched = 0; - - for ( i=strlen(s)-1; i>=0; i-- ) { - if ( (s[i] == ' ') || (s[i] == '\t') ) { - s[i] = '\0'; - munched++; - } else { - return munched; - } - } - - return munched; -} - - -void chomp(char *s) -{ - size_t i; - - if ( !s ) return; - - for ( i=0; i<strlen(s); i++ ) { - if ( (s[i] == '\n') || (s[i] == '\r') ) { - s[i] = '\0'; - return; - } - } -} - - -void progress_bar(int val, int total, const char *text) -{ - double frac; - int n, i; - char s[1024]; - const int width = 50; - - if ( total == 0 ) return; - - if ( !isatty(STDERR_FILENO) ) return; - if ( tcgetpgrp(STDERR_FILENO) != getpgrp() ) return; - - frac = (double)val/total; - n = (int)(frac*width); - - for ( i=0; i<n; i++ ) s[i] = '='; - for ( i=n; i<width; i++ ) s[i] = '.'; - s[width] = '\0'; - - pthread_mutex_lock(&stderr_lock); - fprintf(stderr, "\r%s: |%s|", text, s); - if ( val == total ) fprintf(stderr, "\n"); - pthread_mutex_unlock(&stderr_lock); - - fflush(stdout); -} - - -double random_flat(double max) -{ - return max * (double)random()/RAND_MAX; -} - - -double flat_noise(double expected, double width) -{ - double noise = random_flat(2.0*width); - return expected+noise-width; -} - - -double gaussian_noise(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); - - noise = sqrt(-2.0*log(x1)) * cos(2.0*M_PI*x2); - - return expected + noise*stddev; -} - - -static int fake_poisson_noise(double expected) -{ - double rf = gaussian_noise(expected, sqrt(expected)); - return (int)rf; -} - - -int poisson_noise(double expected) -{ - double L; - int k = 0; - double p = 1.0; - - /* 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); - - L = exp(-expected); - - do { - - double r; - - k++; - r = (double)random()/(double)RAND_MAX; - p *= r; - - } while ( p > L ); - - return k - 1; -} - - -/** - * SECTION:quaternion - * @short_description: Simple quaternion handling - * @title: Quaternion - * @section_id: - * @see_also: - * @include: "utils.h" - * @Image: - * - * There is a simple quaternion structure in CrystFEL. At the moment, it is - * only used when simulating patterns, as an argument to cell_rotate() to - * orient the unit cell. - */ - -/** - * quaternion_modulus: - * @q: A %quaternion - * - * If a quaternion represents a pure rotation, its modulus should be unity. - * - * Returns: the modulus of the given quaternion. - **/ -double quaternion_modulus(struct quaternion q) -{ - return sqrt(q.w*q.w + q.x*q.x + q.y*q.y + q.z*q.z); -} - - -/** - * normalise_quaternion: - * @q: A %quaternion - * - * Rescales the quaternion such that its modulus is unity. - * - * Returns: the normalised version of @q - **/ -struct quaternion normalise_quaternion(struct quaternion q) -{ - double mod; - struct quaternion r; - - mod = quaternion_modulus(q); - - r.w = q.w / mod; - r.x = q.x / mod; - r.y = q.y / mod; - r.z = q.z / mod; - - return r; -} - - -/** - * random_quaternion: - * - * Returns: a randomly generated, normalised, quaternion. - **/ -struct quaternion random_quaternion() -{ - 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 = normalise_quaternion(q); - - return q; -} - - -/** - * quaternion_valid: - * @q: A %quaternion - * - * Checks if the given quaternion is normalised. - * - * This function performs a nasty floating point comparison of the form - * <code>(modulus > 0.999) && (modulus < 1.001)</code>, and so should not be - * relied upon to spot anything other than the most obvious input error. - * - * Returns: 1 if the quaternion is normalised, 0 if not. - **/ -int quaternion_valid(struct quaternion q) -{ - double qmod; - - qmod = quaternion_modulus(q); - - /* Modulus = 1 to within some tolerance? - * Nasty allowance for floating-point accuracy follows... */ - if ( (qmod > 0.999) && (qmod < 1.001) ) return 1; - - return 0; -} - - -/** - * quat_rot - * @q: A vector (in the form of an %rvec) - * @z: A %quaternion - * - * Rotates a vector according to a quaternion. - * - * Returns: A rotated version of @p. - **/ -struct rvec quat_rot(struct rvec q, struct quaternion z) -{ - struct rvec res; - double t01, t02, t03, t11, t12, t13, t22, t23, t33; - - t01 = z.w*z.x; - t02 = z.w*z.y; - t03 = z.w*z.z; - t11 = z.x*z.x; - t12 = z.x*z.y; - t13 = z.x*z.z; - t22 = z.y*z.y; - t23 = z.y*z.z; - t33 = z.z*z.z; - - res.u = (1.0 - 2.0 * (t22 + t33)) * q.u - + (2.0 * (t12 + t03)) * q.v - + (2.0 * (t13 - t02)) * q.w; - - res.v = (2.0 * (t12 - t03)) * q.u - + (1.0 - 2.0 * (t11 + t33)) * q.v - + (2.0 * (t01 + t23)) * q.w; - - res.w = (2.0 * (t02 + t13)) * q.u - + (2.0 * (t23 - t01)) * q.v - + (1.0 - 2.0 * (t11 + t22)) * q.w; - - return res; -} - - -/* Return non-zero if c is in delims */ -static int assplode_isdelim(const char c, const char *delims) -{ - size_t i; - for ( i=0; i<strlen(delims); i++ ) { - if ( c == delims[i] ) return 1; - } - return 0; -} - - -static int assplode_extract(char ***pbits, int n, size_t n_captured, - size_t start, const char *a) -{ - char **bits = *pbits; - bits = realloc(bits, sizeof(char *)*(n+1)); - bits[n] = malloc(n_captured+1); - memcpy(bits[n], a+start, n_captured); - bits[n][n_captured] = '\0'; - n++; - *pbits = bits; - return n; -} - - -/* Split the string 'a' using 'delims' as a zero-terminated list of - * deliminators. - * Store each segment in bits[0...n] where n is the number of segments and is - * the return value. pbits = &bits - * Each segment needs to be freed with free() when finished with. - * The array of bits also needs to be freed with free() when finished with, - * unless n=0 in which case bits==NULL - */ -int assplode(const char *a, const char *delims, char ***pbits, - AssplodeFlag flags) -{ - size_t i, start, n_captured; - int n, last_was_delim; - char **bits; - - n = 0; - i = 0; - n_captured = 0; - start = 0; - last_was_delim = 0; - bits = NULL; - while ( i < strlen(a) ) { - - if ( assplode_isdelim(a[i], delims) ) { - - if ( n_captured > 0 ) { - /* This is a deliminator after a sequence of - * non-deliminator chars */ - n = assplode_extract(&bits, n, n_captured, - start, a); - } - - n_captured = 0; - if ( (flags & ASSPLODE_DUPS) && last_was_delim ) { - n = assplode_extract(&bits, n, 0, start, a); - } - last_was_delim = 1; - - } else { - - if ( n_captured == 0 ) { - /* No characters currently found, so this is the - * start */ - start = i; - } - n_captured++; - last_was_delim = 0; - - } - - i++; - - } - /* Left over characters at the end? */ - if ( n_captured > 0 ) { - n = assplode_extract(&bits, n, n_captured, start, a); - } - - *pbits = bits; - return n; -} - - -char *check_prefix(char *prefix) -{ - int r; - struct stat statbuf; - char *new; - size_t len; - - /* Is "prefix" a directory? */ - r = stat(prefix, &statbuf); - if ( r != 0 ) { - /* "prefix" probably doesn't exist. This is fine - assume - * the user knows what they're doing, and that "prefix" - * suffixed with the actual filename will produce something - * sensible. */ - return prefix; - } - - if ( !S_ISDIR(statbuf.st_mode) ) { - /* Also fine, as above. */ - return prefix; - } - - /* Does the prefix end in a slash? */ - if ( prefix[strlen(prefix)-1] == '/' ) { - /* This looks sensible. */ - return prefix; - } - - STATUS("Your prefix ('%s') is a directory, but doesn't end" - " with a slash. I'm going to add it for you.\n", prefix); - STATUS("If this isn't what you want, run with --no-check-prefix.\n"); - len = strlen(prefix)+2; - new = malloc(len); - snprintf(new, len, "%s/", prefix); - free(prefix); - return new; -} - - -char *safe_basename(const char *in) -{ - int i; - char *cpy; - char *res; - - cpy = strdup(in); - - /* Get rid of any trailing slashes */ - for ( i=strlen(cpy)-1; i>0; i-- ) { - if ( cpy[i] == '/' ) { - cpy[i] = '\0'; - } else { - break; - } - } - - /* Find the base name */ - for ( i=strlen(cpy)-1; i>0; i-- ) { - if ( cpy[i] == '/' ) { - i++; - break; - } - } - - res = strdup(cpy+i); - /* If we didn't find a previous slash, i==0 so res==cpy */ - - free(cpy); - - return res; -} - - -#ifdef GSL_FUDGE -/* Force the linker to bring in CBLAS to make GSL happy */ -void utils_fudge_gslcblas() -{ - STATUS("%p\n", cblas_sgemm); -} -#endif |