indexamajig.c 20.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
 * (c) 2006-2010 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"
33
#include "sfac.h"
Thomas White's avatar
Thomas White committed
34
#include "filters.h"
35
#include "reflections.h"
36
#include "thread-pool.h"
37
#include "beam-parameters.h"
38
#include "geometry.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
54
55
{
	UnitCell *cell;
	int config_cmfilter;
	int config_noisefilter;
	int config_dumpfound;
	int config_verbose;
	int config_nearbragg;
Thomas White's avatar
Thomas White committed
56
	int config_polar;
57
	int config_satcorr;
58
	int config_closer;
59
	int config_insane;
60
	float threshold;
61
	float min_gradient;
Thomas White's avatar
Thomas White committed
62
	struct detector *det;
Thomas White's avatar
Thomas White committed
63
64
	IndexingMethod *indm;
	IndexingPrivate **ipriv;
65
	int peaks;
66
	int cellr;
67
	struct beam_params *beam;
Thomas White's avatar
Thomas White committed
68
	const char *element;
69

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


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

	int n_indexable;
97
98

	char *use_this_one_instead;
99
100
101
};


Thomas White's avatar
Thomas White committed
102
103
104
105
106
107
static void show_help(const char *s)
{
	printf("Syntax: %s [options]\n\n", s);
	printf(
"Process and index FEL diffraction images.\n"
"\n"
108
" -h, --help               Display this help message.\n"
Thomas White's avatar
Thomas White committed
109
"\n"
110
" -i, --input=<filename>   Specify file containing list of images to process.\n"
Thomas White's avatar
Thomas White committed
111
"                           '-' means stdin, which is the default.\n"
112
" -o, --output=<filename>  Write indexed stream to this file. '-' for stdout.\n"
Thomas White's avatar
Thomas White committed
113
"\n"
Thomas White's avatar
Thomas White committed
114
115
116
117
118
119
"     --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"
"                            template : index by template matching\n"
120
" -g. --geometry=<file>    Get detector geometry from file.\n"
121
122
123
" -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"
124
" -p, --pdb=<file>         PDB file from which to get the unit cell to match.\n"
Thomas White's avatar
Thomas White committed
125
"                           Default: 'molecule.pdb'.\n"
126
"     --basename           Remove the directory parts of the filenames.\n"
127
" -x, --prefix=<p>         Prefix filenames from input file with <p>.\n"
128
129
130
131
132
"     --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
133
134
135
"\n"
"\nWith just the above options, this program does not do much of practical use."
"\nYou should also enable some of the following:\n\n"
136
"     --near-bragg         Output a list of reflection intensities to stdout.\n"
Thomas White's avatar
Thomas White committed
137
138
139
140
141
142
143
144
145
146
"                           When pixels with fractional indices within 0.1 of\n"
"                           integer values (the Bragg condition) are found,\n"
"                           the integral of pixels within a ten pixel radius\n"
"                           of the nearest-to-Bragg pixel will be reported as\n"
"                           the intensity.  The centroid of the pixels will\n"
"                           be given as the coordinates, as well as the h,k,l\n"
"                           (integer) indices of the reflection.  If a peak\n"
"                           was located by the initial peak search close to\n"
"                           the \"near Bragg\" location, its coordinates will\n"
"                           be taken as the centre instead.\n"
147
"     --dump-peaks         Write the results of the peak search to stdout.\n"
148
149
"                           The intensities in this list are from the\n"
"                           centroid/integration procedure.\n"
Thomas White's avatar
Thomas White committed
150
151
"\n"
"\nFor more control over the process, you might need:\n\n"
Thomas White's avatar
Thomas White committed
152
153
154
155
156
"     --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"
157
"     --filter-cm          Perform common-mode noise subtraction on images\n"
Thomas White's avatar
Thomas White committed
158
159
"                           before proceeding.  Intensities will be extracted\n"
"                           from the image as it is after this processing.\n"
160
"     --filter-noise       Apply an aggressive noise filter which sets all\n"
Thomas White's avatar
Thomas White committed
161
162
163
"                           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"
164
"     --unpolarized        Don't correct for the polarisation of the X-rays.\n"
165
166
"     --no-sat-corr        Don't correct values of saturated peaks using a\n"
"                           table included in the HDF5 file.\n"
167
"     --threshold=<n>      Only accept peaks above <n> ADU.  Default: 800.\n"
168
169
"     --min-gradient=<n>   Minimum gradient for Zaefferer peak search.\n"
"                           Default: 100,000.\n"
Thomas White's avatar
Thomas White committed
170
171
172
" -e, --image=<element>    Use this image from the HDF5 file.\n"
"                           Example: /data/data0."
"                           Default: The first one found.\n"
Thomas White's avatar
Thomas White committed
173
"\n"
174
175
176
"\nOptions for greater performance or verbosity:\n\n"
"     --verbose            Be verbose about indexing.\n"
" -j <n>                   Run <n> analyses in parallel.  Default 1.\n"
177
178
179
"\n"
"\nOptions you probably won't need:\n\n"
"     --no-check-prefix    Don't attempt to correct the --prefix.\n"
180
181
182
"     --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"
183
184
"     --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
185
);
Thomas White's avatar
Thomas White committed
186
187
188
}


