indexamajig.c 21.6 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
 * (c) 2006-2011 Thomas White <taw@physics.org>
Thomas White's avatar
Thomas White committed
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
 *
 * Part of CrystFEL - crystallography with a FEL
 *
 */


#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
23
#include <hdf5.h>
Thomas White's avatar
Thomas White committed
24
#include <gsl/gsl_errno.h>
25
#include <pthread.h>
26
#include <time.h>
Thomas White's avatar
Thomas White committed
27
28
29

#include "utils.h"
#include "hdf5-file.h"
Thomas White's avatar
Thomas White committed
30
#include "index.h"
31
#include "peaks.h"
32
#include "detector.h"
Thomas White's avatar
Thomas White committed
33
#include "filters.h"
34
#include "thread-pool.h"
35
#include "beam-parameters.h"
36
#include "geometry.h"
Thomas White's avatar
Thomas White committed
37
#include "stream.h"
Thomas White's avatar
Thomas White committed
38
#include "reflist-utils.h"
39

40

Thomas White's avatar
Thomas White committed
41
/* Write statistics at APPROXIMATELY this interval */
42
43
44
#define STATS_EVERY_N_SECONDS (5)


45
46
47
48
49
50
enum {
	PEAK_ZAEF,
	PEAK_HDF5,
};


51
52
/* Information about the indexing process which is common to all patterns */
struct static_index_args
53
54
55
56
57
{
	UnitCell *cell;
	int config_cmfilter;
	int config_noisefilter;
	int config_verbose;
Thomas White's avatar
Thomas White committed
58
	int stream_flags;         /* What goes into the output? */
Thomas White's avatar
Thomas White committed
59
	int config_polar;
60
	int config_satcorr;
61
	int config_closer;
62
	int config_insane;
63
	float threshold;
64
	float min_gradient;
Thomas White's avatar
Thomas White committed
65
	struct detector *det;
Thomas White's avatar
Thomas White committed
66
67
	IndexingMethod *indm;
	IndexingPrivate **ipriv;
Thomas White's avatar
Thomas White committed
68
	int peaks;                /* Peak detection method */
69
	int cellr;
70
	struct beam_params *beam;
Thomas White's avatar
Thomas White committed
71
	const char *element;
72

73
74
75
	/* Output stream */
	pthread_mutex_t *output_mutex;  /* Protects the output stream */
	FILE *ofh;
76
77
78
};


