indexamajig.c 27.8 KB
Newer Older
Thomas White's avatar
Thomas White committed
1
/*
2
 * indexamajig.c
Thomas White's avatar
Thomas White committed
3
 *
Thomas White's avatar
Thomas White committed
4
 * Index patterns, output hkl+intensity etc.
Thomas White's avatar
Thomas White committed
5
 *
Thomas White's avatar
Thomas White committed
6
 * Copyright © 2012 Thomas White <taw@physics.org>
Thomas White's avatar
Thomas White committed
7
 *
Thomas White's avatar
Thomas White committed
8
9
10
11
12
13
14
15
16
17
18
19
20
21
 * 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/>.
Thomas White's avatar
Thomas White committed
22
23
24
25
26
27
28
29
30
31
32
33
34
35
 *
 */


#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <stdarg.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <unistd.h>
#include <getopt.h>
Thomas White's avatar
Thomas White committed
36
#include <hdf5.h>
Thomas White's avatar
Thomas White committed
37
#include <gsl/gsl_errno.h>
38
#include <pthread.h>
39
40

#ifdef HAVE_CLOCK_GETTIME
41
#include <time.h>
42
43
44
#else
#include <sys/time.h>
#endif
Thomas White's avatar
Thomas White committed
45
46
47

#include "utils.h"
#include "hdf5-file.h"
Thomas White's avatar
Thomas White committed
48
#include "index.h"
49
#include "peaks.h"
50
#include "detector.h"
Thomas White's avatar
Thomas White committed
51
#include "filters.h"
52
#include "thread-pool.h"
53
#include "beam-parameters.h"
54
#include "geometry.h"
Thomas White's avatar
Thomas White committed
55
#include "stream.h"
Thomas White's avatar
Thomas White committed
56
#include "reflist-utils.h"
57

58

Thomas White's avatar
Thomas White committed
59
/* Write statistics at APPROXIMATELY this interval */
60
61
62
#define STATS_EVERY_N_SECONDS (5)


63
64
65
66
67
68
enum {
	PEAK_ZAEF,
	PEAK_HDF5,
};


69
70
/* Information about the indexing process which is common to all patterns */
struct static_index_args
71
72
73
74
75
{
	UnitCell *cell;
	int config_cmfilter;
	int config_noisefilter;
	int config_verbose;
Thomas White's avatar
Thomas White committed
76
	int stream_flags;         /* What goes into the output? */
77
	int config_satcorr;
78
	int config_closer;
79
	int config_insane;
80
	int config_bgsub;
81
	float threshold;
82
	float min_gradient;
83
	float min_snr;
84
	double min_int_snr;
Thomas White's avatar
Thomas White committed
85
	struct detector *det;
Thomas White's avatar
Thomas White committed
86
87
	IndexingMethod *indm;
	IndexingPrivate **ipriv;
Thomas White's avatar
Thomas White committed
88
	int peaks;                /* Peak detection method */
89
	int cellr;
90
	float tols[4];
91
	struct beam_params *beam;
Thomas White's avatar
Thomas White committed
92
	const char *element;
93
	const char *hdf5_peak_path;
94

95
96
97
	/* Output stream */
	pthread_mutex_t *output_mutex;  /* Protects the output stream */
	FILE *ofh;
98
	const struct copy_hdf5_field *copyme;
99
100
101
};


102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
/* Information about the indexing process for one pattern */
struct index_args
{
	/* "Input" */
	char *filename;
	struct static_index_args static_args;

	/* "Output" */
	int indexable;
};


/* Information needed to choose the next task and dispatch it */
struct queue_args
{
	FILE *fh;
	char *prefix;
119
	int config_basename;
120
121
	struct static_index_args static_args;

122
	char *use_this_one_instead;
123
124
125
126
127
128

	int n_indexable;
	int n_processed;
	int n_indexable_last_stats;
	int n_processed_last_stats;
	int t_last_stats;
129
130
131
};


