diff options
-rw-r--r-- | libcrystfel/src/crystal.c | 77 | ||||
-rw-r--r-- | libcrystfel/src/crystal.h | 3 | ||||
-rw-r--r-- | libcrystfel/src/dirax.c | 5 | ||||
-rw-r--r-- | libcrystfel/src/image.c | 18 | ||||
-rw-r--r-- | libcrystfel/src/mosflm.c | 2 | ||||
-rw-r--r-- | libcrystfel/src/reax.c | 66 | ||||
-rw-r--r-- | libcrystfel/src/stream.c | 64 | ||||
-rw-r--r-- | libcrystfel/src/stream.h | 5 | ||||
-rw-r--r-- | src/process_hkl.c | 179 |
9 files changed, 256 insertions, 163 deletions
diff --git a/libcrystfel/src/crystal.c b/libcrystfel/src/crystal.c index 3c0649cc..8ba955f6 100644 --- a/libcrystfel/src/crystal.c +++ b/libcrystfel/src/crystal.c @@ -55,7 +55,7 @@ struct _crystal double osf; double profile_radius; int pr_dud; - double diffracting_resolution; + double resolution_limit; /* Integrated (or about-to-be-integrated) reflections */ RefList *reflections; @@ -63,7 +63,7 @@ struct _crystal }; -/************************** Setters and Constructors **************************/ +/************************** Constructor / Destructor **************************/ /** @@ -83,7 +83,7 @@ Crystal *crystal_new() cryst->cell = NULL; cryst->reflections = NULL; - cryst->diffracting_resolution = 0.0; + cryst->resolution_limit = 0.0; cryst->n_saturated = 0; return cryst; @@ -97,11 +97,78 @@ Crystal *crystal_new() * Frees a %Crystal, and all internal resources concerning that crystal. * */ -void crystal_free(UnitCell *cryst) +void crystal_free(Crystal *cryst) { if ( cryst == NULL ) return; - free(crysta); + if ( cryst->cell != NULL ) cell_free(cryst->cell); + if ( cryst->reflections != NULL ) reflist_free(cryst->reflections); + free(cryst); } /********************************** Getters ***********************************/ + + +UnitCell *crystal_get_cell(Crystal *cryst) +{ + return cryst->cell; +} + + +double crystal_get_profile_radius(Crystal *cryst) +{ + return cryst->profile_radius; +} + + +RefList *crystal_get_reflections(Crystal *cryst) +{ + return cryst->reflections; +} + + +double crystal_get_resolution_limit(Crystal *cryst) +{ + return cryst->resolution_limit; +} + + +long long int crystal_get_num_saturated_reflections(Crystal *cryst) +{ + return cryst->n_saturated; +} + + +/********************************** Setters ***********************************/ + + +void crystal_set_cell(Crystal *cryst, UnitCell *cell) +{ + if ( cryst->cell != NULL ) cell_free(cryst->cell); + cryst->cell = cell_new_from_cell(cell); +} + + +void crystal_set_profile_radius(Crystal *cryst, double r) +{ + cryst->profile_radius = r; +} + + +void crystal_set_reflections(Crystal *cryst, RefList *reflist) +{ + if ( cryst->reflections != NULL ) reflist_free(reflist); + cryst->reflections = reflist; +} + + +void crystal_set_resolution_limit(Crystal *cryst, double res) +{ + cryst->resolution_limit = res; +} + + +void crystal_set_num_saturated_reflections(Crystal *cryst, long long int n) +{ + cryst->n_saturated = n; +} diff --git a/libcrystfel/src/crystal.h b/libcrystfel/src/crystal.h index e5e40589..ca215399 100644 --- a/libcrystfel/src/crystal.h +++ b/libcrystfel/src/crystal.h @@ -60,7 +60,8 @@ extern void crystal_set_cell(Crystal *cryst, UnitCell *cell); extern void crystal_set_profile_radius(Crystal *cryst, double r); extern void crystal_set_reflections(Crystal *cryst, RefList *reflist); extern void crystal_set_resolution_limit(Crystal *cryst, double res); -extern void crystal_set_num_saturated_reflections(Crystal *cryst, int n); +extern void crystal_set_num_saturated_reflections(Crystal *cryst, + long long int n); #endif /* CRYSTAL_H */ diff --git a/libcrystfel/src/dirax.c b/libcrystfel/src/dirax.c index 0b2debda..6fbb9663 100644 --- a/libcrystfel/src/dirax.c +++ b/libcrystfel/src/dirax.c @@ -128,17 +128,20 @@ static int check_cell(struct dirax_private *dp, struct image *image, return 0; } - crystal_set_cell(cr, cell); + crystal_set_cell(cr, out); if ( dp->indm & INDEXING_CHECK_PEAKS ) { if ( !peak_sanity_check(image, &cr, 1) ) { crystal_free(cr); /* Frees the cell as well */ + cell_free(out); return 0; } } image_add_crystal(image, cr); + cell_free(out); /* Crystal makes its own copy */ + return 1; } diff --git a/libcrystfel/src/image.c b/libcrystfel/src/image.c index 093e20f2..6e4fd40a 100644 --- a/libcrystfel/src/image.c +++ b/libcrystfel/src/image.c @@ -159,3 +159,21 @@ void image_remove_feature(ImageFeatureList *flist, int idx) { flist->features[idx].valid = 0; } + + +void image_add_crystal(struct image *image, Crystal *cryst) +{ + Crystal **crs; + int n; + + n = image->n_crystals; + crs = realloc(image->crystals, (n+1)*sizeof(Crystal *)); + if ( crs == NULL ) { + ERROR("Failed to allocate memory for crystals.\n"); + return; + } + + crs[n] = cryst; + image->crystals = crs; + image->n_crystals = n+1; +} diff --git a/libcrystfel/src/mosflm.c b/libcrystfel/src/mosflm.c index 8b300012..d1d3c134 100644 --- a/libcrystfel/src/mosflm.c +++ b/libcrystfel/src/mosflm.c @@ -149,7 +149,7 @@ static int check_cell(struct mosflm_private *mp, struct image *image, return 0; } - crystal_set_cell(cr, cell); + crystal_set_cell(cr, out); if ( mp->indm & INDEXING_CHECK_PEAKS ) { if ( !peak_sanity_check(image, &cr, 1) ) { diff --git a/libcrystfel/src/reax.c b/libcrystfel/src/reax.c index 7529d12f..b144e266 100644 --- a/libcrystfel/src/reax.c +++ b/libcrystfel/src/reax.c @@ -103,7 +103,7 @@ struct reax_search struct reax_private { - IndexingPrivate base; + IndexingMethod indm; struct dvec *directions; int n_dir; double angular_inc; @@ -915,9 +915,39 @@ static void add_cell_candidate(struct cell_candidate_list *cl, UnitCell *cnew, } +static int check_cell(struct reax_private *dp, struct image *image, + UnitCell *cell) +{ + UnitCell *out; + Crystal *cr; + + out = cell_new_from_cell(cell); + + cr = crystal_new(); + if ( cr == NULL ) { + ERROR("Failed to allocate crystal.\n"); + return 0; + } + + crystal_set_cell(cr, out); + + if ( dp->indm & INDEXING_CHECK_PEAKS ) { + if ( !peak_sanity_check(image, &cr, 1) ) { + crystal_free(cr); /* Frees the cell as well */ + return 0; + } + } + + image_add_crystal(image, cr); + + return 1; +} + + static void assemble_cells_from_candidates(struct image *image, struct reax_search *s, - UnitCell *cell) + UnitCell *cell, + struct reax_private *p) { int i, j, k; signed int ti, tj, tk; @@ -967,7 +997,9 @@ static void assemble_cells_from_candidates(struct image *image, continue; } - peak_lattice_agreement(image, cnew, &fom); + /* FIXME! */ + //peak_lattice_agreement(image, cnew, &fom); + fom = 1.0; add_cell_candidate(&cl, cnew, fom); } @@ -985,22 +1017,20 @@ static void assemble_cells_from_candidates(struct image *image, cell_get_parameters(cl.cand[i].cell, &a, &b, &c, &al, &be, &ga); cell_get_parameters(cl.cand[i].cell, &aA, &bA, &cA, &alA, &beA, &gaA); - if ( (a - aA) > aA/10.0 ) w = 1; - if ( (b - bA) > bA/10.0 ) w = 1; - if ( (c - cA) > cA/10.0 ) w = 1; - if ( (al - alA) > deg2rad(5.0) ) w = 1; - if ( (be - beA) > deg2rad(5.0) ) w = 1; - if ( (ga - gaA) > deg2rad(5.0) ) w = 1; + if ( fabs(a - aA) > aA/10.0 ) w = 1; + if ( fabs(b - bA) > bA/10.0 ) w = 1; + if ( fabs(c - cA) > cA/10.0 ) w = 1; + if ( fabs(al - alA) > deg2rad(5.0) ) w = 1; + if ( fabs(be - beA) > deg2rad(5.0) ) w = 1; + if ( fabs(ga - gaA) > deg2rad(5.0) ) w = 1; if ( w ) { STATUS("This cell is a long way from that sought:\n"); cell_print(cl.cand[i].cell); } } - image->ncells = cl.n_cand; - assert(image->ncells <= MAX_CELL_CANDIDATES); for ( i=0; i<cl.n_cand; i++ ) { - image->candidate_cells[i] = cl.cand[i].cell; + if ( check_cell(p, image, cl.cand[i].cell) ) break; } free(cl.cand); @@ -1016,7 +1046,6 @@ void reax_index(IndexingPrivate *pp, struct image *image, UnitCell *cell) struct reax_search *s; int i; - assert(pp->indm == INDEXING_REAX); p = (struct reax_private *)pp; fft_in = fftw_malloc(p->nel*sizeof(double)); @@ -1039,7 +1068,7 @@ void reax_index(IndexingPrivate *pp, struct image *image, UnitCell *cell) // fft_in, fft_out, p->plan, smin, smax, // image->det, p); - assemble_cells_from_candidates(image, s, cell); + assemble_cells_from_candidates(image, s, cell, p); for ( i=0; i<s->n_search; i++ ) { free(s->search[i].cand); @@ -1051,7 +1080,9 @@ void reax_index(IndexingPrivate *pp, struct image *image, UnitCell *cell) } -IndexingPrivate *reax_prepare() +IndexingPrivate *reax_prepare(IndexingMethod indm, UnitCell *cell, + const char *filename, struct detector *det, + struct beam_params *beam, float *ltl) { struct reax_private *p; int samp; @@ -1060,8 +1091,6 @@ IndexingPrivate *reax_prepare() p = calloc(1, sizeof(*p)); if ( p == NULL ) return NULL; - p->base.indm = INDEXING_REAX; - p->angular_inc = deg2rad(1.0); /* Reserve memory, over-estimating the number of directions */ @@ -1122,6 +1151,8 @@ IndexingPrivate *reax_prepare() p->r_plan = fftw_plan_dft_2d(p->cw, p->ch, p->r_fft_in, p->r_fft_out, 1, FFTW_MEASURE); + p->indm = indm; + return (IndexingPrivate *)p; } @@ -1130,7 +1161,6 @@ void reax_cleanup(IndexingPrivate *pp) { struct reax_private *p; - assert(pp->indm == INDEXING_REAX); p = (struct reax_private *)pp; free(p->directions); diff --git a/libcrystfel/src/stream.c b/libcrystfel/src/stream.c index 943ee656..31fb8527 100644 --- a/libcrystfel/src/stream.c +++ b/libcrystfel/src/stream.c @@ -69,30 +69,6 @@ struct _stream }; -int count_patterns(FILE *fh) -{ - char *rval; - - int n_total_patterns = 0; - do { - char line[1024]; - - rval = fgets(line, 1023, fh); - if ( rval == NULL ) continue; - chomp(line); - if ( strcmp(line, CHUNK_END_MARKER) == 0 ) n_total_patterns++; - - } while ( rval != NULL ); - - if ( ferror(fh) ) { - ERROR("Read error while counting patterns.\n"); - return 0; - } - - return n_total_patterns; -} - - static int read_peaks(FILE *fh, struct image *image) { char *rval = NULL; @@ -482,28 +458,6 @@ void write_stream_header(FILE *ofh, int argc, char *argv[]) } -int skip_some_files(FILE *fh, int n) -{ - char *rval = NULL; - int n_patterns = 0; - - do { - - char line[1024]; - - if ( n_patterns == n ) return 0; - - rval = fgets(line, 1023, fh); - if ( rval == NULL ) continue; - chomp(line); - if ( strcmp(line, CHUNK_END_MARKER) == 0 ) n_patterns++; - - } while ( rval != NULL ); - - return 1; -} - - Stream *open_stream_for_read(const char *filename) { Stream *st; @@ -591,3 +545,21 @@ int is_stream(const char *filename) return 0; } + + +/** + * rewind_stream: + * @st: A %Stream + * + * Attempts to set the file pointer for @st to the start of the stream, so that + * later calls to read_chunk() will repeat the sequence of chunks from the + * start. + * + * Programs must not assume that this operation always succeeds! + * + * Returns: non-zero if the stream could not be rewound. + */ +int rewind_stream(Stream *st) +{ + return fseek(st->fh, 0, SEEK_SET); +} diff --git a/libcrystfel/src/stream.h b/libcrystfel/src/stream.h index db74583d..9340ca0d 100644 --- a/libcrystfel/src/stream.h +++ b/libcrystfel/src/stream.h @@ -49,9 +49,8 @@ extern void write_chunk(Stream *st, struct image *image, struct hdfile *hdfile, extern int read_chunk(Stream *st, struct image *image); -/* Nasty functions that should be avoided */ -extern int count_patterns(FILE *fh); -extern int skip_some_files(FILE *fh, int n); +extern int rewind_stream(Stream *st); + extern int is_stream(const char *filename); #endif /* STREAM_H */ diff --git a/src/process_hkl.c b/src/process_hkl.c index 8bfbea2b..e9a8daa5 100644 --- a/src/process_hkl.c +++ b/src/process_hkl.c @@ -3,12 +3,12 @@ * * Assemble and process FEL Bragg intensities * - * 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 Lorenzo Galli * * Authors: - * 2009-2012 Thomas White <taw@physics.org> + * 2009-2013 Thomas White <taw@physics.org> * 2011 Andrew Martin <andrew.martin@desy.de> * 2012 Lorenzo Galli <lorenzo.galli@desy.de> * @@ -48,6 +48,8 @@ #include "stream.h" #include "reflist.h" #include "image.h" +#include "crystal.h" +#include "thread-pool.h" static void show_help(const char *s) @@ -170,11 +172,11 @@ static double scale_intensities(RefList *reference, RefList *new, } -static int merge_pattern(RefList *model, struct image *new, RefList *reference, - const SymOpList *sym, - double *hist_vals, signed int hist_h, - signed int hist_k, signed int hist_l, int *hist_n, - int config_nopolar) +static int merge_crystal(RefList *model, struct image *image, Crystal *cr, + RefList *reference, const SymOpList *sym, + double *hist_vals, signed int hist_h, + signed int hist_k, signed int hist_l, int *hist_n, + int config_nopolar) { Reflection *refl; RefListIterator *iter; @@ -184,17 +186,18 @@ static int merge_pattern(RefList *model, struct image *new, RefList *reference, double csx, csy, csz; if ( reference != NULL ) { - scale = scale_intensities(reference, new->reflections, sym); + scale = scale_intensities(reference, + crystal_get_reflections(cr), sym); } else { scale = 1.0; } if ( isnan(scale) ) return 1; - cell_get_reciprocal(new->indexed_cell, &asx, &asy, &asz, - &bsx, &bsy, &bsz, - &csx, &csy, &csz); + cell_get_reciprocal(crystal_get_cell(cr), &asx, &asy, &asz, + &bsx, &bsy, &bsz, + &csx, &csy, &csz); - for ( refl = first_refl(new->reflections, &iter); + for ( refl = first_refl(crystal_get_reflections(cr), &iter); refl != NULL; refl = next_refl(refl, iter) ) { @@ -226,7 +229,7 @@ static int merge_pattern(RefList *model, struct image *new, RefList *reference, yl = h*asy + k*bsy + l*csy; zl = h*asz + k*bsz + l*csz; - ool = 1.0 / new->lambda; + ool = 1.0 / image->lambda; tt = angle_between(0.0, 0.0, 1.0, xl, yl, zl+ool); phi = atan2(yl, xl); pa = pow(sin(phi)*sin(tt), 2.0); @@ -260,63 +263,72 @@ static int merge_pattern(RefList *model, struct image *new, RefList *reference, } -static void merge_all(FILE *fh, RefList *model, RefList *reference, - int config_startafter, int config_stopafter, - const SymOpList *sym, - int n_total_patterns, - double *hist_vals, signed int hist_h, - signed int hist_k, signed int hist_l, - int *hist_i, int config_nopolar, int min_measurements) +static void display_progress(int n_images, int n_crystals, int n_crystals_used) +{ + if ( !isatty(STDERR_FILENO) ) return; + if ( tcgetpgrp(STDERR_FILENO) != getpgrp() ) return; + + pthread_mutex_lock(&stderr_lock); + fprintf(stderr, "\r%i images processed, %i crystals, %i crystals used", + n_images, n_crystals, n_crystals_used); + pthread_mutex_unlock(&stderr_lock); + + fflush(stdout); +} + + + +static int merge_all(Stream *st, RefList *model, RefList *reference, + const SymOpList *sym, + double *hist_vals, signed int hist_h, + signed int hist_k, signed int hist_l, + int *hist_i, int config_nopolar, int min_measurements) { int rval; - int n_patterns = 0; - int n_used = 0; + int n_images = 0; + int n_crystals = 0; + int n_crystals_used = 0; Reflection *refl; RefListIterator *iter; - if ( skip_some_files(fh, config_startafter) ) { - ERROR("Failed to skip first %i files.\n", config_startafter); - return; - } - do { struct image image; + int i; image.det = NULL; /* Get data from next chunk */ - rval = read_chunk(fh, &image); + rval = read_chunk(st, &image); if ( rval ) break; - n_patterns++; + n_images++; - if ( (image.reflections != NULL) && (image.indexed_cell) ) { + for ( i=0; i<image.n_crystals; i++ ) { int r; + Crystal *cr = image.crystals[i]; - r = merge_pattern(model, &image, reference, sym, + n_crystals++; + r = merge_crystal(model, &image, cr, reference, sym, hist_vals, hist_h, hist_k, hist_l, hist_i, config_nopolar); - if ( r == 0 ) n_used++; + if ( r == 0 ) n_crystals_used++; + + crystal_free(cr); } free(image.filename); - reflist_free(image.reflections); image_feature_list_free(image.features); - cell_free(image.indexed_cell); - - progress_bar(n_patterns, n_total_patterns-config_startafter, - "Merging"); - if ( config_stopafter ) { - if ( n_patterns == config_stopafter ) break; - } + display_progress(n_images, n_crystals, n_crystals_used); } while ( rval == 0 ); + if ( rval ) return 1; + for ( refl = first_refl(model, &iter); refl != NULL; refl = next_refl(refl, iter) ) @@ -339,7 +351,7 @@ static void merge_all(FILE *fh, RefList *model, RefList *reference, } - STATUS("%i of the patterns could be used.\n", n_used); + return 0; } @@ -348,14 +360,11 @@ int main(int argc, char *argv[]) int c; char *filename = NULL; char *output = NULL; - FILE *fh; + Stream *st; RefList *model; int config_maxonly = 0; - int config_startafter = 0; - int config_stopafter = 0; int config_sum = 0; int config_scale = 0; - unsigned int n_total_patterns; char *sym_str = NULL; SymOpList *sym; char *histo = NULL; @@ -369,6 +378,7 @@ int main(int argc, char *argv[]) int config_nopolar = 0; char *rval; int min_measurements = 2; + int r; /* Long options */ const struct option longopts[] = { @@ -409,11 +419,11 @@ int main(int argc, char *argv[]) break; case 's' : - config_stopafter = atoi(optarg); + ERROR("The option '--stop-after' no longer works.\n"); break; case 'f' : - config_startafter = atoi(optarg); + ERROR("The option '--start-after' no longer works.\n"); break; case 'y' : @@ -465,25 +475,11 @@ int main(int argc, char *argv[]) free(sym_str); /* Open the data stream */ - if ( strcmp(filename, "-") == 0 ) { - fh = stdin; - } else { - fh = fopen(filename, "r"); - } - free(filename); - if ( fh == NULL ) { - ERROR("Failed to open input file\n"); - return 1; - } - - /* Count the number of patterns in the file */ - n_total_patterns = count_patterns(fh); - if ( n_total_patterns == 0 ) { - ERROR("No patterns to process.\n"); + st = open_stream_for_read(filename); + if ( st == NULL ) { + ERROR("Failed to open stream.\n"); return 1; } - STATUS("There are %i patterns to process\n", n_total_patterns); - rewind(fh); model = reflist_new(); @@ -497,7 +493,8 @@ int main(int argc, char *argv[]) return 1; } - space_for_hist = n_total_patterns * num_equivs(sym, NULL); + /* FIXME: This array must grow as necessary. */ + space_for_hist = 0 * num_equivs(sym, NULL); hist_vals = malloc(space_for_hist * sizeof(double)); free(histo); STATUS("Histogramming %i %i %i -> ", hist_h, hist_k, hist_l); @@ -529,40 +526,46 @@ int main(int argc, char *argv[]) } hist_i = 0; - merge_all(fh, model, NULL, config_startafter, config_stopafter, - sym, n_total_patterns, hist_vals, hist_h, hist_k, hist_l, - &hist_i, config_nopolar, min_measurements); - if ( ferror(fh) ) { - ERROR("Stream read error.\n"); + r = merge_all(st, model, NULL, sym, hist_vals, hist_h, hist_k, hist_l, + &hist_i, config_nopolar, min_measurements); + fprintf(stderr, "\n"); + if ( r ) { + ERROR("Error while reading stream.\n"); return 1; } - rewind(fh); if ( config_scale ) { RefList *reference; - STATUS("Extra pass for scaling...\n"); + if ( rewind_stream(st) ) { - reference = copy_reflist(model); + ERROR("Couldn't rewind stream - scaling cannot be " + "performed.\n"); - reflist_free(model); - model = reflist_new(); + } else { - rewind(fh); + int r; - merge_all(fh, model, reference, - config_startafter, config_stopafter, sym, - n_total_patterns, - hist_vals, hist_h, hist_k, hist_l, &hist_i, - config_nopolar, min_measurements); + STATUS("Extra pass for scaling...\n"); - if ( ferror(fh) ) { - ERROR("Stream read error.\n"); - return 1; - } + reference = copy_reflist(model); + + reflist_free(model); + model = reflist_new(); + + r = merge_all(st, model, reference, sym, + hist_vals, hist_h, hist_k, hist_l, &hist_i, + config_nopolar, min_measurements); + fprintf(stderr, "\n"); + if ( r ) { + ERROR("Error while reading stream.\n"); + return 1; + } + + reflist_free(reference); - reflist_free(reference); + } } @@ -579,7 +582,7 @@ int main(int argc, char *argv[]) write_reflist(output, model); - fclose(fh); + close_stream(st); free(sym); reflist_free(model); |