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

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

80

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


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


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

170
171
	image = malloc(sizeof(*image));

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

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

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

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

Thomas White's avatar
Thomas White committed
214
		image->det->panels = panels;
215

216
217
218
	} else {

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

	}
222

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

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

233
234
	return image;
}
235
236


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

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

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


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

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

293
294
295
296
297
298
	/* 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;

299
	STATUS("Processing '%s'\n", image.filename);
300

301
302
	result.peaks_sane = 0;
	result.hit = 0;
303

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

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

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

	/* Perform 'fine' peak search */
330
	search_peaks(&image, pargs->threshold);
Thomas White's avatar
Thomas White committed
331
332
333
334
335
336

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

337
338
	if ( image_feature_count(image.features) < 5 ) goto done;

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

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

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

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

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

372
373
	simage = get_simage(&image, config_alternate);

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

	/* 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
397
	free(image.flags);
398
	image_feature_list_free(image.features);
Thomas White's avatar
Thomas White committed
399
	free(image.hits);
400
401
	hdfile_close(hdfile);

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

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

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

		switch (c) {
Thomas White's avatar
Thomas White committed
527
		case 'h' :
Thomas White's avatar
Thomas White committed
528
529
530
			show_help(argv[0]);
			return 0;

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

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

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

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

Thomas White's avatar
Thomas White committed
547
		case 'x' :
548
549
550
			prefix = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
551
		case 'j' :
552
553
554
			nthreads = atoi(optarg);
			break;

Thomas White's avatar
Thomas White committed
555
		case 'g' :
Thomas White's avatar
Thomas White committed
556
557
558
			geometry = strdup(optarg);
			break;

559
560
561
562
		case 't' :
			threshold = strtof(optarg, NULL);
			break;

Thomas White's avatar
Thomas White committed
563
		case 0 :
Thomas White's avatar
Thomas White committed
564
565
			break;

Thomas White's avatar
Thomas White committed
566
		default :
Thomas White's avatar
Thomas White committed
567
568
569
570
571
572
573
574
575
576
577
578
579
580
			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
581
		ERROR("Failed to open input file '%s'\n", filename);
Thomas White's avatar
Thomas White committed
582
583
		return 1;
	}
Thomas White's avatar
Thomas White committed
584
	free(filename);
Thomas White's avatar
Thomas White committed
585

586
	if ( intfile != NULL ) {
Thomas White's avatar
Thomas White committed
587
588
589
		ReflItemList *items;
		items = read_reflections(intfile, intensities, NULL, NULL);
		delete_items(items);
590
591
592
593
	} else {
		intensities = NULL;
	}

Thomas White's avatar
Thomas White committed
594
595
596
597
	if ( pdb == NULL ) {
		pdb = strdup("molecule.pdb");
	}

598
	if ( prefix == NULL ) {
Thomas White's avatar
Thomas White committed
599
		prefix = strdup("");
600
601
	}

602
	if ( (nthreads == 0) || (nthreads > MAX_THREADS) ) {
603
604
605
606
		ERROR("Invalid number of threads.\n");
		return 1;
	}

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

Thomas White's avatar
Thomas White committed
625
626
627
628
629
630
631
632
633
634
635
636
	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);

637
	if ( (!config_nomatch) || (indm == INDEXING_TEMPLATE) ) {
638
639
		cell = load_cell_from_pdb(pdb);
		if ( cell == NULL ) {
640
641
642
			ERROR("Couldn't read unit cell (from %s)\n", pdb);
			return 1;
		}
643
644
645
	} else {
		STATUS("No cell needed because --no-match was used.\n");
		cell = NULL;
646
	}
647
	free(pdb);
648

649
650
651
652
653
654
655
656
	/* 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);
657
	ipriv = prepare_indexing(indm, cell, prepare_filename, det);
658
659
660
661
662
	if ( ipriv == NULL ) {
		ERROR("Failed to prepare indexing.\n");
		return 1;
	}

Thomas White's avatar
Thomas White committed
663
	gsl_set_error_handler_off();
Thomas White's avatar
Thomas White committed
664
665
	n_images = 0;
	n_hits = 0;
Thomas White's avatar
Thomas White committed
666
	n_sane = 0;
667

668
669
670
671
672
673
	for ( i=0; i<nthreads; i++ ) {
		worker_args[i] = malloc(sizeof(struct process_args));
		worker_args[i]->filename = malloc(1024);
		worker_active[i] = 0;
	}

674
	/* Start threads off */