Thomas White's avatar
Thomas White committed
132
133
134
135
136
137
static void show_help(const char *s)
{
	printf("Syntax: %s [options]\n\n", s);
	printf(
"Process and index FEL diffraction images.\n"
"\n"
138
" -h, --help               Display this help message.\n"
Thomas White's avatar
Thomas White committed
139
"\n"
140
" -i, --input=<filename>   Specify file containing list of images to process.\n"
Thomas White's avatar
Thomas White committed
141
"                           '-' means stdin, which is the default.\n"
Thomas White's avatar
Thomas White committed
142
143
" -o, --output=<filename>  Write output stream to this file. '-' for stdout.\n"
"                           Default: indexamajig.stream\n"
Thomas White's avatar
Thomas White committed
144
"\n"
Thomas White's avatar
Thomas White committed
145
146
147
148
149
"     --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"
Thomas White's avatar
Thomas White committed
150
"                            reax     : DPS algorithm with known unit cell\n"
151
" -g. --geometry=<file>    Get detector geometry from file.\n"
152
153
154
" -b, --beam=<file>        Get beam parameters from file (provides nominal\n"
"                           wavelength value if no per-shot value is found in\n"
"                           the HDF5 files.\n"
155
" -p, --pdb=<file>         PDB file from which to get the unit cell to match.\n"
Thomas White's avatar
Thomas White committed
156
"                           Default: 'molecule.pdb'.\n"
157
"     --basename           Remove the directory parts of the filenames.\n"
158
" -x, --prefix=<p>         Prefix filenames from input file with <p>.\n"
159
160
161
"     --peaks=<method>     Use 'method' for finding peaks.  Choose from:\n"
"                           zaef  : Use Zaefferer (2000) gradient detection.\n"
"                                    This is the default method.\n"
162
163
164
"                           hdf5  : Get from a table in HDF5 file.\n"
"     --hdf5-peaks=<p>     Find peaks table in HDF5 file here.\n"
"                           Default: /processing/hitfinder/peakinfo\n"
Thomas White's avatar
Thomas White committed
165
166
"\n\n"
"You can control what information is included in the output stream using\n"
Thomas White's avatar
Thomas White committed
167
"' --record=<flag1>,<flag2>,<flag3>' and so on.  Possible flags are:\n\n"
168
169
" integrated        Include a list of reflection intensities, produced by\n"
"                    integrating around predicted peak locations.\n"
Thomas White's avatar
Thomas White committed
170
"\n"
171
172
" peaks             Include peak locations and intensities from the peak\n"
"                    search.\n"
Thomas White's avatar
Thomas White committed
173
"\n"
174
" peaksifindexed    As 'peaks', but only if the pattern could be indexed.\n"
Thomas White's avatar
Thomas White committed
175
"\n"
176
177
" peaksifnotindexed As 'peaks', but only if the pattern could NOT be indexed.\n"
"\n\n"
Thomas White's avatar
Thomas White committed
178
"The default is '--record=integrated'.\n"
Thomas White's avatar
Thomas White committed
179
180
"\n\n"
"For more control over the process, you might need:\n\n"
Thomas White's avatar
Thomas White committed
181
182
183
184
185
"     --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"
186
"                           compare_ab : compare 'a' and 'b' lengths only.\n"
187
188
189
"     --tolerance=<a,b,c,angl>  Set the tolerance for a,b,c axis (in %%)\n"
"                                and for the angles (in deg) when reducing\n"
"                                or comparing (default is 5%% and 1.5deg)\n"
190
"     --filter-cm          Perform common-mode noise subtraction on images\n"
Thomas White's avatar
Thomas White committed
191
192
"                           before proceeding.  Intensities will be extracted\n"
"                           from the image as it is after this processing.\n"
193
"     --filter-noise       Apply an aggressive noise filter which sets all\n"
Thomas White's avatar
Thomas White committed
194
195
196
"                           pixels in each 3x3 region to zero if any of them\n"
"                           have negative values.  Intensity measurement will\n"
"                           be performed on the image as it was before this.\n"
197
198
"     --no-sat-corr        Don't correct values of saturated peaks using a\n"
"                           table included in the HDF5 file.\n"
199
"     --threshold=<n>      Only accept peaks above <n> ADU.  Default: 800.\n"
200
201
"     --min-gradient=<n>   Minimum gradient for Zaefferer peak search.\n"
"                           Default: 100,000.\n"
202
203
"     --min-snr=<n>        Minimum signal-to-noise ratio for peaks.\n"
"                           Default: 5.\n"
204
205
"     --min-integration-snr=<n>  Minimum signal-to-noise ratio for peaks\n"
"                                 during integration. Default: -infinity.\n"
Thomas White's avatar
Thomas White committed
206
" -e, --image=<element>    Use this image from the HDF5 file.\n"
Thomas White's avatar
Thomas White committed
207
"                           Example: /data/data0.\n"
Thomas White's avatar
Thomas White committed
208
"                           Default: The first one found.\n"
Thomas White's avatar
Thomas White committed
209
"\n"
210
211
212
213
"\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"
214
215
216
"\nOptions for greater performance or verbosity:\n\n"
"     --verbose            Be verbose about indexing.\n"
" -j <n>                   Run <n> analyses in parallel.  Default 1.\n"
217
218
219
"\n"
"\nOptions you probably won't need:\n\n"
"     --no-check-prefix    Don't attempt to correct the --prefix.\n"
220
221
222
"     --no-closer-peak     Don't integrate from the location of a nearby peak\n"
"                           instead of the position closest to the reciprocal\n"
"                           lattice point.\n"
223
224
"     --insane             Don't check that the reduced cell accounts for at\n"
"                           least 10%% of the located peaks.\n"
225
226
"     --no-bg-sub          Don't subtract local background estimates from\n"
"                           integrated intensities.\n"
227
"\n"
Thomas White's avatar
Thomas White committed
228
229
"\nYou can tune the CPU affinities for enhanced performance on NUMA machines:\n"
"\n"
Thomas White's avatar
Thomas White committed
230
231
"     --cpus=<n>           Specify number of CPUs.  This is NOT the same as\n"
"                           giving the number of analyses to run in parallel.\n"
232
233
"     --cpugroup=<n>       Batch threads in groups of this size.\n"
"     --cpuoffset=<n>      Start using CPUs at this group number.\n"
Thomas White's avatar
Thomas White committed
234
);
Thomas White's avatar
Thomas White committed
235
236
237
}


