indexamajig.c 24.1 KB
Newer Older
Thomas White's avatar
Thomas White committed
1
/*
2
 * indexamajig.c
Thomas White's avatar
Thomas White committed
3
4
5
 *
 * Find hits, index patterns, output hkl+intensity etc.
 *
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"
Thomas White's avatar
Thomas White committed
38
39


40
41
#define MAX_THREADS (96)

42
43
44
45
46
47
48

enum {
	PEAK_ZAEF,
	PEAK_HDF5,
};


49
50
struct process_args
{
51
	/* Input */
52
	char *filename;
Thomas White's avatar
Thomas White committed
53
	int id;
54
	pthread_mutex_t *gpu_mutex;     /* Protects "gctx" */
55
56
57
58
59
60
61
62
63
64
65
	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
66
	int config_polar;
67
	int config_sanity;
68
	int config_satcorr;
69
	int config_sa;
70
	int config_closer;
71
	float threshold;
Thomas White's avatar
Thomas White committed
72
	struct detector *det;
73
	IndexingMethod indm;
74
	IndexingPrivate *ipriv;
75
76
	const double *intensities;
	struct gpu_context *gctx;
77
	int peaks;
78
79
80
81
82
83

	/* Thread control and output */
	pthread_mutex_t control_mutex;  /* Protects the scary stuff below */
	int start;
	int finish;
	int done;
Thomas White's avatar
Thomas White committed
84
85
	int indexable;
	int sane;
86
87
88
89

	/* Output stream */
	pthread_mutex_t *output_mutex;  /* Protects the output stream */
	FILE *ofh;
90
91
92
};


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


181
static struct image *get_simage(struct image *template, int alternate)
182
{
183
	struct image *image;
184
	struct panel panels[2];
185

186
187
	image = malloc(sizeof(*image));

188
	/* Simulate a diffraction pattern */
189
190
191
	image->twotheta = NULL;
	image->data = NULL;
	image->det = template->det;
192
193
194
	image->flags = NULL;
	image->f0_available = 0;
	image->f0 = 1.0;
195
196

	/* View head-on (unit cell is tilted) */
197
198
199
200
	image->orientation.w = 1.0;
	image->orientation.x = 0.0;
	image->orientation.y = 0.0;
	image->orientation.z = 0.0;
201
202
203

	/* Detector geometry for the simulation
	 * - not necessarily the same as the original. */
204
205
	image->width = 1024;
	image->height = 1024;
Thomas White's avatar
Thomas White committed
206
	image->det->n_panels = 2;
207

208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
	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;
226
		panels[1].cy = 525.0;
227
228
229
		panels[1].clen = 56.7e-2;  /* 56.7 cm */
		panels[1].res = 13333.3;   /* 75 microns/pixel */

Thomas White's avatar
Thomas White committed
230
		image->det->panels = panels;
231

232
233
234
	} else {

		/* Copy pointer to old geometry */
Thomas White's avatar
Thomas White committed
235
		image->det->panels = template->det->panels;
236
237

	}
238

239
	image->lambda = ph_en_to_lambda(eV_to_J(1.8e3));
240
	image->features = template->features;
241
	image->filename = template->filename;
242
	image->indexed_cell = template->indexed_cell;
Thomas White's avatar
Thomas White committed
243
	image->f0 = template->f0;
244

Thomas White's avatar
Thomas White committed
245
246
247
248
	/* Prevent muppetry */
	image->hits = NULL;
	image->n_hits = 0;

249
250
	return image;
}
251
252


253
static void simulate_and_write(struct image *simage, struct gpu_context **gctx,
Thomas White's avatar
Thomas White committed
254
                               const double *intensities, UnitCell *cell)
