From 639fc7e226f63ceb4dd51f6c5cbe1bc05ba96a67 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Sat, 18 Sep 2010 15:42:05 +0200 Subject: cubeit: Add multi-threading --- src/cubeit.c | 495 +++++++++++++++++++++++++++++++---------------------------- 1 file changed, 260 insertions(+), 235 deletions(-) (limited to 'src/cubeit.c') 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 #include #include -#include +#include +#include #include -#include -#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=

Prefix filenames from input file with

.\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 Run 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; ycell, &ax, &ay, &az, &bx, &by, + &bz, &cx, &cy, &cz); - val = vals[xs*ys*z + xs*y + x]; + for ( x=0; x 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; ycontrol_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; ycontrol_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; yfinish ); -#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; ifilename = 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; ifilename, 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; xcontrol_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; i180 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; ifilename != 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; } -- cgit v1.2.3