aboutsummaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
authorThomas White <taw@bitwiz.org.uk>2010-09-18 15:42:05 +0200
committerThomas White <taw@physics.org>2012-02-22 15:26:58 +0100
commit639fc7e226f63ceb4dd51f6c5cbe1bc05ba96a67 (patch)
tree06b0bda8993b43fd9f2cece8f19528204563998a /src
parentd1f773c4d27ad22e8e6c9551383c36739f1bde0f (diff)
cubeit: Add multi-threading
Diffstat (limited to 'src')
-rw-r--r--src/cubeit.c495
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;
}