cubeit: Add multi-threading

1 files changed, 260 insertions, 235 deletions

diff --git a/src/cubeit.c b/src/cubeit.c
index 18b20556..ff6b18c6 100644
--- a/src/cubeit.c
+++ b/src/cubeit.c
@@ -20,18 +20,39 @@
 #include <string.h>
 #include <unistd.h>
 #include <getopt.h>
-#include <errno.h>
+#include <pthread.h>
+#include <sys/time.h>
 #include <assert.h>
-#include <png.h>
 
-#include "image.h"
-#include "cell.h"
+#include "utils.h"
 #include "hdf5-file.h"
 #include "diffraction.h"
-#include "render.h"
 
 
-#define MAX_HITS (1024)
+#define MAX_THREADS (256)
+
+struct process_args
+{
+	char *filename;
+	int id;
+
+	/* Thread control */
+	pthread_mutex_t control_mutex;  /* Protects the scary stuff below */
+	int start;
+	int finish;
+	int done;
+
+	UnitCell *cell;
+	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;
+
+};
 
 
 static void show_help(const char *s)
@@ -48,76 +69,7 @@ static void show_help(const char *s)
 "  -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"
-);
-}
-
-
-static UnitCell *read_orientation_matrix(FILE *fh)
-{
-	float u, v, w;
-	struct rvec as, bs, cs;
-	UnitCell *cell;
-	char line[1024];
-
-	if ( fgets(line, 1023, fh) == NULL ) return NULL;
-	if ( sscanf(line, "astar = %f %f %f", &u, &v, &w) != 3 ) {
-		ERROR("Couldn't read a-star\n");
-		return NULL;
-	}
-	as.u = u*1e9;  as.v = v*1e9;  as.w = w*1e9;
-	if ( fgets(line, 1023, fh) == NULL ) return NULL;
-	if ( sscanf(line, "bstar = %f %f %f", &u, &v, &w) != 3 ) {
-		ERROR("Couldn't read b-star\n");
-		return NULL;
-	}
-	bs.u = u*1e9;  bs.v = v*1e9;  bs.w = w*1e9;
-	if ( fgets(line, 1023, fh) == NULL ) return NULL;
-	if ( sscanf(line, "cstar = %f %f %f", &u, &v, &w) != 3 ) {
-		ERROR("Couldn't read c-star\n");
-		return NULL;
-	}
-	cs.u = u*1e9;  cs.v = v*1e9;  cs.w = w*1e9;
-	cell = cell_new_from_axes(as, bs, cs);
-
-	return cell;
-}
-
-
-static int find_chunk(FILE *fh, UnitCell **cell, char **filename)
-{
-	char line[1024];
-	char *rval = NULL;
-
-	do {
-
-		rval = fgets(line, 1023, fh);
-		if ( rval == NULL ) continue;
-
-		chomp(line);
-
-		if ( strncmp(line, "Reflections from indexing", 25) != 0 ) {
-			continue;
-		}
-
-		*filename = strdup(line+29);
-
-		/* Skip two lines (while checking for errors) */
-		rval = fgets(line, 1023, fh);
-		if ( rval == NULL ) continue;
-		rval = fgets(line, 1023, fh);
-		if ( rval == NULL ) continue;
-
-		*cell = read_orientation_matrix(fh);
-		if ( *cell == NULL ) {
-			STATUS("Got filename but no cell for %s\n", *filename);
-			continue;
-		}
-
-		return 0;
-
-	} while ( rval != NULL );
-
-	return 1;
+"  -j <n>                     Run <n> analyses in parallel.\n");
 }
 
 
@@ -187,109 +139,141 @@ static void interpolate_onto_grid(double *vals, double v,
 }
 
 
