indexamajig.c 23 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
		}
	}

	/* 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 */
436
	cell_free(image.indexed_cell);
437
438
439

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
	int config_checkprefix = 1;
517
	float threshold = 800.0;
Thomas White's avatar
Thomas White committed
518
519
	struct detector *det;
	char *geometry = NULL;
Thomas White's avatar
Thomas White committed
520
521
	IndexingMethod indm;
	char *indm_str = NULL;
522
523
524
	UnitCell *cell;
	double *intensities = NULL;
	char *intfile = NULL;
Thomas White's avatar
Thomas White committed
525
	char *pdb = NULL;
526
	char *prefix = NULL;
527
528
529
	int nthreads = 1;
	pthread_t workers[MAX_THREADS];
	struct process_args *worker_args[MAX_THREADS];
Thomas White's avatar
Thomas White committed
530
	int worker_active[MAX_THREADS];
531
	int i;
532
	pthread_mutex_t output_mutex = PTHREAD_MUTEX_INITIALIZER;
533
	pthread_mutex_t gpu_mutex = PTHREAD_MUTEX_INITIALIZER;
534
535
536
	char prepare_line[1024];
	char prepare_filename[1024];
	IndexingPrivate *ipriv;
Thomas White's avatar
Thomas White committed
537
538
539
540
541

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

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

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

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

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

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

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

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

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

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

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

Thomas White's avatar
Thomas White committed
608
		case 0 :
Thomas White's avatar
Thomas White committed
609
610
			break;

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

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

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

643
	if ( prefix == NULL ) {
Thomas White's avatar
Thomas White committed
644
		prefix = strdup("");
645
	} else {
646
647
648
		if ( config_checkprefix ) {
			prefix = check_prefix(prefix);
		}
649
650
	}

651
	if ( (nthreads == 0) || (nthreads > MAX_THREADS) ) {
652
653
654
655
		ERROR("Invalid number of threads.\n");
		return 1;
	}

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

Thomas White's avatar
Thomas White committed
674
675
676
677
678
679
680
681
682
683
684
685
	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);

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

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

Thomas White's avatar
Thomas White committed
712
	gsl_set_error_handler_off();
Thomas White's avatar
Thomas White committed
713
714
	n_images = 0;
	n_hits = 0;
Thomas White's avatar
Thomas White committed
715
	n_sane = 0;
716

717
718
719
720
721
722
	for ( i=0; i<nthreads; i++ ) {
		worker_args[i] = malloc(sizeof(struct process_args));
		worker_args[i]->filename = malloc(1024);
		worker_active[i] = 0;
	}

723
	/* Start threads off */
724
	for ( i=0; i<nthreads; i++ ) {
Thomas White's avatar
Thomas White committed
725
726

		char line[1024];
727
728
		struct process_args *pargs;
		int r;
Thomas White's avatar
Thomas White committed
729

730
		pargs = worker_args[i];
Thomas White's avatar
Thomas White committed
731

732
733
734
735
736
737
738
		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
739
		chomp(line);
740
		snprintf(pargs->filename, 1023, "%s%s", prefix, line);
741

742
743
		n_images++;

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

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

	}

784
	/* Keep threads busy until the end of the data */
785
786
787
788
789
790
791
792
	do {

		int i;

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

			char line[1024];
			struct process_args *pargs;
793
794
795
796
			int done;

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

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

801
			/* Has the thread finished yet? */
802
			pargs = worker_args[i];
803
804
805
806
			pthread_mutex_lock(&pargs->control_mutex);
			done = pargs->done;
			pthread_mutex_unlock(&pargs->control_mutex);
			if ( !done ) continue;
807

808
809
			/* Results will be processed after checking if
			 * there are any more images to process. */
810

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

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

822
			chomp(line);
823
			snprintf(pargs->filename, 1023, "%s%s", prefix, line);
824
825

			n_images++;
826
827
828
829
830
831
832

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

833
		}
Thomas White's avatar
Thomas White committed
834

Thomas White's avatar
Thomas White committed
835
836
	} while ( rval != NULL );

837
	/* Join threads */
838
	for ( i=0; i<nthreads; i++ ) {
Thomas White's avatar
Thomas White committed
839

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

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

848
849
		/* Wait for it to join */
		pthread_join(workers[i], NULL);
Thomas White's avatar
Thomas White committed
850
851
		worker_active[i] = 0;

852
853
		n_hits += pargs->hit;
		n_sane += pargs->peaks_sane;
Thomas White's avatar
Thomas White committed
854

Thomas White's avatar
Thomas White committed
855
	free:
856
		if ( worker_args[i]->filename != NULL ) {
857
858
			free(worker_args[i]->filename);
		}
859
		free(worker_args[i]);
Thomas White's avatar
Thomas White committed
860

861
862
	}

863
864
	cleanup_indexing(ipriv);

865
	free(prefix);
Thomas White's avatar
Thomas White committed
866
867
	free(det->panels);
	free(det);
868
	cell_free(cell);
Thomas White's avatar
Thomas White committed
869
870
871
	fclose(fh);

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

Thomas White's avatar
Thomas White committed
874
875
876
877
	if ( gctx != NULL ) {
		cleanup_gpu(gctx);
	}

Thomas White's avatar
Thomas White committed
878
879
	return 0;
}