indexamajig.c 26.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
#include "symmetry.h"
41

42
43
44
45
46
47
48

enum {
	PEAK_ZAEF,
	PEAK_HDF5,
};


49
50
/* Information about the indexing process which is common to all patterns */
struct static_index_args
51
{
52
	pthread_mutex_t *gpu_mutex;     /* Protects "gctx" */
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;
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;
69
	float min_gradient;
Thomas White's avatar
Thomas White committed
70
	struct detector *det;
71
	IndexingMethod indm;
72
	IndexingPrivate *ipriv;
73
	const double *intensities;
74
75
	const unsigned char *flags;
	const char *sym;  /* Symmetry of "intensities" and "flags" */
76
	struct gpu_context *gctx;
77
	int gpu_dev;
78
	int peaks;
79
	int cellr;
80
	double nominal_photon_energy;
81

82
83
84
	/* Output stream */
	pthread_mutex_t *output_mutex;  /* Protects the output stream */
	FILE *ofh;
85
86
87
};


88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
/* 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;
110
111

	char *use_this_one_instead;
112
113
114
};


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


216
static struct image *get_simage(struct image *template, int alternate)
217
{
218
	struct image *image;
219
	struct panel panels[2];
220

221
222
	image = malloc(sizeof(*image));

223
	/* Simulate a diffraction pattern */
224
225
226
	image->twotheta = NULL;
	image->data = NULL;
	image->det = template->det;
227
228
229
	image->flags = NULL;
	image->f0_available = 0;
	image->f0 = 1.0;
230
231
232

	/* Detector geometry for the simulation
	 * - not necessarily the same as the original. */
233
234
	image->width = 1024;
	image->height = 1024;
Thomas White's avatar
Thomas White committed
235
	image->det->n_panels = 2;
236

237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
	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;
255
		panels[1].cy = 525.0;
256
257
258
		panels[1].clen = 56.7e-2;  /* 56.7 cm */
		panels[1].res = 13333.3;   /* 75 microns/pixel */

Thomas White's avatar
Thomas White committed
259
		image->det->panels = panels;
260

261
262
263
	} else {

		/* Copy pointer to old geometry */
Thomas White's avatar
Thomas White committed
264
		image->det->panels = template->det->panels;
265
266

	}
267

268
	image->lambda = ph_en_to_lambda(eV_to_J(1.8e3));
269
	image->features = template->features;
270
	image->filename = template->filename;
271
	image->indexed_cell = template->indexed_cell;
Thomas White's avatar
Thomas White committed
272
	image->f0 = template->f0;
273

Thomas White's avatar
Thomas White committed
274
	/* Prevent muppetry */
Thomas White's avatar
Thomas White committed
275
276
	image->cpeaks = NULL;
	image->n_cpeaks = 0;
Thomas White's avatar
Thomas White committed
277

278
279
	return image;
}
280
281


282
static void simulate_and_write(struct image *simage, struct gpu_context **gctx,
283
284
285
                               const double *intensities,
                               const unsigned char *flags, UnitCell *cell,
                               int gpu_dev, const char *sym)
286
{
287
288
289
	/* Set up GPU if necessary.
	 * Unfortunately, setup has to go here since until now we don't know
	 * enough about the situation. */
Thomas White's avatar
Thomas White committed
290
	if ( (gctx != NULL) && (*gctx == NULL) ) {
291
		*gctx = setup_gpu(0, simage, intensities, flags, sym, gpu_dev);
292
293
	}

Thomas White's avatar
Thomas White committed
294
	if ( (gctx != NULL) && (*gctx != NULL) ) {
295
		get_diffraction_gpu(*gctx, simage, 24, 24, 40, cell);
296
	} else {
297
		get_diffraction(simage, 24, 24, 40,
298
299
		                intensities, NULL, flags, cell, 0,
		                GRADIENT_MOSAIC, sym);
300
	}
301
	record_image(simage, 0);
302

303
	hdf5_write("simulated.h5", simage->data, simage->width, simage->height,
304
305
306
307
		   H5T_NATIVE_FLOAT);
}