-static void write_slice(const char *filename, double *vals, int z,
-                        int xs, int ys, int zs, double boost)
+static void process_image(struct process_args *pargs)
 {
-#ifdef HAVE_LIBPNG
-	FILE *fh;
-	png_structp png_ptr;
-	png_infop info_ptr;
-	png_bytep *row_pointers;
+	struct hdfile *hdfile;
+	struct image image;
+	double ax, ay, az;
+	double bx, by, bz;
+	double cx, cy, cz;
 	int x, y;
-	float max = 0.0;
-	int w, h;
+	const signed int ht = 3;
+	const signed int kt = 4;
+	const signed int lt = 5;
 
-	w = xs;
-	h = ys;
+	image.features = NULL;
+	image.data = NULL;
+	image.flags = NULL;
+	image.indexed_cell = NULL;
+	image.id = pargs->id;
+	image.filename = pargs->filename;
+	image.hits = NULL;
+	image.n_hits = 0;
+	image.det = NULL;
+
+	/* View head-on (unit cell is tilted) */
+	image.orientation.w = 1.0;
+	image.orientation.x = 0.0;
+	image.orientation.y = 0.0;
+	image.orientation.z = 0.0;
+
+	STATUS("Processing '%s'\n", pargs->filename);
+
+	hdfile = hdfile_open(pargs->filename);
+	if ( hdfile == NULL ) {
+		return;
+	} else if ( hdfile_set_first_image(hdfile, "/") ) {
+		ERROR("Couldn't select path\n");
+		hdfile_close(hdfile);
+		return;
+	}
 
-	for ( y=0; y<h; y++ ) {
-	for ( x=0; x<w; x++ ) {
+	hdf5_read(hdfile, &image, 1);
 
-		float val;
+	cell_get_cartesian(pargs->cell, &ax, &ay, &az, &bx, &by,
+	                                &bz, &cx, &cy, &cz);
 
-		val = vals[xs*ys*z + xs*y + x];
+	for ( x=0; x<image.width; x++ ) {
+	for ( y=0; y<image.height; y++ ) {
 
-		if ( val > max ) max = val;
+		double hd, kd, ld;
+		signed int h, k, l;
+		double dh, dk, dl;
+		struct rvec q;
 
-	}
-	}
+		q = get_q(&image, x, y, 1, 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 = (signed int)rint(hd);
+		k = (signed int)rint(kd);
+		l = (signed int)rint(ld);
+
+		if ( (h!=ht) || (k!=kt) || (l!=lt) ) continue;
+
+		dh = hd - h;
+		dk = kd - k;
+		dl = ld - l;
+
+		double v = image.data[x+image.width*y];
+
+		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);
 
-	fh = fopen(filename, "wb");
-	if ( !fh ) {
-		ERROR("Couldn't open output file.\n");
-		return;
-	}
-	png_ptr = png_create_write_struct(PNG_LIBPNG_VER_STRING,
-	                                  NULL, NULL, NULL);
-	if ( !png_ptr ) {
-		ERROR("Couldn't create PNG write structure.\n");
-		fclose(fh);
-		return;
-	}
-	info_ptr = png_create_info_struct(png_ptr);
-	if ( !info_ptr ) {
-		png_destroy_write_struct(&png_ptr, (png_infopp)NULL);
-		ERROR("Couldn't create PNG info structure.\n");
-		fclose(fh);
-		return;
 	}
-	if ( setjmp(png_jmpbuf(png_ptr)) ) {
-		png_destroy_write_struct(&png_ptr, &info_ptr);
-		fclose(fh);
-		ERROR("PNG write failed.\n");
-		return;
 	}
-	png_init_io(png_ptr, fh);
 
-	png_set_IHDR(png_ptr, info_ptr, w, h, 8,
-	             PNG_COLOR_TYPE_RGB, PNG_INTERLACE_NONE,
-	             PNG_COMPRESSION_TYPE_DEFAULT, PNG_FILTER_TYPE_DEFAULT);
+	if ( pargs->config_angles ) {
+
+		double asx, asy, asz;
+		double bsx, bsy, bsz;
+		double csx, csy, csz;
+		double ang;
+		int bin;
+
+		cell_get_reciprocal(pargs->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->vals_mutex);
+		pargs->angles[bin]++;
+		pthread_mutex_unlock(&pargs->vals_mutex);
 
-	row_pointers = malloc(h*sizeof(png_bytep *));
+	}
 
-	/* Write the image data */
-	max /= boost;
-	if ( max <= 6 ) { max = 10; }
+	free(image.data);
+	if ( image.flags != NULL ) free(image.flags);
+	hdfile_close(hdfile);
+}
 
-	for ( y=0; y<h; y++ ) {
 
-		row_pointers[y] = malloc(w*3);
+static void *worker_thread(void *pargsv)
+{
+	struct process_args *pargs = pargsv;
+	int finish;
 
-		for ( x=0; x<w; x++ ) {
+	do {
 
-			float r, g, b;
-			float val;
+		int wakeup;
 
-			val = vals[xs*ys*z + xs*y + x];
+		process_image(pargs);
 
-			render_scale(val, max, SCALE_COLOUR, &r, &g, &b);
-			row_pointers[y][3*x] = (png_byte)255*r;
-			row_pointers[y][3*x+1] = (png_byte)255*g;
-			row_pointers[y][3*x+2] = (png_byte)255*b;
+		pthread_mutex_lock(&pargs->control_mutex);
+		pargs->done = 1;
+		pargs->start = 0;
+		pthread_mutex_unlock(&pargs->control_mutex);
 
-		}
-	}
+		/* Go to sleep until told to exit or process next image */
+		do {
 
-	for ( y=0; y<h/2+1; y++ ) {
-		png_bytep scratch;
-		scratch = row_pointers[y];
-		row_pointers[y] = row_pointers[h-y-1];
-		row_pointers[h-y-1] = scratch;
-	}
+			pthread_mutex_lock(&pargs->control_mutex);
+			/* Either of these can result in the thread waking up */
+			wakeup = pargs->start || pargs->finish;
+			finish = pargs->finish;
+			pthread_mutex_unlock(&pargs->control_mutex);
+			usleep(20000);
 
-	png_set_rows(png_ptr, info_ptr, row_pointers);
-	png_write_png(png_ptr, info_ptr, PNG_TRANSFORM_IDENTITY, NULL);
+		} while ( !wakeup );
 
-	png_destroy_write_struct(&png_ptr, &info_ptr);
-	for ( y=0; y<h; y++ ) {
-		free(row_pointers[y]);
-	}
-	free(row_pointers);
-	fclose(fh);
+	} while ( !pargs->finish );
 
-#else
-	STATUS("No PNG support.\n");
-#endif
+	return NULL;
 }
 
 
@@ -297,22 +281,23 @@ int main(int argc, char *argv[])
 {
 	int c;
 	char *infile = NULL;
+	char *geomfile = NULL;
 	FILE *fh;
-	UnitCell *cell;
-	char *filename;
-	unsigned int angles[180];
-	int i;
+	char *rval = NULL;
+	int n_images;
 	char *prefix = NULL;
-	char *geomfile = NULL;
+	int nthreads = 1;
+	pthread_t workers[MAX_THREADS];
+	struct process_args *worker_args[MAX_THREADS];
+	int worker_active[MAX_THREADS];
 	int config_basename = 0;
 	int config_checkprefix = 1;
 	struct detector *det;
-	int config_angles = 0;
-	const signed int ht = 3;
-	const signed int kt = 4;
-	const signed int lt = 5;
-	const int gs = 16;
+	int i;
 	double *vals;
+	const int gs = 16;
+	unsigned int angles[180];
+	int config_angles = 0;
 
 	/* Long options */
 	const struct option longopts[] = {
@@ -326,7 +311,8 @@ int main(int argc, char *argv[])
 	};
 
 	/* Short options */
-	while ((c = getopt_long(argc, argv, "hi:g:x:", longopts, NULL)) != -1) {
+	while ((c = getopt_long(argc, argv, "hi:g:x:j:",
+	                        longopts, NULL)) != -1) {
 
 		switch (c) {
 		case 'h' :
@@ -345,6 +331,10 @@ int main(int argc, char *argv[])
 			prefix = strdup(optarg);
 			break;
 
+		case 'j' :
+			nthreads = atoi(optarg);
+			break;
+
 		case 0 :
 			break;
 
@@ -354,16 +344,20 @@ int main(int argc, char *argv[])
 
 	}
 
-	if ( infile == NULL ) infile = strdup("-");
+
+	if ( infile == NULL ) {
+		infile = strdup("-");
+	}
 	if ( strcmp(infile, "-") == 0 ) {
 		fh = stdin;
 	} else {
 		fh = fopen(infile, "r");
 	}
 	if ( fh == NULL ) {
-		ERROR("Couldn't open input stream '%s'\n", infile);
-		return ENOENT;
+		ERROR("Failed to open input file '%s'\n", infile);
+		return 1;
 	}
+	free(infile);
 
 	if ( prefix == NULL ) {
 		prefix = strdup("");
@@ -386,100 +380,131 @@ int main(int argc, char *argv[])
 	/* Initialise shape transform array */
 	vals = calloc(gs*gs*gs, sizeof(double));
 
-	/* Loop over all successfully indexed patterns */
-	while ( find_chunk(fh, &cell, &filename) == 0 ) {
+	if ( (nthreads == 0) || (nthreads > MAX_THREADS) ) {
+		ERROR("Invalid number of threads.\n");
+		return 1;
+	}
 
-		struct image image;
-		struct hdfile *hdfile;
-		double ang;
-		double ax, ay, az;
-		double bx, by, bz;
-		double cx, cy, cz;
-		unsigned int bin;
-		struct rvec q;
-		int x, y;
+	/* Initialise worker arguments */
+	for ( i=0; i<nthreads; i++ ) {
+
+		worker_args[i] = malloc(sizeof(struct process_args));
+		worker_args[i]->filename = malloc(1024);
+		worker_active[i] = 0;
+		worker_args[i]->xs = gs;
+		worker_args[i]->ys = gs;
+		worker_args[i]->zs = gs;
+		worker_args[i]->config_angles = config_angles;
+		worker_args[i]->vals = vals;
+		worker_args[i]->angles = angles;
+		pthread_mutex_init(&worker_args[i]->control_mutex, NULL);
+		pthread_mutex_init(&worker_args[i]->vals_mutex, NULL);
+		pthread_mutex_init(&worker_args[i]->angles_mutex, NULL);
+
+	}
+
+	n_images = 0;
+
+	/* Start threads off */
+	for ( i=0; i<nthreads; i++ ) {
+
+		char line[1024];
+		struct process_args *pargs;
+		int r;
+
+		pargs = worker_args[i];
+
+		rval = fgets(line, 1023, fh);
+		if ( rval == NULL ) continue;
+		chomp(line);
+		snprintf(pargs->filename, 1023, "%s%s", prefix, line);
 
-		STATUS("Processing '%s'\n", filename);
+		n_images++;
 
-		hdfile = hdfile_open(filename);
-		if ( hdfile == NULL ) {
-			return ENOENT;
-		} else if ( hdfile_set_image(hdfile, "/data/data0") ) {
-			ERROR("Couldn't select path\n");
-			return ENOENT;
+		pthread_mutex_lock(&pargs->control_mutex);
+		pargs->done = 0;
+		pargs->start = 1;
+		pargs->finish = 0;
+		pthread_mutex_unlock(&pargs->control_mutex);
+
+		worker_active[i] = 1;
+		r = pthread_create(&workers[i], NULL, worker_thread, pargs);
+		if ( r != 0 ) {
+			worker_active[i] = 0;
+			ERROR("Couldn't start thread %i\n", i);
 		}
 
-		hdf5_read(hdfile, &image, 0);
+	}
+
+	/* Keep threads busy until the end of the data */
+	do {
 
-		cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by,
-		                         &bz, &cx, &cy, &cz);
+		int i;
 
-		for ( x=0; x<image.width; x++ ) {
-		for ( y=0; y<image.height; y++ ) {
+		for ( i=0; i<nthreads; i++ ) {
 
-			double hd, kd, ld;
-			signed int h, k, l;
-			double dh, dk, dl;
+			char line[1024];
+			struct process_args *pargs;
+			int done;
 
-			q = get_q(&image, x, y, 1, NULL, 1.0/image.lambda);
+			/* Spend time working, not managing threads */
+			usleep(100000);
 
-			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;
+			/* Are we using this thread record at all? */
+			if ( !worker_active[i] ) continue;
 
-			h = (signed int)rint(hd);
-			k = (signed int)rint(kd);
-			l = (signed int)rint(ld);
+			/* Has the thread finished yet? */
+			pargs = worker_args[i];
+			pthread_mutex_lock(&pargs->control_mutex);
+			done = pargs->done;
+			pthread_mutex_unlock(&pargs->control_mutex);
+			if ( !done ) continue;
 
-			if ( !((h==ht) && (k==kt) && (l==lt)) ) continue;
+			/* Get the next filename */
+			rval = fgets(line, 1023, fh);
+			if ( rval == NULL ) break;
+			chomp(line);
+			snprintf(pargs->filename, 1023, "%s%s", prefix, line);
 
-			dh = hd - h;
-			dk = kd - k;
-			dl = ld - l;
+			n_images++;
 
-			double v = image.data[x+image.width*y];
+			STATUS("Done %i images\n", n_images);
 
-			interpolate_onto_grid(vals, v, gs, gs, gs, dh, dk, dl);
+			/* Wake the thread up ... */
+			pthread_mutex_lock(&pargs->control_mutex);
+			pargs->done = 0;
+			pargs->start = 1;
+			pthread_mutex_unlock(&pargs->control_mutex);
 
 		}
-		}
 
-		if ( config_angles ) {
+	} while ( rval != NULL );
 
-			double asx, asy, asz;
-			double bsx, bsy, bsz;
-			double csx, csy, csz;
+	/* Join threads */
+	for ( i=0; i<nthreads; i++ ) {
 
-			cell_get_reciprocal(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);
-			angles[bin]++;
+		if ( !worker_active[i] ) goto free;
 
-		}
 
-		hdfile_close(hdfile);
-		free(image.data);
-		free(image.flags);
+		/* Tell the thread to exit */
+		struct process_args *pargs = worker_args[i];
+		pthread_mutex_lock(&pargs->control_mutex);
+		pargs->finish = 1;
+		pthread_mutex_unlock(&pargs->control_mutex);
 
-	}
+		/* Wait for it to join */
+		pthread_join(workers[i], NULL);
 
-	for ( i=0; i<gs; i++ ) {
-		char line[64];
-		float boost;
-		snprintf(line, 63, "slice-%i.png", i);
-		for ( boost=1; boost<1000; boost+=50 ) {
-			write_slice(line, vals, i, gs, gs, gs, boost);
+	free:
+		if ( worker_args[i]->filename != NULL ) {
+			free(worker_args[i]->filename);
 		}
-	}
 
-	if ( config_angles ) {
-		for ( i=0; i<180; i++ ) {
-			STATUS("%i %i\n", i, angles[i]);
-		}
 	}
 
+	fclose(fh);
+
+	STATUS("There were %i images.\n", n_images);
+
 	return 0;
 }