238
static void process_image(void *pp, int cookie)
239
{
240
	struct index_args *pargs = pp;
241
242
243
244
	struct hdfile *hdfile;
	struct image image;
	float *data_for_measurement;
	size_t data_size;
245
	char *filename = pargs->filename;
246
247
248
249
	UnitCell *cell = pargs->static_args.cell;
	int config_cmfilter = pargs->static_args.config_cmfilter;
	int config_noisefilter = pargs->static_args.config_noisefilter;
	int config_verbose = pargs->static_args.config_verbose;
Thomas White's avatar
Thomas White committed
250
	IndexingMethod *indm = pargs->static_args.indm;
251
	struct beam_params *beam = pargs->static_args.beam;
252
253
254

	image.features = NULL;
	image.data = NULL;
Thomas White's avatar
Thomas White committed
255
	image.flags = NULL;
256
	image.indexed_cell = NULL;
257
	image.id = cookie;
258
	image.filename = filename;
259
	image.det = copy_geom(pargs->static_args.det);
260
	image.copyme = pargs->static_args.copyme;
261
262
263
264
265
266
267
	image.beam = beam;

	if ( beam == NULL ) {
		ERROR("Warning: no beam parameters file.\n");
		ERROR("I'm going to assume 1 ADU per photon, which is almost");
		ERROR(" certainly wrong.  Peak sigmas will be incorrect.\n");
	}
268

Thomas White's avatar
Thomas White committed
269
	pargs->indexable = 0;
270

271
	hdfile = hdfile_open(filename);
Thomas White's avatar
Thomas White committed
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
	if ( hdfile == NULL ) return;

	if ( pargs->static_args.element != NULL ) {

		int r;
		r = hdfile_set_image(hdfile, pargs->static_args.element);
		if ( r ) {
			ERROR("Couldn't select path '%s'\n",
			      pargs->static_args.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;
		}

295
296
	}

297
	hdf5_read(hdfile, &image, pargs->static_args.config_satcorr);
298
299
300
301
302
303
304
305
306
307
308
309
310
311

	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);
		free_detector_geometry(image.det);
		return;
	}

312
313
	if ( image.lambda < 0.0 ) {
		if ( beam != NULL ) {
Thomas White's avatar
Thomas White committed
314
			ERROR("Using nominal photon energy of %.2f eV\n",
315
                              beam->photon_energy);
Thomas White's avatar
Thomas White committed
316
317
			image.lambda = ph_en_to_lambda(
			                          eV_to_J(beam->photon_energy));
318
319
320
321
		} else {
			ERROR("No wavelength in file, so you need to give "
			      "a beam parameters file with -b.\n");
			hdfile_close(hdfile);
Thomas White's avatar
Thomas White committed
322
			free_detector_geometry(image.det);
323
324
325
			return;
		}
	}
326
	fill_in_values(image.det, hdfile);
327
328
329
330
331
332
333
334
335
336
337
338
339

	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 {
340
		memcpy(data_for_measurement, image.data, data_size);
341
342
	}

343
	switch ( pargs->static_args.peaks )
344
345
346
	{
	case PEAK_HDF5 :
		/* Get peaks from HDF5 */
347
348
349
		if ( get_peaks(&image, hdfile,
		               pargs->static_args.hdf5_peak_path) )
		{
350
351
352
353
			ERROR("Failed to get peaks from HDF5 file.\n");
		}
		break;
	case PEAK_ZAEF :
354
		search_peaks(&image, pargs->static_args.threshold,
355
356
		             pargs->static_args.min_gradient,
		             pargs->static_args.min_snr);
357
358
		break;
	}
Thomas White's avatar
Thomas White committed
359
360
361
362
363
364

	/* 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;

365
	/* Calculate orientation matrix (by magic) */
366
367
368
	image.div = beam->divergence;
	image.bw = beam->bandwidth;
	image.profile_radius = 0.0001e9;
369
	index_pattern(&image, cell, indm, pargs->static_args.cellr,
370
		      config_verbose, pargs->static_args.ipriv,
371
		      pargs->static_args.config_insane, pargs->static_args.tols);
372

373
	if ( image.indexed_cell != NULL ) {
Thomas White's avatar
Thomas White committed
374

375
		pargs->indexable = 1;
376

377
378
		image.reflections = find_intersections(&image,
		                                       image.indexed_cell);
379

380
		if ( image.reflections != NULL ) {
381

382
			integrate_reflections(&image,
383
					      pargs->static_args.config_closer,
384
385
					      pargs->static_args.config_bgsub,
					      pargs->static_args.min_int_snr);
386

387
		}
388

389
390
391
392
393
	} else {

		image.reflections = NULL;

	}
Thomas White's avatar
Thomas White committed
394
395

	pthread_mutex_lock(pargs->static_args.output_mutex);
396
	write_chunk(pargs->static_args.ofh, &image, hdfile,
Thomas White's avatar
Thomas White committed
397
398
	            pargs->static_args.stream_flags);
	pthread_mutex_unlock(pargs->static_args.output_mutex);
399
400

	/* Only free cell if found */
401
	cell_free(image.indexed_cell);
402

Thomas White's avatar
Thomas White committed
403
	reflist_free(image.reflections);
404
	free(image.data);
Thomas White's avatar
Thomas White committed
405
	if ( image.flags != NULL ) free(image.flags);
406
407
	image_feature_list_free(image.features);
	hdfile_close(hdfile);
Thomas White's avatar
Thomas White committed
408
	free_detector_geometry(image.det);
409
410
411
}


412
static void *get_image(void *qp)
413
{
414
	char *line;
415
416
417
	struct index_args *pargs;
	char *rval;
	struct queue_args *qargs = qp;
418

419
	/* Initialise new task arguments */
420
421
422
	pargs = malloc(sizeof(struct index_args));
	memcpy(&pargs->static_args, &qargs->static_args,
	       sizeof(struct static_index_args));
423

424
425
426
	/* Get the next filename */
	if ( qargs->use_this_one_instead != NULL ) {

427
		line = qargs->use_this_one_instead;
428
429
430
431
		qargs->use_this_one_instead = NULL;

	} else {

432
		line = malloc(1024*sizeof(char));
433
		rval = fgets(line, 1023, qargs->fh);
Thomas White's avatar
Thomas White committed
434
435
		if ( rval == NULL ) {
			free(pargs);
Thomas White's avatar
Thomas White committed
436
			free(line);
Thomas White's avatar
Thomas White committed
437
438
			return NULL;
		}
439
440
441
		chomp(line);

	}
442

443
444
445
446
447
448
449
450
451
452
453
454
455
	if ( qargs->config_basename ) {
		char *tmp;
		tmp = safe_basename(line);
		free(line);
		line = tmp;
	}

	pargs->filename = malloc(strlen(qargs->prefix)+strlen(line)+1);

	snprintf(pargs->filename, 1023, "%s%s", qargs->prefix, line);

	free(line);

456
457
	return pargs;
}
458
459


460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
#ifdef HAVE_CLOCK_GETTIME

static time_t get_monotonic_seconds()
{
	struct timespec tp;
	clock_gettime(CLOCK_MONOTONIC, &tp);
	return tp.tv_sec;
}

#else

/* Fallback version of the above.  The time according to gettimeofday() is not
 * monotonic, so measuring intervals based on it will screw up if there's a
 * timezone change (e.g. daylight savings) while the program is running. */
static time_t get_monotonic_seconds()
{
	struct timeval tp;
	gettimeofday(&tp, NULL);
	return tp.tv_sec;
}

#endif

483
484
485
486
static void finalise_image(void *qp, void *pp)
{
	struct queue_args *qargs = qp;
	struct index_args *pargs = pp;
487
	time_t monotonic_seconds;
488

489
	qargs->n_indexable += pargs->indexable;
490
491
	qargs->n_processed++;

492
493
	monotonic_seconds = get_monotonic_seconds();
	if ( monotonic_seconds >= qargs->t_last_stats+STATS_EVERY_N_SECONDS ) {
494
495

		STATUS("%i out of %i indexed so far,"
Thomas White's avatar
Thomas White committed
496
		       " %i out of %i since the last message.\n",
497
498
		       qargs->n_indexable, qargs->n_processed,
		       qargs->n_indexable - qargs->n_indexable_last_stats,
Thomas White's avatar
Thomas White committed
499
		       qargs->n_processed - qargs->n_processed_last_stats);
500
501
502

		qargs->n_processed_last_stats = qargs->n_processed;
		qargs->n_indexable_last_stats = qargs->n_indexable;
503
		qargs->t_last_stats = monotonic_seconds;
504
505

	}
506

507
508
	free(pargs->filename);
	free(pargs);
509
510
511
}


512
513
514
static int parse_cell_reduction(const char *scellr, int *err,
                                int *reduction_needs_cell)
{
Thomas White's avatar
Thomas White committed
515
	*err = 0;
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
	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;
	}
}


Thomas White's avatar
Thomas White committed
536
537
538
539
int main(int argc, char *argv[])
{
	int c;
	char *filename = NULL;
540
	char *outfile = NULL;
Thomas White's avatar
Thomas White committed
541
	FILE *fh;
542
	FILE *ofh;
543
	char *rval = NULL;
Thomas White's avatar
Thomas White committed
544
	int n_images;
545
	int config_noindex = 0;
546
547
	int config_cmfilter = 0;
	int config_noisefilter = 0;
548
	int config_verbose = 0;
549
	int config_satcorr = 1;
550
	int config_checkprefix = 1;
551
	int config_closer = 1;
552
	int config_insane = 0;
553
	int config_bgsub = 1;
554
	int config_basename = 0;
555
	float threshold = 800.0;
556
	float min_gradient = 100000.0;
Thomas White's avatar
Thomas White committed
557
	float min_snr = 5.0;
558
	double min_int_snr = -INFINITY;
Thomas White's avatar
Thomas White committed
559
560
	struct detector *det;
	char *geometry = NULL;
Thomas White's avatar
Thomas White committed
561
562
563
564
	IndexingMethod *indm;
	IndexingPrivate **ipriv;
	int indexer_needs_cell;
	int reduction_needs_cell;
Thomas White's avatar
Thomas White committed
565
	char *indm_str = NULL;
566
	UnitCell *cell;
Thomas White's avatar
Thomas White committed
567
	char *pdb = NULL;
568
	char *prefix = NULL;
569
	char *speaks = NULL;
570
	char *scellr = NULL;
571
	char *toler = NULL;
Thomas White's avatar
Thomas White committed
572
	float tols[4] = {5.0, 5.0, 5.0, 1.5}; /* a,b,c,angles (%,%,%,deg) */
573
	int cellr;
574
	int peaks;
575
	int nthreads = 1;
576
	pthread_mutex_t output_mutex = PTHREAD_MUTEX_INITIALIZER;
577
	char *prepare_line;
578
	char prepare_filename[1024];
579
	struct queue_args qargs;
580
	struct beam_params *beam = NULL;
Thomas White's avatar
Thomas White committed
581
	char *element = NULL;
582
	double nominal_photon_energy;
Thomas White's avatar
Thomas White committed
583
	int stream_flags = STREAM_INTEGRATED;
584
585
586
587
	int cpu_num = 0;
	int cpu_groupsize = 1;
	int cpu_offset = 0;
	char *endptr;
588
	char *hdf5_peak_path = NULL;
589
590
591
592
593
594
595
	struct copy_hdf5_field *copyme;

	copyme = new_copy_hdf5_field_list();
	if ( copyme == NULL ) {
		ERROR("Couldn't allocate HDF5 field list.\n");
		return 1;
	}
Thomas White's avatar
Thomas White committed
596
597
598
599
600

	/* Long options */
	const struct option longopts[] = {
		{"help",               0, NULL,               'h'},
		{"input",              1, NULL,               'i'},
601
		{"output",             1, NULL,               'o'},
602
		{"no-index",           0, &config_noindex,     1},
603
		{"peaks",              1, NULL,                2},
604
		{"cell-reduction",     1, NULL,                3},
605
		{"tolerance",          1, NULL,               13},
Thomas White's avatar
Thomas White committed
606
		{"indexing",           1, NULL,               'z'},
Thomas White's avatar
Thomas White committed
607
		{"geometry",           1, NULL,               'g'},
608
		{"beam",               1, NULL,               'b'},
609
610
		{"filter-cm",          0, &config_cmfilter,    1},
		{"filter-noise",       0, &config_noisefilter, 1},
611
		{"verbose",            0, &config_verbose,     1},
Thomas White's avatar
Thomas White committed
612
		{"pdb",                1, NULL,               'p'},
613
		{"prefix",             1, NULL,               'x'},
614
615
		{"no-sat-corr",        0, &config_satcorr,     0},
		{"sat-corr",           0, &config_satcorr,     1}, /* Compat */
616
		{"threshold",          1, NULL,               't'},
617
		{"min-gradient",       1, NULL,                4},
618
		{"min-snr",            1, NULL,               11},
619
		{"min-integration-snr",1, NULL,               12},
620
		{"no-check-prefix",    0, &config_checkprefix, 0},
621
		{"no-closer-peak",     0, &config_closer,      0},
622
		{"insane",             0, &config_insane,      1},
Thomas White's avatar
Thomas White committed
623
		{"image",              1, NULL,               'e'},
624
		{"basename",           0, &config_basename,    1},
Thomas White's avatar
Thomas White committed
625
		{"record",             1, NULL,                5},
626
627
		{"cpus",               1, NULL,                6},
		{"cpugroup",           1, NULL,                7},
628
		{"cpuoffset",          1, NULL,                8},
Thomas White's avatar
Thomas White committed
629
		{"bg-sub",             0, &config_bgsub,       1}, /* Compat */
630
		{"no-bg-sub",          0, &config_bgsub,       0},
631
		{"hdf5-peaks",         1, NULL,                9},
632
		{"copy-hdf5-field",    1, NULL,               10},
Thomas White's avatar
Thomas White committed
633
634
635
636
		{0, 0, NULL, 0}
	};

	/* Short options */
Thomas White's avatar
Thomas White committed
637
	while ((c = getopt_long(argc, argv, "hi:wp:j:x:g:t:o:b:e:",
Thomas White's avatar
Thomas White committed
638
	                        longopts, NULL)) != -1) {
Thomas White's avatar
Thomas White committed
639
640

		switch (c) {
Thomas White's avatar
Thomas White committed
641
		case 'h' :
Thomas White's avatar
Thomas White committed
642
643
644
			show_help(argv[0]);
			return 0;

Thomas White's avatar
Thomas White committed
645
		case 'i' :
Thomas White's avatar
Thomas White committed
646
647
648
			filename = strdup(optarg);
			break;

649
650
651
652
		case 'o' :
			outfile = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
653
		case 'z' :
Thomas White's avatar
Thomas White committed
654
655
656
			indm_str = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
657
		case 'p' :
Thomas White's avatar
Thomas White committed
658
659
660
			pdb = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
661
		case 'x' :
662
663
664
			prefix = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
665
		case 'j' :
666
667
668
			nthreads = atoi(optarg);
			break;

Thomas White's avatar
Thomas White committed
669
		case 'g' :
Thomas White's avatar
Thomas White committed
670
671
672
			geometry = strdup(optarg);
			break;

673
674
675
676
		case 't' :
			threshold = strtof(optarg, NULL);
			break;

677
678
679
680
681
682
683
684
685
		case 'b' :
			beam = get_beam_parameters(optarg);
			if ( beam == NULL ) {
				ERROR("Failed to load beam parameters"
				      " from '%s'\n", optarg);
				return 1;
			}
			break;

686
687
688
689
		case 2 :
			speaks = strdup(optarg);
			break;

690
691
692
693
		case 3 :
			scellr = strdup(optarg);
			break;

694
695
696
697
		case 13 :
			toler = strdup(optarg);
			break;

698
699
700
701
		case 4 :
			min_gradient = strtof(optarg, NULL);
			break;

702
703
704
705
		case 11 :
			min_snr = strtof(optarg, NULL);
			break;

706
707
708
709
		case 12 :
			min_int_snr = strtof(optarg, NULL);
			break;

Thomas White's avatar
Thomas White committed
710
711
712
713
		case 'e' :
			element = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
714
715
716
717
718
		case 5 :
			stream_flags = parse_stream_flags(optarg);
			if ( stream_flags < 0 ) return 1;
			break;

719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
		case 6 :
			cpu_num = strtol(optarg, &endptr, 10);
			if ( !( (optarg[0] != '\0') && (endptr[0] == '\0') ) ) {
				ERROR("Invalid number of CPUs ('%s')\n",
				      optarg);
				return 1;
			}
			break;

		case 7 :
			cpu_groupsize = strtol(optarg, &endptr, 10);
			if ( !( (optarg[0] != '\0') && (endptr[0] == '\0') ) ) {
				ERROR("Invalid CPU group size ('%s')\n",
				      optarg);
				return 1;
			}
			if ( cpu_groupsize < 1 ) {
				ERROR("CPU group size cannot be"
				      " less than 1.\n");
				return 1;
			}
			break;

		case 8 :
			cpu_offset = strtol(optarg, &endptr, 10);
			if ( !( (optarg[0] != '\0') && (endptr[0] == '\0') ) ) {
				ERROR("Invalid CPU offset ('%s')\n",
				      optarg);
				return 1;
			}
			if ( cpu_offset < 0 ) {
				ERROR("CPU offset must be positive.\n");
				return 1;
			}
			break;

755
756
757
758
		case 9 :
			hdf5_peak_path = strdup(optarg);
			break;

759
760
761
762
		case 10 :
			add_copy_hdf5_field(copyme, optarg);
			break;

Thomas White's avatar
Thomas White committed
763
		case 0 :
Thomas White's avatar
Thomas White committed
764
765
			break;

Thomas White's avatar
Thomas White committed
766
		default :
Thomas White's avatar
Thomas White committed
767
768
769
770
771
			return 1;
		}

	}

772
773
774
775
776
777
778
	if ( (cpu_num > 0) && (cpu_num % cpu_groupsize != 0) ) {
		ERROR("Number of CPUs must be divisible by"
		      " the CPU group size.\n");
		return 1;
	}


Thomas White's avatar
Thomas White committed
779
780
781
782
783
784
785
786
787
	if ( filename == NULL ) {
		filename = strdup("-");
	}
	if ( strcmp(filename, "-") == 0 ) {
		fh = stdin;
	} else {
		fh = fopen(filename, "r");
	}
	if ( fh == NULL ) {
Thomas White's avatar
Thomas White committed
788
		ERROR("Failed to open input file '%s'\n", filename);
Thomas White's avatar
Thomas White committed
789
790
		return 1;
	}
Thomas White's avatar
Thomas White committed
791
	free(filename);
Thomas White's avatar
Thomas White committed
792

793
794
795
796
797
798
799
800
801
802
803
804
805
806
	if ( outfile == NULL ) {
		outfile = strdup("-");
	}
	if ( strcmp(outfile, "-") == 0 ) {
		ofh = stdout;
	} else {
		ofh = fopen(outfile, "w");
	}
	if ( ofh == NULL ) {
		ERROR("Failed to open output file '%s'\n", outfile);
		return 1;
	}
	free(outfile);

807
808
809
810
	if ( hdf5_peak_path == NULL ) {
		hdf5_peak_path = strdup("/processing/hitfinder/peakinfo");
	}

811
812
813
814
815
816
817
818
819
820
821
822
823
	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;
	} else if ( strcmp(speaks, "hdf5") == 0 ) {
		peaks = PEAK_HDF5;
	} else {
		ERROR("Unrecognised peak detection method '%s'\n", speaks);
		return 1;
	}
