diff options
Diffstat (limited to 'src/calibrate_detector.c')
| -rw-r--r-- | src/calibrate_detector.c | 479 |
1 files changed, 479 insertions, 0 deletions
diff --git a/src/calibrate_detector.c b/src/calibrate_detector.c new file mode 100644 index 00000000..3dc57e4b --- /dev/null +++ b/src/calibrate_detector.c @@ -0,0 +1,479 @@ +/* + * calibrate-detector.c + * + * Detector calibration + * + * (c) 2006-2010 Thomas White <taw@physics.org> + * + * Part of CrystFEL - crystallography with a FEL + * + */ + + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#define _GNU_SOURCE 1 +#include <stdarg.h> +#include <stdlib.h> +#include <stdio.h> +#include <string.h> +#include <unistd.h> +#include <getopt.h> +#include <pthread.h> +#include <sys/time.h> + +#include "utils.h" +#include "hdf5-file.h" +#include "filters.h" +#include "peaks.h" + + +#define INTEGRATION_RADIUS (10) +#define MAX_THREADS (96) + + +typedef enum +{ + SUM_THRESHOLD, + SUM_PEAKS +} SumMethod; + +struct process_args +{ + char *filename; + int id; + int config_cmfilter; + int config_noisefilter; + double *sum; + int w; + int h; + SumMethod sum_method; +}; + + +static void show_help(const char *s) +{ + printf("Syntax: %s [options]\n\n", s); + printf( +"Calibrate detector geometry from FEL diffraction images.\n" +"\n" +" -h, --help Display this help message.\n" +"\n" +" -i, --input=<filename> Specify file containing list of images to process.\n" +" '-' means stdin, which is the default.\n" +" -o, --output=<filename> Output filename for summed image in HDF5 format.\n" +" Default: summed.h5.\n" +"\n" +" -s, --sum=<method> Use this method for summation. Choose from:\n" +" peaks : sum 10px radius circles around peaks.\n" +" threshold : sum thresholded images.\n" +" --filter-cm Perform common-mode noise subtraction on images\n" +" before proceeding.\n" +" --filter-noise Apply an aggressive noise filter which sets all\n" +" pixels in each 3x3 region to zero if any of them\n" +" have negative values.\n" +" -j <n> Run <n> analyses in parallel. Default 1.\n"); +} + + +static void sum_peaks(struct image *image, double *sum) +{ + int x, y, i; + int w = image->width; + const int lim = INTEGRATION_RADIUS * INTEGRATION_RADIUS; + + search_peaks(image); + + for ( i=0; i<image_feature_count(image->features); i++ ) { + + struct imagefeature *f = image_get_feature(image->features, i); + int xp, yp; + + /* This is not an error. */ + if ( f == NULL ) continue; + + xp = f->x; + yp = f->y; + + for ( x=-INTEGRATION_RADIUS; x<+INTEGRATION_RADIUS; x++ ) { + for ( y=-INTEGRATION_RADIUS; y<+INTEGRATION_RADIUS; y++ ) { + + /* Circular mask */ + if ( x*x + y*y > lim ) continue; + + if ( ((x+xp)>=image->width) || ((x+xp)<0) ) continue; + if ( ((y+yp)>=image->height) || ((y+yp)<0) ) continue; + + float val = image->data[(x+xp)+w*(y+yp)]; + sum[(x+xp)+w*(y+yp)] += val; + + } + } + + } +} + + +static void sum_threshold(struct image *image, double *sum) +{ + int x, y; + + for ( x=0; x<image->width; x++ ) { + for ( y=0; y<image->height; y++ ) { + float val = image->data[x+image->width*y]; + if ( val > 400.0 ) { + sum[x+image->width*y] += val; + } + } + } +} + + +static void *process_image(void *pargsv) +{ + struct process_args *pargs = pargsv; + 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 = 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 NULL; + } else if ( hdfile_set_first_image(hdfile, "/") ) { + ERROR("Couldn't select path\n"); + hdfile_close(hdfile); + return NULL; + } + + hdf5_read(hdfile, &image); + + if ( pargs->config_cmfilter ) { + filter_cm(&image); + } + + if ( pargs->config_noisefilter ) { + filter_noise(&image, NULL); + } + + if ( (pargs->w != image.width) || (pargs->h != image.height) ) { + ERROR("Wrong image size.\n"); + goto out; + } + + switch ( pargs->sum_method ) { + + case SUM_THRESHOLD : + sum_threshold(&image, pargs->sum); + break; + + case SUM_PEAKS : + sum_peaks(&image, pargs->sum); + break; + + } + +out: + free(image.data); + if ( image.flags != NULL ) free(image.flags); + hdfile_close(hdfile); + + return NULL; +} + + +static void dump_to_file(struct process_args *worker_args[], int nthreads, + int w, int h, int n) +{ + int i; + double *total; + char outfile[256]; + + total = calloc(w*h, sizeof(double)); + + /* Add the individual sums to the 0th sum */ + for ( i=0; i<nthreads; i++ ) { + + int x, y; + + for ( x=0; x<w; x++ ) { + for ( y=0; y<h; y++ ) { + double val = worker_args[i]->sum[x+w*y]; + total[x+w*y] += val; + } + } + + } + + snprintf(outfile, 255, "sum-%i.h5", n); + + hdf5_write(outfile, total, w, h, H5T_NATIVE_DOUBLE); +} + + +int main(int argc, char *argv[]) +{ + int c; + char *filename = NULL; + char *outfile = NULL; + FILE *fh; + char *rval = NULL; + int n_images; + int config_cmfilter = 0; + int config_noisefilter = 0; + char *prefix = NULL; + char *sum_str = NULL; + SumMethod sum; + int nthreads = 1; + pthread_t workers[MAX_THREADS]; + struct process_args *worker_args[MAX_THREADS]; + int worker_active[MAX_THREADS]; + int i; + const int w = 1024; /* FIXME! */ + const int h = 1024; /* FIXME! */ + + /* Long options */ + const struct option longopts[] = { + {"help", 0, NULL, 'h'}, + {"input", 1, NULL, 'i'}, + {"output", 1, NULL, 'o'}, + {"filter-cm", 0, &config_cmfilter, 1}, + {"filter-noise", 0, &config_noisefilter, 1}, + {"prefix", 1, NULL, 'x'}, + {"sum", 1, NULL, 's'}, + {0, 0, NULL, 0} + }; + + /* Short options */ + while ((c = getopt_long(argc, argv, "hi:x:j:o:", + longopts, NULL)) != -1) { + + switch (c) { + case 'h' : + show_help(argv[0]); + return 0; + + case 'i' : + filename = strdup(optarg); + break; + + case 'o' : + outfile = strdup(optarg); + break; + + case 'x' : + prefix = strdup(optarg); + break; + + case 'j' : + nthreads = atoi(optarg); + break; + + case 's' : + sum_str = strdup(optarg); + break; + + case 0 : + break; + + default : + return 1; + } + + } + + if ( filename == NULL ) { + filename = strdup("-"); + } + if ( strcmp(filename, "-") == 0 ) { + fh = stdin; + } else { + fh = fopen(filename, "r"); + } + if ( fh == NULL ) { + ERROR("Failed to open input file '%s'\n", filename); + return 1; + } + free(filename); + + if ( sum_str == NULL ) { + STATUS("You didn't specify a summation method, so I'm using" + " the 'peaks' method, which gives the best results.\n"); + sum = SUM_PEAKS; + } else if ( strcmp(sum_str, "peaks") == 0 ) { + sum = SUM_PEAKS; + } else if ( strcmp(sum_str, "threshold") == 0) { + sum = SUM_THRESHOLD; + } else { + ERROR("Unrecognised summation method '%s'\n", sum_str); + return 1; + } + free(sum_str); + + if ( prefix == NULL ) { + prefix = strdup(""); + } + + if ( outfile == NULL ) { + outfile = strdup("summed.h5"); + } + + if ( (nthreads == 0) || (nthreads > MAX_THREADS) ) { + ERROR("Invalid number of threads.\n"); + return 1; + } + + /* Initialise worker arguments */ + for ( i=0; i<nthreads; i++ ) { + + worker_args[i] = malloc(sizeof(struct process_args)); + worker_args[i]->filename = malloc(1024); + worker_args[i]->sum = calloc(w*h, sizeof(double)); + worker_active[i] = 0; + + worker_args[i]->w = w; + worker_args[i]->h = h; + worker_args[i]->sum_method = sum; + + } + + n_images = 0; + + /* Initially, fire off the full number of threads */ + 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); + + n_images++; + + pargs->config_cmfilter = config_cmfilter; + pargs->config_noisefilter = config_noisefilter; + + worker_active[i] = 1; + r = pthread_create(&workers[i], NULL, process_image, pargs); + if ( r != 0 ) { + worker_active[i] = 0; + ERROR("Couldn't start thread %i\n", i); + } + + } + + /* Start new threads as old ones finish */ + do { + + int i; + + for ( i=0; i<nthreads; i++ ) { + + char line[1024]; + int r; + struct timespec t; + struct timeval tv; + struct process_args *pargs; + + if ( !worker_active[i] ) continue; + + pargs = worker_args[i]; + + gettimeofday(&tv, NULL); + t.tv_sec = tv.tv_sec; + t.tv_nsec = tv.tv_usec * 1000 + 20000; + + r = pthread_timedjoin_np(workers[i], NULL, &t); + if ( r != 0 ) continue; /* Not ready yet */ + + worker_active[i] = 0; + + rval = fgets(line, 1023, fh); + if ( rval == NULL ) break; + chomp(line); + snprintf(pargs->filename, 1023, "%s%s", prefix, line); + + worker_active[i] = 1; + r = pthread_create(&workers[i], NULL, process_image, pargs); + if ( r != 0 ) { + worker_active[i] = 0; + ERROR("Couldn't start thread %i\n", i); + } + + n_images++; + STATUS("Done %i images\n", n_images); + + if ( n_images % 1000 == 0 ) { + dump_to_file(worker_args, nthreads, w, h, n_images); + } + } + + } while ( rval != NULL ); + + /* Catch all remaining threads */ + for ( i=0; i<nthreads; i++ ) { + + if ( !worker_active[i] ) goto free; + + pthread_join(workers[i], NULL); + + worker_active[i] = 0; + + free: + if ( worker_args[i]->filename != NULL ) { + free(worker_args[i]->filename); + } + + } + + /* Add the individual sums to the 0th sum */ + for ( i=1; i<nthreads; i++ ) { + + int x, y; + + for ( x=0; x<w; x++ ) { + for ( y=0; y<h; y++ ) { + double val = worker_args[i]->sum[x+w*y]; + worker_args[0]->sum[x+w*y] += val; + } + } + free(worker_args[i]->sum); + free(worker_args[i]); + + } + + hdf5_write(outfile, worker_args[0]->sum, w, h, H5T_NATIVE_DOUBLE); + + free(worker_args[0]->sum); + free(worker_args[0]); + free(prefix); + free(outfile); + fclose(fh); + + STATUS("There were %i images.\n", n_images); + + return 0; +} |