From 29cca07716b48f9e433087f5dbb202165b1897e1 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Tue, 5 Feb 2013 11:36:44 +0100 Subject: WIP on bringing programs up to date --- libcrystfel/src/crystal.c | 77 ++++++++++++++++++++++++++++++++++++++++++++--- libcrystfel/src/crystal.h | 3 +- libcrystfel/src/dirax.c | 5 ++- libcrystfel/src/image.c | 18 +++++++++++ libcrystfel/src/mosflm.c | 2 +- libcrystfel/src/reax.c | 66 +++++++++++++++++++++++++++++----------- libcrystfel/src/stream.c | 64 +++++++++++---------------------------- libcrystfel/src/stream.h | 5 ++- 8 files changed, 165 insertions(+), 75 deletions(-) (limited to 'libcrystfel') 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; icandidate_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; in_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 */ -- cgit v1.2.3