indexamajig.c 20.5 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
 *
 * Part of CrystFEL - crystallography with a FEL
 *
 */


#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

17
#define _GNU_SOURCE 1
Thomas White's avatar
Thomas White committed
18
19
20
21
22
23
#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
24
#include <hdf5.h>
Thomas White's avatar
Thomas White committed
25
#include <gsl/gsl_errno.h>
26
#include <pthread.h>
Thomas White's avatar
Thomas White committed
27
#include <sys/time.h>
Thomas White's avatar
Thomas White committed
28
29
30

#include "utils.h"
#include "hdf5-file.h"
Thomas White's avatar
Thomas White committed
31
#include "index.h"
32
#include "peaks.h"
33
#include "diffraction.h"
Thomas White's avatar
Thomas White committed
34
#include "diffraction-gpu.h"
35
#include "detector.h"
36
#include "sfac.h"
Thomas White's avatar
Thomas White committed
37
#include "filters.h"
38
#include "reflections.h"
Thomas White's avatar
Thomas White committed
39
40


41
42
43
44
#define MAX_THREADS (96)

struct process_args
{
45
	/* Input */
46
	char *filename;
Thomas White's avatar
Thomas White committed
47
	int id;
48
49
	pthread_mutex_t *output_mutex;  /* Protects stdout */
	pthread_mutex_t *gpu_mutex;     /* Protects "gctx" */
50
51
52
53
54
55
56
57
58
59
60
	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
61
	int config_polar;
62
	int config_sanity;
63
	int config_satcorr;
64
	int config_sa;
Thomas White's avatar
Thomas White committed
65
	struct detector *det;
66
67
68
	IndexingMethod indm;
	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"
Thomas White's avatar
Thomas White committed
101
"  -g. --geometry=<file>   Get detector geometry from file.\n"
Thomas White's avatar
Thomas White committed
102
103
"\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"
104
"      --near-bragg        Output a list of reflection intensities to stdout.\n"
Thomas White's avatar
Thomas White committed
105
106
107
108
109
110
111
112
113
114
"                           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"
115
"      --simulate          Simulate the diffraction pattern using the indexed\n"
Thomas White's avatar
Thomas White committed
116
117
118
119
"                           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"
120
121
122
123
124
125
126
"      --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"
127
128
129
130
131
"      --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"
132
133
"                           The intensities in this list are from the\n"
"                           centroid/integration procedure.\n"
134
135
136
"      --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"
137
"      --unpolarized       Don't correct for the polarisation of the X-rays.\n"
138
139
"      --check-sanity      Check that indexed locations approximately correspond\n"
"                           with detected peaks.\n"
140
141
"      --sat-corr          Correct values of saturated peaks using a table\n"
"                           included in the HDF5 file.\n"
142
143
"      --no-sa             Don't correct for the differing solid angles of\n"
"                           the pixels.\n"
Thomas White's avatar
Thomas White committed
144
145
146
147
148
"\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"
"\n\nControl of model and data input:\n\n"
Thomas White's avatar
Thomas White committed
149
150
151
"     --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"
152
" -x, --prefix=<p>         Prefix filenames from input file with 'p'.\n"
Thomas White's avatar
Thomas White committed
153
);
Thomas White's avatar
Thomas White committed
154
155
156
}


157
static struct image *get_simage(struct image *template, int alternate)
158
{
159
	struct image *image;
160
	struct panel panels[2];
161

162
163
	image = malloc(sizeof(*image));

164
	/* Simulate a diffraction pattern */
165
166
167
	image->twotheta = NULL;
	image->data = NULL;
	image->det = template->det;
168
169
170
	image->flags = NULL;
	image->f0_available = 0;
	image->f0 = 1.0;
171
172

	/* View head-on (unit cell is tilted) */
173
174
175
176
	image->orientation.w = 1.0;
	image->orientation.x = 0.0;
	image->orientation.y = 0.0;
	image->orientation.z = 0.0;
177
178
179

	/* Detector geometry for the simulation
	 * - not necessarily the same as the original. */
180
181
	image->width = 1024;
	image->height = 1024;
Thomas White's avatar
Thomas White committed
182
	image->det->n_panels = 2;
183

184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
	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;
202
		panels[1].cy = 525.0;
203
204
205
		panels[1].clen = 56.7e-2;  /* 56.7 cm */
		panels[1].res = 13333.3;   /* 75 microns/pixel */

Thomas White's avatar
Thomas White committed
206
		image->det->panels = panels;
207

208
209
210
	} else {

		/* Copy pointer to old geometry */
Thomas White's avatar
Thomas White committed
211
		image->det->panels = template->det->panels;
212
213

	}