189
static void process_image(void *pp, int cookie)
190
{
191
	struct index_args *pargs = pp;
192
193
194
195
	struct hdfile *hdfile;
	struct image image;
	float *data_for_measurement;
	size_t data_size;
196
	char *filename = pargs->filename;
197
198
199
200
201
202
203
	UnitCell *cell = pargs->static_args.cell;
	int config_cmfilter = pargs->static_args.config_cmfilter;
	int config_noisefilter = pargs->static_args.config_noisefilter;
	int config_dumpfound = pargs->static_args.config_dumpfound;
	int config_verbose = pargs->static_args.config_verbose;
	int config_nearbragg = pargs->static_args.config_nearbragg;
	int config_polar = pargs->static_args.config_polar;
Thomas White's avatar
Thomas White committed
204
	IndexingMethod *indm = pargs->static_args.indm;
205
	const struct beam_params *beam = pargs->static_args.beam;
206
207
208
209

	image.features = NULL;
	image.data = NULL;
	image.indexed_cell = NULL;
210
	image.id = cookie;
211
	image.filename = filename;
212
	image.det = pargs->static_args.det;
213

214
	STATUS("Processing '%s'\n", image.filename);
215

Thomas White's avatar
Thomas White committed
216
	pargs->indexable = 0;
217

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

242
243
	}

244
	hdf5_read(hdfile, &image, pargs->static_args.config_satcorr,
245
	          beam->photon_energy);
246
247
248
249
250
251
252
253
254
255
256
257
258

	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 {
259
		memcpy(data_for_measurement, image.data, data_size);
260
261
	}

262
	switch ( pargs->static_args.peaks )
263
264
265
266
267
268
269
270
271
	{
	case PEAK_HDF5 :
		/* Get peaks from HDF5 */
		if ( get_peaks(&image, hdfile) ) {
			ERROR("Failed to get peaks from HDF5 file.\n");
			return;
		}
		break;
	case PEAK_ZAEF :
272
273
		search_peaks(&image, pargs->static_args.threshold,
		             pargs->static_args.min_gradient);
274
275
		break;
	}
Thomas White's avatar
Thomas White committed
276
277
278
279
280
281

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

282
	if ( config_dumpfound ) {
283
284
		dump_peaks(&image, pargs->static_args.ofh,
		           pargs->static_args.output_mutex);
285
	}
286

287
288
	/* Not indexing? Then there's nothing left to do. */
	if ( indm == NULL ) goto done;
289
290

	/* Calculate orientation matrix (by magic) */
291
	index_pattern(&image, cell, indm, pargs->static_args.cellr,
292
293
		      config_verbose, pargs->static_args.ipriv,
		      pargs->static_args.config_insane);
294
295
296

	/* No cell at this point?  Then we're done. */
	if ( image.indexed_cell == NULL ) goto done;
Thomas White's avatar
Thomas White committed
297
	pargs->indexable = 1;
298
299
300

	/* Measure intensities if requested */
	if ( config_nearbragg ) {
301
302
303

		RefList *reflections;

304
305
		image.div = beam->divergence;
		image.bw = beam->bandwidth;
Thomas White's avatar
Thomas White committed
306
		image.profile_radius = 0.0001e9;
307

Thomas White's avatar
Thomas White committed
308
309
310
311
		//reflections = find_intersections(&image, image.indexed_cell,
		//                                 0);
		reflections = find_projected_peaks(&image, image.indexed_cell,
		                                   0, 0.1);
312
313

		output_intensities(&image, image.indexed_cell, reflections,
314
		                   pargs->static_args.output_mutex,
Thomas White's avatar
Thomas White committed
315
		                   config_polar,
316
		                   pargs->static_args.config_closer,
317
318
319
320
		                   pargs->static_args.ofh);

		reflist_free(reflections);

321
322
323
	}

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

