diff options
Diffstat (limited to 'src/cubeit.c')
-rw-r--r-- | src/cubeit.c | 624 |
1 files changed, 0 insertions, 624 deletions
diff --git a/src/cubeit.c b/src/cubeit.c deleted file mode 100644 index 212be426..00000000 --- a/src/cubeit.c +++ /dev/null @@ -1,624 +0,0 @@ -/* - * cubeit.c - * - * "Full integration" of diffraction data - * - * (c) 2006-2011 Thomas White <taw@physics.org> - * - * Part of CrystFEL - crystallography with a FEL - * - */ - - -#ifdef HAVE_CONFIG_H -#include <config.h> -#endif - -#include <stdarg.h> -#include <stdlib.h> -#include <stdio.h> -#include <string.h> -#include <unistd.h> -#include <getopt.h> -#include <assert.h> -#include <png.h> -#include <fenv.h> -#include <pthread.h> -#include <libgen.h> -#include <cairo.h> - -#include "utils.h" -#include "hdf5-file.h" -#include "diffraction.h" -#include "render.h" -#include "symmetry.h" -#include "stream.h" -#include "thread-pool.h" - - -struct static_sum_args -{ - pthread_mutex_t *vals_mutex; /* Protects "vals" */ - double *vals; - int xs; - int ys; - int zs; - int config_angles; - pthread_mutex_t *angles_mutex; /* Protects "angles" */ - unsigned int *angles; - - pthread_mutex_t *cell_mutex; /* Protects "angles" */ - double *as; - double *bs; - double *cs; - double *als; - double *bes; - double *gas; - - struct detector *det; - char *element; - signed int ht; - signed int kt; - signed int lt; - SymOpList *sym; -}; - - -struct sum_args -{ - char *filename; - UnitCell *cell; - - struct static_sum_args static_args; -}; - - -struct queue_args -{ - FILE *fh; - char *prefix; - int config_basename; - - struct static_sum_args static_args; -}; - - -static void show_help(const char *s) -{ - printf("Syntax: %s [options]\n\n", s); - printf( -"'Full integration' of diffraction data.\n" -"\n" -" -h, --help Display this help message.\n" -"\n" -" -i, --input=<filename> Specify the name of the input stream.\n" -" Can be '-' for stdin.\n" -" -g. --geometry=<file> Get detector geometry from file.\n" -" -x, --prefix=<p> Prefix filenames from input file with <p>.\n" -" --basename Remove the directory parts of the filenames.\n" -" --no-check-prefix Don't attempt to correct the --prefix.\n" -" -j <n> Run <n> analyses in parallel.\n" -" -e, --image=<element> Use this image from the HDF5 file.\n" -" Example: /data/data0.\n" -" Default: The first one found.\n"); -} - - -static void add_to_mean(UnitCell *cell, double *ast, double *bst, double *cst, - double *alst, double *best, double *gast) -{ - double asx, asy, asz; - double bsx, bsy, bsz; - double csx, csy, csz; - - cell_get_reciprocal(cell, &asx, &asy, &asz, &bsx, &bsy, &bsz, - &csx, &csy, &csz); - *ast += modulus(asx, asy, asz); - *bst += modulus(bsx, bsy, bsz); - *cst += modulus(csx, csy, csz); - *alst += angle_between(bsx, bsy, bsz, csx, csy, csz); - *best += angle_between(asx, asy, asz, csx, csy, csz); - *gast += angle_between(asx, asy, asz, bsx, bsy, bsz); -} - - -static void interpolate_linear(double *vals, double v, - int xs, int ys, int zs, - int xv, int yv, double zv) -{ - int k; - double val1, val2; - float f, c; - - c = (zv+0.5)*(float)zs; - c -= 0.5; - k = floor(c); - f = c - (float)k; - assert(f >= 0.0); - assert(k+1 <= zs); - - val1 = v*(1.0-f); - val2 = v*f; - - /* Intensity may belong to the next reflection along */ - if ( k >= 0 ) { - vals[xs*ys*k + xs*yv + xv] += val1; - } - if ( k+1 < zs ) { - vals[xs*ys*(k+1) + xs*yv + xv] += val2; - } -} - - -static void interpolate_bilinear(double *vals, double v, - int xs, int ys, int zs, - int xv, double yv, double zv) -{ - int k; - double val1, val2; - float f, c; - - c = (yv+0.5)*(float)ys; - c -= 0.5; - k = floor(c); - f = c - (float)k; - assert(f >= 0.0); - assert(k+1 <= ys); - - val1 = v*(1.0-f); - val2 = v*f; - - /* Intensity may partially belong to the next reflection along */ - if ( k >= 0 ) { - interpolate_linear(vals, val1, xs, ys, zs, xv, k, zv); - } - if ( k+1 < ys ) { - interpolate_linear(vals, val2, xs, ys, zs, xv, k+1, zv); - } -} - - -static void interpolate_onto_grid(double *vals, double v, - int xs, int ys, int zs, - double xv, double yv, double zv) -{ - int k; - double val1, val2; - float f, c; - - c = (xv+0.5)*(float)xs; - c -= 0.5; - k = floor(c); - f = c - (float)k; - assert(f >= 0.0); - assert(k+1 <= xs); - - val1 = v*(1.0-f); - val2 = v*f; - - /* Intensity may partially belong to the next reflection along */ - if ( k >= 0 ) { - interpolate_bilinear(vals, val1, xs, ys, zs, k, yv, zv); - } - if ( k+1 < xs ) { - interpolate_bilinear(vals, val2, xs, ys, zs, k+1, yv, zv); - } -} - - -static void sum_image(void *pg, int cookie) -{ - struct sum_args *apargs = pg; - struct static_sum_args *pargs = &apargs->static_args; - struct hdfile *hdfile; - struct image image; - double ax, ay, az; - double bx, by, bz; - double cx, cy, cz; - int fs, ss; - - image.features = NULL; - image.data = NULL; - image.flags = NULL; - image.indexed_cell = NULL; - image.filename = apargs->filename; - image.det = pargs->det; - - STATUS("Processing '%s'\n", apargs->filename); - - hdfile = hdfile_open(apargs->filename); - if ( hdfile == NULL ) return; - if ( pargs->element != NULL ) { - - int r; - r = hdfile_set_image(hdfile, pargs->element); - if ( r ) { - ERROR("Couldn't select path '%s'\n", - pargs->element); - hdfile_close(hdfile); - return; - } - - } else { - - int r; - r = hdfile_set_first_image(hdfile, "/"); - if ( r ) { - ERROR("Couldn't select first path\n"); - hdfile_close(hdfile); - return; - } - - } - - hdf5_read(hdfile, &image, 1); - - cell_get_cartesian(apargs->cell, &ax, &ay, &az, &bx, &by, - &bz, &cx, &cy, &cz); - - fesetround(1); /* Round towards nearest */ - for ( fs=0; fs<image.width; fs++ ) { - for ( ss=0; ss<image.height; ss++ ) { - - double hd, kd, ld; - signed int h, k, l; - double dh, dk, dl; - struct rvec q; - //signed int ha, ka, la; - - q = get_q(&image, fs, ss, NULL, 1.0/image.lambda); - - hd = q.u * ax + q.v * ay + q.w * az; - kd = q.u * bx + q.v * by + q.w * bz; - ld = q.u * cx + q.v * cy + q.w * cz; - - h = lrint(hd); - k = lrint(kd); - l = lrint(ld); - - /* Uncomment this to use symmetry if you only want to look - * at a particular reflection. But it's really really slow. - * And wrong. To get useful information from symmetry - * averaging, the pattern must be transformed by the - * appropriate symmetry operator(s) to bring it into - * alignment. */ - //get_asymm(pargs->sym, h, k, l, &ha, &ka, &la); - //if ( (ha!=pargs->ht) || (ka!=pargs->kt) || (la!=pargs->lt) ) { - // continue; - //} - - dh = hd - h; - dk = kd - k; - dl = ld - l; - - double v = image.data[fs+image.width*ss]; - - pthread_mutex_lock(pargs->vals_mutex); - interpolate_onto_grid(pargs->vals, v, pargs->xs, pargs->ys, - pargs->zs, dh, dk, dl); - pthread_mutex_unlock(pargs->vals_mutex); - - } - } - - if ( pargs->config_angles ) { - - double asx, asy, asz; - double bsx, bsy, bsz; - double csx, csy, csz; - double ang; - int bin; - - cell_get_reciprocal(apargs->cell, &asx, &asy, &asz, - &bsx, &bsy, &bsz, - &csx, &csy, &csz); - ang = angle_between(csx, csy, csz, 0.0, 0.0, 1.0); - ang = rad2deg(ang); /* 0->180 deg */ - bin = rint(ang); - pthread_mutex_lock(pargs->angles_mutex); - pargs->angles[bin]++; - pthread_mutex_unlock(pargs->angles_mutex); - - } - - pthread_mutex_lock(pargs->cell_mutex); - add_to_mean(apargs->cell, pargs->as, pargs->bs, pargs->cs, - pargs->als, pargs->bes, pargs->gas); - pthread_mutex_unlock(pargs->cell_mutex); - - free(image.data); - cell_free(apargs->cell); - if ( image.flags != NULL ) free(image.flags); - hdfile_close(hdfile); - - free(apargs->filename); - free(apargs); -} - - -static void write_slice(const char *filename, double *vals, int z, - int xs, int ys, int zs, double boost, - double as, double bs, double ang) -{ - int x, y, zf; - float max = 0.0; - int zoom = 16; - double s = zoom * 30.0 / 1e9; - cairo_surface_t *surface; - cairo_t *c; - int w, h; - double xl, yl, xli; - - w = xs; h = ys; - - /* Find maximum value */ - for ( zf=0; zf<zs; zf++ ) { - for ( y=0; y<h; y++ ) { - for ( x=0; x<w; x++ ) { - float val = vals[xs*ys*zf + xs*y + x]; - if ( val > max ) max = val; - } - } - } - max /= boost; - - xl = s*as; - yl = s*bs*sin(ang); - xli = s*bs*cos(ang); - - surface = cairo_image_surface_create(CAIRO_FORMAT_ARGB32, - xs*xl + ys*xli, ys*yl); - - c = cairo_create(surface); - - cairo_scale(c, 1.0, -1.0); - cairo_translate(c, 0.0, -ys*yl); - - cairo_rectangle(c, 0.0, 0.0, xs*xl + ys*xli, ys*yl); - cairo_set_source_rgb(c, 1.0, 1.0, 1.0); - cairo_fill(c); - - for ( y=0; y<h; y++ ) { - for ( x=0; x<w; x++ ) { - - double r, g, b; - double val; - - val = vals[xs*ys*z + xs*y + x]; - - render_scale(val, max, SCALE_COLOUR, &r, &g, &b); - - cairo_new_path(c); - cairo_move_to(c, x*xl+y*xli, y*yl); - cairo_line_to(c, (x+1)*xl+y*xli, y*yl); - cairo_line_to(c, (x+1)*xl+(y+1)*xli, (y+1)*yl); - cairo_line_to(c, x*xl+(y+1)*xli, (y+1)*yl); - cairo_set_source_rgb(c, r, g, b); - cairo_fill(c); - cairo_stroke(c); - - } - } - - cairo_surface_write_to_png(surface, filename); - cairo_surface_destroy(surface); -} - - -static void *get_image(void *qp) -{ - struct sum_args *pargs; - struct queue_args *qargs = qp; - struct image image; - - /* Get the next filename */ - if ( read_chunk(qargs->fh, &image) == 1 ) { - if ( ferror(qargs->fh) ) { - ERROR("Stream read error.\n"); - } - return NULL; - } - - /* Won't be needing these */ - image_feature_list_free(image.features); - image.features = NULL; - reflist_free(image.reflections); - - - pargs = malloc(sizeof(struct sum_args)); - - if ( qargs->config_basename ) { - char *tmp; - tmp = safe_basename(image.filename); - free(image.filename); - image.filename = tmp; - } - - memcpy(&pargs->static_args, &qargs->static_args, - sizeof(struct static_sum_args)); - - pargs->cell = image.indexed_cell; - pargs->filename = malloc(1024); - snprintf(pargs->filename, 1023, "%s%s", qargs->prefix, image.filename); - free(image.filename); - - return pargs; -} - - -int main(int argc, char *argv[]) -{ - int c; - char *infile = NULL; - char *geomfile = NULL; - FILE *fh; - int n_images; - char *prefix = NULL; - int nthreads = 1; - int config_basename = 0; - int config_checkprefix = 1; - struct detector *det; - int i; - double *vals; - const int gs = 16; - unsigned int angles[180]; - int config_angles = 0; - signed int ht, kt, lt; - SymOpList *sym; - struct queue_args qargs; - pthread_mutex_t vals_mutex = PTHREAD_MUTEX_INITIALIZER; - pthread_mutex_t angles_mutex = PTHREAD_MUTEX_INITIALIZER; - pthread_mutex_t cell_mutex = PTHREAD_MUTEX_INITIALIZER; - double as; - double bs; - double cs; - double als; - double bes; - double gas; - char *element = NULL; - - /* Long options */ - const struct option longopts[] = { - {"help", 0, NULL, 'h'}, - {"input", 1, NULL, 'i'}, - {"geometry", 1, NULL, 'g'}, - {"prefix", 1, NULL, 'x'}, - {"basename", 0, &config_basename, 1}, - {"no-check-prefix", 0, &config_checkprefix, 0}, - {"image", 1, NULL, 'e'}, - {0, 0, NULL, 0} - }; - - /* Short options */ - while ((c = getopt_long(argc, argv, "hi:g:x:j:", - longopts, NULL)) != -1) { - - switch (c) { - case 'h' : - show_help(argv[0]); - return 0; - - case 'i' : - infile = strdup(optarg); - break; - - case 'g' : - geomfile = strdup(optarg); - break; - - case 'x' : - prefix = strdup(optarg); - break; - - case 'j' : - nthreads = atoi(optarg); - break; - - case 'e' : - element = strdup(optarg); - break; - - case 0 : - break; - - default : - return 1; - } - - } - - if ( infile == NULL ) { - infile = strdup("-"); - } - if ( strcmp(infile, "-") == 0 ) { - fh = stdin; - } else { - fh = fopen(infile, "r"); - } - if ( fh == NULL ) { - ERROR("Failed to open input file '%s'\n", infile); - return 1; - } - free(infile); - - if ( prefix == NULL ) { - prefix = strdup(""); - } else { - if ( config_checkprefix ) { - prefix = check_prefix(prefix); - } - } - - det = get_detector_geometry(geomfile); - if ( det == NULL ) { - ERROR("Failed to read detector geometry from '%s'\n", geomfile); - return 1; - } - free(geomfile); - - sym = get_pointgroup("6/mmm"); /* FIXME: Should be on command line */ - - /* Initialise histogram */ - for ( i=0; i<180; i++ ) angles[i] = 0; - - /* Initialise shape transform array */ - vals = calloc(gs*gs*gs, sizeof(double)); - - if ( nthreads == 0 ) { - ERROR("Invalid number of threads.\n"); - return 1; - } - - /* FIXME: Get indices on command line (or elsewhere) */ - get_asymm(sym, 3, 4, 5, &ht, &kt, <); - - as = 0.0; bs = 0.0; cs = 0.0; als = 0.0; bes = 0.0; gas = 0.0; - - qargs.fh = fh; - qargs.prefix = prefix; - qargs.config_basename = config_basename; - qargs.static_args.xs = gs; - qargs.static_args.ys = gs; - qargs.static_args.zs = gs; - qargs.static_args.config_angles = config_angles; - qargs.static_args.vals = vals; - qargs.static_args.angles = angles; - qargs.static_args.det = det; - qargs.static_args.vals_mutex = &vals_mutex; - qargs.static_args.angles_mutex = &angles_mutex; - qargs.static_args.ht = ht; - qargs.static_args.kt = kt; - qargs.static_args.lt = lt; - qargs.static_args.sym = sym; - qargs.static_args.cell_mutex = &cell_mutex; - qargs.static_args.as = &as; - qargs.static_args.bs = &bs; - qargs.static_args.cs = &cs; - qargs.static_args.als = &als; - qargs.static_args.bes = &bes; - qargs.static_args.gas = &gas; - qargs.static_args.element = element; - - n_images = run_threads(nthreads, sum_image, get_image, NULL, &qargs, 0, - 0, 0, 0); - - fclose(fh); - - for ( i=0; i<gs; i++ ) { - char line[64]; - float boost = 1.0; - snprintf(line, 63, "slice-%i.png", i); - write_slice(line, vals, i, gs, gs, gs, boost, - as/n_images, bs/n_images, gas/n_images); - } - - if ( config_angles ) { - for ( i=0; i<180; i++ ) { - STATUS("%i %i\n", i, angles[i]); - } - } - - STATUS("There were %i images.\n", n_images); - - return 0; -} |