indexamajig.c 20.4 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 ) {
Thomas White's avatar
Thomas White committed
620
		ERROR("Couldn't read unit cell (from %s)\n", pdb);
621
622
		return 1;
	}
623
	free(pdb);
624

Thomas White's avatar
Thomas White committed
625
	gsl_set_error_handler_off();
Thomas White's avatar
Thomas White committed
626
627
	n_images = 0;
	n_hits = 0;
Thomas White's avatar
Thomas White committed
628
	n_sane = 0;
629

630
631
632
633
634
635
	for ( i=0; i<nthreads; i++ ) {
		worker_args[i] = malloc(sizeof(struct process_args));
		worker_args[i]->filename = malloc(1024);
		worker_active[i] = 0;
	}

636
	/* Start threads off */
637
	for ( i=0; i<nthreads; i++ ) {
Thomas White's avatar
Thomas White committed
638
639

		char line[1024];
640
641
		struct process_args *pargs;
		int r;
Thomas White's avatar
Thomas White committed
642

643
		pargs = worker_args[i];
Thomas White's avatar
Thomas White committed
644

Thomas White's avatar
Thomas White committed
645
		rval = fgets(line, 1023, fh);
Thomas White's avatar
Thomas White committed
646
		if ( rval == NULL ) continue;
Thomas White's avatar
Thomas White committed
647
		chomp(line);
648
		snprintf(pargs->filename, 1023, "%s%s", prefix, line);
649

650
651
		n_images++;

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

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

	}

690
	/* Keep threads busy until the end of the data */
691
692
693
694
695
696
697
698
	do {

		int i;

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

			char line[1024];
			struct process_args *pargs;
699
700
701
702
			int done;

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

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

707
			/* Has the thread finished yet? */
708
			pargs = worker_args[i];
709
710
711
712
			pthread_mutex_lock(&pargs->control_mutex);
			done = pargs->done;
			pthread_mutex_unlock(&pargs->control_mutex);
			if ( !done ) continue;
713

714
715
716
			/* Record the result */
			n_hits += pargs->hit;
			n_sane += pargs->peaks_sane;
717

718
			/* Get next filename */
719
720
721
			rval = fgets(line, 1023, fh);
			if ( rval == NULL ) break;
			chomp(line);
722
			snprintf(pargs->filename, 1023, "%s%s", prefix, line);
723
724

			n_images++;
725
726
727
728
729
730
731

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

732
		}
Thomas White's avatar
Thomas White committed
733

Thomas White's avatar
Thomas White committed
734
735
	} while ( rval != NULL );

736
	/* Join threads */
737
	for ( i=0; i<nthreads; i++ ) {
Thomas White's avatar
Thomas White committed
738

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

741
742
		/* Tell the thread to exit */
		struct process_args *pargs = worker_args[i];
Thomas White's avatar
Thomas White committed
743
		pthread_mutex_lock(&pargs->control_mutex);
744
		pargs->finish = 1;
Thomas White's avatar
Thomas White committed
745
		pthread_mutex_unlock(&pargs->control_mutex);
Thomas White's avatar
Thomas White committed
746

747
748
		/* Wait for it to join */
		pthread_join(workers[i], NULL);
Thomas White's avatar
Thomas White committed
749
750
		worker_active[i] = 0;

751
752
		n_hits += pargs->hit;
		n_sane += pargs->peaks_sane;
Thomas White's avatar
Thomas White committed
753

Thomas White's avatar
Thomas White committed
754
	free:
755
		if ( worker_args[i]->filename != NULL ) {
756
757
			free(worker_args[i]->filename);
		}
758
		free(worker_args[i]);
Thomas White's avatar
Thomas White committed
759

760
761
762
	}

	free(prefix);
Thomas White's avatar
Thomas White committed
763
764
	free(det->panels);
	free(det);
765
	free(cell);
Thomas White's avatar
Thomas White committed
766
767
768
	fclose(fh);

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

Thomas White's avatar
Thomas White committed
771
772
773
774
	if ( gctx != NULL ) {
		cleanup_gpu(gctx);
	}

Thomas White's avatar
Thomas White committed
775
776
	return 0;
}