308
static void process_image(void *pp, int cookie)
309
{
310
	struct index_args *pargs = pp;
311
312
313
314
315
	struct hdfile *hdfile;
	struct image image;
	struct image *simage;
	float *data_for_measurement;
	size_t data_size;
316
	char *filename = pargs->filename;
317
318
319
320
321
322
323
324
325
326
327
328
329
	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;
330
	const unsigned char *flags = pargs->static_args.flags;
331
	struct gpu_context *gctx = pargs->static_args.gctx;
332
	const char *sym = pargs->static_args.sym;
333
334
335
336

	image.features = NULL;
	image.data = NULL;
	image.indexed_cell = NULL;
337
	image.id = cookie;
338
	image.filename = filename;
Thomas White's avatar
Thomas White committed
339
340
	image.cpeaks = NULL;
	image.n_cpeaks = 0;
341
	image.det = pargs->static_args.det;
342

343
	STATUS("Processing '%s'\n", image.filename);
344

Thomas White's avatar
Thomas White committed
345
346
	pargs->sane = 0;
	pargs->indexable = 0;
347

348
349
	hdfile = hdfile_open(filename);
	if ( hdfile == NULL ) {
Thomas White's avatar
Thomas White committed
350
		return;
351
352
	} else if ( hdfile_set_first_image(hdfile, "/") ) {
		ERROR("Couldn't select path\n");
Thomas White's avatar
Thomas White committed
353
		return;
354
355
	}

356
357
	hdf5_read(hdfile, &image, pargs->static_args.config_satcorr,
	          pargs->static_args.nominal_photon_energy);
358
359
360
361
362
363
364
365
366
367
368
369
370

	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 {
371
		memcpy(data_for_measurement, image.data, data_size);
372
373
	}

374
	switch ( pargs->static_args.peaks )
375
376
377
378
379
380
381
382
383
	{
	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 :
384
385
		search_peaks(&image, pargs->static_args.threshold,
		             pargs->static_args.min_gradient);
386
387
		break;
	}
Thomas White's avatar
Thomas White committed
388
389
390
391
392
393

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

394
	if ( config_dumpfound ) {
395
396
		dump_peaks(&image, pargs->static_args.ofh,
		           pargs->static_args.output_mutex);
397
	}
398
399
400
401
402
403
404
405
406

	/* 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) ) {
407
		index_pattern(&image, cell, indm, pargs->static_args.cellr,
408
		              config_verbose, pargs->static_args.ipriv);
409
410
411
412
	}

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

Thomas White's avatar
Thomas White committed
415
	/* Sanity check */
416
	if ( pargs->static_args.config_sanity
417
	  && !peak_sanity_check(&image, image.indexed_cell, 0, 0.1) ) {
Thomas White's avatar
Thomas White committed
418
419
420
		STATUS("Failed peak sanity check.\n");
		goto done;
	} else {
Thomas White's avatar
Thomas White committed
421
		pargs->sane = 1;
Thomas White's avatar
Thomas White committed
422
423
	}

424
425
	/* Measure intensities if requested */
	if ( config_nearbragg ) {
426
		output_intensities(&image, image.indexed_cell,
427
428
429
430
		                   pargs->static_args.output_mutex,
		                   config_polar, pargs->static_args.config_sa,
		                   pargs->static_args.config_closer,
		                   pargs->static_args.ofh, 0, 0.1);
431
432
	}

433
434
	simage = get_simage(&image, config_alternate);

435
436
437
	/* Simulate if requested */
	if ( config_simulate ) {
		if ( config_gpu ) {
438
			pthread_mutex_lock(pargs->static_args.gpu_mutex);
439
			simulate_and_write(simage, &gctx, intensities, flags,
440
			                   image.indexed_cell,
441
			                   pargs->static_args.gpu_dev, sym);
442
			pthread_mutex_unlock(pargs->static_args.gpu_mutex);
443
		} else {
444
445
			simulate_and_write(simage, NULL, intensities, flags,
			                   image.indexed_cell, 0, sym);
446
447
448
449
450
451
452
453
454
		}
	}

	/* 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 */
455
	cell_free(image.indexed_cell);
456
457
458

done:
	free(image.data);
Thomas White's avatar
Thomas White committed
459
	free(image.flags);
460
	image_feature_list_free(image.features);
Thomas White's avatar
Thomas White committed
461
	free(image.cpeaks);
462
463
464
465
	hdfile_close(hdfile);
}