255
{
256
	/* Set up GPU if necessary */
Thomas White's avatar
Thomas White committed
257
	if ( (gctx != NULL) && (*gctx == NULL) ) {
Thomas White's avatar
Thomas White committed
258
		*gctx = setup_gpu(0, simage, intensities);
259
260
	}

Thomas White's avatar
Thomas White committed
261
	if ( (gctx != NULL) && (*gctx != NULL) ) {
262
		get_diffraction_gpu(*gctx, simage, 24, 24, 40, cell);
263
	} else {
264
		get_diffraction(simage, 24, 24, 40,
Thomas White's avatar
Thomas White committed
265
		                intensities, NULL, cell, 0,
266
		                GRADIENT_MOSAIC);
267
	}
268
	record_image(simage, 0);
269

270
	hdf5_write("simulated.h5", simage->data, simage->width, simage->height,
271
272
273
274
		   H5T_NATIVE_FLOAT);
}


Thomas White's avatar
Thomas White committed
275
static void process_image(struct process_args *pargs)
276
277
278
279
280
281
{
	struct hdfile *hdfile;
	struct image image;
	struct image *simage;
	float *data_for_measurement;
	size_t data_size;
282
283
284
285
286
287
288
289
290
291
292
293
	const char *filename = pargs->filename;
	UnitCell *cell = pargs->cell;
	int config_cmfilter = pargs->config_cmfilter;
	int config_noisefilter = pargs->config_noisefilter;
	int config_writedrx = pargs->config_writedrx;
	int config_dumpfound = pargs->config_dumpfound;
	int config_verbose = pargs->config_verbose;
	int config_alternate  = pargs->config_alternate;
	int config_nearbragg = pargs->config_nearbragg;
	int config_gpu = pargs->config_gpu;
	int config_simulate = pargs->config_simulate;
	int config_nomatch = pargs->config_nomatch;
Thomas White's avatar
Thomas White committed
294
	int config_polar = pargs->config_polar;
295
296
297
	IndexingMethod indm = pargs->indm;
	const double *intensities = pargs->intensities;
	struct gpu_context *gctx = pargs->gctx;
298
299
300
301

	image.features = NULL;
	image.data = NULL;
	image.indexed_cell = NULL;
Thomas White's avatar
Thomas White committed
302
	image.id = pargs->id;
303
	image.filename = filename;
Thomas White's avatar
Thomas White committed
304
305
	image.hits = NULL;
	image.n_hits = 0;
Thomas White's avatar
Thomas White committed
306
	image.det = pargs->det;
307

308
309
310
311
312
313
	/* 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;

314
	STATUS("Processing '%s'\n", image.filename);
315

Thomas White's avatar
Thomas White committed
316
317
	pargs->sane = 0;
	pargs->indexable = 0;
318

319
320
	hdfile = hdfile_open(filename);
	if ( hdfile == NULL ) {
Thomas White's avatar
Thomas White committed
321
		return;
322
323
	} else if ( hdfile_set_first_image(hdfile, "/") ) {
		ERROR("Couldn't select path\n");
Thomas White's avatar
Thomas White committed
324
		return;
325
326
	}

327
	hdf5_read(hdfile, &image, pargs->config_satcorr);
328
329
330
331
332
333
334
335
336
337
338
339
340

	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 {
341
		memcpy(data_for_measurement, image.data, data_size);
342
343
	}

344
345
346
347
348
349
350
351
352
353
354
355
356
	switch ( pargs->peaks )
	{
	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 :
		search_peaks(&image, pargs->threshold);
		break;
	}
Thomas White's avatar
Thomas White committed
357
358
359
360
361
362

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

363
364
365
	if ( config_dumpfound ) {
		dump_peaks(&image, pargs->ofh, pargs->output_mutex);
	}
366
367
368
369
370
371
372
373
374
375

	/* 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,
376
		              config_verbose, pargs->ipriv);
377
378
379
380
	}

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

Thomas White's avatar
Thomas White committed
383
	/* Sanity check */
384
	if ( pargs->config_sanity
385
	  && !peak_sanity_check(&image, image.indexed_cell, 0, 0.1) ) {
Thomas White's avatar
Thomas White committed
386
387
388
		STATUS("Failed peak sanity check.\n");
		goto done;
	} else {
Thomas White's avatar
Thomas White committed
389
		pargs->sane = 1;
Thomas White's avatar
Thomas White committed
390
391
	}

392
393
	/* Measure intensities if requested */
	if ( config_nearbragg ) {
394
		output_intensities(&image, image.indexed_cell,
Thomas White's avatar
Thomas White committed
395
		                   pargs->output_mutex, config_polar,
396
		                   pargs->config_sa, pargs->config_closer,
397
		                   pargs->ofh, 0, 0.1);
398
399
	}

400
401
	simage = get_simage(&image, config_alternate);

402
403
404
	/* Simulate if requested */
	if ( config_simulate ) {
		if ( config_gpu ) {
405
			pthread_mutex_lock(pargs->gpu_mutex);
406
			simulate_and_write(simage, &gctx, intensities,
Thomas White's avatar
Thomas White committed
407
			                   image.indexed_cell);
408
			pthread_mutex_unlock(pargs->gpu_mutex);
409
410
		} else {
			simulate_and_write(simage, NULL, intensities,
Thomas White's avatar
Thomas White committed
411
			                   image.indexed_cell);
412
413
414
415
416
417
418
419
420
		}
	}

	/* 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 */
421
	cell_free(image.indexed_cell);
422
423
424

done:
	free(image.data);
Thomas White's avatar
Thomas White committed
425
	free(image.flags);
426
	image_feature_list_free(image.features);
Thomas White's avatar
Thomas White committed
427
	free(image.hits);
428
429
430
431
	hdfile_close(hdfile);
}