214

215
	image->lambda = ph_en_to_lambda(eV_to_J(1.8e3));
216
	image->features = template->features;
217
	image->filename = template->filename;
218
	image->indexed_cell = template->indexed_cell;
Thomas White's avatar
Thomas White committed
219
	image->f0 = template->f0;
220

Thomas White's avatar
Thomas White committed
221
222
223
224
	/* Prevent muppetry */
	image->hits = NULL;
	image->n_hits = 0;

225
226
	return image;
}
227
228


229
static void simulate_and_write(struct image *simage, struct gpu_context **gctx,
Thomas White's avatar
Thomas White committed
230
                               const double *intensities, UnitCell *cell)
231
{
232
	/* Set up GPU if necessary */
Thomas White's avatar
Thomas White committed
233
	if ( (gctx != NULL) && (*gctx == NULL) ) {
Thomas White's avatar
Thomas White committed
234
		*gctx = setup_gpu(0, simage, intensities);
235
236
	}

Thomas White's avatar
Thomas White committed
237
	if ( (gctx != NULL) && (*gctx != NULL) ) {
238
		get_diffraction_gpu(*gctx, simage, 24, 24, 40, cell);
239
	} else {
240
		get_diffraction(simage, 24, 24, 40,
Thomas White's avatar
Thomas White committed
241
		                intensities, NULL, cell, 0,
242
		                GRADIENT_MOSAIC);
243
	}
244
	record_image(simage, 0);
245

246
	hdf5_write("simulated.h5", simage->data, simage->width, simage->height,
247
248
249
250
		   H5T_NATIVE_FLOAT);
}


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

	image.features = NULL;
	image.data = NULL;
	image.indexed_cell = NULL;
Thomas White's avatar
Thomas White committed
279
	image.id = pargs->id;
280
	image.filename = filename;
Thomas White's avatar
Thomas White committed
281
282
	image.hits = NULL;
	image.n_hits = 0;
Thomas White's avatar
Thomas White committed
283
	image.det = pargs->det;
284

285
286
287
288
289
290
	/* 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;

291
	STATUS("Processing '%s'\n", image.filename);
292

293
294
	result.peaks_sane = 0;
	result.hit = 0;
295

296
297
	hdfile = hdfile_open(filename);
	if ( hdfile == NULL ) {
298
		return result;
299
300
	} else if ( hdfile_set_first_image(hdfile, "/") ) {
		ERROR("Couldn't select path\n");
301
		return result;
302
303
	}

304
	hdf5_read(hdfile, &image, pargs->config_satcorr);
305
306
307
308
309
310
311
312
313
314
315
316
317

	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 {
318
		memcpy(data_for_measurement, image.data, data_size);
319
320
321
322
	}

	/* Perform 'fine' peak search */
	search_peaks(&image);
Thomas White's avatar
Thomas White committed
323
324
325
326
327
328

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

329
330
	if ( image_feature_count(image.features) < 5 ) goto done;

331
	if ( config_dumpfound ) dump_peaks(&image, pargs->output_mutex);
332
333
334
335
336
337
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,
		              config_verbose);
	}

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

Thomas White's avatar
Thomas White committed
348
	/* Sanity check */
349
350
	if ( pargs->config_sanity
	  && !peak_sanity_check(&image, image.indexed_cell) ) {
Thomas White's avatar
Thomas White committed
351
352
353
		STATUS("Failed peak sanity check.\n");
		goto done;
	} else {
354
		result.peaks_sane = 1;
Thomas White's avatar
Thomas White committed
355
356
	}

357
358
	/* Measure intensities if requested */
	if ( config_nearbragg ) {
359
		output_intensities(&image, image.indexed_cell,
Thomas White's avatar
Thomas White committed
360
		                   pargs->output_mutex, config_polar,
361
		                   pargs->config_sa);
362
363
	}

