indexamajig.c 20.7 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>
Thomas White's avatar
Thomas White committed
26
#include <sys/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 "reflections.h"
35
#include "thread-pool.h"
36
#include "beam-parameters.h"
37
#include "geometry.h"
Thomas White's avatar
Thomas White committed
38
#include "stream.h"
39

40
41
42
43
44
45
46

enum {
	PEAK_ZAEF,
	PEAK_HDF5,
};


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

69
70
71
	/* Output stream */
	pthread_mutex_t *output_mutex;  /* Protects the output stream */
	FILE *ofh;
72
73
74
};


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

	int n_indexable;
96
97

	char *use_this_one_instead;
98
99
100
};


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


187
static void process_image(void *pp, int cookie)
188
{
189
	struct index_args *pargs = pp;
190
191
192
193
	struct hdfile *hdfile;
	struct image image;
	float *data_for_measurement;
	size_t data_size;
194
	char *filename = pargs->filename;
195
196
197
198
199
	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
200
	IndexingMethod *indm = pargs->static_args.indm;
201
	const struct beam_params *beam = pargs->static_args.beam;
202
203
204
205

	image.features = NULL;
	image.data = NULL;
	image.indexed_cell = NULL;
206
	image.id = cookie;
207
	image.filename = filename;
208
	image.det = copy_geom(pargs->static_args.det);
209

210
	STATUS("Processing '%s'\n", image.filename);
211

Thomas White's avatar
Thomas White committed
212
	pargs->indexable = 0;
213

214
	hdfile = hdfile_open(filename);
Thomas White's avatar
Thomas White committed
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
	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;
		}

238
239
	}

240
241
242
243
244
245
246
247
248
249
250
	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;
		}
	}
251
	fill_in_values(image.det, hdfile);
252
253
254
255
256
257
258
259
260
261
262
263
264

	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 {
265
		memcpy(data_for_measurement, image.data, data_size);
266
267
	}

268
	switch ( pargs->static_args.peaks )
269
270
271
272
273
274
275
276
	{
	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 :
277
278
		search_peaks(&image, pargs->static_args.threshold,
		             pargs->static_args.min_gradient);
279
280
		break;
	}
Thomas White's avatar
Thomas White committed
281
282
283
284
285
286

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

287
	/* Calculate orientation matrix (by magic) */
288
	index_pattern(&image, cell, indm, pargs->static_args.cellr,
289
290
		      config_verbose, pargs->static_args.ipriv,
		      pargs->static_args.config_insane);
291
292

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

	/* Measure intensities if requested */
296

Thomas White's avatar
Thomas White committed
297
	/* Do EITHER: */
298

Thomas White's avatar
Thomas White committed
299
300
301
302
303
	//image.div = beam->divergence;
	//image.bw = beam->bandwidth;
	//image.profile_radius = 0.0001e9;
	//image.reflections = find_intersections(&image,
	//                                       image.indexed_cell, 0);
304

Thomas White's avatar
Thomas White committed
305
306
307
	image.reflections = find_projected_peaks(&image,
	                                         image.indexed_cell,
	                                         0, 0.1);
308

Thomas White's avatar
Thomas White committed
309
310
	integrate_reflections(&image, config_polar,
	                      pargs->static_args.config_closer);
311

Thomas White's avatar
Thomas White committed
312
	/* OR */
313

Thomas White's avatar
Thomas White committed
314
315
316
317
318
319
320
	//image.reflections = integrate_pixels(&image, 0, 0.1,
	//                                     config_polar);

	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);
321
322

	/* Only free cell if found */
323
	cell_free(image.indexed_cell);
324
325

	free(image.data);
Thomas White's avatar
Thomas White committed
326
	free(image.flags);
327
328
	image_feature_list_free(image.features);
	hdfile_close(hdfile);
Thomas White's avatar
Thomas White committed
329
	free_detector_geometry(image.det);
330
331
332
}


