aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2013-02-05 16:21:40 +0100
committerThomas White <taw@physics.org>2013-02-05 16:21:40 +0100
commit2ce85a95d86e350b785d206405e97d7317672188 (patch)
tree5cecce161733acf1a276e5a0b55e5c22f5ffb20d
parent29cca07716b48f9e433087f5dbb202165b1897e1 (diff)
Indexing pipeline - "done"!
-rw-r--r--libcrystfel/src/index.c5
-rw-r--r--libcrystfel/src/index.h2
-rw-r--r--libcrystfel/src/stream.c24
-rw-r--r--libcrystfel/src/stream.h2
-rw-r--r--src/im-sandbox.c129
-rw-r--r--src/im-sandbox.h27
-rw-r--r--src/indexamajig.c128
7 files changed, 147 insertions, 170 deletions
diff --git a/libcrystfel/src/index.c b/libcrystfel/src/index.c
index 2370e26e..821be39a 100644
--- a/libcrystfel/src/index.c
+++ b/libcrystfel/src/index.c
@@ -241,7 +241,7 @@ static IndexingMethod set_axes(IndexingMethod a)
}
-IndexingMethod *build_indexer_list(const char *str, int *need_cell)
+IndexingMethod *build_indexer_list(const char *str)
{
int n, i;
char **methods;
@@ -249,9 +249,6 @@ IndexingMethod *build_indexer_list(const char *str, int *need_cell)
int tmp;
int nmeth = 0;
- if ( need_cell == NULL ) need_cell = &tmp;
- *need_cell = 0;
-
n = assplode(str, ",-", &methods, ASSPLODE_NONE);
list = malloc((n+1)*sizeof(IndexingMethod));
diff --git a/libcrystfel/src/index.h b/libcrystfel/src/index.h
index e0dcb8e0..093e4ad6 100644
--- a/libcrystfel/src/index.h
+++ b/libcrystfel/src/index.h
@@ -82,7 +82,7 @@ typedef enum {
**/
typedef void *IndexingPrivate;
-extern IndexingMethod *build_indexer_list(const char *str, int *need_cell);
+extern IndexingMethod *build_indexer_list(const char *str);
extern IndexingPrivate **prepare_indexing(IndexingMethod *indm, UnitCell *cell,
const char *filename,
diff --git a/libcrystfel/src/stream.c b/libcrystfel/src/stream.c
index 31fb8527..113dc359 100644
--- a/libcrystfel/src/stream.c
+++ b/libcrystfel/src/stream.c
@@ -38,6 +38,9 @@
#include <stdio.h>
#include <string.h>
#include <assert.h>
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <fcntl.h>
#include "cell.h"
#include "utils.h"
@@ -236,6 +239,8 @@ void write_chunk(Stream *st, struct image *i, struct hdfile *hdfile,
}
fprintf(st->fh, CHUNK_END_MARKER"\n\n");
+
+ fflush(st->fh);
}
@@ -497,14 +502,14 @@ Stream *open_stream_for_read(const char *filename)
}
-Stream *open_stream_for_write(const char *filename)
+Stream *open_stream_fd_for_write(int fd)
{
Stream *st;
st = malloc(sizeof(struct _stream));
if ( st == NULL ) return NULL;
- st->fh = fopen(filename, "w");
+ st->fh = fdopen(fd, "w");
if ( st->fh == NULL ) {
free(st);
return NULL;
@@ -520,6 +525,21 @@ Stream *open_stream_for_write(const char *filename)
}
+
+Stream *open_stream_for_write(const char *filename)
+{
+ int fd;
+
+ fd = open(filename, O_CREAT | O_TRUNC | O_WRONLY, S_IRUSR | S_IWUSR);
+ if ( fd == -1 ) {
+ ERROR("Failed to open stream.\n");
+ return NULL;
+ }
+
+ return open_stream_fd_for_write(fd);
+}
+
+
void close_stream(Stream *st)
{
fclose(st->fh);
diff --git a/libcrystfel/src/stream.h b/libcrystfel/src/stream.h
index 9340ca0d..4d561cd2 100644
--- a/libcrystfel/src/stream.h
+++ b/libcrystfel/src/stream.h
@@ -42,6 +42,8 @@ typedef struct _stream Stream;
extern Stream *open_stream_for_read(const char *filename);
extern Stream *open_stream_for_write(const char *filename);
+extern Stream *open_stream_fd_for_write(int fd);
+
extern void close_stream(Stream *st);
extern void write_chunk(Stream *st, struct image *image, struct hdfile *hdfile,
diff --git a/src/im-sandbox.c b/src/im-sandbox.c
index c8cceb2b..825f2d6a 100644
--- a/src/im-sandbox.c
+++ b/src/im-sandbox.c
@@ -73,14 +73,27 @@
#define STATS_EVERY_N_SECONDS (5)
+/* Information about the indexing process for one pattern */
+struct pattern_args
+{
+ /* "Input" */
+ char *filename;
+
+ /* "Output" */
+ int n_crystals;
+};
+
+
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;
@@ -164,31 +177,29 @@ 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,
+ struct pattern_args *pargs, Stream *st,
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;
+ int i;
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;
+ image.crystals = NULL;
+ image.n_crystals = 0;
hdfile = hdfile_open(image.filename);
if ( hdfile == NULL ) return;
@@ -291,28 +302,34 @@ static void process_image(const struct index_args *iargs,
free(image.data);
image.data = data_for_measurement;
- /* Calculate orientation matrix (by magic) */
+ /* Index the pattern */
+ index_pattern(&image, indm, iargs->ipriv);
+
+ /* Default beam parameters */
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);
+ /* Integrate each crystal's diffraction spots */
+ for ( i=0; i<image.n_crystals; i++ ) {
- if ( image.indexed_cell != NULL ) {
+ RefList *reflections;
- pargs->indexable = 1;
+ /* Set default crystal parameter(s) */
+ crystal_set_profile_radius(image.crystals[i],
+ image.beam->profile_radius);
if ( iargs->integrate_found ) {
- image.reflections = select_intersections(&image,
- image.indexed_cell);
+ reflections = select_intersections(&image,
+ image.crystals[i]);
} else {
- image.reflections = find_intersections(&image,
- image.indexed_cell);
+ reflections = find_intersections(&image,
+ image.crystals[i]);
}
- if (image.reflections != NULL) {
+ crystal_set_reflections(image.crystals[i], reflections);
+
+ if ( reflections != NULL ) {
+
integrate_reflections(&image,
iargs->config_closer,
iargs->config_bgsub,
@@ -323,18 +340,15 @@ static void process_image(const struct index_args *iargs,
iargs->integrate_saturated);
}
- } else {
- image.reflections = NULL;
}
- write_chunk(ofh, &image, hdfile, iargs->stream_flags);
- fprintf(ofh, "END\n");
- fflush(ofh);
+ write_chunk(st, &image, hdfile, iargs->include_peaks,
+ iargs->include_reflections);
- /* Only free cell if found */
- cell_free(image.indexed_cell);
+ for ( i=0; i<image.n_crystals; i++ ) {
+ crystal_free(image.crystals[i]);
+ }
- reflist_free(image.reflections);
free(image.data);
if ( image.flags != NULL ) free(image.flags);
image_feature_list_free(image.features);
@@ -343,7 +357,7 @@ static void process_image(const struct index_args *iargs,
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 +403,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 +420,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);
@@ -577,7 +591,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;
@@ -614,10 +628,10 @@ 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(sb->stream_pipe_write[slot]);
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]);
@@ -697,7 +711,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)
+ Stream *st)
{
int i;
int allDone;
@@ -712,13 +726,14 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix,
return;
}
- sb->n_indexable = 0;
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);
@@ -895,8 +910,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 +950,19 @@ 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,
- sb->n_processed - sb->n_processed_last_stats);
+ STATUS("Total so far: %i images processed, "
+ "%i had crystals, %i crystals overall. "
+ "Since the last message: %i images processed,"
+ "%i had crystals, %i crystals overall.\n",
+ sb->n_processed, sb->n_hadcrystals,
+ sb->n_crystals,
+ sb->n_processed - sb->n_processed_last_stats,
+ sb->n_hadcrystals - sb->n_hadcrystals_last_stats,
+ sb->n_crystals - sb->n_crystals_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;
}
@@ -975,9 +1001,8 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix,
free(sb->pids);
free(sb->running);
- if ( ofh != stdout ) fclose(ofh);
-
- 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);
}
diff --git a/src/im-sandbox.h b/src/im-sandbox.h
index 872f84ad..176bfde0 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,6 +31,8 @@
*
*/
+#include "stream.h"
+
enum {
PEAK_ZAEF,
PEAK_HDF5,
@@ -44,10 +46,8 @@ struct index_args
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;
@@ -57,7 +57,6 @@ struct index_args
IndexingMethod *indm;
IndexingPrivate **ipriv;
int peaks; /* Peak detection method */
- int cellr;
float tols[4];
struct beam_params *beam;
char *element;
@@ -70,21 +69,11 @@ struct index_args
int use_saturated;
int no_revalidate;
int integrate_found;
+ int include_peaks;
+ int include_reflections;
};
-
-
-/* Information about the indexing process for one pattern */
-struct pattern_args
-{
- /* "Input" */
- char *filename;
-
- /* "Output" */
- int indexable;
-};
-
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, Stream *st);
diff --git a/src/indexamajig.c b/src/indexamajig.c
index 01377de3..3076ea86 100644
--- a/src/indexamajig.c
+++ b/src/indexamajig.c
@@ -174,37 +174,13 @@ static void show_help(const char *s)
}
-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;
+ Stream *st;
char *rval = NULL;
int config_noindex = 0;
int config_cmfilter = 0;
@@ -213,7 +189,6 @@ int main(int argc, char *argv[])
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;
@@ -224,17 +199,13 @@ int main(int argc, char *argv[])
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;
@@ -243,7 +214,6 @@ int main(int argc, char *argv[])
struct index_args iargs;
struct beam_params *beam = NULL;
char *element = NULL;
- int stream_flags = STREAM_INTEGRATED;
char *hdf5_peak_path = NULL;
struct copy_hdf5_field *copyme;
char *intrad = NULL;
@@ -255,6 +225,10 @@ int main(int argc, char *argv[])
int no_revalidate = 0;
int integrate_found = 0;
+ /* For ease of upgrading from old method */
+ char *scellr = NULL;
+ int config_insane = 0;
+
copyme = new_copy_hdf5_field_list();
if ( copyme == NULL ) {
ERROR("Couldn't allocate HDF5 field list.\n");
@@ -290,7 +264,6 @@ int main(int argc, char *argv[])
{"peaks", 1, NULL, 2},
{"cell-reduction", 1, NULL, 3},
{"min-gradient", 1, NULL, 4},
- {"record", 1, NULL, 5},
{"cpus", 1, NULL, 6},
{"cpugroup", 1, NULL, 7},
{"cpuoffset", 1, NULL, 8},
@@ -301,6 +274,8 @@ int main(int argc, char *argv[])
{"tolerance", 1, NULL, 13},
{"int-radius", 1, NULL, 14},
+ /* FIXME: Add '--no-peaks' and '--no-reflections' */
+
{"integrate-saturated",0, &integrate_saturated,1},
{"use-saturated", 0, &use_saturated, 1},
{"no-revalidate", 0, &no_revalidate, 1},
@@ -376,8 +351,7 @@ int main(int argc, char *argv[])
break;
case 5 :
- stream_flags = parse_stream_flags(optarg);
- if ( stream_flags < 0 ) return 1;
+ ERROR("The option '--record' is no longer used.\n");
break;
case 6 :
@@ -425,6 +399,12 @@ int main(int argc, char *argv[])
}
+ if ( scellr != NULL ) {
+ ERROR("Old-style indexing method specification.\n");
+ /* FIXME: Translate to new style */
+ return 1;
+ }
+
if ( filename == NULL ) {
filename = strdup("-");
}
@@ -439,17 +419,6 @@ 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");
}
@@ -482,53 +451,24 @@ 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",
@@ -586,7 +526,12 @@ int main(int argc, char *argv[])
cell = NULL;
}
- write_stream_header(ofh, argc, argv);
+ st = open_stream_for_write(outfile);
+ if ( st == NULL ) {
+ 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);
@@ -609,7 +554,7 @@ int main(int argc, char *argv[])
/* Prepare the indexer */
if ( indm != NULL ) {
ipriv = prepare_indexing(indm, cell, prepare_filename,
- det, beam);
+ det, beam, tols);
if ( ipriv == NULL ) {
ERROR("Failed to prepare indexing.\n");
return 1;
@@ -627,9 +572,7 @@ int main(int argc, char *argv[])
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];
@@ -644,7 +587,6 @@ int main(int argc, char *argv[])
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;
@@ -654,9 +596,11 @@ int main(int argc, char *argv[])
iargs.integrate_saturated = integrate_saturated;
iargs.no_revalidate = no_revalidate;
iargs.integrate_found = integrate_found;
+ iargs.include_peaks = 1;
+ iargs.include_reflections = 1; /* FIXME! */
create_sandbox(&iargs, n_proc, prefix, config_basename, fh,
- use_this_one_instead, ofh);
+ use_this_one_instead, st);
free(prefix);