aboutsummaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/hdf5-file.c2
-rw-r--r--src/image.h2
-rw-r--r--src/index.c2
-rw-r--r--src/indexamajig.c91
-rw-r--r--src/pattern_sim.c20
-rw-r--r--src/peaks.c134
-rw-r--r--src/peaks.h13
-rw-r--r--src/reintegrate.c11
-rw-r--r--src/stream.c79
-rw-r--r--src/stream.h3
10 files changed, 165 insertions, 192 deletions
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 <taw@physics.org>
+ * (c) 2006-2011 Thomas White <taw@physics.org>
*
* 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=<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;
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; 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;
}
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; 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);
+ } 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 */
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-15 21:40:19 +0000