/* * facetron.c * * Profile fitting for coherent nanocrystallography * * (c) 2006-2010 Thomas White * * Part of CrystFEL - crystallography with a FEL * */ #ifdef HAVE_CONFIG_H #include #endif #include #include #include #include #include #include #include #include #include #include "utils.h" #include "hdf5-file.h" #include "symmetry.h" #include "reflections.h" #include "stream.h" #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; struct detector *det; const char *sym; }; static void show_help(const char *s) { printf("Syntax: %s [options]\n\n", s); printf( "Post-refinement and profile fitting for coherent nanocrystallography.\n" "\n" " -h, --help Display this help message.\n" "\n" " -i, --input= Specify the name of the input 'stream'.\n" " (must be a file, not e.g. stdin)\n" " -o, --output= Output filename. Default: facetron.hkl.\n" " -g. --geometry= Get detector geometry from file.\n" " -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" " -j Run analyses in parallel.\n"); } static void process_image(struct process_args *pargs) { struct hdfile *hdfile; struct image image; 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 = pargs->det; /* 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; //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; //} //hdf5_read(hdfile, &image, 1); free(image.data); cell_free(pargs->cell); if ( image.flags != NULL ) free(image.flags); //hdfile_close(hdfile); } static void *worker_thread(void *pargsv) { struct process_args *pargs = pargsv; int finish; do { int wakeup; process_image(pargs); 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 { 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); } while ( !wakeup ); } while ( !pargs->finish ); return NULL; } static void optimise_all(int nthreads, struct detector *det, const char *sym, FILE *fh, int config_basename, const char *prefix, int n_total_patterns) { pthread_t workers[MAX_THREADS]; struct process_args *worker_args[MAX_THREADS]; int worker_active[MAX_THREADS]; int i; int rval; int n_done = 0; /* Initialise worker arguments */ for ( i=0; ifilename = malloc(1024); worker_active[i] = 0; worker_args[i]->det = det; pthread_mutex_init(&worker_args[i]->control_mutex, NULL); worker_args[i]->sym = sym; } /* Start threads off */ for ( i=0; ifilename, 1023, "%s%s", prefix, filename); pargs->cell = cell; free(filename); 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); } } /* Keep threads busy until the end of the data */ do { int i; rval = 0; for ( i=0; icontrol_mutex); done = pargs->done; pthread_mutex_unlock(&pargs->control_mutex); if ( !done ) continue; n_done++; progress_bar(n_done, n_total_patterns, "Refining"); /* Get the next filename */ rval = find_chunk(fh, &cell, &filename); if ( rval == 1 ) break; if ( config_basename ) { char *tmp; tmp = strdup(basename(filename)); free(filename); filename = tmp; } snprintf(pargs->filename, 1023, "%s%s", prefix, filename); pargs->cell = cell; free(filename); /* Wake the thread up ... */ pthread_mutex_lock(&pargs->control_mutex); pargs->done = 0; pargs->start = 1; pthread_mutex_unlock(&pargs->control_mutex); } } while ( rval == 0 ); /* Join threads */ for ( i=0; icontrol_mutex); pargs->finish = 1; pthread_mutex_unlock(&pargs->control_mutex); /* Wait for it to join */ pthread_join(workers[i], NULL); free: if ( worker_args[i]->filename != NULL ) { free(worker_args[i]->filename); } } } int main(int argc, char *argv[]) { int c; char *infile = NULL; char *outfile = NULL; char *geomfile = NULL; char *prefix = NULL; char *sym = NULL; FILE *fh; int nthreads = 1; int config_basename = 0; int config_checkprefix = 1; struct detector *det; double *i_full; ReflItemList *obs; int i; int n_total_patterns; /* Long options */ const struct option longopts[] = { {"help", 0, NULL, 'h'}, {"input", 1, NULL, 'i'}, {"output", 1, NULL, 'o'}, {"geometry", 1, NULL, 'g'}, {"prefix", 1, NULL, 'x'}, {"basename", 0, &config_basename, 1}, {"no-check-prefix", 0, &config_checkprefix, 0}, {"symmetry", 1, NULL, 'y'}, {0, 0, NULL, 0} }; /* Short options */ while ((c = getopt_long(argc, argv, "hi:g:x:j:y:o:", 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 'y' : sym = strdup(optarg); break; case 'o' : outfile = strdup(optarg); break; case 0 : break; default : return 1; } } /* Sanitise input filename and open */ 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); /* Sanitise output filename */ if ( outfile == NULL ) { outfile = strdup("facetron.hkl"); } /* Sanitise prefix */ if ( prefix == NULL ) { prefix = strdup(""); } else { if ( config_checkprefix ) { prefix = check_prefix(prefix); } } /* Get detector geometry */ det = get_detector_geometry(geomfile); if ( det == NULL ) { ERROR("Failed to read detector geometry from '%s'\n", geomfile); return 1; } free(geomfile); /* Prepare for iteration */ i_full = new_list_intensity(); obs = new_items(); n_total_patterns = count_patterns(fh); STATUS("There are %i patterns to process\n", n_total_patterns); /* Iterate */ for ( i=0; i<10; i++ ) { STATUS("Post refinement iteration %i of 10\n", i+1); /* Refine the geometry of all patterns to get the best fit */ rewind(fh); optimise_all(nthreads, det, sym, fh, config_basename, prefix, n_total_patterns); /* Re-estimate all the full intensities */ } /* Output results */ write_reflections(outfile, obs, i_full, NULL, NULL, NULL); /* Clean up */ free(i_full); delete_items(obs); fclose(fh); free(sym); free(prefix); free(outfile); return 0; }