indexamajig.c 24.3 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 "diffraction.h"
Thomas White's avatar
Thomas White committed
33
#include "diffraction-gpu.h"
34
#include "detector.h"
35
#include "sfac.h"
Thomas White's avatar
Thomas White committed
36
#include "filters.h"
37
#include "reflections.h"
38
#include "thread-pool.h"
39
#include "beam-parameters.h"
40

41
42
43
44
45
46
47

enum {
	PEAK_ZAEF,
	PEAK_HDF5,
};


48
49
/* Information about the indexing process which is common to all patterns */
struct static_index_args
50
{
51
	pthread_mutex_t *gpu_mutex;     /* Protects "gctx" */
52
53
54
55
56
57
58
59
60
61
	UnitCell *cell;
	int config_cmfilter;
	int config_noisefilter;
	int config_writedrx;
	int config_dumpfound;
	int config_verbose;
	int config_alternate;
	int config_nearbragg;
	int config_gpu;
	int config_simulate;
Thomas White's avatar
Thomas White committed
62
	int config_polar;
63
	int config_sanity;
64
	int config_satcorr;
65
	int config_sa;
66
	int config_closer;
67
	float threshold;
Thomas White's avatar
Thomas White committed
68
	struct detector *det;
69
	IndexingMethod indm;
70
	IndexingPrivate *ipriv;
71
72
	const double *intensities;
	struct gpu_context *gctx;
73
	int peaks;
74
	int cellr;
75
	double nominal_photon_energy;
76

77
78
79
	/* Output stream */
	pthread_mutex_t *output_mutex;  /* Protects the output stream */
	FILE *ofh;
80
81
82
};


83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
/* Information about the indexing process for one pattern */
struct index_args
{
	/* "Input" */
	char *filename;
	struct static_index_args static_args;

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


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

	int n_indexable;
	int n_sane;
105
106

	char *use_this_one_instead;
107
108
109
};


Thomas White's avatar
Thomas White committed
110
111
112
113
114
115
static void show_help(const char *s)
{
	printf("Syntax: %s [options]\n\n", s);
	printf(
"Process and index FEL diffraction images.\n"
"\n"
116
" -h, --help               Display this help message.\n"
Thomas White's avatar
Thomas White committed
117
"\n"
118
" -i, --input=<filename>   Specify file containing list of images to process.\n"
Thomas White's avatar
Thomas White committed
119
"                           '-' means stdin, which is the default.\n"
120
" -o, --output=<filename>  Write indexed stream to this file. '-' for stdout.\n"
Thomas White's avatar
Thomas White committed
121
"\n"
122
"     --indexing=<method>  Use 'method' for indexing.  Choose from:\n"
Thomas White's avatar
Thomas White committed
123
"                           none     : no indexing (default)\n"
Thomas White's avatar
Thomas White committed
124
"                           dirax    : invoke DirAx\n"
125
"                           template : index by template matching\n"
126
" -g. --geometry=<file>    Get detector geometry from file.\n"
127
128
129
" -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"
130
" -p, --pdb=<file>         PDB file from which to get the unit cell to match.\n"
Thomas White's avatar
Thomas White committed
131
"                           Default: 'molecule.pdb'.\n"
132
" -x, --prefix=<p>         Prefix filenames from input file with <p>.\n"
133
134
135
136
137
"     --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
138
139
140
"\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"
141
"     --near-bragg         Output a list of reflection intensities to stdout.\n"
Thomas White's avatar
Thomas White committed
142
143
144
145
146
147
148
149
150
151
"                           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"
152
"     --simulate           Simulate the diffraction pattern using the indexed\n"
Thomas White's avatar
Thomas White committed
153
154
155
156
"                           unit cell.  The simulated pattern will be saved\n"
"                           as \"simulated.h5\".  You can TRY to combine this\n"
"                           with \"-j <n>\" with n greater than 1, but it's\n"
"                           not a good idea.\n"
157
"     --write-drx          Write 'xfel.drx' for visualisation of reciprocal\n"
158
159
160
"                           space.  Implied by any indexing method other than\n"
"                           'none'.  Beware: the units in this file are\n"
"                           reciprocal Angstroms.\n"
161
"     --dump-peaks         Write the results of the peak search to stdout.\n"
162
163
"                           The intensities in this list are from the\n"
"                           centroid/integration procedure.\n"
Thomas White's avatar
Thomas White committed
164
165
"\n"
"\nFor more control over the process, you might need:\n\n"
Thomas White's avatar
Thomas White committed
166
167
168
169
170
"     --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"
171
"     --check-sanity       Check that indexed locations approximately correspond\n"
172
"                           with detected peaks.\n"
173
"     --filter-cm          Perform common-mode noise subtraction on images\n"
Thomas White's avatar
Thomas White committed
174
175
"                           before proceeding.  Intensities will be extracted\n"
"                           from the image as it is after this processing.\n"
176
"     --filter-noise       Apply an aggressive noise filter which sets all\n"
Thomas White's avatar
Thomas White committed
177
178
179
"                           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"
180
"     --unpolarized        Don't correct for the polarisation of the X-rays.\n"
181
182
"     --no-sat-corr        Don't correct values of saturated peaks using a\n"
"                           table included in the HDF5 file.\n"
183
"     --no-sa              Don't correct for the differing solid angles of\n"
184
"                           the pixels.\n"
185
"     --threshold=<n>      Only accept peaks above <n> ADU.  Default: 800.\n"
Thomas White's avatar
Thomas White committed
186
"\n"
187
"\nIf you used --simulate, you may also want:\n\n"
Thomas White's avatar
Thomas White committed
188
189
"     --intensities=<file> Specify file containing reflection intensities\n"
"                           to use when simulating.\n"
190
191
192
193
194
"\n"
"\nOptions for greater performance or verbosity:\n\n"
"     --verbose            Be verbose about indexing.\n"
"     --gpu                Use the GPU to speed up the simulation.\n"
" -j <n>                   Run <n> analyses in parallel.  Default 1.\n"
195
196
197
"\n"
"\nOptions you probably won't need:\n\n"
"     --no-check-prefix    Don't attempt to correct the --prefix.\n"
198
199
200
"     --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"
Thomas White's avatar
Thomas White committed
201
);
Thomas White's avatar
Thomas White committed
202
203
204
}