432
433
434
435
436
437
438
439
440
static void *worker_thread(void *pargsv)
{
	struct process_args *pargs = pargsv;
	int finish;

	do {

		int wakeup;

Thomas White's avatar
Thomas White committed
441
		process_image(pargs);
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465

		pthread_mutex_lock(&pargs->control_mutex);
		pargs->done = 1;
		pargs->start = 0;
		pthread_mutex_unlock(&pargs->control_mutex);

		/* Go to sleep until told to exit or process next image */
		do {

			pthread_mutex_lock(&pargs->control_mutex);
			/* Either of these can result in the thread waking up */
			wakeup = pargs->start || pargs->finish;
			finish = pargs->finish;
			pthread_mutex_unlock(&pargs->control_mutex);
			usleep(20000);

		} while ( !wakeup );

	} while ( !pargs->finish );

	return NULL;
}


Thomas White's avatar
Thomas White committed
466
467
468
int main(int argc, char *argv[])
{
	int c;
Thomas White's avatar
Thomas White committed
469
	struct gpu_context *gctx = NULL;
Thomas White's avatar
Thomas White committed
470
	char *filename = NULL;
471
	char *outfile = NULL;
Thomas White's avatar
Thomas White committed
472
	FILE *fh;
473
	FILE *ofh;
474
	char *rval = NULL;
Thomas White's avatar
Thomas White committed
475
	int n_images;
Thomas White's avatar
Thomas White committed
476
	int n_indexable;
Thomas White's avatar
Thomas White committed
477
	int n_sane;
478
	int config_noindex = 0;
Thomas White's avatar
Thomas White committed
479
	int config_dumpfound = 0;
480
	int config_nearbragg = 0;
Thomas White's avatar
Thomas White committed
481
	int config_writedrx = 0;
482
	int config_simulate = 0;
483
484
	int config_cmfilter = 0;
	int config_noisefilter = 0;
485
	int config_nomatch = 0;
Thomas White's avatar
Thomas White committed
486
	int config_gpu = 0;
487
	int config_verbose = 0;
488
	int config_alternate = 0;
Thomas White's avatar
Thomas White committed
489
	int config_polar = 1;
490
	int config_sanity = 0;
491
	int config_satcorr = 1;
492
	int config_sa = 1;
493
	int config_checkprefix = 1;
494
	int config_closer = 1;
495
	float threshold = 800.0;
Thomas White's avatar
Thomas White committed
496
497
	struct detector *det;
	char *geometry = NULL;
Thomas White's avatar
Thomas White committed
498
499
	IndexingMethod indm;
	char *indm_str = NULL;
500
501
502
	UnitCell *cell;
	double *intensities = NULL;
	char *intfile = NULL;
Thomas White's avatar
Thomas White committed
503
	char *pdb = NULL;
504
	char *prefix = NULL;
505
506
	char *speaks = NULL;
	int peaks;
507
508
509
	int nthreads = 1;
	pthread_t workers[MAX_THREADS];
	struct process_args *worker_args[MAX_THREADS];
Thomas White's avatar
Thomas White committed
510
	int worker_active[MAX_THREADS];
511
	int i;
512
	pthread_mutex_t output_mutex = PTHREAD_MUTEX_INITIALIZER;
513
	pthread_mutex_t gpu_mutex = PTHREAD_MUTEX_INITIALIZER;
514
515
516
	char prepare_line[1024];
	char prepare_filename[1024];
	IndexingPrivate *ipriv;
Thomas White's avatar
Thomas White committed
517
518
519
520
521

	/* Long options */
	const struct option longopts[] = {
		{"help",               0, NULL,               'h'},
		{"input",              1, NULL,               'i'},
522
		{"output",             1, NULL,               'o'},
Thomas White's avatar
Thomas White committed
523
		{"gpu",                0, &config_gpu,         1},
524
		{"no-index",           0, &config_noindex,     1},
Thomas White's avatar
Thomas White committed
525
		{"dump-peaks",         0, &config_dumpfound,   1},
526
		{"peaks",              1, NULL,                2},
527
		{"near-bragg",         0, &config_nearbragg,   1},
Thomas White's avatar
Thomas White committed
528
529
		{"write-drx",          0, &config_writedrx,    1},
		{"indexing",           1, NULL,               'z'},
Thomas White's avatar
Thomas White committed
530
		{"geometry",           1, NULL,               'g'},
531
		{"simulate",           0, &config_simulate,    1},
532
533
		{"filter-cm",          0, &config_cmfilter,    1},
		{"filter-noise",       0, &config_noisefilter, 1},
534
		{"no-match",           0, &config_nomatch,     1},
535
		{"verbose",            0, &config_verbose,     1},
536
		{"alternate",          0, &config_alternate,   1},
537
		{"intensities",        1, NULL,               'q'},
Thomas White's avatar
Thomas White committed
538
		{"pdb",                1, NULL,               'p'},
539
		{"prefix",             1, NULL,               'x'},
Thomas White's avatar
Thomas White committed
540
		{"unpolarized",        0, &config_polar,       0},
541
		{"check-sanity",       0, &config_sanity,      1},
542
543
		{"no-sat-corr",        0, &config_satcorr,     0},
		{"sat-corr",           0, &config_satcorr,     1}, /* Compat */
544
		{"no-sa",              0, &config_sa,          0},
545
		{"threshold",          1, NULL,               't'},
546
		{"no-check-prefix",    0, &config_checkprefix, 0},
547
		{"no-closer-peak",     0, &config_closer,      0},
Thomas White's avatar
Thomas White committed
548
549
550
551
		{0, 0, NULL, 0}
	};

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

		switch (c) {
Thomas White's avatar
Thomas White committed
556
		case 'h' :
Thomas White's avatar
Thomas White committed
557
558
559
			show_help(argv[0]);
			return 0;

Thomas White's avatar
Thomas White committed
560
		case 'i' :
Thomas White's avatar
Thomas White committed
561
562
563
			filename = strdup(optarg);
			break;

564
565
566
567
		case 'o' :
			outfile = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
568
		case 'z' :
Thomas White's avatar
Thomas White committed
569
570
571
			indm_str = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
572
		case 'q' :
573
574
575
			intfile = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
576
		case 'p' :
Thomas White's avatar
Thomas White committed
577
578
579
			pdb = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
580
		case 'x' :
581
582
583
			prefix = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
584
		case 'j' :
585
586
587
			nthreads = atoi(optarg);
			break;

Thomas White's avatar
Thomas White committed
588
		case 'g' :
Thomas White's avatar
Thomas White committed
589
590
591
			geometry = strdup(optarg);
			break;

592
593
594
595
		case 't' :
			threshold = strtof(optarg, NULL);
			break;

596
597
598
599
		case 2 :
			speaks = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
600
		case 0 :
Thomas White's avatar
Thomas White committed
601
602
			break;

Thomas White's avatar
Thomas White committed
603
		default :
Thomas White's avatar
Thomas White committed
604
605
606
607
608
609
610
611
612
613
614
615
616
617
			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
618
		ERROR("Failed to open input file '%s'\n", filename);
Thomas White's avatar
Thomas White committed
619
620
		return 1;
	}
Thomas White's avatar
Thomas White committed
621
	free(filename);
Thomas White's avatar
Thomas White committed
622

623
624
625
626
627
628
629
630
631
632
633
634
635
636
	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);

