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

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

81

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


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


172
static struct image *get_simage(struct image *template, int alternate)
173
{
174
	struct image *image;
175
	struct panel panels[2];
176

177
178
	image = malloc(sizeof(*image));

179
	/* Simulate a diffraction pattern */
180
181
182
	image->twotheta = NULL;
	image->data = NULL;
	image->det = template->det;
183
184
185
	image->flags = NULL;
	image->f0_available = 0;
	image->f0 = 1.0;
186
187

	/* View head-on (unit cell is tilted) */
188
189
190
191
	image->orientation.w = 1.0;
	image->orientation.x = 0.0;
	image->orientation.y = 0.0;
	image->orientation.z = 0.0;
192
193
194

	/* Detector geometry for the simulation
	 * - not necessarily the same as the original. */
195
196
	image->width = 1024;
	image->height = 1024;
Thomas White's avatar
Thomas White committed
197
	image->det->n_panels = 2;
198

199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
	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;
217
		panels[1].cy = 525.0;
218
219
220
		panels[1].clen = 56.7e-2;  /* 56.7 cm */
		panels[1].res = 13333.3;   /* 75 microns/pixel */

Thomas White's avatar
Thomas White committed
221
		image->det->panels = panels;
222

223
224
225
	} else {

		/* Copy pointer to old geometry */
Thomas White's avatar
Thomas White committed
226
		image->det->panels = template->det->panels;
227
228

	}
229

230
	image->lambda = ph_en_to_lambda(eV_to_J(1.8e3));
231
	image->features = template->features;
232
	image->filename = template->filename;
233
	image->indexed_cell = template->indexed_cell;
Thomas White's avatar
Thomas White committed
234
	image->f0 = template->f0;
235

Thomas White's avatar
Thomas White committed
236
237
238
239
	/* Prevent muppetry */
	image->hits = NULL;
	image->n_hits = 0;

240
241
	return image;
}
242
243


244
static void simulate_and_write(struct image *simage, struct gpu_context **gctx,
Thomas White's avatar
Thomas White committed
245
                               const double *intensities, UnitCell *cell)
246
{
247
	/* Set up GPU if necessary */
Thomas White's avatar
Thomas White committed
248
	if ( (gctx != NULL) && (*gctx == NULL) ) {
Thomas White's avatar
Thomas White committed
249
		*gctx = setup_gpu(0, simage, intensities);
250
251
	}

Thomas White's avatar
Thomas White committed
252
	if ( (gctx != NULL) && (*gctx != NULL) ) {
253
		get_diffraction_gpu(*gctx, simage, 24, 24, 40, cell);
254
	} else {
255
		get_diffraction(simage, 24, 24, 40,
Thomas White's avatar
Thomas White committed
256
		                intensities, NULL, cell, 0,
257
		                GRADIENT_MOSAIC);
258
	}
259
	record_image(simage, 0);
260

261
	hdf5_write("simulated.h5", simage->data, simage->width, simage->height,
262
263
264
265
		   H5T_NATIVE_FLOAT);
}


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

	image.features = NULL;
	image.data = NULL;
	image.indexed_cell = NULL;
Thomas White's avatar
Thomas White committed
294
	image.id = pargs->id;
295
	image.filename = filename;
Thomas White's avatar
Thomas White committed
296
297
	image.hits = NULL;
	image.n_hits = 0;
Thomas White's avatar
Thomas White committed
298
	image.det = pargs->det;
299

300
301
302
303
304
305
	/* 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;

306
	STATUS("Processing '%s'\n", image.filename);
307

308
309
	result.peaks_sane = 0;
	result.hit = 0;
310

311
312
	hdfile = hdfile_open(filename);
	if ( hdfile == NULL ) {
313
		return result;
314
315
	} else if ( hdfile_set_first_image(hdfile, "/") ) {
		ERROR("Couldn't select path\n");
316
		return result;
317
318
	}

319
	hdf5_read(hdfile, &image, pargs->config_satcorr);
320
321
322
323
324
325
326
327
328
329
330
331
332

	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 {
333
		memcpy(data_for_measurement, image.data, data_size);
334
335
336
	}

	/* Perform 'fine' peak search */
337
	search_peaks(&image, pargs->threshold);
Thomas White's avatar
Thomas White committed
338
339
340
341
342
343

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

344
345
	if ( image_feature_count(image.features) < 5 ) goto done;