205
static struct image *get_simage(struct image *template, int alternate)
206
{
207
	struct image *image;
208
	struct panel panels[2];
209

210
211
	image = malloc(sizeof(*image));

212
	/* Simulate a diffraction pattern */
213
214
215
	image->twotheta = NULL;
	image->data = NULL;
	image->det = template->det;
216
217
218
	image->flags = NULL;
	image->f0_available = 0;
	image->f0 = 1.0;
219
220
221

	/* Detector geometry for the simulation
	 * - not necessarily the same as the original. */
222
223
	image->width = 1024;
	image->height = 1024;
Thomas White's avatar
Thomas White committed
224
	image->det->n_panels = 2;
225

226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
	if ( alternate ) {

		/* Upper */
		panels[0].min_x = 0;
		panels[0].max_x = 1023;
		panels[0].min_y = 512;
		panels[0].max_y = 1023;
		panels[0].cx = 523.6;
		panels[0].cy = 502.5;
		panels[0].clen = 56.4e-2;  /* 56.4 cm */
		panels[0].res = 13333.3;   /* 75 microns/pixel */

		/* Lower */
		panels[1].min_x = 0;
		panels[1].max_x = 1023;
		panels[1].min_y = 0;
		panels[1].max_y = 511;
		panels[1].cx = 520.8;
244
		panels[1].cy = 525.0;
245
246
247
		panels[1].clen = 56.7e-2;  /* 56.7 cm */
		panels[1].res = 13333.3;   /* 75 microns/pixel */

Thomas White's avatar
Thomas White committed
248
		image->det->panels = panels;
249

250
251
252
	} else {

		/* Copy pointer to old geometry */
Thomas White's avatar
Thomas White committed
253
		image->det->panels = template->det->panels;
254
255

	}
256

257
	image->lambda = ph_en_to_lambda(eV_to_J(1.8e3));
258
	image->features = template->features;
259
	image->filename = template->filename;
260
	image->indexed_cell = template->indexed_cell;
Thomas White's avatar
Thomas White committed
261
	image->f0 = template->f0;
262

Thomas White's avatar
Thomas White committed
263
	/* Prevent muppetry */
Thomas White's avatar
Thomas White committed
264
265
	image->cpeaks = NULL;
	image->n_cpeaks = 0;
Thomas White's avatar
Thomas White committed
266

267
268
	return image;
}
269
270