637
638
639
640
641
642
643
644
645
646
647
648
649
650
	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;
	}

651
	if ( intfile != NULL ) {
Thomas White's avatar
Thomas White committed
652
		ReflItemList *items;
653
654
		items = read_reflections(intfile, intensities,
		                         NULL, NULL, NULL);
Thomas White's avatar
Thomas White committed
655
		delete_items(items);
656
657
658
659
	} else {
		intensities = NULL;
	}

Thomas White's avatar
Thomas White committed
660
661
662
663
	if ( pdb == NULL ) {
		pdb = strdup("molecule.pdb");
	}

664
	if ( prefix == NULL ) {
Thomas White's avatar
Thomas White committed
665
		prefix = strdup("");
666
	} else {
667
668
669
		if ( config_checkprefix ) {
			prefix = check_prefix(prefix);
		}
670
671
	}

672
	if ( (nthreads == 0) || (nthreads > MAX_THREADS) ) {
673
674
675
676
		ERROR("Invalid number of threads.\n");
		return 1;
	}

Thomas White's avatar
Thomas White committed
677
678
	if ( indm_str == NULL ) {
		STATUS("You didn't specify an indexing method, so I won't"
679
680
681
		       " 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
682
683
684
685
686
		indm = INDEXING_NONE;
	} else if ( strcmp(indm_str, "none") == 0 ) {
		indm = INDEXING_NONE;
	} else if ( strcmp(indm_str, "dirax") == 0) {
		indm = INDEXING_DIRAX;
687
688
	} else if ( strcmp(indm_str, "template") == 0) {
		indm = INDEXING_TEMPLATE;
Thomas White's avatar
Thomas White committed
689
690
691
692
	} else {
		ERROR("Unrecognised indexing method '%s'\n", indm_str);
		return 1;
	}
693
	free(indm_str);
Thomas White's avatar
Thomas White committed
694

Thomas White's avatar
Thomas White committed
695
696
697
698
699
700
701
702
703
704
705
706
	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);

