diff options
author | Thomas White <taw@physics.org> | 2010-10-08 14:44:21 +0200 |
---|---|---|
committer | Thomas White <taw@physics.org> | 2012-02-22 15:27:01 +0100 |
commit | 59b806107c0bd532d4f315422547f2aaede460cd (patch) | |
tree | a4fccf08bd46ea747a4e145613d58b1f60c96aff | |
parent | 8ed37ffb1dbdaeddb557216a5c50561bda8c40ca (diff) |
facetron: Lots of framework stuff
-rw-r--r-- | src/facetron.c | 206 | ||||
-rw-r--r-- | src/image.h | 19 |
2 files changed, 138 insertions, 87 deletions
diff --git a/src/facetron.c b/src/facetron.c index 9d05cb05..f7b801e6 100644 --- a/src/facetron.c +++ b/src/facetron.c @@ -29,14 +29,14 @@ #include "symmetry.h" #include "reflections.h" #include "stream.h" +#include "geometry.h" #define MAX_THREADS (256) struct process_args { - char *filename; - int id; + struct image *image; /* Thread control */ pthread_mutex_t control_mutex; /* Protects the scary stuff below */ @@ -44,9 +44,10 @@ struct process_args int finish; int done; - UnitCell *cell; - struct detector *det; + /* Analysis parameters */ const char *sym; + ReflItemList *obs; + double *i_full; }; @@ -65,30 +66,28 @@ 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" +" -y, --symmetry=<sym> Merge according to symmetry <sym>.\n" +" -n, --iterations=<n> Run <n> cycles of post-refinement.\n" +"\n" " -j <n> Run <n> analyses in parallel.\n"); } -static void process_image(struct process_args *pargs) +static void refine_image(struct image *image, ReflItemList *obs, double *i_full) { - 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; + //struct hdfile *hdfile; + + image->features = NULL; + image->data = NULL; + image->flags = NULL; + image->hits = NULL; + image->n_hits = 0; /* 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; + 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 ) { @@ -101,10 +100,8 @@ static void process_image(struct process_args *pargs) //hdf5_read(hdfile, &image, 1); - - free(image.data); - cell_free(pargs->cell); - if ( image.flags != NULL ) free(image.flags); + free(image->data); + if ( image->flags != NULL ) free(image->flags); //hdfile_close(hdfile); } @@ -118,11 +115,10 @@ static void *worker_thread(void *pargsv) int wakeup; - process_image(pargs); + refine_image(pargs->image, pargs->obs, pargs->i_full); 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 */ @@ -143,24 +139,22 @@ static void *worker_thread(void *pargsv) } -static void optimise_all(int nthreads, struct detector *det, const char *sym, - FILE *fh, int config_basename, const char *prefix, - int n_total_patterns) +static void refine_all(struct image *images, int n_total_patterns, + struct detector *det, const char *sym, + ReflItemList *obs, double *i_full, int nthreads) { 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; + int n_started = 0; - /* Initialise worker arguments */ + /* Initialise worker arguments with the unchanging data */ 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]->det = det; pthread_mutex_init(&worker_args[i]->control_mutex, NULL); worker_args[i]->sym = sym; @@ -171,25 +165,9 @@ static void optimise_all(int nthreads, struct detector *det, const char *sym, struct process_args *pargs; int r; - int rval; - char *filename; - UnitCell *cell; pargs = worker_args[i]; - - /* 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); + pargs->image = &images[n_started++]; pthread_mutex_lock(&pargs->control_mutex); pargs->done = 0; @@ -210,14 +188,11 @@ static void optimise_all(int nthreads, struct detector *det, const char *sym, do { int i; - rval = 0; for ( i=0; i<nthreads; i++ ) { struct process_args *pargs; int done; - char *filename; - UnitCell *cell; /* Spend time working, not managing threads */ usleep(100000); @@ -235,19 +210,9 @@ static void optimise_all(int nthreads, struct detector *det, const char *sym, 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); + /* Start work on the next pattern */ + if ( n_started == n_total_patterns ) break; + pargs->image = &images[n_started++]; /* Wake the thread up ... */ pthread_mutex_lock(&pargs->control_mutex); @@ -257,12 +222,12 @@ static void optimise_all(int nthreads, struct detector *det, const char *sym, } - } while ( rval == 0 ); + } while ( n_started < n_total_patterns ); /* Join threads */ for ( i=0; i<nthreads; i++ ) { - if ( !worker_active[i] ) goto free; + if ( !worker_active[i] ) continue; /* Tell the thread to exit */ struct process_args *pargs = worker_args[i]; @@ -273,11 +238,42 @@ static void optimise_all(int nthreads, struct detector *det, const char *sym, /* Wait for it to join */ pthread_join(workers[i], NULL); - free: - if ( worker_args[i]->filename != NULL ) { - free(worker_args[i]->filename); + n_done++; + progress_bar(n_done, n_total_patterns, "Refining"); + + } +} + + +static double partiality(struct image *image, double x, double y) +{ + return 1.0; +} + + +static void estimate_full(struct image *images, int n_total_patterns, + struct detector *det, const char *sym, + ReflItemList *obs, double *i_full) +{ + int i; + + for ( i=0; i<n_total_patterns; i++ ) { + + struct reflhit *spots; + struct image *image = &images[i]; + int j, n; + + /* Get the "spots" appearing in this pattern */ + spots = find_intersections(image, image->indexed_cell, + image->div, image->bw, &n, 0); + + /* For each spot, estimate the partiality */ + for ( j=0; j<n; j++ ) { + double p = partiality(image, spots[j].x, spots[j].y); } + progress_bar(i, n_total_patterns-1, "Integrating"); + } } @@ -299,6 +295,8 @@ int main(int argc, char *argv[]) ReflItemList *obs; int i; int n_total_patterns; + struct image *images; + int n_iter = 10; /* Long options */ const struct option longopts[] = { @@ -310,6 +308,7 @@ int main(int argc, char *argv[]) {"basename", 0, &config_basename, 1}, {"no-check-prefix", 0, &config_checkprefix, 0}, {"symmetry", 1, NULL, 'y'}, + {"iterations", 1, NULL, 'n'}, {0, 0, NULL, 0} }; @@ -347,6 +346,10 @@ int main(int argc, char *argv[]) outfile = strdup(optarg); break; + case 'n' : + n_iter = atoi(optarg); + break; + case 0 : break; @@ -400,18 +403,63 @@ int main(int argc, char *argv[]) n_total_patterns = count_patterns(fh); STATUS("There are %i patterns to process\n", n_total_patterns); + images = malloc(n_total_patterns * sizeof(struct image)); + if ( images == NULL ) { + ERROR("Couldn't allocate memory for images.\n"); + return 1; + } + + /* Fill in what we know about the images so far */ + rewind(fh); + for ( i=0; i<n_total_patterns; i++ ) { + + UnitCell *cell; + char *filename; + char fnamereal[1024]; + + if ( find_chunk(fh, &cell, &filename) == 1 ) { + ERROR("Couldn't get all of the filenames and cells" + " from the input stream.\n"); + return 1; + } + + images[i].indexed_cell = cell; + + /* Mangle the filename now */ + if ( config_basename ) { + char *tmp; + tmp = strdup(basename(filename)); + free(filename); + filename = tmp; + } + snprintf(fnamereal, 1023, "%s%s", prefix, filename); + + images[i].filename = fnamereal; + images[i].div = 0.5e-3; + images[i].bw = 0.001; + + free(filename); + + progress_bar(i, n_total_patterns-1, "Loading pattern data"); + + } + fclose(fh); + free(prefix); + + /* Make initial estimates */ + estimate_full(images, n_total_patterns, det, sym, obs, i_full); + /* Iterate */ - for ( i=0; i<10; i++ ) { + for ( i=0; i<n_iter; i++ ) { - STATUS("Post refinement iteration %i of 10\n", i+1); + STATUS("Post refinement iteration %i of %i\n", i+1, n_iter); /* 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); + refine_all(images, n_total_patterns, det, sym, obs, i_full, + nthreads); /* Re-estimate all the full intensities */ - + estimate_full(images, n_total_patterns, det, sym, obs, i_full); } @@ -421,9 +469,7 @@ int main(int argc, char *argv[]) /* Clean up */ free(i_full); delete_items(obs); - fclose(fh); free(sym); - free(prefix); free(outfile); return 0; diff --git a/src/image.h b/src/image.h index 613b0989..decbbab8 100644 --- a/src/image.h +++ b/src/image.h @@ -82,21 +82,26 @@ struct image { struct reflhit *hits; int n_hits; - int id; + int id; /* ID number of the thread + * handling this image */ + /* Information about the crystal */ struct quaternion orientation; + double m; /* Mosaicity in radians */ - /* Wavelength must always be given */ - double lambda; /* Wavelength in m */ - /* Incident intensity (if unknown, put 1.0) */ - double f0; - int f0_available; + /* Information about the radiation */ + double lambda; /* Wavelength in m */ + double div; /* Divergence in radians */ + double bw; /* Bandwidth as a fraction */ + double f0; /* Incident intensity */ + int f0_available; /* 0 if f0 wasn't available + * from the input. */ int width; int height; - ImageFeatureList *features; /* "Experimental" features */ + ImageFeatureList *features; /* DirAx auto-indexing low-level stuff */ int dirax_pty; |