271
static void simulate_and_write(struct image *simage, struct gpu_context **gctx,
Thomas White's avatar
Thomas White committed
272
                               const double *intensities, UnitCell *cell)
273
{
274
	/* Set up GPU if necessary */
Thomas White's avatar
Thomas White committed
275
	if ( (gctx != NULL) && (*gctx == NULL) ) {
Thomas White's avatar
Thomas White committed
276
		*gctx = setup_gpu(0, simage, intensities);
277
278
	}

Thomas White's avatar
Thomas White committed
279
	if ( (gctx != NULL) && (*gctx != NULL) ) {
280
		get_diffraction_gpu(*gctx, simage, 24, 24, 40, cell);
281
	} else {
282
		get_diffraction(simage, 24, 24, 40,
Thomas White's avatar
Thomas White committed
283
		                intensities, NULL, cell, 0,
284
		                GRADIENT_MOSAIC);
285
	}
286
	record_image(simage, 0);
287

288
	hdf5_write("simulated.h5", simage->data, simage->width, simage->height,
289
290
291
292
		   H5T_NATIVE_FLOAT);
}


293
static void process_image(void *pp, int cookie)
294
{
295
	struct index_args *pargs = pp;
296
297
298
299
300
	struct hdfile *hdfile;
	struct image image;
	struct image *simage;
	float *data_for_measurement;
	size_t data_size;
301
	char *filename = pargs->filename;
302
303
304
305
306
307
308
309
310
311
312
313
314
315
	UnitCell *cell = pargs->static_args.cell;
	int config_cmfilter = pargs->static_args.config_cmfilter;
	int config_noisefilter = pargs->static_args.config_noisefilter;
	int config_writedrx = pargs->static_args.config_writedrx;
	int config_dumpfound = pargs->static_args.config_dumpfound;
	int config_verbose = pargs->static_args.config_verbose;
	int config_alternate  = pargs->static_args.config_alternate;
	int config_nearbragg = pargs->static_args.config_nearbragg;
	int config_gpu = pargs->static_args.config_gpu;
	int config_simulate = pargs->static_args.config_simulate;
	int config_polar = pargs->static_args.config_polar;
	IndexingMethod indm = pargs->static_args.indm;
	const double *intensities = pargs->static_args.intensities;
	struct gpu_context *gctx = pargs->static_args.gctx;
316
317
318
319

	image.features = NULL;
	image.data = NULL;
	image.indexed_cell = NULL;
320
	image.id = cookie;
321
	image.filename = filename;
Thomas White's avatar
Thomas White committed
322
323
	image.cpeaks = NULL;
	image.n_cpeaks = 0;
324
	image.det = pargs->static_args.det;
325

326
	STATUS("Processing '%s'\n", image.filename);
327

Thomas White's avatar
Thomas White committed
328
329
	pargs->sane = 0;
	pargs->indexable = 0;
330

331
332
	hdfile = hdfile_open(filename);
	if ( hdfile == NULL ) {
Thomas White's avatar
Thomas White committed
333
		return;
334
335
	} else if ( hdfile_set_first_image(hdfile, "/") ) {
		ERROR("Couldn't select path\n");
Thomas White's avatar
Thomas White committed
336
		return;
337
338
	}

339
340
	hdf5_read(hdfile, &image, pargs->static_args.config_satcorr,
	          pargs->static_args.nominal_photon_energy);
341
342
343
344
345
346
347
348
349
350
351
352
353

	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 {
354
		memcpy(data_for_measurement, image.data, data_size);
355
356
	}

357
	switch ( pargs->static_args.peaks )
358
359
360
361
362
363
364
365
366
	{
	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 :
367
		search_peaks(&image, pargs->static_args.threshold);
368
369
		break;
	}
Thomas White's avatar
Thomas White committed
370
371
372
373
374
375

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

376
	if ( config_dumpfound ) {
377
378
		dump_peaks(&image, pargs->static_args.ofh,
		           pargs->static_args.output_mutex);
379
	}
380
381
382
383
384
385
386
387
388

	/* Not indexing nor writing xfel.drx?
	 * Then there's nothing left to do. */
	if ( (!config_writedrx) && (indm == INDEXING_NONE) ) {
		goto done;
	}

	/* Calculate orientation matrix (by magic) */
	if ( config_writedrx || (indm != INDEXING_NONE) ) {
389
		index_pattern(&image, cell, indm, pargs->static_args.cellr,
390
		              config_verbose, pargs->static_args.ipriv);
391
392
393
394
	}

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

Thomas White's avatar
Thomas White committed
397
	/* Sanity check */
398
	if ( pargs->static_args.config_sanity
399
	  && !peak_sanity_check(&image, image.indexed_cell, 0, 0.1) ) {
Thomas White's avatar
Thomas White committed
400
401
402
		STATUS("Failed peak sanity check.\n");
		goto done;
	} else {
Thomas White's avatar
Thomas White committed
403
		pargs->sane = 1;
Thomas White's avatar
Thomas White committed
404
405
	}

406
407
	/* Measure intensities if requested */
	if ( config_nearbragg ) {
408
		output_intensities(&image, image.indexed_cell,
409
410
411
412
		                   pargs->static_args.output_mutex,
		                   config_polar, pargs->static_args.config_sa,
		                   pargs->static_args.config_closer,
		                   pargs->static_args.ofh, 0, 0.1);
413
414
	}

415
416
	simage = get_simage(&image, config_alternate);

417
418
419
	/* Simulate if requested */
	if ( config_simulate ) {
		if ( config_gpu ) {
420
			pthread_mutex_lock(pargs->static_args.gpu_mutex);
421
			simulate_and_write(simage, &gctx, intensities,
Thomas White's avatar
Thomas White committed
422
			                   image.indexed_cell);
423
			pthread_mutex_unlock(pargs->static_args.gpu_mutex);
424
425
		} else {
			simulate_and_write(simage, NULL, intensities,
Thomas White's avatar
Thomas White committed
426
			                   image.indexed_cell);
427
428
429
430
431
432
433
434
435
		}
	}

	/* Finished with alternate image */
	if ( simage->twotheta != NULL ) free(simage->twotheta);
	if ( simage->data != NULL ) free(simage->data);
	free(simage);

	/* Only free cell if found */
436
	cell_free(image.indexed_cell);
437
438
439

done:
	free(image.data);
Thomas White's avatar
Thomas White committed
440
	free(image.flags);
441
	image_feature_list_free(image.features);
Thomas White's avatar
Thomas White committed
442
	free(image.cpeaks);
443
444
445
446
	hdfile_close(hdfile);
}