79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
/* 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;
96
	int config_basename;
97
98
	struct static_index_args static_args;

99
	char *use_this_one_instead;
100
101
102
103
104
105

	int n_indexable;
	int n_processed;
	int n_indexable_last_stats;
	int n_processed_last_stats;
	int t_last_stats;
106
107
108
};


Thomas White's avatar
Thomas White committed
109
110
111
112
113
114
static void show_help(const char *s)
{
	printf("Syntax: %s [options]\n\n", s);
	printf(
"Process and index FEL diffraction images.\n"
"\n"
115
" -h, --help               Display this help message.\n"
Thomas White's avatar
Thomas White committed
116
"\n"
117
" -i, --input=<filename>   Specify file containing list of images to process.\n"
Thomas White's avatar
Thomas White committed
118
"                           '-' means stdin, which is the default.\n"
Thomas White's avatar
Thomas White committed
119
120
" -o, --output=<filename>  Write output stream to this file. '-' for stdout.\n"
"                           Default: indexamajig.stream\n"
Thomas White's avatar
Thomas White committed
121
"\n"
Thomas White's avatar
Thomas White committed
122
123
124
125
126
"     --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"
127
" -g. --geometry=<file>    Get detector geometry from file.\n"
128
129
130
" -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"
131
" -p, --pdb=<file>         PDB file from which to get the unit cell to match.\n"
Thomas White's avatar
Thomas White committed
132
"                           Default: 'molecule.pdb'.\n"
133
"     --basename           Remove the directory parts of the filenames.\n"
134
" -x, --prefix=<p>         Prefix filenames from input file with <p>.\n"
135
136
137
138
139
"     --peaks=<method>     Use 'method' for finding peaks.  Choose from:\n"
"                           zaef  : Use Zaefferer (2000) gradient detection.\n"
"                                    This is the default method.\n"
"                           hdf5  : Get from /processing/hitfinder/peakinfo\n"
"                                    in the HDF5 file.\n"
Thomas White's avatar
Thomas White committed
140
141
"\n\n"
"You can control what information is included in the output stream using\n"
Thomas White's avatar
Thomas White committed
142
"' --record=<flags>'.  Possible flags are:\n\n"
143
144
145
" pixels            Include a list of sums of pixel values within the\n"
"                    integration domain, correcting for individual pixel\n"
"                    solid angles.\n"
Thomas White's avatar
Thomas White committed
146
"\n"
147
148
" integrated        Include a list of reflection intensities, produced by\n"
"                    integrating around predicted peak locations.\n"
Thomas White's avatar
Thomas White committed
149
"\n"
150
151
" peaks             Include peak locations and intensities from the peak\n"
"                    search.\n"
Thomas White's avatar
Thomas White committed
152
"\n"
153
" peaksifindexed    As 'peaks', but only if the pattern could be indexed.\n"
Thomas White's avatar
Thomas White committed
154
"\n"
155
156
" peaksifnotindexed As 'peaks', but only if the pattern could NOT be indexed.\n"
"\n\n"
Thomas White's avatar
Thomas White committed
157
158
"The default is '--record=integrated'.  The flags 'pixels' and 'integrated'\n"
"are mutually exclusive, as are the flags 'peaks' and 'peaksifindexed'.\n"
Thomas White's avatar
Thomas White committed
159
160
"\n\n"
"For more control over the process, you might need:\n\n"
Thomas White's avatar
Thomas White committed
161
162
163
164
165
"     --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"
166
"     --filter-cm          Perform common-mode noise subtraction on images\n"
Thomas White's avatar
Thomas White committed
167
168
"                           before proceeding.  Intensities will be extracted\n"
"                           from the image as it is after this processing.\n"
169
"     --filter-noise       Apply an aggressive noise filter which sets all\n"
Thomas White's avatar
Thomas White committed
170
171
172
"                           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"
173
"     --unpolarized        Don't correct for the polarisation of the X-rays.\n"
174
175
"     --no-sat-corr        Don't correct values of saturated peaks using a\n"
"                           table included in the HDF5 file.\n"
176
"     --threshold=<n>      Only accept peaks above <n> ADU.  Default: 800.\n"
177
178
"     --min-gradient=<n>   Minimum gradient for Zaefferer peak search.\n"
"                           Default: 100,000.\n"
Thomas White's avatar
Thomas White committed
179
" -e, --image=<element>    Use this image from the HDF5 file.\n"
Thomas White's avatar
Thomas White committed
180
"                           Example: /data/data0.\n"
Thomas White's avatar
Thomas White committed
181
"                           Default: The first one found.\n"
Thomas White's avatar
Thomas White committed
182
"\n"
183
184
185
"\nOptions for greater performance or verbosity:\n\n"
"     --verbose            Be verbose about indexing.\n"
" -j <n>                   Run <n> analyses in parallel.  Default 1.\n"
186
187
188
"\n"
"\nOptions you probably won't need:\n\n"
"     --no-check-prefix    Don't attempt to correct the --prefix.\n"
189
190
191
"     --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"
192
193
"     --insane             Don't check that the reduced cell accounts for at\n"
"                           least 10%% of the located peaks.\n"
Thomas White's avatar
Thomas White committed
194
);
Thomas White's avatar
Thomas White committed
195
196
197
}


198
static void process_image(void *pp, int cookie)
199
{
200
	struct index_args *pargs = pp;
201
202
203
204
	struct hdfile *hdfile;
	struct image image;
	float *data_for_measurement;
	size_t data_size;
205
	char *filename = pargs->filename;
206
207
208
209
210
	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;
	int config_polar = pargs->static_args.config_polar;
Thomas White's avatar
Thomas White committed
211
	IndexingMethod *indm = pargs->static_args.indm;
212
	const struct beam_params *beam = pargs->static_args.beam;
213
214
215
216

	image.features = NULL;
	image.data = NULL;
	image.indexed_cell = NULL;
217
	image.id = cookie;
218
	image.filename = filename;
219
	image.det = copy_geom(pargs->static_args.det);
220

Thomas White's avatar
Thomas White committed
221
	pargs->indexable = 0;
222

223
	hdfile = hdfile_open(filename);
Thomas White's avatar
Thomas White committed
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
	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;
		}

247
248
	}

249
250
251
252
253
254
255
256
257
258
259
	hdf5_read(hdfile, &image, pargs->static_args.config_satcorr);
	if ( image.lambda < 0.0 ) {
		if ( beam != NULL ) {
			image.lambda = beam->photon_energy;
		} else {
			ERROR("No wavelength in file, so you need to give "
			      "a beam parameters file with -b.\n");
			hdfile_close(hdfile);
			return;
		}
	}
260
	fill_in_values(image.det, hdfile);
261
262
263
264
265
266
267
268
269
270
271
272
273

	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 {
274
		memcpy(data_for_measurement, image.data, data_size);
275
276
	}

277
	switch ( pargs->static_args.peaks )
278
279
280
281
282
283
284
285
	{
	case PEAK_HDF5 :
		/* Get peaks from HDF5 */
		if ( get_peaks(&image, hdfile) ) {
			ERROR("Failed to get peaks from HDF5 file.\n");
		}
		break;
	case PEAK_ZAEF :