675
	for ( i=0; i<nthreads; i++ ) {
Thomas White's avatar
Thomas White committed
676
677

		char line[1024];
678
679
		struct process_args *pargs;
		int r;
Thomas White's avatar
Thomas White committed
680

681
		pargs = worker_args[i];
Thomas White's avatar
Thomas White committed
682

683
684
685
686
687
688
689
		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
690
		chomp(line);
691
		snprintf(pargs->filename, 1023, "%s%s", prefix, line);
692

693
694
		n_images++;

695
		pargs->output_mutex = &output_mutex;
696
		pargs->gpu_mutex = &gpu_mutex;
697
		pthread_mutex_init(&pargs->control_mutex, NULL);
698
699
700
701
702
703
704
705
706
707
		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
708
		pargs->config_polar = config_polar;
709
		pargs->config_sanity = config_sanity;
710
		pargs->config_satcorr = config_satcorr;
711
		pargs->config_sa = config_sa;
712
		pargs->cell = cell;
Thomas White's avatar
Thomas White committed
713
		pargs->det = det;
714
		pargs->ipriv = ipriv;
715
716
717
		pargs->indm = indm;
		pargs->intensities = intensities;
		pargs->gctx = gctx;
718
		pargs->threshold = threshold;
Thomas White's avatar
Thomas White committed
719
		pargs->id = i;
720
721
722
723
724
		pthread_mutex_lock(&pargs->control_mutex);
		pargs->done = 0;
		pargs->start = 1;
		pargs->finish = 0;
		pthread_mutex_unlock(&pargs->control_mutex);
725

Thomas White's avatar
Thomas White committed
726
		worker_active[i] = 1;
727
		r = pthread_create(&workers[i], NULL, worker_thread, pargs);
728
		if ( r != 0 ) {
Thomas White's avatar
Thomas White committed
729
			worker_active[i] = 0;
730
731
732
733
734
			ERROR("Couldn't start thread %i\n", i);
		}

	}

735
	/* Keep threads busy until the end of the data */
736
737
738
739
740
741
742
743
	do {

		int i;

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

			char line[1024];
			struct process_args *pargs;
744
745
746
747
			int done;

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

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

752
			/* Has the thread finished yet? */
753
			pargs = worker_args[i];
754
755
756
757
			pthread_mutex_lock(&pargs->control_mutex);
			done = pargs->done;
			pthread_mutex_unlock(&pargs->control_mutex);
			if ( !done ) continue;
758

759
760
			/* Results will be processed after checking if
			 * there are any more images to process. */
761

762
			/* Get next filename */
763
			rval = fgets(line, 1023, fh);
764
765
766
			/* In this case, the result of the last file
			 * file will be processed when the thread is
			 * joined. */
767
			if ( rval == NULL ) break;
768
769
770
771
772

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

773
			chomp(line);
774
			snprintf(pargs->filename, 1023, "%s%s", prefix, line);
775
776

			n_images++;
777
778
779
780
781
782
783

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

784
		}
Thomas White's avatar
Thomas White committed
785

Thomas White's avatar
Thomas White committed
786
787
	} while ( rval != NULL );

788
	/* Join threads */
789
	for ( i=0; i<nthreads; i++ ) {
Thomas White's avatar
Thomas White committed
790

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

793
794
		/* Tell the thread to exit */
		struct process_args *pargs = worker_args[i];
Thomas White's avatar
Thomas White committed
795
		pthread_mutex_lock(&pargs->control_mutex);
796
		pargs->finish = 1;
Thomas White's avatar
Thomas White committed
797
		pthread_mutex_unlock(&pargs->control_mutex);
Thomas White's avatar
Thomas White committed
798

799
800
		/* Wait for it to join */
		pthread_join(workers[i], NULL);
Thomas White's avatar
Thomas White committed
801
802
		worker_active[i] = 0;

803
804
		n_hits += pargs->hit;
		n_sane += pargs->peaks_sane;
Thomas White's avatar
Thomas White committed
805

Thomas White's avatar
Thomas White committed
806
	free:
807
		if ( worker_args[i]->filename != NULL ) {
808
809
			free(worker_args[i]->filename);
		}
810
		free(worker_args[i]);
Thomas White's avatar
Thomas White committed
811

812
813
814
	}

	free(prefix);
Thomas White's avatar
Thomas White committed
815
816
	free(det->panels);
	free(det);
817
	free(cell);
Thomas White's avatar
Thomas White committed
818
819
820
	fclose(fh);

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

Thomas White's avatar
Thomas White committed
823
824
825
826
	if ( gctx != NULL ) {
		cleanup_gpu(gctx);
	}

Thomas White's avatar
Thomas White committed
827
828
	return 0;
}