364
365
	simage = get_simage(&image, config_alternate);

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

	/* 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
389
	free(image.flags);
390
	image_feature_list_free(image.features);
Thomas White's avatar
Thomas White committed
391
	free(image.hits);
392
393
	hdfile_close(hdfile);

394
	if ( image.indexed_cell == NULL ) {
395
		result.hit = 0;
396
	} else {
397
		result.hit = 1;
398
399
	}
	return result;
400
401
402
}


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

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

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

		switch (c) {
Thomas White's avatar
Thomas White committed
514
		case 'h' :
Thomas White's avatar
Thomas White committed
515
516
517
			show_help(argv[0]);
			return 0;

Thomas White's avatar
Thomas White committed
518
		case 'i' :
Thomas White's avatar
Thomas White committed
519
520
521
			filename = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
522
		case 'z' :
Thomas White's avatar
Thomas White committed
523
524
525
			indm_str = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
526
		case 'q' :
527
528
529
			intfile = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
530
		case 'p' :
Thomas White's avatar
Thomas White committed
531
532
533
			pdb = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
534
		case 'x' :
535
536
537
			prefix = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
538
		case 'j' :
539
540
541
			nthreads = atoi(optarg);
			break;

Thomas White's avatar
Thomas White committed
542
		case 'g' :
Thomas White's avatar
Thomas White committed
543
544
545
			geometry = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
546
		case 0 :
Thomas White's avatar
Thomas White committed
547
548
			break;

Thomas White's avatar
Thomas White committed
549
		default :
Thomas White's avatar
Thomas White committed
550
551
552
553
554
555
556
557
558
559
560
561
562
563
			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
564
		ERROR("Failed to open input file '%s'\n", filename);
Thomas White's avatar
Thomas White committed
565
566
		return 1;
	}
Thomas White's avatar
Thomas White committed
567
	free(filename);
Thomas White's avatar
Thomas White committed
568

569
	if ( intfile != NULL ) {
Thomas White's avatar
Thomas White committed
570
571
572
		ReflItemList *items;
		items = read_reflections(intfile, intensities, NULL, NULL);
		delete_items(items);
573
574
575
576
	} else {
		intensities = NULL;
	}

Thomas White's avatar
Thomas White committed
577
578
579
580
	if ( pdb == NULL ) {
		pdb = strdup("molecule.pdb");
	}

581
	if ( prefix == NULL ) {
Thomas White's avatar
Thomas White committed
582
		prefix = strdup("");
583
584
	}

585
	if ( (nthreads == 0) || (nthreads > MAX_THREADS) ) {
586
587
588
589
		ERROR("Invalid number of threads.\n");
		return 1;
	}

Thomas White's avatar
Thomas White committed
590
591
	if ( indm_str == NULL ) {
		STATUS("You didn't specify an indexing method, so I won't"
592
593
594
		       " 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
595
596
597
598
599
600
601
602
603
		indm = INDEXING_NONE;
	} else if ( strcmp(indm_str, "none") == 0 ) {
		indm = INDEXING_NONE;
	} else if ( strcmp(indm_str, "dirax") == 0) {
		indm = INDEXING_DIRAX;
	} else {
		ERROR("Unrecognised indexing method '%s'\n", indm_str);
		return 1;
	}
604
	free(indm_str);
Thomas White's avatar
Thomas White committed
605

Thomas White's avatar
Thomas White committed
606
607
608
609
610
611
612
613
614
615
616
617
	if ( geometry == NULL ) {
		ERROR("You need to specify a geometry file with --geometry\n");
		return 1;
	}

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

Thomas White's avatar
Thomas White committed
618
	cell = load_cell_from_pdb(pdb);
619
	if ( cell == NULL ) {
620
621
622
623
624
625
626
		if ( config_nomatch ) {
			STATUS("Couldn't read unit cell, but no matching is"
			       " requested, so that's OK.\n");
		} else {
			ERROR("Couldn't read unit cell (from %s)\n", pdb);
			return 1;
		}
627
	}
628
	free(pdb);
629

Thomas White's avatar
Thomas White committed
630
	gsl_set_error_handler_off();
Thomas White's avatar
Thomas White committed
631
632
	n_images = 0;
	n_hits = 0;
Thomas White's avatar
Thomas White committed
633
	n_sane = 0;
634

635
636
637
638
639
640
	for ( i=0; i<nthreads; i++ ) {
		worker_args[i] = malloc(sizeof(struct process_args));
		worker_args[i]->filename = malloc(1024);
		worker_active[i] = 0;
	}

641
	/* Start threads off */
