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;
Thomas White's avatar
Thomas White committed
71
72
	IndexingMethod *indm;
	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
	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;
Thomas White's avatar
Thomas White committed
328
	IndexingMethod *indm = pargs->static_args.indm;
329
	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
548
549
	IndexingMethod *indm;
	IndexingPrivate **ipriv;
	int indexer_needs_cell;
	int reduction_needs_cell;
Thomas White's avatar
Thomas White committed
550
	char *indm_str = NULL;
551
552
	UnitCell *cell;
	double *intensities = NULL;
553
	unsigned char *flags;
554
	char *intfile = NULL;
Thomas White's avatar
Thomas White committed
555
	char *pdb = NULL;
556
	char *prefix = NULL;
557
	char *speaks = NULL;
558
559
	char *scellr = NULL;
	int cellr;
560
	int peaks;
561
562
	int nthreads = 1;
	int i;
563
	pthread_mutex_t output_mutex = PTHREAD_MUTEX_INITIALIZER;
564
	pthread_mutex_t gpu_mutex = PTHREAD_MUTEX_INITIALIZER;
565
566
	char prepare_line[1024];
	char prepare_filename[1024];
567
	struct queue_args qargs;
568
569
	struct beam_params *beam = NULL;
	double nominal_photon_energy;
570
	int gpu_dev = -1;
571
	char *sym = NULL;
Thomas White's avatar
Thomas White committed
572
573
574
575
576

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

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

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

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

623
624
625
626
		case 'o' :
			outfile = strdup(optarg);
			break;

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

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

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

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

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

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

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

655
656
657
658
		case 'y' :
			sym = strdup(optarg);
			break;

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

668
669
670
671
		case 2 :
			speaks = strdup(optarg);
			break;

672
673
674
675
		case 3 :
			scellr = strdup(optarg);
			break;

676
677
678
679
		case 4 :
			min_gradient = strtof(optarg, NULL);
			break;

680
681
682
683
		case 5 :
			gpu_dev = atoi(optarg);
			break;

Thomas White's avatar
Thomas White committed
684
		case 0 :
Thomas White's avatar
Thomas White committed
685
686
			break;

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

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

721
722
	if ( sym == NULL ) sym = strdup("1");

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

738
	if ( intfile != NULL ) {
739

Thomas White's avatar
Thomas White committed
740
		ReflItemList *items;
741
742
		int i;

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

		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
752
753
754
755
756
757
		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
758
		delete_items(items);
759

760
	} else {
761

762
		intensities = NULL;
763
764
		flags = NULL;

765
766
	}

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

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

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

Thomas White's avatar
Thomas White committed
784
785
	if ( indm_str == NULL ) {
		STATUS("You didn't specify an indexing method, so I won't"
786
787
788
		       " 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
789
		indm = NULL;
Thomas White's avatar
Thomas White committed
790
	} else {
Thomas White's avatar
Thomas White committed
791
792
793
794
795
796
		indm = build_indexer_list(indm_str, &indexer_needs_cell);
		if ( indm == NULL ) {
			ERROR("Invalid indexer list '%s'\n", indm_str);
			return 1;
		}
		free(indm_str);
Thomas White's avatar
Thomas White committed
797
798
	}

Thomas White's avatar
Thomas White committed
799
	reduction_needs_cell = 0;
800
801
802
803
	if ( scellr == NULL ) {
		STATUS("You didn't specify a cell reduction method, so I'm"
		       " going to use 'reduce'.\n");
		cellr = CELLR_REDUCE;
Thomas White's avatar
Thomas White committed
804
		reduction_needs_cell = 1;
805
806
807
808
	} else if ( strcmp(scellr, "none") == 0 ) {
		cellr = CELLR_NONE;
	} else if ( strcmp(scellr, "reduce") == 0) {
		cellr = CELLR_REDUCE;
Thomas White's avatar
Thomas White committed
809
		reduction_needs_cell = 1;
810
811
	} else if ( strcmp(scellr, "compare") == 0) {
		cellr = CELLR_COMPARE;
Thomas White's avatar
Thomas White committed
812
		reduction_needs_cell = 1;
813
814
815
816
817
818
	} else {
		ERROR("Unrecognised cell reduction method '%s'\n", scellr);
		return 1;
	}
	free(scellr);  /* free(NULL) is OK. */

Thomas White's avatar
Thomas White committed
819
820
821
822
823
824
825
826
827
828
829
830
	if ( geometry == NULL ) {
		ERROR("You need to specify a geometry file with --geometry\n");
		return 1;
	}

	det = get_detector_geometry(geometry);
	if ( det == NULL ) {
		ERROR("Failed to read detector geometry from '%s'\n", geometry);
		return 1;
	}
	free(geometry);

Thomas White's avatar
Thomas White committed
831
	if ( reduction_needs_cell || indexer_needs_cell ) {
832
833
		cell = load_cell_from_pdb(pdb);
		if ( cell == NULL ) {
834
835
836
			ERROR("Couldn't read unit cell (from %s)\n", pdb);
			return 1;
		}
837
	} else {
Thomas White's avatar
Thomas White committed
838
839
		STATUS("No cell needed for these choices of indexing"
		       " and reduction.\n");
840
		cell = NULL;
841
	}
842
	free(pdb);
843

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

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

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

	/* Prepare the indexer */
Thomas White's avatar
Thomas White committed
871
872
873
874
875
876
877
878
879
	if ( indm != NULL ) {
		ipriv = prepare_indexing(indm, cell, prepare_filename, det,
		                         nominal_photon_energy);
		if ( ipriv == NULL ) {
			ERROR("Failed to prepare indexing.\n");
			return 1;
		}
	} else {
		ipriv = NULL;
880
881
	}

Thomas White's avatar
Thomas White committed
882
	gsl_set_error_handler_off();
Thomas White's avatar
Thomas White committed
883

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

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

924
925
	cleanup_indexing(ipriv);

926
	free(prefix);
Thomas White's avatar
Thomas White committed
927
928
	free(det->panels);
	free(det);
929
	cell_free(cell);
930
	if ( fh != stdout ) fclose(fh);
931
	free(sym);
Thomas White's avatar
Thomas White committed
932

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

Thomas White's avatar
Thomas White committed
936
937
938
939
	if ( gctx != NULL ) {
		cleanup_gpu(gctx);
	}

Thomas White's avatar
Thomas White committed
940
941
	return 0;
}