indexamajig.c 22.8 KB
Newer Older
Thomas White's avatar
Thomas White committed
1
/*
2
 * indexamajig.c
Thomas White's avatar
Thomas White committed
3
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
42
43
#define MAX_THREADS (96)

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

	/* Thread control and output */
	pthread_mutex_t control_mutex;  /* Protects the scary stuff below */
	int start;
	int finish;
	int done;
	int hit;
	int peaks_sane;
79
80
};

81

82
83
84
struct process_result
{
	int hit;
Thomas White's avatar
Thomas White committed
85
	int peaks_sane;
86
87
88
};


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


172
static struct image *get_simage(struct image *template, int alternate)
173
{
174
	struct image *image;
175
	struct panel panels[2];
176

177
178
	image = malloc(sizeof(*image));

179
	/* Simulate a diffraction pattern */
180
181
182
	image->twotheta = NULL;
	image->data = NULL;
	image->det = template->det;
183
184
185
	image->flags = NULL;
	image->f0_available = 0;
	image->f0 = 1.0;
186
187

	/* View head-on (unit cell is tilted) */
188
189
190
191
	image->orientation.w = 1.0;
	image->orientation.x = 0.0;
	image->orientation.y = 0.0;
	image->orientation.z = 0.0;
192
193
194

	/* Detector geometry for the simulation
	 * - not necessarily the same as the original. */
195
196
	image->width = 1024;
	image->height = 1024;
Thomas White's avatar
Thomas White committed
197
	image->det->n_panels = 2;
198

199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
	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;
217
		panels[1].cy = 525.0;
218
219
220
		panels[1].clen = 56.7e-2;  /* 56.7 cm */
		panels[1].res = 13333.3;   /* 75 microns/pixel */

Thomas White's avatar
Thomas White committed
221
		image->det->panels = panels;
222

223
224
225
	} else {

		/* Copy pointer to old geometry */
Thomas White's avatar
Thomas White committed
226
		image->det->panels = template->det->panels;
227
228

	}
229

230
	image->lambda = ph_en_to_lambda(eV_to_J(1.8e3));
231
	image->features = template->features;
232
	image->filename = template->filename;
233
	image->indexed_cell = template->indexed_cell;
Thomas White's avatar
Thomas White committed
234
	image->f0 = template->f0;
235

Thomas White's avatar
Thomas White committed
236
237
238
239
	/* Prevent muppetry */
	image->hits = NULL;
	image->n_hits = 0;

240
241
	return image;
}
242
243


244
static void simulate_and_write(struct image *simage, struct gpu_context **gctx,
Thomas White's avatar
Thomas White committed
245
                               const double *intensities, UnitCell *cell)
246
{
247
	/* Set up GPU if necessary */
Thomas White's avatar
Thomas White committed
248
	if ( (gctx != NULL) && (*gctx == NULL) ) {
Thomas White's avatar
Thomas White committed
249
		*gctx = setup_gpu(0, simage, intensities);
250
251
	}

Thomas White's avatar
Thomas White committed
252
	if ( (gctx != NULL) && (*gctx != NULL) ) {
253
		get_diffraction_gpu(*gctx, simage, 24, 24, 40, cell);
254
	} else {
255
		get_diffraction(simage, 24, 24, 40,
Thomas White's avatar
Thomas White committed
256
		                intensities, NULL, cell, 0,
257
		                GRADIENT_MOSAIC);
258
	}
259
	record_image(simage, 0);
260

261
	hdf5_write("simulated.h5", simage->data, simage->width, simage->height,
262
263
264
265
		   H5T_NATIVE_FLOAT);
}


