indexamajig.c 22.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 "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

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
	pthread_mutex_t *gpu_mutex;     /* Protects "gctx" */
51
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;
	int config_nomatch;
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

75
76
77
	/* Output stream */
	pthread_mutex_t *output_mutex;  /* Protects the output stream */
	FILE *ofh;
78
79
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
105
/* 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;
};


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


196
static struct image *get_simage(struct image *template, int alternate)
197
{
198
	struct image *image;
199
	struct panel panels[2];
200

201
202
	image = malloc(sizeof(*image));

203
	/* Simulate a diffraction pattern */
204
205
206
	image->twotheta = NULL;
	image->data = NULL;
	image->det = template->det;
207
208
209
	image->flags = NULL;
	image->f0_available = 0;
	image->f0 = 1.0;
210
211

	/* View head-on (unit cell is tilted) */
212
213
214
215
	image->orientation.w = 1.0;
	image->orientation.x = 0.0;
	image->orientation.y = 0.0;
	image->orientation.z = 0.0;
216
217
218

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

223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
	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;
241
		panels[1].cy = 525.0;
242
243
244
		panels[1].clen = 56.7e-2;  /* 56.7 cm */
		panels[1].res = 13333.3;   /* 75 microns/pixel */

Thomas White's avatar
Thomas White committed
245
		image->det->panels = panels;
246

247
248
249
	} else {

		/* Copy pointer to old geometry */
Thomas White's avatar
Thomas White committed
250
		image->det->panels = template->det->panels;
251
252

	}
253

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

Thomas White's avatar
Thomas White committed
260
	/* Prevent muppetry */
Thomas White's avatar
Thomas White committed
261
262
	image->cpeaks = NULL;
	image->n_cpeaks = 0;
Thomas White's avatar
Thomas White committed
263

264
265
	return image;
}
266
267


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

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

285
	hdf5_write("simulated.h5", simage->data, simage->width, simage->height,
286
287
288
289
		   H5T_NATIVE_FLOAT);
}


290
static void process_image(void *pp, int cookie)
291
{
292
	struct index_args *pargs = pp;
293
294
295
296
297
	struct hdfile *hdfile;
	struct image image;
	struct image *simage;
	float *data_for_measurement;
	size_t data_size;
298
	char *filename = pargs->filename;
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
	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;
314
315
316
317

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

324
325
326
327
328
329
	/* 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;

330
	STATUS("Processing '%s'\n", image.filename);
331

Thomas White's avatar
Thomas White committed
332
333
	pargs->sane = 0;
	pargs->indexable = 0;
334

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

343
	hdf5_read(hdfile, &image, pargs->static_args.config_satcorr);
344
345
346
347
348
349
350
351
352
353
354
355
356

	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 {
357
		memcpy(data_for_measurement, image.data, data_size);
358
359
	}

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

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

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

	/* 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,
393
		              config_verbose, pargs->static_args.ipriv);
394
395
396
397
	}

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

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

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

418
419
	simage = get_simage(&image, config_alternate);

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

	/* 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 */
439
	cell_free(image.indexed_cell);
440
441
442

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


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

457
458
459
	/* Get the next filename */
	rval = fgets(line, 1023, qargs->fh);
	if ( rval == NULL ) return NULL;
460

461
	pargs = malloc(sizeof(struct index_args));
462

463
464
	memcpy(&pargs->static_args, &qargs->static_args,
	       sizeof(struct static_index_args));
465

466
467
468
	chomp(line);
	pargs->filename = malloc(strlen(qargs->prefix) + strlen(line) + 1);
	snprintf(pargs->filename, 1023, "%s%s", qargs->prefix, line);
469

470
471
	return pargs;
}
472
473