707
	if ( (!config_nomatch) || (indm == INDEXING_TEMPLATE) ) {
708
709
		cell = load_cell_from_pdb(pdb);
		if ( cell == NULL ) {
710
711
712
			ERROR("Couldn't read unit cell (from %s)\n", pdb);
			return 1;
		}
713
714
715
	} else {
		STATUS("No cell needed because --no-match was used.\n");
		cell = NULL;
716
	}
717
	free(pdb);
718

719
720
721
722
723
	/* Start by writing the entire command line to stdout */
	fprintf(ofh, "Command line:");
	for ( i=0; i<argc; i++ ) {
		fprintf(ofh, " %s", argv[i]);
	}
724
725
	fprintf(ofh, "\n");
	fflush(ofh);
726

727
728
729
730
731
732
733
734
	/* Get first filename and use it to set up the indexing */
	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);
735
	ipriv = prepare_indexing(indm, cell, prepare_filename, det);
736
737
738
739
740
	if ( ipriv == NULL ) {
		ERROR("Failed to prepare indexing.\n");
		return 1;
	}

Thomas White's avatar
Thomas White committed
741
	gsl_set_error_handler_off();
Thomas White's avatar
Thomas White committed
742
	n_images = 0;
Thomas White's avatar
Thomas White committed
743
	n_indexable = 0;
