diff options
| -rw-r--r-- | Makefile.am | 5 | ||||
| -rw-r--r-- | Makefile.in | 7 | ||||
| -rw-r--r-- | src/detector.c | 23 | ||||
| -rw-r--r-- | src/detector.h | 4 | ||||
| -rw-r--r-- | src/indexamajig.c | 1 | ||||
| -rw-r--r-- | src/powder_plot.c | 913 | ||||
| -rw-r--r-- | src/stream.c | 21 | ||||
| -rw-r--r-- | src/stream.h | 2 |
8 files changed, 924 insertions, 52 deletions
diff --git a/Makefile.am b/Makefile.am index e1ff7441..00c9fd8b 100644 --- a/Makefile.am +++ b/Makefile.am @@ -94,7 +94,8 @@ src_check_hkl_SOURCES = src/check_hkl.c src/cell.c src/utils.c \ src_powder_plot_SOURCES = src/powder_plot.c src/cell.c src/utils.c src/image.c \ src/hdf5-file.c src/detector.c src/thread-pool.c \ - src/beam-parameters.c + src/reflist-utils.c src/beam-parameters.c src/stream.c \ + src/reflist.c src/symmetry.c src_render_hkl_SOURCES = src/render_hkl.c src/cell.c src/reflist-utils.c \ src/utils.c src/povray.c src/symmetry.c src/render.c \ @@ -118,6 +119,8 @@ src_partialator_SOURCES = src/partialator.c src/cell.c src/hdf5-file.c \ src/symmetry.c src/post-refinement.c \ src/hrs-scaling.c src/reflist.c +src_stream_SOURCES = src/reflist-utils.c + if BUILD_CUBEIT src_cubeit_SOURCES = src/cubeit.c src/cell.c src/hdf5-file.c src/utils.c \ src/detector.c src/render.c src/filters.c src/image.c \ diff --git a/Makefile.in b/Makefile.in index 1efd5f8c..584d4ace 100644 --- a/Makefile.in +++ b/Makefile.in @@ -207,7 +207,9 @@ src_pattern_sim_DEPENDENCIES = $(top_builddir)/lib/libgnu.a am_src_powder_plot_OBJECTS = src/powder_plot.$(OBJEXT) \ src/cell.$(OBJEXT) src/utils.$(OBJEXT) src/image.$(OBJEXT) \ src/hdf5-file.$(OBJEXT) src/detector.$(OBJEXT) \ - src/thread-pool.$(OBJEXT) src/beam-parameters.$(OBJEXT) + src/thread-pool.$(OBJEXT) src/reflist-utils.$(OBJEXT) \ + src/beam-parameters.$(OBJEXT) src/stream.$(OBJEXT) \ + src/reflist.$(OBJEXT) src/symmetry.$(OBJEXT) src_powder_plot_OBJECTS = $(am_src_powder_plot_OBJECTS) src_powder_plot_LDADD = $(LDADD) src_powder_plot_DEPENDENCIES = $(top_builddir)/lib/libgnu.a @@ -669,7 +671,8 @@ src_check_hkl_SOURCES = src/check_hkl.c src/cell.c src/utils.c \ src_powder_plot_SOURCES = src/powder_plot.c src/cell.c src/utils.c src/image.c \ src/hdf5-file.c src/detector.c src/thread-pool.c \ - src/beam-parameters.c + src/reflist-utils.c src/beam-parameters.c src/stream.c \ + src/reflist.c src/symmetry.c src_render_hkl_SOURCES = src/render_hkl.c src/cell.c src/reflist-utils.c \ src/utils.c src/povray.c src/symmetry.c src/render.c \ diff --git a/src/detector.c b/src/detector.c index 5c1db570..e201c454 100644 --- a/src/detector.c +++ b/src/detector.c @@ -919,6 +919,29 @@ double largest_q(struct image *image) return qmax; } +double smallest_q(struct image *image) +{ + int fs, ss; + double ttm = +INFINITY; + double qmin = +INFINITY; + for ( fs=0; fs<image->width; fs++ ) { + for ( ss=0; ss<image->height; ss++ ) { + + struct rvec q; + double tt; + + q = get_q(image, fs, ss, &tt, 1.0/image->lambda); + + if ( tt < ttm ) { + qmin = modulus(q.u, q.v, q.w); + ttm = tt; + } + + } + } + + return qmin; +} void get_pixel_extents(struct detector *det, double *min_x, double *min_y, diff --git a/src/detector.h b/src/detector.h index fce60ff6..7073a47a 100644 --- a/src/detector.h +++ b/src/detector.h @@ -118,9 +118,9 @@ extern int reverse_2d_mapping(double x, double y, double *pfs, double *pss, extern double largest_q(struct image *image); -extern int write_detector_geometry(const char *filename, struct detector *det); - +extern double smallest_q(struct image *image); +extern int write_detector_geometry(const char *filename, struct detector *det); #endif /* DETECTOR_H */ diff --git a/src/indexamajig.c b/src/indexamajig.c index 569a41bd..4985ccbe 100644 --- a/src/indexamajig.c +++ b/src/indexamajig.c @@ -43,6 +43,7 @@ #include "reflist-utils.h" + /* Write statistics at APPROXIMATELY this interval */ #define STATS_EVERY_N_SECONDS (5) diff --git a/src/powder_plot.c b/src/powder_plot.c index c6d2f045..515c8c4e 100644 --- a/src/powder_plot.c +++ b/src/powder_plot.c @@ -3,7 +3,8 @@ * * Plot powder patterns * - * (c) 2006-2010 Thomas White <taw@physics.org> + * (c) 2011 Andrew Aquila <andrew.aquila@cfel.de> + * (c) 2006-2010 Thomas White <taw@physics.org> * * Part of CrystFEL - crystallography with a FEL * @@ -21,64 +22,611 @@ #include <unistd.h> #include <getopt.h> +#include "stream.h" +#include "reflist.h" #include "utils.h" #include "image.h" #include "detector.h" #include "index.h" #include "hdf5-file.h" #include "beam-parameters.h" +#include "reflist-utils.h" +#include "symmetry.h" +struct bin_stats { + unsigned int N; + double total; + double mean; + double std_dev; + double q_min; + double q_max; + double q_value; + double fract; +}; + +struct histogram_info { + double q_max; + double q_min; + double q_delta; + unsigned int histsize; + int spacing; //linear, q^2, & equal volume +}; + +enum { + PLOT_PEAKS, + PLOT_HKL, + PLOT_REFL, + PLOT_H5, + PLOT_D +}; + +enum { + FILE_STREAM, + FILE_HKL, + FILE_H5 +}; + +enum { + LINEAR, + q2, + VOLUME +}; + +static int find_q_bin_index(double q, struct histogram_info *info, + struct bin_stats *hist) +{ + /* bisection search alg. find q_bin index of order Log(n) time */ + int mid; + int min = 0; + int max = (*info).histsize-1; + if (q < hist[min].q_max) {return min;} + if (q > hist[max].q_min) {return max;} + do { + mid = (min + max) / 2; + if (q < hist[mid].q_min) { + max = mid; + } else if (q > hist[mid].q_max){ + min = mid ; + } else { + return mid; + } + } while(max - min > 1); + return mid; +} + +static int add_peak_to_histogram(double fs, double ss, double intensity, + struct image *image, struct histogram_info *info, + struct bin_stats *hist) +{ + /* used for h5 files, peak/feature list and stream positions*/ + /* See Knuth TAOCP vol 2, 3rd ed, pg 232 for running variance*/ + struct rvec r; + double q, delta; + int i; + + r = get_q(image, fs, ss, NULL, 1.0/ (*image).lambda); + q = modulus(r.u, r.v, r.w); + + /* ignore q value if outside of range */ + if ((q<(*info).q_min) || (q>(*info).q_max)) { + return 1; + } + i = find_q_bin_index(q, info, hist); + + delta = intensity - hist[i].mean; + hist[i].N++; + hist[i].total += intensity; + hist[i].mean = hist[i].mean + delta /hist[i].N; + hist[i].std_dev = hist[i].std_dev + + (delta *(intensity - hist[i].mean)); + return 0; +} + +static int add_d_to_histogram(double q, double intensity, + struct histogram_info *info, struct bin_stats *hist) +{ + /* used for d and hkl of stream files where redundancy = 1 */ + double delta; + int i; + + /* ignore q value if outside of range */ + if ((q<(*info).q_min) || (q>(*info).q_max)) { + return 1; + } + i = find_q_bin_index(q, info, hist); + + delta = intensity - hist[i].mean; + hist[i].N++; + hist[i].total += intensity; + hist[i].mean = hist[i].mean + delta /hist[i].N; + hist[i].std_dev = hist[i].std_dev + + (delta *(intensity - hist[i].mean)); + return 0; +} + +static int add_hkl_to_histogram(double q, double intensity, + int redundancy, int q_scaling, struct histogram_info *info, + struct bin_stats *hist) +{ + int i = 0; + + /* ignore q value if outside of range */ + if ((q<(*info).q_min) || (q>(*info).q_max)) { + return 1; + } + + /* The accounting is the intensity of the reflection times the + The number of occurance of that reflection smeered out over the + surface area which is 4 Pi q^2 the 4 pi is left out as it is a + common consant and the total is in arbiturary units. + */ + for (i=0; i<redundancy; i++) { + if (q_scaling) { + add_d_to_histogram(q, intensity/(q*q), info, hist); + } else { + add_d_to_histogram(q, intensity, info, hist); + } + } + return 0; +} + + +static int histogram_setup(struct histogram_info *info, + struct bin_stats *histdata) +{ + int i; + double x; + /* setup histogram*/ + if ((*info).spacing == LINEAR) {x = 1.0;} + else if ((*info).spacing == q2) {x = 2.0;} + else {x = 3.0;} + for ( i=0; i<(*info).histsize; i++) { + histdata[i].N = 0; + histdata[i].total = 0.0; + histdata[i].mean = 0.0; + histdata[i].std_dev = 0.0; + histdata[i].fract = 0.0; + histdata[i].q_min = pow(( i *(*info).q_delta) + + pow((*info).q_min, x), 1.0/x); + histdata[i].q_max = pow(((i+1.0)*(*info).q_delta) + + pow((*info).q_min, x), 1.0/x); + histdata[i].q_value = pow(((i+0.5)*(*info).q_delta) + + pow((*info).q_min, x), 1.0/x); + } + return 0; +} + +static int ring_fraction_calc(struct histogram_info *info, + struct bin_stats *hist, struct image *image) +{ + struct panel *p; + struct rvec r; + int i,fs,ss; + double q, q_fs, q_ss; + double pi = 3.14159265; + + /* check if detector geometry is loaded */ + if ( ((*image).det == NULL) || ((*image).lambda < 0.0)) { + return 1; + } + + /* loop over all pixels */ + for ( ss=0; ss<(*image).height; ss++ ) { + for ( fs=0; fs<(*image).width; fs++ ) { + r = get_q(image, fs, ss, NULL, 1.0/ (*image).lambda); + q = modulus(r.u, r.v, r.w); + + /* if pixel is valid (not a bad pixel and not out of range) */ + if ( (q>(*info).q_min) && (q<(*info).q_max) && + (in_bad_region((*image).det,fs,ss) == 0) ) { + + /* using the panel to determine the pixel size + takes care of edge pixels correctly */ + p = find_panel((*image).det, fs, ss); + r = get_q_for_panel(p, (fs+1)-(double)p->min_fs, + ss-(double)p->min_ss, + NULL, 1.0/(*image).lambda); + q_fs = modulus(r.u, r.v, r.w); + r = get_q_for_panel(p, fs-(double)p->min_fs, + (ss+1) -(double)p->min_ss, + NULL, 1.0/(*image).lambda); + q_ss = modulus(r.u, r.v, r.w); + i = find_q_bin_index(q, info, hist); + hist[i].fract = hist[i].fract + + fabs((q_fs - q) * (q_ss - q)); //approx pixel size + } + } + } + /* divide measured area by ring area */ + for ( i=0; i<(*info).histsize; i++) { + hist[i].fract = hist[i].fract / + (pi * (pow(hist[i].q_max,2.0) - pow(hist[i].q_min,2.0))); + } + return 0; +} + +static unsigned int process_h5(struct image *image, struct histogram_info *info, + struct bin_stats *histdata, unsigned int processing_total) +{ + int fs, ss; + double intensity; + for ( ss=0; ss<(*image).height; ss++ ) { + for ( fs=0; fs<(*image).width; fs++ ) { + intensity = (*image).data[fs + (*image).width*ss]; + progress_bar(fs + (*image).width*ss, processing_total, + "Processing"); + if (in_bad_region((*image).det,fs,ss) == 0) { + add_peak_to_histogram(fs, ss, intensity, + image, info, histdata); + } + } + } + return 0; +} + +static unsigned int process_hkl(struct image *image, char *sym, UnitCell *cell, + struct histogram_info *info, struct bin_stats *histdata, + int q_scaling,unsigned int processing_total, int use_redundency) +{ + Reflection *refl; + RefListIterator *iter; + unsigned int i = 0; + unsigned int n_peaks = 0; + int h,k,l,redundancy; + double q, intensity; + for ( refl = first_refl((*image).reflections, &iter); + refl != NULL; + refl = next_refl(refl, iter) ) { + i++; + progress_bar(i, processing_total, "Processing"); + get_indices(refl, &h, &k, &l); + intensity = get_intensity(refl); + if (use_redundency == 1) { + redundancy = get_redundancy(refl); + } else { + redundancy = num_equivs(h, k, l, sym); + } + /*note: resolution returns 1/d/2 in m */ + q = 2.0 * resolution(cell, h, k, l); + //sigma = get_esd_intensity(refl); + add_hkl_to_histogram(q, intensity, redundancy, + q_scaling, info, histdata); + n_peaks+=redundancy; + } + return n_peaks; +} + +static unsigned int process_stream_reflection(FILE *fh, struct image *image, + struct histogram_info *info, struct bin_stats *histdata, + unsigned int processing_total, unsigned int *n_patterns) +{ + int rval; + unsigned int i = 0; + unsigned int n_peaks = 0; + Reflection *refl; + RefListIterator *iter; + double intensity, fs_double, ss_double; + do { + /* Get data from next chunk */ + rval = read_chunk(fh, image); + if ( rval ) continue; + i++; + progress_bar(i, processing_total, "Processing"); + // check if the pattern indexed, if so use those peaks + if ((*image).reflections != NULL) { + (*n_patterns)++; //inc number of patterns used + for ( refl = first_refl((*image).reflections, &iter); + refl != NULL; + refl = next_refl(refl, iter) ) { + /* note added fs_double as fs is an int */ + intensity = get_intensity(refl); + get_detector_pos(refl, &fs_double, &ss_double); + if(add_peak_to_histogram(fs_double, ss_double, + intensity, image, info, histdata) == 0) { + n_peaks++; + } + } + } + free((*image).filename); + reflist_free((*image).reflections); + image_feature_list_free((*image).features); + cell_free((*image).indexed_cell); + } while (rval == 0); + return n_peaks; +} + +static unsigned int process_stream_d(FILE *fh, struct image *image, + struct histogram_info *info, struct bin_stats *histdata, + unsigned int processing_total, unsigned int *n_patterns) +{ + int h, k, l, rval; + unsigned int i = 0; + unsigned int n_peaks = 0; + Reflection *refl; + RefListIterator *iter; + double intensity, q; + do { + /* Get data from next chunk */ + rval = read_chunk(fh, image); + if ( rval ) continue; + i++; + progress_bar(i, processing_total, "Processing"); + // check if the pattern indexed, if so use those peaks + if ((*image).reflections != NULL) { + (*n_patterns)++; //inc number of patterns used + for ( refl = first_refl((*image).reflections, &iter); + refl != NULL; + refl = next_refl(refl, iter) ) { + get_indices(refl, &h, &k, &l); + intensity = get_intensity(refl); + /*resolution returns 1/d/2 in nm */ + q = 2.0 * resolution((*image).indexed_cell, h, k, l); + if (add_d_to_histogram(q, intensity, info, + histdata) == 0) { + n_peaks++; + } + } + } + free((*image).filename); + reflist_free((*image).reflections); + image_feature_list_free((*image).features); + cell_free((*image).indexed_cell); + } while (rval == 0); + return n_peaks; +} + + +static unsigned int process_stream_hkl(FILE *fh, struct image *image, + struct histogram_info *info, struct bin_stats *histdata, UnitCell *cell, + unsigned int processing_total, unsigned int *n_patterns) +{ + int h, k, l, rval; + unsigned int i = 0; + unsigned int n_peaks = 0; + Reflection *refl; + RefListIterator *iter; + double intensity, q; + do { + /* Get data from next chunk */ + rval = read_chunk(fh, image); + if ( rval ) continue; + i++; + progress_bar(i, processing_total, "Processing"); + if ((*image).reflections != NULL) { + n_patterns++; //inc number of patterns used + for ( refl = first_refl((*image).reflections, &iter); + refl != NULL; + refl = next_refl(refl, iter) ) { + get_indices(refl, &h, &k, &l); + intensity = get_intensity(refl); + q = 2.0 * resolution(cell, h, k, l); + if(add_d_to_histogram(q, intensity, info, + histdata) == 0) { + n_peaks++; + } + } + } + free((*image).filename); + reflist_free((*image).reflections); + image_feature_list_free((*image).features); + cell_free((*image).indexed_cell); + } while (rval == 0); + return n_peaks; +} + +static unsigned int process_stream_peaks(FILE *fh, struct image *image, + struct histogram_info *info, struct bin_stats *histdata, + unsigned int processing_total, unsigned int *n_patterns, int only_indexed) +{ + struct imagefeature *f; + int rval; + unsigned int i = 0; + unsigned int n_peaks = 0; + unsigned int j; + do { + /* Get data from next chunk */ + rval = read_chunk(fh, image); + if ( rval ) continue; + i++; + progress_bar(i, processing_total, "Processing"); + if ((*image).features != NULL) { + if ((only_indexed == 0) || + ((only_indexed == 1) && + ((*image).reflections != NULL))) { + (*n_patterns)++; + for (j = 0; + j<image_feature_count((*image).features); + j++) { + f = image_get_feature((*image).features, j); + if (in_bad_region((*image).det, + (*f).fs,(*f).ss) == 0) { + if (add_peak_to_histogram((*f).fs, + (*f).ss, (*f).intensity, image, + info, histdata) == 0) { + n_peaks++; + } + } + } + } + } + free((*image).filename); + reflist_free((*image).reflections); + image_feature_list_free((*image).features); + cell_free((*image).indexed_cell); + } while (rval == 0); + return n_peaks; +} + +static unsigned int process_stream_h5(FILE *fh, struct image *image, + struct histogram_info *info, struct bin_stats *histdata, int config_satcorr, + int only_indexed, unsigned int processing_total, unsigned int *n_patterns) +{ + int fs, ss, rval; + double intensity; + unsigned int i = 0; + unsigned int n_peaks = 0; + struct hdfile *hdfile = NULL; + do { + /* Get data from next chunk */ + rval = read_chunk(fh, image); + if ( rval ) continue; + i++; + progress_bar(i, processing_total, "Processing"); + if ((only_indexed == 0) || + ((only_indexed == 1) && + ((*image).reflections != NULL))) { + hdfile = hdfile_open((*image).filename); + hdfile_set_image(hdfile, "/data/data"); + hdf5_read(hdfile, image, config_satcorr); + hdfile_close(hdfile); + n_patterns++; + for ( ss=0; ss<(*image).height; ss++ ) { + for ( fs=0; fs<(*image).width; fs++ ) { + intensity = (*image).data[fs + (*image).width*ss]; + if (in_bad_region((*image).det,fs,ss) == 0) { + add_peak_to_histogram(fs, ss, intensity, + image, info, histdata); + } + } + } + } + free((*image).data); + free((*image).filename); + reflist_free((*image).reflections); + image_feature_list_free((*image).features); + cell_free((*image).indexed_cell); + } while (rval == 0); + return n_peaks; +} static void show_help(const char *s) { - printf("Syntax: %s [options] <file.h5>\n\n", s); + printf("Syntax: %s [options]\n\n", s); printf( "Plot a powder pattern as a 1D graph using the detector geometry.\n" "\n" -" -h, --help Display this help message.\n" -" -g. --geometry=<file> Get detector geometry from file.\n" -" -i, --input=<file> Input filename.\n" -" -b, --beam=<file> Get beam parameters (specifically, wavelength)" -" from file if non present in the HDF5 file.\n" +" -h, --help Display this help message.\n" +" -i, --input=<file> Input filename. (*.stream, *.hkl, or *.h5)\n" +" -o, --output=<file> Output filename. (default stdout)\n" +" -g. --geometry=<file> Get detector geometry from file.\n" +" -b, --beam=<file> Get beam parameters (wavelength) from file.\n" +" -p, --pdb=<file> Get unit cell from PDB file. (.hkl files only)\n" +" -y, --symmetry=<sym> The symmetry of crystal (.hkl files only)\n" +" -s, --bins=n Makes histogram wiht n bins (default is 100).\n" +" --spacing=<type> Use 'type' to select the q spacing.\n" +" Choose from:\n" +" linear : linear (default)\n" +" q2 : even spacing in Wilson plots\n" +" volume : constant volume spacing\n" +" --q-max=n The max q to be considered in plot.\n" +" --q-min=n The min q to be considered in plot.\n" +" -d, --data=<type> Use to select the kind of stream data in histogram.\n" +" Choose from:\n" +" reflection : uses peak positons from indexed\n" +" reflection \n" +" hkl : uses the hkl list from indexed\n" +" reflections (requires pdb file)\n" +" d : uses the 1/d list from indexed\n" +" reflections (default)\n" +" peaks : uses all peaks found from peak\n" +" search\n" +" h5 : all points in h5, excluding bad\n" +" regions\n" +" --no-sat-corr Don't correct values of saturated peaks using a\n" +" table included in the HDF5 file.\n" +" --only-indexed Use wiht -data=peaks or h5 if you want to use the\n" +" peak list of only indexed patterns\n" +" --no-q-scaling Use with .hkl files if you want to not scale the\n" +" powder by 1/q^2\n" +" --ring-corr Use if you want to scale the powder plot to\n" +" correct for the fractional area sampled of the\n" +" powder ring\n" +" --use-redundency Use with .hkl files if you want to use the number\n" +" of measurements and not the number of symetrical\n" +" equivelent reflections as the number of time a\n" +" reflection occurs in the powder\n" "\n"); } - int main(int argc, char *argv[]) { - int c; - struct image image; - int x, y; - struct hdfile *hdfile; + FILE *fh = NULL; + UnitCell *cell = NULL; + struct image image; + struct hdfile *hdfile = NULL; + struct bin_stats *histdata = NULL; + struct histogram_info hist_info; + + hist_info.histsize = 100; //default settings + hist_info.q_min = -1.0; + hist_info.q_max = -1.0; + hist_info.spacing = LINEAR; + image.lambda = -1.0; + image.beam = NULL; + image.det = NULL; + + unsigned int n_patterns = 0; + unsigned int n_peaks = 0; + + int c, rval, file_type, data_type; + int config_satcorr = 1; //true by default + int need_geometry = 0; //false + int need_beam = 0; //false + int need_pdb = 0; //false + int only_indexed = 0; //false + int q_scaling = 1; //true by default + int ring_corr = 0; //false + int use_redundency = 0; //false + unsigned int i, processing_total; + char *filename = NULL; char *geometry = NULL; - char *beamf = NULL; - struct beam_params *beam = NULL; + char *beamf = NULL; + char *pdb = NULL; + char *output = NULL; + char *datatype = NULL; + char *sym = NULL; /* Long options */ const struct option longopts[] = { {"help", 0, NULL, 'h'}, {"input", 1, NULL, 'i'}, + {"output", 1, NULL, 'o'}, {"geometry", 1, NULL, 'g'}, {"beam", 1, NULL, 'b'}, + {"pdb", 1, NULL, 'p'}, + {"symmetry", 1, NULL, 'y'}, + {"bins", 1, NULL, 's'}, + {"q-max", 1, NULL, 1 }, + {"q-min", 1, NULL, 2 }, + {"spacing", 1, NULL, 3 }, + {"no-sat-corr", 0, &config_satcorr, 0 }, + {"sat-corr", 0, &config_satcorr, 1 }, + {"only-indexed", 0, &only_indexed, 1 }, + {"no-q-scaling", 0, &q_scaling, 0 }, + {"ring-corr", 0, &ring_corr, 1 }, + {"use-redundency", 0, &use_redundency, 1 }, + {"data", 1, NULL, 'd'}, {0, 0, NULL, 0} }; /* Short options */ - while ((c = getopt_long(argc, argv, "hi:g:", longopts, NULL)) != -1) { + while ((c = getopt_long(argc, argv, "hi:o:g:b:p:s:d:y:", + longopts, NULL)) != -1) { switch (c) { case 'h' : show_help(argv[0]); return 0; - case 0 : - break; - case 'i' : filename = strdup(optarg); break; + case 'o' : + output = strdup(optarg); + break; + case 'g' : geometry = strdup(optarg); break; @@ -87,65 +635,338 @@ int main(int argc, char *argv[]) beamf = strdup(optarg); break; + case 'p' : + pdb = strdup(optarg); + break; + + case 'y' : + sym = strdup(optarg); + break; + + case 's' : + hist_info.histsize = atoi(optarg); + break; + + case 1 : + hist_info.q_max = atof(optarg); + break; + + case 2 : + hist_info.q_min = atof(optarg); + break; + + case 3 : + if (strcmp(optarg, "linear") == 0 ) { + hist_info.spacing = LINEAR; + } else if (strcmp(optarg, "q2") == 0 ) { + hist_info.spacing = q2; + } else if (strcmp(optarg, "volume") == 0) { + hist_info.spacing = VOLUME; + } else { + ERROR("Failed to read spacing plot type: '%s'\n", + optarg); + return 1; + } + break; + + case 'd' : + datatype = strdup(optarg); + break; + + case 0 : + break; + default : return 1; } - } + /* Process input file type */ if ( filename == NULL ) { ERROR("You must specify the input filename with -i\n"); return 1; + } else if ( is_stream(filename) == 1 ) { //open .stream file + file_type = FILE_STREAM; + fh = fopen(filename, "r"); + if ( fh == NULL ) { + ERROR("Failed to open input file\n"); + return 1; + } + rval = read_chunk(fh, &image); //read chunk to get wavelength + rewind(fh); + processing_total = count_patterns(fh); + rewind(fh); + + } else if ( H5Fis_hdf5(filename) > 0) { //open .h5 file + file_type = FILE_H5; + need_geometry = 1; //true + hdfile = hdfile_open(filename); + hdfile_set_image(hdfile, "/data/data"); + hdf5_read(hdfile, &image, config_satcorr); + hdfile_close(hdfile); + processing_total = image.width * image.height; + } else { + image.reflections = read_reflections(filename); + if (image.reflections != NULL) { //open .hkl file + file_type = FILE_HKL; + need_pdb = 1; + processing_total = num_reflections(image.reflections); + } else { + ERROR("Input file: %s is an invalid type.\n",filename); + return 1; + } } + free(filename); - if ( geometry == NULL ) { - ERROR("You need to specify a geometry file with --geometry\n"); + /*Process data read type */ + if ( datatype == NULL ) { + data_type = PLOT_D; + if ((hist_info.q_min < 0.0) || (hist_info.q_max < 0.0)) { + need_geometry = 1; + } + } else if ( strcmp(datatype, "reflection") == 0 ) { + data_type = PLOT_REFL; + need_geometry = 1; + } else if ( strcmp(datatype, "hkl") == 0 ) { + data_type = PLOT_HKL; + need_pdb = 1; + if ((hist_info.q_min <= 0.0) || (hist_info.q_max <= 0.0)) { + need_geometry = 1; + } + } else if ( strcmp(datatype, "d") == 0 ) { + data_type = PLOT_D; + if ((hist_info.q_min <= 0.0) || (hist_info.q_max <= 0.0)) { + need_geometry = 1; + } + } else if ( strcmp(datatype, "peaks") == 0 ) { + data_type = PLOT_PEAKS; + need_geometry = 1; + } + else if ( strcmp(datatype, "h5") == 0 ) { + data_type = PLOT_H5; + need_geometry = 1; + } else { + ERROR("Failed to read data plot type: '%s'\n", datatype); return 1; } - if ( beamf != NULL ) { - beam = get_beam_parameters(beamf); - free(beamf); + /* logic checks */ + if ((need_geometry) && (image.lambda < 0.0)) { + need_beam = 1; } - - image.det = get_detector_geometry(geometry); - if ( image.det == NULL ) { - ERROR("Failed to read detector geometry from '%s'\n", geometry); + if ( hist_info.histsize <= 0 ) { + ERROR("You need to specify a histogram with more then 0 " + "bins\n"); + return 1; + } + if (hist_info.q_min > hist_info.q_max) { + ERROR("the minimum q value of: %e " + "is greator then your max q value of: %e\n", + hist_info.q_min, hist_info.q_max); return 1; } + + /*get geometry, beam and pdb files and parameters as needed */ + if (need_geometry) { + if (geometry == NULL) { + ERROR("You need to specify a geometry file with " + "--geometry\n"); + return 1; + } else { + image.det = get_detector_geometry(geometry); + if ( image.det == NULL ) { + ERROR("Failed to read detector geometry " + "from '%s'\n", geometry); + return 1; + } + image.width = image.det->max_fs; + image.height = image.det->max_ss; + } + } free(geometry); + + if (need_beam) { + if ( beamf == NULL ) { + ERROR("No wavelength in file, so you need to specify " + "a beam parameters file with --beam\n"); + return 1; + } else { + image.beam = get_beam_parameters(beamf); + if ( image.beam == NULL ) { + ERROR("Failed to read beam from '%s'\n", + beamf); + return 1; + } + image.lambda = ph_en_to_lambda(eV_to_J( + image.beam->photon_energy)); + } + } + free(beamf); - hdfile = hdfile_open(filename); - hdfile_set_image(hdfile, "/data/data"); - hdf5_read(hdfile, &image, 1); - if ( image.lambda < 0.0 ) { - if ( beam != NULL ) { - image.lambda = beam->photon_energy; + if (need_pdb) { + if (pdb == NULL) { + ERROR("You need to specify a pdb file with --pdb.\n"); + return 1; } else { - ERROR("No wavelength in file, so you need to give " - "a beam parameters file with -b.\n"); + cell = load_cell_from_pdb(pdb); + if ( cell == NULL ) { + ERROR("Couldn't read unit cell (from %s)\n", + pdb); + return 1; + } + } + } + free(pdb); + + if ( sym == NULL ) { + sym = strdup("1"); + } + + /* setup histogram info*/ + if (hist_info.q_min <= 0 ) { + hist_info.q_min = smallest_q(&image); + } + if (hist_info.q_max <= 0 ) { + hist_info.q_max = largest_q(&image); + } + if (hist_info.spacing == LINEAR) { + hist_info.q_delta = (hist_info.q_max - hist_info.q_min)/ + hist_info.histsize; + } else if (hist_info.spacing == q2) { + hist_info.q_delta = (pow(hist_info.q_max, 2.0) - + pow(hist_info.q_min, 2.0)) / + hist_info.histsize; + } else { //by default must be in VOLUME + hist_info.q_delta = (pow(hist_info.q_max, 3.0) - + pow(hist_info.q_min, 3.0)) / + hist_info.histsize; + } + /* setup histogram data */ + histdata = malloc((hist_info.histsize) * sizeof(struct bin_stats)); + histogram_setup(&hist_info, histdata); + + /* setup ring scaling */ + if (ring_corr == 1) { + if (ring_fraction_calc(&hist_info, histdata, &image) == 1) { + ERROR("Detector is not loaded could not correct for" + " finiate ring measurement. Do not use --ring-corr\n"); return 1; } } - for ( x=0; x<image.width; x++ ) { - for ( y=0; y<image.height; y++ ) { + /* process reflections based on file type and data type*/ + switch (file_type) { + case FILE_H5 : + n_patterns++; + n_peaks = process_h5(&image, &hist_info, histdata, processing_total); + free(image.data); + break; - double q; - int intensity; - struct rvec r; + case FILE_HKL : + n_patterns++; //inc number of patterns used + n_peaks = process_hkl(&image, sym, cell, &hist_info, histdata, + q_scaling, processing_total, use_redundency); + break; - r = get_q(&image, x, y, NULL, 1.0/image.lambda); - q = modulus(r.u, r.v, r.w); + case FILE_STREAM : + switch (data_type) { + case PLOT_REFL : + n_peaks = process_stream_reflection(fh, &image, + &hist_info, histdata, processing_total, + &n_patterns); + break; + case PLOT_D : + n_peaks = process_stream_d(fh, &image, &hist_info, + histdata, processing_total, + &n_patterns); + break; + case PLOT_HKL : + n_peaks = process_stream_hkl(fh, &image, &hist_info, + histdata, cell, processing_total, + &n_patterns); + break; + case PLOT_PEAKS : + n_peaks = process_stream_peaks(fh, &image, &hist_info, + histdata, processing_total, + &n_patterns, only_indexed); + break; + case PLOT_H5 : + n_peaks = process_stream_h5(fh, &image, &hist_info, + histdata, config_satcorr, only_indexed, + processing_total, &n_patterns); + break; + default : + break; + } + fclose(fh); + break; - intensity = image.data[x + image.width*y]; + default : + break; - printf("%5i\t%5i\t%7.3f\t%7i\n", x, y, q/1.0e9, intensity); + } + + /* sqrt the variance to get the std_dev */ + + for(i=0; i<hist_info.histsize; i++) { + if (histdata[i].N > 1) { + histdata[i].std_dev = sqrt(histdata[i].std_dev/ + (histdata[i].N-1)); + } + } + if (ring_corr == 1) { + for(i=0; i<hist_info.histsize; i++) { + histdata[i].N = histdata[i].N / + histdata[i].fract; + histdata[i].total = histdata[i].total / + histdata[i].fract; + histdata[i].mean = histdata[i].mean / + histdata[i].fract; + histdata[i].std_dev = histdata[i].total / + histdata[i].fract; + } + } + /* print out the results */ + if (output != NULL) { + fh = fopen(output, "w"); + if ( fh == NULL ) { + ERROR("Failed to open output file\n"); + return 1; + } + fprintf(fh,"I read %i patterns with %i peaks\n", + n_patterns, n_peaks); + fprintf(fh,"q\tN\ttotal\tmean\tstd dev\t std dev of mean\n"); } + else { + printf("I read %i patterns with %i peaks\n", + n_patterns, n_peaks); + printf("q\tN\ttotal\tmean\tstd dev\t std dev of mean\n"); + } + for(i=0; i<hist_info.histsize; i++) { + if (output != NULL) { + fprintf(fh,"%5e\t%i\t%5e\t%5e\t%5e\t%5e\n", + histdata[i].q_min, histdata[i].N, histdata[i].total, + histdata[i].mean, histdata[i].std_dev, + histdata[i].std_dev/sqrt(histdata[i].N)); + } + else { + printf("%5e\t%i\t%5e\t%5e\t%5e\t%5e\n", + histdata[i].q_min, histdata[i].N, histdata[i].total, + histdata[i].mean, histdata[i].std_dev, + histdata[i].std_dev/sqrt(histdata[i].N)); + } } - free(beamf); + if ( cell != NULL ) cell_free(cell); + if ( image.det != NULL ) free(image.det); + if ( image.beam != NULL ) free(image.beam); + if (output != NULL) { + fclose(fh); + free(output); + } + free(histdata); + free(sym); return 0; } diff --git a/src/stream.c b/src/stream.c index 99344404..e3bc3c87 100644 --- a/src/stream.c +++ b/src/stream.c @@ -168,7 +168,7 @@ static int read_peaks(FILE *fh, struct image *image) first = 0; if ( r == 4 ) { image_add_feature(image->features, x, y, - image, 1.0, NULL); + image, intensity, NULL); } } while ( rval != NULL ); @@ -450,3 +450,22 @@ int skip_some_files(FILE *fh, int n) return 1; } + +int is_stream(const char *filename) { + FILE *fh; + char line[1024]; + char *rval = NULL; + fh = fopen(filename, "r"); + rval = fgets(line, 1023, fh); + fclose(fh); + if ( rval == NULL ) { + return -1; + } + if ( strncmp(line, "CrystFEL stream format 2.0", 26) == 0 ) { + return 1; + } + else { + return 0; + } + return -1; +} diff --git a/src/stream.h b/src/stream.h index 02a01c48..e4e71728 100644 --- a/src/stream.h +++ b/src/stream.h @@ -43,4 +43,6 @@ extern int read_chunk(FILE *fh, struct image *image); extern int skip_some_files(FILE *fh, int n); +extern int is_stream(const char *filename); + #endif /* STREAM_H */ |