indexamajig.c 21.6 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;
Thomas White's avatar
Thomas White committed
64
	struct detector *det;
65
	IndexingMethod indm;
66
	IndexingPrivate *ipriv;
67
68
	const double *intensities;
	struct gpu_context *gctx;
69
70
71
72
73
74
75
76

	/* 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;
77
78
};

79

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


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


163
static struct image *get_simage(struct image *template, int alternate)
164
{
165
	struct image *image;
166
	struct panel panels[2];
167

168
169
	image = malloc(sizeof(*image));

170
	/* Simulate a diffraction pattern */
171
172
173
	image->twotheta = NULL;
	image->data = NULL;
	image->det = template->det;
174
175
176
	image->flags = NULL;
	image->f0_available = 0;
	image->f0 = 1.0;
177
178

	/* View head-on (unit cell is tilted) */
179
180
181
182
	image->orientation.w = 1.0;
	image->orientation.x = 0.0;
	image->orientation.y = 0.0;
	image->orientation.z = 0.0;
183
184
185

	/* Detector geometry for the simulation
	 * - not necessarily the same as the original. */
186
187
	image->width = 1024;
	image->height = 1024;
Thomas White's avatar
Thomas White committed
188
	image->det->n_panels = 2;
189

190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
	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;
208
		panels[1].cy = 525.0;
209
210
211
		panels[1].clen = 56.7e-2;  /* 56.7 cm */
		panels[1].res = 13333.3;   /* 75 microns/pixel */

Thomas White's avatar
Thomas White committed
212
		image->det->panels = panels;
213

214
215
216
	} else {

		/* Copy pointer to old geometry */
Thomas White's avatar
Thomas White committed
217
		image->det->panels = template->det->panels;
218
219

	}
220

221
	image->lambda = ph_en_to_lambda(eV_to_J(1.8e3));
222
	image->features = template->features;
223
	image->filename = template->filename;
224
	image->indexed_cell = template->indexed_cell;
Thomas White's avatar
Thomas White committed
225
	image->f0 = template->f0;
226

Thomas White's avatar
Thomas White committed
227
228
229
230
	/* Prevent muppetry */
	image->hits = NULL;
	image->n_hits = 0;

231
232
	return image;
}
233
234


235
static void simulate_and_write(struct image *simage, struct gpu_context **gctx,
Thomas White's avatar
Thomas White committed
236
                               const double *intensities, UnitCell *cell)
237
{
238
	/* Set up GPU if necessary */
Thomas White's avatar
Thomas White committed
239
	if ( (gctx != NULL) && (*gctx == NULL) ) {
Thomas White's avatar
Thomas White committed
240
		*gctx = setup_gpu(0, simage, intensities);
241
242
	}

Thomas White's avatar
Thomas White committed
243
	if ( (gctx != NULL) && (*gctx != NULL) ) {
244
		get_diffraction_gpu(*gctx, simage, 24, 24, 40, cell);
245
	} else {
246
		get_diffraction(simage, 24, 24, 40,
Thomas White's avatar
Thomas White committed
247
		                intensities, NULL, cell, 0,
248
		                GRADIENT_MOSAIC);
249
	}
250
	record_image(simage, 0);
251

252
	hdf5_write("simulated.h5", simage->data, simage->width, simage->height,
253
254
255
256
		   H5T_NATIVE_FLOAT);
}