286
287
		search_peaks(&image, pargs->static_args.threshold,
		             pargs->static_args.min_gradient);
288
289
		break;
	}
Thomas White's avatar
Thomas White committed
290
291
292
293
294
295

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

296
	/* Calculate orientation matrix (by magic) */
297
	index_pattern(&image, cell, indm, pargs->static_args.cellr,
298
299
		      config_verbose, pargs->static_args.ipriv,
		      pargs->static_args.config_insane);
300
301

	/* No cell at this point?  Then we're done. */
Thomas White's avatar
Thomas White committed
302
	if ( image.indexed_cell != NULL ) pargs->indexable = 1;
303
304

	/* Measure intensities if requested */
305

Thomas White's avatar
Thomas White committed
306
	/* Do EITHER: */
307

Thomas White's avatar
Thomas White committed
308
309
310
311
312
	//image.div = beam->divergence;
	//image.bw = beam->bandwidth;
	//image.profile_radius = 0.0001e9;
	//image.reflections = find_intersections(&image,
	//                                       image.indexed_cell, 0);
313

314
315
316
317
	if ( image.indexed_cell != NULL ) {
		image.reflections = find_projected_peaks(&image,
			                                 image.indexed_cell,
			                                 0, 0.1);
318

319
320
		integrate_reflections(&image, config_polar,
			              pargs->static_args.config_closer);
321

322
		/* OR */
323

324
325
326
327
328
329
330
331
		//image.reflections = integrate_pixels(&image, 0, 0.1,
		//                                     config_polar);

	} else {

		image.reflections = NULL;

	}
Thomas White's avatar
Thomas White committed
332
333
334
335
336

	pthread_mutex_lock(pargs->static_args.output_mutex);
	write_chunk(pargs->static_args.ofh, &image,
	            pargs->static_args.stream_flags);
	pthread_mutex_unlock(pargs->static_args.output_mutex);
337
338

	/* Only free cell if found */
339
	cell_free(image.indexed_cell);
340
341

	free(image.data);
Thomas White's avatar
Thomas White committed
342
	free(image.flags);
343
344
	image_feature_list_free(image.features);
	hdfile_close(hdfile);
Thomas White's avatar
Thomas White committed
345
	free_detector_geometry(image.det);
346
347
348
}


349
static void *get_image(void *qp)
350
{
351
	char *line;
352
353
354
	struct index_args *pargs;
	char *rval;
	struct queue_args *qargs = qp;
355

356
	/* Initialise new task arguments */
357
358
359
	pargs = malloc(sizeof(struct index_args));
	memcpy(&pargs->static_args, &qargs->static_args,
	       sizeof(struct static_index_args));
360

361
362
363
	/* Get the next filename */
	if ( qargs->use_this_one_instead != NULL ) {

364
		line = qargs->use_this_one_instead;
365
366
367
368
		qargs->use_this_one_instead = NULL;

	} else {

369
		line = malloc(1024*sizeof(char));
370
		rval = fgets(line, 1023, qargs->fh);
Thomas White's avatar
Thomas White committed
371
372
373
374
		if ( rval == NULL ) {
			free(pargs);
			return NULL;
		}
375
376
377
		chomp(line);

	}
378

379
380
381
382
383
384
385
386
387
388
389
390
391
	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);

392
393
	return pargs;
}
394
395


