indexamajig.c 22.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>
27
28
#include <sys/types.h>
#include <sys/stat.h>
Thomas White's avatar
Thomas White committed
29
30
31

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


42
43
44
45
#define MAX_THREADS (96)

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

	/* 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;
80
81
};

82

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


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


170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
static char *check_prefix(char *prefix)
{
	int r;
	struct stat statbuf;
	char *new;
	size_t len;

	/* Is "prefix" a directory? */
	r = stat(prefix, &statbuf);
	if ( r != 0 ) {
		/* "prefix" probably doesn't exist.  This is fine - assume
		 * the user knows what they're doing, and that "prefix"
		 * suffixed with the actual filename will produce something
		 * sensible. */
		return prefix;
	}

	if ( !S_ISDIR(statbuf.st_mode) ) {
		/* Also fine, as above. */
		return prefix;
	}

	/* Does the prefix end in a slash? */
	if ( prefix[strlen(prefix)-1] == '/' ) {
		/* This looks sensible. */
		return prefix;
	}

	STATUS("Your prefix ('%s') is a directory, but doesn't end"
	       " with a slash.  I'm going to add it for you.\n", prefix);
	len = strlen(prefix)+2;
	new = malloc(len);
	snprintf(new, len, "%s/", prefix);
	free(prefix);
	return new;
}


208
static struct image *get_simage(struct image *template, int alternate)
209
{
210
	struct image *image;
211
	struct panel panels[2];
212

213
214
	image = malloc(sizeof(*image));

215
	/* Simulate a diffraction pattern */
216
217
218
	image->twotheta = NULL;
	image->data = NULL;
	image->det = template->det;
219
220
221
	image->flags = NULL;
	image->f0_available = 0;
	image->f0 = 1.0;
222
223

	/* View head-on (unit cell is tilted) */
224
225
226
227
	image->orientation.w = 1.0;
	image->orientation.x = 0.0;
	image->orientation.y = 0.0;
	image->orientation.z = 0.0;
228
229
230

	/* Detector geometry for the simulation
	 * - not necessarily the same as the original. */
231
232
	image->width = 1024;
	image->height = 1024;
Thomas White's avatar
Thomas White committed
233
	image->det->n_panels = 2;
234

235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
	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;
253
		panels[1].cy = 525.0;
254
255
256
		panels[1].clen = 56.7e-2;  /* 56.7 cm */
		panels[1].res = 13333.3;   /* 75 microns/pixel */

Thomas White's avatar
Thomas White committed
257
		image->det->panels = panels;
258

259
260
261
	} else {

		/* Copy pointer to old geometry */
Thomas White's avatar
Thomas White committed
262
		image->det->panels = template->det->panels;
263
264

	}
265

266
	image->lambda = ph_en_to_lambda(eV_to_J(1.8e3));
267
	image->features = template->features;
268
	image->filename = template->filename;
269
	image->indexed_cell = template->indexed_cell;
Thomas White's avatar
Thomas White committed
270
	image->f0 = template->f0;
271

Thomas White's avatar
Thomas White committed
272
273
274
275
	/* Prevent muppetry */
	image->hits = NULL;
	image->n_hits = 0;

276
277
	return image;
}
278
279


280
static void simulate_and_write(struct image *simage, struct gpu_context **gctx,
Thomas White's avatar
Thomas White committed
281
                               const double *intensities, UnitCell *cell)
282
{
283
	/* Set up GPU if necessary */
Thomas White's avatar
Thomas White committed
284
	if ( (gctx != NULL) && (*gctx == NULL) ) {
Thomas White's avatar
Thomas White committed
285
		*gctx = setup_gpu(0, simage, intensities);
286
287
	}

Thomas White's avatar
Thomas White committed
288
	if ( (gctx != NULL) && (*gctx != NULL) ) {
289
		get_diffraction_gpu(*gctx, simage, 24, 24, 40, cell);
290
	} else {
291
		get_diffraction(simage, 24, 24, 40,
Thomas White's avatar
Thomas White committed
292
		                intensities, NULL, cell, 0,
293
		                GRADIENT_MOSAIC);
294
	}
295
	record_image(simage, 0);
296

297
	hdf5_write("simulated.h5", simage->data, simage->width, simage->height,
298
299
300
301
		   H5T_NATIVE_FLOAT);
}


