aboutsummaryrefslogtreecommitdiff
path: root/contrib/cubeit.c
diff options
context:
space:
mode:
Diffstat (limited to 'contrib/cubeit.c')
-rw-r--r--contrib/cubeit.c624
1 files changed, 624 insertions, 0 deletions
diff --git a/contrib/cubeit.c b/contrib/cubeit.c
new file mode 100644
index 00000000..33dce012
--- /dev/null
+++ b/contrib/cubeit.c
@@ -0,0 +1,624 @@
+/*
+ * 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 "../src/utils.h"
+#include "../src/hdf5-file.h"
+#include "../src/diffraction.h"
+#include "../src/render.h"
+#include "../src/symmetry.h"
+#include "../src/stream.h"
+#include "../src/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, &lt);
+
+ 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;
+}
generated by cgit v1.2.3 (git 2.25.1) at 2024-07-02 13:25:48 +0000