396
397
398
399
static void finalise_image(void *qp, void *pp)
{
	struct queue_args *qargs = qp;
	struct index_args *pargs = pp;
400
	struct timespec tp;
401

402
	qargs->n_indexable += pargs->indexable;
403
404
405
	qargs->n_processed++;

	clock_gettime(CLOCK_REALTIME, &tp);
Thomas White's avatar
Thomas White committed
406
	if ( tp.tv_sec >= qargs->t_last_stats+STATS_EVERY_N_SECONDS ) {
407
408

		STATUS("%i out of %i indexed so far,"
Thomas White's avatar
Thomas White committed
409
		       " %i out of %i since the last message.\n",
410
411
		       qargs->n_indexable, qargs->n_processed,
		       qargs->n_indexable - qargs->n_indexable_last_stats,
Thomas White's avatar
Thomas White committed
412
		       qargs->n_processed - qargs->n_processed_last_stats);
413
414
415
416
417
418

		qargs->n_processed_last_stats = qargs->n_processed;
		qargs->n_indexable_last_stats = qargs->n_indexable;
		qargs->t_last_stats = tp.tv_sec;

	}
419

420
421
	free(pargs->filename);
	free(pargs);
422
423
424
}


Thomas White's avatar
Thomas White committed
425
426
427
428
int main(int argc, char *argv[])
{
	int c;
	char *filename = NULL;
429
	char *outfile = NULL;
Thomas White's avatar
Thomas White committed
430
	FILE *fh;
431
	FILE *ofh;
432
	char *rval = NULL;
Thomas White's avatar
Thomas White committed
433
	int n_images;
434
	int config_noindex = 0;
435
436
	int config_cmfilter = 0;
	int config_noisefilter = 0;
437
	int config_verbose = 0;
Thomas White's avatar
Thomas White committed
438
	int config_polar = 1;
439
	int config_satcorr = 1;
440
	int config_checkprefix = 1;
441
	int config_closer = 1;
442
	int config_insane = 0;
443
	int config_basename = 0;
444
	float threshold = 800.0;
445
	float min_gradient = 100000.0;
Thomas White's avatar
Thomas White committed
446
447
	struct detector *det;
	char *geometry = NULL;
Thomas White's avatar
Thomas White committed
448
449
450
451
	IndexingMethod *indm;
	IndexingPrivate **ipriv;
	int indexer_needs_cell;
	int reduction_needs_cell;
Thomas White's avatar
Thomas White committed
452
	char *indm_str = NULL;
453
	UnitCell *cell;
Thomas White's avatar
Thomas White committed
454
	char *pdb = NULL;
455
	char *prefix = NULL;
456
	char *speaks = NULL;
457
458
	char *scellr = NULL;
	int cellr;
459
	int peaks;
460
461
	int nthreads = 1;
	int i;
462
	pthread_mutex_t output_mutex = PTHREAD_MUTEX_INITIALIZER;
463
	char *prepare_line;
464
	char prepare_filename[1024];
465
	struct queue_args qargs;
466
	struct beam_params *beam = NULL;
Thomas White's avatar
Thomas White committed
467
	char *element = NULL;
468
	double nominal_photon_energy;
Thomas White's avatar
Thomas White committed
469
	int stream_flags = STREAM_INTEGRATED;
470
	struct timespec tp;
Thomas White's avatar
Thomas White committed
471
472
473
474
475

	/* Long options */
	const struct option longopts[] = {
		{"help",               0, NULL,               'h'},
		{"input",              1, NULL,               'i'},
476
		{"output",             1, NULL,               'o'},
477
		{"no-index",           0, &config_noindex,     1},
478
		{"peaks",              1, NULL,                2},
479
		{"cell-reduction",     1, NULL,                3},
Thomas White's avatar
Thomas White committed
480
		{"indexing",           1, NULL,               'z'},
Thomas White's avatar
Thomas White committed
481
		{"geometry",           1, NULL,               'g'},
482
		{"beam",               1, NULL,               'b'},
483
484
		{"filter-cm",          0, &config_cmfilter,    1},
		{"filter-noise",       0, &config_noisefilter, 1},
485
		{"verbose",            0, &config_verbose,     1},
Thomas White's avatar
Thomas White committed
486
		{"pdb",                1, NULL,               'p'},
487
		{"prefix",             1, NULL,               'x'},
Thomas White's avatar
Thomas White committed
488
		{"unpolarized",        0, &config_polar,       0},
489
490
		{"no-sat-corr",        0, &config_satcorr,     0},
		{"sat-corr",           0, &config_satcorr,     1}, /* Compat */
491
		{"threshold",          1, NULL,               't'},
492
		{"min-gradient",       1, NULL,                4},
493
		{"no-check-prefix",    0, &config_checkprefix, 0},
494
		{"no-closer-peak",     0, &config_closer,      0},
495
		{"insane",             0, &config_insane,      1},
Thomas White's avatar
Thomas White committed
496
		{"image",              1, NULL,               'e'},
497
		{"basename",           0, &config_basename,    1},
Thomas White's avatar
Thomas White committed
498
		{"record",             1, NULL,                5},
Thomas White's avatar
Thomas White committed
499
500
501
502
		{0, 0, NULL, 0}
	};

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

		switch (c) {
Thomas White's avatar
Thomas White committed
507
		case 'h' :
Thomas White's avatar
Thomas White committed
508
509
510
			show_help(argv[0]);
			return 0;

Thomas White's avatar
Thomas White committed
511
		case 'i' :
Thomas White's avatar
Thomas White committed
512
513
514
			filename = strdup(optarg);
			break;

515
516
517
518
		case 'o' :
			outfile = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
519
		case 'z' :
Thomas White's avatar
Thomas White committed
520
521
522
			indm_str = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
523
		case 'p' :
Thomas White's avatar
Thomas White committed
524
525
526
			pdb = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
527
		case 'x' :
528
529
530
			prefix = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
531
		case 'j' :
532
533
534
			nthreads = atoi(optarg);
			break;

Thomas White's avatar
Thomas White committed
535
		case 'g' :
Thomas White's avatar
Thomas White committed
536
537
538
			geometry = strdup(optarg);
			break;

539
540
541
542
		case 't' :
			threshold = strtof(optarg, NULL);
			break;

543
544
545
546
547
548
549
550
551
		case 'b' :
			beam = get_beam_parameters(optarg);
			if ( beam == NULL ) {
				ERROR("Failed to load beam parameters"
				      " from '%s'\n", optarg);
				return 1;
			}
			break;

552
553
554
555
		case 2 :
			speaks = strdup(optarg);
			break;

556
557
558
559
		case 3 :
			scellr = strdup(optarg);
			break;

560
561
562
563
		case 4 :
			min_gradient = strtof(optarg, NULL);
			break;

Thomas White's avatar
Thomas White committed
564
565
566
567
		case 'e' :
			element = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
568
569
570
571
572
		case 5 :
			stream_flags = parse_stream_flags(optarg);
			if ( stream_flags < 0 ) return 1;
			break;

Thomas White's avatar
Thomas White committed
573
		case 0 :
Thomas White's avatar
Thomas White committed
574
575
			break;

Thomas White's avatar
Thomas White committed
576
		default :
Thomas White's avatar
Thomas White committed
577
578
579
580
581
582
583
584
585
586
587
588
589
590
			return 1;
		}

	}

	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