302
static struct process_result process_image(struct process_args *pargs)
303
304
305
306
307
{
	struct hdfile *hdfile;
	struct image image;
	struct image *simage;
	float *data_for_measurement;
308
	struct process_result result;
309
	size_t data_size;
310
311
312
313
314
315
316
317
318
319
320
321
	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
322
	int config_polar = pargs->config_polar;
323
324
325
	IndexingMethod indm = pargs->indm;
	const double *intensities = pargs->intensities;
	struct gpu_context *gctx = pargs->gctx;
326
327
328
329

	image.features = NULL;
	image.data = NULL;
	image.indexed_cell = NULL;
Thomas White's avatar
Thomas White committed
330
	image.id = pargs->id;
331
	image.filename = filename;
Thomas White's avatar
Thomas White committed
332
333
	image.hits = NULL;
	image.n_hits = 0;
Thomas White's avatar
Thomas White committed
334
	image.det = pargs->det;
335

336
337
338
339
340
341
	/* 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;

342
	STATUS("Processing '%s'\n", image.filename);
343

344
345
	result.peaks_sane = 0;
	result.hit = 0;
346

347
348
	hdfile = hdfile_open(filename);
	if ( hdfile == NULL ) {
349
		return result;
350
351
	} else if ( hdfile_set_first_image(hdfile, "/") ) {
		ERROR("Couldn't select path\n");
352
		return result;
353
354
	}

355
	hdf5_read(hdfile, &image, pargs->config_satcorr);
356
357
358
359
360
361
362
363
364
365
366
367
368

	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 {
369
		memcpy(data_for_measurement, image.data, data_size);
370
371
372
	}

	/* Perform 'fine' peak search */
373
	search_peaks(&image, pargs->threshold);
Thomas White's avatar
Thomas White committed
374
375
376
377
378
379

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

380
381
	if ( image_feature_count(image.features) < 5 ) goto done;

382
	if ( config_dumpfound ) dump_peaks(&image, pargs->output_mutex);
383
384
385
386
387
388
389
390
391
392

	/* 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,
393
		              config_verbose, pargs->ipriv);
394
395
396
397
398
	}

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

Thomas White's avatar
Thomas White committed
399
	/* Sanity check */
400
401
	if ( pargs->config_sanity
	  && !peak_sanity_check(&image, image.indexed_cell) ) {
Thomas White's avatar
Thomas White committed
402
403
404
		STATUS("Failed peak sanity check.\n");
		goto done;
	} else {
405
		result.peaks_sane = 1;
Thomas White's avatar
Thomas White committed
406
407
	}

408
409
	/* Measure intensities if requested */
	if ( config_nearbragg ) {
410
		output_intensities(&image, image.indexed_cell,
Thomas White's avatar
Thomas White committed
411
		                   pargs->output_mutex, config_polar,
412
		                   pargs->config_sa);
413
414
	}

415
416
	simage = get_simage(&image, config_alternate);

417
418
419
	/* Simulate if requested */
	if ( config_simulate ) {
		if ( config_gpu ) {
420
			pthread_mutex_lock(pargs->gpu_mutex);
421
			simulate_and_write(simage, &gctx, intensities,
Thomas White's avatar
Thomas White committed
422
			                   image.indexed_cell);
423
			pthread_mutex_unlock(pargs->gpu_mutex);
424
425
		} else {
			simulate_and_write(simage, NULL, intensities,
Thomas White's avatar
Thomas White committed
426
			                   image.indexed_cell);
427
428
429
430
431
432
433
434
435
436
437
438
439
		}
	}

	/* 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
440
	free(image.flags);
441
	image_feature_list_free(image.features);
Thomas White's avatar
Thomas White committed
442
	free(image.hits);
443
444
	hdfile_close(hdfile);

445
	if ( image.indexed_cell == NULL ) {
446
		result.hit = 0;
447
	} else {
448
		result.hit = 1;
449
450
	}
	return result;
451
452
453
}


454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
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
491
492
493
int main(int argc, char *argv[])
{
	int c;
Thomas White's avatar
Thomas White committed
494
	struct gpu_context *gctx = NULL;
Thomas White's avatar
Thomas White committed
495
496
	char *filename = NULL;
	FILE *fh;
497
	char *rval = NULL;
Thomas White's avatar
Thomas White committed
498
499
	int n_images;
	int n_hits;
Thomas White's avatar
Thomas White committed
500
	int n_sane;
501
	int config_noindex = 0;
Thomas White's avatar
Thomas White committed
502
	int config_dumpfound = 0;
503
	int config_nearbragg = 0;
Thomas White's avatar
Thomas White committed
504
	int config_writedrx = 0;
505
	int config_simulate = 0;
506
507
	int config_cmfilter = 0;
	int config_noisefilter = 0;
508
	int config_nomatch = 0;
Thomas White's avatar
Thomas White committed
509
	int config_gpu = 0;
510
	int config_verbose = 0;
511
	int config_alternate = 0;
Thomas White's avatar
Thomas White committed
512
	int config_polar = 1;
513
	int config_sanity = 0;
514
	int config_satcorr = 0;
515
	int config_sa = 1;
516
	float threshold = 800.0;
Thomas White's avatar
Thomas White committed
517
518
	struct detector *det;
	char *geometry = NULL;
Thomas White's avatar
Thomas White committed
519
520
	IndexingMethod indm;
	char *indm_str = NULL;
521
522
523
	UnitCell *cell;
	double *intensities = NULL;
	char *intfile = NULL;
Thomas White's avatar
Thomas White committed
524
	char *pdb = NULL;
525
	char *prefix = NULL;
526
527
528
	int nthreads = 1;
	pthread_t workers[MAX_THREADS];
	struct process_args *worker_args[MAX_THREADS];
Thomas White's avatar
Thomas White committed
529
	int worker_active[MAX_THREADS];
530
	int i;
531
	pthread_mutex_t output_mutex = PTHREAD_MUTEX_INITIALIZER;
532
	pthread_mutex_t gpu_mutex = PTHREAD_MUTEX_INITIALIZER;
533
534
535
	char prepare_line[1024];
	char prepare_filename[1024];
	IndexingPrivate *ipriv;
Thomas White's avatar
Thomas White committed
536
537
538
539
540

	/* Long options */
	const struct option longopts[] = {
		{"help",               0, NULL,               'h'},
		{"input",              1, NULL,               'i'},
Thomas White's avatar
Thomas White committed
541
		{"gpu",                0, &config_gpu,         1},
542
		{"no-index",           0, &config_noindex,     1},
Thomas White's avatar
Thomas White committed
543
		{"dump-peaks",         0, &config_dumpfound,   1},
544
		{"near-bragg",         0, &config_nearbragg,   1},
Thomas White's avatar
Thomas White committed
545
546
		{"write-drx",          0, &config_writedrx,    1},
		{"indexing",           1, NULL,               'z'},
Thomas White's avatar
Thomas White committed
547
		{"geometry",           1, NULL,               'g'},
548
		{"simulate",           0, &config_simulate,    1},
549
550
		{"filter-cm",          0, &config_cmfilter,    1},
		{"filter-noise",       0, &config_noisefilter, 1},
551
		{"no-match",           0, &config_nomatch,     1},
552
		{"verbose",            0, &config_verbose,     1},
553
		{"alternate",          0, &config_alternate,   1},
554
		{"intensities",        1, NULL,               'q'},
Thomas White's avatar
Thomas White committed
555
		{"pdb",                1, NULL,               'p'},
556
		{"prefix",             1, NULL,               'x'},
Thomas White's avatar
Thomas White committed
557
		{"unpolarized",        0, &config_polar,       0},
558
		{"check-sanity",       0, &config_sanity,      1},
559
		{"sat-corr",           0, &config_satcorr,     1},
560
		{"no-sa",              0, &config_sa,          0},
561
		{"threshold",          1, NULL,               't'},
Thomas White's avatar
Thomas White committed
562
563
564
565
		{0, 0, NULL, 0}
	};

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

		switch (c) {
Thomas White's avatar
Thomas White committed
570
		case 'h' :
Thomas White's avatar
Thomas White committed
571
572
573
			show_help(argv[0]);
			return 0;

Thomas White's avatar
Thomas White committed
574
		case 'i' :
Thomas White's avatar
Thomas White committed
575
576
577
			filename = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
578
		case 'z' :
Thomas White's avatar
Thomas White committed
579
580
581
			indm_str = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
582
		case 'q' :
583
584
585
			intfile = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
586
		case 'p' :
Thomas White's avatar
Thomas White committed
587
588
589
			pdb = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
590
		case 'x' :
591
592
593
			prefix = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
594
		case 'j' :
595
596
597
			nthreads = atoi(optarg);
			break;

Thomas White's avatar
Thomas White committed
598
		case 'g' :
Thomas White's avatar
Thomas White committed
599
600
601
			geometry = strdup(optarg);
			break;

602
603
604
605
		case 't' :
			threshold = strtof(optarg, NULL);
			break;

Thomas White's avatar
Thomas White committed
606
		case 0 :
Thomas White's avatar
Thomas White committed
607
608
			break;

Thomas White's avatar
Thomas White committed
609
		default :
Thomas White's avatar
Thomas White committed
610
611
612
613
614
615
616
617
618
619
620
621
622
623
			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
624
		ERROR("Failed to open input file '%s'\n", filename);
Thomas White's avatar
Thomas White committed
625
626
		return 1;
	}
Thomas White's avatar
Thomas White committed
627
	free(filename);
Thomas White's avatar
Thomas White committed
628

629
	if ( intfile != NULL ) {
Thomas White's avatar
Thomas White committed
630
631
632
		ReflItemList *items;
		items = read_reflections(intfile, intensities, NULL, NULL);
		delete_items(items);
633
634
635
636
	} else {
		intensities = NULL;
	}

Thomas White's avatar
Thomas White committed
637
638
639
640
	if ( pdb == NULL ) {
		pdb = strdup("molecule.pdb");
	}

641
	if ( prefix == NULL ) {
Thomas White's avatar
Thomas White committed
642
		prefix = strdup("");
643
644
	} else {
		prefix = check_prefix(prefix);
645
646
	}

647
	if ( (nthreads == 0) || (nthreads > MAX_THREADS) ) {
648
649
650
651
		ERROR("Invalid number of threads.\n");
		return 1;
	}

Thomas White's avatar
Thomas White committed
652
653
	if ( indm_str == NULL ) {
		STATUS("You didn't specify an indexing method, so I won't"
654
655
656
		       " 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
657
658
659
660
661
		indm = INDEXING_NONE;
	} else if ( strcmp(indm_str, "none") == 0 ) {
		indm = INDEXING_NONE;
	} else if ( strcmp(indm_str, "dirax") == 0) {
		indm = INDEXING_DIRAX;
662
663
	} else if ( strcmp(indm_str, "template") == 0) {
		indm = INDEXING_TEMPLATE;
Thomas White's avatar
Thomas White committed
664
665
666
667
	} else {
		ERROR("Unrecognised indexing method '%s'\n", indm_str);
		return 1;
	}
668
	free(indm_str);
Thomas White's avatar
Thomas White committed
669

Thomas White's avatar
Thomas White committed
670
671
672
673
674
675
676
677
678
679
680
681
	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);

682
	if ( (!config_nomatch) || (indm == INDEXING_TEMPLATE) ) {
683
684
		cell = load_cell_from_pdb(pdb);
		if ( cell == NULL ) {
685
686
687
			ERROR("Couldn't read unit cell (from %s)\n", pdb);
			return 1;
		}
688
689
690
	} else {
		STATUS("No cell needed because --no-match was used.\n");
		cell = NULL;
691
	}
692
	free(pdb);
693

694
695
696
697
698
699
700
701
	/* 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);
702
	ipriv = prepare_indexing(indm, cell, prepare_filename, det);
703
704
705
706
707
	if ( ipriv == NULL ) {
		ERROR("Failed to prepare indexing.\n");
		return 1;
	}

Thomas White's avatar
Thomas White committed
708
	gsl_set_error_handler_off();
Thomas White's avatar
Thomas White committed
709
710
	n_images = 0;
	n_hits = 0;
Thomas White's avatar
Thomas White committed
711
	n_sane = 0;
712

713
714
715
716
717
718
	for ( i=0; i<nthreads; i++ ) {
		worker_args[i] = malloc(sizeof(struct process_args));
		worker_args[i]->filename = malloc(1024);
		worker_active[i] = 0;
	}

719
	/* Start threads off */
720
	for ( i=0; i<nthreads; i++ ) {
Thomas White's avatar
Thomas White committed
721
722

		char line[1024];
723
724
		struct process_args *pargs;
		int r;
Thomas White's avatar
Thomas White committed
725

726
		pargs = worker_args[i];
Thomas White's avatar
Thomas White committed
727

728
729
730
731
732
733
734
		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
735
		chomp(line);
736
		snprintf(pargs->filename, 1023, "%s%s", prefix, line);
737

738
739
		n_images++;

740
		pargs->output_mutex = &output_mutex;
741
		pargs->gpu_mutex = &gpu_mutex;
742
		pthread_mutex_init(&pargs->control_mutex, NULL);
743
744
745
746
747
748
749
750
751
752
		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
753
		pargs->config_polar = config_polar;
754
		pargs->config_sanity = config_sanity;
755
		pargs->config_satcorr = config_satcorr;
756
		pargs->config_sa = config_sa;
757
		pargs->cell = cell;
Thomas White's avatar
Thomas White committed
758
		pargs->det = det;
759
		pargs->ipriv = ipriv;
760
761
762
		pargs->indm = indm;
		pargs->intensities = intensities;
		pargs->gctx = gctx;
763
		pargs->threshold = threshold;
Thomas White's avatar
Thomas White committed
764
		pargs->id = i;
765
766
767
768
769
		pthread_mutex_lock(&pargs->control_mutex);
		pargs->done = 0;
		pargs->start = 1;
		pargs->finish = 0;
		pthread_mutex_unlock(&pargs->control_mutex);
770

Thomas White's avatar
Thomas White committed
771
		worker_active[i] = 1;
772
		r = pthread_create(&workers[i], NULL, worker_thread, pargs);
773
		if ( r != 0 ) {
Thomas White's avatar
Thomas White committed
774
			worker_active[i] = 0;
775
776
777
778
779
			ERROR("Couldn't start thread %i\n", i);
		}

	}

780
	/* Keep threads busy until the end of the data */
781
782
783
784
785
786
787
788
	do {

		int i;

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

			char line[1024];
			struct process_args *pargs;
789
790
791
792
			int done;

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

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

797
			/* Has the thread finished yet? */
798
			pargs = worker_args[i];
799
800
801
802
			pthread_mutex_lock(&pargs->control_mutex);
			done = pargs->done;
			pthread_mutex_unlock(&pargs->control_mutex);
			if ( !done ) continue;
803

804
805
			/* Results will be processed after checking if
			 * there are any more images to process. */
806

807
			/* Get next filename */
808
			rval = fgets(line, 1023, fh);
809
810
811
			/* In this case, the result of the last file
			 * file will be processed when the thread is
			 * joined. */
812
			if ( rval == NULL ) break;
813
814
815
816
817

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

818
			chomp(line);
819
			snprintf(pargs->filename, 1023, "%s%s", prefix, line);
820
821

			n_images++;
822
823
824
825
826
827
828

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

829
		}
Thomas White's avatar
Thomas White committed
830

Thomas White's avatar
Thomas White committed
831
832
	} while ( rval != NULL );

