From 3f7f7b74dba93b7ac105071927664d7b1f2c31de Mon Sep 17 00:00:00 2001
From: Thomas White <taw@physics.org>
Date: Thu, 6 Oct 2011 14:33:10 +0200
Subject: Move "cubeit" to contrib

---
 src/cubeit.c | 624 -----------------------------------------------------------
 1 file changed, 624 deletions(-)
 delete mode 100644 src/cubeit.c

(limited to 'src/cubeit.c')

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, &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;
-}
-- 
cgit v1.2.3