591
		ERROR("Failed to open input file '%s'\n", filename);
Thomas White's avatar
Thomas White committed
592
593
		return 1;
	}
Thomas White's avatar
Thomas White committed
594
	free(filename);
Thomas White's avatar
Thomas White committed
595

596
597
598
599
600
601
602
603
604
605
606
607
608
609
	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);

610
611
612
613
614
615
616
617
618
619
620
621
622
	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
623
	free(speaks);
624

Thomas White's avatar
Thomas White committed
625
626
627
628
	if ( pdb == NULL ) {
		pdb = strdup("molecule.pdb");
	}

629
	if ( prefix == NULL ) {
Thomas White's avatar
Thomas White committed
630
		prefix = strdup("");
631
	} else {
632
633
634
		if ( config_checkprefix ) {
			prefix = check_prefix(prefix);
		}
635
636
	}

637
	if ( nthreads == 0 ) {
638
639
640
641
		ERROR("Invalid number of threads.\n");
		return 1;
	}

642
643
	if ( (indm_str == NULL) ||
	     ((indm_str != NULL) && (strcmp(indm_str, "none") == 0)) ) {
644
		STATUS("Not indexing anything.\n");
645
		indexer_needs_cell = 0;
646
647
648
		reduction_needs_cell = 0;
		indm = NULL;
		cellr = CELLR_NONE;
Thomas White's avatar
Thomas White committed
649
	} else {
650
651
652
653
654
655
656
657
658
659
660
661
662
663
		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);
Thomas White's avatar
Thomas White committed
664
		}
Thomas White's avatar
Thomas White committed
665

