From 38ab00ccb6869318f1a12a96443fa6c23de1f6e7 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 11 Mar 2011 17:51:03 +0100 Subject: First part of stream rework --- src/hdf5-file.c | 2 +- src/image.h | 2 +- src/index.c | 2 + src/indexamajig.c | 91 +++++++++++++++--------------------- src/pattern_sim.c | 20 +++++--- src/peaks.c | 134 ++++++++---------------------------------------------- src/peaks.h | 13 ++---- src/reintegrate.c | 11 ++--- src/stream.c | 79 ++++++++++++++++++++++++++++++++ src/stream.h | 3 ++ 10 files changed, 165 insertions(+), 192 deletions(-) (limited to 'src') diff --git a/src/hdf5-file.c b/src/hdf5-file.c index 1b957b4c..20141522 100644 --- a/src/hdf5-file.c +++ b/src/hdf5-file.c @@ -3,7 +3,7 @@ * * Read/write HDF5 data files * - * (c) 2006-2010 Thomas White + * (c) 2006-2011 Thomas White * * Part of CrystFEL - crystallography with a FEL * diff --git a/src/image.h b/src/image.h index e0c190f2..69c36de8 100644 --- a/src/image.h +++ b/src/image.h @@ -85,7 +85,7 @@ struct image { int width; int height; - /* Reflections (used for scaling ONLY) */ + /* Integrated (or about-to-be-integrated) reflections */ RefList *reflections; /* Detected peaks */ diff --git a/src/index.c b/src/index.c index 554c099d..ddb25eb4 100644 --- a/src/index.c +++ b/src/index.c @@ -128,6 +128,8 @@ void index_pattern(struct image *image, UnitCell *cell, IndexingMethod *indm, int i; int n = 0; + if ( indm == NULL ) return; + map_all_peaks(image); while ( indm[n] != INDEXING_NONE ) { diff --git a/src/indexamajig.c b/src/indexamajig.c index 51461c91..72491450 100644 --- a/src/indexamajig.c +++ b/src/indexamajig.c @@ -36,6 +36,7 @@ #include "thread-pool.h" #include "beam-parameters.h" #include "geometry.h" +#include "stream.h" enum { @@ -50,9 +51,8 @@ struct static_index_args UnitCell *cell; int config_cmfilter; int config_noisefilter; - int config_dumpfound; int config_verbose; - int config_nearbragg; + int stream_flags; int config_polar; int config_satcorr; int config_closer; @@ -109,7 +109,8 @@ static void show_help(const char *s) "\n" " -i, --input= Specify file containing list of images to process.\n" " '-' means stdin, which is the default.\n" -" -o, --output= Write indexed stream to this file. '-' for stdout.\n" +" -o, --output= Write output stream to this file. '-' for stdout.\n" +" Default: indexamajig.stream\n" "\n" " --indexing= Use 'methods' for indexing. Provide one or more\n" " methods separated by commas. Choose from:\n" @@ -129,25 +130,19 @@ static void show_help(const char *s) " This is the default method.\n" " hdf5 : Get from /processing/hitfinder/peakinfo\n" " in the HDF5 file.\n" -"\n" -"\nWith just the above options, this program does not do much of practical use." -"\nYou should also enable some of the following:\n\n" -" --near-bragg Output a list of reflection intensities to stdout.\n" -" When pixels with fractional indices within 0.1 of\n" -" integer values (the Bragg condition) are found,\n" -" the integral of pixels within a ten pixel radius\n" -" of the nearest-to-Bragg pixel will be reported as\n" -" the intensity. The centroid of the pixels will\n" -" be given as the coordinates, as well as the h,k,l\n" -" (integer) indices of the reflection. If a peak\n" -" was located by the initial peak search close to\n" -" the \"near Bragg\" location, its coordinates will\n" -" be taken as the centre instead.\n" -" --dump-peaks Write the results of the peak search to stdout.\n" -" The intensities in this list are from the\n" -" centroid/integration procedure.\n" -"\n" -"\nFor more control over the process, you might need:\n\n" +"\n\n" +"You can control what information is included in the output stream using\n" +"' --record='. Possible flags are:\n" +"pixels Include a list of sums of pixel values within the\n" +" integration domain, correcting for individual pixel\n" +" solid angles.\n" +"integrated Include a list of reflection intensities, produced by\n" +" integrating around predicted peak locations.\n" +" The flags 'pixels' and 'integrated' are mutually exclusive.\n" +"peaks Include peak locations and intensities from the peak search.\n" +"peaksifindexed Include peak locations only if the pattern could be indexed.\n" +"\n\n" +"For more control over the process, you might need:\n\n" " --cell-reduction= Use as the cell reduction method. Choose from:\n" " none : no matching, just use the raw cell.\n" " reduce : full cell reduction.\n" @@ -196,9 +191,7 @@ static void process_image(void *pp, int cookie) UnitCell *cell = pargs->static_args.cell; int config_cmfilter = pargs->static_args.config_cmfilter; int config_noisefilter = pargs->static_args.config_noisefilter; - int config_dumpfound = pargs->static_args.config_dumpfound; int config_verbose = pargs->static_args.config_verbose; - int config_nearbragg = pargs->static_args.config_nearbragg; int config_polar = pargs->static_args.config_polar; IndexingMethod *indm = pargs->static_args.indm; const struct beam_params *beam = pargs->static_args.beam; @@ -274,7 +267,6 @@ static void process_image(void *pp, int cookie) /* Get peaks from HDF5 */ if ( get_peaks(&image, hdfile) ) { ERROR("Failed to get peaks from HDF5 file.\n"); - return; } break; case PEAK_ZAEF : @@ -288,51 +280,44 @@ static void process_image(void *pp, int cookie) free(image.data); image.data = data_for_measurement; - if ( config_dumpfound ) { - dump_peaks(&image, pargs->static_args.ofh, - pargs->static_args.output_mutex); - } - - /* Not indexing? Then there's nothing left to do. */ - if ( indm == NULL ) goto done; - /* Calculate orientation matrix (by magic) */ index_pattern(&image, cell, indm, pargs->static_args.cellr, config_verbose, pargs->static_args.ipriv, pargs->static_args.config_insane); /* No cell at this point? Then we're done. */ - if ( image.indexed_cell == NULL ) goto done; - pargs->indexable = 1; + if ( image.indexed_cell != NULL ) pargs->indexable = 1; /* Measure intensities if requested */ - if ( config_nearbragg ) { - RefList *reflections; + /* Do EITHER: */ - image.div = beam->divergence; - image.bw = beam->bandwidth; - image.profile_radius = 0.0001e9; + //image.div = beam->divergence; + //image.bw = beam->bandwidth; + //image.profile_radius = 0.0001e9; + //image.reflections = find_intersections(&image, + // image.indexed_cell, 0); - //reflections = find_intersections(&image, image.indexed_cell, - // 0); - reflections = find_projected_peaks(&image, image.indexed_cell, - 0, 0.1); + image.reflections = find_projected_peaks(&image, + image.indexed_cell, + 0, 0.1); - output_intensities(&image, image.indexed_cell, reflections, - pargs->static_args.output_mutex, - config_polar, - pargs->static_args.config_closer, - pargs->static_args.ofh); + integrate_reflections(&image, config_polar, + pargs->static_args.config_closer); - reflist_free(reflections); + /* OR */ - } + //image.reflections = integrate_pixels(&image, 0, 0.1, + // config_polar); + + pthread_mutex_lock(pargs->static_args.output_mutex); + write_chunk(pargs->static_args.ofh, &image, + pargs->static_args.stream_flags); + pthread_mutex_unlock(pargs->static_args.output_mutex); /* Only free cell if found */ cell_free(image.indexed_cell); -done: free(image.data); free(image.flags); image_feature_list_free(image.features); @@ -737,9 +722,7 @@ int main(int argc, char *argv[]) qargs.static_args.cell = cell; qargs.static_args.config_cmfilter = config_cmfilter; qargs.static_args.config_noisefilter = config_noisefilter; - qargs.static_args.config_dumpfound = config_dumpfound; qargs.static_args.config_verbose = config_verbose; - qargs.static_args.config_nearbragg = config_nearbragg; qargs.static_args.config_polar = config_polar; qargs.static_args.config_satcorr = config_satcorr; qargs.static_args.config_closer = config_closer; diff --git a/src/pattern_sim.c b/src/pattern_sim.c index c823c3a3..948458bc 100644 --- a/src/pattern_sim.c +++ b/src/pattern_sim.c @@ -551,15 +551,23 @@ int main(int argc, char *argv[]) if ( config_nearbragg ) { - RefList *reflections; + /* Do EITHER: */ - reflections = find_projected_peaks(&image, cell, - 0, 0.1); + //image.div = beam->divergence; + //image.bw = beam->bandwidth; + //image.profile_radius = 0.0001e9; + //image.reflections = find_intersections(&image, + // image.indexed_cell, 0); - output_intensities(&image, cell, reflections, NULL, - 0, 0, stdout); + image.reflections = find_projected_peaks(&image, + image.indexed_cell, 0, 0.1); - reflist_free(reflections); + integrate_reflections(&image, 0, 0); + + /* OR */ + + //image.reflections = integrate_pixels(&image, 0, 0.1, + // config_polar); } diff --git a/src/peaks.c b/src/peaks.c index 4bf56ce9..69afc20f 100644 --- a/src/peaks.c +++ b/src/peaks.c @@ -416,39 +416,6 @@ void search_peaks(struct image *image, float threshold, float min_gradient) } -void dump_peaks(struct image *image, FILE *ofh, pthread_mutex_t *mutex) -{ - int i; - - /* Get exclusive access to the output stream if necessary */ - if ( mutex != NULL ) pthread_mutex_lock(mutex); - - fprintf(ofh, "Peaks from peak search in %s\n", image->filename); - fprintf(ofh, " x/px y/px (1/d)/nm^-1 Intensity\n"); - - for ( i=0; ifeatures); i++ ) { - - struct imagefeature *f; - struct rvec r; - double q; - - f = image_get_feature(image->features, i); - if ( f == NULL ) continue; - - r = get_q(image, f->x, f->y, NULL, 1.0/image->lambda); - q = modulus(r.u, r.v, r.w); - - fprintf(ofh, "%8.3f %8.3f %8.3f %12.3f\n", - f->x, f->y, q/1.0e9, f->intensity); - - } - - fprintf(ofh, "\n"); - - if ( mutex != NULL ) pthread_mutex_unlock(mutex); -} - - RefList *find_projected_peaks(struct image *image, UnitCell *cell, int circular_domain, double domain_r) { @@ -600,63 +567,13 @@ int peak_sanity_check(struct image *image, UnitCell *cell, } -static void output_header(FILE *ofh, UnitCell *cell, struct image *image) -{ - double asx, asy, asz; - double bsx, bsy, bsz; - double csx, csy, csz; - double a, b, c, al, be, ga; - - fprintf(ofh, "Reflections from indexing in %s\n", image->filename); - - cell_get_parameters(cell, &a, &b, &c, &al, &be, &ga); - fprintf(ofh, "Cell parameters %7.5f %7.5f %7.5f nm, %7.5f %7.5f %7.5f deg\n", - a*1.0e9, b*1.0e9, c*1.0e9, - rad2deg(al), rad2deg(be), rad2deg(ga)); - - cell_get_reciprocal(cell, &asx, &asy, &asz, - &bsx, &bsy, &bsz, - &csx, &csy, &csz); - fprintf(ofh, "astar = %+9.7f %+9.7f %+9.7f nm^-1\n", - asx/1e9, asy/1e9, asz/1e9); - fprintf(ofh, "bstar = %+9.7f %+9.7f %+9.7f nm^-1\n", - bsx/1e9, bsy/1e9, bsz/1e9); - fprintf(ofh, "cstar = %+9.7f %+9.7f %+9.7f nm^-1\n", - csx/1e9, csy/1e9, csz/1e9); - - if ( image->f0_available ) { - fprintf(ofh, "f0 = %7.5f (arbitrary gas detector units)\n", - image->f0); - } else { - fprintf(ofh, "f0 = invalid\n"); - } - - fprintf(ofh, "photon_energy_eV = %f\n", - J_to_eV(ph_lambda_to_en(image->lambda))); - -} - - -void output_intensities(struct image *image, UnitCell *cell, - RefList *reflections, pthread_mutex_t *mutex, int polar, - int use_closer, FILE *ofh) +/* Integrate the list of predicted reflections in "image" */ +void integrate_reflections(struct image *image, int polar, int use_closer) { - double asx, asy, asz; - double bsx, bsy, bsz; - double csx, csy, csz; Reflection *refl; RefListIterator *iter; - /* Get exclusive access to the output stream if necessary */ - if ( mutex != NULL ) pthread_mutex_lock(mutex); - - output_header(ofh, cell, image); - - cell_get_reciprocal(cell, &asx, &asy, &asz, - &bsx, &bsy, &bsz, - &csx, &csy, &csz); - - for ( refl = first_refl(reflections, &iter); + for ( refl = first_refl(image->reflections, &iter); refl != NULL; refl = next_refl(refl, iter) ) { @@ -666,7 +583,6 @@ void output_intensities(struct image *image, UnitCell *cell, double bg, max; struct panel *p; double px, py; - signed int h, k, l; get_detector_pos(refl, &px, &py); p = find_panel(image->det, px, py); @@ -734,23 +650,14 @@ void output_intensities(struct image *image, UnitCell *cell, } - /* Write h,k,l, integrated intensity and centroid coordinates */ - get_indices(refl, &h, &k, &l); - fprintf(ofh, "%3i %3i %3i %6f (at %5.2f,%5.2f) max=%6f bg=%6f\n", - h, k, l, intensity, x, y, max, bg); + set_int(refl, intensity); } - - /* Blank line at end */ - fprintf(ofh, "\n"); - - if ( mutex != NULL ) pthread_mutex_unlock(mutex); } -void output_pixels(struct image *image, UnitCell *cell, - pthread_mutex_t *mutex, int do_polar, - FILE *ofh, int circular_domain, double domain_r) +RefList *integrate_pixels(struct image *image, int circular_domain, + double domain_r, int do_polar) { int i; double ax, ay, az; @@ -762,26 +669,25 @@ void output_pixels(struct image *image, UnitCell *cell, double *xmom; double *ymom; ReflItemList *obs; - - /* Get exclusive access to the output stream if necessary */ - if ( mutex != NULL ) pthread_mutex_lock(mutex); - - output_header(ofh, cell, image); + RefList *reflections; obs = new_items(); intensities = new_list_intensity(); xmom = new_list_intensity(); ymom = new_list_intensity(); + reflections = reflist_new(); /* "Borrow" direction values to get reciprocal lengths */ - cell_get_reciprocal(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz); + cell_get_reciprocal(image->indexed_cell, &ax, &ay, &az, + &bx, &by, &bz, + &cx, &cy, &cz); aslen = modulus(ax, ay, az); bslen = modulus(bx, by, bz); cslen = modulus(cx, cy, cz); - cell_get_cartesian(cell, &ax, &ay, &az, - &bx, &by, &bz, - &cx, &cy, &cz); + cell_get_cartesian(image->indexed_cell, &ax, &ay, &az, + &bx, &by, &bz, + &cx, &cy, &cz); /* For each pixel */ fesetround(1); /* Round towards nearest */ for ( fs=0; fswidth; fs++ ) { @@ -890,6 +796,7 @@ void output_pixels(struct image *image, UnitCell *cell, struct refl_item *it; double intensity, xmomv, ymomv; double xp, yp; + Reflection *refl; it = get_item(obs, i); intensity = lookup_intensity(intensities, it->h, it->k, it->l); @@ -899,19 +806,16 @@ void output_pixels(struct image *image, UnitCell *cell, xp = xmomv / (double)intensity; yp = ymomv / (double)intensity; - fprintf(ofh, "%3i %3i %3i %6f (at %5.2f,%5.2f)\n", - it->h, it->k, it->l, intensity, xp, yp); + refl = add_refl(reflections, it->h, it->k, it->l); + set_int(refl, intensity); + set_detector_pos(refl, 0.0, xp, yp); } - fprintf(ofh, "No peak statistics, because output_pixels() was used.\n"); - /* Blank line at end */ - fprintf(ofh, "\n"); - free(xmom); free(ymom); free(intensities); delete_items(obs); - if ( mutex != NULL ) pthread_mutex_unlock(mutex); + return reflections; } diff --git a/src/peaks.h b/src/peaks.h index 2acb249b..7a6557f1 100644 --- a/src/peaks.h +++ b/src/peaks.h @@ -23,21 +23,16 @@ extern void search_peaks(struct image *image, float threshold, float min_gradient); -extern void dump_peaks(struct image *image, FILE *ofh, pthread_mutex_t *mutex); -extern void output_intensities(struct image *image, UnitCell *cell, - RefList *reflections, - pthread_mutex_t *mutex, int polar, - int use_closer, FILE *ofh); - -extern void output_pixels(struct image *image, UnitCell *cell, - pthread_mutex_t *mutex, int do_polar, - FILE *ofh, int circular_domain, double domain_r); +extern void integrate_reflections(struct image *image, + int polar, int use_closer); extern int peak_sanity_check(struct image *image, UnitCell *cell, int circular_domain, double domain_r); + extern RefList *find_projected_peaks(struct image *image, UnitCell *cell, int circular_domain, double domain_r); + extern int integrate_peak(struct image *image, int xp, int yp, double *xc, double *yc, double *intensity, double *pbg, double *pmax, diff --git a/src/reintegrate.c b/src/reintegrate.c index 67e7b347..fcffd499 100644 --- a/src/reintegrate.c +++ b/src/reintegrate.c @@ -140,15 +140,14 @@ static void process_image(void *pg, int cookie) reflections = find_projected_peaks(&image, image.indexed_cell, 0, 0.1); - output_intensities(&image, image.indexed_cell, reflections, - pargs->output_mutex, - pargs->config_polar, - pargs->config_closer, - pargs->ofh); - reflist_free(reflections); } + pthread_mutex_lock(pargs->output_mutex); + write_chunk(pargs->ofh, &image, pargs->stream_flags); + pthread_mutex_unlock(pargs->output_mutex); + + free(image.data); if ( image.flags != NULL ) free(image.flags); hdfile_close(hdfile); diff --git a/src/stream.c b/src/stream.c index ffee4a6c..fbe0abad 100644 --- a/src/stream.c +++ b/src/stream.c @@ -20,6 +20,7 @@ #include "cell.h" #include "utils.h" +#include "image.h" int count_patterns(FILE *fh) @@ -73,6 +74,84 @@ static UnitCell *read_orientation_matrix(FILE *fh) } +static void write_reflections(struct image *image, FILE *ofh) +{ +} + + +static void write_peaks(struct image *image, FILE *ofh) +{ + int i; + + fprintf(ofh, "Peaks from peak search in %s\n", image->filename); + fprintf(ofh, " x/px y/px (1/d)/nm^-1 Intensity\n"); + + for ( i=0; ifeatures); i++ ) { + + struct imagefeature *f; + struct rvec r; + double q; + + f = image_get_feature(image->features, i); + if ( f == NULL ) continue; + + r = get_q(image, f->x, f->y, NULL, 1.0/image->lambda); + q = modulus(r.u, r.v, r.w); + + fprintf(ofh, "%8.3f %8.3f %8.3f %12.3f\n", + f->x, f->y, q/1.0e9, f->intensity); + + } + + fprintf(ofh, "\n"); +} + + + +void write_chunk(FILE *ofh, struct image *image, int flags) +{ + double asx, asy, asz; + double bsx, bsy, bsz; + double csx, csy, csz; + double a, b, c, al, be, ga; + + fprintf(ofh, "----- Begin chunk -----\n"); + + fprintf(ofh, "Image filename: %s\n", image->filename); + + cell_get_parameters(image->indexed_cell, &a, &b, &c, &al, &be, &ga); + fprintf(ofh, "Cell parameters %7.5f %7.5f %7.5f nm," + " %7.5f %7.5f %7.5f deg\n", + a*1.0e9, b*1.0e9, c*1.0e9, + rad2deg(al), rad2deg(be), rad2deg(ga)); + + cell_get_reciprocal(image->indexed_cell, &asx, &asy, &asz, + &bsx, &bsy, &bsz, + &csx, &csy, &csz); + fprintf(ofh, "astar = %+9.7f %+9.7f %+9.7f nm^-1\n", + asx/1e9, asy/1e9, asz/1e9); + fprintf(ofh, "bstar = %+9.7f %+9.7f %+9.7f nm^-1\n", + bsx/1e9, bsy/1e9, bsz/1e9); + fprintf(ofh, "cstar = %+9.7f %+9.7f %+9.7f nm^-1\n", + csx/1e9, csy/1e9, csz/1e9); + + if ( image->f0_available ) { + fprintf(ofh, "I0 = %7.5f (arbitrary gas detector units)\n", + image->f0); + } else { + fprintf(ofh, "I0 = invalid\n"); + } + + fprintf(ofh, "photon_energy_eV = %f\n", + J_to_eV(ph_lambda_to_en(image->lambda))); + + write_peaks(image, ofh); + write_reflections(image, ofh); + + fprintf(ofh, "----- End chunk -----\n\n"); +} + + int find_chunk(FILE *fh, UnitCell **cell, char **filename, double *ev) { char line[1024]; diff --git a/src/stream.h b/src/stream.h index 5865aeca..b5f05b19 100644 --- a/src/stream.h +++ b/src/stream.h @@ -17,8 +17,11 @@ #endif +struct image; + extern int count_patterns(FILE *fh); extern int find_chunk(FILE *fh, UnitCell **cell, char **filename, double *ev); +extern void write_chunk(FILE *ofh, struct image *image, int flags); #endif /* STREAM_H */ -- cgit v1.2.3