Commit d004d970 authored by Thomas White's avatar Thomas White
Browse files

New stream writing

parent f27607b8
......@@ -45,7 +45,7 @@ src_process_hkl_SOURCES = src/process_hkl.c src/sfac.c src/statistics.c \
src/cell.c src/utils.c src/reflections.c \
src/symmetry.c src/stream.c src/beam-parameters.c \
src/thread-pool.c src/image.c src/detector.c \
src/hdf5-file.c
src/hdf5-file.c src/reflist.c
src_indexamajig_SOURCES = src/indexamajig.c src/hdf5-file.c src/utils.c \
src/cell.c src/image.c src/peaks.c src/index.c \
......@@ -104,11 +104,13 @@ src_partialator_SOURCES = src/partialator.c src/cell.c src/hdf5-file.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 \
src/symmetry.c src/stream.c src/thread-pool.c
src/symmetry.c src/stream.c src/thread-pool.c src/reflist.c
endif
src_estimate_background_SOURCES = src/estimate_background.c src/stream.c \
src/utils.c src/cell.c src/thread-pool.c
src/utils.c src/cell.c src/thread-pool.c \
src/image.c src/detector.c src/hdf5-file.c \
src/reflist.c
tests_list_check_SOURCES = tests/list_check.c src/reflist.c src/thread-pool.c \
src/utils.c
......
......@@ -106,21 +106,25 @@ src_compare_hkl_LDADD = $(LDADD)
src_compare_hkl_DEPENDENCIES = $(top_builddir)/lib/libgnu.a
am__src_cubeit_SOURCES_DIST = 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 src/symmetry.c src/stream.c src/thread-pool.c
src/image.c src/symmetry.c src/stream.c src/thread-pool.c \
src/reflist.c
@BUILD_CUBEIT_TRUE@am_src_cubeit_OBJECTS = src/cubeit.$(OBJEXT) \
@BUILD_CUBEIT_TRUE@ src/cell.$(OBJEXT) src/hdf5-file.$(OBJEXT) \
@BUILD_CUBEIT_TRUE@ src/utils.$(OBJEXT) src/detector.$(OBJEXT) \
@BUILD_CUBEIT_TRUE@ src/render.$(OBJEXT) src/filters.$(OBJEXT) \
@BUILD_CUBEIT_TRUE@ src/image.$(OBJEXT) src/symmetry.$(OBJEXT) \
@BUILD_CUBEIT_TRUE@ src/stream.$(OBJEXT) \
@BUILD_CUBEIT_TRUE@ src/thread-pool.$(OBJEXT)
@BUILD_CUBEIT_TRUE@ src/thread-pool.$(OBJEXT) \
@BUILD_CUBEIT_TRUE@ src/reflist.$(OBJEXT)
src_cubeit_OBJECTS = $(am_src_cubeit_OBJECTS)
src_cubeit_LDADD = $(LDADD)
src_cubeit_DEPENDENCIES = $(top_builddir)/lib/libgnu.a
am_src_estimate_background_OBJECTS = \
src/estimate_background.$(OBJEXT) src/stream.$(OBJEXT) \
src/utils.$(OBJEXT) src/cell.$(OBJEXT) \
src/thread-pool.$(OBJEXT)
src/thread-pool.$(OBJEXT) src/image.$(OBJEXT) \
src/detector.$(OBJEXT) src/hdf5-file.$(OBJEXT) \
src/reflist.$(OBJEXT)
src_estimate_background_OBJECTS = \
$(am_src_estimate_background_OBJECTS)
src_estimate_background_LDADD = $(LDADD)
......@@ -205,7 +209,7 @@ am_src_process_hkl_OBJECTS = src/process_hkl.$(OBJEXT) \
src/symmetry.$(OBJEXT) src/stream.$(OBJEXT) \
src/beam-parameters.$(OBJEXT) src/thread-pool.$(OBJEXT) \
src/image.$(OBJEXT) src/detector.$(OBJEXT) \
src/hdf5-file.$(OBJEXT)
src/hdf5-file.$(OBJEXT) src/reflist.$(OBJEXT)
src_process_hkl_OBJECTS = $(am_src_process_hkl_OBJECTS)
src_process_hkl_LDADD = $(LDADD)
src_process_hkl_DEPENDENCIES = $(top_builddir)/lib/libgnu.a
......@@ -611,7 +615,7 @@ src_process_hkl_SOURCES = src/process_hkl.c src/sfac.c src/statistics.c \
src/cell.c src/utils.c src/reflections.c \
src/symmetry.c src/stream.c src/beam-parameters.c \
src/thread-pool.c src/image.c src/detector.c \
src/hdf5-file.c
src/hdf5-file.c src/reflist.c
src_indexamajig_SOURCES = src/indexamajig.c src/hdf5-file.c \
src/utils.c src/cell.c src/image.c src/peaks.c src/index.c \
......@@ -661,10 +665,12 @@ src_partialator_SOURCES = src/partialator.c src/cell.c src/hdf5-file.c \
@BUILD_CUBEIT_TRUE@src_cubeit_SOURCES = src/cubeit.c src/cell.c src/hdf5-file.c src/utils.c \
@BUILD_CUBEIT_TRUE@ src/detector.c src/render.c src/filters.c src/image.c \
@BUILD_CUBEIT_TRUE@ src/symmetry.c src/stream.c src/thread-pool.c
@BUILD_CUBEIT_TRUE@ src/symmetry.c src/stream.c src/thread-pool.c src/reflist.c
src_estimate_background_SOURCES = src/estimate_background.c src/stream.c \
src/utils.c src/cell.c src/thread-pool.c
src/utils.c src/cell.c src/thread-pool.c \
src/image.c src/detector.c src/hdf5-file.c \
src/reflist.c
tests_list_check_SOURCES = tests/list_check.c src/reflist.c src/thread-pool.c \
src/utils.c
......@@ -843,6 +849,8 @@ src/filters.$(OBJEXT): src/$(am__dirstamp) \
src/$(DEPDIR)/$(am__dirstamp)
src/stream.$(OBJEXT): src/$(am__dirstamp) \
src/$(DEPDIR)/$(am__dirstamp)
src/reflist.$(OBJEXT): src/$(am__dirstamp) \
src/$(DEPDIR)/$(am__dirstamp)
src/cubeit$(EXEEXT): $(src_cubeit_OBJECTS) $(src_cubeit_DEPENDENCIES) src/$(am__dirstamp)
@rm -f src/cubeit$(EXEEXT)
$(AM_V_CCLD)$(LINK) $(src_cubeit_OBJECTS) $(src_cubeit_LDADD) $(LIBS)
......@@ -876,8 +884,6 @@ src/mosflm.$(OBJEXT): src/$(am__dirstamp) \
src/$(DEPDIR)/$(am__dirstamp)
src/geometry.$(OBJEXT): src/$(am__dirstamp) \
src/$(DEPDIR)/$(am__dirstamp)
src/reflist.$(OBJEXT): src/$(am__dirstamp) \
src/$(DEPDIR)/$(am__dirstamp)
src/diffraction-gpu.$(OBJEXT): src/$(am__dirstamp) \
src/$(DEPDIR)/$(am__dirstamp)
src/cl-utils.$(OBJEXT): src/$(am__dirstamp) \
......
......@@ -61,12 +61,7 @@ In addition, there is also:
- estimate_background, for calculating signal to noise ratios from
the indexed peaks.
- calibrate_detector, which splits
- reintegrate, which is like "indexamajig" but without the indexing
step, instead getting the orientation matrix from the
output of a previous run of either "indexamajig" or
"reintegrate".
- calibrate_detector, which does nothing useful at the moment.
Included at no extra cost are:
......
......@@ -292,21 +292,21 @@ static double get_wavelength(struct hdfile *f)
}
static double get_f0(struct hdfile *f)
static double get_i0(struct hdfile *f)
{
herr_t r;
hid_t dh;
double f0;
double i0;
dh = H5Dopen2(f->fh, "/LCLS/f_11_ENRC", H5P_DEFAULT);
if ( dh < 0 ) return -1.0;
r = H5Dread(dh, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL,
H5P_DEFAULT, &f0);
H5P_DEFAULT, &i0);
H5Dclose(dh);
if ( r < 0 ) return -1.0;
return f0;
return i0;
}
......@@ -434,13 +434,13 @@ int hdf5_read(struct hdfile *f, struct image *image, int satcorr)
/* Read wavelength from file */
image->lambda = get_wavelength(f);
image->f0 = get_f0(f);
if ( image->f0 < 0.0 ) {
image->i0 = get_i0(f);
if ( image->i0 < 0.0 ) {
ERROR("Couldn't read incident intensity - using 1.0.\n");
image->f0 = 1.0;
image->f0_available = 0;
image->i0 = 1.0;
image->i0_available = 0;
} else {
image->f0_available = 1;
image->i0_available = 1;
}
if ( satcorr ) debodge_saturation(f, image);
......
......@@ -76,8 +76,8 @@ struct image {
double lambda; /* Wavelength in m */
double div; /* Divergence in radians */
double bw; /* Bandwidth as a fraction */
double f0; /* Incident intensity */
int f0_available; /* 0 if f0 wasn't available
double i0; /* Incident intensity */
int i0_available; /* 0 if f0 wasn't available
* from the input. */
double osf; /* Overall scaling factor */
double profile_radius; /* Radius of reflection */
......
......@@ -52,7 +52,7 @@ struct static_index_args
int config_cmfilter;
int config_noisefilter;
int config_verbose;
int stream_flags;
int stream_flags; /* What goes into the output? */
int config_polar;
int config_satcorr;
int config_closer;
......@@ -62,7 +62,7 @@ struct static_index_args
struct detector *det;
IndexingMethod *indm;
IndexingPrivate **ipriv;
int peaks;
int peaks; /* Peak detection method */
int cellr;
struct beam_params *beam;
const char *element;
......@@ -132,15 +132,20 @@ static void show_help(const char *s)
" in the HDF5 file.\n"
"\n\n"
"You can control what information is included in the output stream using\n"
"' --record=<flags>'. 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"
"' --record=<flags>'. Possible flags are:\n\n"
" pixels Include a list of sums of pixel values within the\n"
" integration domain, correcting for individual pixel\n"
" solid angles.\n"
"\n"
" integrated Include a list of reflection intensities, produced by\n"
" integrating around predicted peak locations.\n"
"\n"
" peaks Include peak locations and intensities from the peak search.\n"
"\n"
" peaksifindexed As 'peaks', but only if the pattern could be indexed.\n\n"
"\n"
"The default is '--record=integrated'. The flags 'pixels' and 'integrated'\n"
"are mutually exclusive, as are the flags 'peaks' and 'peaksifindexed'.\n"
"\n\n"
"For more control over the process, you might need:\n\n"
" --cell-reduction=<m> Use <m> as the cell reduction method. Choose from:\n"
......@@ -162,7 +167,7 @@ static void show_help(const char *s)
" --min-gradient=<n> Minimum gradient for Zaefferer peak search.\n"
" Default: 100,000.\n"
" -e, --image=<element> Use this image from the HDF5 file.\n"
" Example: /data/data0."
" Example: /data/data0.\n"
" Default: The first one found.\n"
"\n"
"\nOptions for greater performance or verbosity:\n\n"
......@@ -431,6 +436,7 @@ int main(int argc, char *argv[])
struct beam_params *beam = NULL;
char *element = NULL;
double nominal_photon_energy;
int stream_flags = STREAM_INTEGRATED;
/* Long options */
const struct option longopts[] = {
......@@ -460,6 +466,7 @@ int main(int argc, char *argv[])
{"insane", 0, &config_insane, 1},
{"image", 1, NULL, 'e'},
{"basename", 0, &config_basename, 1},
{"record", 1, NULL, 5},
{0, 0, NULL, 0}
};
......@@ -529,6 +536,11 @@ int main(int argc, char *argv[])
element = strdup(optarg);
break;
case 5 :
stream_flags = parse_stream_flags(optarg);
if ( stream_flags < 0 ) return 1;
break;
case 0 :
break;
......@@ -673,6 +685,7 @@ int main(int argc, char *argv[])
free(pdb);
/* Start by writing the entire command line to stdout */
fprintf(ofh, "CrystFEL stream format 2.0\n");
fprintf(ofh, "Command line:");
for ( i=0; i<argc; i++ ) {
fprintf(ofh, " %s", argv[i]);
......
......@@ -472,8 +472,8 @@ int main(int argc, char *argv[])
image.filename = NULL;
image.features = NULL;
image.flags = NULL;
image.f0 = 1.0;
image.f0_available = 1;
image.i0 = 1.0;
image.i0_available = 1;
powder = calloc(image.width*image.height, sizeof(*powder));
......
......@@ -44,6 +44,7 @@ struct _refldata {
/* Intensity */
double intensity;
double esd_i;
};
......
......@@ -21,6 +21,63 @@
#include "cell.h"
#include "utils.h"
#include "image.h"
#include "stream.h"
static void exclusive(const char *a, const char *b)
{
ERROR("The stream options '%s' and '%s' are mutually exclusive.\n",
a, b);
}
int parse_stream_flags(const char *a)
{
int n, i;
int ret = STREAM_NONE;
char **flags;
n = assplode(a, ",", &flags, ASSPLODE_NONE);
for ( i=0; i<n; i++ ) {
if ( strcmp(flags[i], "pixels") == 0) {
if ( ret & STREAM_INTEGRATED ) {
exclusive("pixels", "integrated");
return -1;
}
ret |= STREAM_PIXELS;
} else if ( strcmp(flags[i], "integrated") == 0) {
if ( ret & STREAM_PIXELS ) {
exclusive("pixels", "integrated");
return -1;
}
ret |= STREAM_INTEGRATED;
} else if ( strcmp(flags[i], "peaks") == 0) {
if ( ret & STREAM_PEAKS_IF_INDEXED ) {
exclusive("peaks", "peaksifindexed");
return -1;
}
ret |= STREAM_PEAKS;
} else if ( strcmp(flags[i], "peaksifindexed") == 0) {
if ( ret & STREAM_PEAKS ) {
exclusive("peaks", "peaksifindexed");
return -1;
}
ret |= STREAM_PEAKS_IF_INDEXED;
} else {
ERROR("Unrecognised stream flag '%s'\n", flags[i]);
return 0;
}
free(flags[i]);
}
free(flags);
return ret;
}
int count_patterns(FILE *fh)
......@@ -76,6 +133,31 @@ static UnitCell *read_orientation_matrix(FILE *fh)
static void write_reflections(struct image *image, FILE *ofh)
{
Reflection *refl;
RefListIterator *iter;
fprintf(ofh, "Reflections measured after indexing\n");
/* FIXME: Unify this with write_reflections() over in reflections.c */
fprintf(ofh, " h k l I phase sigma(I) "
" 1/d(nm^-1) counts\n");
for ( refl = first_refl(image->reflections, &iter);
refl != NULL;
refl = next_refl(refl, iter) ) {
signed int h, k, l;
double intensity, esd_i, s;
get_indices(refl, &h, &k, &l);
intensity = get_intensity(refl);
esd_i = 0.0; /* FIXME! */
s = 0.0; /* FIXME! */
/* h, k, l, I, sigma(I), s */
fprintf(ofh, "%3i %3i %3i %10.2f %s %10.2f %10.2f %7i\n",
h, k, l, intensity, " -", esd_i, s/1.0e9, 1);
}
}
......@@ -83,7 +165,7 @@ static void write_peaks(struct image *image, FILE *ofh)
{
int i;
fprintf(ofh, "Peaks from peak search in %s\n", image->filename);
fprintf(ofh, "Peaks from peak search\n");
fprintf(ofh, " x/px y/px (1/d)/nm^-1 Intensity\n");
for ( i=0; i<image_feature_count(image->features); i++ ) {
......@@ -102,13 +184,10 @@ static void write_peaks(struct image *image, FILE *ofh)
f->x, f->y, q/1.0e9, f->intensity);
}
fprintf(ofh, "\n");
}
void write_chunk(FILE *ofh, struct image *image, int flags)
void write_chunk(FILE *ofh, struct image *i, int f)
{
double asx, asy, asz;
double bsx, bsy, bsz;
......@@ -117,36 +196,51 @@ void write_chunk(FILE *ofh, struct image *image, int flags)
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);
fprintf(ofh, "Image filename: %s\n", i->filename);
if ( i->indexed_cell != NULL ) {
cell_get_parameters(i->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(i->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);
} else {
fprintf(ofh, "No unit cell from indexing.\n");
}
if ( i->i0_available ) {
fprintf(ofh, "I0 = %7.5f (arbitrary units)\n", i->i0);
} else {
fprintf(ofh, "I0 = invalid\n");
}
fprintf(ofh, "photon_energy_eV = %f\n",
J_to_eV(ph_lambda_to_en(image->lambda)));
J_to_eV(ph_lambda_to_en(i->lambda)));
if ( (f & STREAM_PEAKS)
|| ((f & STREAM_PEAKS_IF_INDEXED) && (i->indexed_cell != NULL)) ) {
fprintf(ofh, "\n");
write_peaks(i, ofh);
}
write_peaks(image, ofh);
write_reflections(image, ofh);
if ( (f & STREAM_PIXELS) || (f & STREAM_INTEGRATED) ) {
fprintf(ofh, "\n");
write_reflections(i, ofh);
}
fprintf(ofh, "----- End chunk -----\n\n");
}
......
......@@ -19,9 +19,23 @@
struct image;
/* Possible options dictating what goes into the output stream */
enum
{
STREAM_NONE = 0,
STREAM_INTEGRATED = 1<<0,
STREAM_PIXELS = 1<<1,
STREAM_PEAKS = 1<<2,
STREAM_PEAKS_IF_INDEXED = 1<<3,
};
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);
extern int parse_stream_flags(const char *a);
#endif /* STREAM_H */
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment