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},
Thomas White's avatar
Thomas White committed
603
		{"indexing",           1, NULL,               'z'},
Thomas White's avatar
Thomas White committed
604
		{"geometry",           1, NULL,               'g'},
605
		{"beam",               1, NULL,               'b'},
606
607
		{"filter-cm",          0, &config_cmfilter,    1},
		{"filter-noise",       0, &config_noisefilter, 1},
608
		{"verbose",            0, &config_verbose,     1},
Thomas White's avatar
Thomas White committed
609
		{"pdb",                1, NULL,               'p'},
610
		{"prefix",             1, NULL,               'x'},
611
612
		{"no-sat-corr",        0, &config_satcorr,     0},
		{"sat-corr",           0, &config_satcorr,     1}, /* Compat */
613
		{"threshold",          1, NULL,               't'},
614
		{"no-check-prefix",    0, &config_checkprefix, 0},
615
		{"no-closer-peak",     0, &config_closer,      0},
616
		{"insane",             0, &config_insane,      1},
Thomas White's avatar
Thomas White committed
617
		{"image",              1, NULL,               'e'},
618
		{"basename",           0, &config_basename,    1},
Thomas White's avatar
Thomas White committed
619
620
621
622
623
624
		{"bg-sub",             0, &config_bgsub,       1}, /* Compat */
		{"no-bg-sub",          0, &config_bgsub,       0},

		{"peaks",              1, NULL,                2},
		{"cell-reduction",     1, NULL,                3},
		{"min-gradient",       1, NULL,                4},
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},
629
		{"hdf5-peaks",         1, NULL,                9},
630
		{"copy-hdf5-field",    1, NULL,               10},
Thomas White's avatar
Thomas White committed
631
632
633
		{"min-snr",            1, NULL,               11},
		{"min-integration-snr",1, NULL,               12},
		{"tolerance",          1, NULL,               13},
Thomas White's avatar
Thomas White committed
634
635
636
637
		{0, 0, NULL, 0}
	};

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

			case 'h' :
Thomas White's avatar
Thomas White committed
644
645
646
			show_help(argv[0]);
			return 0;

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

Thomas White's avatar
Thomas White committed
651
			case 'o' :
652
653
654
			outfile = strdup(optarg);
			break;

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

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

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

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

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

Thomas White's avatar
Thomas White committed
675
			case 't' :
676
677
678
			threshold = strtof(optarg, NULL);
			break;

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

Thomas White's avatar
Thomas White committed
688
689
			case 'e' :
			element = strdup(optarg);
690
691
			break;

Thomas White's avatar
Thomas White committed
692
693
			case 2 :
			speaks = strdup(optarg);
694
695
			break;

Thomas White's avatar
Thomas White committed
696
697
			case 3 :
			scellr = strdup(optarg);
698
699
			break;

Thomas White's avatar
Thomas White committed
700
			case 4 :
701
702
703
			min_gradient = strtof(optarg, NULL);
			break;

Thomas White's avatar
Thomas White committed
704
			case 5 :
Thomas White's avatar
Thomas White committed
705
706
707
708
			stream_flags = parse_stream_flags(optarg);
			if ( stream_flags < 0 ) return 1;
			break;

Thomas White's avatar
Thomas White committed
709
			case 6 :
710
711
712
713
714
715
716
717
			cpu_num = strtol(optarg, &endptr, 10);
			if ( !( (optarg[0] != '\0') && (endptr[0] == '\0') ) ) {
				ERROR("Invalid number of CPUs ('%s')\n",
				      optarg);
				return 1;
			}
			break;

Thomas White's avatar
Thomas White committed
718
			case 7 :
719
720
721
722
723
724
725
726
727
728
729
730
731
			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;

Thomas White's avatar
Thomas White committed
732
			case 8 :
733
734
735
736
737
738
739
740
741
742
743
744
			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;

Thomas White's avatar
Thomas White committed
745
			case 9 :
746
747
748
			hdf5_peak_path = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
749
			case 10 :
750
751
752
			add_copy_hdf5_field(copyme, optarg);
			break;

Thomas White's avatar
Thomas White committed
753
754
755
756
757
758
			case 11 :
			min_snr = strtof(optarg, NULL);
			break;

			case 12 :
			min_int_snr = strtof(optarg, NULL);
Thomas White's avatar
Thomas White committed
759
760
			break;

Thomas White's avatar
Thomas White committed
761
762
763
764
765
766
767
768
			case 13 :
			toler = strdup(optarg);
			break;

			case 0 :
			break;

			default :
Thomas White's avatar
Thomas White committed
769
			return 1;
Thomas White's avatar
Thomas White committed
770

Thomas White's avatar
Thomas White committed
771
772
773
774
		}

	}

775
776
777
778
779
780
781
	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
782
783
784
785
786
787
788
789
790
	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
791
		ERROR("Failed to open input file '%s'\n", filename);
Thomas White's avatar
Thomas White committed
792
793
		return 1;
	}
Thomas White's avatar
Thomas White committed
794
	free(filename);
Thomas White's avatar
Thomas White committed
795

796
797
798
799
800
801
802
803
804
805
806
807
808
809
	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);

810
811
812
813
	if ( hdf5_peak_path == NULL ) {
		hdf5_peak_path = strdup("/processing/hitfinder/peakinfo");
	}

814
815
816
817
818
819
820
821
822
823
824
825
826
	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
827
	free(speaks);
828

Thomas White's avatar
Thomas White committed
829
830
831
832
	if ( pdb == NULL ) {
		pdb = strdup("molecule.pdb");
	}

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

841
	if ( nthreads == 0 ) {
842
843
844
845
		ERROR("Invalid number of threads.\n");
		return 1;
	}

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

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

890
891
892
	/* No indexing -> no reduction */
	if ( indm == NULL ) reduction_needs_cell = 0;

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

Thomas White's avatar
Thomas White committed
903
904
905
906
907
908
909
910
911
912
913
914
	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
915
	if ( reduction_needs_cell || indexer_needs_cell ) {
916
917
		cell = load_cell_from_pdb(pdb);
		if ( cell == NULL ) {
918
919
920
			ERROR("Couldn't read unit cell (from %s)\n", pdb);
			return 1;
		}
921
	} else {
Thomas White's avatar
Thomas White committed
922
923
		STATUS("No cell needed for these choices of indexing"
		       " and reduction.\n");
924
		cell = NULL;
925
	}
926
	free(pdb);
927

Thomas White's avatar
Thomas White committed
928
	write_stream_header(ofh, argc, argv);
929

930
931
932
933
934
935
936
937
	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;
	}

938
	/* Get first filename and use it to set up the indexing */
939
	prepare_line = malloc(1024*sizeof(char));
940
941
942
943
	rval = fgets(prepare_line, 1023, fh);
	if ( rval == NULL ) {
		ERROR("Failed to get filename to prepare indexing.\n");
		return 1;
944
	}
945
	chomp(prepare_line);
946
947
948
949
950
951
	if ( config_basename ) {
		char *tmp;
		tmp = safe_basename(prepare_line);
		free(prepare_line);
		prepare_line = tmp;
	}
952
953
954
955
	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
956
957
958
959
960
961
962
963
964
	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;
965
966
	}

Thomas White's avatar
Thomas White committed
967
	gsl_set_error_handler_off();
Thomas White's avatar
Thomas White committed
968

969
970
971
972
973
974
	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;
975
	qargs.static_args.config_insane = config_insane;
976
	qargs.static_args.config_bgsub = config_bgsub;
977
	qargs.static_args.cellr = cellr;
978
979
980
981
	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];
982
	qargs.static_args.threshold = threshold;
983
	qargs.static_args.min_gradient = min_gradient;
984
	qargs.static_args.min_snr = min_snr;
Thomas White's avatar
Thomas White committed
985
	qargs.static_args.min_int_snr = min_int_snr;
986
987
988
989
990
991
	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;
992
	qargs.static_args.beam = beam;
Thomas White's avatar
Thomas White committed
993
	qargs.static_args.element = element;
994
	qargs.static_args.stream_flags = stream_flags;
995
	qargs.static_args.hdf5_peak_path = hdf5_peak_path;
996
	qargs.static_args.copyme = copyme;
997
998
999

	qargs.fh = fh;
	qargs.prefix = prefix;
1000
	qargs.config_basename = config_basename;