333
static void *get_image(void *qp)
334
{
335
	char *line;
336
337
338
	struct index_args *pargs;
	char *rval;
	struct queue_args *qargs = qp;
339

340
	/* Initialise new task arguments */
341
342
343
	pargs = malloc(sizeof(struct index_args));
	memcpy(&pargs->static_args, &qargs->static_args,
	       sizeof(struct static_index_args));
344

345
346
347
	/* Get the next filename */
	if ( qargs->use_this_one_instead != NULL ) {

348
		line = qargs->use_this_one_instead;
349
350
351
352
		qargs->use_this_one_instead = NULL;

	} else {

353
		line = malloc(1024*sizeof(char));
354
		rval = fgets(line, 1023, qargs->fh);
Thomas White's avatar
Thomas White committed
355
356
357
358
		if ( rval == NULL ) {
			free(pargs);
			return NULL;
		}
359
360
361
		chomp(line);

	}
362

363
364
365
366
367
368
369
370
371
372
373
374
375
	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);

376
377
	return pargs;
}
378
379


380
381
382
383
static void finalise_image(void *qp, void *pp)
{
	struct queue_args *qargs = qp;
	struct index_args *pargs = pp;
384

385
	qargs->n_indexable += pargs->indexable;
386

387
388
	free(pargs->filename);
	free(pargs);
389
390
391
}