447
static void *get_image(void *qp)
448
{
449
450
451
452
	char line[1024];
	struct index_args *pargs;
	char *rval;
	struct queue_args *qargs = qp;
453

454
	/* Initialise new task arguments */
455
456
457
	pargs = malloc(sizeof(struct index_args));
	memcpy(&pargs->static_args, &qargs->static_args,
	       sizeof(struct static_index_args));
458

459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
	/* Get the next filename */
	if ( qargs->use_this_one_instead != NULL ) {

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

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

		qargs->use_this_one_instead = NULL;

	} else {

		rval = fgets(line, 1023, qargs->fh);
		if ( rval == NULL ) return NULL;
		chomp(line);
		pargs->filename = malloc(strlen(qargs->prefix)+strlen(line)+1);
		snprintf(pargs->filename, 1023, "%s%s", qargs->prefix, line);

	}
479

480
481
	return pargs;
}
482
483


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

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

492
493
	free(pargs->filename);
	free(pargs);
494
495
496
}


Thomas White's avatar
Thomas White committed
497
498
499
int main(int argc, char *argv[])
{
	int c;
Thomas White's avatar
Thomas White committed
500
	struct gpu_context *gctx = NULL;
Thomas White's avatar
Thomas White committed
501
	char *filename = NULL;
502
	char *outfile = NULL;
Thomas White's avatar
Thomas White committed
503
	FILE *fh;
504
	FILE *ofh;
505
	char *rval = NULL;
Thomas White's avatar
Thomas White committed
506
	int n_images;
507
	int config_noindex = 0;
Thomas White's avatar
Thomas White committed
508
	int config_dumpfound = 0;
509
	int config_nearbragg = 0;
Thomas White's avatar
Thomas White committed
510
	int config_writedrx = 0;
511
	int config_simulate = 0;
512
513
	int config_cmfilter = 0;
	int config_noisefilter = 0;
Thomas White's avatar
Thomas White committed
514
	int config_gpu = 0;
515
	int config_verbose = 0;
516
	int config_alternate = 0;
Thomas White's avatar
Thomas White committed
517
	int config_polar = 1;
518
	int config_sanity = 0;
519
	int config_satcorr = 1;
520
	int config_sa = 1;
521
	int config_checkprefix = 1;
522
	int config_closer = 1;
523
	float threshold = 800.0;
Thomas White's avatar
Thomas White committed
524
525
	struct detector *det;
	char *geometry = NULL;
Thomas White's avatar
Thomas White committed
526
527
	IndexingMethod indm;
	char *indm_str = NULL;
528
529
530
	UnitCell *cell;
	double *intensities = NULL;
	char *intfile = NULL;
Thomas White's avatar
Thomas White committed
531
	char *pdb = NULL;
532
	char *prefix = NULL;
533
	char *speaks = NULL;
534
535
	char *scellr = NULL;
	int cellr;
536
	int peaks;
537
538
	int nthreads = 1;
	int i;
539
	pthread_mutex_t output_mutex = PTHREAD_MUTEX_INITIALIZER;
540
	pthread_mutex_t gpu_mutex = PTHREAD_MUTEX_INITIALIZER;
541
542
543
	char prepare_line[1024];
	char prepare_filename[1024];
	IndexingPrivate *ipriv;
544
	struct queue_args qargs;
545
546
	struct beam_params *beam = NULL;
	double nominal_photon_energy;
Thomas White's avatar
Thomas White committed
547
548
549
550
551

	/* Long options */
	const struct option longopts[] = {
		{"help",               0, NULL,               'h'},
		{"input",              1, NULL,               'i'},
552
		{"output",             1, NULL,               'o'},
Thomas White's avatar
Thomas White committed
553
		{"gpu",                0, &config_gpu,         1},
554
		{"no-index",           0, &config_noindex,     1},
Thomas White's avatar
Thomas White committed
555
		{"dump-peaks",         0, &config_dumpfound,   1},
556
		{"peaks",              1, NULL,                2},
557
		{"cell-reduction",     1, NULL,                3},
558
		{"near-bragg",         0, &config_nearbragg,   1},
Thomas White's avatar
Thomas White committed
559
560
		{"write-drx",          0, &config_writedrx,    1},
		{"indexing",           1, NULL,               'z'},
Thomas White's avatar
Thomas White committed
561
		{"geometry",           1, NULL,               'g'},
562
		{"beam",               1, NULL,               'b'},
563
		{"simulate",           0, &config_simulate,    1},
564
565
		{"filter-cm",          0, &config_cmfilter,    1},
		{"filter-noise",       0, &config_noisefilter, 1},
566
		{"verbose",            0, &config_verbose,     1},
567
		{"alternate",          0, &config_alternate,   1},
568
		{"intensities",        1, NULL,               'q'},
Thomas White's avatar
Thomas White committed
569
		{"pdb",                1, NULL,               'p'},
570
		{"prefix",             1, NULL,               'x'},
Thomas White's avatar
Thomas White committed
571
		{"unpolarized",        0, &config_polar,       0},
572
		{"check-sanity",       0, &config_sanity,      1},
573
574
		{"no-sat-corr",        0, &config_satcorr,     0},
		{"sat-corr",           0, &config_satcorr,     1}, /* Compat */
575
		{"no-sa",              0, &config_sa,          0},
576
		{"threshold",          1, NULL,               't'},
577
		{"no-check-prefix",    0, &config_checkprefix, 0},
578
		{"no-closer-peak",     0, &config_closer,      0},
Thomas White's avatar
Thomas White committed
579
580
581
582
		{0, 0, NULL, 0}
	};

	/* Short options */
583
	while ((c = getopt_long(argc, argv, "hi:wp:j:x:g:t:o:b:",
Thomas White's avatar
Thomas White committed
584
	                        longopts, NULL)) != -1) {
Thomas White's avatar
Thomas White committed
585
586

		switch (c) {
Thomas White's avatar
Thomas White committed
587
		case 'h' :
Thomas White's avatar
Thomas White committed
588
589
590
			show_help(argv[0]);
			return 0;

Thomas White's avatar
Thomas White committed
591
		case 'i' :
Thomas White's avatar
Thomas White committed
592
593
594
			filename = strdup(optarg);
			break;

595
596
597
598
		case 'o' :
			outfile = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
599
		case 'z' :
Thomas White's avatar
Thomas White committed
600
601
602
			indm_str = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
603
		case 'q' :
604
605
606
			intfile = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
607
		case 'p' :
Thomas White's avatar
Thomas White committed
608
609
610
			pdb = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
611
		case 'x' :
612
613
614
			prefix = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
615
		case 'j' :
616
617
618
			nthreads = atoi(optarg);
			break;

Thomas White's avatar
Thomas White committed
619
		case 'g' :
Thomas White's avatar
Thomas White committed
620
621
622
			geometry = strdup(optarg);
			break;

623
624
625
626
		case 't' :
			threshold = strtof(optarg, NULL);
			break;

627
628
629
630
631
632
633
634
635
		case 'b' :
			beam = get_beam_parameters(optarg);
			if ( beam == NULL ) {
				ERROR("Failed to load beam parameters"
				      " from '%s'\n", optarg);
				return 1;
			}
			break;

636
637
638
639
		case 2 :
			speaks = strdup(optarg);
			break;

640
641
642
643
		case 3 :
			scellr = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
644
		case 0 :
Thomas White's avatar
Thomas White committed
645
646
			break;

Thomas White's avatar
Thomas White committed
647
		default :
Thomas White's avatar
Thomas White committed
648
649
650
651
652
653
654
655
656
657
658
659
660
661
			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
662
		ERROR("Failed to open input file '%s'\n", filename);
Thomas White's avatar
Thomas White committed
663
664
		return 1;
	}
Thomas White's avatar
Thomas White committed
665
	free(filename);
Thomas White's avatar
Thomas White committed
666

667
668
669
670
671
672
673
674
675
676
677
678
679
680
	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);

681
682
683
684
685
686
687
688
689
690
691
692
693
	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
694
	free(speaks);
695

696
	if ( intfile != NULL ) {
Thomas White's avatar
Thomas White committed
697
		ReflItemList *items;
698
699
		items = read_reflections(intfile, intensities,
		                         NULL, NULL, NULL);
Thomas White's avatar
Thomas White committed
700
		delete_items(items);
701
702
703
704
	} else {
		intensities = NULL;
	}

Thomas White's avatar
Thomas White committed
705
706
707
708
	if ( pdb == NULL ) {
		pdb = strdup("molecule.pdb");
	}

709
	if ( prefix == NULL ) {
Thomas White's avatar
Thomas White committed
710
		prefix = strdup("");
711
	} else {
712
713
714
		if ( config_checkprefix ) {
			prefix = check_prefix(prefix);
		}
715
716
	}

717
	if ( nthreads == 0 ) {
718
719
720
721
		ERROR("Invalid number of threads.\n");
		return 1;
	}

Thomas White's avatar
Thomas White committed
722
723
	if ( indm_str == NULL ) {
		STATUS("You didn't specify an indexing method, so I won't"
724
725
726
		       " try to index anything.\n"
		       "If that isn't what you wanted, re-run with"
		       " --indexing=<method>.\n");
Thomas White's avatar
Thomas White committed
727
728
729
730
731
		indm = INDEXING_NONE;
	} else if ( strcmp(indm_str, "none") == 0 ) {
		indm = INDEXING_NONE;
	} else if ( strcmp(indm_str, "dirax") == 0) {
		indm = INDEXING_DIRAX;
732
733
	} else if ( strcmp(indm_str, "template") == 0) {
		indm = INDEXING_TEMPLATE;
Thomas White's avatar
Thomas White committed
734
735
736
737
	} else {
		ERROR("Unrecognised indexing method '%s'\n", indm_str);
		return 1;
	}
738
	free(indm_str);
Thomas White's avatar
Thomas White committed
739

740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
	if ( scellr == NULL ) {
		STATUS("You didn't specify a cell reduction method, so I'm"
		       " going to use 'reduce'.\n");
		cellr = CELLR_REDUCE;
	} else if ( strcmp(scellr, "none") == 0 ) {
		cellr = CELLR_NONE;
	} else if ( strcmp(scellr, "reduce") == 0) {
		cellr = CELLR_REDUCE;
	} else if ( strcmp(scellr, "compare") == 0) {
		cellr = CELLR_COMPARE;
	} else {
		ERROR("Unrecognised cell reduction method '%s'\n", scellr);
		return 1;
	}
	free(scellr);  /* free(NULL) is OK. */

Thomas White's avatar
Thomas White committed
756
757
758
759
760
761
762
763
764
765
766
767
	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);