466
static void *get_image(void *qp)
467
{
468
469
470
471
	char line[1024];
	struct index_args *pargs;
	char *rval;
	struct queue_args *qargs = qp;
472

473
	/* Initialise new task arguments */
474
475
476
	pargs = malloc(sizeof(struct index_args));
	memcpy(&pargs->static_args, &qargs->static_args,
	       sizeof(struct static_index_args));
477

478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
	/* 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);

	}
498

499
500
	return pargs;
}
501
502


503
504
505
506
static void finalise_image(void *qp, void *pp)
{
	struct queue_args *qargs = qp;
	struct index_args *pargs = pp;
507

508
509
	qargs->n_indexable += pargs->indexable;
	qargs->n_sane += pargs->sane;
510

511
512
	free(pargs->filename);
	free(pargs);
513
514
515
}


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

	/* Long options */
	const struct option longopts[] = {
		{"help",               0, NULL,               'h'},
		{"input",              1, NULL,               'i'},
575
		{"output",             1, NULL,               'o'},
Thomas White's avatar
Thomas White committed
576
		{"gpu",                0, &config_gpu,         1},
577
		{"no-index",           0, &config_noindex,     1},
Thomas White's avatar
Thomas White committed
578
		{"dump-peaks",         0, &config_dumpfound,   1},
579
		{"peaks",              1, NULL,                2},
580
		{"cell-reduction",     1, NULL,                3},
581
		{"near-bragg",         0, &config_nearbragg,   1},
Thomas White's avatar
Thomas White committed
582
583
		{"write-drx",          0, &config_writedrx,    1},
		{"indexing",           1, NULL,               'z'},
Thomas White's avatar
Thomas White committed
584
		{"geometry",           1, NULL,               'g'},
585
		{"beam",               1, NULL,               'b'},
586
		{"simulate",           0, &config_simulate,    1},
587
588
		{"filter-cm",          0, &config_cmfilter,    1},
		{"filter-noise",       0, &config_noisefilter, 1},
589
		{"verbose",            0, &config_verbose,     1},
590
		{"alternate",          0, &config_alternate,   1},
591
		{"intensities",        1, NULL,               'q'},
592
		{"symmetry",           1, NULL,               'y'},
Thomas White's avatar
Thomas White committed
593
		{"pdb",                1, NULL,               'p'},
594
		{"prefix",             1, NULL,               'x'},
Thomas White's avatar
Thomas White committed
595
		{"unpolarized",        0, &config_polar,       0},
596
		{"check-sanity",       0, &config_sanity,      1},
597
598
		{"no-sat-corr",        0, &config_satcorr,     0},
		{"sat-corr",           0, &config_satcorr,     1}, /* Compat */
599
		{"no-sa",              0, &config_sa,          0},
600
		{"threshold",          1, NULL,               't'},
601
		{"min-gradient",       1, NULL,                4},
602
		{"no-check-prefix",    0, &config_checkprefix, 0},
603
		{"no-closer-peak",     0, &config_closer,      0},
604
		{"gpu-dev",            1, NULL,                5},
Thomas White's avatar
Thomas White committed
605
606
607
608
		{0, 0, NULL, 0}
	};

	/* Short options */
609
	while ((c = getopt_long(argc, argv, "hi:wp:j:x:g:t:o:b:y:",
Thomas White's avatar
Thomas White committed
610
	                        longopts, NULL)) != -1) {
Thomas White's avatar
Thomas White committed
611
612

		switch (c) {
Thomas White's avatar
Thomas White committed
613
		case 'h' :
Thomas White's avatar
Thomas White committed
614
615
616
			show_help(argv[0]);
			return 0;

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

621
622
623
624
		case 'o' :
			outfile = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
625
		case 'z' :
Thomas White's avatar
Thomas White committed
626
627
628
			indm_str = strdup(optarg);
			break;

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

Thomas White's avatar
Thomas White committed
633
		case 'p' :
Thomas White's avatar
Thomas White committed
634
635
636
			pdb = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
637
		case 'x' :
638
639
640
			prefix = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
641
		case 'j' :
642
643
644
			nthreads = atoi(optarg);
			break;

Thomas White's avatar
Thomas White committed
645
		case 'g' :
Thomas White's avatar
Thomas White committed
646
647
648
			geometry = strdup(optarg);
			break;

649
650
651
652
		case 't' :
			threshold = strtof(optarg, NULL);
			break;

653
654
655
656
		case 'y' :
			sym = strdup(optarg);
			break;

657
658
659
660
661
662
663
664
665
		case 'b' :
			beam = get_beam_parameters(optarg);
			if ( beam == NULL ) {
				ERROR("Failed to load beam parameters"
				      " from '%s'\n", optarg);
				return 1;
			}
			break;

666
667
668
669
		case 2 :
			speaks = strdup(optarg);
			break;

670
671
672
673
		case 3 :
			scellr = strdup(optarg);
			break;

674
675
676
677
		case 4 :
			min_gradient = strtof(optarg, NULL);
			break;

678
679
680
681
		case 5 :
			gpu_dev = atoi(optarg);
			break;

Thomas White's avatar
Thomas White committed
682
		case 0 :
Thomas White's avatar
Thomas White committed
683
684
			break;

Thomas White's avatar
Thomas White committed
685
		default :
Thomas White's avatar
Thomas White committed
686
687
688
689
690
691
692
693
694
695
696
697
698
699
			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
700
		ERROR("Failed to open input file '%s'\n", filename);
Thomas White's avatar
Thomas White committed
701
702
		return 1;
	}
Thomas White's avatar
Thomas White committed
703
	free(filename);
Thomas White's avatar
Thomas White committed
704

705
706
707
708
709
710
711
712
713
714
715
716
717
718
	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);

719
720
	if ( sym == NULL ) sym = strdup("1");

721
722
723
724
725
726
727
728
729
730
731
732
733
	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
734
	free(speaks);
735

736
	if ( intfile != NULL ) {
737

Thomas White's avatar
Thomas White committed
738
		ReflItemList *items;
739
740
		int i;

741
742
		items = read_reflections(intfile, intensities,
		                         NULL, NULL, NULL);
743
744
745
746
747
748
749

		flags = new_list_flag();
		for ( i=0; i<num_items(items); i++ ) {
			struct refl_item *it = get_item(items, i);
			set_flag(flags, it->h, it->k, it->l, 1);
		}

Thomas White's avatar
Thomas White committed
750
751
752
753
754
755
		if ( check_symmetry(items, sym) ) {
			ERROR("The input reflection list does not appear to"
			      " have symmetry %s\n", sym);
			return 1;
		}

Thomas White's avatar
Thomas White committed
756
		delete_items(items);
757

758
	} else {
759

760
		intensities = NULL;
761
762
		flags = NULL;

763
764
	}

Thomas White's avatar
Thomas White committed
765
766
767
768
	if ( pdb == NULL ) {
		pdb = strdup("molecule.pdb");
	}

769
	if ( prefix == NULL ) {
Thomas White's avatar
Thomas White committed
770
		prefix = strdup("");
771
	} else {
772
773
774
		if ( config_checkprefix ) {
			prefix = check_prefix(prefix);
		}
775
776
	}

777
	if ( nthreads == 0 ) {
778
779
780
781
		ERROR("Invalid number of threads.\n");
		return 1;
	}

Thomas White's avatar
Thomas White committed
782
783
	if ( indm_str == NULL ) {
		STATUS("You didn't specify an indexing method, so I won't"
784
785
786
		       " 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
787
788
789
790
791
		indm = INDEXING_NONE;
	} else if ( strcmp(indm_str, "none") == 0 ) {
		indm = INDEXING_NONE;
	} else if ( strcmp(indm_str, "dirax") == 0) {
		indm = INDEXING_DIRAX;
792
793
	} else if ( strcmp(indm_str, "mosflm") == 0) {
		indm = INDEXING_MOSFLM;
794
795
	} else if ( strcmp(indm_str, "template") == 0) {
		indm = INDEXING_TEMPLATE;
Thomas White's avatar
Thomas White committed
796
797
798
799
	} else {
		ERROR("Unrecognised indexing method '%s'\n", indm_str);
		return 1;
	}
800
	free(indm_str);
Thomas White's avatar
Thomas White committed
801

802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
	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
818
819
820
821
822
823
824
825
826
827
828
829
	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);

830
	if ( (cellr != CELLR_NONE) || (indm == INDEXING_TEMPLATE) ) {
831
832
		cell = load_cell_from_pdb(pdb);
		if ( cell == NULL ) {
833
834
835
			ERROR("Couldn't read unit cell (from %s)\n", pdb);
			return 1;
		}
836
837
838
	} else {
		STATUS("No cell needed because --no-match was used.\n");
		cell = NULL;
839
	}
840
	free(pdb);
841

842
843
844
845
846
	/* Start by writing the entire command line to stdout */
	fprintf(ofh, "Command line:");
	for ( i=0; i<argc; i++ ) {
		fprintf(ofh, " %s", argv[i]);
	}
847
848
	fprintf(ofh, "\n");
	fflush(ofh);
849

850
851
852
853
854
855
856
857
	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;
	}