done:
	free(image.data);
Thomas White's avatar
Thomas White committed
328
	free(image.flags);
329
330
331
332
333
	image_feature_list_free(image.features);
	hdfile_close(hdfile);
}


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

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

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

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

	} else {

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

	}
363

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

377
378
	return pargs;
}
379
380


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

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

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


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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Thomas White's avatar
Thomas White committed
537
		case 0 :
Thomas White's avatar
Thomas White committed
538
539
			break;

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

560
561
562
563
564
565
566
567
568
569
570
571
572
573
	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);

574
575
576
577
578
579
580
581
582
583
584
585
586
	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
587
	free(speaks);
588

Thomas White's avatar
Thomas White committed
589
590
591
592
	if ( pdb == NULL ) {
		pdb = strdup("molecule.pdb");
	}

593
	if ( prefix == NULL ) {
Thomas White's avatar
Thomas White committed
594
		prefix = strdup("");
595
	} else {
596
597
598
		if ( config_checkprefix ) {
			prefix = check_prefix(prefix);
		}
599
600
	}

601
	if ( nthreads == 0 ) {
602
603
604
605
		ERROR("Invalid number of threads.\n");
		return 1;
	}

606
607
	if ( strcmp(indm_str, "none") == 0 ) {
		STATUS("Not indexing anything.\n");
608
		indexer_needs_cell = 0;
609
610
611
		reduction_needs_cell = 0;
		indm = NULL;
		cellr = CELLR_NONE;
Thomas White's avatar
Thomas White committed
612
	} else {
613
614
615
616
617
618
619
620
621
622
623
624
625
626
		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
627
		}
Thomas White's avatar
Thomas White committed
628

629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
		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. */
649
650
	}

651
652
653
	/* No indexing -> no reduction */
	if ( indm == NULL ) reduction_needs_cell = 0;

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

679
680
681
682
683
	/* Start by writing the entire command line to stdout */
	fprintf(ofh, "Command line:");
	for ( i=0; i<argc; i++ ) {
		fprintf(ofh, " %s", argv[i]);
	}
684
685
	fprintf(ofh, "\n");
	fflush(ofh);
686

687
688
689
690
691
692
693
694
	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;
	}

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

Thomas White's avatar
Thomas White committed
724
	gsl_set_error_handler_off();
Thomas White's avatar
Thomas White committed
725

726
727
728
729
730
731
732
733
734
	qargs.static_args.cell = cell;
	qargs.static_args.config_cmfilter = config_cmfilter;
	qargs.static_args.config_noisefilter = config_noisefilter;
	qargs.static_args.config_dumpfound = config_dumpfound;
	qargs.static_args.config_verbose = config_verbose;
	qargs.static_args.config_nearbragg = config_nearbragg;
	qargs.static_args.config_polar = config_polar;
	qargs.static_args.config_satcorr = config_satcorr;
	qargs.static_args.config_closer = config_closer;
735
	qargs.static_args.config_insane = config_insane;
736
	qargs.static_args.cellr = cellr;
737
	qargs.static_args.threshold = threshold;
738
	qargs.static_args.min_gradient = min_gradient;
739
740
741
742
743
744
	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;
745
	qargs.static_args.beam = beam;
Thomas White's avatar
Thomas White committed
746
	qargs.static_args.element = element;
747
748
749

	qargs.fh = fh;
	qargs.prefix = prefix;
750
	qargs.config_basename = config_basename;
751
752
753
754
	qargs.n_indexable = 0;

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

756
757
	cleanup_indexing(ipriv);

Thomas White's avatar
Thomas White committed
758
759
	free(indm);
	free(ipriv);
760
	free(prefix);
Thomas White's avatar
Thomas White committed
761
762
	free(det->panels);
	free(det);
Thomas White's avatar
Thomas White committed
763
	free(element);
764
	cell_free(cell);
Thomas White's avatar
Thomas White committed
765
766
	if ( fh != stdin ) fclose(fh);
	if ( ofh != stdout ) fclose(ofh);
Thomas White's avatar
Thomas White committed
767

768
769
	STATUS("There were %i images, of which %i could be indexed.\n",
	        n_images, qargs.n_indexable);
Thomas White's avatar
Thomas White committed
770
771
772

	return 0;
}