Thomas White's avatar
Thomas White committed
824
	free(speaks);
825

Thomas White's avatar
Thomas White committed
826
827
828
829
	if ( pdb == NULL ) {
		pdb = strdup("molecule.pdb");
	}

830
	if ( prefix == NULL ) {
Thomas White's avatar
Thomas White committed
831
		prefix = strdup("");
832
	} else {
833
834
835
		if ( config_checkprefix ) {
			prefix = check_prefix(prefix);
		}
836
837
	}

838
	if ( nthreads == 0 ) {
839
840
841
842
		ERROR("Invalid number of threads.\n");
		return 1;
	}

843
844
	if ( (indm_str == NULL) ||
	     ((indm_str != NULL) && (strcmp(indm_str, "none") == 0)) ) {
845
		STATUS("Not indexing anything.\n");
846
		indexer_needs_cell = 0;
847
848
849
		reduction_needs_cell = 0;
		indm = NULL;
		cellr = CELLR_NONE;
Thomas White's avatar
Thomas White committed
850
	} else {
851
852
853
854
855
856
857
858
		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 {
Thomas White's avatar
Thomas White committed
859
860
			indm = build_indexer_list(indm_str,
			                          &indexer_needs_cell);
861
862
863
864
865
			if ( indm == NULL ) {
				ERROR("Invalid indexer list '%s'\n", indm_str);
				return 1;
			}
			free(indm_str);
Thomas White's avatar
Thomas White committed
866
		}
Thomas White's avatar
Thomas White committed
867

868
869
870
871
872
873
874
		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 {
875
876
877
878
879
880
881
882
883
			int err;
			cellr = parse_cell_reduction(scellr, &err,
			                             &reduction_needs_cell);
			if ( err ) {
				ERROR("Unrecognised cell reduction '%s'\n",
			              scellr);
				return 1;
			}
			free(scellr);
884
		}
885
886
	}