642
	for ( i=0; i<nthreads; i++ ) {
Thomas White's avatar
Thomas White committed
643
644

		char line[1024];
645
646
		struct process_args *pargs;
		int r;
Thomas White's avatar
Thomas White committed
647

648
		pargs = worker_args[i];
Thomas White's avatar
Thomas White committed
649

Thomas White's avatar
Thomas White committed
650
		rval = fgets(line, 1023, fh);
Thomas White's avatar
Thomas White committed
651
		if ( rval == NULL ) continue;
Thomas White's avatar
Thomas White committed
652
		chomp(line);
653
		snprintf(pargs->filename, 1023, "%s%s", prefix, line);
654

655
656
		n_images++;

657
		pargs->output_mutex = &output_mutex;
658
		pargs->gpu_mutex = &gpu_mutex;
659
		pthread_mutex_init(&pargs->control_mutex, NULL);
660
661
662
663
664
665
666
667
668
669
		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
670
		pargs->config_polar = config_polar;
671
		pargs->config_sanity = config_sanity;
672
		pargs->config_satcorr = config_satcorr;
673
		pargs->config_sa = config_sa;
674
		pargs->cell = cell;
Thomas White's avatar
Thomas White committed
675
		pargs->det = det;
676
677
678
		pargs->indm = indm;
		pargs->intensities = intensities;
		pargs->gctx = gctx;
Thomas White's avatar
Thomas White committed
679
		pargs->id = i;
680
681
682
683
684
		pthread_mutex_lock(&pargs->control_mutex);
		pargs->done = 0;
		pargs->start = 1;
		pargs->finish = 0;
		pthread_mutex_unlock(&pargs->control_mutex);
685

Thomas White's avatar
Thomas White committed
686
		worker_active[i] = 1;
687
		r = pthread_create(&workers[i], NULL, worker_thread, pargs);
688
		if ( r != 0 ) {
Thomas White's avatar
Thomas White committed
689
			worker_active[i] = 0;
690
691
692
693
694
			ERROR("Couldn't start thread %i\n", i);
		}

	}

695
	/* Keep threads busy until the end of the data */
696
697
698
699
700
701
702
703
	do {

		int i;

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

			char line[1024];
			struct process_args *pargs;
704
705
706
707
			int done;

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

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

712
			/* Has the thread finished yet? */
713
			pargs = worker_args[i];
714
715
716
717
			pthread_mutex_lock(&pargs->control_mutex);
			done = pargs->done;
			pthread_mutex_unlock(&pargs->control_mutex);
			if ( !done ) continue;
718

719
720
721
			/* Record the result */
			n_hits += pargs->hit;
			n_sane += pargs->peaks_sane;
722

723
			/* Get next filename */
724
725
726
			rval = fgets(line, 1023, fh);
			if ( rval == NULL ) break;
			chomp(line);
727
			snprintf(pargs->filename, 1023, "%s%s", prefix, line);
728
729

			n_images++;
730
731
732
733
734
735
736

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

737
		}
Thomas White's avatar
Thomas White committed
738

Thomas White's avatar
Thomas White committed
739
740
	} while ( rval != NULL );

741
	/* Join threads */
742
	for ( i=0; i<nthreads; i++ ) {
Thomas White's avatar
Thomas White committed
743

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

746
747
		/* Tell the thread to exit */
		struct process_args *pargs = worker_args[i];
Thomas White's avatar
Thomas White committed
748
		pthread_mutex_lock(&pargs->control_mutex);
749
		pargs->finish = 1;
Thomas White's avatar
Thomas White committed
750
		pthread_mutex_unlock(&pargs->control_mutex);
Thomas White's avatar
Thomas White committed
751

752
753
		/* Wait for it to join */
		pthread_join(workers[i], NULL);
Thomas White's avatar
Thomas White committed
754
755
		worker_active[i] = 0;

756
757
		n_hits += pargs->hit;
		n_sane += pargs->peaks_sane;
Thomas White's avatar
Thomas White committed
758

Thomas White's avatar
Thomas White committed
759
	free:
760
		if ( worker_args[i]->filename != NULL ) {
761
762
			free(worker_args[i]->filename);
		}
763
		free(worker_args[i]);
Thomas White's avatar
Thomas White committed
764

765
766
767
	}

	free(prefix);
Thomas White's avatar
Thomas White committed
768
769
	free(det->panels);
	free(det);
770
	free(cell);
Thomas White's avatar
Thomas White committed
771
772
773
	fclose(fh);

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

Thomas White's avatar
Thomas White committed
776
777
778
779
	if ( gctx != NULL ) {
		cleanup_gpu(gctx);
	}

Thomas White's avatar
Thomas White committed
780
781
	return 0;
}