346
	if ( config_dumpfound ) dump_peaks(&image, pargs->output_mutex);
347
348
349
350
351
352
353
354
355
356

	/* 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,
357
		              config_verbose, pargs->ipriv);
358
359
360
361
362
	}

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

Thomas White's avatar
Thomas White committed
363
	/* Sanity check */
364
	if ( pargs->config_sanity
365
	  && !peak_sanity_check(&image, image.indexed_cell, 0, 0.1) ) {
Thomas White's avatar
Thomas White committed
366
367
368
		STATUS("Failed peak sanity check.\n");
		goto done;
	} else {
369
		result.peaks_sane = 1;
Thomas White's avatar
Thomas White committed
370
371
	}

372
373
	/* Measure intensities if requested */
	if ( config_nearbragg ) {
374
		output_intensities(&image, image.indexed_cell,
Thomas White's avatar
Thomas White committed
375
		                   pargs->output_mutex, config_polar,
376
377
		                   pargs->config_sa, pargs->config_closer,
		                   0, 0.1);
378
379
	}

380
381
	simage = get_simage(&image, config_alternate);

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

	/* 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 */
401
	cell_free(image.indexed_cell);
402
403
404

done:
	free(image.data);
Thomas White's avatar
Thomas White committed
405
	free(image.flags);
406
	image_feature_list_free(image.features);
Thomas White's avatar
Thomas White committed
407
	free(image.hits);
408
409
	hdfile_close(hdfile);

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

	/* Long options */
	const struct option longopts[] = {
		{"help",               0, NULL,               'h'},
		{"input",              1, NULL,               'i'},
Thomas White's avatar
Thomas White committed
508
		{"gpu",                0, &config_gpu,         1},
509
		{"no-index",           0, &config_noindex,     1},
Thomas White's avatar
Thomas White committed
510
		{"dump-peaks",         0, &config_dumpfound,   1},
511
		{"near-bragg",         0, &config_nearbragg,   1},
Thomas White's avatar
Thomas White committed
512
513
		{"write-drx",          0, &config_writedrx,    1},
		{"indexing",           1, NULL,               'z'},
Thomas White's avatar
Thomas White committed
514
		{"geometry",           1, NULL,               'g'},
515
		{"simulate",           0, &config_simulate,    1},
516
517
		{"filter-cm",          0, &config_cmfilter,    1},
		{"filter-noise",       0, &config_noisefilter, 1},
518
		{"no-match",           0, &config_nomatch,     1},
519
		{"verbose",            0, &config_verbose,     1},
520
		{"alternate",          0, &config_alternate,   1},
521
		{"intensities",        1, NULL,               'q'},
Thomas White's avatar
Thomas White committed
522
		{"pdb",                1, NULL,               'p'},
523
		{"prefix",             1, NULL,               'x'},
Thomas White's avatar
Thomas White committed
524
		{"unpolarized",        0, &config_polar,       0},
525
		{"check-sanity",       0, &config_sanity,      1},
526
		{"sat-corr",           0, &config_satcorr,     1},
527
		{"no-sa",              0, &config_sa,          0},
528
		{"threshold",          1, NULL,               't'},
529
		{"no-check-prefix",    0, &config_checkprefix, 0},
530
		{"no-closer-peak",     0, &config_closer,      0},
Thomas White's avatar
Thomas White committed
531
532
533
534
		{0, 0, NULL, 0}
	};

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

		switch (c) {
Thomas White's avatar
Thomas White committed
539
		case 'h' :
Thomas White's avatar
Thomas White committed
540
541
542
			show_help(argv[0]);
			return 0;

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

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

Thomas White's avatar
Thomas White committed
551
		case 'q' :
552
553
554
			intfile = strdup(optarg);
			break;

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

Thomas White's avatar
Thomas White committed
559
		case 'x' :
560
561
562
			prefix = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
563
		case 'j' :
564
565
566
			nthreads = atoi(optarg);
			break;

Thomas White's avatar
Thomas White committed
567
		case 'g' :
Thomas White's avatar
Thomas White committed
568
569
570
			geometry = strdup(optarg);
			break;

571
572
573
574
		case 't' :
			threshold = strtof(optarg, NULL);
			break;

Thomas White's avatar
Thomas White committed
575
		case 0 :
Thomas White's avatar
Thomas White committed
576
577
			break;

Thomas White's avatar
Thomas White committed
578
		default :
Thomas White's avatar
Thomas White committed
579
580
581
582
583
			return 1;
		}

	}