266
static struct process_result process_image(struct process_args *pargs)
267
268
269
270
271
{
	struct hdfile *hdfile;
	struct image image;
	struct image *simage;
	float *data_for_measurement;
272
	struct process_result result;
273
	size_t data_size;
274
275
276
277
278
279
280
281
282
283
284
285
	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
286
	int config_polar = pargs->config_polar;
287
288
289
	IndexingMethod indm = pargs->indm;
	const double *intensities = pargs->intensities;
	struct gpu_context *gctx = pargs->gctx;
290
291
292
293

	image.features = NULL;
	image.data = NULL;
	image.indexed_cell = NULL;
Thomas White's avatar
Thomas White committed
294
	image.id = pargs->id;
295
	image.filename = filename;
Thomas White's avatar
Thomas White committed
296
297
	image.hits = NULL;
	image.n_hits = 0;
Thomas White's avatar
Thomas White committed
298
	image.det = pargs->det;
299

300
301
302
303
304
305
	/* 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;

306
	STATUS("Processing '%s'\n", image.filename);
307

308
309
	result.peaks_sane = 0;
	result.hit = 0;
310

311
312
	hdfile = hdfile_open(filename);
	if ( hdfile == NULL ) {
313
		return result;
314
315
	} else if ( hdfile_set_first_image(hdfile, "/") ) {
		ERROR("Couldn't select path\n");
316
		return result;
317
318
	}

319
	hdf5_read(hdfile, &image, pargs->config_satcorr);
320
321
322
323
324
325
326
327
328
329
330
331
332

	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 {
333
		memcpy(data_for_measurement, image.data, data_size);
334
335
336
	}

	/* Perform 'fine' peak search */
337
	search_peaks(&image, pargs->threshold);
Thomas White's avatar
Thomas White committed
338
339
340
341
342
343

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

344
345
	if ( image_feature_count(image.features) < 5 ) goto done;

346
	if ( config_dumpfound ) dump_peaks(&image, pargs->output_mutex);
347
348
349
350
351
352
353
354
355
356

	/* 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,
357
		              config_verbose, pargs->ipriv);
358
359
360
361
362
	}

	/* No cell at this point?  Then we're done. */
	if ( image.indexed_cell == NULL ) goto done;

Thomas White's avatar
Thomas White committed
363
	/* Sanity check */
364
	if ( pargs->config_sanity
365
	  && !peak_sanity_check(&image, image.indexed_cell, 0, 0.1) ) {
Thomas White's avatar
Thomas White committed
366
367
368
		STATUS("Failed peak sanity check.\n");
		goto done;
	} else {
369
		result.peaks_sane = 1;
Thomas White's avatar
Thomas White committed
370
371
	}

372
373
	/* Measure intensities if requested */
	if ( config_nearbragg ) {
374
		output_intensities(&image, image.indexed_cell,
Thomas White's avatar
Thomas White committed
375
		                   pargs->output_mutex, config_polar,
376
377
		                   pargs->config_sa, pargs->config_closer,
		                   0, 0.1);
378
379
	}

380
381
	simage = get_simage(&image, config_alternate);

382
383
384
	/* Simulate if requested */
	if ( config_simulate ) {
		if ( config_gpu ) {
385
			pthread_mutex_lock(pargs->gpu_mutex);
386
			simulate_and_write(simage, &gctx, intensities,
Thomas White's avatar
Thomas White committed
387
			                   image.indexed_cell);
388
			pthread_mutex_unlock(pargs->gpu_mutex);
389
390
		} else {
			simulate_and_write(simage, NULL, intensities,
Thomas White's avatar
Thomas White committed
391
			                   image.indexed_cell);
392
393
394
395
396
397
398
399
400
		}
	}

	/* 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 */
401
	cell_free(image.indexed_cell);
402
403
404

done:
	free(image.data);
Thomas White's avatar
Thomas White committed
405
	free(image.flags);
406
	image_feature_list_free(image.features);
Thomas White's avatar
Thomas White committed
407
	free(image.hits);
408
409
	hdfile_close(hdfile);

410
	if ( image.indexed_cell == NULL ) {
411
		result.hit = 0;
412
	} else {
413
		result.hit = 1;
414
415
	}
	return result;
416
417
418
}


