indexamajig.c 24 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
62
	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;
	int config_nomatch;
Thomas White's avatar
Thomas White committed
63
	int config_polar;
64
	int config_sanity;
65
	int config_satcorr;
66
	int config_sa;
67
	int config_closer;
68
	float threshold;
Thomas White's avatar
Thomas White committed
69
	struct detector *det;
70
	IndexingMethod indm;
71
	IndexingPrivate *ipriv;
72
73
	const double *intensities;
	struct gpu_context *gctx;
74
	int peaks;
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"
166
"     --no-match           Don't attempt to match the indexed cell to the\n"
167
168
"                           model, just proceed with the one generated by the\n"
"                           auto-indexing procedure.\n"
169
"     --check-sanity       Check that indexed locations approximately correspond\n"
170
"                           with detected peaks.\n"
171
"     --filter-cm          Perform common-mode noise subtraction on images\n"
Thomas White's avatar
Thomas White committed
172
173
"                           before proceeding.  Intensities will be extracted\n"
"                           from the image as it is after this processing.\n"
174
"     --filter-noise       Apply an aggressive noise filter which sets all\n"
Thomas White's avatar
Thomas White committed
175
176
177
"                           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"
178
"     --unpolarized        Don't correct for the polarisation of the X-rays.\n"
179
180
"     --no-sat-corr        Don't correct values of saturated peaks using a\n"
"                           table included in the HDF5 file.\n"
181
"     --no-sa              Don't correct for the differing solid angles of\n"
182
"                           the pixels.\n"
183
"     --threshold=<n>      Only accept peaks above <n> ADU.  Default: 800.\n"
Thomas White's avatar
Thomas White committed
184
"\n"
185
"\nIf you used --simulate, you may also want:\n\n"
Thomas White's avatar
Thomas White committed
186
187
"     --intensities=<file> Specify file containing reflection intensities\n"
"                           to use when simulating.\n"
188
189
190
191
192
"\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"
193
194
195
"\n"
"\nOptions you probably won't need:\n\n"
"     --no-check-prefix    Don't attempt to correct the --prefix.\n"
196
197
198
"     --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
199
);
Thomas White's avatar
Thomas White committed
200
201
202
}


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

208
209
	image = malloc(sizeof(*image));

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

	/* View head-on (unit cell is tilted) */
219
220
221
222
	image->orientation.w = 1.0;
	image->orientation.x = 0.0;
	image->orientation.y = 0.0;
	image->orientation.z = 0.0;
223
224
225

	/* Detector geometry for the simulation
	 * - not necessarily the same as the original. */
226
227
	image->width = 1024;
	image->height = 1024;
Thomas White's avatar
Thomas White committed
228
	image->det->n_panels = 2;
229

230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
	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;
248
		panels[1].cy = 525.0;
249
250
251
		panels[1].clen = 56.7e-2;  /* 56.7 cm */
		panels[1].res = 13333.3;   /* 75 microns/pixel */

Thomas White's avatar
Thomas White committed
252
		image->det->panels = panels;
253

254
255
256
	} else {

		/* Copy pointer to old geometry */
Thomas White's avatar
Thomas White committed
257
		image->det->panels = template->det->panels;
258
259

	}
260

261
	image->lambda = ph_en_to_lambda(eV_to_J(1.8e3));
262
	image->features = template->features;
263
	image->filename = template->filename;
264
	image->indexed_cell = template->indexed_cell;
Thomas White's avatar
Thomas White committed
265
	image->f0 = template->f0;
266

Thomas White's avatar
Thomas White committed
267
	/* Prevent muppetry */
Thomas White's avatar
Thomas White committed
268
269
	image->cpeaks = NULL;
	image->n_cpeaks = 0;
Thomas White's avatar
Thomas White committed
270

271
272
	return image;
}
273
274


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

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

292
	hdf5_write("simulated.h5", simage->data, simage->width, simage->height,
293
294
295
296
		   H5T_NATIVE_FLOAT);
}