474
475
476
477
static void finalise_image(void *qp, void *pp)
{
	struct queue_args *qargs = qp;
	struct index_args *pargs = pp;
478

479
480
	qargs->n_indexable += pargs->indexable;
	qargs->n_sane += pargs->sane;
481

482
483
	free(pargs->filename);
	free(pargs);
484
485
486
}


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

	/* Long options */
	const struct option longopts[] = {
		{"help",               0, NULL,               'h'},
		{"input",              1, NULL,               'i'},
539
		{"output",             1, NULL,               'o'},
Thomas White's avatar
Thomas White committed
540
		{"gpu",                0, &config_gpu,         1},
541
		{"no-index",           0, &config_noindex,     1},
Thomas White's avatar
Thomas White committed
542
		{"dump-peaks",         0, &config_dumpfound,   1},
543
		{"peaks",              1, NULL,                2},
544
		{"near-bragg",         0, &config_nearbragg,   1},
Thomas White's avatar
Thomas White committed
545
546
		{"write-drx",          0, &config_writedrx,    1},
		{"indexing",           1, NULL,               'z'},
Thomas White's avatar
Thomas White committed
547
		{"geometry",           1, NULL,               'g'},
548
		{"simulate",           0, &config_simulate,    1},
549
550
		{"filter-cm",          0, &config_cmfilter,    1},
		{"filter-noise",       0, &config_noisefilter, 1},
551
		{"no-match",           0, &config_nomatch,     1},
552
		{"verbose",            0, &config_verbose,     1},
553
		{"alternate",          0, &config_alternate,   1},
554
		{"intensities",        1, NULL,               'q'},
Thomas White's avatar
Thomas White committed
555
		{"pdb",                1, NULL,               'p'},
556
		{"prefix",             1, NULL,               'x'},
Thomas White's avatar
Thomas White committed
557
		{"unpolarized",        0, &config_polar,       0},
558
		{"check-sanity",       0, &config_sanity,      1},
559
560
		{"no-sat-corr",        0, &config_satcorr,     0},
		{"sat-corr",           0, &config_satcorr,     1}, /* Compat */
561
		{"no-sa",              0, &config_sa,          0},
562
		{"threshold",          1, NULL,               't'},
563
		{"no-check-prefix",    0, &config_checkprefix, 0},
564
		{"no-closer-peak",     0, &config_closer,      0},
Thomas White's avatar
Thomas White committed
565
566
567
568
		{0, 0, NULL, 0}
	};

	/* Short options */
569
	while ((c = getopt_long(argc, argv, "hi:wp:j:x:g:t:o:",
Thomas White's avatar
Thomas White committed
570
	                        longopts, NULL)) != -1) {
Thomas White's avatar
Thomas White committed
571
572

		switch (c) {
Thomas White's avatar
Thomas White committed
573
		case 'h' :
Thomas White's avatar
Thomas White committed
574
575
576
			show_help(argv[0]);
			return 0;

Thomas White's avatar
Thomas White committed
577
		case 'i' :
Thomas White's avatar
Thomas White committed
578
579
580
			filename = strdup(optarg);
			break;

581
582
583
584
		case 'o' :
			outfile = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
585
		case 'z' :
Thomas White's avatar
Thomas White committed
586
587
588
			indm_str = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
589
		case 'q' :
590
591
592
			intfile = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
593
		case 'p' :
Thomas White's avatar
Thomas White committed
594
595
596
			pdb = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
597
		case 'x' :
598
599
600
			prefix = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
601
		case 'j' :
602
603
604
			nthreads = atoi(optarg);
			break;

Thomas White's avatar
Thomas White committed
605
		case 'g' :
Thomas White's avatar
Thomas White committed
606
607
608
			geometry = strdup(optarg);
			break;

609
610
611
612
		case 't' :
			threshold = strtof(optarg, NULL);
			break;

613
614
615
616
		case 2 :
			speaks = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
617
		case 0 :
Thomas White's avatar
Thomas White committed
618
619
			break;

Thomas White's avatar
Thomas White committed
620
		default :
Thomas White's avatar
Thomas White committed
621
622
623
624
625
626
627
628
629
630
631
632
633
634
			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
635
		ERROR("Failed to open input file '%s'\n", filename);
Thomas White's avatar
Thomas White committed
636
637
		return 1;
	}
Thomas White's avatar
Thomas White committed
638
	free(filename);
Thomas White's avatar
Thomas White committed
639

640
641
642
643
644
645
646
647
648
649
650
651
652
653
	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);

654
655
656
657
658
659
660
661
662
663
664
665
666
	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
667
	free(speaks);
668