768
	if ( (cellr != CELLR_NONE) || (indm == INDEXING_TEMPLATE) ) {
769
770
		cell = load_cell_from_pdb(pdb);
		if ( cell == NULL ) {
771
772
773
			ERROR("Couldn't read unit cell (from %s)\n", pdb);
			return 1;
		}
774
775
776
	} else {
		STATUS("No cell needed because --no-match was used.\n");
		cell = NULL;
777
	}
778
	free(pdb);
779

780
781
782
783
784
	/* Start by writing the entire command line to stdout */
	fprintf(ofh, "Command line:");
	for ( i=0; i<argc; i++ ) {
		fprintf(ofh, " %s", argv[i]);
	}
785
786
	fprintf(ofh, "\n");
	fflush(ofh);
787

788
789
790
791
792
793
794
795
	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;
	}

796
	/* Get first filename and use it to set up the indexing */
797
798
799
800
	rval = fgets(prepare_line, 1023, fh);
	if ( rval == NULL ) {
		ERROR("Failed to get filename to prepare indexing.\n");
		return 1;
801
	}
802
803
804
805
806
	chomp(prepare_line);
	snprintf(prepare_filename, 1023, "%s%s", prefix, prepare_line);
	qargs.use_this_one_instead = prepare_line;

	/* Prepare the indexer */