297
static void process_image(void *pp, int cookie)
298
{
299
	struct index_args *pargs = pp;
300
301
302
303
304
	struct hdfile *hdfile;
	struct image image;
	struct image *simage;
	float *data_for_measurement;
	size_t data_size;
305
	char *filename = pargs->filename;
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
	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_nomatch = pargs->static_args.config_nomatch;
	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;
321
322
323
324

	image.features = NULL;
	image.data = NULL;
	image.indexed_cell = NULL;
325
	image.id = cookie;
326
	image.filename = filename;
Thomas White's avatar
Thomas White committed
327
328
	image.cpeaks = NULL;
	image.n_cpeaks = 0;
329
	image.det = pargs->static_args.det;
330

331
332
333
334
335
336
	/* View head-on (unit cell is tilted) */
	image.orientation.w = 1.0;
	image.orientation.x = 0.0;
	image.orientation.y = 0.0;
	image.orientation.z = 0.0;

337
	STATUS("Processing '%s'\n", image.filename);
338

Thomas White's avatar
Thomas White committed
339
340
	pargs->sane = 0;
	pargs->indexable = 0;
341

342
343
	hdfile = hdfile_open(filename);
	if ( hdfile == NULL ) {
Thomas White's avatar
Thomas White committed
344
		return;
345
346
	} else if ( hdfile_set_first_image(hdfile, "/") ) {
		ERROR("Couldn't select path\n");
Thomas White's avatar
Thomas White committed
347
		return;
348
349
	}

350
351
	hdf5_read(hdfile, &image, pargs->static_args.config_satcorr,
	          pargs->static_args.nominal_photon_energy);
352
353
354
355
356
357
358
359
360
361
362
363
364

	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 {
365
		memcpy(data_for_measurement, image.data, data_size);
366
367
	}

368
	switch ( pargs->static_args.peaks )
369
370
371
372
373
374
375
376
377
	{
	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 :
378
		search_peaks(&image, pargs->static_args.threshold);
379
380
		break;
	}
Thomas White's avatar
Thomas White committed
381
382
383
384
385
386

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

387
	if ( config_dumpfound ) {
388
389
		dump_peaks(&image, pargs->static_args.ofh,
		           pargs->static_args.output_mutex);
390
	}
391
392
393
394
395
396
397
398
399
400

	/* 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) ) {
		index_pattern(&image, cell, indm, config_nomatch,
401
		              config_verbose, pargs->static_args.ipriv);
402
403
404
405
	}

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

Thomas White's avatar
Thomas White committed
408
	/* Sanity check */
409
	if ( pargs->static_args.config_sanity
410
	  && !peak_sanity_check(&image, image.indexed_cell, 0, 0.1) ) {
Thomas White's avatar
Thomas White committed
411
412
413
		STATUS("Failed peak sanity check.\n");
		goto done;
	} else {
Thomas White's avatar
Thomas White committed
414
		pargs->sane = 1;
Thomas White's avatar
Thomas White committed
415
416
	}

417
418
	/* Measure intensities if requested */
	if ( config_nearbragg ) {
419
		output_intensities(&image, image.indexed_cell,
420
421
422
423
		                   pargs->static_args.output_mutex,
		                   config_polar, pargs->static_args.config_sa,
		                   pargs->static_args.config_closer,
		                   pargs->static_args.ofh, 0, 0.1);
424
425
	}

426
427
	simage = get_simage(&image, config_alternate);

428
429
430
	/* Simulate if requested */
	if ( config_simulate ) {
		if ( config_gpu ) {
431
			pthread_mutex_lock(pargs->static_args.gpu_mutex);
432
			simulate_and_write(simage, &gctx, intensities,
Thomas White's avatar
Thomas White committed
433
			                   image.indexed_cell);
434
			pthread_mutex_unlock(pargs->static_args.gpu_mutex);
435
436
		} else {
			simulate_and_write(simage, NULL, intensities,
Thomas White's avatar
Thomas White committed
437
			                   image.indexed_cell);
438
439
440
441
442
443
444
445
446
		}
	}

	/* 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 */
447
	cell_free(image.indexed_cell);
448
449
450

done:
	free(image.data);
Thomas White's avatar
Thomas White committed
451
	free(image.flags);
452
	image_feature_list_free(image.features);
Thomas White's avatar
Thomas White committed
453
	free(image.cpeaks);
454
455
456
457
	hdfile_close(hdfile);
}