833
	/* Join threads */
834
	for ( i=0; i<nthreads; i++ ) {
Thomas White's avatar
Thomas White committed
835

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

838
839
		/* Tell the thread to exit */
		struct process_args *pargs = worker_args[i];
Thomas White's avatar
Thomas White committed
840
		pthread_mutex_lock(&pargs->control_mutex);
841
		pargs->finish = 1;
Thomas White's avatar
Thomas White committed
842
		pthread_mutex_unlock(&pargs->control_mutex);
Thomas White's avatar
Thomas White committed
843

844
845
		/* Wait for it to join */
		pthread_join(workers[i], NULL);
Thomas White's avatar
Thomas White committed
846
847
		worker_active[i] = 0;

848
849
		n_hits += pargs->hit;
		n_sane += pargs->peaks_sane;
Thomas White's avatar
Thomas White committed
850

Thomas White's avatar
Thomas White committed
851
	free:
852
		if ( worker_args[i]->filename != NULL ) {
853
854
			free(worker_args[i]->filename);
		}
855
		free(worker_args[i]);
Thomas White's avatar
Thomas White committed
856

857
858
859
	}

	free(prefix);
Thomas White's avatar
Thomas White committed
860
861
	free(det->panels);
	free(det);
862
	free(cell);
Thomas White's avatar
Thomas White committed
863
864
865
	fclose(fh);

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

Thomas White's avatar
Thomas White committed
868
869
870
871
	if ( gctx != NULL ) {
		cleanup_gpu(gctx);
	}

Thomas White's avatar
Thomas White committed
872
873
	return 0;
}