807
808
	ipriv = prepare_indexing(indm, cell, prepare_filename, det,
	                         nominal_photon_energy);
809
810
811
812
813
	if ( ipriv == NULL ) {
		ERROR("Failed to prepare indexing.\n");
		return 1;
	}

Thomas White's avatar
Thomas White committed
814
	gsl_set_error_handler_off();
Thomas White's avatar
Thomas White committed
815

816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
	qargs.static_args.gpu_mutex = &gpu_mutex;
	qargs.static_args.cell = cell;
	qargs.static_args.config_cmfilter = config_cmfilter;
	qargs.static_args.config_noisefilter = config_noisefilter;
	qargs.static_args.config_writedrx = config_writedrx;
	qargs.static_args.config_dumpfound = config_dumpfound;
	qargs.static_args.config_verbose = config_verbose;
	qargs.static_args.config_alternate = config_alternate;
	qargs.static_args.config_nearbragg = config_nearbragg;
	qargs.static_args.config_gpu = config_gpu;
	qargs.static_args.config_simulate = config_simulate;
	qargs.static_args.config_polar = config_polar;
	qargs.static_args.config_sanity = config_sanity;
	qargs.static_args.config_satcorr = config_satcorr;
	qargs.static_args.config_sa = config_sa;
	qargs.static_args.config_closer = config_closer;