584
585
586
587
588
589
590
	/* Start by writing the entire command line to stdout */
	printf("Command line:");
	for ( i=0; i<argc; i++ ) {
		printf(" %s", argv[i]);
	}
	printf("\n");

Thomas White's avatar
Thomas White committed
591
592
593
594
595
596
597
598
599
	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
600
		ERROR("Failed to open input file '%s'\n", filename);
Thomas White's avatar
Thomas White committed
601
602
		return 1;
	}
Thomas White's avatar
Thomas White committed
603
	free(filename);
Thomas White's avatar
Thomas White committed
604

605
	if ( intfile != NULL ) {
Thomas White's avatar
Thomas White committed
606
		ReflItemList *items;
607
608
		items = read_reflections(intfile, intensities,
		                         NULL, NULL, NULL);
Thomas White's avatar
Thomas White committed
609
		delete_items(items);
610
611
612
613
	} else {
		intensities = NULL;
	}

Thomas White's avatar
Thomas White committed
614
615
616
617
	if ( pdb == NULL ) {
		pdb = strdup("molecule.pdb");
	}

618
	if ( prefix == NULL ) {
Thomas White's avatar
Thomas White committed
619
		prefix = strdup("");
620
	} else {
621
622
623
		if ( config_checkprefix ) {
			prefix = check_prefix(prefix);
		}
624
625
	}

626
	if ( (nthreads == 0) || (nthreads > MAX_THREADS) ) {
627
628
629
630
		ERROR("Invalid number of threads.\n");
		return 1;
	}

Thomas White's avatar
Thomas White committed
631
632
	if ( indm_str == NULL ) {
		STATUS("You didn't specify an indexing method, so I won't"
633
634
635
		       " 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
636
637
638
639
640
		indm = INDEXING_NONE;
	} else if ( strcmp(indm_str, "none") == 0 ) {
		indm = INDEXING_NONE;
	} else if ( strcmp(indm_str, "dirax") == 0) {
		indm = INDEXING_DIRAX;
641
642
	} else if ( strcmp(indm_str, "template") == 0) {
		indm = INDEXING_TEMPLATE;
Thomas White's avatar
Thomas White committed
643
644
645
646
	} else {
		ERROR("Unrecognised indexing method '%s'\n", indm_str);
		return 1;
	}
647
	free(indm_str);
Thomas White's avatar
Thomas White committed
648

Thomas White's avatar
Thomas White committed
649
650
651
652
653
654
655
656
657
658
659
660
	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);

661
	if ( (!config_nomatch) || (indm == INDEXING_TEMPLATE) ) {
662
663
		cell = load_cell_from_pdb(pdb);
		if ( cell == NULL ) {
664
665
666
			ERROR("Couldn't read unit cell (from %s)\n", pdb);
			return 1;
		}
667
668
669
	} else {
		STATUS("No cell needed because --no-match was used.\n");
		cell = NULL;
670
	}
671
	free(pdb);
672

673
674
675
676
677
678
679
680
	/* 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);
681
	ipriv = prepare_indexing(indm, cell, prepare_filename, det);
682
683
684
685
686
	if ( ipriv == NULL ) {
		ERROR("Failed to prepare indexing.\n");
		return 1;
	}

Thomas White's avatar
Thomas White committed
687
	gsl_set_error_handler_off();
Thomas White's avatar
Thomas White committed
688
689
	n_images = 0;
	n_hits = 0;
Thomas White's avatar
Thomas White committed
690
	n_sane = 0;
691

692
693
694
695
696
697
	for ( i=0; i<nthreads; i++ ) {
		worker_args[i] = malloc(sizeof(struct process_args));
		worker_args[i]->filename = malloc(1024);
		worker_active[i] = 0;
	}

698
	/* Start threads off */