887
888
889
	/* No indexing -> no reduction */
	if ( indm == NULL ) reduction_needs_cell = 0;

890
891
	if ( toler != NULL ) {
		int ttt;
Thomas White's avatar
Thomas White committed
892
893
		ttt = sscanf(toler, "%f,%f,%f,%f",
		             &tols[0], &tols[1], &tols[2], &tols[3] );
894
895
896
897
898
899
		if ( ttt != 4 ) {
			ERROR("Invalid parameters for '--tolerance'\n");
			return 1;
		}
	}

Thomas White's avatar
Thomas White committed
900
901
902
903
904
905
906
907
908
909
910
911
	if ( geometry == NULL ) {
		ERROR("You need to specify a geometry file with --geometry\n");
		return 1;
	}

	det = get_detector_geometry(geometry);
	if ( det == NULL ) {
		ERROR("Failed to read detector geometry from '%s'\n", geometry);
		return 1;
	}
	free(geometry);

Thomas White's avatar
Thomas White committed
912
	if ( reduction_needs_cell || indexer_needs_cell ) {
913
914
		cell = load_cell_from_pdb(pdb);
		if ( cell == NULL ) {
915
916
917
			ERROR("Couldn't read unit cell (from %s)\n", pdb);
			return 1;
		}
918
	} else {
Thomas White's avatar
Thomas White committed
919
920
		STATUS("No cell needed for these choices of indexing"
		       " and reduction.\n");
921
		cell = NULL;
922
	}