Thomas White's avatar
Thomas White committed
392
393
394
395
int main(int argc, char *argv[])
{
	int c;
	char *filename = NULL;
396
	char *outfile = NULL;
Thomas White's avatar
Thomas White committed
397
	FILE *fh;
398
	FILE *ofh;
399
	char *rval = NULL;
Thomas White's avatar
Thomas White committed
400
	int n_images;
401
	int config_noindex = 0;
Thomas White's avatar
Thomas White committed
402
	int config_dumpfound = 0;
403
	int config_nearbragg = 0;
404
405
	int config_cmfilter = 0;
	int config_noisefilter = 0;
406
	int config_verbose = 0;
Thomas White's avatar
Thomas White committed
407
	int config_polar = 1;
408
	int config_satcorr = 1;
409
	int config_checkprefix = 1;
410
	int config_closer = 1;
411
	int config_insane = 0;
412
	int config_basename = 0;
413
	float threshold = 800.0;
414
	float min_gradient = 100000.0;
Thomas White's avatar
Thomas White committed
415
416
	struct detector *det;
	char *geometry = NULL;
Thomas White's avatar
Thomas White committed
417
418
419
420
	IndexingMethod *indm;
	IndexingPrivate **ipriv;
	int indexer_needs_cell;
	int reduction_needs_cell;
Thomas White's avatar
Thomas White committed
421
	char *indm_str = NULL;
422
	UnitCell *cell;
Thomas White's avatar
Thomas White committed
423
	char *pdb = NULL;
424
	char *prefix = NULL;
425
	char *speaks = NULL;
426
427
	char *scellr = NULL;
	int cellr;
428
	int peaks;
429
430
	int nthreads = 1;
	int i;
431
	pthread_mutex_t output_mutex = PTHREAD_MUTEX_INITIALIZER;
432
	char *prepare_line;
433
	char prepare_filename[1024];
434
	struct queue_args qargs;
435
	struct beam_params *beam = NULL;
Thomas White's avatar
Thomas White committed
436
	char *element = NULL;
437
	double nominal_photon_energy;
Thomas White's avatar
Thomas White committed
438
	int stream_flags = STREAM_INTEGRATED;
Thomas White's avatar
Thomas White committed
439
440
441
442
443

	/* Long options */
	const struct option longopts[] = {
		{"help",               0, NULL,               'h'},
		{"input",              1, NULL,               'i'},
444
		{"output",             1, NULL,               'o'},
445
		{"no-index",           0, &config_noindex,     1},
Thomas White's avatar
Thomas White committed
446
		{"dump-peaks",         0, &config_dumpfound,   1},
447
		{"peaks",              1, NULL,                2},
448
		{"cell-reduction",     1, NULL,                3},
449
		{"near-bragg",         0, &config_nearbragg,   1},
Thomas White's avatar
Thomas White committed
450
		{"indexing",           1, NULL,               'z'},
Thomas White's avatar
Thomas White committed
451
		{"geometry",           1, NULL,               'g'},
452
		{"beam",               1, NULL,               'b'},
453
454
		{"filter-cm",          0, &config_cmfilter,    1},
		{"filter-noise",       0, &config_noisefilter, 1},
455
		{"verbose",            0, &config_verbose,     1},
Thomas White's avatar
Thomas White committed
456
		{"pdb",                1, NULL,               'p'},
457
		{"prefix",             1, NULL,               'x'},
Thomas White's avatar
Thomas White committed
458
		{"unpolarized",        0, &config_polar,       0},
459
460
		{"no-sat-corr",        0, &config_satcorr,     0},
		{"sat-corr",           0, &config_satcorr,     1}, /* Compat */
461
		{"threshold",          1, NULL,               't'},
462
		{"min-gradient",       1, NULL,                4},
463
		{"no-check-prefix",    0, &config_checkprefix, 0},
464
		{"no-closer-peak",     0, &config_closer,      0},
465
		{"insane",             0, &config_insane,      1},
Thomas White's avatar
Thomas White committed
466
		{"image",              1, NULL,               'e'},
467
		{"basename",           0, &config_basename,    1},
Thomas White's avatar
Thomas White committed
468
		{"record",             1, NULL,                5},
Thomas White's avatar
Thomas White committed
469
470
471
472
		{0, 0, NULL, 0}
	};

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

		switch (c) {
Thomas White's avatar
Thomas White committed
477
		case 'h' :
Thomas White's avatar
Thomas White committed
478
479
480
			show_help(argv[0]);
			return 0;

Thomas White's avatar
Thomas White committed
481
		case 'i' :
Thomas White's avatar
Thomas White committed
482
483
484
			filename = strdup(optarg);
			break;

485
486
487
488
		case 'o' :
			outfile = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
489
		case 'z' :
Thomas White's avatar
Thomas White committed
490
491
492
			indm_str = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
493
		case 'p' :
Thomas White's avatar
Thomas White committed
494
495
496
			pdb = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
497
		case 'x' :
498
499
500
			prefix = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
501
		case 'j' :
502
503
504
			nthreads = atoi(optarg);
			break;

Thomas White's avatar
Thomas White committed
505
		case 'g' :
Thomas White's avatar
Thomas White committed
506
507
508
			geometry = strdup(optarg);
			break;

509
510
511
512
		case 't' :
			threshold = strtof(optarg, NULL);
			break;

513
514
515
516
517
518
519
520
521
		case 'b' :
			beam = get_beam_parameters(optarg);
			if ( beam == NULL ) {
				ERROR("Failed to load beam parameters"
				      " from '%s'\n", optarg);
				return 1;
			}
			break;

522
523
524
525
		case 2 :
			speaks = strdup(optarg);
			break;

526
527
528
529
		case 3 :
			scellr = strdup(optarg);
			break;

530
531
532
533
		case 4 :
			min_gradient = strtof(optarg, NULL);
			break;

Thomas White's avatar
Thomas White committed
534
535
536
537
		case 'e' :
			element = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
538
539
540
541
542
		case 5 :
			stream_flags = parse_stream_flags(optarg);
			if ( stream_flags < 0 ) return 1;
			break;

Thomas White's avatar
Thomas White committed
543
		case 0 :
Thomas White's avatar
Thomas White committed
544
545
			break;

Thomas White's avatar
Thomas White committed
546
		default :
Thomas White's avatar
Thomas White committed
547
548
549
550
551
552
553
554
555
556
557
558
559
560
			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
561
		ERROR("Failed to open input file '%s'\n", filename);
Thomas White's avatar
Thomas White committed
562
563
		return 1;
	}
Thomas White's avatar
Thomas White committed
564
	free(filename);
Thomas White's avatar
Thomas White committed
565

566
567
568
569
570
571
572
573
574
575
576
577
578
579
	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);

580
581
582
583
584
585
586
587
588
589
590
591
592
	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
593
	free(speaks);
594

Thomas White's avatar
Thomas White committed
595
596
597
598
	if ( pdb == NULL ) {
		pdb = strdup("molecule.pdb");
	}

599
	if ( prefix == NULL ) {
Thomas White's avatar
Thomas White committed
600
		prefix = strdup("");
601
	} else {
602
603
604
		if ( config_checkprefix ) {
			prefix = check_prefix(prefix);
		}
605
606
	}

