diff options
author | Thomas White <taw@physics.org> | 2010-10-08 17:15:32 +0200 |
---|---|---|
committer | Thomas White <taw@physics.org> | 2012-02-22 15:27:02 +0100 |
commit | f7be30c4e727ad517bd3bf1ee3dde914c01910ea (patch) | |
tree | d441fc525ec184cf434f507eacbc42cec5611ebc | |
parent | 234aedb229f26074322d912066499b26fa079805 (diff) |
facetron: Multithread the integration, as well
-rw-r--r-- | src/facetron.c | 190 |
1 files changed, 143 insertions, 47 deletions
diff --git a/src/facetron.c b/src/facetron.c index 486cde6d..55acb645 100644 --- a/src/facetron.c +++ b/src/facetron.c @@ -30,6 +30,7 @@ #include "reflections.h" #include "stream.h" #include "geometry.h" +#include "peaks.h" #define MAX_THREADS (256) @@ -44,10 +45,15 @@ struct process_args int finish; int done; + /* Analysis routine */ + void (*func)(struct process_args *); + /* Analysis parameters */ const char *sym; + pthread_mutex_t *list_lock; /* Protects 'obs', 'i_full' and 'cts' */ ReflItemList *obs; double *i_full; + unsigned int *cts; }; @@ -73,36 +79,93 @@ static void show_help(const char *s) } -static void refine_image(struct image *image, ReflItemList *obs, double *i_full) +static void refine_image(struct process_args *pargs) { - //struct hdfile *hdfile; + /* Do, er, something. */ +} + + +static double partiality(struct image *image, + signed int h, signed int k, signed int l) +{ + return 1.0; +} + + +static void integrate_image(struct process_args *pargs) +{ + struct reflhit *spots; + int j, n; + struct hdfile *hdfile; + struct image *image = pargs->image; - 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; - - //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); + hdfile = hdfile_open(image->filename); + if ( hdfile == NULL ) { + ERROR("Couldn't open '%s'\n", image->filename); + return; + } else if ( hdfile_set_image(hdfile, "/data/data0") ) { + ERROR("Couldn't select path\n"); + hdfile_close(hdfile); + return; + } + + if ( hdf5_read(hdfile, pargs->image, 0) ) { + ERROR("Couldn't read '%s'\n", image->filename); + hdfile_close(hdfile); + return; + } + +goto skip; + + /* Figure out which spots should appear in this pattern, + * using a large divergence and bandwidth to avoid missing + * reflection tails. */ + spots = find_intersections(image, image->indexed_cell, + image->div, image->bw, &n, 0); + + /* For each reflection, estimate the partiality */ + for ( j=0; j<n; j++ ) { + + signed int h, k, l; + float i_partial; + double p; + float xc, yc; + + h = spots[j].h; + k = spots[j].k; + l = spots[j].l; + + /* Calculated partiality of this spot in this pattern */ + p = partiality(image, h, k, l); + + /* Don't attempt to use spots with very small + * partialities, since it won't be accurate. */ + if ( p < 0.1 ) continue; + + /* Actual measurement of this reflection from this + * pattern? */ + /* FIXME: Coordinates aren't whole numbers */ + if ( integrate_peak(image, spots[j].x, spots[j].y, + &xc, &yc, &i_partial, 1, 1) ) continue; + + pthread_mutex_lock(pargs->list_lock); + integrate_intensity(pargs->i_full, h, k, l, i_partial); + integrate_count(pargs->cts, h, k, l, 1); + if ( !find_item(pargs->obs, h, k, l) ) { + add_item(pargs->obs, h, k, l); + } + pthread_mutex_unlock(pargs->list_lock); + + } + +skip: free(image->data); if ( image->flags != NULL ) free(image->flags); - //hdfile_close(hdfile); + hdfile_close(hdfile); + //free(spots); } @@ -115,7 +178,7 @@ static void *worker_thread(void *pargsv) int wakeup; - refine_image(pargs->image, pargs->obs, pargs->i_full); + pargs->func(pargs); pthread_mutex_lock(&pargs->control_mutex); pargs->done = 1; @@ -139,12 +202,15 @@ static void *worker_thread(void *pargsv) } -static void refine_all(struct image *images, int n_total_patterns, - struct detector *det, const char *sym, - ReflItemList *obs, double *i_full, int nthreads) +static void munch_threads(struct image *images, int n_total_patterns, + struct detector *det, const char *sym, + ReflItemList *obs, double *i_full, unsigned int *cts, + int nthreads, void (*func)(struct process_args *), + const char *text) { pthread_t workers[MAX_THREADS]; struct process_args *worker_args[MAX_THREADS]; + pthread_mutex_t list_lock = PTHREAD_MUTEX_INITIALIZER; int worker_active[MAX_THREADS]; int i; int n_done = 0; @@ -157,6 +223,11 @@ static void refine_all(struct image *images, int n_total_patterns, worker_active[i] = 0; pthread_mutex_init(&worker_args[i]->control_mutex, NULL); worker_args[i]->sym = sym; + worker_args[i]->obs = obs; + worker_args[i]->i_full = i_full; + worker_args[i]->cts = cts; + worker_args[i]->list_lock = &list_lock; + worker_args[i]->func = func; } @@ -211,7 +282,7 @@ static void refine_all(struct image *images, int n_total_patterns, pargs->done = 0; n_done++; - progress_bar(n_done, n_total_patterns, "Refining"); + progress_bar(n_done, n_total_patterns, text); /* If there are no more patterns, "done" will remain * zero, so the last pattern will not be re-counted. */ @@ -243,43 +314,53 @@ static void refine_all(struct image *images, int n_total_patterns, if ( pargs->done ) { n_done++; - progress_bar(n_done, n_total_patterns, "Refining"); + progress_bar(n_done, n_total_patterns, text); } /* else this thread was not busy */ } + + for ( i=0; i<nthreads; i++ ) { + free(worker_args[i]); + } } -static double partiality(struct image *image, double x, double y) +static void refine_all(struct image *images, int n_total_patterns, + struct detector *det, const char *sym, + ReflItemList *obs, double *i_full, int nthreads) { - return 1.0; + munch_threads(images, n_total_patterns, det, sym, obs, i_full, NULL, + nthreads, refine_image, "Refining"); } static void estimate_full(struct image *images, int n_total_patterns, struct detector *det, const char *sym, - ReflItemList *obs, double *i_full) + ReflItemList *obs, double *i_full, int nthreads) { int i; + unsigned int *cts; - for ( i=0; i<n_total_patterns; i++ ) { + cts = new_list_count(); + clear_items(obs); - struct reflhit *spots; - struct image *image = &images[i]; - int j, n; + munch_threads(images, n_total_patterns, det, sym, obs, i_full, cts, + nthreads, integrate_image, "Integrating"); - /* Get the "spots" appearing in this pattern */ - spots = find_intersections(image, image->indexed_cell, - image->div, image->bw, &n, 0); + /* Divide the totals to get the means */ + for ( i=0; i<num_items(obs); i++ ) { - /* For each spot, estimate the partiality */ - for ( j=0; j<n; j++ ) { - double p = partiality(image, spots[j].x, spots[j].y); - } + struct refl_item *it; + double total; - progress_bar(i, n_total_patterns-1, "Integrating"); + it = get_item(obs, i); + total = lookup_intensity(i_full, it->h, it->k, it->l); + total /= lookup_count(cts, it->h, it->k, it->l); + set_intensity(i_full, it->h, it->k, it->l, total); } + + free(cts); } @@ -420,7 +501,7 @@ int main(int argc, char *argv[]) UnitCell *cell; char *filename; - char fnamereal[1024]; + char *fnamereal; if ( find_chunk(fh, &cell, &filename) == 1 ) { ERROR("Couldn't get all of the filenames and cells" @@ -437,11 +518,17 @@ int main(int argc, char *argv[]) free(filename); filename = tmp; } + fnamereal = malloc(1024); snprintf(fnamereal, 1023, "%s%s", prefix, filename); images[i].filename = fnamereal; images[i].div = 0.5e-3; images[i].bw = 0.001; + images[i].orientation.w = 1.0; + images[i].orientation.x = 0.0; + images[i].orientation.y = 0.0; + images[i].orientation.z = 0.0; + images[i].det = det; free(filename); @@ -452,7 +539,8 @@ int main(int argc, char *argv[]) free(prefix); /* Make initial estimates */ - estimate_full(images, n_total_patterns, det, sym, obs, i_full); + estimate_full(images, n_total_patterns, det, sym, obs, i_full, + nthreads); /* Iterate */ for ( i=0; i<n_iter; i++ ) { @@ -464,7 +552,8 @@ int main(int argc, char *argv[]) nthreads); /* Re-estimate all the full intensities */ - estimate_full(images, n_total_patterns, det, sym, obs, i_full); + estimate_full(images, n_total_patterns, det, sym, obs, i_full, + nthreads); } @@ -476,6 +565,13 @@ int main(int argc, char *argv[]) delete_items(obs); free(sym); free(outfile); + free(det->panels); + free(det); + for ( i=0; i<n_total_patterns; i++ ) { + cell_free(images[i].indexed_cell); + free(images[i].filename); + } + free(images); return 0; } |