666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
		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 if ( strcmp(scellr, "none") == 0 ) {
			cellr = CELLR_NONE;
		} else if ( strcmp(scellr, "reduce") == 0) {
			cellr = CELLR_REDUCE;
			reduction_needs_cell = 1;
		} else if ( strcmp(scellr, "compare") == 0) {
			cellr = CELLR_COMPARE;
			reduction_needs_cell = 1;
		} else {
			ERROR("Unrecognised cell reduction method '%s'\n",
			      scellr);
			return 1;
		}
		free(scellr);  /* free(NULL) is OK. */
686
687
	}

688
689
690
	/* No indexing -> no reduction */
	if ( indm == NULL ) reduction_needs_cell = 0;

Thomas White's avatar
Thomas White committed
691
692
693
694
695
696
697
698
699
700
701
702
	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
703
	if ( reduction_needs_cell || indexer_needs_cell ) {
704
705
		cell = load_cell_from_pdb(pdb);
		if ( cell == NULL ) {
706
707
708
			ERROR("Couldn't read unit cell (from %s)\n", pdb);
			return 1;
		}
709
	} else {
Thomas White's avatar
Thomas White committed
710
711
		STATUS("No cell needed for these choices of indexing"
		       " and reduction.\n");
712
		cell = NULL;
713
	}
714
	free(pdb);
715

716
	/* Start by writing the entire command line to stdout */
Thomas White's avatar
Thomas White committed
717
	fprintf(ofh, "CrystFEL stream format 2.0\n");
718
719
720
721
	fprintf(ofh, "Command line:");
	for ( i=0; i<argc; i++ ) {
		fprintf(ofh, " %s", argv[i]);
	}
722
723
	fprintf(ofh, "\n");
	fflush(ofh);
724

725
726
727
728
729
730
731
732
	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;
	}

733
	/* Get first filename and use it to set up the indexing */
734
	prepare_line = malloc(1024*sizeof(char));
735
736
737
738
	rval = fgets(prepare_line, 1023, fh);
	if ( rval == NULL ) {
		ERROR("Failed to get filename to prepare indexing.\n");
		return 1;
739
	}
740
	chomp(prepare_line);
741
742
743
744
745
746
	if ( config_basename ) {
		char *tmp;
		tmp = safe_basename(prepare_line);
		free(prepare_line);
		prepare_line = tmp;
	}
747
748
749
750
	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
751
752
753
754
755
756
757
758
759
	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;
760
761
	}

Thomas White's avatar
Thomas White committed
762
	gsl_set_error_handler_off();
Thomas White's avatar
Thomas White committed
763

764
765
766
767
768
769
770
	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_polar = config_polar;
	qargs.static_args.config_satcorr = config_satcorr;
	qargs.static_args.config_closer = config_closer;
771
	qargs.static_args.config_insane = config_insane;
772
	qargs.static_args.cellr = cellr;
773
	qargs.static_args.threshold = threshold;
774
	qargs.static_args.min_gradient = min_gradient;
775
776
777
778
779
780
	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;
781
	qargs.static_args.beam = beam;
Thomas White's avatar
Thomas White committed
782
	qargs.static_args.element = element;
783
	qargs.static_args.stream_flags = stream_flags;
784
785
786

	qargs.fh = fh;
	qargs.prefix = prefix;
787
	qargs.config_basename = config_basename;
788
	qargs.n_indexable = 0;
789
790
791
792
793
794
	qargs.n_processed = 0;
	qargs.n_indexable_last_stats = 0;
	qargs.n_processed_last_stats = 0;
	clock_gettime(CLOCK_REALTIME, &tp);
	qargs.t_last_stats = tp.tv_sec;

795
796
797

	n_images = run_threads(nthreads, process_image, get_image,
	                       finalise_image, &qargs, 0);
798

799
800
	cleanup_indexing(ipriv);

Thomas White's avatar
Thomas White committed
801
802
	free(indm);
	free(ipriv);
803
	free(prefix);
Thomas White's avatar
Thomas White committed
804
805
	free(det->panels);
	free(det);
Thomas White's avatar
Thomas White committed
806
	free(element);
807
	cell_free(cell);
Thomas White's avatar
Thomas White committed
808
809
	if ( fh != stdin ) fclose(fh);
	if ( ofh != stdout ) fclose(ofh);
Thomas White's avatar
Thomas White committed
810

811
812
	STATUS("There were %i images, of which %i could be indexed.\n",
	        n_images, qargs.n_indexable);
Thomas White's avatar
Thomas White committed
813
814
815

	return 0;
}