858
	/* Get first filename and use it to set up the indexing */
859
860
861
862
	rval = fgets(prepare_line, 1023, fh);
	if ( rval == NULL ) {
		ERROR("Failed to get filename to prepare indexing.\n");
		return 1;
863
	}
864
865
866
867
868
	chomp(prepare_line);
	snprintf(prepare_filename, 1023, "%s%s", prefix, prepare_line);
	qargs.use_this_one_instead = prepare_line;

	/* Prepare the indexer */
869
870
	ipriv = prepare_indexing(indm, cell, prepare_filename, det,
	                         nominal_photon_energy);
871
872
873
874
875
	if ( ipriv == NULL ) {
		ERROR("Failed to prepare indexing.\n");
		return 1;
	}

Thomas White's avatar
Thomas White committed
876
	gsl_set_error_handler_off();
Thomas White's avatar
Thomas White committed
877

878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
	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;
894
	qargs.static_args.cellr = cellr;
895
	qargs.static_args.threshold = threshold;
896
	qargs.static_args.min_gradient = min_gradient;
897
898
899
900
	qargs.static_args.det = det;
	qargs.static_args.indm = indm;
	qargs.static_args.ipriv = ipriv;
	qargs.static_args.intensities = intensities;
901
902
	qargs.static_args.flags = flags;
	qargs.static_args.sym = sym;
903
	qargs.static_args.gctx = gctx;
904
	qargs.static_args.gpu_dev = gpu_dev;
905
906
907
	qargs.static_args.peaks = peaks;
	qargs.static_args.output_mutex = &output_mutex;
	qargs.static_args.ofh = ofh;
908
	qargs.static_args.nominal_photon_energy = nominal_photon_energy;
909
910
911
912
913
914
915
916

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

918
919
	cleanup_indexing(ipriv);

920
	free(prefix);
Thomas White's avatar
Thomas White committed
921
922
	free(det->panels);
	free(det);
923
	cell_free(cell);
924
	if ( fh != stdout ) fclose(fh);
925
	free(sym);
Thomas White's avatar
Thomas White committed
926

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

Thomas White's avatar
Thomas White committed
930
931
932
933
	if ( gctx != NULL ) {
		cleanup_gpu(gctx);
	}

Thomas White's avatar
Thomas White committed
934
935
	return 0;
}