923
	free(pdb);
924

Thomas White's avatar
Thomas White committed
925
	write_stream_header(ofh, argc, argv);
926

927
928
929
930
931
932
933
934
	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");
		nominal_photon_energy = 2000.0;
	}

935
	/* Get first filename and use it to set up the indexing */
936
	prepare_line = malloc(1024*sizeof(char));
937
938
939
940
	rval = fgets(prepare_line, 1023, fh);
	if ( rval == NULL ) {
		ERROR("Failed to get filename to prepare indexing.\n");
		return 1;
941
	}
942
	chomp(prepare_line);
943
944
945
946
947
948
	if ( config_basename ) {
		char *tmp;
		tmp = safe_basename(prepare_line);
		free(prepare_line);
		prepare_line = tmp;
	}
949
950
951
952
	snprintf(prepare_filename, 1023, "%s%s", prefix, prepare_line);
	qargs.use_this_one_instead = prepare_line;

	/* Prepare the indexer */
Thomas White's avatar
Thomas White committed
953
954
955
956
957
958
959
960
961
	if ( indm != NULL ) {
		ipriv = prepare_indexing(indm, cell, prepare_filename, det,
		                         nominal_photon_energy);
		if ( ipriv == NULL ) {
			ERROR("Failed to prepare indexing.\n");
			return 1;
		}
	} else {
		ipriv = NULL;
962
963
	}

Thomas White's avatar
Thomas White committed
964
	gsl_set_error_handler_off();
Thomas White's avatar
Thomas White committed
965

966
967
968
969
970
971
	qargs.static_args.cell = cell;
	qargs.static_args.config_cmfilter = config_cmfilter;
	qargs.static_args.config_noisefilter = config_noisefilter;
	qargs.static_args.config_verbose = config_verbose;
	qargs.static_args.config_satcorr = config_satcorr;
	qargs.static_args.config_closer = config_closer;
972
	qargs.static_args.config_insane = config_insane;
973
	qargs.static_args.config_bgsub = config_bgsub;
974
	qargs.static_args.cellr = cellr;
975
976
977
978
	qargs.static_args.tols[0] = tols[0];
	qargs.static_args.tols[1] = tols[1];
	qargs.static_args.tols[2] = tols[2];
	qargs.static_args.tols[3] = tols[3];
979
	qargs.static_args.threshold = threshold;
980
	qargs.static_args.min_gradient = min_gradient;
981
	qargs.static_args.min_snr = min_snr;
Thomas White's avatar
Thomas White committed
982
	qargs.static_args.min_int_snr = min_int_snr;
983
984
985
986
987
988
	qargs.static_args.det = det;
	qargs.static_args.indm = indm;
	qargs.static_args.ipriv = ipriv;
	qargs.static_args.peaks = peaks;
	qargs.static_args.output_mutex = &output_mutex;
	qargs.static_args.ofh = ofh;
989
	qargs.static_args.beam = beam;
Thomas White's avatar
Thomas White committed
990
	qargs.static_args.element = element;
991
	qargs.static_args.stream_flags = stream_flags;
992
	qargs.static_args.hdf5_peak_path = hdf5_peak_path;
993
	qargs.static_args.copyme = copyme;
994
995
996

	qargs.fh = fh;
	qargs.prefix = prefix;
997
	qargs.config_basename = config_basename;
998
	qargs.n_indexable = 0;
999
1000
	qargs.n_processed = 0;
	qargs.n_indexable_last_stats = 0;