419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
static void *worker_thread(void *pargsv)
{
	struct process_args *pargs = pargsv;
	int finish;

	do {

		struct process_result result;
		int wakeup;

		result = process_image(pargs);

		pthread_mutex_lock(&pargs->control_mutex);
		pargs->hit = result.hit;
		pargs->peaks_sane = result.peaks_sane;
		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
456
457
458
int main(int argc, char *argv[])
{
	int c;
Thomas White's avatar
Thomas White committed
459
	struct gpu_context *gctx = NULL;
Thomas White's avatar
Thomas White committed
460
461
	char *filename = NULL;
	FILE *fh;
462
	char *rval = NULL;
Thomas White's avatar
Thomas White committed
463
464
	int n_images;
	int n_hits;
Thomas White's avatar
Thomas White committed
465
	int n_sane;
466
	int config_noindex = 0;
Thomas White's avatar
Thomas White committed
467
	int config_dumpfound = 0;
468
	int config_nearbragg = 0;
Thomas White's avatar
Thomas White committed
469
	int config_writedrx = 0;
470
	int config_simulate = 0;
471
472
	int config_cmfilter = 0;
	int config_noisefilter = 0;
473
	int config_nomatch = 0;
Thomas White's avatar
Thomas White committed
474
	int config_gpu = 0;
475
	int config_verbose = 0;
476
	int config_alternate = 0;
Thomas White's avatar
Thomas White committed
477
	int config_polar = 1;
478
	int config_sanity = 0;
479
	int config_satcorr = 1;
480
	int config_sa = 1;
481
	int config_checkprefix = 1;
482
	int config_closer = 1;
483
	float threshold = 800.0;
Thomas White's avatar
Thomas White committed
484
485
	struct detector *det;
	char *geometry = NULL;
Thomas White's avatar
Thomas White committed
486
487
	IndexingMethod indm;
	char *indm_str = NULL;
488
489
490
	UnitCell *cell;
	double *intensities = NULL;
	char *intfile = NULL;
Thomas White's avatar
Thomas White committed
491
	char *pdb = NULL;
492
	char *prefix = NULL;
493
494
495
	int nthreads = 1;
	pthread_t workers[MAX_THREADS];
	struct process_args *worker_args[MAX_THREADS];
Thomas White's avatar
Thomas White committed
496
	int worker_active[MAX_THREADS];
497
	int i;
498
	pthread_mutex_t output_mutex = PTHREAD_MUTEX_INITIALIZER;
499
	pthread_mutex_t gpu_mutex = PTHREAD_MUTEX_INITIALIZER;
500
501
502
	char prepare_line[1024];
	char prepare_filename[1024];
	IndexingPrivate *ipriv;
Thomas White's avatar
Thomas White committed
503
504
505
506
507

	/* Long options */
	const struct option longopts[] = {
		{"help",               0, NULL,               'h'},
		{"input",              1, NULL,               'i'},
Thomas White's avatar
Thomas White committed
508
		{"gpu",                0, &config_gpu,         1},
509
		{"no-index",           0, &config_noindex,     1},
Thomas White's avatar
Thomas White committed
510
		{"dump-peaks",         0, &config_dumpfound,   1},
511
		{"near-bragg",         0, &config_nearbragg,   1},
Thomas White's avatar
Thomas White committed
512
513
		{"write-drx",          0, &config_writedrx,    1},
		{"indexing",           1, NULL,               'z'},
Thomas White's avatar
Thomas White committed
514
		{"geometry",           1, NULL,               'g'},
515
		{"simulate",           0, &config_simulate,    1},
516
517
		{"filter-cm",          0, &config_cmfilter,    1},
		{"filter-noise",       0, &config_noisefilter, 1},
518
		{"no-match",           0, &config_nomatch,     1},
519
		{"verbose",            0, &config_verbose,     1},
520
		{"alternate",          0, &config_alternate,   1},
521
		{"intensities",        1, NULL,               'q'},
Thomas White's avatar
Thomas White committed
522
		{"pdb",                1, NULL,               'p'},
523
		{"prefix",             1, NULL,               'x'},
Thomas White's avatar
Thomas White committed
524
		{"unpolarized",        0, &config_polar,       0},
525
		{"check-sanity",       0, &config_sanity,      1},
526
527
		{"no-sat-corr",        0, &config_satcorr,     0},
		{"sat-corr",           0, &config_satcorr,     1}, /* Compat */
528
		{"no-sa",              0, &config_sa,          0},
529
		{"threshold",          1, NULL,               't'},
530
		{"no-check-prefix",    0, &config_checkprefix, 0},
531
		{"no-closer-peak",     0, &config_closer,      0},
Thomas White's avatar
Thomas White committed
532
533
534
535
		{0, 0, NULL, 0}
	};

	/* Short options */
536
	while ((c = getopt_long(argc, argv, "hi:wp:j:x:g:t:",
Thomas White's avatar
Thomas White committed
537
	                        longopts, NULL)) != -1) {
Thomas White's avatar
Thomas White committed
538
539

		switch (c) {
Thomas White's avatar
Thomas White committed
540
		case 'h' :
Thomas White's avatar
Thomas White committed
541
542
543
			show_help(argv[0]);
			return 0;

Thomas White's avatar
Thomas White committed
544
		case 'i' :
Thomas White's avatar
Thomas White committed
545
546
547
			filename = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
548
		case 'z' :
Thomas White's avatar
Thomas White committed
549
550
551
			indm_str = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
552
		case 'q' :
553
554
555
			intfile = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
556
		case 'p' :
Thomas White's avatar
Thomas White committed
557
558
559
			pdb = strdup(optarg);
			break;

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

Thomas White's avatar
Thomas White committed
564
		case 'j' :
565
566
567
			nthreads = atoi(optarg);
			break;

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

572
573
574
575
		case 't' :
			threshold = strtof(optarg, NULL);
			break;

Thomas White's avatar
Thomas White committed
576
		case 0 :
Thomas White's avatar
Thomas White committed
577
578
			break;

Thomas White's avatar
Thomas White committed
579
		default :
Thomas White's avatar
Thomas White committed
580
581
582
583
584
			return 1;
		}

	}

585
586
587
588
589
590
591
	/* Start by writing the entire command line to stdout */
	printf("Command line:");
	for ( i=0; i<argc; i++ ) {
		printf(" %s", argv[i]);
	}
	printf("\n");

Thomas White's avatar
Thomas White committed
592
593
594
595
596
597
598
599
600
	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
601
		ERROR("Failed to open input file '%s'\n", filename);
Thomas White's avatar
Thomas White committed
602
603
		return 1;
	}
Thomas White's avatar
Thomas White committed
604
	free(filename);
Thomas White's avatar
Thomas White committed
605

606
	if ( intfile != NULL ) {
Thomas White's avatar
Thomas White committed
607
		ReflItemList *items;
608
609
		items = read_reflections(intfile, intensities,
		                         NULL, NULL, NULL);
Thomas White's avatar
Thomas White committed
610
		delete_items(items);
611
612
613
614
	} else {
		intensities = NULL;
	}

Thomas White's avatar
Thomas White committed
615
616
617
618
	if ( pdb == NULL ) {
		pdb = strdup("molecule.pdb");
	}

619
	if ( prefix == NULL ) {
Thomas White's avatar
Thomas White committed
620
		prefix = strdup("");
621
	} else {
622
623
624
		if ( config_checkprefix ) {
			prefix = check_prefix(prefix);
		}
625
626
	}

627
	if ( (nthreads == 0) || (nthreads > MAX_THREADS) ) {
628
629
630
631
		ERROR("Invalid number of threads.\n");
		return 1;
	}

Thomas White's avatar
Thomas White committed
632
633
	if ( indm_str == NULL ) {
		STATUS("You didn't specify an indexing method, so I won't"
634
635
636
		       " 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
637
638
639
640
641
		indm = INDEXING_NONE;
	} else if ( strcmp(indm_str, "none") == 0 ) {
		indm = INDEXING_NONE;
	} else if ( strcmp(indm_str, "dirax") == 0) {
		indm = INDEXING_DIRAX;
642
643
	} else if ( strcmp(indm_str, "template") == 0) {
		indm = INDEXING_TEMPLATE;
Thomas White's avatar
Thomas White committed
644
645
646
647
	} else {
		ERROR("Unrecognised indexing method '%s'\n", indm_str);
		return 1;
	}
648
	free(indm_str);
Thomas White's avatar
Thomas White committed
649

Thomas White's avatar
Thomas White committed
650
651
652
653
654
655
656
657
658
659
660
661
	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);

662
	if ( (!config_nomatch) || (indm == INDEXING_TEMPLATE) ) {
663
664
		cell = load_cell_from_pdb(pdb);
		if ( cell == NULL ) {
665
666
667
			ERROR("Couldn't read unit cell (from %s)\n", pdb);
			return 1;
		}
668
669
670
	} else {
		STATUS("No cell needed because --no-match was used.\n");
		cell = NULL;
671
	}
672
	free(pdb);
673

674
675
676
677
678
679
680
681
	/* 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);
682
	ipriv = prepare_indexing(indm, cell, prepare_filename, det);
683
684
685
686
687
	if ( ipriv == NULL ) {
		ERROR("Failed to prepare indexing.\n");
		return 1;
	}

Thomas White's avatar
Thomas White committed
688
	gsl_set_error_handler_off();
Thomas White's avatar
Thomas White committed
689
690
	n_images = 0;
	n_hits = 0;
Thomas White's avatar
Thomas White committed
691
	n_sane = 0;
692

693
694
695
696
697
698
	for ( i=0; i<nthreads; i++ ) {
		worker_args[i] = malloc(sizeof(struct process_args));
		worker_args[i]->filename = malloc(1024);
		worker_active[i] = 0;
	}

699
	/* Start threads off */
700
	for ( i=0; i<nthreads; i++ ) {
Thomas White's avatar
Thomas White committed
701
702

		char line[1024];
703
704
		struct process_args *pargs;
		int r;
Thomas White's avatar
Thomas White committed
705

706
		pargs = worker_args[i];
Thomas White's avatar
Thomas White committed
707

708
709
710
711
712
713
714
		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
715
		chomp(line);
716
		snprintf(pargs->filename, 1023, "%s%s", prefix, line);
717

718
719
		n_images++;

720
		pargs->output_mutex = &output_mutex;
721
		pargs->gpu_mutex = &gpu_mutex;
722
		pthread_mutex_init(&pargs->control_mutex, NULL);
723
724
725
726
727
728
729
730
731
732
		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
733
		pargs->config_polar = config_polar;
734
		pargs->config_sanity = config_sanity;
735
		pargs->config_satcorr = config_satcorr;
736
		pargs->config_sa = config_sa;
737
		pargs->config_closer = config_closer;
738
		pargs->cell = cell;
Thomas White's avatar
Thomas White committed
739
		pargs->det = det;
740
		pargs->ipriv = ipriv;
741
742
743
		pargs->indm = indm;
		pargs->intensities = intensities;
		pargs->gctx = gctx;
744
		pargs->threshold = threshold;
Thomas White's avatar
Thomas White committed
745
		pargs->id = i;
746
747
748
749
750
		pthread_mutex_lock(&pargs->control_mutex);
		pargs->done = 0;
		pargs->start = 1;
		pargs->finish = 0;
		pthread_mutex_unlock(&pargs->control_mutex);
751

Thomas White's avatar
Thomas White committed
752
		worker_active[i] = 1;
753
		r = pthread_create(&workers[i], NULL, worker_thread, pargs);
754
		if ( r != 0 ) {
Thomas White's avatar
Thomas White committed
755
			worker_active[i] = 0;
756
757
758
759
760
			ERROR("Couldn't start thread %i\n", i);
		}

	}

761
	/* Keep threads busy until the end of the data */
762
763
764
765
766
767
768
769
	do {

		int i;

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

			char line[1024];
			struct process_args *pargs;
770
771
772
773
			int done;

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

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

778
			/* Has the thread finished yet? */
779
			pargs = worker_args[i];
780
781
782
783
			pthread_mutex_lock(&pargs->control_mutex);
			done = pargs->done;
			pthread_mutex_unlock(&pargs->control_mutex);
			if ( !done ) continue;
784

785
786
			/* Results will be processed after checking if
			 * there are any more images to process. */
787

788
			/* Get next filename */
789
			rval = fgets(line, 1023, fh);
790
791
792
			/* In this case, the result of the last file
			 * file will be processed when the thread is
			 * joined. */
793
			if ( rval == NULL ) break;
794
795
796
797
798

			/* Record the result */
			n_hits += pargs->hit;
			n_sane += pargs->peaks_sane;

799
			chomp(line);
800
			snprintf(pargs->filename, 1023, "%s%s", prefix, line);
801
802

			n_images++;
803
804
805
806
807
808
809

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

810
		}
Thomas White's avatar
Thomas White committed
811

Thomas White's avatar
Thomas White committed
812
813
	} while ( rval != NULL );

814
	/* Join threads */
815
	for ( i=0; i<nthreads; i++ ) {
Thomas White's avatar
Thomas White committed
816

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

819
820
		/* Tell the thread to exit */
		struct process_args *pargs = worker_args[i];
Thomas White's avatar
Thomas White committed
821
		pthread_mutex_lock(&pargs->control_mutex);
822
		pargs->finish = 1;
Thomas White's avatar
Thomas White committed
823
		pthread_mutex_unlock(&pargs->control_mutex);
Thomas White's avatar
Thomas White committed
824

825
826
		/* Wait for it to join */
		pthread_join(workers[i], NULL);
Thomas White's avatar
Thomas White committed
827
828
		worker_active[i] = 0;

829
830
		n_hits += pargs->hit;
		n_sane += pargs->peaks_sane;
Thomas White's avatar
Thomas White committed
831

Thomas White's avatar
Thomas White committed
832
	free:
833
		if ( worker_args[i]->filename != NULL ) {
834
835
			free(worker_args[i]->filename);
		}
836
		free(worker_args[i]);
Thomas White's avatar
Thomas White committed
837

838
839
	}

840
841
	cleanup_indexing(ipriv);

842
	free(prefix);
Thomas White's avatar
Thomas White committed
843
844
	free(det->panels);
	free(det);
845
	cell_free(cell);
Thomas White's avatar
Thomas White committed
846
847
848
	fclose(fh);

	STATUS("There were %i images.\n", n_images);
Thomas White's avatar
Thomas White committed
849
	STATUS("%i hits were found, of which %i were sane.\n", n_hits, n_sane);
Thomas White's avatar
Thomas White committed
850

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

Thomas White's avatar
Thomas White committed
855
856
	return 0;
}