832
	qargs.static_args.cellr = cellr;
833
834
835
836
837
838
839
840
841
	qargs.static_args.threshold = threshold;
	qargs.static_args.det = det;
	qargs.static_args.indm = indm;
	qargs.static_args.ipriv = ipriv;
	qargs.static_args.intensities = intensities;
	qargs.static_args.gctx = gctx;
	qargs.static_args.peaks = peaks;
	qargs.static_args.output_mutex = &output_mutex;
	qargs.static_args.ofh = ofh;
842
	qargs.static_args.nominal_photon_energy = nominal_photon_energy;
843
844
845
846
847
848
849
850

	qargs.fh = fh;
	qargs.prefix = prefix;
	qargs.n_indexable = 0;
	qargs.n_sane = 0;

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

852
853
	cleanup_indexing(ipriv);

854
	free(prefix);
Thomas White's avatar
Thomas White committed
855
856
	free(det->panels);
	free(det);
857
	cell_free(cell);
858
	if ( fh != stdout ) fclose(fh);
Thomas White's avatar
Thomas White committed
859

Thomas White's avatar
Thomas White committed
860
	STATUS("There were %i images.  %i could be indexed, of which %i"
861
	       " looked sane.\n", n_images, qargs.n_indexable, qargs.n_sane);
Thomas White's avatar
Thomas White committed
862

Thomas White's avatar
Thomas White committed
863
864
865
866
	if ( gctx != NULL ) {
		cleanup_gpu(gctx);
	}

Thomas White's avatar
Thomas White committed
867
868
	return 0;
}