257
static struct process_result process_image(struct process_args *pargs)
258
259
260
261
262
{
	struct hdfile *hdfile;
	struct image image;
	struct image *simage;
	float *data_for_measurement;
263
	struct process_result result;
264
	size_t data_size;
265
266
267
268
269
270
271
272
273
274
275
276
	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
277
	int config_polar = pargs->config_polar;
278
279
280
	IndexingMethod indm = pargs->indm;
	const double *intensities = pargs->intensities;
	struct gpu_context *gctx = pargs->gctx;
281
282
283
284

	image.features = NULL;
	image.data = NULL;
	image.indexed_cell = NULL;
Thomas White's avatar
Thomas White committed
285
	image.id = pargs->id;
286
	image.filename = filename;
Thomas White's avatar
Thomas White committed
287
288
	image.hits = NULL;
	image.n_hits = 0;
Thomas White's avatar
Thomas White committed
289
	image.det = pargs->det;
290

291
292
293
294
295
296
	/* 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;

297
	STATUS("Processing '%s'\n", image.filename);
298

299
300
	result.peaks_sane = 0;
	result.hit = 0;
301

302
303
	hdfile = hdfile_open(filename);
	if ( hdfile == NULL ) {
304
		return result;
305
306
	} else if ( hdfile_set_first_image(hdfile, "/") ) {
		ERROR("Couldn't select path\n");
307
		return result;
308
309
	}

310
	hdf5_read(hdfile, &image, pargs->config_satcorr);
311
312
313
314
315
316
317
318
319
320
321
322
323

	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 {
324
		memcpy(data_for_measurement, image.data, data_size);
325
326
327
328
	}

	/* Perform 'fine' peak search */
	search_peaks(&image);
Thomas White's avatar
Thomas White committed
329
330
331
332
333
334

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

335
336
	if ( image_feature_count(image.features) < 5 ) goto done;

337
	if ( config_dumpfound ) dump_peaks(&image, pargs->output_mutex);
338
339
340
341
342
343
344
345
346
347

	/* 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,
348
		              config_verbose, pargs->ipriv);
349
350
351
352
353
	}

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

Thomas White's avatar
Thomas White committed
354
	/* Sanity check */
355
356
	if ( pargs->config_sanity
	  && !peak_sanity_check(&image, image.indexed_cell) ) {
Thomas White's avatar
Thomas White committed
357
358
359
		STATUS("Failed peak sanity check.\n");
		goto done;
	} else {
360
		result.peaks_sane = 1;
Thomas White's avatar
Thomas White committed
361
362
	}

363
364
	/* Measure intensities if requested */
	if ( config_nearbragg ) {
365
		output_intensities(&image, image.indexed_cell,
Thomas White's avatar
Thomas White committed
366
		                   pargs->output_mutex, config_polar,
367
		                   pargs->config_sa);
368
369
	}

370
371
	simage = get_simage(&image, config_alternate);

372
373
374
	/* Simulate if requested */
	if ( config_simulate ) {
		if ( config_gpu ) {
375
			pthread_mutex_lock(pargs->gpu_mutex);
376
			simulate_and_write(simage, &gctx, intensities,
Thomas White's avatar
Thomas White committed
377
			                   image.indexed_cell);
378
			pthread_mutex_unlock(pargs->gpu_mutex);
379
380
		} else {
			simulate_and_write(simage, NULL, intensities,
Thomas White's avatar
Thomas White committed
381
			                   image.indexed_cell);
382
383
384
385
386
387
388
389
390
391
392
393
394
		}
	}

	/* 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 */
	free(image.indexed_cell);

done:
	free(image.data);
Thomas White's avatar
Thomas White committed
395
	free(image.flags);
396
	image_feature_list_free(image.features);
Thomas White's avatar
Thomas White committed
397
	free(image.hits);
398
399
	hdfile_close(hdfile);

400
	if ( image.indexed_cell == NULL ) {
401
		result.hit = 0;
402
	} else {
403
		result.hit = 1;
404
405
	}
	return result;
406
407
408
}