607
	if ( nthreads == 0 ) {
608
609
610
611
		ERROR("Invalid number of threads.\n");
		return 1;
	}

612
613
	if ( (indm_str == NULL) ||
	     ((indm_str != NULL) && (strcmp(indm_str, "none") == 0)) ) {
614
		STATUS("Not indexing anything.\n");
615
		indexer_needs_cell = 0;
616
617
618
		reduction_needs_cell = 0;
		indm = NULL;
		cellr = CELLR_NONE;
Thomas White's avatar
Thomas White committed
619
	} else {
620
621
622
623
624
625
626
627
628
629
630
631
632
633
		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
634
		}
Thomas White's avatar
Thomas White committed
635

636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
		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. */
656
657
	}

658
659
660
	/* No indexing -> no reduction */
	if ( indm == NULL ) reduction_needs_cell = 0;

Thomas White's avatar
Thomas White committed
661
662
663
664
665
666
667
668
669
670
671
672
	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
673
	if ( reduction_needs_cell || indexer_needs_cell ) {
674
675
		cell = load_cell_from_pdb(pdb);
		if ( cell == NULL ) {
676
677
678
			ERROR("Couldn't read unit cell (from %s)\n", pdb);
			return 1;
		}
679
	} else {
Thomas White's avatar
Thomas White committed
680
681
		STATUS("No cell needed for these choices of indexing"
		       " and reduction.\n");
682
		cell = NULL;
683
	}
684
	free(pdb);
685

686
	/* Start by writing the entire command line to stdout */
Thomas White's avatar
Thomas White committed
687
	fprintf(ofh, "CrystFEL stream format 2.0\n");
688
689
690
691
	fprintf(ofh, "Command line:");
	for ( i=0; i<argc; i++ ) {
		fprintf(ofh, " %s", argv[i]);
	}
692
693
	fprintf(ofh, "\n");
	fflush(ofh);
694

695
696
697
698
699
700
701
702
	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;
	}

703
	/* Get first filename and use it to set up the indexing */
704
	prepare_line = malloc(1024*sizeof(char));
705
706
707
708
	rval = fgets(prepare_line, 1023, fh);
	if ( rval == NULL ) {
		ERROR("Failed to get filename to prepare indexing.\n");
		return 1;
709
	}
710
	chomp(prepare_line);
711
712
713
714
715
716
	if ( config_basename ) {
		char *tmp;
		tmp = safe_basename(prepare_line);
		free(prepare_line);
		prepare_line = tmp;
	}
717
718
719
720
	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
721
722
723
724
725
726
727
728
729
	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;
730
731
	}

Thomas White's avatar
Thomas White committed
732
	gsl_set_error_handler_off();
Thomas White's avatar
Thomas White committed
733

734
735
736
737
738
739
740
	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;
741
	qargs.static_args.config_insane = config_insane;
742
	qargs.static_args.cellr = cellr;
743
	qargs.static_args.threshold = threshold;
744
	qargs.static_args.min_gradient = min_gradient;
745
746
747
748
749
750
	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;
751
	qargs.static_args.beam = beam;
Thomas White's avatar
Thomas White committed
752
	qargs.static_args.element = element;
753
	qargs.static_args.stream_flags = stream_flags;
754
755
756

	qargs.fh = fh;
	qargs.prefix = prefix;
757
	qargs.config_basename = config_basename;
758
759
760
761
	qargs.n_indexable = 0;

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

763
764
	cleanup_indexing(ipriv);

Thomas White's avatar
Thomas White committed
765
766
	free(indm);
	free(ipriv);
767
	free(prefix);
Thomas White's avatar
Thomas White committed
768
769
	free(det->panels);
	free(det);
Thomas White's avatar
Thomas White committed
770
	free(element);
771
	cell_free(cell);
Thomas White's avatar
Thomas White committed
772
773
	if ( fh != stdin ) fclose(fh);
	if ( ofh != stdout ) fclose(ofh);
Thomas White's avatar
Thomas White committed
774

775
776
	STATUS("There were %i images, of which %i could be indexed.\n",
	        n_images, qargs.n_indexable);
Thomas White's avatar
Thomas White committed
777
778
779

	return 0;
}