Thomas White's avatar
Thomas White committed
744
	n_sane = 0;
745

746
747
748
	for ( i=0; i<nthreads; i++ ) {
		worker_args[i] = malloc(sizeof(struct process_args));
		worker_args[i]->filename = malloc(1024);
749
		worker_args[i]->ofh = ofh;
750
		worker_args[i]->peaks = peaks;
751
752
753
		worker_active[i] = 0;
	}

754
	/* Start threads off */
755
	for ( i=0; i<nthreads; i++ ) {
Thomas White's avatar
Thomas White committed
756
757

		char line[1024];
758
759
		struct process_args *pargs;
		int r;
Thomas White's avatar
Thomas White committed
760

761
		pargs = worker_args[i];
Thomas White's avatar
Thomas White committed
762

763
764
765
766
767
768
769
		if ( strlen(prepare_line) > 0 ) {
			strcpy(line, prepare_line);
			prepare_line[0] = '\0';
		} else {
			rval = fgets(line, 1023, fh);
			if ( rval == NULL ) continue;
		}
Thomas White's avatar
Thomas White committed
770
		chomp(line);
771
		snprintf(pargs->filename, 1023, "%s%s", prefix, line);
772

773
774
		n_images++;

775
		pargs->output_mutex = &output_mutex;
776
		pargs->gpu_mutex = &gpu_mutex;
777
		pthread_mutex_init(&pargs->control_mutex, NULL);
778
779
780
781
782
783
784
785
786
787
		pargs->config_cmfilter = config_cmfilter;
		pargs->config_noisefilter = config_noisefilter;
		pargs->config_writedrx = config_writedrx;
		pargs->config_dumpfound = config_dumpfound;
		pargs->config_verbose = config_verbose;
		pargs->config_alternate = config_alternate;
		pargs->config_nearbragg = config_nearbragg;
		pargs->config_gpu = config_gpu;
		pargs->config_simulate = config_simulate;
		pargs->config_nomatch = config_nomatch;
Thomas White's avatar
Thomas White committed
788
		pargs->config_polar = config_polar;
789
		pargs->config_sanity = config_sanity;
790
		pargs->config_satcorr = config_satcorr;
791
		pargs->config_sa = config_sa;
792
		pargs->config_closer = config_closer;
793
		pargs->cell = cell;
Thomas White's avatar
Thomas White committed
794
		pargs->det = det;
795
		pargs->ipriv = ipriv;
796
797
798
		pargs->indm = indm;
		pargs->intensities = intensities;
		pargs->gctx = gctx;
799
		pargs->threshold = threshold;
Thomas White's avatar
Thomas White committed
800
		pargs->id = i;
801
802
803
804
805
		pthread_mutex_lock(&pargs->control_mutex);
		pargs->done = 0;
		pargs->start = 1;
		pargs->finish = 0;
		pthread_mutex_unlock(&pargs->control_mutex);
806

Thomas White's avatar
Thomas White committed
807
		worker_active[i] = 1;
808
		r = pthread_create(&workers[i], NULL, worker_thread, pargs);
809
		if ( r != 0 ) {
Thomas White's avatar
Thomas White committed
810
			worker_active[i] = 0;
811
812
813
814
815
			ERROR("Couldn't start thread %i\n", i);
		}

	}

