aboutsummaryrefslogtreecommitdiff
path: root/src/indexamajig.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/indexamajig.c')
-rw-r--r--src/indexamajig.c261
1 files changed, 126 insertions, 135 deletions
diff --git a/src/indexamajig.c b/src/indexamajig.c
index bfeb2960..4b12194c 100644
--- a/src/indexamajig.c
+++ b/src/indexamajig.c
@@ -3,13 +3,13 @@
*
* Index patterns, output hkl+intensity etc.
*
- * Copyright © 2012-2020 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2021 Deutsches Elektronen-Synchrotron DESY,
* a research centre of the Helmholtz Association.
* Copyright © 2012 Richard Kirian
* Copyright © 2012 Lorenzo Galli
*
* Authors:
- * 2010-2019 Thomas White <taw@physics.org>
+ * 2010-2021 Thomas White <taw@physics.org>
* 2011 Richard Kirian
* 2012 Lorenzo Galli
* 2012 Chunhong Yoon
@@ -44,48 +44,26 @@
#include <string.h>
#include <unistd.h>
#include <argp.h>
-#include <hdf5.h>
#include <gsl/gsl_errno.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <fcntl.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 "geometry.h"
-#include "stream.h"
-#include "reflist-utils.h"
-#include "cell-utils.h"
-#include "integration.h"
-#include "taketwo.h"
-#include "im-sandbox.h"
-#include "image.h"
-
-
-static void add_geom_beam_stuff_to_field_list(struct imagefile_field_list *copyme,
- struct detector *det,
- struct beam_params *beam)
-{
- int i;
-
- for ( i=0; i<det->n_panels; i++ ) {
-
- struct panel *p = &det->panels[i];
+#include <utils.h>
+#include <index.h>
+#include <peaks.h>
+#include <filters.h>
+#include <thread-pool.h>
+#include <geometry.h>
+#include <stream.h>
+#include <reflist-utils.h>
+#include <cell-utils.h>
+#include <integration.h>
+#include <image.h>
+#include <datatemplate.h>
- if ( p->clen_from != NULL ) {
- add_imagefile_field(copyme, p->clen_from);
- }
- }
-
- if ( beam->photon_energy_from != NULL ) {
- add_imagefile_field(copyme, beam->photon_energy_from);
- }
-}
+#include "im-sandbox.h"
+#include "version.h"
struct indexamajig_arguments
@@ -99,20 +77,21 @@ struct indexamajig_arguments
int check_prefix;
int n_proc;
char *cellfile;
- char *spectrum_fn;
char *indm_str;
int basename;
int zmq;
int no_image_data;
+ int no_mask_data;
int serial_start;
char *temp_location;
- char *command_line_peak_path;
int if_refine;
int if_checkcell;
int if_peaks;
int if_multi;
int if_retry;
int profile; /* Whether to do wall-clock time profiling */
+ char **copy_headers;
+ int n_copy_headers;
TakeTwoOptions **taketwo_opts_ptr;
FelixOptions **felix_opts_ptr;
@@ -123,8 +102,23 @@ struct indexamajig_arguments
static void show_version(FILE *fh, struct argp_state *state)
{
- printf("CrystFEL: " CRYSTFEL_VERSIONSTRING "\n");
- printf(CRYSTFEL_BOILERPLATE"\n");
+ printf("CrystFEL: %s\n", crystfel_version_string());
+ printf("%s\n", crystfel_licence_string());
+}
+
+
+static void add_copy_header(struct indexamajig_arguments *args,
+ const char *header)
+{
+ char **new_copy_headers = realloc(args->copy_headers,
+ (args->n_copy_headers+1)*sizeof(char *));
+ if ( new_copy_headers == NULL ) {
+ ERROR("Failed to add copy header '%s'\n", header);
+ return;
+ }
+
+ args->copy_headers = new_copy_headers;
+ args->copy_headers[args->n_copy_headers++] = strdup(header);
}
@@ -212,9 +206,11 @@ static error_t parse_arg(int key, char *arg, struct argp_state *state)
break;
case 209 :
- args->spectrum_fn = strdup(arg);
- ERROR("WARNING: Prediction using arbitrary spectrum does not "
- "yet work in a useful way.\n");
+ ERROR("--spectrum-filename is no longer used.\n");
+ return 1;
+
+ case 210 :
+ args->no_mask_data = 1;
break;
/* ---------- Peak search ---------- */
@@ -224,21 +220,8 @@ static error_t parse_arg(int key, char *arg, struct argp_state *state)
break;
case 301 :
- if ( strcmp(arg, "zaef") == 0 ) {
- args->iargs.peaks = PEAK_ZAEF;
- } else if ( strcmp(arg, "peakfinder8") == 0 ) {
- args->iargs.peaks = PEAK_PEAKFINDER8;
- } else if ( strcmp(arg, "hdf5") == 0 ) {
- args->iargs.peaks = PEAK_HDF5;
- } else if ( strcmp(arg, "cxi") == 0 ) {
- args->iargs.peaks = PEAK_CXI;
- } else if ( strcmp(arg, "peakfinder9") == 0 ) {
- args->iargs.peaks = PEAK_PEAKFINDER9;
- } else if ( strcmp(arg, "msgpack") == 0 ) {
- args->iargs.peaks = PEAK_MSGPACK;
- } else if ( strcmp(arg, "none") == 0 ) {
- args->iargs.peaks = PEAK_NONE;
- } else {
+ args->iargs.peaks = parse_peaksearch(arg);
+ if ( args->iargs.peaks == PEAK_ERROR ) {
ERROR("Unrecognised peak detection method '%s'\n", arg);
return EINVAL;
}
@@ -259,12 +242,11 @@ static error_t parse_arg(int key, char *arg, struct argp_state *state)
ERROR("Invalid value for --min-peaks\n");
return EINVAL;
}
- (*(args->pinkindexer_opts_ptr))->min_peaks = args->iargs.min_peaks;
break;
case 304 :
- free(args->command_line_peak_path);
- args->command_line_peak_path = strdup(arg);
+ ERROR("The option --hdf5-peak-path is no longer used.\n");
+ ERROR("Set the peak path in the geometry file.\n");
break;
case 305 :
@@ -395,6 +377,8 @@ static error_t parse_arg(int key, char *arg, struct argp_state *state)
break;
case 401 :
+ /* Values in 'tols' are in frac (not %) and rad
+ * Conversion happens a few lines below */
r = sscanf(arg, "%f,%f,%f,%f,%f,%f",
&args->iargs.tols[0], &args->iargs.tols[1], &args->iargs.tols[2],
&args->iargs.tols[3], &args->iargs.tols[4], &args->iargs.tols[5]);
@@ -466,6 +450,21 @@ static error_t parse_arg(int key, char *arg, struct argp_state *state)
* better than the user expected. */
break;
+ case 413 :
+ if (sscanf(arg, "%f", &args->iargs.wavelength_estimate) != 1)
+ {
+ ERROR("Invalid value for --wavelength-estimate\n");
+ return EINVAL;
+ }
+ break;
+
+ case 414 :
+ if (sscanf(arg, "%d", &args->iargs.n_threads) != 1)
+ {
+ ERROR("Invalid value for --max-indexer-threads\n");
+ return EINVAL;
+ }
+ break;
/* ---------- Integration ---------- */
case 501 :
@@ -559,15 +558,17 @@ static error_t parse_arg(int key, char *arg, struct argp_state *state)
break;
case 602 :
- add_imagefile_field(args->iargs.copyme, arg);
+ add_copy_header(args, arg);
break;
case 603 :
- args->iargs.stream_peaks = 0;
+ args->iargs.stream_flags = CLEAR_BIT(args->iargs.stream_flags,
+ STREAM_PEAKS);
break;
case 604 :
- args->iargs.stream_refls = 0;
+ args->iargs.stream_flags = CLEAR_BIT(args->iargs.stream_flags,
+ STREAM_REFLECTIONS);
break;
case 605 :
@@ -594,13 +595,13 @@ int main(int argc, char *argv[])
char *tmpdir; /* e.g. /tmp/indexamajig.12345 */
char *rn; /* e.g. /home/taw/indexing */
int r;
- struct beam_params beam;
char *zmq_address = NULL;
int timeout = 240;
TakeTwoOptions *taketwo_opts = NULL;
FelixOptions *felix_opts = NULL;
XGandalfOptions *xgandalf_opts = NULL;
PinkIndexerOptions *pinkindexer_opts = NULL;
+ double wl_from_dt;
/* Defaults for "top level" arguments */
args.filename = NULL;
@@ -611,18 +612,18 @@ int main(int argc, char *argv[])
args.check_prefix = 1;
args.n_proc = 1;
args.cellfile = NULL;
- args.spectrum_fn = NULL;
args.indm_str = NULL;
args.basename = 0;
args.zmq = 0;
args.serial_start = 1;
- args.command_line_peak_path = NULL;
args.if_peaks = 1;
args.if_multi = 0;
args.if_retry = 1;
args.if_refine = 1;
args.if_checkcell = 1;
args.profile = 0;
+ args.copy_headers = NULL;
+ args.n_copy_headers = 0;
args.taketwo_opts_ptr = &taketwo_opts;
args.felix_opts_ptr = &felix_opts;
args.xgandalf_opts_ptr = &xgandalf_opts;
@@ -632,12 +633,12 @@ int main(int argc, char *argv[])
args.iargs.cell = NULL;
args.iargs.noisefilter = 0;
args.iargs.median_filter = 0;
- args.iargs.tols[0] = 0.05;
- args.iargs.tols[1] = 0.05;
- args.iargs.tols[2] = 0.05;
- args.iargs.tols[3] = deg2rad(1.5);
- args.iargs.tols[4] = deg2rad(1.5);
- args.iargs.tols[5] = deg2rad(1.5);
+ args.iargs.tols[0] = 0.05; /* frac (not %) */
+ args.iargs.tols[1] = 0.05; /* frac (not %) */
+ args.iargs.tols[2] = 0.05; /* frac (not %) */
+ args.iargs.tols[3] = deg2rad(1.5); /* radians */
+ args.iargs.tols[4] = deg2rad(1.5); /* radians */
+ args.iargs.tols[5] = deg2rad(1.5); /* radians */
args.iargs.threshold = 800.0;
args.iargs.min_sq_gradient = 100000.0;
args.iargs.min_snr = 5.0;
@@ -651,12 +652,9 @@ int main(int argc, char *argv[])
args.iargs.min_sig = 11.0;
args.iargs.min_peak_over_neighbour = -INFINITY;
args.iargs.check_hdf5_snr = 0;
- args.iargs.det = NULL;
+ args.iargs.dtempl = NULL;
args.iargs.peaks = PEAK_ZAEF;
- args.iargs.beam = &beam;
- args.iargs.hdf5_peak_path = NULL;
args.iargs.half_pixel_shift = 1;
- args.iargs.copyme = NULL;
args.iargs.pk_inn = -1.0;
args.iargs.pk_mid = -1.0;
args.iargs.pk_out = -1.0;
@@ -665,18 +663,12 @@ int main(int argc, char *argv[])
args.iargs.ir_out = -1.0;
args.iargs.use_saturated = 1;
args.iargs.no_revalidate = 0;
- args.iargs.stream_peaks = 1;
- args.iargs.stream_refls = 1;
+ args.iargs.stream_flags = STREAM_PEAKS | STREAM_REFLECTIONS;
args.iargs.stream_nonhits = 1;
args.iargs.int_diag = INTDIAG_NONE;
args.iargs.min_peaks = 0;
args.iargs.overpredict = 0;
args.iargs.wait_for_file = 0;
- args.iargs.copyme = new_imagefile_field_list();
- if ( args.iargs.copyme == NULL ) {
- ERROR("Couldn't allocate HDF5 field list.\n");
- return 1;
- }
args.iargs.ipriv = NULL; /* No default */
args.iargs.int_meth = integration_method("rings-nocen-nosat-nograd", NULL);
args.iargs.push_res = +INFINITY;
@@ -684,6 +676,9 @@ int main(int argc, char *argv[])
args.iargs.fix_profile_r = -1.0;
args.iargs.fix_divergence = -1.0;
args.iargs.no_image_data = 0;
+ args.iargs.no_mask_data = 0;
+ args.iargs.wavelength_estimate = NAN;
+ args.iargs.n_threads = 1;
argp_program_version_hook = show_version;
@@ -716,15 +711,16 @@ int main(int argc, char *argv[])
"processing"},
{"zmq-msgpack", 207, NULL, OPTION_NO_USAGE, "Receive data in MessagePack format "
"over ZMQ"},
- {"no-image-data", 208, NULL, OPTION_NO_USAGE, "Do not load image data (from ZMQ)"},
+ {"no-image-data", 208, NULL, OPTION_NO_USAGE, "Do not load image data"},
{"spectrum-file", 209, "fn", OPTION_NO_USAGE | OPTION_HIDDEN,
"File containing radiation spectrum"},
+ {"no-mask-data", 210, NULL, OPTION_NO_USAGE, "Do not load mask data"},
{NULL, 0, 0, OPTION_DOC, "Peak search options:", 3},
{"peaks", 301, "method", 0, "Peak search method. Default: zaef"},
{"peak-radius", 302, "r1,r2,r3", OPTION_NO_USAGE, "Radii for peak search"},
{"min-peaks", 303, "n", OPTION_NO_USAGE, "Minimum number of peaks for indexing"},
- {"hdf5-peaks", 304, "p", OPTION_NO_USAGE, "Location of peak table in HDF5 file"},
+ {"hdf5-peaks", 304, "p", OPTION_HIDDEN, "Location of peak table in HDF5 file"},
{"median-filter", 305, "n", OPTION_NO_USAGE, "Apply median filter to image data"},
{"filter-noise", 306, NULL, OPTION_NO_USAGE, "Apply noise filter to image data"},
{"threshold", 't', "adu", OPTION_NO_USAGE, "Threshold for peak detection "
@@ -780,6 +776,10 @@ int main(int argc, char *argv[])
"accounted for by the indexing solution"},
{"check-peaks", 411, NULL, OPTION_HIDDEN, NULL},
{"no-cell-combinations", 412, NULL, OPTION_HIDDEN, NULL},
+ {"wavelength-estimate", 413, "metres", 0,
+ "Estimate of the incident radiation wavelength, in metres."},
+ {"max-indexer-threads", 414, "n", 0,
+ "Maximum number of threads allowed for indexing engines."},
{NULL, 0, 0, OPTION_DOC, "Integration options:", 5},
{"integration", 501, "method", OPTION_NO_USAGE, "Integration method"},
@@ -797,8 +797,9 @@ int main(int argc, char *argv[])
{NULL, 0, 0, OPTION_DOC, "Output options:", 6},
{"no-non-hits-in-stream", 601, NULL, OPTION_NO_USAGE, "Don't include non-hits in "
"stream (see --min-peaks)"},
- {"copy-hdf5-field", 602, "f", OPTION_NO_USAGE, "Put the value of this HDF5 field "
- "into the stream"},
+ {"copy-hdf5-field", 602, "f", OPTION_HIDDEN, NULL},
+ {"copy-header", 602, "f", OPTION_NO_USAGE, "Put the value of this image header "
+ "field into the stream"},
{"no-peaks-in-stream", 603, NULL, OPTION_NO_USAGE, "Don't put peak search results "
"in stream"},
{"no-refls-in-stream", 604, NULL, OPTION_NO_USAGE, "Don't put integration results "
@@ -861,31 +862,18 @@ int main(int argc, char *argv[])
return 1;
}
- /* Load detector geometry */
- args.iargs.det = get_detector_geometry_2(args.geom_filename,
- args.iargs.beam,
- &args.iargs.hdf5_peak_path);
- if ( args.iargs.det == NULL ) {
- ERROR("Failed to read detector geometry from '%s'\n",
- args.geom_filename);
+ /* Load data template (new API) */
+ args.iargs.dtempl = data_template_new_from_file(args.geom_filename);
+ if ( args.iargs.dtempl == NULL ) {
+ ERROR("Failed to read detector geometry from '%s'"
+ " (for new API)\n", args.geom_filename);
return 1;
}
- add_geom_beam_stuff_to_field_list(args.iargs.copyme, args.iargs.det,
- args.iargs.beam);
-
- /* If no peak path from geometry file, use these (but see later) */
- if ( args.iargs.hdf5_peak_path == NULL ) {
- if ( args.iargs.peaks == PEAK_HDF5 ) {
- args.iargs.hdf5_peak_path = strdup("/processing/hitfinder/peakinfo");
- } else if ( args.iargs.peaks == PEAK_CXI ) {
- args.iargs.hdf5_peak_path = strdup("/entry_1/result_1");
- }
- }
- /* If an HDF5 peak path was given on the command line, use it */
- if ( args.command_line_peak_path != NULL ) {
- free(args.iargs.hdf5_peak_path);
- args.iargs.hdf5_peak_path = args.command_line_peak_path;
+ /* Add any headers we need to copy */
+ for ( r=0; r<args.n_copy_headers; r++ ) {
+ data_template_add_copy_header(args.iargs.dtempl,
+ args.copy_headers[r]);
}
/* If no integration radii were given, apply the defaults */
@@ -917,18 +905,6 @@ int main(int argc, char *argv[])
args.iargs.cell = NULL;
}
- /* Load spectrum from file if given */
- if ( args.spectrum_fn != NULL ) {
- args.iargs.spectrum = spectrum_load(args.spectrum_fn);
- if ( args.iargs.spectrum == NULL ) {
- ERROR("Couldn't read spectrum (from %s)\n", args.spectrum_fn);
- return 1;
- }
- free(args.spectrum_fn);
- } else {
- args.iargs.spectrum = NULL;
- }
-
tmpdir = create_tempdir(args.temp_location);
if ( tmpdir == NULL ) return 1;
@@ -955,6 +931,16 @@ int main(int argc, char *argv[])
}
+ wl_from_dt = data_template_get_wavelength_if_possible(args.iargs.dtempl);
+ if ( !isnan(wl_from_dt) ) {
+ if ( !isnan(args.iargs.wavelength_estimate) ) {
+ ERROR("WARNING: Ignoring your value for "
+ "--wavelength-estimate because the geometry file "
+ "already contains a static value.\n");
+ }
+ args.iargs.wavelength_estimate = wl_from_dt;
+ }
+
/* Prepare the indexing system */
if ( args.indm_str == NULL ) {
@@ -1001,9 +987,12 @@ int main(int argc, char *argv[])
flags |= INDEXING_RETRY;
}
- args.iargs.ipriv = setup_indexing(args.indm_str, args.iargs.cell,
- args.iargs.det, args.iargs.beam,
- args.iargs.tols, flags,
+ args.iargs.ipriv = setup_indexing(args.indm_str,
+ args.iargs.cell,
+ args.iargs.tols,
+ flags,
+ args.iargs.wavelength_estimate,
+ args.iargs.n_threads,
taketwo_opts,
xgandalf_opts,
pinkindexer_opts,
@@ -1039,13 +1028,18 @@ int main(int argc, char *argv[])
free(rn);
/* Open output stream */
- st = open_stream_for_write_4(args.outfile, args.geom_filename,
- args.iargs.cell, argc, argv,
- args.indm_str);
+ st = stream_open_for_write(args.outfile, args.iargs.dtempl);
if ( st == NULL ) {
ERROR("Failed to open stream '%s'\n", args.outfile);
return 1;
}
+
+ /* Write audit info */
+ stream_write_commandline_args(st, argc, argv);
+ stream_write_geometry_file(st, args.geom_filename);
+ stream_write_target_cell(st, args.iargs.cell);
+ stream_write_indexing_methods(st, args.indm_str);
+
free(args.outfile);
free(args.indm_str);
@@ -1069,15 +1063,12 @@ int main(int argc, char *argv[])
fh, st, tmpdir, args.serial_start, zmq_address,
timeout, args.profile);
- free_imagefile_field_list(args.iargs.copyme);
cell_free(args.iargs.cell);
- free(args.iargs.beam->photon_energy_from);
free(args.prefix);
free(args.temp_location);
free(tmpdir);
- free_detector_geometry(args.iargs.det);
- free(args.iargs.hdf5_peak_path);
- close_stream(st);
+ data_template_free(args.iargs.dtempl);
+ stream_close(st);
cleanup_indexing(args.iargs.ipriv);
return r;
generated by cgit v1.2.3 (git 2.25.1) at 2024-08-18 14:16:04 +0000