458
static void *get_image(void *qp)
459
{
460
461
462
463
	char line[1024];
	struct index_args *pargs;
	char *rval;
	struct queue_args *qargs = qp;
464

465
	/* Initialise new task arguments */
466
467
468
	pargs = malloc(sizeof(struct index_args));
	memcpy(&pargs->static_args, &qargs->static_args,
	       sizeof(struct static_index_args));
469

470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
	/* 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);

	}
490

491
492
	return pargs;
}
493
494


495
496
497
498
static void finalise_image(void *qp, void *pp)
{
	struct queue_args *qargs = qp;
	struct index_args *pargs = pp;
499

500
501
	qargs->n_indexable += pargs->indexable;
	qargs->n_sane += pargs->sane;
502

503
504
	free(pargs->filename);
	free(pargs);
505
506
507
}


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

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

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

		switch (c) {
Thomas White's avatar
Thomas White committed
597
		case 'h' :
Thomas White's avatar
Thomas White committed
598
599
600
			show_help(argv[0]);
			return 0;

Thomas White's avatar
Thomas White committed
601
		case 'i' :
Thomas White's avatar
Thomas White committed
602
603
604
			filename = strdup(optarg);
			break;

605
606
607
608
		case 'o' :
			outfile = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
609
		case 'z' :
Thomas White's avatar
Thomas White committed
610
611
612
			indm_str = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
613
		case 'q' :
614
615
616
			intfile = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
617
		case 'p' :
Thomas White's avatar
Thomas White committed
618
619
620
			pdb = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
621
		case 'x' :
622
623
624
			prefix = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
625
		case 'j' :
626
627
628
			nthreads = atoi(optarg);
			break;

Thomas White's avatar
Thomas White committed
629
		case 'g' :
Thomas White's avatar
Thomas White committed
630
631
632
			geometry = strdup(optarg);
			break;

633
634
635
636
		case 't' :
			threshold = strtof(optarg, NULL);
			break;

637
638
639
640
641
642
643
644
645
		case 'b' :
			beam = get_beam_parameters(optarg);
			if ( beam == NULL ) {
				ERROR("Failed to load beam parameters"
				      " from '%s'\n", optarg);
				return 1;
			}
			break;

646
647
648
649
		case 2 :
			speaks = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
650
		case 0 :
Thomas White's avatar
Thomas White committed
651
652
			break;

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

673
674
675
676
677
678
679
680
681
682
683
684
685
686
	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);

687
688
689
690
691
692
693
694
695
696
697
698
699
	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
700
	free(speaks);
701

702
	if ( intfile != NULL ) {
Thomas White's avatar
Thomas White committed
703
		ReflItemList *items;
704
705
		items = read_reflections(intfile, intensities,
		                         NULL, NULL, NULL);
Thomas White's avatar
Thomas White committed
706
		delete_items(items);
707
708
709
710
	} else {
		intensities = NULL;
	}

Thomas White's avatar
Thomas White committed
711
712
713
714
	if ( pdb == NULL ) {
		pdb = strdup("molecule.pdb");
	}

715
	if ( prefix == NULL ) {
Thomas White's avatar
Thomas White committed
716
		prefix = strdup("");
717
	} else {
718
719
720
		if ( config_checkprefix ) {
			prefix = check_prefix(prefix);
		}
721
722
	}

723
	if ( nthreads == 0 ) {
724
725
726
727
		ERROR("Invalid number of threads.\n");
		return 1;
	}

Thomas White's avatar
Thomas White committed
728
729
	if ( indm_str == NULL ) {
		STATUS("You didn't specify an indexing method, so I won't"
730
731
732
		       " 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
733
734
735
736
737
		indm = INDEXING_NONE;
	} else if ( strcmp(indm_str, "none") == 0 ) {
		indm = INDEXING_NONE;
	} else if ( strcmp(indm_str, "dirax") == 0) {
		indm = INDEXING_DIRAX;
738
739
	} else if ( strcmp(indm_str, "template") == 0) {
		indm = INDEXING_TEMPLATE;
Thomas White's avatar
Thomas White committed
740
741
742
743
	} else {
		ERROR("Unrecognised indexing method '%s'\n", indm_str);
		return 1;
	}
744
	free(indm_str);
Thomas White's avatar
Thomas White committed
745

Thomas White's avatar
Thomas White committed
746
747
748
749
750
751
752
753
754
755
756
757
	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);

758
	if ( (!config_nomatch) || (indm == INDEXING_TEMPLATE) ) {
759
760
		cell = load_cell_from_pdb(pdb);
		if ( cell == NULL ) {
761
762
763
			ERROR("Couldn't read unit cell (from %s)\n", pdb);
			return 1;
		}
764
765
766
	} else {
		STATUS("No cell needed because --no-match was used.\n");
		cell = NULL;
767
	}
768
	free(pdb);
769

770
771
772
773
774
	/* Start by writing the entire command line to stdout */
	fprintf(ofh, "Command line:");
	for ( i=0; i<argc; i++ ) {
		fprintf(ofh, " %s", argv[i]);
	}
775
776
	fprintf(ofh, "\n");
	fflush(ofh);
777

778
779
780
781
782
783
784
785
	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;
	}

786
	/* Get first filename and use it to set up the indexing */
787
788
789
790
	rval = fgets(prepare_line, 1023, fh);
	if ( rval == NULL ) {
		ERROR("Failed to get filename to prepare indexing.\n");
		return 1;
791
	}
792
793
794
795
796
	chomp(prepare_line);
	snprintf(prepare_filename, 1023, "%s%s", prefix, prepare_line);
	qargs.use_this_one_instead = prepare_line;

	/* Prepare the indexer */
797
798
	ipriv = prepare_indexing(indm, cell, prepare_filename, det,
	                         nominal_photon_energy);
799
800
801
802
803
	if ( ipriv == NULL ) {
		ERROR("Failed to prepare indexing.\n");
		return 1;
	}

Thomas White's avatar
Thomas White committed
804
	gsl_set_error_handler_off();
Thomas White's avatar
Thomas White committed
805

806
807
808
809
810
811
812
813
814
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_nomatch = config_nomatch;
	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;
	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;
832
	qargs.static_args.nominal_photon_energy = nominal_photon_energy;
833
834
835
836
837
838
839
840

	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);
841

842
843
	cleanup_indexing(ipriv);

844
	free(prefix);
Thomas White's avatar
Thomas White committed
845
846
	free(det->panels);
	free(det);
847
	cell_free(cell);
848
	if ( fh != stdout ) fclose(fh);
Thomas White's avatar
Thomas White committed
849

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

Thomas White's avatar
Thomas White committed
853
854
855
856
	if ( gctx != NULL ) {
		cleanup_gpu(gctx);
	}

Thomas White's avatar
Thomas White committed
857
858
	return 0;
}