Skip to content
GitLab
Commit 38ab00cc authored by Thomas White's avatar Thomas White
Browse files

First part of stream rework

parent 5c198c0e
Hide whitespace changes
Inline Side-by-side
Makefile.am
View file @ 38ab00cc
......@@ -44,7 +44,8 @@ endif
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/thread-pool.c src/image.c src/detector.c \
src/hdf5-file.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 \
......@@ -52,7 +53,7 @@ src_indexamajig_SOURCES = src/indexamajig.c src/hdf5-file.c src/utils.c \
src/sfac.c src/dirax.c src/mosflm.c \
src/reflections.c src/symmetry.c \
src/geometry.c src/thread-pool.c \
src/beam-parameters.c src/reflist.c
src/beam-parameters.c src/reflist.c src/stream.c
if HAVE_OPENCL
src_indexamajig_SOURCES += src/diffraction-gpu.c src/cl-utils.c
endif
......
Makefile.in
View file @ 38ab00cc
......@@ -150,7 +150,8 @@ am__src_indexamajig_SOURCES_DIST = src/indexamajig.c src/hdf5-file.c \
src/filters.c src/diffraction.c src/detector.c src/sfac.c \
src/dirax.c src/mosflm.c src/reflections.c src/symmetry.c \
src/geometry.c src/thread-pool.c src/beam-parameters.c \
src/reflist.c src/diffraction-gpu.c src/cl-utils.c
src/reflist.c src/stream.c src/diffraction-gpu.c \
src/cl-utils.c
@HAVE_OPENCL_TRUE@am__objects_1 = src/diffraction-gpu.$(OBJEXT) \
@HAVE_OPENCL_TRUE@ src/cl-utils.$(OBJEXT)
am_src_indexamajig_OBJECTS = src/indexamajig.$(OBJEXT) \
......@@ -161,7 +162,7 @@ am_src_indexamajig_OBJECTS = src/indexamajig.$(OBJEXT) \
src/mosflm.$(OBJEXT) src/reflections.$(OBJEXT) \
src/symmetry.$(OBJEXT) src/geometry.$(OBJEXT) \
src/thread-pool.$(OBJEXT) src/beam-parameters.$(OBJEXT) \
src/reflist.$(OBJEXT) $(am__objects_1)
src/reflist.$(OBJEXT) src/stream.$(OBJEXT) $(am__objects_1)
src_indexamajig_OBJECTS = $(am_src_indexamajig_OBJECTS)
src_indexamajig_LDADD = $(LDADD)
src_indexamajig_DEPENDENCIES = $(top_builddir)/lib/libgnu.a
......@@ -202,7 +203,9 @@ am_src_process_hkl_OBJECTS = src/process_hkl.$(OBJEXT) \
src/sfac.$(OBJEXT) src/statistics.$(OBJEXT) src/cell.$(OBJEXT) \
src/utils.$(OBJEXT) src/reflections.$(OBJEXT) \
src/symmetry.$(OBJEXT) src/stream.$(OBJEXT) \
src/beam-parameters.$(OBJEXT) src/thread-pool.$(OBJEXT)
src/beam-parameters.$(OBJEXT) src/thread-pool.$(OBJEXT) \
src/image.$(OBJEXT) src/detector.$(OBJEXT) \
src/hdf5-file.$(OBJEXT)
src_process_hkl_OBJECTS = $(am_src_process_hkl_OBJECTS)
src_process_hkl_LDADD = $(LDADD)
src_process_hkl_DEPENDENCIES = $(top_builddir)/lib/libgnu.a
......@@ -618,14 +621,15 @@ src_pattern_sim_SOURCES = src/pattern_sim.c src/diffraction.c \
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/thread-pool.c src/image.c src/detector.c \
src/hdf5-file.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 \
src/filters.c src/diffraction.c src/detector.c src/sfac.c \
src/dirax.c src/mosflm.c src/reflections.c src/symmetry.c \
src/geometry.c src/thread-pool.c src/beam-parameters.c \
src/reflist.c $(am__append_6)
src/reflist.c src/stream.c $(am__append_6)
@BUILD_HDFSEE_TRUE@src_hdfsee_SOURCES = src/hdfsee.c src/dw-hdfsee.c src/render.c \
@BUILD_HDFSEE_TRUE@ src/hdf5-file.c src/utils.c src/image.c src/filters.c \
@BUILD_HDFSEE_TRUE@ src/thread-pool.c src/detector.c
......
src/hdf5-file.c
View file @ 38ab00cc
......@@ -3,7 +3,7 @@
*
* Read/write HDF5 data files
*
* (c) 2006-2010 Thomas White <taw@physics.org>
* (c) 2006-2011 Thomas White <taw@physics.org>
*
* Part of CrystFEL - crystallography with a FEL
*
......
src/image.h
View file @ 38ab00cc
......@@ -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 */
......
src/index.c
View file @ 38ab00cc
......@@ -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 ) {
......
src/indexamajig.c
View file @ 38ab00cc
......@@ -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=<filename> Specify file containing list of images to process.\n"
" '-' means stdin, which is the default.\n"
" -o, --output=<filename> Write indexed stream to this file. '-' for stdout.\n"
" -o, --output=<filename> Write output stream to this file. '-' for stdout.\n"
" Default: indexamajig.stream\n"
"\n"
" --indexing=<methods> 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=<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"
"\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"
" 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;
......
src/pattern_sim.c
View file @ 38ab00cc
......@@ -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);
}
......
src/peaks.c
View file @ 38ab00cc
......@@ -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; i<image_feature_count(image->features); 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; fs<image->width; 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;
}
src/peaks.h
View file @ 38ab00cc
......@@ -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,
......
src/reintegrate.c
View file @ 38ab00cc
......@@ -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);
......
src/stream.c
View file @ 38ab00cc
......@@ -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; i<image_feature_count(image->features); 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);