409
410
411
412
413
414
415
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
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
446
447
448
int main(int argc, char *argv[])
{
	int c;
Thomas White's avatar
Thomas White committed
449
	struct gpu_context *gctx = NULL;
Thomas White's avatar
Thomas White committed
450
451
	char *filename = NULL;
	FILE *fh;
452
	char *rval = NULL;
Thomas White's avatar
Thomas White committed
453
454
	int n_images;
	int n_hits;
Thomas White's avatar
Thomas White committed
455
	int n_sane;
456
	int config_noindex = 0;
Thomas White's avatar
Thomas White committed
457
	int config_dumpfound = 0;
458
	int config_nearbragg = 0;
Thomas White's avatar
Thomas White committed
459
	int config_writedrx = 0;
460
	int config_simulate = 0;
461
462
	int config_cmfilter = 0;
	int config_noisefilter = 0;
463
	int config_nomatch = 0;
Thomas White's avatar
Thomas White committed
464
	int config_gpu = 0;
465
	int config_verbose = 0;
466
	int config_alternate = 0;
Thomas White's avatar
Thomas White committed
467
	int config_polar = 1;
468
	int config_sanity = 0;
469
	int config_satcorr = 0;
470
	int config_sa = 1;
Thomas White's avatar
Thomas White committed
471
472
	struct detector *det;
	char *geometry = NULL;
Thomas White's avatar
Thomas White committed
473
474
	IndexingMethod indm;
	char *indm_str = NULL;
475
476
477
	UnitCell *cell;
	double *intensities = NULL;
	char *intfile = NULL;
Thomas White's avatar
Thomas White committed
478
	char *pdb = NULL;
479
	char *prefix = NULL;
480
481
482
	int nthreads = 1;
	pthread_t workers[MAX_THREADS];
	struct process_args *worker_args[MAX_THREADS];
Thomas White's avatar
Thomas White committed
483
	int worker_active[MAX_THREADS];
484
	int i;
485
	pthread_mutex_t output_mutex = PTHREAD_MUTEX_INITIALIZER;
486
	pthread_mutex_t gpu_mutex = PTHREAD_MUTEX_INITIALIZER;
487
488
489
	char prepare_line[1024];
	char prepare_filename[1024];
	IndexingPrivate *ipriv;
Thomas White's avatar
Thomas White committed
490
491
492
493
494

	/* Long options */
	const struct option longopts[] = {
		{"help",               0, NULL,               'h'},
		{"input",              1, NULL,               'i'},
Thomas White's avatar
Thomas White committed
495
		{"gpu",                0, &config_gpu,         1},
496
		{"no-index",           0, &config_noindex,     1},
Thomas White's avatar
Thomas White committed
497
		{"dump-peaks",         0, &config_dumpfound,   1},
498
		{"near-bragg",         0, &config_nearbragg,   1},
Thomas White's avatar
Thomas White committed
499
500
		{"write-drx",          0, &config_writedrx,    1},
		{"indexing",           1, NULL,               'z'},
Thomas White's avatar
Thomas White committed
501
		{"geometry",           1, NULL,               'g'},
502
		{"simulate",           0, &config_simulate,    1},
503
504
		{"filter-cm",          0, &config_cmfilter,    1},
		{"filter-noise",       0, &config_noisefilter, 1},
505
		{"no-match",           0, &config_nomatch,     1},
506
		{"verbose",            0, &config_verbose,     1},
507
		{"alternate",          0, &config_alternate,   1},
508
		{"intensities",        1, NULL,               'q'},
Thomas White's avatar
Thomas White committed
509
		{"pdb",                1, NULL,               'p'},
510
		{"prefix",             1, NULL,               'x'},
Thomas White's avatar
Thomas White committed
511
		{"unpolarized",        0, &config_polar,       0},
512
		{"check-sanity",       0, &config_sanity,      1},
513
		{"sat-corr",           0, &config_satcorr,     1},
514
		{"no-sa",              0, &config_sa,          0},
Thomas White's avatar
Thomas White committed
515
516
517
518
		{0, 0, NULL, 0}
	};

	/* Short options */
Thomas White's avatar
Thomas White committed
519
520
	while ((c = getopt_long(argc, argv, "hi:wp:j:x:g:",
	                        longopts, NULL)) != -1) {
Thomas White's avatar
Thomas White committed
521
522

		switch (c) {
Thomas White's avatar
Thomas White committed
523
		case 'h' :
Thomas White's avatar
Thomas White committed
524
525
526
			show_help(argv[0]);
			return 0;

Thomas White's avatar
Thomas White committed
527
		case 'i' :
Thomas White's avatar
Thomas White committed
528
529
530
			filename = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
531
		case 'z' :
Thomas White's avatar
Thomas White committed
532
533
534
			indm_str = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
535
		case 'q' :
536
537
538
			intfile = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
539
		case 'p' :
Thomas White's avatar
Thomas White committed
540
541
542
			pdb = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
543
		case 'x' :
544
545
546
			prefix = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
547
		case 'j' :
548
549
550
			nthreads = atoi(optarg);
			break;

Thomas White's avatar
Thomas White committed
551
		case 'g' :
Thomas White's avatar
Thomas White committed
552
553
554
			geometry = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
555
		case 0 :
Thomas White's avatar
Thomas White committed
556
557
			break;

Thomas White's avatar
Thomas White committed
558
		default :
Thomas White's avatar
Thomas White committed
559
560
561
562
563
564
565
566
567
568
569
570
571
572
			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
573
		ERROR("Failed to open input file '%s'\n", filename);
Thomas White's avatar
Thomas White committed
574
575
		return 1;
	}
Thomas White's avatar
Thomas White committed
576
	free(filename);
Thomas White's avatar
Thomas White committed
577

578
	if ( intfile != NULL ) {
Thomas White's avatar
Thomas White committed
579
580
581
		ReflItemList *items;
		items = read_reflections(intfile, intensities, NULL, NULL);
		delete_items(items);
582
583
584
585
	} else {
		intensities = NULL;
	}

Thomas White's avatar
Thomas White committed
586
587
588
589
	if ( pdb == NULL ) {
		pdb = strdup("molecule.pdb");
	}

590
	if ( prefix == NULL ) {
Thomas White's avatar
Thomas White committed
591
		prefix = strdup("");
592
593
	}

594
	if ( (nthreads == 0) || (nthreads > MAX_THREADS) ) {
595
596
597
598
		ERROR("Invalid number of threads.\n");
		return 1;
	}

Thomas White's avatar
Thomas White committed
599
600
	if ( indm_str == NULL ) {
		STATUS("You didn't specify an indexing method, so I won't"
601
602
603
		       " 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
604
605
606
607
608
		indm = INDEXING_NONE;
	} else if ( strcmp(indm_str, "none") == 0 ) {
		indm = INDEXING_NONE;
	} else if ( strcmp(indm_str, "dirax") == 0) {
		indm = INDEXING_DIRAX;
609
610
	} else if ( strcmp(indm_str, "template") == 0) {
		indm = INDEXING_TEMPLATE;
Thomas White's avatar
Thomas White committed
611
612
613
614
	} else {
		ERROR("Unrecognised indexing method '%s'\n", indm_str);
		return 1;
	}
615
	free(indm_str);
Thomas White's avatar
Thomas White committed
616

Thomas White's avatar
Thomas White committed
617
618
619
620
621
622
623
624
625
626
627
628
	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);

629
	if ( (!config_nomatch) || (indm == INDEXING_TEMPLATE) ) {
630
631
		cell = load_cell_from_pdb(pdb);
		if ( cell == NULL ) {
632
633
634
			ERROR("Couldn't read unit cell (from %s)\n", pdb);
			return 1;
		}
635
636
637
	} else {
		STATUS("No cell needed because --no-match was used.\n");
		cell = NULL;
638
	}
639
	free(pdb);
640

641
642
643
644
645
646
647
648
	/* 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);
649
	ipriv = prepare_indexing(indm, cell, prepare_filename, det);
650
651
652
653
654
	if ( ipriv == NULL ) {
		ERROR("Failed to prepare indexing.\n");
		return 1;
	}

Thomas White's avatar
Thomas White committed
655
	gsl_set_error_handler_off();
Thomas White's avatar
Thomas White committed
656
657
	n_images = 0;
	n_hits = 0;
Thomas White's avatar
Thomas White committed
658
	n_sane = 0;
659

660
661
662
663
664
665
	for ( i=0; i<nthreads; i++ ) {
		worker_args[i] = malloc(sizeof(struct process_args));
		worker_args[i]->filename = malloc(1024);
		worker_active[i] = 0;
	}

666
	/* Start threads off */
667
	for ( i=0; i<nthreads; i++ ) {
Thomas White's avatar
Thomas White committed
668
669

		char line[1024];
670
671
		struct process_args *pargs;
		int r;
Thomas White's avatar
Thomas White committed
672

673
		pargs = worker_args[i];
Thomas White's avatar
Thomas White committed
674

675
676
677
678
679
680
681
		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
682
		chomp(line);
683
		snprintf(pargs->filename, 1023, "%s%s", prefix, line);
684

685
686
		n_images++;

687
		pargs->output_mutex = &output_mutex;
688
		pargs->gpu_mutex = &gpu_mutex;
689
		pthread_mutex_init(&pargs->control_mutex, NULL);
690
691
692
693
694
695
696
697
698
699
		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
700
		pargs->config_polar = config_polar;
701
		pargs->config_sanity = config_sanity;
702
		pargs->config_satcorr = config_satcorr;
703
		pargs->config_sa = config_sa;
704
		pargs->cell = cell;
Thomas White's avatar
Thomas White committed
705
		pargs->det = det;
706
		pargs->ipriv = ipriv;
707
708
709
		pargs->indm = indm;
		pargs->intensities = intensities;
		pargs->gctx = gctx;
Thomas White's avatar
Thomas White committed
710
		pargs->id = i;
711
712
713
714
715
		pthread_mutex_lock(&pargs->control_mutex);
		pargs->done = 0;
		pargs->start = 1;
		pargs->finish = 0;
		pthread_mutex_unlock(&pargs->control_mutex);
716

Thomas White's avatar
Thomas White committed
717
		worker_active[i] = 1;
718
		r = pthread_create(&workers[i], NULL, worker_thread, pargs);
719
		if ( r != 0 ) {
Thomas White's avatar
Thomas White committed
720
			worker_active[i] = 0;
721
722
723
724
725
			ERROR("Couldn't start thread %i\n", i);
		}

	}

726
	/* Keep threads busy until the end of the data */
727
728
729
730
731
732
733
734
	do {

		int i;

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

			char line[1024];
			struct process_args *pargs;
735
736
737
738
			int done;

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

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

743
			/* Has the thread finished yet? */
744
			pargs = worker_args[i];
745
746
747
748
			pthread_mutex_lock(&pargs->control_mutex);
			done = pargs->done;
			pthread_mutex_unlock(&pargs->control_mutex);
			if ( !done ) continue;
749

750
751
			/* Results will be processed after checking if
			 * there are any more images to process. */
752

753
			/* Get next filename */
754
			rval = fgets(line, 1023, fh);
755
756
757
			/* In this case, the result of the last file
			 * file will be processed when the thread is
			 * joined. */
758
			if ( rval == NULL ) break;
759
760
761
762
763

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

764
			chomp(line);
765
			snprintf(pargs->filename, 1023, "%s%s", prefix, line);
766
767

			n_images++;
768
769
770
771
772
773
774

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

775
		}
Thomas White's avatar
Thomas White committed
776

Thomas White's avatar
Thomas White committed
777
778
	} while ( rval != NULL );

779
	/* Join threads */
780
	for ( i=0; i<nthreads; i++ ) {
Thomas White's avatar
Thomas White committed
781

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

784
785
		/* Tell the thread to exit */
		struct process_args *pargs = worker_args[i];
Thomas White's avatar
Thomas White committed
786
		pthread_mutex_lock(&pargs->control_mutex);
787
		pargs->finish = 1;
Thomas White's avatar
Thomas White committed
788
		pthread_mutex_unlock(&pargs->control_mutex);
Thomas White's avatar
Thomas White committed
789

790
791
		/* Wait for it to join */
		pthread_join(workers[i], NULL);
Thomas White's avatar
Thomas White committed
792
793
		worker_active[i] = 0;

794
795
		n_hits += pargs->hit;
		n_sane += pargs->peaks_sane;
Thomas White's avatar
Thomas White committed
796

Thomas White's avatar
Thomas White committed
797
	free:
798
		if ( worker_args[i]->filename != NULL ) {
799
800
			free(worker_args[i]->filename);
		}
801
		free(worker_args[i]);
Thomas White's avatar
Thomas White committed
802

803
804
805
	}

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

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

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

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