aboutsummaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/dw-hdfsee.c1
-rw-r--r--src/hdfsee.c6
-rw-r--r--src/hrs-scaling.c77
-rw-r--r--src/hrs-scaling.h5
-rw-r--r--src/im-sandbox.c594
-rw-r--r--src/im-sandbox.h66
-rw-r--r--src/indexamajig.c383
-rw-r--r--src/partial_sim.c69
-rw-r--r--src/partialator.c235
-rw-r--r--src/post-refinement.c64
-rw-r--r--src/post-refinement.h5
-rw-r--r--src/process_hkl.c208
-rw-r--r--src/process_image.c231
-rw-r--r--src/process_image.h94
-rw-r--r--src/scaling-report.c83
-rw-r--r--src/scaling-report.h10
16 files changed, 1173 insertions, 958 deletions
diff --git a/src/dw-hdfsee.c b/src/dw-hdfsee.c
index 6ef12e7a..c4eabaf2 100644
--- a/src/dw-hdfsee.c
+++ b/src/dw-hdfsee.c
@@ -1819,7 +1819,6 @@ DisplayWindow *displaywindow_open(const char *filename, const char *peaks,
dw->hdfile = hdfile_open(filename);
if ( dw->hdfile == NULL ) {
- ERROR("Couldn't open file '%s'\n", filename);
free(dw);
return NULL;
} else {
diff --git a/src/hdfsee.c b/src/hdfsee.c
index 620b684b..36472e46 100644
--- a/src/hdfsee.c
+++ b/src/hdfsee.c
@@ -266,11 +266,7 @@ int main(int argc, char *argv[])
ring_radii,
n_rings,
ring_size);
- if ( main_window_list[i] == NULL ) {
- ERROR("Couldn't open display window\n");
- } else {
- main_n_windows++;
- }
+ if ( main_window_list[i] != NULL ) main_n_windows++;
}
if ( main_n_windows == 0 ) return 0;
diff --git a/src/hrs-scaling.c b/src/hrs-scaling.c
index f86ae0bb..0498a532 100644
--- a/src/hrs-scaling.c
+++ b/src/hrs-scaling.c
@@ -57,7 +57,7 @@
struct scale_queue_args
{
RefList *reference;
- struct image *images;
+ Crystal **crystals;
int n_started;
double max_shift;
};
@@ -65,7 +65,7 @@ struct scale_queue_args
struct scale_worker_args
{
- struct image *image;
+ Crystal *crystal;
double shift;
RefList *reference;
};
@@ -79,7 +79,7 @@ static void *create_scale_job(void *vqargs)
wargs = malloc(sizeof(struct scale_worker_args));
wargs->reference = qargs->reference;
- wargs->image = &qargs->images[qargs->n_started++];
+ wargs->crystal = qargs->crystals[qargs->n_started++];
return wargs;
}
@@ -88,20 +88,21 @@ static void *create_scale_job(void *vqargs)
static void run_scale_job(void *vwargs, int cookie)
{
struct scale_worker_args *wargs = vwargs;
- struct image *image = wargs->image;
+ Crystal *cr = wargs->crystal;
RefList *reference = wargs->reference;
Reflection *refl;
RefListIterator *iter;
double num = 0.0;
double den = 0.0;
double corr;
+ const double osf = crystal_get_osf(cr);
- if ( image->pr_dud ) {
+ if ( crystal_get_user_flag(cr) ) {
wargs->shift = 0.0;
return;
}
- for ( refl = first_refl(image->reflections, &iter);
+ for ( refl = first_refl(crystal_get_reflections(cr), &iter);
refl != NULL;
refl = next_refl(refl, iter) )
{
@@ -128,8 +129,8 @@ static void run_scale_job(void *vwargs, int cookie)
Ihl = get_intensity(refl) / get_partiality(refl);
esd = get_esd_intensity(refl) / get_partiality(refl);
- num += Ih * (Ihl/image->osf) / pow(esd/image->osf, 2.0);
- den += pow(Ih, 2.0)/pow(esd/image->osf, 2.0);
+ num += Ih * (Ihl/osf) / pow(esd/osf, 2.0);
+ den += pow(Ih, 2.0)/pow(esd/osf, 2.0);
}
@@ -138,7 +139,7 @@ static void run_scale_job(void *vwargs, int cookie)
corr = num / den;
if ( !isnan(corr) && !isinf(corr) ) {
- image->osf *= corr;
+ crystal_set_osf(cr, osf*corr);
}
wargs->shift = fabs(corr-1.0);
@@ -155,7 +156,7 @@ static void finalise_scale_job(void *vqargs, void *vwargs)
}
-static double iterate_scale(struct image *images, int n, RefList *reference,
+static double iterate_scale(Crystal **crystals, int n, RefList *reference,
int n_threads)
{
struct scale_queue_args qargs;
@@ -164,7 +165,7 @@ static double iterate_scale(struct image *images, int n, RefList *reference,
qargs.reference = reference;
qargs.n_started = 0;
- qargs.images = images;
+ qargs.crystals = crystals;
qargs.max_shift = 0.0;
run_threads(n_threads, run_scale_job, create_scale_job,
@@ -178,14 +179,14 @@ struct merge_queue_args
{
RefList *full;
pthread_mutex_t full_lock;
- struct image *images;
+ Crystal **crystals;
int n_started;
};
struct merge_worker_args
{
- struct image *image;
+ Crystal *crystal;
RefList *full;
pthread_mutex_t *full_lock;
};
@@ -200,7 +201,7 @@ static void *create_merge_job(void *vqargs)
wargs->full = qargs->full;
wargs->full_lock = &qargs->full_lock;
- wargs->image = &qargs->images[qargs->n_started++];
+ wargs->crystal = qargs->crystals[qargs->n_started++];
return wargs;
}
@@ -209,17 +210,17 @@ static void *create_merge_job(void *vqargs)
static void run_merge_job(void *vwargs, int cookie)
{
struct merge_worker_args *wargs = vwargs;
- struct image *image = wargs->image;
+ Crystal *cr = wargs->crystal;
RefList *full = wargs->full;
Reflection *refl;
RefListIterator *iter;
double G;
- if ( image->pr_dud ) return;
+ if ( crystal_get_user_flag(cr)) return;
- G = image->osf;
+ G = crystal_get_osf(cr);
- for ( refl = first_refl(image->reflections, &iter);
+ for ( refl = first_refl(crystal_get_reflections(cr), &iter);
refl != NULL;
refl = next_refl(refl, iter) )
{
@@ -272,7 +273,7 @@ static void finalise_merge_job(void *vqargs, void *vwargs)
}
-static RefList *lsq_intensities(struct image *images, int n, int n_threads)
+static RefList *lsq_intensities(Crystal **crystals, int n, int n_threads)
{
RefList *full;
struct merge_queue_args qargs;
@@ -283,7 +284,7 @@ static RefList *lsq_intensities(struct image *images, int n, int n_threads)
qargs.full = full;
qargs.n_started = 0;
- qargs.images = images;
+ qargs.crystals = crystals;
pthread_mutex_init(&qargs.full_lock, NULL);
run_threads(n_threads, run_merge_job, create_merge_job,
@@ -308,14 +309,14 @@ static RefList *lsq_intensities(struct image *images, int n, int n_threads)
struct esd_queue_args
{
RefList *full;
- struct image *images;
+ Crystal **crystals;
int n_started;
};
struct esd_worker_args
{
- struct image *image;
+ Crystal *crystal;
RefList *full;
};
@@ -328,7 +329,7 @@ static void *create_esd_job(void *vqargs)
wargs = malloc(sizeof(struct esd_worker_args));
wargs->full = qargs->full;
- wargs->image = &qargs->images[qargs->n_started++];
+ wargs->crystal = qargs->crystals[qargs->n_started++];
return wargs;
}
@@ -337,17 +338,17 @@ static void *create_esd_job(void *vqargs)
static void run_esd_job(void *vwargs, int cookie)
{
struct esd_worker_args *wargs = vwargs;
- struct image *image = wargs->image;
+ Crystal *cr = wargs->crystal;
RefList *full = wargs->full;
Reflection *refl;
RefListIterator *iter;
double G;
- if ( image->pr_dud ) return;
+ if ( crystal_get_user_flag(cr) ) return;
- G = image->osf;
+ G = crystal_get_osf(cr);
- for ( refl = first_refl(image->reflections, &iter);
+ for ( refl = first_refl(crystal_get_reflections(cr), &iter);
refl != NULL;
refl = next_refl(refl, iter) )
{
@@ -383,7 +384,7 @@ static void finalise_esd_job(void *vqargs, void *vwargs)
}
-static void calculate_esds(struct image *images, int n, RefList *full,
+static void calculate_esds(Crystal **crystals, int n, RefList *full,
int n_threads, int min_red)
{
struct esd_queue_args qargs;
@@ -392,7 +393,7 @@ static void calculate_esds(struct image *images, int n, RefList *full,
qargs.full = full;
qargs.n_started = 0;
- qargs.images = images;
+ qargs.crystals = crystals;
for ( refl = first_refl(full, &iter);
refl != NULL;
@@ -423,7 +424,7 @@ static void calculate_esds(struct image *images, int n, RefList *full,
/* Scale the stack of images */
-RefList *scale_intensities(struct image *images, int n, RefList *gref,
+RefList *scale_intensities(Crystal **crystals, int n, RefList *gref,
int n_threads, int noscale)
{
int i;
@@ -431,17 +432,17 @@ RefList *scale_intensities(struct image *images, int n, RefList *gref,
RefList *full = NULL;
const int min_redundancy = 3;
- for ( i=0; i<n; i++ ) images[i].osf = 1.0;
+ for ( i=0; i<n; i++ ) crystal_set_osf(crystals[i], 1.0);
if ( noscale ) {
- full = lsq_intensities(images, n, n_threads);
- calculate_esds(images, n, full, n_threads, min_redundancy);
+ full = lsq_intensities(crystals, n, n_threads);
+ calculate_esds(crystals, n, full, n_threads, min_redundancy);
return full;
}
/* No reference -> create an initial list to refine against */
if ( gref == NULL ) {
- full = lsq_intensities(images, n, n_threads);
+ full = lsq_intensities(crystals, n, n_threads);
}
/* Iterate */
@@ -457,14 +458,14 @@ RefList *scale_intensities(struct image *images, int n, RefList *gref,
reference = full;
}
- max_corr = iterate_scale(images, n, reference, n_threads);
+ max_corr = iterate_scale(crystals, n, reference, n_threads);
//STATUS("Scaling iteration %2i: max correction = %5.2f\n",
// i+1, max_corr);
/* No reference -> generate list for next iteration */
if ( gref == NULL ) {
reflist_free(full);
- full = lsq_intensities(images, n, n_threads);
+ full = lsq_intensities(crystals, n, n_threads);
}
//show_scale_factors(images, n);
@@ -474,10 +475,10 @@ RefList *scale_intensities(struct image *images, int n, RefList *gref,
} while ( (max_corr > 0.01) && (i < MAX_CYCLES) );
if ( gref != NULL ) {
- full = lsq_intensities(images, n, n_threads);
+ full = lsq_intensities(crystals, n, n_threads);
} /* else we already did it */
- calculate_esds(images, n, full, n_threads, min_redundancy);
+ calculate_esds(crystals, n, full, n_threads, min_redundancy);
return full;
}
diff --git a/src/hrs-scaling.h b/src/hrs-scaling.h
index 5425b1ff..80940347 100644
--- a/src/hrs-scaling.h
+++ b/src/hrs-scaling.h
@@ -35,9 +35,10 @@
#endif
-#include "image.h"
+#include "crystal.h"
+#include "reflist.h"
-extern RefList *scale_intensities(struct image *images, int n,
+extern RefList *scale_intensities(Crystal **crystals, int n,
RefList *reference, int n_threads,
int noscale);
diff --git a/src/im-sandbox.c b/src/im-sandbox.c
index c8cceb2b..03c553b7 100644
--- a/src/im-sandbox.c
+++ b/src/im-sandbox.c
@@ -3,13 +3,13 @@
*
* Sandbox for indexing
*
- * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY,
- * a research centre of the Helmholtz Association.
+ * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY,
+ * a research centre of the Helmholtz Association.
* Copyright © 2012 Richard Kirian
* Copyright © 2012 Lorenzo Galli
*
* Authors:
- * 2010-2012 Thomas White <taw@physics.org>
+ * 2010-2013 Thomas White <taw@physics.org>
* 2011 Richard Kirian
* 2012 Lorenzo Galli
* 2012 Chunhong Yoon
@@ -41,12 +41,11 @@
#include <string.h>
#include <unistd.h>
#include <getopt.h>
-#include <hdf5.h>
-#include <gsl/gsl_errno.h>
#include <pthread.h>
#include <sys/wait.h>
#include <fcntl.h>
#include <signal.h>
+#include <sys/stat.h>
#ifdef HAVE_CLOCK_GETTIME
#include <time.h>
@@ -54,48 +53,56 @@
#include <sys/time.h>
#endif
-#include "utils.h"
-#include "hdf5-file.h"
-#include "index.h"
-#include "peaks.h"
-#include "detector.h"
-#include "filters.h"
-#include "thread-pool.h"
-#include "beam-parameters.h"
-#include "geometry.h"
-#include "stream.h"
-#include "reflist-utils.h"
-
#include "im-sandbox.h"
+#include "process_image.h"
/* Write statistics at APPROXIMATELY this interval */
#define STATS_EVERY_N_SECONDS (5)
+struct sb_reader
+{
+ pthread_mutex_t lock;
+ int done;
+
+ /* If a worker process dies unexpectedly (e.g. if it segfaults), then
+ * the pipe for its output can still stay open for a little while while
+ * its buffer empties. The number of pipes being read from is therefore
+ * not necessarily the same as the number of worker processes. */
+ int n_read;
+ FILE **fhs;
+ int *fds;
+
+ /* Final output fd */
+ int ofd;
+};
+
+
struct sandbox
{
pthread_mutex_t lock;
- int n_indexable;
int n_processed;
- int n_indexable_last_stats;
+ int n_hadcrystals;
+ int n_crystals;
int n_processed_last_stats;
+ int n_hadcrystals_last_stats;
+ int n_crystals_last_stats;
int t_last_stats;
struct index_args *iargs;
int n_proc;
pid_t *pids;
- FILE *ofh;
- FILE **fhs;
int *running;
FILE **result_fhs;
int *filename_pipes;
- int *stream_pipe_read;
int *stream_pipe_write;
char **last_filename;
+
+ struct sb_reader *reader;
};
@@ -120,6 +127,7 @@ static char *get_pattern(FILE *fh, char **use_this_one_instead,
{
char *line;
char *filename;
+ size_t len;
do {
@@ -153,7 +161,13 @@ static char *get_pattern(FILE *fh, char **use_this_one_instead,
line = tmp;
}
- filename = malloc(strlen(prefix)+strlen(line)+1);
+ len = strlen(prefix)+strlen(line)+1;
+
+ /* Round the length of the buffer, too keep Valgrind quiet when it gets
+ * given to write() a bit later on */
+ len += 4 - (len % 4);
+
+ filename = malloc(len);
snprintf(filename, 1023, "%s%s", prefix, line);
@@ -163,187 +177,8 @@ static char *get_pattern(FILE *fh, char **use_this_one_instead,
}
-static void process_image(const struct index_args *iargs,
- struct pattern_args *pargs, FILE *ofh,
- int cookie)
-{
- float *data_for_measurement;
- size_t data_size;
- UnitCell *cell = iargs->cell;
- int config_cmfilter = iargs->config_cmfilter;
- int config_noisefilter = iargs->config_noisefilter;
- int config_verbose = iargs->config_verbose;
- IndexingMethod *indm = iargs->indm;
- int check;
- struct hdfile *hdfile;
- struct image image;
-
- image.features = NULL;
- image.data = NULL;
- image.flags = NULL;
- image.indexed_cell = NULL;
- image.copyme = iargs->copyme;
- image.reflections = NULL;
- image.n_saturated = 0;
- image.id = cookie;
- image.filename = pargs->filename;
- image.beam = iargs->beam;
- image.det = iargs->det;
-
- hdfile = hdfile_open(image.filename);
- if ( hdfile == NULL ) return;
-
- if ( iargs->element != NULL ) {
-
- int r;
- r = hdfile_set_image(hdfile, iargs->element);
- if ( r ) {
- ERROR("Couldn't select path '%s'\n", iargs->element);
- hdfile_close(hdfile);
- return;
- }
-
- } else {
-
- int r;
- r = hdfile_set_first_image(hdfile, "/");
- if ( r ) {
- ERROR("Couldn't select first path\n");
- hdfile_close(hdfile);
- return;
- }
-
- }
-
- check = hdf5_read(hdfile, &image, 1);
- if ( check ) {
- hdfile_close(hdfile);
- return;
- }
-
- if ( (image.width != image.det->max_fs + 1 )
- || (image.height != image.det->max_ss + 1))
- {
- ERROR("Image size doesn't match geometry size"
- " - rejecting image.\n");
- ERROR("Image size: %i,%i. Geometry size: %i,%i\n",
- image.width, image.height,
- image.det->max_fs + 1, image.det->max_ss + 1);
- hdfile_close(hdfile);
- return;
- }
-
- fill_in_values(image.det, hdfile);
- fill_in_beam_parameters(image.beam, hdfile);
-
- image.lambda = ph_en_to_lambda(eV_to_J(image.beam->photon_energy));
-
- if ( (image.beam->photon_energy < 0.0) || (image.lambda > 1000) ) {
- /* Error message covers a silly value in the beam file or in
- * the HDF5 file. */
- ERROR("Nonsensical wavelength (%e m or %e eV) value for %s.\n",
- image.lambda, image.beam->photon_energy, image.filename);
- hdfile_close(hdfile);
- return;
- }
-
- if ( config_cmfilter ) {
- filter_cm(&image);
- }
-
- /* Take snapshot of image after CM subtraction but before
- * the aggressive noise filter. */
- data_size = image.width * image.height * sizeof(float);
- data_for_measurement = malloc(data_size);
-
- if ( config_noisefilter ) {
- filter_noise(&image, data_for_measurement);
- } else {
- memcpy(data_for_measurement, image.data, data_size);
- }
-
- switch ( iargs->peaks ) {
-
- case PEAK_HDF5:
- // Get peaks from HDF5
- if (get_peaks(&image, hdfile,
- iargs->hdf5_peak_path)) {
- ERROR("Failed to get peaks from HDF5 file.\n");
- }
- if ( !iargs->no_revalidate ) {
- validate_peaks(&image, iargs->min_int_snr,
- iargs->ir_inn, iargs->ir_mid,
- iargs->ir_out, iargs->use_saturated);
- }
- break;
-
- case PEAK_ZAEF:
- search_peaks(&image, iargs->threshold,
- iargs->min_gradient, iargs->min_snr,
- iargs->ir_inn, iargs->ir_mid, iargs->ir_out,
- iargs->use_saturated);
- break;
-
- }
-
- /* Get rid of noise-filtered version at this point
- * - it was strictly for the purposes of peak detection. */
- free(image.data);
- image.data = data_for_measurement;
-
- /* Calculate orientation matrix (by magic) */
- image.div = image.beam->divergence;
- image.bw = image.beam->bandwidth;
- image.profile_radius = image.beam->profile_radius;
-
- index_pattern(&image, cell, indm, iargs->cellr,
- config_verbose, iargs->ipriv,
- iargs->config_insane, iargs->tols);
-
- if ( image.indexed_cell != NULL ) {
-
- pargs->indexable = 1;
-
- if ( iargs->integrate_found ) {
- image.reflections = select_intersections(&image,
- image.indexed_cell);
- } else {
- image.reflections = find_intersections(&image,
- image.indexed_cell);
- }
-
- if (image.reflections != NULL) {
- integrate_reflections(&image,
- iargs->config_closer,
- iargs->config_bgsub,
- iargs->min_int_snr,
- iargs->ir_inn,
- iargs->ir_mid,
- iargs->ir_out,
- iargs->integrate_saturated);
- }
-
- } else {
- image.reflections = NULL;
- }
-
- write_chunk(ofh, &image, hdfile, iargs->stream_flags);
- fprintf(ofh, "END\n");
- fflush(ofh);
-
- /* Only free cell if found */
- cell_free(image.indexed_cell);
-
- reflist_free(image.reflections);
- free(image.data);
- if ( image.flags != NULL ) free(image.flags);
- image_feature_list_free(image.features);
- hdfile_close(hdfile);
-}
-
-
static void run_work(const struct index_args *iargs,
- int filename_pipe, int results_pipe, FILE *ofh,
+ int filename_pipe, int results_pipe, Stream *st,
int cookie)
{
int allDone = 0;
@@ -389,12 +224,12 @@ static void run_work(const struct index_args *iargs,
} else {
pargs.filename = line;
- pargs.indexable = 0;
+ pargs.n_crystals = 0;
- process_image(iargs, &pargs, ofh, cookie);
+ process_image(iargs, &pargs, st, cookie);
/* Request another image */
- c = sprintf(buf, "%i\n", pargs.indexable);
+ c = sprintf(buf, "%i\n", pargs.n_crystals);
w = write(results_pipe, buf, c);
if ( w < 0 ) {
ERROR("write P0\n");
@@ -406,7 +241,7 @@ static void run_work(const struct index_args *iargs,
}
- cleanup_indexing(iargs->ipriv);
+ cleanup_indexing(iargs->indm, iargs->ipriv);
free(iargs->indm);
free(iargs->ipriv);
free_detector_geometry(iargs->det);
@@ -442,11 +277,19 @@ static time_t get_monotonic_seconds()
#endif
+size_t vol = 0;
+
+
+static ssize_t lwrite(int fd, const char *a)
+{
+ size_t l = strlen(a);
+ return write(fd, a, l);
+}
+
-static int pump_chunk(FILE *fh, FILE *ofh)
+static int pump_chunk(FILE *fh, int ofd)
{
int chunk_started = 0;
- int chunk_finished = 0;
do {
@@ -457,65 +300,121 @@ static int pump_chunk(FILE *fh, FILE *ofh)
if ( rval == NULL ) {
if ( feof(fh) ) {
- /* Process died */
+ /* Whoops, connection lost */
if ( chunk_started ) {
ERROR("EOF during chunk!\n");
- fprintf(ofh, "Chunk is unfinished!\n");
- }
+ lwrite(ofd, "Unfinished chunk!\n");
+ lwrite(ofd, CHUNK_END_MARKER"\n");
+ } /* else normal end of output */
return 1;
- } else {
- ERROR("fgets() failed: %s\n", strerror(errno));
}
- chunk_finished = 1;
- continue;
- }
+ ERROR("fgets() failed: %s\n", strerror(errno));
+ return 1;
- if ( strcmp(line, "END\n") == 0 ) {
- chunk_finished = 1;
- } else {
- chunk_started = 1;
- fprintf(ofh, "%s", line);
}
- } while ( !chunk_finished );
+ if ( strcmp(line, "FLUSH\n") == 0 ) break;
+ lwrite(ofd, line);
+
+ if ( strcmp(line, CHUNK_END_MARKER"\n") == 0 ) break;
+ if ( strcmp(line, CHUNK_START_MARKER"\n") == 0 ) break;
+ } while ( 1 );
return 0;
}
-static void *run_reader(void *sbv)
+/* Add an fd to the list of pipes to be read from */
+static void add_pipe(struct sb_reader *rd, int fd)
{
- struct sandbox *sb = sbv;
- int done = 0;
+ int *fds_new;
+ FILE **fhs_new;
+ int slot;
- while ( !done ) {
+ pthread_mutex_lock(&rd->lock);
+
+ fds_new = realloc(rd->fds, (rd->n_read+1)*sizeof(int));
+ if ( fds_new == NULL ) {
+ ERROR("Failed to allocate memory for new pipe.\n");
+ return;
+ }
+
+ fhs_new = realloc(rd->fhs, (rd->n_read+1)*sizeof(FILE *));
+ if ( fhs_new == NULL ) {
+ ERROR("Failed to allocate memory for new FH.\n");
+ return;
+ }
+
+ rd->fds = fds_new;
+ rd->fhs = fhs_new;
+ slot = rd->n_read;
+
+ rd->fds[slot] = fd;
+
+ rd->fhs[slot] = fdopen(fd, "r");
+ if ( rd->fhs[slot] == NULL ) {
+ ERROR("Couldn't fdopen() stream!\n");
+ return;
+ }
+
+ rd->n_read++;
+
+ pthread_mutex_unlock(&rd->lock);
+}
+
+
+/* Assumes that the caller is already holding rd->lock! */
+static void remove_pipe(struct sb_reader *rd, int d)
+{
+ int i;
+
+ for ( i=d; i<rd->n_read; i++ ) {
+ if ( i < rd->n_read-1 ) {
+ rd->fds[i] = rd->fds[i+1];
+ rd->fhs[i] = rd->fhs[i+1];
+ } /* else don't bother */
+ }
+
+ rd->n_read--;
+
+ /* We don't bother shrinking the arrays */
+}
+
+
+static void *run_reader(void *rdv)
+{
+ struct sb_reader *rd = rdv;
+
+ while ( 1 ) {
int r, i;
struct timeval tv;
fd_set fds;
int fdmax;
- tv.tv_sec = 5;
+ /* Exit when:
+ * - No fhs left open to read from
+ * AND - Main thread says "done" */
+ if ( (rd->n_read == 0) && rd->done ) break;
+
+ tv.tv_sec = 1;
tv.tv_usec = 0;
FD_ZERO(&fds);
fdmax = 0;
- lock_sandbox(sb);
- for ( i=0; i<sb->n_proc; i++ ) {
+ pthread_mutex_lock(&rd->lock);
+ for ( i=0; i<rd->n_read; i++ ) {
int fd;
- if ( !sb->running[i] ) continue;
-
- fd = sb->stream_pipe_read[i];
+ fd = rd->fds[i];
FD_SET(fd, &fds);
if ( fd > fdmax ) fdmax = fd;
}
-
- unlock_sandbox(sb);
+ pthread_mutex_unlock(&rd->lock);
r = select(fdmax+1, &fds, NULL, NULL, &tv);
@@ -526,38 +425,75 @@ static void *run_reader(void *sbv)
continue;
}
- if ( r == 0 ) continue; /* Nothing this time. Try again */
+ pthread_mutex_lock(&rd->lock);
+ for ( i=0; i<rd->n_read; i++ ) {
- lock_sandbox(sb);
- for ( i=0; i<sb->n_proc; i++ ) {
+ if ( !FD_ISSET(rd->fds[i], &fds) ) {
+ continue;
+ }
- if ( !sb->running[i] ) continue;
+ /* If the chunk cannot be read, assume the connection
+ * is broken and that the process will die soon. */
+ if ( pump_chunk(rd->fhs[i], rd->ofd) ) {
+ /* remove_pipe() assumes that the caller is
+ * holding rd->lock ! */
+ remove_pipe(rd, i);
+ }
- if ( !FD_ISSET(sb->stream_pipe_read[i], &fds) ) continue;
+ }
+ pthread_mutex_unlock(&rd->lock);
+
+ }
+
+ return NULL;
+}
- if ( pump_chunk(sb->fhs[i], sb->ofh) ) {
- sb->running[i] = 0;
- }
+static int create_temporary_folder(signed int n)
+{
+ int r;
+ char tmp[64];
+ struct stat s;
+
+ if ( n < 0 ) {
+ snprintf(tmp, 63, "indexamajig.%i", getpid());
+ } else {
+ snprintf(tmp, 63, "worker.%i", n);
+ }
+
+ if ( stat(tmp, &s) == -1 ) {
+ if ( errno != ENOENT ) {
+ ERROR("Failed to stat temporary folder.\n");
+ return 1;
}
- done = 1;
- for ( i=0; i<sb->n_proc; i++ ) {
- if ( sb->running[i] ) done = 0;
+ r = mkdir(tmp, S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
+ if ( r ) {
+ ERROR("Failed to create temporary folder: %s\n",
+ strerror(errno));
+ return 1;
}
- unlock_sandbox(sb);
}
- return NULL;
+ r = chdir(tmp);
+ if ( r ) {
+ ERROR("Failed to chdir to temporary folder: %s\n",
+ strerror(errno));
+ return 1;
+ }
+
+ return 0;
}
-static void start_worker_process(struct sandbox *sb, int slot)
+static void start_worker_process(struct sandbox *sb, int slot,
+ int argc, char *argv[])
{
pid_t p;
int filename_pipe[2];
int result_pipe[2];
+ int stream_pipe[2];
if ( pipe(filename_pipe) == - 1 ) {
ERROR("pipe() failed!\n");
@@ -569,6 +505,11 @@ static void start_worker_process(struct sandbox *sb, int slot)
return;
}
+ if ( pipe(stream_pipe) == - 1 ) {
+ ERROR("pipe() failed!\n");
+ return;
+ }
+
p = fork();
if ( p == -1 ) {
ERROR("fork() failed!\n");
@@ -577,7 +518,7 @@ static void start_worker_process(struct sandbox *sb, int slot)
if ( p == 0 ) {
- FILE *sfh;
+ Stream *st;
int j;
struct sigaction sa;
int r;
@@ -592,6 +533,8 @@ static void start_worker_process(struct sandbox *sb, int slot)
return;
}
+ create_temporary_folder(slot);
+
/* Free resources which will not be needed by worker */
for ( j=0; j<sb->n_proc; j++ ) {
if ( (j != slot) && (sb->running[j]) ) {
@@ -600,7 +543,9 @@ static void start_worker_process(struct sandbox *sb, int slot)
}
for ( j=0; j<sb->n_proc; j++ ) {
if ( (j != slot) && (sb->running[j]) ) {
- fclose(sb->result_fhs[j]);
+ if ( sb->result_fhs[j] != NULL ) {
+ fclose(sb->result_fhs[j]);
+ }
close(sb->filename_pipes[j]);
}
}
@@ -614,10 +559,12 @@ static void start_worker_process(struct sandbox *sb, int slot)
close(filename_pipe[1]);
close(result_pipe[0]);
- sfh = fdopen(sb->stream_pipe_write[slot], "w");
+ st = open_stream_fd_for_write(stream_pipe[1]);
+ write_command(st, argc, argv);
+ write_line(st, "FLUSH");
run_work(sb->iargs, filename_pipe[0], result_pipe[1],
- sfh, slot);
- fclose(sfh);
+ st, slot);
+ close_stream(st);
//close(filename_pipe[0]);
close(result_pipe[1]);
@@ -630,14 +577,11 @@ static void start_worker_process(struct sandbox *sb, int slot)
* and the 'read' end of the result pipe. */
sb->pids[slot] = p;
sb->running[slot] = 1;
+ add_pipe(sb->reader, stream_pipe[0]);
close(filename_pipe[0]);
close(result_pipe[1]);
+ close(stream_pipe[1]);
sb->filename_pipes[slot] = filename_pipe[1];
- sb->fhs[slot] = fdopen(sb->stream_pipe_read[slot], "r");
- if ( sb->fhs[slot] == NULL ) {
- ERROR("Couldn't fdopen() stream!\n");
- return;
- }
sb->result_fhs[slot] = fdopen(result_pipe[0], "r");
if ( sb->result_fhs[slot] == NULL ) {
@@ -685,7 +629,7 @@ static void handle_zombie(struct sandbox *sb)
STATUS("Last filename was: %s\n",
sb->last_filename[i]);
sb->n_processed++;
- start_worker_process(sb, i);
+ start_worker_process(sb, i, 0, NULL);
}
}
@@ -697,7 +641,7 @@ static void handle_zombie(struct sandbox *sb)
void create_sandbox(struct index_args *iargs, int n_proc, char *prefix,
int config_basename, FILE *fh, char *use_this_one_instead,
- FILE *ofh)
+ int ofd, int argc, char *argv[])
{
int i;
int allDone;
@@ -712,48 +656,41 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix,
return;
}
- sb->n_indexable = 0;
+ sb->reader = calloc(1, sizeof(struct sb_reader));
+ if ( sb->reader == NULL ) {
+ ERROR("Couldn't allocate memory for SB reader.\n");
+ free(sb);
+ return;
+ }
+
+ pthread_mutex_init(&sb->lock, NULL);
+ pthread_mutex_init(&sb->reader->lock, NULL);
+
sb->n_processed = 0;
- sb->n_indexable_last_stats = 0;
+ sb->n_hadcrystals = 0;
+ sb->n_crystals = 0;
sb->n_processed_last_stats = 0;
+ sb->n_hadcrystals_last_stats = 0;
+ sb->n_crystals_last_stats = 0;
sb->t_last_stats = get_monotonic_seconds();
sb->n_proc = n_proc;
- sb->ofh = ofh;
sb->iargs = iargs;
- pthread_mutex_init(&sb->lock, NULL);
+ sb->reader->fds = NULL;
+ sb->reader->fhs = NULL;
+ sb->reader->ofd = ofd;
- sb->stream_pipe_read = calloc(n_proc, sizeof(int));
sb->stream_pipe_write = calloc(n_proc, sizeof(int));
- if ( sb->stream_pipe_read == NULL ) {
- ERROR("Couldn't allocate memory for pipes.\n");
- return;
- }
if ( sb->stream_pipe_write == NULL ) {
ERROR("Couldn't allocate memory for pipes.\n");
return;
}
- for ( i=0; i<n_proc; i++ ) {
-
- int stream_pipe[2];
-
- if ( pipe(stream_pipe) == - 1 ) {
- ERROR("pipe() failed!\n");
- return;
- }
-
- sb->stream_pipe_read[i] = stream_pipe[0];
- sb->stream_pipe_write[i] = stream_pipe[1];
-
- }
-
lock_sandbox(sb);
sb->filename_pipes = calloc(n_proc, sizeof(int));
sb->result_fhs = calloc(n_proc, sizeof(FILE *));
sb->pids = calloc(n_proc, sizeof(pid_t));
sb->running = calloc(n_proc, sizeof(int));
- sb->fhs = calloc(sb->n_proc, sizeof(FILE *));
if ( sb->filename_pipes == NULL ) {
ERROR("Couldn't allocate memory for pipes.\n");
return;
@@ -776,17 +713,8 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix,
ERROR("Couldn't allocate memory for last filename list.\n");
return;
}
- if ( sb->fhs == NULL ) {
- ERROR("Couldn't allocate memory for file handles!\n");
- return;
- }
unlock_sandbox(sb);
- if ( pthread_create(&reader_thread, NULL, run_reader, (void *)sb) ) {
- ERROR("Failed to create reader thread.\n");
- return;
- }
-
if ( pipe(signal_pipe) == -1 ) {
ERROR("Failed to create signal pipe.\n");
return;
@@ -802,13 +730,23 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix,
return;
}
+ if ( create_temporary_folder(-1) ) return;
+
/* Fork the right number of times */
lock_sandbox(sb);
for ( i=0; i<n_proc; i++ ) {
- start_worker_process(sb, i);
+ start_worker_process(sb, i, argc, argv);
}
unlock_sandbox(sb);
+ /* Start reader thread after forking, so that things are definitely
+ * "running" */
+ if ( pthread_create(&reader_thread, NULL, run_reader,
+ (void *)sb->reader) ) {
+ ERROR("Failed to create reader thread.\n");
+ return;
+ }
+
allDone = 0;
while ( !allDone ) {
@@ -828,9 +766,7 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix,
int fd;
- if ( !sb->running[i] ) {
- continue;
- }
+ if ( sb->result_fhs[i] == NULL) continue;
fd = fileno(sb->result_fhs[i]);
FD_SET(fd, &fds);
@@ -844,14 +780,11 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix,
r = select(fdmax+1, &fds, NULL, NULL, &tv);
if ( r == -1 ) {
- if ( errno != EINTR ) {
- ERROR("select() failed: %s\n", strerror(errno));
- }
- continue;
+ if ( errno == EINTR ) continue;
+ ERROR("select() failed: %s\n", strerror(errno));
+ break;
}
- if ( r == 0 ) continue; /* No progress this time. Try again */
-
if ( FD_ISSET(signal_pipe[0], &fds) ) {
char d;
@@ -869,14 +802,10 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix,
int fd;
char *eptr;
- if ( !sb->running[i] ) {
- continue;
- }
+ if ( sb->result_fhs[i] == NULL ) continue;
fd = fileno(sb->result_fhs[i]);
- if ( !FD_ISSET(fd, &fds) ) {
- continue;
- }
+ if ( !FD_ISSET(fd, &fds) ) continue;
rval = fgets(results, 1024, sb->result_fhs[i]);
if ( rval == NULL ) {
@@ -884,6 +813,7 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix,
ERROR("fgets() failed: %s\n",
strerror(errno));
}
+ sb->result_fhs[i] = NULL;
continue;
}
@@ -895,8 +825,14 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix,
ERROR("Invalid result '%s'\n", results);
}
} else {
- sb->n_indexable += atoi(results);
+
+ int nc = atoi(results);
+ sb->n_crystals += nc;
+ if ( nc > 0 ) {
+ sb->n_hadcrystals++;
+ }
sb->n_processed++;
+
}
/* Send next filename */
@@ -929,14 +865,16 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix,
tNow = get_monotonic_seconds();
if ( tNow >= sb->t_last_stats+STATS_EVERY_N_SECONDS ) {
- STATUS("%i out of %i indexed so far,"
- " %i out of %i since the last message.\n",
- sb->n_indexable, sb->n_processed,
- sb->n_indexable - sb->n_indexable_last_stats,
+ STATUS("%4i indexable out of %4i processed, "
+ "%4i crystals so far. "
+ "%4i images processed since the last message.\n",
+ sb->n_hadcrystals, sb->n_processed,
+ sb->n_crystals,
sb->n_processed - sb->n_processed_last_stats);
- sb->n_indexable_last_stats = sb->n_indexable;
sb->n_processed_last_stats = sb->n_processed;
+ sb->n_hadcrystals_last_stats = sb->n_hadcrystals;
+ sb->n_crystals_last_stats = sb->n_crystals;
sb->t_last_stats = tNow;
}
@@ -947,14 +885,18 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix,
for ( i=0; i<n_proc; i++ ) {
if ( sb->running[i] ) allDone = 0;
}
-
unlock_sandbox(sb);
}
fclose(fh);
- pthread_mutex_destroy(&sb->lock);
+ /* Indicate to the reader thread that we are done */
+ pthread_mutex_lock(&sb->reader->lock);
+ sb->reader->done = 1;
+ pthread_mutex_unlock(&sb->reader->lock);
+
+ pthread_join(reader_thread, NULL);
for ( i=0; i<n_proc; i++ ) {
int status;
@@ -963,21 +905,19 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix,
for ( i=0; i<n_proc; i++ ) {
close(sb->filename_pipes[i]);
- fclose(sb->result_fhs[i]);
+ if ( sb->result_fhs[i] != NULL ) fclose(sb->result_fhs[i]);
}
- for ( i=0; i<sb->n_proc; i++ ) {
- fclose(sb->fhs[i]);
- }
- free(sb->fhs);
+ free(sb->running);
free(sb->filename_pipes);
free(sb->result_fhs);
free(sb->pids);
- free(sb->running);
- if ( ofh != stdout ) fclose(ofh);
+ pthread_mutex_destroy(&sb->lock);
- STATUS("There were %i images, of which %i could be indexed.\n",
- sb->n_processed, sb->n_indexable);
+ STATUS("Final:"
+ " %i images processed, %i had crystals, %i crystals overall.\n",
+ sb->n_processed, sb->n_hadcrystals, sb->n_crystals);
+ free(sb);
}
diff --git a/src/im-sandbox.h b/src/im-sandbox.h
index 872f84ad..9248b226 100644
--- a/src/im-sandbox.h
+++ b/src/im-sandbox.h
@@ -3,13 +3,13 @@
*
* Sandbox for indexing
*
- * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY,
- * a research centre of the Helmholtz Association.
+ * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY,
+ * a research centre of the Helmholtz Association.
* Copyright © 2012 Richard Kirian
* Copyright © 2012 Lorenzo Galli
*
* Authors:
- * 2010-2012 Thomas White <taw@physics.org>
+ * 2010-2013 Thomas White <taw@physics.org>
* 2011 Richard Kirian
* 2012 Lorenzo Galli
* 2012 Chunhong Yoon
@@ -31,60 +31,12 @@
*
*/
-enum {
- PEAK_ZAEF,
- PEAK_HDF5,
-};
-
-
-/* Information about the indexing process which is common to all patterns */
-struct index_args
-{
- UnitCell *cell;
- int config_cmfilter;
- int config_noisefilter;
- int config_verbose;
- int stream_flags; /* What goes into the output? */
- int config_satcorr;
- int config_closer;
- int config_insane;
- int config_bgsub;
- float threshold;
- float min_gradient;
- float min_snr;
- double min_int_snr;
- struct detector *det;
- IndexingMethod *indm;
- IndexingPrivate **ipriv;
- int peaks; /* Peak detection method */
- int cellr;
- float tols[4];
- struct beam_params *beam;
- char *element;
- char *hdf5_peak_path;
- double ir_inn;
- double ir_mid;
- double ir_out;
- struct copy_hdf5_field *copyme;
- int integrate_saturated;
- int use_saturated;
- int no_revalidate;
- int integrate_found;
-};
-
-
-
-
-/* Information about the indexing process for one pattern */
-struct pattern_args
-{
- /* "Input" */
- char *filename;
-
- /* "Output" */
- int indexable;
-};
+#include "index.h"
+#include "stream.h"
+#include "cell.h"
+#include "process_image.h"
extern void create_sandbox(struct index_args *iargs, int n_proc, char *prefix,
int config_basename, FILE *fh,
- char *use_this_one_instead, FILE *ofh);
+ char *use_this_one_instead, int streamfd,
+ int argc, char *argv[]);
diff --git a/src/indexamajig.c b/src/indexamajig.c
index accfa788..9b868b9f 100644
--- a/src/indexamajig.c
+++ b/src/indexamajig.c
@@ -3,13 +3,13 @@
*
* Index patterns, output hkl+intensity etc.
*
- * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY,
- * a research centre of the Helmholtz Association.
+ * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY,
+ * a research centre of the Helmholtz Association.
* Copyright © 2012 Richard Kirian
* Copyright © 2012 Lorenzo Galli
*
* Authors:
- * 2010-2012 Thomas White <taw@physics.org>
+ * 2010-2013 Thomas White <taw@physics.org>
* 2011 Richard Kirian
* 2012 Lorenzo Galli
* 2012 Chunhong Yoon
@@ -44,6 +44,9 @@
#include <getopt.h>
#include <hdf5.h>
#include <gsl/gsl_errno.h>
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <fcntl.h>
#ifdef HAVE_CLOCK_GETTIME
#include <time.h>
@@ -81,11 +84,8 @@ static void show_help(const char *s)
" Default: indexamajig.stream\n"
"\n"
" --indexing=<methods> Use 'methods' for indexing. Provide one or more\n"
-" methods separated by commas. Choose from:\n"
-" none : no indexing (default)\n"
-" dirax : invoke DirAx\n"
-" mosflm : invoke MOSFLM (DPS)\n"
-" reax : DPS algorithm with known unit cell\n"
+" methods separated by commas.\n"
+" See 'man indexamajig' for details.\n"
" -g. --geometry=<file> Get detector geometry from file.\n"
" -b, --beam=<file> Get beam parameters from file (provides nominal\n"
" wavelength value if no per-shot value is found in\n"
@@ -101,28 +101,8 @@ static void show_help(const char *s)
" --hdf5-peaks=<p> Find peaks table in HDF5 file here.\n"
" Default: /processing/hitfinder/peakinfo\n"
"\n\n"
-"You can control what information is included in the output stream using\n"
-"' --record=<flag1>,<flag2>,<flag3>' and so on. Possible flags are:\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\n"
-" search.\n"
-"\n"
-" peaksifindexed As 'peaks', but only if the pattern could be indexed.\n"
-"\n"
-" peaksifnotindexed As 'peaks', but only if the pattern could NOT be indexed.\n"
-"\n\n"
-"The default is '--record=integrated'.\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"
-" compare : match by at most changing the order of\n"
-" the indices.\n"
-" compare_ab : compare 'a' and 'b' lengths only.\n"
-" --tolerance=<tol> Set the tolerances for cell reduction.\n"
+" --tolerance=<tol> Set the tolerances for cell comparison.\n"
" Default: 5,5,5,1.5.\n"
" --filter-cm Perform common-mode noise subtraction on images\n"
" before proceeding. Intensities will be extracted\n"
@@ -134,7 +114,7 @@ static void show_help(const char *s)
" --no-sat-corr Don't correct values of saturated peaks using a\n"
" table included in the HDF5 file.\n"
" --threshold=<n> Only accept peaks above <n> ADU. Default: 800.\n"
-" --min-gradient=<n> Minimum gradient for Zaefferer peak search.\n"
+" --min-gradient=<n> Minimum squared gradient for Zaefferer peak search.\n"
" Default: 100,000.\n"
" --min-snr=<n> Minimum signal-to-noise ratio for peaks.\n"
" Default: 5.\n"
@@ -144,21 +124,20 @@ static void show_help(const char *s)
"-e, --image=<element> Use this image from the HDF5 file.\n"
" Example: /data/data0.\n"
" Default: The first one found.\n"
+" --res-cutoff Estimate the resolution limit of each pattern, and\n"
+" don't integrate reflections further out.\n"
"\n"
"\nFor time-resolved stuff, you might want to use:\n\n"
" --copy-hdf5-field <f> Copy the value of field <f> into the stream. You\n"
" can use this option as many times as you need.\n"
"\n"
-"\nOptions for greater performance or verbosity:\n\n"
-" --verbose Be verbose about indexing.\n"
+"\nOptions for greater performance:\n\n"
" -j <n> Run <n> analyses in parallel. Default 1.\n"
"\n"
"\nOptions you probably won't need:\n\n"
" --no-check-prefix Don't attempt to correct the --prefix.\n"
" --closer-peak Don't integrate from the location of a nearby peak\n"
" instead of the predicted spot. Don't use.\n"
-" --insane Don't check that the reduced cell accounts for at\n"
-" least 10%% of the located peaks.\n"
" --no-bg-sub Don't subtract local background estimates from\n"
" integrated intensities.\n"
" --use-saturated During the initial peak search, don't reject\n"
@@ -170,124 +149,110 @@ static void show_help(const char *s)
" --integrate-found Skip the spot prediction step, and just integrate\n"
" the intensities of the spots found by the initial\n"
" peak search.\n"
+" --no-peaks-in-stream Do not record peak search results in the stream.\n"
+" --no-refls-in-stream Do not record integrated reflections in the stream.\n"
);
}
-static int parse_cell_reduction(const char *scellr, int *err,
- int *reduction_needs_cell)
-{
- *err = 0;
- if ( strcmp(scellr, "none") == 0 ) {
- *reduction_needs_cell = 0;
- return CELLR_NONE;
- } else if ( strcmp(scellr, "reduce") == 0) {
- *reduction_needs_cell = 1;
- return CELLR_REDUCE;
- } else if ( strcmp(scellr, "compare") == 0) {
- *reduction_needs_cell = 1;
- return CELLR_COMPARE;
- } else if ( strcmp(scellr, "compare_ab") == 0) {
- *reduction_needs_cell = 1;
- return CELLR_COMPARE_AB;
- } else {
- *err = 1;
- *reduction_needs_cell = 0;
- return CELLR_NONE;
- }
-}
-
-
int main(int argc, char *argv[])
{
int c;
char *filename = NULL;
char *outfile = NULL;
FILE *fh;
- FILE *ofh;
+ int ofd;
char *rval = NULL;
- int config_noindex = 0;
- int config_cmfilter = 0;
- int config_noisefilter = 0;
- int config_verbose = 0;
- int config_satcorr = 1;
int config_checkprefix = 1;
- int config_closer = 0;
- int config_insane = 0;
- int config_bgsub = 1;
int config_basename = 0;
- float threshold = 800.0;
- float min_gradient = 100000.0;
- float min_snr = 5.0;
- double min_int_snr = -INFINITY;
- struct detector *det;
- char *geometry = NULL;
IndexingMethod *indm;
IndexingPrivate **ipriv;
- int indexer_needs_cell;
- int reduction_needs_cell;
char *indm_str = NULL;
- UnitCell *cell;
char *pdb = NULL;
char *prefix = NULL;
char *speaks = NULL;
- char *scellr = NULL;
char *toler = NULL;
- float tols[4] = {5.0, 5.0, 5.0, 1.5}; /* a,b,c,angles (%,%,%,deg) */
- int cellr;
- int peaks;
int n_proc = 1;
char *prepare_line;
char prepare_filename[1024];
char *use_this_one_instead;
struct index_args iargs;
- struct beam_params *beam = NULL;
- char *element = NULL;
- double nominal_photon_energy;
- int stream_flags = STREAM_INTEGRATED;
- char *hdf5_peak_path = NULL;
- struct copy_hdf5_field *copyme;
char *intrad = NULL;
- float ir_inn = 4.0;
- float ir_mid = 5.0;
- float ir_out = 7.0;
- int integrate_saturated = 0;
- int use_saturated = 0;
- int no_revalidate = 0;
- int integrate_found = 0;
-
- copyme = new_copy_hdf5_field_list();
- if ( copyme == NULL ) {
+
+ /* Defaults */
+ iargs.cell = NULL;
+ iargs.cmfilter = 0;
+ iargs.noisefilter = 0;
+ iargs.satcorr = 1;
+ iargs.closer = 0;
+ iargs.bgsub = 1;
+ iargs.tols[0] = 5.0;
+ iargs.tols[1] = 5.0;
+ iargs.tols[2] = 5.0;
+ iargs.tols[3] = 1.5;
+ iargs.threshold = 800.0;
+ iargs.min_gradient = 100000.0;
+ iargs.min_snr = 5.0;
+ iargs.min_int_snr = -INFINITY;
+ iargs.det = NULL;
+ iargs.peaks = PEAK_ZAEF;
+ iargs.beam = NULL;
+ iargs.element = NULL;
+ iargs.hdf5_peak_path = strdup("/processing/hitfinder/peakinfo");
+ iargs.copyme = NULL;
+ iargs.ir_inn = 4.0;
+ iargs.ir_mid = 5.0;
+ iargs.ir_out = 7.0;
+ iargs.use_saturated = 0;
+ iargs.integrate_saturated = 0;
+ iargs.no_revalidate = 0;
+ iargs.integrate_found = 0;
+ iargs.stream_peaks = 1;
+ iargs.stream_refls = 1;
+ iargs.res_cutoff = 0;
+ iargs.copyme = new_copy_hdf5_field_list();
+ if ( iargs.copyme == NULL ) {
ERROR("Couldn't allocate HDF5 field list.\n");
return 1;
}
+ iargs.indm = NULL; /* No default */
+ iargs.ipriv = NULL; /* No default */
/* Long options */
const struct option longopts[] = {
+
+ /* Options with long and short versions */
{"help", 0, NULL, 'h'},
{"input", 1, NULL, 'i'},
{"output", 1, NULL, 'o'},
- {"no-index", 0, &config_noindex, 1},
{"indexing", 1, NULL, 'z'},
{"geometry", 1, NULL, 'g'},
{"beam", 1, NULL, 'b'},
- {"filter-cm", 0, &config_cmfilter, 1},
- {"filter-noise", 0, &config_noisefilter, 1},
- {"verbose", 0, &config_verbose, 1},
{"pdb", 1, NULL, 'p'},
{"prefix", 1, NULL, 'x'},
- {"no-sat-corr", 0, &config_satcorr, 0},
- {"sat-corr", 0, &config_satcorr, 1}, /* Compat */
{"threshold", 1, NULL, 't'},
- {"no-check-prefix", 0, &config_checkprefix, 0},
- {"no-closer-peak", 0, &config_closer, 0},
- {"closer-peak", 0, &config_closer, 1},
- {"insane", 0, &config_insane, 1},
{"image", 1, NULL, 'e'},
- {"basename", 0, &config_basename, 1},
- {"bg-sub", 0, &config_bgsub, 1}, /* Compat */
- {"no-bg-sub", 0, &config_bgsub, 0},
+ /* Long-only options with no arguments */
+ {"filter-cm", 0, &iargs.cmfilter, 1},
+ {"filter-noise", 0, &iargs.noisefilter, 1},
+ {"no-sat-corr", 0, &iargs.satcorr, 0},
+ {"sat-corr", 0, &iargs.satcorr, 1},
+ {"no-check-prefix", 0, &config_checkprefix, 0},
+ {"no-closer-peak", 0, &iargs.closer, 0},
+ {"closer-peak", 0, &iargs.closer, 1},
+ {"basename", 0, &config_basename, 1},
+ {"bg-sub", 0, &iargs.bgsub, 1},
+ {"no-bg-sub", 0, &iargs.bgsub, 0},
+ {"no-peaks-in-stream", 0, &iargs.stream_peaks, 0},
+ {"no-refls-in-stream", 0, &iargs.stream_refls, 0},
+ {"res-cutoff", 0, &iargs.res_cutoff, 1},
+ {"integrate-saturated",0, &iargs.integrate_saturated,1},
+ {"use-saturated", 0, &iargs.use_saturated, 1},
+ {"no-revalidate", 0, &iargs.no_revalidate, 1},
+ {"integrate-found", 0, &iargs.integrate_found, 1},
+
+ /* Long-only options with arguments */
{"peaks", 1, NULL, 2},
{"cell-reduction", 1, NULL, 3},
{"min-gradient", 1, NULL, 4},
@@ -302,10 +267,6 @@ int main(int argc, char *argv[])
{"tolerance", 1, NULL, 13},
{"int-radius", 1, NULL, 14},
- {"integrate-saturated",0, &integrate_saturated,1},
- {"use-saturated", 0, &use_saturated, 1},
- {"no-revalidate", 0, &no_revalidate, 1},
- {"integrate-found", 0, &integrate_found, 1},
{0, 0, NULL, 0}
};
@@ -344,16 +305,21 @@ int main(int argc, char *argv[])
break;
case 'g' :
- geometry = strdup(optarg);
+ iargs.det = get_detector_geometry(optarg);
+ if ( iargs.det == NULL ) {
+ ERROR("Failed to read detector geometry from "
+ "'%s'\n", optarg);
+ return 1;
+ }
break;
case 't' :
- threshold = strtof(optarg, NULL);
+ iargs.threshold = strtof(optarg, NULL);
break;
case 'b' :
- beam = get_beam_parameters(optarg);
- if ( beam == NULL ) {
+ iargs.beam = get_beam_parameters(optarg);
+ if ( iargs.beam == NULL ) {
ERROR("Failed to load beam parameters"
" from '%s'\n", optarg);
return 1;
@@ -361,7 +327,7 @@ int main(int argc, char *argv[])
break;
case 'e' :
- element = strdup(optarg);
+ iargs.element = strdup(optarg);
break;
case 2 :
@@ -369,17 +335,24 @@ int main(int argc, char *argv[])
break;
case 3 :
- scellr = strdup(optarg);
- break;
+ ERROR("The option '--cell-reduction' is no longer "
+ "used.\n"
+ "The complete indexing behaviour is now "
+ "controlled using '--indexing'.\n"
+ "See 'man indexamajig' for details of the "
+ "available methods.\n");
+ return 1;
case 4 :
- min_gradient = strtof(optarg, NULL);
+ iargs.min_gradient = strtof(optarg, NULL);
break;
case 5 :
- stream_flags = parse_stream_flags(optarg);
- if ( stream_flags < 0 ) return 1;
- break;
+ ERROR("The option '--record' is no longer used.\n"
+ "Use '--no-peaks-in-stream' and"
+ "'--no-refls-in-stream' if you need to control"
+ "the contents of the stream.\n");
+ return 1;
case 6 :
case 7 :
@@ -389,19 +362,20 @@ int main(int argc, char *argv[])
break;
case 9 :
- hdf5_peak_path = strdup(optarg);
+ free(iargs.hdf5_peak_path);
+ iargs.hdf5_peak_path = strdup(optarg);
break;
case 10 :
- add_copy_hdf5_field(copyme, optarg);
+ add_copy_hdf5_field(iargs.copyme, optarg);
break;
case 11 :
- min_snr = strtof(optarg, NULL);
+ iargs.min_snr = strtof(optarg, NULL);
break;
case 12 :
- min_int_snr = strtof(optarg, NULL);
+ iargs.min_int_snr = strtof(optarg, NULL);
break;
case 13 :
@@ -440,30 +414,15 @@ int main(int argc, char *argv[])
}
free(filename);
- if ( outfile == NULL ) {
- ofh = stdout;
- } else {
- ofh = fopen(outfile, "w");
- if ( ofh == NULL ) {
- ERROR("Failed to open output file '%s'\n", outfile);
- return 1;
- }
- free(outfile);
- }
-
- if ( hdf5_peak_path == NULL ) {
- hdf5_peak_path = strdup("/processing/hitfinder/peakinfo");
- }
-
if ( speaks == NULL ) {
speaks = strdup("zaef");
STATUS("You didn't specify a peak detection method.\n");
STATUS("I'm using 'zaef' for you.\n");
}
if ( strcmp(speaks, "zaef") == 0 ) {
- peaks = PEAK_ZAEF;
+ iargs.peaks = PEAK_ZAEF;
} else if ( strcmp(speaks, "hdf5") == 0 ) {
- peaks = PEAK_HDF5;
+ iargs.peaks = PEAK_HDF5;
} else {
ERROR("Unrecognised peak detection method '%s'\n", speaks);
return 1;
@@ -483,57 +442,29 @@ int main(int argc, char *argv[])
return 1;
}
- if ( (indm_str == NULL) ||
- ((indm_str != NULL) && (strcmp(indm_str, "none") == 0)) ) {
- STATUS("Not indexing anything.\n");
- indexer_needs_cell = 0;
- reduction_needs_cell = 0;
+ if ( indm_str == NULL ) {
+
+ STATUS("You didn't specify an indexing method, so I won't try "
+ " to index anything.\n"
+ "If that isn't what you wanted, re-run with"
+ " --indexing=<methods>.\n");
indm = NULL;
- cellr = CELLR_NONE;
+
} else {
- if ( indm_str == NULL ) {
- STATUS("You didn't specify an indexing method, so I "
- " won't try to index anything.\n"
- "If that isn't what you wanted, re-run with"
- " --indexing=<method>.\n");
- indm = NULL;
- indexer_needs_cell = 0;
- } else {
- indm = build_indexer_list(indm_str,
- &indexer_needs_cell);
- if ( indm == NULL ) {
- ERROR("Invalid indexer list '%s'\n", indm_str);
- return 1;
- }
- free(indm_str);
- }
- reduction_needs_cell = 0;
- if ( scellr == NULL ) {
- STATUS("You didn't specify a cell reduction method, so"
- " I'm going to use 'reduce'.\n");
- cellr = CELLR_REDUCE;
- reduction_needs_cell = 1;
- } else {
- int err;
- cellr = parse_cell_reduction(scellr, &err,
- &reduction_needs_cell);
- if ( err ) {
- ERROR("Unrecognised cell reduction '%s'\n",
- scellr);
- return 1;
- }
- free(scellr);
+ indm = build_indexer_list(indm_str);
+ if ( indm == NULL ) {
+ ERROR("Invalid indexer list '%s'\n", indm_str);
+ return 1;
}
+ free(indm_str);
}
- /* No indexing -> no reduction */
- if ( indm == NULL ) reduction_needs_cell = 0;
-
if ( toler != NULL ) {
int ttt;
ttt = sscanf(toler, "%f,%f,%f,%f",
- &tols[0], &tols[1], &tols[2], &tols[3] );
+ &iargs.tols[0], &iargs.tols[1],
+ &iargs.tols[2], &iargs.tols[3]);
if ( ttt != 4 ) {
ERROR("Invalid parameters for '--tolerance'\n");
return 1;
@@ -543,7 +474,8 @@ int main(int argc, char *argv[])
if ( intrad != NULL ) {
int r;
- r = sscanf(intrad, "%f,%f,%f", &ir_inn, &ir_mid, &ir_out);
+ r = sscanf(intrad, "%f,%f,%f",
+ &iargs.ir_inn, &iargs.ir_mid, &iargs.ir_out);
if ( r != 3 ) {
ERROR("Invalid parameters for '--int-radius'\n");
return 1;
@@ -555,43 +487,38 @@ int main(int argc, char *argv[])
" probably not appropriate for your patterns.\n");
}
- if ( geometry == NULL ) {
- ERROR("You need to specify a geometry file with --geometry\n");
+ if ( iargs.det == NULL ) {
+ ERROR("You need to provide a geometry file (please read the"
+ " manual for more details).\n");
return 1;
}
- det = get_detector_geometry(geometry);
- if ( det == NULL ) {
- ERROR("Failed to read detector geometry from '%s'\n", geometry);
+ if ( iargs.beam == NULL ) {
+ ERROR("You need to provide a beam parameters file (please read"
+ " the manual for more details).\n");
return 1;
}
- free(geometry);
if ( pdb != NULL ) {
- cell = load_cell_from_pdb(pdb);
- if ( cell == NULL ) {
+ iargs.cell = load_cell_from_pdb(pdb);
+ if ( iargs.cell == NULL ) {
ERROR("Couldn't read unit cell (from %s)\n", pdb);
return 1;
}
free(pdb);
- cell_print(cell);
+ STATUS("This is what I understood your unit cell to be:\n");
+ cell_print(iargs.cell);
} else {
STATUS("No unit cell given.\n");
- cell = NULL;
+ iargs.cell = NULL;
}
- write_stream_header(ofh, argc, argv);
-
- if ( beam != NULL ) {
- nominal_photon_energy = beam->photon_energy;
- } else {
- STATUS("No beam parameters file was given, so I'm taking the"
- " nominal photon energy to be 2 keV.\n");
- ERROR("I'm also going to assume 1 ADU per photon, which is");
- ERROR(" almost certainly wrong. Peak sigmas will be"
- " incorrect.\n");
- nominal_photon_energy = 2000.0;
+ ofd = open(outfile, O_CREAT | O_TRUNC | O_WRONLY, S_IRUSR | S_IWUSR);
+ if ( ofd == -1 ) {
+ ERROR("Failed to open stream '%s'\n", outfile);
+ return 1;
}
+ free(outfile);
/* Get first filename and use it to set up the indexing */
prepare_line = malloc(1024);
@@ -613,8 +540,8 @@ int main(int argc, char *argv[])
/* Prepare the indexer */
if ( indm != NULL ) {
- ipriv = prepare_indexing(indm, cell, prepare_filename, det,
- nominal_photon_energy);
+ ipriv = prepare_indexing(indm, iargs.cell, prepare_filename,
+ iargs.det, iargs.beam, iargs.tols);
if ( ipriv == NULL ) {
ERROR("Failed to prepare indexing.\n");
return 1;
@@ -625,43 +552,11 @@ int main(int argc, char *argv[])
gsl_set_error_handler_off();
- /* Static worker args */
- iargs.cell = cell;
- iargs.config_cmfilter = config_cmfilter;
- iargs.config_noisefilter = config_noisefilter;
- iargs.config_verbose = config_verbose;
- iargs.config_satcorr = config_satcorr;
- iargs.config_closer = config_closer;
- iargs.config_insane = config_insane;
- iargs.config_bgsub = config_bgsub;
- iargs.cellr = cellr;
- iargs.tols[0] = tols[0];
- iargs.tols[1] = tols[1];
- iargs.tols[2] = tols[2];
- iargs.tols[3] = tols[3];
- iargs.threshold = threshold;
- iargs.min_gradient = min_gradient;
- iargs.min_snr = min_snr;
- iargs.min_int_snr = min_int_snr;
- iargs.det = det;
iargs.indm = indm;
iargs.ipriv = ipriv;
- iargs.peaks = peaks;
- iargs.beam = beam;
- iargs.element = element;
- iargs.stream_flags = stream_flags;
- iargs.hdf5_peak_path = hdf5_peak_path;
- iargs.copyme = copyme;
- iargs.ir_inn = ir_inn;
- iargs.ir_mid = ir_mid;
- iargs.ir_out = ir_out;
- iargs.use_saturated = use_saturated;
- iargs.integrate_saturated = integrate_saturated;
- iargs.no_revalidate = no_revalidate;
- iargs.integrate_found = integrate_found;
create_sandbox(&iargs, n_proc, prefix, config_basename, fh,
- use_this_one_instead, ofh);
+ use_this_one_instead, ofd, argc, argv);
free(prefix);
diff --git a/src/partial_sim.c b/src/partial_sim.c
index facf14ef..a4a2f875 100644
--- a/src/partial_sim.c
+++ b/src/partial_sim.c
@@ -54,11 +54,12 @@
#define NBINS 50
-static void mess_up_cell(UnitCell *cell, double cnoise)
+static void mess_up_cell(Crystal *cr, double cnoise)
{
double ax, ay, az;
double bx, by, bz;
double cx, cy, cz;
+ UnitCell *cell = crystal_get_cell(cr);
//STATUS("Real:\n");
//cell_print(cell);
@@ -82,17 +83,18 @@ static void mess_up_cell(UnitCell *cell, double cnoise)
/* For each reflection in "partial", fill in what the intensity would be
* according to "full" */
-static void calculate_partials(RefList *partial, double osf,
+static void calculate_partials(Crystal *cr,
RefList *full, const SymOpList *sym,
int random_intensities,
pthread_mutex_t *full_lock,
unsigned long int *n_ref, double *p_hist,
- double *p_max, double max_q, UnitCell *cell)
+ double *p_max, double max_q)
{
Reflection *refl;
RefListIterator *iter;
+ double res;
- for ( refl = first_refl(partial, &iter);
+ for ( refl = first_refl(crystal_get_reflections(cr), &iter);
refl != NULL;
refl = next_refl(refl, iter) )
{
@@ -137,16 +139,17 @@ static void calculate_partials(RefList *partial, double osf,
}
}
- Ip = osf * p * If;
+ Ip = crystal_get_osf(cr) * p * If;
- bin = NBINS*2.0*resolution(cell, h, k, l) / max_q;
+ res = resolution(crystal_get_cell(cr), h, k, l);
+ bin = NBINS*2.0*res/max_q;
if ( (bin < NBINS) && (bin>=0) ) {
p_hist[bin] += p;
n_ref[bin]++;
if ( p > p_max[bin] ) p_max[bin] = p;
} else {
STATUS("Reflection out of histogram range: %e %i %f\n",
- resolution(cell, h, k, l), bin, p);
+ res, bin, p);
}
Ip = gaussian_noise(Ip, 100.0);
@@ -206,13 +209,14 @@ struct queue_args
unsigned long int n_ref[NBINS];
double p_max[NBINS];
- FILE *stream;
+ Stream *stream;
};
struct worker_args
{
struct queue_args *qargs;
+ Crystal *crystal;
struct image image;
/* Histogram for this image */
@@ -243,21 +247,31 @@ static void *create_job(void *vqargs)
static void run_job(void *vwargs, int cookie)
{
- double osf;
struct quaternion orientation;
struct worker_args *wargs = vwargs;
struct queue_args *qargs = wargs->qargs;
int i;
+ Crystal *cr;
+ RefList *reflections;
- osf = gaussian_noise(1.0, 0.3);
+ cr = crystal_new();
+ if ( cr == NULL ) {
+ ERROR("Failed to create crystal.\n");
+ return;
+ }
+ wargs->crystal = cr;
+ crystal_set_image(cr, &wargs->image);
+
+ crystal_set_osf(cr, gaussian_noise(1.0, 0.3));
+ crystal_set_profile_radius(cr, wargs->image.beam->profile_radius);
/* Set up a random orientation */
orientation = random_quaternion();
- wargs->image.indexed_cell = cell_rotate(qargs->cell, orientation);
+ crystal_set_cell(cr, cell_rotate(qargs->cell, orientation));
snprintf(wargs->image.filename, 255, "dummy.h5");
- wargs->image.reflections = find_intersections(&wargs->image,
- wargs->image.indexed_cell);
+ reflections = find_intersections(&wargs->image, cr);
+ crystal_set_reflections(cr, reflections);
for ( i=0; i<NBINS; i++ ) {
wargs->n_ref[i] = 0;
@@ -265,14 +279,16 @@ static void run_job(void *vwargs, int cookie)
wargs->p_max[i] = 0.0;
}
- calculate_partials(wargs->image.reflections, osf, qargs->full,
+ calculate_partials(cr, qargs->full,
qargs->sym, qargs->random_intensities,
&qargs->full_lock,
wargs->n_ref, wargs->p_hist, wargs->p_max,
- qargs->max_q, wargs->image.indexed_cell);
+ qargs->max_q);
/* Give a slightly incorrect cell in the stream */
- mess_up_cell(wargs->image.indexed_cell, qargs->cnoise);
+ mess_up_cell(cr, qargs->cnoise);
+
+ image_add_crystal(&wargs->image, cr);
}
@@ -282,7 +298,7 @@ static void finalise_job(void *vqargs, void *vwargs)
struct queue_args *qargs = vqargs;
int i;
- write_chunk(qargs->stream, &wargs->image, NULL, STREAM_INTEGRATED);
+ write_chunk(qargs->stream, &wargs->image, NULL, 0, 1);
for ( i=0; i<NBINS; i++ ) {
qargs->n_ref[i] += wargs->n_ref[i];
@@ -295,8 +311,7 @@ static void finalise_job(void *vqargs, void *vwargs)
qargs->n_done++;
progress_bar(qargs->n_done, qargs->n_to_do, "Simulating");
- reflist_free(wargs->image.reflections);
- cell_free(wargs->image.indexed_cell);
+ crystal_free(wargs->crystal);
free(wargs);
}
@@ -315,7 +330,7 @@ int main(int argc, char *argv[])
char *sym_str = NULL;
SymOpList *sym;
UnitCell *cell = NULL;
- FILE *ofh;
+ Stream *stream;
int n = 2;
int random_intensities = 0;
char *save_file = NULL;
@@ -496,13 +511,12 @@ int main(int argc, char *argv[])
ERROR("You must give a filename for the output.\n");
return 1;
}
- ofh = fopen(output_file, "w");
- if ( ofh == NULL ) {
+ stream = open_stream_for_write(output_file);
+ if ( stream == NULL ) {
ERROR("Couldn't open output file '%s'\n", output_file);
return 1;
}
free(output_file);
- write_stream_header(ofh, argc, argv);
image.det = det;
image.width = det->max_fs;
@@ -511,9 +525,12 @@ int main(int argc, char *argv[])
image.lambda = ph_en_to_lambda(eV_to_J(beam->photon_energy));
image.div = beam->divergence;
image.bw = beam->bandwidth;
- image.profile_radius = beam->profile_radius;
+ image.beam = beam;
image.filename = malloc(256);
image.copyme = NULL;
+ image.crystals = NULL;
+ image.n_crystals = 0;
+ image.indexed_by = INDEXING_NONE;
if ( random_intensities ) {
full = reflist_new();
@@ -528,7 +545,7 @@ int main(int argc, char *argv[])
qargs.random_intensities = random_intensities;
qargs.cell = cell;
qargs.template_image = &image;
- qargs.stream = ofh;
+ qargs.stream = stream;
qargs.cnoise = cnoise;
qargs.max_q = largest_q(&image);
@@ -574,7 +591,7 @@ int main(int argc, char *argv[])
}
- fclose(ofh);
+ close_stream(stream);
cell_free(cell);
free_detector_geometry(det);
free(beam);
diff --git a/src/partialator.c b/src/partialator.c
index c2f6c299..a4af3c18 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -3,11 +3,11 @@
*
* Scaling and post refinement for coherent nanocrystallography
*
- * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY,
- * a research centre of the Helmholtz Association.
+ * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY,
+ * a research centre of the Helmholtz Association.
*
* Authors:
- * 2010-2012 Thomas White <taw@physics.org>
+ * 2010-2013 Thomas White <taw@physics.org>
*
* This file is part of CrystFEL.
*
@@ -87,16 +87,16 @@ static void show_help(const char *s)
struct refine_args
{
RefList *full;
- struct image *image;
+ Crystal *crystal;
};
struct queue_args
{
- int n;
+ int n_started;
int n_done;
- int n_total_patterns;
- struct image *images;
+ Crystal **crystals;
+ int n_crystals;
struct refine_args task_defaults;
};
@@ -104,10 +104,9 @@ struct queue_args
static void refine_image(void *task, int id)
{
struct refine_args *pargs = task;
- struct image *image = pargs->image;
- image->id = id;
+ Crystal *cr = pargs->crystal;
- pr_refine(image, pargs->full);
+ pr_refine(cr, pargs->full);
}
@@ -119,9 +118,9 @@ static void *get_image(void *vqargs)
task = malloc(sizeof(struct refine_args));
memcpy(task, &qargs->task_defaults, sizeof(struct refine_args));
- task->image = &qargs->images[qargs->n];
+ task->crystal = qargs->crystals[qargs->n_started];
- qargs->n++;
+ qargs->n_started++;
return task;
}
@@ -133,12 +132,12 @@ static void done_image(void *vqargs, void *task)
qargs->n_done++;
- progress_bar(qargs->n_done, qargs->n_total_patterns, "Refining");
+ progress_bar(qargs->n_done, qargs->n_crystals, "Refining");
free(task);
}
-static void refine_all(struct image *images, int n_total_patterns,
+static void refine_all(Crystal **crystals, int n_crystals,
struct detector *det,
RefList *full, int nthreads)
{
@@ -146,19 +145,19 @@ static void refine_all(struct image *images, int n_total_patterns,
struct queue_args qargs;
task_defaults.full = full;
- task_defaults.image = NULL;
+ task_defaults.crystal = NULL;
qargs.task_defaults = task_defaults;
- qargs.n = 0;
+ qargs.n_started = 0;
qargs.n_done = 0;
- qargs.n_total_patterns = n_total_patterns;
- qargs.images = images;
+ qargs.n_crystals = n_crystals;
+ qargs.crystals = crystals;
/* Don't have threads which are doing nothing */
- if ( n_total_patterns < nthreads ) nthreads = n_total_patterns;
+ if ( n_crystals < nthreads ) nthreads = n_crystals;
run_threads(nthreads, refine_image, get_image, done_image,
- &qargs, n_total_patterns, 0, 0, 0);
+ &qargs, n_crystals, 0, 0, 0);
}
@@ -201,13 +200,14 @@ static int select_scalable_reflections(RefList *list, RefList *reference)
}
-static void select_reflections_for_refinement(struct image *images, int n,
+static void select_reflections_for_refinement(Crystal **crystals, int n,
RefList *full, int have_reference)
{
int i;
for ( i=0; i<n; i++ ) {
+ RefList *reflist;
Reflection *refl;
RefListIterator *iter;
int n_acc = 0;
@@ -216,9 +216,10 @@ static void select_reflections_for_refinement(struct image *images, int n,
int n_fewmatch = 0;
int n_ref = 0;
- if ( images[i].pr_dud ) continue;
+ if ( crystal_get_user_flag(crystals[i]) ) continue;
- for ( refl = first_refl(images[i].reflections, &iter);
+ reflist = crystal_get_reflections(crystals[i]);
+ for ( refl = first_refl(reflist, &iter);
refl != NULL;
refl = next_refl(refl, iter) )
{
@@ -272,6 +273,20 @@ static void select_reflections_for_refinement(struct image *images, int n,
}
+static void display_progress(int n_images, int n_crystals)
+{
+ if ( !isatty(STDERR_FILENO) ) return;
+ if ( tcgetpgrp(STDERR_FILENO) != getpgrp() ) return;
+
+ pthread_mutex_lock(&stderr_lock);
+ fprintf(stderr, "\r%i images loaded, %i crystals.",
+ n_images, n_crystals);
+ pthread_mutex_unlock(&stderr_lock);
+
+ fflush(stdout);
+}
+
+
int main(int argc, char *argv[])
{
int c;
@@ -280,16 +295,15 @@ int main(int argc, char *argv[])
char *geomfile = NULL;
char *sym_str = NULL;
SymOpList *sym;
- FILE *fh;
int nthreads = 1;
struct detector *det;
int i;
- int n_total_patterns;
struct image *images;
int n_iter = 10;
struct beam_params *beam = NULL;
RefList *full;
- int n_usable_patterns = 0;
+ int n_images = 0;
+ int n_crystals = 0;
int nobs;
char *reference_file = NULL;
RefList *reference = NULL;
@@ -298,6 +312,8 @@ int main(int argc, char *argv[])
char cmdline[1024];
SRContext *sr;
int noscale = 0;
+ Stream *st;
+ Crystal **crystals;
/* Long options */
const struct option longopts[] = {
@@ -386,19 +402,15 @@ int main(int argc, char *argv[])
return 1;
}
- /* Sanitise input filename and open */
if ( infile == NULL ) {
infile = strdup("-");
}
- if ( strcmp(infile, "-") == 0 ) {
- fh = stdin;
- } else {
- fh = fopen(infile, "r");
- }
- if ( fh == NULL ) {
- ERROR("Failed to open input file '%s'\n", infile);
+ st = open_stream_for_read(infile);
+ if ( st == NULL ) {
+ ERROR("Failed to open input stream '%s'\n", infile);
return 1;
}
+ /* Don't free "infile", because it's needed for the scaling report */
/* Sanitise output filename */
if ( outfile == NULL ) {
@@ -438,86 +450,116 @@ int main(int argc, char *argv[])
}
- n_total_patterns = count_patterns(fh);
- if ( n_total_patterns == 0 ) {
- ERROR("No patterns to process.\n");
- return 1;
- }
- STATUS("There are %i patterns to process\n", n_total_patterns);
-
gsl_set_error_handler_off();
- images = malloc(n_total_patterns * sizeof(struct image));
- if ( images == NULL ) {
- ERROR("Couldn't allocate memory for images.\n");
- return 1;
- }
-
/* Fill in what we know about the images so far */
- rewind(fh);
- nobs = 0;
- for ( i=0; i<n_total_patterns; i++ ) {
+ n_images = 0;
+ n_crystals = 0;
+ images = NULL;
+ crystals = NULL;
+ do {
RefList *as;
+ int i;
+ struct image *images_new;
struct image *cur;
- cur = &images[n_usable_patterns];
+ images_new = realloc(images, (n_images+1)*sizeof(struct image));
+ if ( images_new == NULL ) {
+ ERROR("Failed to allocate memory for image list.\n");
+ return 1;
+ }
+ images = images_new;
+ cur = &images[n_images];
cur->det = det;
-
- if ( read_chunk(fh, cur) != 0 ) {
- /* Should not happen, because we counted the patterns
- * earlier. */
- ERROR("Failed to read chunk from the input stream.\n");
- return 1;
+ if ( read_chunk(st, cur) != 0 ) {
+ break;
}
/* Won't be needing this, if it exists */
image_feature_list_free(cur->features);
cur->features = NULL;
-
- /* "n_usable_patterns" will not be incremented in this case */
- if ( cur->indexed_cell == NULL ) continue;
-
- /* Fill in initial estimates of stuff */
cur->div = beam->divergence;
cur->bw = beam->bandwidth;
cur->width = det->max_fs;
cur->height = det->max_ss;
- cur->osf = 1.0;
- cur->profile_radius = beam->profile_radius;
- cur->pr_dud = 0;
-
- /* Muppet proofing */
cur->data = NULL;
cur->flags = NULL;
cur->beam = NULL;
- /* This is the raw list of reflections */
- as = asymmetric_indices(cur->reflections, sym);
- reflist_free(cur->reflections);
- cur->reflections = as;
+ n_images++;
+
+ for ( i=0; i<cur->n_crystals; i++ ) {
+
+ Crystal *cr;
+ Crystal **crystals_new;
+ RefList *cr_refl;
+
+ crystals_new = realloc(crystals,
+ (n_crystals+1)*sizeof(Crystal *));
+ if ( crystals_new == NULL ) {
+ ERROR("Failed to allocate memory for crystal "
+ "list.\n");
+ return 1;
+ }
+ crystals = crystals_new;
+ crystals[n_crystals] = cur->crystals[i];
+ cr = crystals[n_crystals];
+
+ /* Image pointer will change due to later reallocs */
+ crystal_set_image(cr, NULL);
+
+ /* Fill in initial estimates of stuff */
+ crystal_set_osf(cr, 1.0);
+ crystal_set_profile_radius(cr, beam->profile_radius);
+ crystal_set_user_flag(cr, 0);
+
+ /* This is the raw list of reflections */
+ cr_refl = crystal_get_reflections(cr);
+ as = asymmetric_indices(cr_refl, sym);
+ crystal_set_reflections(cr, as);
+ reflist_free(cr_refl);
- update_partialities(cur);
+ n_crystals++;
- nobs += select_scalable_reflections(cur->reflections,
- reference);
+ }
+
+ display_progress(n_images, n_crystals);
+
+ } while ( 1 );
+
+ close_stream(st);
+
+ /* Fill in image pointers */
+ nobs = 0;
+ for ( i=0; i<n_images; i++ ) {
+ int j;
+ for ( j=0; j<images[i].n_crystals; j++ ) {
+
+ Crystal *cryst;
+ RefList *as;
- progress_bar(i, n_total_patterns-1, "Loading pattern data");
- n_usable_patterns++;
+ cryst = images[i].crystals[j];
+ crystal_set_image(cryst, &images[i]);
+ /* Now it's safe to do the following */
+ update_partialities(cryst);
+ as = crystal_get_reflections(cryst);
+ nobs += select_scalable_reflections(as, reference);
+
+ }
}
- fclose(fh);
/* Make initial estimates */
- STATUS("Performing initial scaling.\n");
+ STATUS("\nPerforming initial scaling.\n");
if ( noscale ) STATUS("Scale factors fixed at 1.\n");
- full = scale_intensities(images, n_usable_patterns, reference,
+ full = scale_intensities(crystals, n_crystals, reference,
nthreads, noscale);
- sr = sr_titlepage(images, n_usable_patterns, "scaling-report.pdf",
+ sr = sr_titlepage(crystals, n_crystals, "scaling-report.pdf",
infile, cmdline);
- sr_iteration(sr, 0, images, n_usable_patterns, full);
+ sr_iteration(sr, 0, crystals, n_crystals, full);
/* Iterate */
for ( i=0; i<n_iter; i++ ) {
@@ -534,54 +576,55 @@ int main(int argc, char *argv[])
}
/* Refine the geometry of all patterns to get the best fit */
- select_reflections_for_refinement(images, n_usable_patterns,
+ select_reflections_for_refinement(crystals, n_crystals,
comp, have_reference);
- refine_all(images, n_usable_patterns, det, comp, nthreads);
+ refine_all(crystals, n_crystals, det, comp, nthreads);
nobs = 0;
- for ( j=0; j<n_usable_patterns; j++ ) {
+ for ( j=0; j<n_crystals; j++ ) {
- struct image *cur = &images[j];
+ Crystal *cr = crystals[j];
+ RefList *rf = crystal_get_reflections(cr);
- nobs += select_scalable_reflections(cur->reflections,
- reference);
+ nobs += select_scalable_reflections(rf, reference);
}
/* Re-estimate all the full intensities */
reflist_free(full);
- full = scale_intensities(images, n_usable_patterns,
+ full = scale_intensities(crystals, n_crystals,
reference, nthreads, noscale);
- sr_iteration(sr, i+1, images, n_usable_patterns, full);
+ sr_iteration(sr, i+1, crystals, n_crystals, full);
}
sr_finish(sr);
n_dud = 0;
- for ( i=0; i<n_usable_patterns; i++ ) {
- if ( images[i].pr_dud ) n_dud++;
+ for ( i=0; i<n_crystals; i++ ) {
+ if ( crystal_get_user_flag(crystals[i]) ) n_dud++;
}
- STATUS("%i images could not be refined on the last cycle.\n", n_dud);
+ STATUS("%i crystals could not be refined on the last cycle.\n", n_dud);
/* Output results */
write_reflist(outfile, full);
/* Clean up */
- for ( i=0; i<n_usable_patterns; i++ ) {
- reflist_free(images[i].reflections);
+ for ( i=0; i<n_crystals; i++ ) {
+ reflist_free(crystal_get_reflections(crystals[i]));
+ crystal_free(crystals[i]);
}
reflist_free(full);
free(sym);
free(outfile);
free_detector_geometry(det);
free(beam);
+ free(crystals);
if ( reference != NULL ) {
reflist_free(reference);
}
- for ( i=0; i<n_usable_patterns; i++ ) {
- cell_free(images[i].indexed_cell);
+ for ( i=0; i<n_images; i++ ) {
free(images[i].filename);
}
free(images);
diff --git a/src/post-refinement.c b/src/post-refinement.c
index 7410d931..1439b148 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -3,11 +3,11 @@
*
* Post refinement
*
- * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY,
- * a research centre of the Helmholtz Association.
+ * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY,
+ * a research centre of the Helmholtz Association.
*
* Authors:
- * 2010-2012 Thomas White <taw@physics.org>
+ * 2010-2013 Thomas White <taw@physics.org>
*
* This file is part of CrystFEL.
*
@@ -88,7 +88,7 @@ static double partiality_rgradient(double r, double profile_radius)
/* Return the gradient of parameter 'k' given the current status of 'image'. */
-double gradient(struct image *image, int k, Reflection *refl, double r)
+double gradient(Crystal *cr, int k, Reflection *refl)
{
double ds, azix, aziy;
double ttlow, tthigh, tt;
@@ -103,17 +103,19 @@ double gradient(struct image *image, int k, Reflection *refl, double r)
int clamp_low, clamp_high;
double klow, khigh;
double gr;
+ struct image *image = crystal_get_image(cr);
+ double r = crystal_get_profile_radius(cr);
get_symmetric_indices(refl, &hs, &ks, &ls);
- cell_get_reciprocal(image->indexed_cell, &asx, &asy, &asz,
- &bsx, &bsy, &bsz,
- &csx, &csy, &csz);
+ cell_get_reciprocal(crystal_get_cell(cr), &asx, &asy, &asz,
+ &bsx, &bsy, &bsz,
+ &csx, &csy, &csz);
xl = hs*asx + ks*bsx + ls*csx;
yl = hs*asy + ks*bsy + ls*csy;
zl = hs*asz + ks*bsz + ls*csz;
- ds = 2.0 * resolution(image->indexed_cell, hs, ks, ls);
+ ds = 2.0 * resolution(crystal_get_cell(cr), hs, ks, ls);
get_partial(refl, &r1, &r2, &p, &clamp_low, &clamp_high);
klow = 1.0/(image->lambda - image->lambda*image->bw/2.0);
@@ -237,8 +239,11 @@ static void apply_cell_shift(UnitCell *cell, int k, double shift)
/* Apply the given shift to the 'k'th parameter of 'image'. */
-static void apply_shift(struct image *image, int k, double shift)
+static void apply_shift(Crystal *cr, int k, double shift)
{
+ double t;
+ struct image *image = crystal_get_image(cr);
+
switch ( k ) {
case REF_DIV :
@@ -251,7 +256,9 @@ static void apply_shift(struct image *image, int k, double shift)
break;
case REF_R :
- image->profile_radius += shift;
+ t = crystal_get_profile_radius(cr);
+ t += shift;
+ crystal_set_profile_radius(cr, t);
break;
case REF_ASX :
@@ -263,7 +270,7 @@ static void apply_shift(struct image *image, int k, double shift)
case REF_CSX :
case REF_CSY :
case REF_CSZ :
- apply_cell_shift(image->indexed_cell, k, shift);
+ apply_cell_shift(crystal_get_cell(cr), k, shift);
break;
default :
@@ -357,7 +364,7 @@ static gsl_vector *solve_svd(gsl_vector *v, gsl_matrix *M)
/* Perform one cycle of post refinement on 'image' against 'full' */
-static double pr_iterate(struct image *image, const RefList *full)
+static double pr_iterate(Crystal *cr, const RefList *full)
{
gsl_matrix *M;
gsl_vector *v;
@@ -369,7 +376,7 @@ static double pr_iterate(struct image *image, const RefList *full)
double max_shift;
int nref = 0;
- reflections = image->reflections;
+ reflections = crystal_get_reflections(cr);
M = gsl_matrix_calloc(NUM_PARAMS, NUM_PARAMS);
v = gsl_vector_calloc(NUM_PARAMS);
@@ -387,6 +394,7 @@ static double pr_iterate(struct image *image, const RefList *full)
double p;
Reflection *match;
double gradients[NUM_PARAMS];
+ const double osf = crystal_get_osf(cr);
if ( !get_refinable(refl) ) continue;
@@ -398,7 +406,7 @@ static double pr_iterate(struct image *image, const RefList *full)
" is it marked as refinable?\n", ha, ka, la);
continue;
}
- I_full = image->osf * get_intensity(match);
+ I_full = osf * get_intensity(match);
/* Actual measurement of this reflection from this pattern? */
I_partial = get_intensity(refl);
@@ -412,7 +420,7 @@ static double pr_iterate(struct image *image, const RefList *full)
/* Calculate all gradients for this reflection */
for ( k=0; k<NUM_PARAMS; k++ ) {
double gr;
- gr = gradient(image, k, refl, image->profile_radius);
+ gr = gradient(cr, k, refl);
gradients[k] = gr;
}
@@ -429,7 +437,7 @@ static double pr_iterate(struct image *image, const RefList *full)
if ( g > k ) continue;
M_c = gradients[g] * gradients[k];
- M_c *= w * pow(image->osf * I_full, 2.0);
+ M_c *= w * pow(osf * I_full, 2.0);
M_curr = gsl_matrix_get(M, k, g);
gsl_matrix_set(M, k, g, M_curr + M_c);
@@ -450,7 +458,7 @@ static double pr_iterate(struct image *image, const RefList *full)
//STATUS("%i reflections went into the equations.\n", nref);
if ( nref == 0 ) {
- image->pr_dud = 1;
+ crystal_set_user_flag(cr, 1);
gsl_matrix_free(M);
gsl_vector_free(v);
return 0.0;
@@ -462,7 +470,7 @@ static double pr_iterate(struct image *image, const RefList *full)
for ( param=0; param<NUM_PARAMS; param++ ) {
double shift = gsl_vector_get(shifts, param);
- apply_shift(image, param, shift);
+ apply_shift(cr, param, shift);
//STATUS("Shift %i: %e\n", param, shift);
if ( fabs(shift) > max_shift ) max_shift = fabs(shift);
}
@@ -480,7 +488,7 @@ static double pr_iterate(struct image *image, const RefList *full)
}
-static double guide_dev(struct image *image, const RefList *full)
+static double guide_dev(Crystal *cr, const RefList *full)
{
double dev = 0.0;
@@ -488,7 +496,7 @@ static double guide_dev(struct image *image, const RefList *full)
Reflection *refl;
RefListIterator *iter;
- for ( refl = first_refl(image->reflections, &iter);
+ for ( refl = first_refl(crystal_get_reflections(cr), &iter);
refl != NULL;
refl = next_refl(refl, iter) ) {
@@ -508,7 +516,7 @@ static double guide_dev(struct image *image, const RefList *full)
* scale_intensities() might not yet have been called, so the
* full version may not have been calculated yet. */
- G = image->osf;
+ G = crystal_get_osf(cr);
p = get_partiality(refl);
I_partial = get_intensity(refl);
I_full = get_intensity(full_version);
@@ -524,21 +532,21 @@ static double guide_dev(struct image *image, const RefList *full)
}
-void pr_refine(struct image *image, const RefList *full)
+void pr_refine(Crystal *cr, const RefList *full)
{
double max_shift, dev;
int i;
const int verbose = 0;
if ( verbose ) {
- dev = guide_dev(image, full);
+ dev = guide_dev(cr, full);
STATUS("\n"); /* Deal with progress bar */
STATUS("Before iteration: dev = %10.5e\n",
dev);
}
i = 0;
- image->pr_dud = 0;
+ crystal_set_user_flag(cr, 0);
do {
double asx, asy, asz;
@@ -546,15 +554,15 @@ void pr_refine(struct image *image, const RefList *full)
double csx, csy, csz;
double dev;
- cell_get_reciprocal(image->indexed_cell, &asx, &asy, &asz,
+ cell_get_reciprocal(crystal_get_cell(cr), &asx, &asy, &asz,
&bsx, &bsy, &bsz, &csx, &csy, &csz);
- max_shift = pr_iterate(image, full);
+ max_shift = pr_iterate(cr, full);
- update_partialities(image);
+ update_partialities(cr);
if ( verbose ) {
- dev = guide_dev(image, full);
+ dev = guide_dev(cr, full);
STATUS("PR Iteration %2i: max shift = %10.2f"
" dev = %10.5e\n", i+1, max_shift, dev);
}
diff --git a/src/post-refinement.h b/src/post-refinement.h
index 2223dcdf..fe171882 100644
--- a/src/post-refinement.h
+++ b/src/post-refinement.h
@@ -39,6 +39,7 @@
#include "image.h"
#include "utils.h"
+#include "crystal.h"
/* Refineable parameters */
@@ -58,10 +59,10 @@ enum {
};
-extern void pr_refine(struct image *image, const RefList *full);
+extern void pr_refine(Crystal *cr, const RefList *full);
/* Exported so it can be poked by tests/pr_gradient_check */
-extern double gradient(struct image *image, int k, Reflection *refl, double r);
+extern double gradient(Crystal *cr, int k, Reflection *refl);
#endif /* POST_REFINEMENT_H */
diff --git a/src/process_hkl.c b/src/process_hkl.c
index 8bfbea2b..8705ce0f 100644
--- a/src/process_hkl.c
+++ b/src/process_hkl.c
@@ -3,12 +3,12 @@
*
* Assemble and process FEL Bragg intensities
*
- * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY,
- * a research centre of the Helmholtz Association.
+ * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY,
+ * a research centre of the Helmholtz Association.
* Copyright © 2012 Lorenzo Galli
*
* Authors:
- * 2009-2012 Thomas White <taw@physics.org>
+ * 2009-2013 Thomas White <taw@physics.org>
* 2011 Andrew Martin <andrew.martin@desy.de>
* 2012 Lorenzo Galli <lorenzo.galli@desy.de>
*
@@ -48,6 +48,8 @@
#include "stream.h"
#include "reflist.h"
#include "image.h"
+#include "crystal.h"
+#include "thread-pool.h"
static void show_help(const char *s)
@@ -62,8 +64,8 @@ static void show_help(const char *s)
" Default: processed.hkl).\n"
" -y, --symmetry=<sym> Merge according to point group <sym>.\n"
"\n"
-" --start-after=<n> Skip n patterns at the start of the stream.\n"
-" --stop-after=<n> Stop after processing n patterns.\n"
+" --start-after=<n> Skip <n> crystals at the start of the stream.\n"
+" --stop-after=<n> Stop after merging <n> crystals.\n"
" -g, --histogram=<h,k,l> Calculate the histogram of measurements for this\n"
" reflection.\n"
" -z, --hist-parameters Set the range for the histogram and the number of\n"
@@ -170,11 +172,11 @@ static double scale_intensities(RefList *reference, RefList *new,
}
-static int merge_pattern(RefList *model, struct image *new, RefList *reference,
- const SymOpList *sym,
- double *hist_vals, signed int hist_h,
- signed int hist_k, signed int hist_l, int *hist_n,
- int config_nopolar)
+static int merge_crystal(RefList *model, struct image *image, Crystal *cr,
+ RefList *reference, const SymOpList *sym,
+ double *hist_vals, signed int hist_h,
+ signed int hist_k, signed int hist_l, int *hist_n,
+ int config_nopolar)
{
Reflection *refl;
RefListIterator *iter;
@@ -184,17 +186,18 @@ static int merge_pattern(RefList *model, struct image *new, RefList *reference,
double csx, csy, csz;
if ( reference != NULL ) {
- scale = scale_intensities(reference, new->reflections, sym);
+ scale = scale_intensities(reference,
+ crystal_get_reflections(cr), sym);
} else {
scale = 1.0;
}
if ( isnan(scale) ) return 1;
- cell_get_reciprocal(new->indexed_cell, &asx, &asy, &asz,
- &bsx, &bsy, &bsz,
- &csx, &csy, &csz);
+ cell_get_reciprocal(crystal_get_cell(cr), &asx, &asy, &asz,
+ &bsx, &bsy, &bsz,
+ &csx, &csy, &csz);
- for ( refl = first_refl(new->reflections, &iter);
+ for ( refl = first_refl(crystal_get_reflections(cr), &iter);
refl != NULL;
refl = next_refl(refl, iter) )
{
@@ -226,7 +229,7 @@ static int merge_pattern(RefList *model, struct image *new, RefList *reference,
yl = h*asy + k*bsy + l*csy;
zl = h*asz + k*bsz + l*csz;
- ool = 1.0 / new->lambda;
+ ool = 1.0 / image->lambda;
tt = angle_between(0.0, 0.0, 1.0, xl, yl, zl+ool);
phi = atan2(yl, xl);
pa = pow(sin(phi)*sin(tt), 2.0);
@@ -260,60 +263,76 @@ static int merge_pattern(RefList *model, struct image *new, RefList *reference,
}
-static void merge_all(FILE *fh, RefList *model, RefList *reference,
- int config_startafter, int config_stopafter,
- const SymOpList *sym,
- int n_total_patterns,
- double *hist_vals, signed int hist_h,
- signed int hist_k, signed int hist_l,
- int *hist_i, int config_nopolar, int min_measurements)
+static void display_progress(int n_images, int n_crystals, int n_crystals_used)
+{
+ if ( !isatty(STDERR_FILENO) ) return;
+ if ( tcgetpgrp(STDERR_FILENO) != getpgrp() ) return;
+
+ pthread_mutex_lock(&stderr_lock);
+ fprintf(stderr, "\r%i images processed, %i crystals, %i crystals used.",
+ n_images, n_crystals, n_crystals_used);
+ pthread_mutex_unlock(&stderr_lock);
+
+ fflush(stdout);
+}
+
+
+
+static int merge_all(Stream *st, RefList *model, RefList *reference,
+ const SymOpList *sym,
+ double *hist_vals, signed int hist_h,
+ signed int hist_k, signed int hist_l,
+ int *hist_i, int config_nopolar, int min_measurements,
+ int start_after, int stop_after)
{
int rval;
- int n_patterns = 0;
- int n_used = 0;
+ int n_images = 0;
+ int n_crystals = 0;
+ int n_crystals_used = 0;
Reflection *refl;
RefListIterator *iter;
-
- if ( skip_some_files(fh, config_startafter) ) {
- ERROR("Failed to skip first %i files.\n", config_startafter);
- return;
- }
+ int n_crystals_seen = 0;
do {
struct image image;
+ int i;
image.det = NULL;
/* Get data from next chunk */
- rval = read_chunk(fh, &image);
+ rval = read_chunk(st, &image);
if ( rval ) break;
- n_patterns++;
+ n_images++;
- if ( (image.reflections != NULL) && (image.indexed_cell) ) {
+ for ( i=0; i<image.n_crystals; i++ ) {
int r;
+ Crystal *cr = image.crystals[i];
- r = merge_pattern(model, &image, reference, sym,
+ n_crystals_seen++;
+ if ( n_crystals_seen <= start_after ) continue;
+
+ n_crystals++;
+ r = merge_crystal(model, &image, cr, reference, sym,
hist_vals, hist_h, hist_k, hist_l,
hist_i, config_nopolar);
- if ( r == 0 ) n_used++;
+ if ( r == 0 ) n_crystals_used++;
+
+ crystal_free(cr);
+
+ if ( n_crystals_used == stop_after ) break;
}
free(image.filename);
- reflist_free(image.reflections);
image_feature_list_free(image.features);
- cell_free(image.indexed_cell);
- progress_bar(n_patterns, n_total_patterns-config_startafter,
- "Merging");
+ display_progress(n_images, n_crystals_seen, n_crystals_used);
- if ( config_stopafter ) {
- if ( n_patterns == config_stopafter ) break;
- }
+ if ( (stop_after>0) && (n_crystals_used == stop_after) ) break;
} while ( rval == 0 );
@@ -339,7 +358,7 @@ static void merge_all(FILE *fh, RefList *model, RefList *reference,
}
- STATUS("%i of the patterns could be used.\n", n_used);
+ return 0;
}
@@ -348,14 +367,11 @@ int main(int argc, char *argv[])
int c;
char *filename = NULL;
char *output = NULL;
- FILE *fh;
+ Stream *st;
RefList *model;
int config_maxonly = 0;
- int config_startafter = 0;
- int config_stopafter = 0;
int config_sum = 0;
int config_scale = 0;
- unsigned int n_total_patterns;
char *sym_str = NULL;
SymOpList *sym;
char *histo = NULL;
@@ -369,6 +385,9 @@ int main(int argc, char *argv[])
int config_nopolar = 0;
char *rval;
int min_measurements = 2;
+ int r;
+ int start_after = 0;
+ int stop_after = 0;
/* Long options */
const struct option longopts[] = {
@@ -377,8 +396,8 @@ int main(int argc, char *argv[])
{"output", 1, NULL, 'o'},
{"max-only", 0, &config_maxonly, 1},
{"output-every", 1, NULL, 'e'},
- {"stop-after", 1, NULL, 's'},
- {"start-after", 1, NULL, 'f'},
+ {"start-after", 1, NULL, 's'},
+ {"stop-after", 1, NULL, 'f'},
{"sum", 0, &config_sum, 1},
{"scale", 0, &config_scale, 1},
{"no-polarisation", 0, &config_nopolar, 1},
@@ -391,7 +410,7 @@ int main(int argc, char *argv[])
};
/* Short options */
- while ((c = getopt_long(argc, argv, "hi:e:o:y:g:f:b:z:",
+ while ((c = getopt_long(argc, argv, "hi:e:o:y:g:s:f:z:",
longopts, NULL)) != -1) {
switch (c) {
@@ -409,11 +428,21 @@ int main(int argc, char *argv[])
break;
case 's' :
- config_stopafter = atoi(optarg);
+ errno = 0;
+ start_after = strtod(optarg, &rval);
+ if ( *rval != '\0' ) {
+ ERROR("Invalid value for --start-after.\n");
+ return 1;
+ }
break;
case 'f' :
- config_startafter = atoi(optarg);
+ errno = 0;
+ stop_after = strtod(optarg, &rval);
+ if ( *rval != '\0' ) {
+ ERROR("Invalid value for --stop-after.\n");
+ return 1;
+ }
break;
case 'y' :
@@ -465,25 +494,11 @@ int main(int argc, char *argv[])
free(sym_str);
/* Open the data stream */
- if ( strcmp(filename, "-") == 0 ) {
- fh = stdin;
- } else {
- fh = fopen(filename, "r");
- }
- free(filename);
- if ( fh == NULL ) {
- ERROR("Failed to open input file\n");
- return 1;
- }
-
- /* Count the number of patterns in the file */
- n_total_patterns = count_patterns(fh);
- if ( n_total_patterns == 0 ) {
- ERROR("No patterns to process.\n");
+ st = open_stream_for_read(filename);
+ if ( st == NULL ) {
+ ERROR("Failed to open stream.\n");
return 1;
}
- STATUS("There are %i patterns to process\n", n_total_patterns);
- rewind(fh);
model = reflist_new();
@@ -497,7 +512,8 @@ int main(int argc, char *argv[])
return 1;
}
- space_for_hist = n_total_patterns * num_equivs(sym, NULL);
+ /* FIXME: This array must grow as necessary. */
+ space_for_hist = 0 * num_equivs(sym, NULL);
hist_vals = malloc(space_for_hist * sizeof(double));
free(histo);
STATUS("Histogramming %i %i %i -> ", hist_h, hist_k, hist_l);
@@ -529,40 +545,48 @@ int main(int argc, char *argv[])
}
hist_i = 0;
- merge_all(fh, model, NULL, config_startafter, config_stopafter,
- sym, n_total_patterns, hist_vals, hist_h, hist_k, hist_l,
- &hist_i, config_nopolar, min_measurements);
- if ( ferror(fh) ) {
- ERROR("Stream read error.\n");
+ r = merge_all(st, model, NULL, sym, hist_vals, hist_h, hist_k, hist_l,
+ &hist_i, config_nopolar, min_measurements,
+ start_after, stop_after);
+ fprintf(stderr, "\n");
+ if ( r ) {
+ ERROR("Error while reading stream.\n");
return 1;
}
- rewind(fh);
if ( config_scale ) {
RefList *reference;
- STATUS("Extra pass for scaling...\n");
+ if ( rewind_stream(st) ) {
- reference = copy_reflist(model);
+ ERROR("Couldn't rewind stream - scaling cannot be "
+ "performed.\n");
- reflist_free(model);
- model = reflist_new();
+ } else {
- rewind(fh);
+ int r;
- merge_all(fh, model, reference,
- config_startafter, config_stopafter, sym,
- n_total_patterns,
- hist_vals, hist_h, hist_k, hist_l, &hist_i,
- config_nopolar, min_measurements);
+ STATUS("Extra pass for scaling...\n");
- if ( ferror(fh) ) {
- ERROR("Stream read error.\n");
- return 1;
- }
+ reference = copy_reflist(model);
+
+ reflist_free(model);
+ model = reflist_new();
+
+ r = merge_all(st, model, reference, sym,
+ hist_vals, hist_h, hist_k, hist_l, &hist_i,
+ config_nopolar, min_measurements,
+ start_after, stop_after);
+ fprintf(stderr, "\n");
+ if ( r ) {
+ ERROR("Error while reading stream.\n");
+ return 1;
+ }
- reflist_free(reference);
+ reflist_free(reference);
+
+ }
}
@@ -579,7 +603,7 @@ int main(int argc, char *argv[])
write_reflist(output, model);
- fclose(fh);
+ close_stream(st);
free(sym);
reflist_free(model);
diff --git a/src/process_image.c b/src/process_image.c
new file mode 100644
index 00000000..010054c8
--- /dev/null
+++ b/src/process_image.c
@@ -0,0 +1,231 @@
+/*
+ * process_image.c
+ *
+ * The processing pipeline for one image
+ *
+ * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY,
+ * a research centre of the Helmholtz Association.
+ *
+ * Authors:
+ * 2010-2013 Thomas White <taw@physics.org>
+ *
+ * This file is part of CrystFEL.
+ *
+ * CrystFEL is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * CrystFEL is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with CrystFEL. If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdlib.h>
+#include <hdf5.h>
+#include <gsl/gsl_errno.h>
+
+#include "utils.h"
+#include "hdf5-file.h"
+#include "index.h"
+#include "peaks.h"
+#include "detector.h"
+#include "filters.h"
+#include "thread-pool.h"
+#include "beam-parameters.h"
+#include "geometry.h"
+#include "stream.h"
+#include "reflist-utils.h"
+#include "process_image.h"
+
+
+void process_image(const struct index_args *iargs, struct pattern_args *pargs,
+ Stream *st, int cookie)
+{
+ float *data_for_measurement;
+ size_t data_size;
+ int check;
+ struct hdfile *hdfile;
+ struct image image;
+ int i;
+ char filename[1024];
+
+ /* Prefix to jump out of temporary folder */
+ if ( pargs->filename[0] != '/' ) {
+ snprintf(filename, 1023, "../../%s", pargs->filename);
+ } else {
+ snprintf(filename, 1023, "%s", pargs->filename);
+ }
+
+ image.features = NULL;
+ image.data = NULL;
+ image.flags = NULL;
+ image.copyme = iargs->copyme;
+ image.id = cookie;
+ image.filename = pargs->filename; /* Relative to top level */
+ image.beam = iargs->beam;
+ image.det = iargs->det;
+ image.crystals = NULL;
+ image.n_crystals = 0;
+
+ hdfile = hdfile_open(filename); /* Relative to temporary folder */
+ if ( hdfile == NULL ) return;
+
+ if ( iargs->element != NULL ) {
+
+ int r;
+ r = hdfile_set_image(hdfile, iargs->element);
+ if ( r ) {
+ ERROR("Couldn't select path '%s'\n", iargs->element);
+ hdfile_close(hdfile);
+ return;
+ }
+
+ } else {
+
+ int r;
+ r = hdfile_set_first_image(hdfile, "/");
+ if ( r ) {
+ ERROR("Couldn't select first path\n");
+ hdfile_close(hdfile);
+ return;
+ }
+
+ }
+
+ check = hdf5_read(hdfile, &image, 1);
+ if ( check ) {
+ hdfile_close(hdfile);
+ return;
+ }
+
+ if ( (image.width != image.det->max_fs + 1 )
+ || (image.height != image.det->max_ss + 1))
+ {
+ ERROR("Image size doesn't match geometry size"
+ " - rejecting image.\n");
+ ERROR("Image size: %i,%i. Geometry size: %i,%i\n",
+ image.width, image.height,
+ image.det->max_fs + 1, image.det->max_ss + 1);
+ hdfile_close(hdfile);
+ return;
+ }
+
+ fill_in_values(image.det, hdfile);
+ fill_in_beam_parameters(image.beam, hdfile);
+
+ image.lambda = ph_en_to_lambda(eV_to_J(image.beam->photon_energy));
+
+ if ( (image.beam->photon_energy < 0.0) || (image.lambda > 1000) ) {
+ /* Error message covers a silly value in the beam file or in
+ * the HDF5 file. */
+ ERROR("Nonsensical wavelength (%e m or %e eV) value for %s.\n",
+ image.lambda, image.beam->photon_energy, image.filename);
+ hdfile_close(hdfile);
+ return;
+ }
+
+ if ( iargs->cmfilter ) filter_cm(&image);
+
+ /* Take snapshot of image after CM subtraction but before
+ * the aggressive noise filter. */
+ data_size = image.width * image.height * sizeof(float);
+ data_for_measurement = malloc(data_size);
+
+ if ( iargs->noisefilter ) {
+ filter_noise(&image, data_for_measurement);
+ } else {
+ memcpy(data_for_measurement, image.data, data_size);
+ }
+
+ switch ( iargs->peaks ) {
+
+ case PEAK_HDF5:
+ // Get peaks from HDF5
+ if (get_peaks(&image, hdfile,
+ iargs->hdf5_peak_path)) {
+ ERROR("Failed to get peaks from HDF5 file.\n");
+ }
+ if ( !iargs->no_revalidate ) {
+ validate_peaks(&image, iargs->min_int_snr,
+ iargs->ir_inn, iargs->ir_mid,
+ iargs->ir_out, iargs->use_saturated);
+ }
+ break;
+
+ case PEAK_ZAEF:
+ search_peaks(&image, iargs->threshold,
+ iargs->min_gradient, iargs->min_snr,
+ iargs->ir_inn, iargs->ir_mid, iargs->ir_out,
+ iargs->use_saturated);
+ break;
+
+ }
+
+ /* Get rid of noise-filtered version at this point
+ * - it was strictly for the purposes of peak detection. */
+ free(image.data);
+ image.data = data_for_measurement;
+
+ /* Index the pattern */
+ index_pattern(&image, iargs->indm, iargs->ipriv);
+
+ pargs->n_crystals = image.n_crystals;
+
+ /* Default beam parameters */
+ image.div = image.beam->divergence;
+ image.bw = image.beam->bandwidth;
+
+ /* Integrate each crystal's diffraction spots */
+ for ( i=0; i<image.n_crystals; i++ ) {
+
+ RefList *reflections;
+
+ /* Set default crystal parameter(s) */
+ crystal_set_profile_radius(image.crystals[i],
+ image.beam->profile_radius);
+
+ if ( iargs->integrate_found ) {
+ reflections = select_intersections(&image,
+ image.crystals[i]);
+ } else {
+ reflections = find_intersections(&image,
+ image.crystals[i]);
+ }
+
+ crystal_set_reflections(image.crystals[i], reflections);
+
+ }
+
+ /* Integrate all the crystals at once - need all the crystals so that
+ * overlaps can be detected. */
+ integrate_reflections(&image, iargs->closer,
+ iargs->bgsub,
+ iargs->min_int_snr,
+ iargs->ir_inn,
+ iargs->ir_mid,
+ iargs->ir_out,
+ iargs->integrate_saturated,
+ iargs->res_cutoff);
+
+ write_chunk(st, &image, hdfile,
+ iargs->stream_peaks, iargs->stream_refls);
+
+ for ( i=0; i<image.n_crystals; i++ ) {
+ crystal_free(image.crystals[i]);
+ }
+
+ free(image.data);
+ if ( image.flags != NULL ) free(image.flags);
+ image_feature_list_free(image.features);
+ hdfile_close(hdfile);
+}
diff --git a/src/process_image.h b/src/process_image.h
new file mode 100644
index 00000000..f2f034d9
--- /dev/null
+++ b/src/process_image.h
@@ -0,0 +1,94 @@
+/*
+ * process_image.h
+ *
+ * The processing pipeline for one image
+ *
+ * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY,
+ * a research centre of the Helmholtz Association.
+ *
+ * Authors:
+ * 2010-2013 Thomas White <taw@physics.org>
+ *
+ * This file is part of CrystFEL.
+ *
+ * CrystFEL is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * CrystFEL is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with CrystFEL. If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+#ifndef PROCESS_IMAGE_H
+#define PROCESS_IMAGE_H
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+
+enum {
+ PEAK_ZAEF,
+ PEAK_HDF5,
+};
+
+
+/* Information about the indexing process which is common to all patterns */
+struct index_args
+{
+ UnitCell *cell;
+ int cmfilter;
+ int noisefilter;
+ int satcorr;
+ int closer;
+ int bgsub;
+ float threshold;
+ float min_gradient;
+ float min_snr;
+ double min_int_snr;
+ struct detector *det;
+ IndexingMethod *indm;
+ IndexingPrivate **ipriv;
+ int peaks; /* Peak detection method */
+ float tols[4];
+ struct beam_params *beam;
+ char *element;
+ char *hdf5_peak_path;
+ float ir_inn;
+ float ir_mid;
+ float ir_out;
+ struct copy_hdf5_field *copyme;
+ int integrate_saturated;
+ int use_saturated;
+ int no_revalidate;
+ int integrate_found;
+ int stream_peaks;
+ int stream_refls;
+ int res_cutoff;
+};
+
+
+/* Information about the indexing process for one pattern */
+struct pattern_args
+{
+ /* "Input" */
+ char *filename;
+
+ /* "Output" */
+ int n_crystals;
+};
+
+
+extern void process_image(const struct index_args *iargs,
+ struct pattern_args *pargs, Stream *st,
+ int cookie);
+
+
+#endif /* PROCESS_IMAGEs_H */
diff --git a/src/scaling-report.c b/src/scaling-report.c
index ba425b96..a02a1796 100644
--- a/src/scaling-report.c
+++ b/src/scaling-report.c
@@ -3,11 +3,11 @@
*
* Write a nice PDF of scaling parameters
*
- * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY,
- * a research centre of the Helmholtz Association.
+ * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY,
+ * a research centre of the Helmholtz Association.
*
* Authors:
- * 2010-2012 Thomas White <taw@physics.org>
+ * 2010-2013 Thomas White <taw@physics.org>
*
* This file is part of CrystFEL.
*
@@ -183,7 +183,7 @@ static void plot_point(cairo_t *cr, double g_width, double g_height,
}
-static void partiality_graph(cairo_t *cr, const struct image *images, int n,
+static void partiality_graph(cairo_t *cr, Crystal **crystals, int n,
RefList *full)
{
const double g_width = 200.0;
@@ -219,7 +219,7 @@ static void partiality_graph(cairo_t *cr, const struct image *images, int n,
num_nondud = 0;
for ( i=0; i<n; i++ ) {
- if ( images[i].pr_dud ) continue;
+ if ( crystal_get_user_flag(crystals[i]) ) continue;
num_nondud++;
}
@@ -229,10 +229,13 @@ static void partiality_graph(cairo_t *cr, const struct image *images, int n,
Reflection *refl;
RefListIterator *iter;
+ Crystal *cryst;
- if ( images[i].pr_dud ) continue;
+ cryst = crystals[i];
- for ( refl = first_refl(images[i].reflections, &iter);
+ if ( crystal_get_user_flag(cryst) ) continue;
+
+ for ( refl = first_refl(crystal_get_reflections(cryst), &iter);
refl != NULL;
refl = next_refl(refl, iter) )
{
@@ -249,7 +252,7 @@ static void partiality_graph(cairo_t *cr, const struct image *images, int n,
if ( get_redundancy(f) < 2 ) continue;
Ipart = get_intensity(refl);
- Ifull = images[i].osf * get_intensity(f);
+ Ifull = crystal_get_osf(cryst) * get_intensity(f);
//if ( Ifull < 10 ) continue; /* FIXME: Ugh */
@@ -267,7 +270,7 @@ static void partiality_graph(cairo_t *cr, const struct image *images, int n,
double esd_pobs, esd_Ip, esd_If;
esd_Ip = get_esd_intensity(refl);
esd_If = get_esd_intensity(f);
- esd_If *= images[i].osf;
+ esd_If *= crystal_get_osf(cryst);
esd_pobs = pow(esd_Ip/Ipart, 2.0);
esd_pobs += pow(esd_If/Ifull, 2.0);
esd_pobs = sqrt(esd_pobs);
@@ -313,8 +316,8 @@ static void partiality_graph(cairo_t *cr, const struct image *images, int n,
}
-static void partiality_histogram(cairo_t *cr, const struct image *images,
- int n, RefList *full, int calc, int backwards)
+static void partiality_histogram(cairo_t *cr, Crystal **crystals, int n,
+ RefList *full, int calc, int backwards)
{
int f_max;
int i, b;
@@ -344,10 +347,11 @@ static void partiality_histogram(cairo_t *cr, const struct image *images,
Reflection *refl;
RefListIterator *iter;
+ Crystal *cryst = crystals[i];
- if ( images[i].pr_dud ) continue;
+ if ( crystal_get_user_flag(cryst) ) continue;
- for ( refl = first_refl(images[i].reflections, &iter);
+ for ( refl = first_refl(crystal_get_reflections(cryst), &iter);
refl != NULL;
refl = next_refl(refl, iter) )
{
@@ -364,7 +368,7 @@ static void partiality_histogram(cairo_t *cr, const struct image *images,
Ipart = get_intensity(refl);
Ifull = get_intensity(f);
- pobs = Ipart/(images[i].osf*Ifull);
+ pobs = Ipart/(crystal_get_osf(cryst)*Ifull);
pcalc = get_partiality(refl);
if ( calc ) {
@@ -426,8 +430,8 @@ static void partiality_histogram(cairo_t *cr, const struct image *images,
}
-static void scale_factor_histogram(cairo_t *cr, const struct image *images,
- int n, const char *title)
+static void scale_factor_histogram(cairo_t *cr, Crystal **crystals, int n,
+ const char *title)
{
int f_max;
int i, b;
@@ -451,8 +455,8 @@ static void scale_factor_histogram(cairo_t *cr, const struct image *images,
osf_max = 0.0;
for ( i=0; i<n; i++ ) {
- double osf = images[i].osf;
- if ( images[i].pr_dud ) continue;
+ double osf = crystal_get_osf(crystals[i]);
+ if ( crystal_get_user_flag(crystals[i]) ) continue;
if ( osf > osf_max ) osf_max = osf;
}
osf_max = ceil(osf_max+osf_max/10000.0);
@@ -473,9 +477,9 @@ static void scale_factor_histogram(cairo_t *cr, const struct image *images,
for ( i=0; i<n; i++ ) {
- double osf = images[i].osf;
+ double osf = crystal_get_osf(crystals[i]);
- if ( images[i].pr_dud ) continue;
+ if ( crystal_get_user_flag(crystals[i]) ) continue;
for ( b=0; b<nbins; b++ ) {
if ( (osf >= osf_low[b])
@@ -544,8 +548,8 @@ static void scale_factor_histogram(cairo_t *cr, const struct image *images,
}
-static void intensity_histogram(cairo_t *cr, const struct image *images,
- int n, RefList *full,
+static void intensity_histogram(cairo_t *cr, Crystal **crystals, int n,
+ RefList *full,
signed int h, signed int k, signed int l)
{
int f_max;
@@ -582,12 +586,13 @@ static void intensity_histogram(cairo_t *cr, const struct image *images,
for ( i=0; i<n; i++ ) {
double osf;
+ Crystal *cryst = crystals[i];
- if ( images[i].pr_dud ) continue;
+ if ( crystal_get_user_flag(cryst) ) continue;
- osf = images[i].osf;
+ osf = crystal_get_osf(cryst);
- for ( f = find_refl(images[i].reflections, h, k, l);
+ for ( f = find_refl(crystal_get_reflections(cryst), h, k, l);
f != NULL;
f = next_found_refl(f) )
{
@@ -616,12 +621,13 @@ static void intensity_histogram(cairo_t *cr, const struct image *images,
for ( i=0; i<n; i++ ) {
double osf;
+ Crystal *cryst = crystals[i];
- if ( images[i].pr_dud ) continue;
+ if ( crystal_get_user_flag(cryst) ) continue;
- osf = images[i].osf;
+ osf = crystal_get_osf(cryst);
- for ( f = find_refl(images[i].reflections, h, k, l);
+ for ( f = find_refl(crystal_get_reflections(cryst), h, k, l);
f != NULL;
f = next_found_refl(f) )
{
@@ -727,6 +733,13 @@ static void find_most_sampled_reflections(RefList *list, int n, signed int *h,
Reflection *refl;
RefListIterator *iter;
int *samples;
+ int j;
+
+ for ( j=0; j<n; j++ ) {
+ h[j] = 0;
+ k[j] = 0;
+ l[j] = 0;
+ }
samples = calloc(n, sizeof(int));
@@ -771,7 +784,7 @@ static void find_most_sampled_reflections(RefList *list, int n, signed int *h,
-SRContext *sr_titlepage(struct image *images, int n,
+SRContext *sr_titlepage(Crystal **crystals, int n,
const char *filename, const char *stream_filename,
const char *cmdline)
{
@@ -805,7 +818,7 @@ SRContext *sr_titlepage(struct image *images, int n,
}
-void sr_iteration(SRContext *sr, int iteration, struct image *images, int n,
+void sr_iteration(SRContext *sr, int iteration, Crystal **crystals, int n,
RefList *full)
{
int i;
@@ -822,7 +835,7 @@ void sr_iteration(SRContext *sr, int iteration, struct image *images, int n,
cairo_save(sr->cr);
cairo_translate(sr->cr, 480.0, 350.0);
- scale_factor_histogram(sr->cr, images, n,
+ scale_factor_histogram(sr->cr, crystals, n,
"Distribution of overall scale factors");
cairo_restore(sr->cr);
@@ -830,7 +843,7 @@ void sr_iteration(SRContext *sr, int iteration, struct image *images, int n,
cairo_save(sr->cr);
cairo_translate(sr->cr, 70.0, 330.0);
- partiality_graph(sr->cr, images, n, full);
+ partiality_graph(sr->cr, crystals, n, full);
cairo_save(sr->cr);
cairo_move_to(sr->cr, 0.0, 0.0);
@@ -841,7 +854,7 @@ void sr_iteration(SRContext *sr, int iteration, struct image *images, int n,
cairo_stroke(sr->cr);
cairo_set_dash(sr->cr, NULL, 0, 0.0);
cairo_translate(sr->cr, 0.0, -150.0);
- partiality_histogram(sr->cr, images, n, full, 1, 0);
+ partiality_histogram(sr->cr, crystals, n, full, 1, 0);
cairo_restore(sr->cr);
cairo_save(sr->cr);
@@ -854,7 +867,7 @@ void sr_iteration(SRContext *sr, int iteration, struct image *images, int n,
cairo_set_dash(sr->cr, NULL, 0, 0.0);
cairo_translate(sr->cr, 230.0, 200.0);
cairo_rotate(sr->cr, -M_PI_2);
- partiality_histogram(sr->cr, images, n, full, 0, 1);
+ partiality_histogram(sr->cr, crystals, n, full, 0, 1);
cairo_restore(sr->cr);
cairo_restore(sr->cr);
@@ -873,7 +886,7 @@ void sr_iteration(SRContext *sr, int iteration, struct image *images, int n,
cairo_save(sr->cr);
cairo_translate(sr->cr, 400.0+140.0*x, 60.0+80.0*y);
- intensity_histogram(sr->cr, images, n, full,
+ intensity_histogram(sr->cr, crystals, n, full,
sr->ms_h[i], sr->ms_k[i], sr->ms_l[i]);
cairo_restore(sr->cr);
diff --git a/src/scaling-report.h b/src/scaling-report.h
index b9ac3fb7..5b153377 100644
--- a/src/scaling-report.h
+++ b/src/scaling-report.h
@@ -40,25 +40,25 @@ typedef struct _srcontext SRContext; /* Opaque */
#ifdef HAVE_CAIRO
-extern SRContext *sr_titlepage(struct image *images, int n,
+extern SRContext *sr_titlepage(Crystal **crystals, int n,
const char *filename,
const char *stream_filename,
const char *cmdline);
-extern void sr_iteration(SRContext *sr, int iteration, struct image *images,
- int n, RefList *full);
+extern void sr_iteration(SRContext *sr, int iteration,
+ Crystal **crystals, int n, RefList *full);
extern void sr_finish(SRContext *sr);
#else
-SRContext *sr_titlepage(struct image *images, int n, const char *filename,
+SRContext *sr_titlepage(Crystal **crystals, int n, const char *filename,
const char *stream_filename, const char *cmdline)
{
return NULL;
}
-void sr_iteration(SRContext *sr, int iteration, struct image *images, int n,
+void sr_iteration(SRContext *sr, int iteration, Crystal **crystals, int n,
RefList *full)
{
}
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-16 16:04:54 +0000