669
	if ( intfile != NULL ) {
Thomas White's avatar
Thomas White committed
670
		ReflItemList *items;
671
672
		items = read_reflections(intfile, intensities,
		                         NULL, NULL, NULL);
Thomas White's avatar
Thomas White committed
673
		delete_items(items);
674
675
676
677
	} else {
		intensities = NULL;
	}

Thomas White's avatar
Thomas White committed
678
679
680
681
	if ( pdb == NULL ) {
		pdb = strdup("molecule.pdb");
	}

682
	if ( prefix == NULL ) {
Thomas White's avatar
Thomas White committed
683
		prefix = strdup("");
684
	} else {
685
686
687
		if ( config_checkprefix ) {
			prefix = check_prefix(prefix);
		}
688
689
	}

690
	if ( nthreads == 0 ) {
691
692
693
694
		ERROR("Invalid number of threads.\n");
		return 1;
	}

Thomas White's avatar
Thomas White committed
695
696
	if ( indm_str == NULL ) {
		STATUS("You didn't specify an indexing method, so I won't"
697
698
699
		       " 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
700
701
702
703
704
		indm = INDEXING_NONE;
	} else if ( strcmp(indm_str, "none") == 0 ) {
		indm = INDEXING_NONE;
	} else if ( strcmp(indm_str, "dirax") == 0) {
		indm = INDEXING_DIRAX;
705
706
	} else if ( strcmp(indm_str, "template") == 0) {
		indm = INDEXING_TEMPLATE;
Thomas White's avatar
Thomas White committed
707
708
709
710
	} else {
		ERROR("Unrecognised indexing method '%s'\n", indm_str);
		return 1;
	}
711
	free(indm_str);
Thomas White's avatar
Thomas White committed
712

Thomas White's avatar
Thomas White committed
713
714
715
716
717
718
719
720
721
722
723
724
	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);

725
	if ( (!config_nomatch) || (indm == INDEXING_TEMPLATE) ) {
726
727
		cell = load_cell_from_pdb(pdb);
		if ( cell == NULL ) {
728
729
730
			ERROR("Couldn't read unit cell (from %s)\n", pdb);
			return 1;
		}
731
732
733
	} else {
		STATUS("No cell needed because --no-match was used.\n");
		cell = NULL;
734
	}
735
	free(pdb);
736

737
738
739
740
741
	/* Start by writing the entire command line to stdout */
	fprintf(ofh, "Command line:");
	for ( i=0; i<argc; i++ ) {
		fprintf(ofh, " %s", argv[i]);
	}
742
743
	fprintf(ofh, "\n");
	fflush(ofh);
744

745
	/* Get first filename and use it to set up the indexing */
746
747
748
749
750
751
752
753
754
755
756
757
758
759
	if ( fh != stdin ) {
		rval = fgets(prepare_line, 1023, fh);
		if ( rval == NULL ) {
			ERROR("Failed to get filename to prepare indexing.\n");
			return 1;
		}
		chomp(prepare_line);
		snprintf(prepare_filename, 1023, "%s%s", prefix, prepare_line);
		rewind(fh);
	} else {
		STATUS("Reading input filenames from stdin, so can't (yet)");
		STATUS(" use the first file for preparing the indexing.\n");
		STATUS("Stuff might break.\n");
		prepare_filename[0] = '\0';
760
	}
761
	ipriv = prepare_indexing(indm, cell, prepare_filename, det);
762
763
764
765
766
	if ( ipriv == NULL ) {
		ERROR("Failed to prepare indexing.\n");
		return 1;
	}

Thomas White's avatar
Thomas White committed
767
	gsl_set_error_handler_off();
Thomas White's avatar
Thomas White committed
768

769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
	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;

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

804
805
	cleanup_indexing(ipriv);

806
	free(prefix);
Thomas White's avatar
Thomas White committed
807
808
	free(det->panels);
	free(det);
809
	cell_free(cell);
810
	if ( fh != stdout ) fclose(fh);
Thomas White's avatar
Thomas White committed
811

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

Thomas White's avatar
Thomas White committed
815
816
817
818
	if ( gctx != NULL ) {
		cleanup_gpu(gctx);
	}

Thomas White's avatar
Thomas White committed
819
820
	return 0;
}