699
	for ( i=0; i<nthreads; i++ ) {
Thomas White's avatar
Thomas White committed
700
701

		char line[1024];
702
703
		struct process_args *pargs;
		int r;
Thomas White's avatar
Thomas White committed
704

705
		pargs = worker_args[i];
Thomas White's avatar
Thomas White committed
706

707
708
709
710
711
712
713
		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
714
		chomp(line);
715
		snprintf(pargs->filename, 1023, "%s%s", prefix, line);
716

717
718
		n_images++;

719
		pargs->output_mutex = &output_mutex;
720
		pargs->gpu_mutex = &gpu_mutex;
721
		pthread_mutex_init(&pargs->control_mutex, NULL);
722
723
724
725
726
727
728
729
730
731
		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
732
		pargs->config_polar = config_polar;
733
		pargs->config_sanity = config_sanity;
734
		pargs->config_satcorr = config_satcorr;
735
		pargs->config_sa = config_sa;
736
		pargs->config_closer = config_closer;
737
		pargs->cell = cell;
Thomas White's avatar
Thomas White committed
738
		pargs->det = det;
739
		pargs->ipriv = ipriv;
740
741
742
		pargs->indm = indm;
		pargs->intensities = intensities;
		pargs->gctx = gctx;
743
		pargs->threshold = threshold;
Thomas White's avatar
Thomas White committed
744
		pargs->id = i;
745
746
747
748
749
		pthread_mutex_lock(&pargs->control_mutex);
		pargs->done = 0;
		pargs->start = 1;
		pargs->finish = 0;
		pthread_mutex_unlock(&pargs->control_mutex);
750

Thomas White's avatar
Thomas White committed
751
		worker_active[i] = 1;
752
		r = pthread_create(&workers[i], NULL, worker_thread, pargs);
753
		if ( r != 0 ) {
Thomas White's avatar
Thomas White committed
754
			worker_active[i] = 0;
755
756
757
758
759
			ERROR("Couldn't start thread %i\n", i);
		}

	}

760
	/* Keep threads busy until the end of the data */
761
762
763
764
765
766
767
768
	do {

		int i;

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

			char line[1024];
			struct process_args *pargs;
769
770
771
772
			int done;

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

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

777
			/* Has the thread finished yet? */
778
			pargs = worker_args[i];
779
780
781
782
			pthread_mutex_lock(&pargs->control_mutex);
			done = pargs->done;
			pthread_mutex_unlock(&pargs->control_mutex);
			if ( !done ) continue;
783

784
785
			/* Results will be processed after checking if
			 * there are any more images to process. */
786

787
			/* Get next filename */
788
			rval = fgets(line, 1023, fh);
789
790
791
			/* In this case, the result of the last file
			 * file will be processed when the thread is
			 * joined. */
792
			if ( rval == NULL ) break;
793
794
795
796
797

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

798
			chomp(line);
799
			snprintf(pargs->filename, 1023, "%s%s", prefix, line);
800
801

			n_images++;
802
803
804
805
806
807
808

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

809
		}
Thomas White's avatar
Thomas White committed
810

Thomas White's avatar
Thomas White committed
811
812
	} while ( rval != NULL );

813
	/* Join threads */
814
	for ( i=0; i<nthreads; i++ ) {
Thomas White's avatar
Thomas White committed
815

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

818
819
		/* Tell the thread to exit */
		struct process_args *pargs = worker_args[i];
Thomas White's avatar
Thomas White committed
820
		pthread_mutex_lock(&pargs->control_mutex);
821
		pargs->finish = 1;
Thomas White's avatar
Thomas White committed
822
		pthread_mutex_unlock(&pargs->control_mutex);
Thomas White's avatar
Thomas White committed
823

824
825
		/* Wait for it to join */
		pthread_join(workers[i], NULL);
Thomas White's avatar
Thomas White committed
826
827
		worker_active[i] = 0;

828
829
		n_hits += pargs->hit;
		n_sane += pargs->peaks_sane;
Thomas White's avatar
Thomas White committed
830

Thomas White's avatar
Thomas White committed
831
	free:
832
		if ( worker_args[i]->filename != NULL ) {
833
834
			free(worker_args[i]->filename);
		}
835
		free(worker_args[i]);
Thomas White's avatar
Thomas White committed
836

837
838
	}

839
840
	cleanup_indexing(ipriv);

841
	free(prefix);
Thomas White's avatar
Thomas White committed
842
843
	free(det->panels);
	free(det);
844
	cell_free(cell);
Thomas White's avatar
Thomas White committed
845
846
847
	fclose(fh);

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

Thomas White's avatar
Thomas White committed
850
851
852
853
	if ( gctx != NULL ) {
		cleanup_gpu(gctx);
	}

Thomas White's avatar
Thomas White committed
854
855
	return 0;
}