816
	/* Keep threads busy until the end of the data */
817
818
819
820
821
822
823
824
	do {

		int i;

		for ( i=0; i<nthreads; i++ ) {

			char line[1024];
			struct process_args *pargs;
825
826
827
828
			int done;

			/* Spend CPU time indexing, not checking results */
			usleep(100000);
829

830
			/* Are we using this thread record at all? */
Thomas White's avatar
Thomas White committed
831
832
			if ( !worker_active[i] ) continue;

833
			/* Has the thread finished yet? */
834
			pargs = worker_args[i];
835
836
837
838
			pthread_mutex_lock(&pargs->control_mutex);
			done = pargs->done;
			pthread_mutex_unlock(&pargs->control_mutex);
			if ( !done ) continue;
839

840
841
			/* Results will be processed after checking if
			 * there are any more images to process. */
842

843
			/* Get next filename */
844
			rval = fgets(line, 1023, fh);
845
846
847
			/* In this case, the result of the last file
			 * file will be processed when the thread is
			 * joined. */
848
			if ( rval == NULL ) break;
849
850

			/* Record the result */
Thomas White's avatar
Thomas White committed
851
852
			n_indexable += pargs->indexable;
			n_sane += pargs->sane;
853

854
			chomp(line);
855
			snprintf(pargs->filename, 1023, "%s%s", prefix, line);
856
857

			n_images++;
858
859
860
861
862
863
864

			/* Wake the thread up ... */
			pthread_mutex_lock(&pargs->control_mutex);
			pargs->done = 0;
			pargs->start = 1;
			pthread_mutex_unlock(&pargs->control_mutex);

865
		}
Thomas White's avatar
Thomas White committed
866

Thomas White's avatar
Thomas White committed
867
868
	} while ( rval != NULL );

869
	/* Join threads */
870
	for ( i=0; i<nthreads; i++ ) {
Thomas White's avatar
Thomas White committed
871

Thomas White's avatar
Thomas White committed
872
		if ( !worker_active[i] ) goto free;
Thomas White's avatar
Thomas White committed
873

874
875
		/* Tell the thread to exit */
		struct process_args *pargs = worker_args[i];
Thomas White's avatar
Thomas White committed
876
		pthread_mutex_lock(&pargs->control_mutex);
877
		pargs->finish = 1;
Thomas White's avatar
Thomas White committed
878
		pthread_mutex_unlock(&pargs->control_mutex);
Thomas White's avatar
Thomas White committed
879

880
881
		/* Wait for it to join */
		pthread_join(workers[i], NULL);
Thomas White's avatar
Thomas White committed
882
883
		worker_active[i] = 0;

Thomas White's avatar
Thomas White committed
884
885
		n_indexable += pargs->indexable;
		n_sane += pargs->sane;
Thomas White's avatar
Thomas White committed
886

Thomas White's avatar
Thomas White committed
887
	free:
888
		if ( worker_args[i]->filename != NULL ) {
889
890
			free(worker_args[i]->filename);
		}
891
		free(worker_args[i]);
Thomas White's avatar
Thomas White committed
892

893
894
	}

895
896
	cleanup_indexing(ipriv);

897
	free(prefix);
Thomas White's avatar
Thomas White committed
898
899
	free(det->panels);
	free(det);
900
	cell_free(cell);
Thomas White's avatar
Thomas White committed
901
902
	fclose(fh);

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

Thomas White's avatar
Thomas White committed
906
907
908
909
	if ( gctx != NULL ) {
		cleanup_gpu(gctx);
	}

Thomas White's avatar
Thomas White committed
910
911
	return 0;
}