diff options
author | Thomas White <taw@physics.org> | 2013-02-05 16:21:40 +0100 |
---|---|---|
committer | Thomas White <taw@physics.org> | 2013-02-05 16:21:40 +0100 |
commit | 2ce85a95d86e350b785d206405e97d7317672188 (patch) | |
tree | 5cecce161733acf1a276e5a0b55e5c22f5ffb20d | |
parent | 29cca07716b48f9e433087f5dbb202165b1897e1 (diff) |
Indexing pipeline - "done"!
-rw-r--r-- | libcrystfel/src/index.c | 5 | ||||
-rw-r--r-- | libcrystfel/src/index.h | 2 | ||||
-rw-r--r-- | libcrystfel/src/stream.c | 24 | ||||
-rw-r--r-- | libcrystfel/src/stream.h | 2 | ||||
-rw-r--r-- | src/im-sandbox.c | 129 | ||||
-rw-r--r-- | src/im-sandbox.h | 27 | ||||
-rw-r--r-- | src/indexamajig.c | 128 |
7 files changed, 147 insertions, 170 deletions
diff --git a/libcrystfel/src/index.c b/libcrystfel/src/index.c index 2370e26e..821be39a 100644 --- a/libcrystfel/src/index.c +++ b/libcrystfel/src/index.c @@ -241,7 +241,7 @@ static IndexingMethod set_axes(IndexingMethod a) } -IndexingMethod *build_indexer_list(const char *str, int *need_cell) +IndexingMethod *build_indexer_list(const char *str) { int n, i; char **methods; @@ -249,9 +249,6 @@ IndexingMethod *build_indexer_list(const char *str, int *need_cell) int tmp; int nmeth = 0; - if ( need_cell == NULL ) need_cell = &tmp; - *need_cell = 0; - n = assplode(str, ",-", &methods, ASSPLODE_NONE); list = malloc((n+1)*sizeof(IndexingMethod)); diff --git a/libcrystfel/src/index.h b/libcrystfel/src/index.h index e0dcb8e0..093e4ad6 100644 --- a/libcrystfel/src/index.h +++ b/libcrystfel/src/index.h @@ -82,7 +82,7 @@ typedef enum { **/ typedef void *IndexingPrivate; -extern IndexingMethod *build_indexer_list(const char *str, int *need_cell); +extern IndexingMethod *build_indexer_list(const char *str); extern IndexingPrivate **prepare_indexing(IndexingMethod *indm, UnitCell *cell, const char *filename, diff --git a/libcrystfel/src/stream.c b/libcrystfel/src/stream.c index 31fb8527..113dc359 100644 --- a/libcrystfel/src/stream.c +++ b/libcrystfel/src/stream.c @@ -38,6 +38,9 @@ #include <stdio.h> #include <string.h> #include <assert.h> +#include <sys/types.h> +#include <sys/stat.h> +#include <fcntl.h> #include "cell.h" #include "utils.h" @@ -236,6 +239,8 @@ void write_chunk(Stream *st, struct image *i, struct hdfile *hdfile, } fprintf(st->fh, CHUNK_END_MARKER"\n\n"); + + fflush(st->fh); } @@ -497,14 +502,14 @@ Stream *open_stream_for_read(const char *filename) } -Stream *open_stream_for_write(const char *filename) +Stream *open_stream_fd_for_write(int fd) { Stream *st; st = malloc(sizeof(struct _stream)); if ( st == NULL ) return NULL; - st->fh = fopen(filename, "w"); + st->fh = fdopen(fd, "w"); if ( st->fh == NULL ) { free(st); return NULL; @@ -520,6 +525,21 @@ Stream *open_stream_for_write(const char *filename) } + +Stream *open_stream_for_write(const char *filename) +{ + int fd; + + fd = open(filename, O_CREAT | O_TRUNC | O_WRONLY, S_IRUSR | S_IWUSR); + if ( fd == -1 ) { + ERROR("Failed to open stream.\n"); + return NULL; + } + + return open_stream_fd_for_write(fd); +} + + void close_stream(Stream *st) { fclose(st->fh); diff --git a/libcrystfel/src/stream.h b/libcrystfel/src/stream.h index 9340ca0d..4d561cd2 100644 --- a/libcrystfel/src/stream.h +++ b/libcrystfel/src/stream.h @@ -42,6 +42,8 @@ typedef struct _stream Stream; extern Stream *open_stream_for_read(const char *filename); extern Stream *open_stream_for_write(const char *filename); +extern Stream *open_stream_fd_for_write(int fd); + extern void close_stream(Stream *st); extern void write_chunk(Stream *st, struct image *image, struct hdfile *hdfile, diff --git a/src/im-sandbox.c b/src/im-sandbox.c index c8cceb2b..825f2d6a 100644 --- a/src/im-sandbox.c +++ b/src/im-sandbox.c @@ -73,14 +73,27 @@ #define STATS_EVERY_N_SECONDS (5) +/* Information about the indexing process for one pattern */ +struct pattern_args +{ + /* "Input" */ + char *filename; + + /* "Output" */ + int n_crystals; +}; + + struct sandbox { pthread_mutex_t lock; - int n_indexable; int n_processed; - int n_indexable_last_stats; + int n_hadcrystals; + int n_crystals; int n_processed_last_stats; + int n_hadcrystals_last_stats; + int n_crystals_last_stats; int t_last_stats; struct index_args *iargs; @@ -164,31 +177,29 @@ static char *get_pattern(FILE *fh, char **use_this_one_instead, static void process_image(const struct index_args *iargs, - struct pattern_args *pargs, FILE *ofh, + struct pattern_args *pargs, Stream *st, int cookie) { float *data_for_measurement; size_t data_size; - UnitCell *cell = iargs->cell; int config_cmfilter = iargs->config_cmfilter; int config_noisefilter = iargs->config_noisefilter; - int config_verbose = iargs->config_verbose; IndexingMethod *indm = iargs->indm; int check; struct hdfile *hdfile; struct image image; + int i; image.features = NULL; image.data = NULL; image.flags = NULL; - image.indexed_cell = NULL; image.copyme = iargs->copyme; - image.reflections = NULL; - image.n_saturated = 0; image.id = cookie; image.filename = pargs->filename; image.beam = iargs->beam; image.det = iargs->det; + image.crystals = NULL; + image.n_crystals = 0; hdfile = hdfile_open(image.filename); if ( hdfile == NULL ) return; @@ -291,28 +302,34 @@ static void process_image(const struct index_args *iargs, free(image.data); image.data = data_for_measurement; - /* Calculate orientation matrix (by magic) */ + /* Index the pattern */ + index_pattern(&image, indm, iargs->ipriv); + + /* Default beam parameters */ image.div = image.beam->divergence; image.bw = image.beam->bandwidth; - image.profile_radius = image.beam->profile_radius; - index_pattern(&image, cell, indm, iargs->cellr, - config_verbose, iargs->ipriv, - iargs->config_insane, iargs->tols); + /* Integrate each crystal's diffraction spots */ + for ( i=0; i<image.n_crystals; i++ ) { - if ( image.indexed_cell != NULL ) { + RefList *reflections; - pargs->indexable = 1; + /* Set default crystal parameter(s) */ + crystal_set_profile_radius(image.crystals[i], + image.beam->profile_radius); if ( iargs->integrate_found ) { - image.reflections = select_intersections(&image, - image.indexed_cell); + reflections = select_intersections(&image, + image.crystals[i]); } else { - image.reflections = find_intersections(&image, - image.indexed_cell); + reflections = find_intersections(&image, + image.crystals[i]); } - if (image.reflections != NULL) { + crystal_set_reflections(image.crystals[i], reflections); + + if ( reflections != NULL ) { + integrate_reflections(&image, iargs->config_closer, iargs->config_bgsub, @@ -323,18 +340,15 @@ static void process_image(const struct index_args *iargs, iargs->integrate_saturated); } - } else { - image.reflections = NULL; } - write_chunk(ofh, &image, hdfile, iargs->stream_flags); - fprintf(ofh, "END\n"); - fflush(ofh); + write_chunk(st, &image, hdfile, iargs->include_peaks, + iargs->include_reflections); - /* Only free cell if found */ - cell_free(image.indexed_cell); + for ( i=0; i<image.n_crystals; i++ ) { + crystal_free(image.crystals[i]); + } - reflist_free(image.reflections); free(image.data); if ( image.flags != NULL ) free(image.flags); image_feature_list_free(image.features); @@ -343,7 +357,7 @@ static void process_image(const struct index_args *iargs, static void run_work(const struct index_args *iargs, - int filename_pipe, int results_pipe, FILE *ofh, + int filename_pipe, int results_pipe, Stream *st, int cookie) { int allDone = 0; @@ -389,12 +403,12 @@ static void run_work(const struct index_args *iargs, } else { pargs.filename = line; - pargs.indexable = 0; + pargs.n_crystals = 0; - process_image(iargs, &pargs, ofh, cookie); + process_image(iargs, &pargs, st, cookie); /* Request another image */ - c = sprintf(buf, "%i\n", pargs.indexable); + c = sprintf(buf, "%i\n", pargs.n_crystals); w = write(results_pipe, buf, c); if ( w < 0 ) { ERROR("write P0\n"); @@ -406,7 +420,7 @@ static void run_work(const struct index_args *iargs, } - cleanup_indexing(iargs->ipriv); + cleanup_indexing(iargs->indm, iargs->ipriv); free(iargs->indm); free(iargs->ipriv); free_detector_geometry(iargs->det); @@ -577,7 +591,7 @@ static void start_worker_process(struct sandbox *sb, int slot) if ( p == 0 ) { - FILE *sfh; + Stream *st; int j; struct sigaction sa; int r; @@ -614,10 +628,10 @@ static void start_worker_process(struct sandbox *sb, int slot) close(filename_pipe[1]); close(result_pipe[0]); - sfh = fdopen(sb->stream_pipe_write[slot], "w"); + st = open_stream_fd_for_write(sb->stream_pipe_write[slot]); run_work(sb->iargs, filename_pipe[0], result_pipe[1], - sfh, slot); - fclose(sfh); + st, slot); + close_stream(st); //close(filename_pipe[0]); close(result_pipe[1]); @@ -697,7 +711,7 @@ static void handle_zombie(struct sandbox *sb) void create_sandbox(struct index_args *iargs, int n_proc, char *prefix, int config_basename, FILE *fh, char *use_this_one_instead, - FILE *ofh) + Stream *st) { int i; int allDone; @@ -712,13 +726,14 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix, return; } - sb->n_indexable = 0; sb->n_processed = 0; - sb->n_indexable_last_stats = 0; + sb->n_hadcrystals = 0; + sb->n_crystals = 0; sb->n_processed_last_stats = 0; + sb->n_hadcrystals_last_stats = 0; + sb->n_crystals_last_stats = 0; sb->t_last_stats = get_monotonic_seconds(); sb->n_proc = n_proc; - sb->ofh = ofh; sb->iargs = iargs; pthread_mutex_init(&sb->lock, NULL); @@ -895,8 +910,14 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix, ERROR("Invalid result '%s'\n", results); } } else { - sb->n_indexable += atoi(results); + + int nc = atoi(results); + sb->n_crystals += nc; + if ( nc > 0 ) { + sb->n_hadcrystals++; + } sb->n_processed++; + } /* Send next filename */ @@ -929,14 +950,19 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix, tNow = get_monotonic_seconds(); if ( tNow >= sb->t_last_stats+STATS_EVERY_N_SECONDS ) { - STATUS("%i out of %i indexed so far," - " %i out of %i since the last message.\n", - sb->n_indexable, sb->n_processed, - sb->n_indexable - sb->n_indexable_last_stats, - sb->n_processed - sb->n_processed_last_stats); + STATUS("Total so far: %i images processed, " + "%i had crystals, %i crystals overall. " + "Since the last message: %i images processed," + "%i had crystals, %i crystals overall.\n", + sb->n_processed, sb->n_hadcrystals, + sb->n_crystals, + sb->n_processed - sb->n_processed_last_stats, + sb->n_hadcrystals - sb->n_hadcrystals_last_stats, + sb->n_crystals - sb->n_crystals_last_stats); - sb->n_indexable_last_stats = sb->n_indexable; sb->n_processed_last_stats = sb->n_processed; + sb->n_hadcrystals_last_stats = sb->n_hadcrystals; + sb->n_crystals_last_stats = sb->n_crystals; sb->t_last_stats = tNow; } @@ -975,9 +1001,8 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix, free(sb->pids); free(sb->running); - if ( ofh != stdout ) fclose(ofh); - - STATUS("There were %i images, of which %i could be indexed.\n", - sb->n_processed, sb->n_indexable); + STATUS("Final:" + " %i images processed, %i had crystals, %i crystals overall.\n", + sb->n_processed, sb->n_hadcrystals, sb->n_crystals); } diff --git a/src/im-sandbox.h b/src/im-sandbox.h index 872f84ad..176bfde0 100644 --- a/src/im-sandbox.h +++ b/src/im-sandbox.h @@ -3,13 +3,13 @@ * * Sandbox for indexing * - * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY, - * a research centre of the Helmholtz Association. + * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. * Copyright © 2012 Richard Kirian * Copyright © 2012 Lorenzo Galli * * Authors: - * 2010-2012 Thomas White <taw@physics.org> + * 2010-2013 Thomas White <taw@physics.org> * 2011 Richard Kirian * 2012 Lorenzo Galli * 2012 Chunhong Yoon @@ -31,6 +31,8 @@ * */ +#include "stream.h" + enum { PEAK_ZAEF, PEAK_HDF5, @@ -44,10 +46,8 @@ struct index_args int config_cmfilter; int config_noisefilter; int config_verbose; - int stream_flags; /* What goes into the output? */ int config_satcorr; int config_closer; - int config_insane; int config_bgsub; float threshold; float min_gradient; @@ -57,7 +57,6 @@ struct index_args IndexingMethod *indm; IndexingPrivate **ipriv; int peaks; /* Peak detection method */ - int cellr; float tols[4]; struct beam_params *beam; char *element; @@ -70,21 +69,11 @@ struct index_args int use_saturated; int no_revalidate; int integrate_found; + int include_peaks; + int include_reflections; }; - - -/* Information about the indexing process for one pattern */ -struct pattern_args -{ - /* "Input" */ - char *filename; - - /* "Output" */ - int indexable; -}; - extern void create_sandbox(struct index_args *iargs, int n_proc, char *prefix, int config_basename, FILE *fh, - char *use_this_one_instead, FILE *ofh); + char *use_this_one_instead, Stream *st); diff --git a/src/indexamajig.c b/src/indexamajig.c index 01377de3..3076ea86 100644 --- a/src/indexamajig.c +++ b/src/indexamajig.c @@ -174,37 +174,13 @@ static void show_help(const char *s) } -static int parse_cell_reduction(const char *scellr, int *err, - int *reduction_needs_cell) -{ - *err = 0; - if ( strcmp(scellr, "none") == 0 ) { - *reduction_needs_cell = 0; - return CELLR_NONE; - } else if ( strcmp(scellr, "reduce") == 0) { - *reduction_needs_cell = 1; - return CELLR_REDUCE; - } else if ( strcmp(scellr, "compare") == 0) { - *reduction_needs_cell = 1; - return CELLR_COMPARE; - } else if ( strcmp(scellr, "compare_ab") == 0) { - *reduction_needs_cell = 1; - return CELLR_COMPARE_AB; - } else { - *err = 1; - *reduction_needs_cell = 0; - return CELLR_NONE; - } -} - - int main(int argc, char *argv[]) { int c; char *filename = NULL; char *outfile = NULL; FILE *fh; - FILE *ofh; + Stream *st; char *rval = NULL; int config_noindex = 0; int config_cmfilter = 0; @@ -213,7 +189,6 @@ int main(int argc, char *argv[]) int config_satcorr = 1; int config_checkprefix = 1; int config_closer = 0; - int config_insane = 0; int config_bgsub = 1; int config_basename = 0; float threshold = 800.0; @@ -224,17 +199,13 @@ int main(int argc, char *argv[]) char *geometry = NULL; IndexingMethod *indm; IndexingPrivate **ipriv; - int indexer_needs_cell; - int reduction_needs_cell; char *indm_str = NULL; UnitCell *cell; char *pdb = NULL; char *prefix = NULL; char *speaks = NULL; - char *scellr = NULL; char *toler = NULL; float tols[4] = {5.0, 5.0, 5.0, 1.5}; /* a,b,c,angles (%,%,%,deg) */ - int cellr; int peaks; int n_proc = 1; char *prepare_line; @@ -243,7 +214,6 @@ int main(int argc, char *argv[]) struct index_args iargs; struct beam_params *beam = NULL; char *element = NULL; - int stream_flags = STREAM_INTEGRATED; char *hdf5_peak_path = NULL; struct copy_hdf5_field *copyme; char *intrad = NULL; @@ -255,6 +225,10 @@ int main(int argc, char *argv[]) int no_revalidate = 0; int integrate_found = 0; + /* For ease of upgrading from old method */ + char *scellr = NULL; + int config_insane = 0; + copyme = new_copy_hdf5_field_list(); if ( copyme == NULL ) { ERROR("Couldn't allocate HDF5 field list.\n"); @@ -290,7 +264,6 @@ int main(int argc, char *argv[]) {"peaks", 1, NULL, 2}, {"cell-reduction", 1, NULL, 3}, {"min-gradient", 1, NULL, 4}, - {"record", 1, NULL, 5}, {"cpus", 1, NULL, 6}, {"cpugroup", 1, NULL, 7}, {"cpuoffset", 1, NULL, 8}, @@ -301,6 +274,8 @@ int main(int argc, char *argv[]) {"tolerance", 1, NULL, 13}, {"int-radius", 1, NULL, 14}, + /* FIXME: Add '--no-peaks' and '--no-reflections' */ + {"integrate-saturated",0, &integrate_saturated,1}, {"use-saturated", 0, &use_saturated, 1}, {"no-revalidate", 0, &no_revalidate, 1}, @@ -376,8 +351,7 @@ int main(int argc, char *argv[]) break; case 5 : - stream_flags = parse_stream_flags(optarg); - if ( stream_flags < 0 ) return 1; + ERROR("The option '--record' is no longer used.\n"); break; case 6 : @@ -425,6 +399,12 @@ int main(int argc, char *argv[]) } + if ( scellr != NULL ) { + ERROR("Old-style indexing method specification.\n"); + /* FIXME: Translate to new style */ + return 1; + } + if ( filename == NULL ) { filename = strdup("-"); } @@ -439,17 +419,6 @@ int main(int argc, char *argv[]) } free(filename); - if ( outfile == NULL ) { - ofh = stdout; - } else { - ofh = fopen(outfile, "w"); - if ( ofh == NULL ) { - ERROR("Failed to open output file '%s'\n", outfile); - return 1; - } - free(outfile); - } - if ( hdf5_peak_path == NULL ) { hdf5_peak_path = strdup("/processing/hitfinder/peakinfo"); } @@ -482,53 +451,24 @@ int main(int argc, char *argv[]) return 1; } - if ( (indm_str == NULL) || - ((indm_str != NULL) && (strcmp(indm_str, "none") == 0)) ) { - STATUS("Not indexing anything.\n"); - indexer_needs_cell = 0; - reduction_needs_cell = 0; + if ( indm_str == NULL ) { + + STATUS("You didn't specify an indexing method, so I won't try " + " to index anything.\n" + "If that isn't what you wanted, re-run with" + " --indexing=<methods>.\n"); indm = NULL; - cellr = CELLR_NONE; + } else { - if ( indm_str == NULL ) { - STATUS("You didn't specify an indexing method, so I " - " won't try to index anything.\n" - "If that isn't what you wanted, re-run with" - " --indexing=<method>.\n"); - indm = NULL; - indexer_needs_cell = 0; - } else { - indm = build_indexer_list(indm_str, - &indexer_needs_cell); - if ( indm == NULL ) { - ERROR("Invalid indexer list '%s'\n", indm_str); - return 1; - } - free(indm_str); - } - reduction_needs_cell = 0; - if ( scellr == NULL ) { - STATUS("You didn't specify a cell reduction method, so" - " I'm going to use 'reduce'.\n"); - cellr = CELLR_REDUCE; - reduction_needs_cell = 1; - } else { - int err; - cellr = parse_cell_reduction(scellr, &err, - &reduction_needs_cell); - if ( err ) { - ERROR("Unrecognised cell reduction '%s'\n", - scellr); - return 1; - } - free(scellr); + indm = build_indexer_list(indm_str); + if ( indm == NULL ) { + ERROR("Invalid indexer list '%s'\n", indm_str); + return 1; } + free(indm_str); } - /* No indexing -> no reduction */ - if ( indm == NULL ) reduction_needs_cell = 0; - if ( toler != NULL ) { int ttt; ttt = sscanf(toler, "%f,%f,%f,%f", @@ -586,7 +526,12 @@ int main(int argc, char *argv[]) cell = NULL; } - write_stream_header(ofh, argc, argv); + st = open_stream_for_write(outfile); + if ( st == NULL ) { + ERROR("Failed to open stream '%s'\n", outfile); + return 1; + } + free(outfile); /* Get first filename and use it to set up the indexing */ prepare_line = malloc(1024); @@ -609,7 +554,7 @@ int main(int argc, char *argv[]) /* Prepare the indexer */ if ( indm != NULL ) { ipriv = prepare_indexing(indm, cell, prepare_filename, - det, beam); + det, beam, tols); if ( ipriv == NULL ) { ERROR("Failed to prepare indexing.\n"); return 1; @@ -627,9 +572,7 @@ int main(int argc, char *argv[]) iargs.config_verbose = config_verbose; iargs.config_satcorr = config_satcorr; iargs.config_closer = config_closer; - iargs.config_insane = config_insane; iargs.config_bgsub = config_bgsub; - iargs.cellr = cellr; iargs.tols[0] = tols[0]; iargs.tols[1] = tols[1]; iargs.tols[2] = tols[2]; @@ -644,7 +587,6 @@ int main(int argc, char *argv[]) iargs.peaks = peaks; iargs.beam = beam; iargs.element = element; - iargs.stream_flags = stream_flags; iargs.hdf5_peak_path = hdf5_peak_path; iargs.copyme = copyme; iargs.ir_inn = ir_inn; @@ -654,9 +596,11 @@ int main(int argc, char *argv[]) iargs.integrate_saturated = integrate_saturated; iargs.no_revalidate = no_revalidate; iargs.integrate_found = integrate_found; + iargs.include_peaks = 1; + iargs.include_reflections = 1; /* FIXME! */ create_sandbox(&iargs, n_proc, prefix, config_basename, fh, - use_this_one_instead, ofh); + use_this_one_instead, st); free(prefix); |