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

	/* Thread control and output */
	pthread_mutex_t control_mutex;  /* Protects the scary stuff below */
	int start;
	int finish;
	int done;
Thomas White's avatar
Thomas White committed
76
77
	int indexable;
	int sane;
78
79
80
81

	/* Output stream */
	pthread_mutex_t *output_mutex;  /* Protects the output stream */
	FILE *ofh;
82
83
84
};


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


168
static struct image *get_simage(struct image *template, int alternate)
169
{
170
	struct image *image;
171
	struct panel panels[2];
172

173
174
	image = malloc(sizeof(*image));

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

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

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

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

Thomas White's avatar
Thomas White committed
217
		image->det->panels = panels;
218

219
220
221
	} else {

		/* Copy pointer to old geometry */
Thomas White's avatar
Thomas White committed
222
		image->det->panels = template->det->panels;
223
224

	}
225

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

Thomas White's avatar
Thomas White committed
232
233
234
235
	/* Prevent muppetry */
	image->hits = NULL;
	image->n_hits = 0;

236
237
	return image;
}
238
239


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

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

257
	hdf5_write("simulated.h5", simage->data, simage->width, simage->height,
258
259
260
261
		   H5T_NATIVE_FLOAT);
}


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

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

295
296
297
298
299
300
	/* 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;

301
	STATUS("Processing '%s'\n", image.filename);
302

Thomas White's avatar
Thomas White committed
303
304
	pargs->sane = 0;
	pargs->indexable = 0;
305

306
307
	hdfile = hdfile_open(filename);
	if ( hdfile == NULL ) {
Thomas White's avatar
Thomas White committed
308
		return;
309
310
	} else if ( hdfile_set_first_image(hdfile, "/") ) {
		ERROR("Couldn't select path\n");
Thomas White's avatar
Thomas White committed
311
		return;
312
313
	}

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

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

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

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

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

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

	/* No cell at this point?  Then we're done. */
	if ( image.indexed_cell == NULL ) goto done;
Thomas White's avatar
Thomas White committed
357
	pargs->indexable = 1;
358

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

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

376
377
	simage = get_simage(&image, config_alternate);

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

	/* Finished with alternate image */
	if ( simage->twotheta != NULL ) free(simage->twotheta);
	if ( simage->data != NULL ) free(simage->data);
	free(simage);

	/* Only free cell if found */
397
	cell_free(image.indexed_cell);
398
399
400

done:
	free(image.data);
Thomas White's avatar
Thomas White committed
401
	free(image.flags);
402
	image_feature_list_free(image.features);
Thomas White's avatar
Thomas White committed
403
	free(image.hits);
404
405
406
407
	hdfile_close(hdfile);
}


408
409
410
411
412
413
414
415
416
static void *worker_thread(void *pargsv)
{
	struct process_args *pargs = pargsv;
	int finish;

	do {

		int wakeup;

Thomas White's avatar
Thomas White committed
417
		process_image(pargs);
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441

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

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

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

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

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

537
538
539
540
		case 'o' :
			outfile = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
541
		case 'z' :
Thomas White's avatar
Thomas White committed
542
543
544
			indm_str = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
545
		case 'q' :
546
547
548
			intfile = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
549
		case 'p' :
Thomas White's avatar
Thomas White committed
550
551
552
			pdb = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
553
		case 'x' :
554
555
556
			prefix = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
557
		case 'j' :
558
559
560
			nthreads = atoi(optarg);
			break;

Thomas White's avatar
Thomas White committed
561
		case 'g' :
Thomas White's avatar
Thomas White committed
562
563
564
			geometry = strdup(optarg);
			break;

565
566
567
568
		case 't' :
			threshold = strtof(optarg, NULL);
			break;

Thomas White's avatar
Thomas White committed
569
		case 0 :
Thomas White's avatar
Thomas White committed
570
571
			break;

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

592
593
594
595
596
597
598
599
600
601
602
603
604
605
	if ( outfile == NULL ) {
		outfile = strdup("-");
	}
	if ( strcmp(outfile, "-") == 0 ) {
		ofh = stdout;
	} else {
		ofh = fopen(outfile, "w");
	}
	if ( ofh == NULL ) {
		ERROR("Failed to open output file '%s'\n", outfile);
		return 1;
	}
	free(outfile);

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

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

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

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

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

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

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

674
675
676
677
678
	/* Start by writing the entire command line to stdout */
	fprintf(ofh, "Command line:");
	for ( i=0; i<argc; i++ ) {
		fprintf(ofh, " %s", argv[i]);
	}
679
680
	fprintf(ofh, "\n");
	fflush(ofh);
681

682
683
684
685
686
687
688
689
	/* 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);
690
	ipriv = prepare_indexing(indm, cell, prepare_filename, det);
691
692
693
694
695
	if ( ipriv == NULL ) {
		ERROR("Failed to prepare indexing.\n");
		return 1;
	}

Thomas White's avatar
Thomas White committed
696
	gsl_set_error_handler_off();
Thomas White's avatar
Thomas White committed
697
	n_images = 0;
Thomas White's avatar
Thomas White committed
698
	n_indexable = 0;
Thomas White's avatar
Thomas White committed
699
	n_sane = 0;
700

701
702
703
	for ( i=0; i<nthreads; i++ ) {
		worker_args[i] = malloc(sizeof(struct process_args));
		worker_args[i]->filename = malloc(1024);
704
		worker_args[i]->ofh = ofh;
705
706
707
		worker_active[i] = 0;
	}

708
	/* Start threads off */
709
	for ( i=0; i<nthreads; i++ ) {
Thomas White's avatar
Thomas White committed
710
711

		char line[1024];
712
713
		struct process_args *pargs;
		int r;
Thomas White's avatar
Thomas White committed
714

715
		pargs = worker_args[i];
Thomas White's avatar
Thomas White committed
716

717
718
719
720
721
722
723
		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
724
		chomp(line);
725
		snprintf(pargs->filename, 1023, "%s%s", prefix, line);
726

727
728
		n_images++;

729
		pargs->output_mutex = &output_mutex;
730
		pargs->gpu_mutex = &gpu_mutex;
731
		pthread_mutex_init(&pargs->control_mutex, NULL);
732
733
734
735
736
737
738
739
740
741
		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
742
		pargs->config_polar = config_polar;
743
		pargs->config_sanity = config_sanity;
744
		pargs->config_satcorr = config_satcorr;
745
		pargs->config_sa = config_sa;
746
		pargs->config_closer = config_closer;
747
		pargs->cell = cell;
Thomas White's avatar
Thomas White committed
748
		pargs->det = det;
749
		pargs->ipriv = ipriv;
750
751
752
		pargs->indm = indm;
		pargs->intensities = intensities;
		pargs->gctx = gctx;
753
		pargs->threshold = threshold;
Thomas White's avatar
Thomas White committed
754
		pargs->id = i;
755
756
757
758
759
		pthread_mutex_lock(&pargs->control_mutex);
		pargs->done = 0;
		pargs->start = 1;
		pargs->finish = 0;
		pthread_mutex_unlock(&pargs->control_mutex);
760

Thomas White's avatar
Thomas White committed
761
		worker_active[i] = 1;
762
		r = pthread_create(&workers[i], NULL, worker_thread, pargs);
763
		if ( r != 0 ) {
Thomas White's avatar
Thomas White committed
764
			worker_active[i] = 0;
765
766
767
768
769
			ERROR("Couldn't start thread %i\n", i);
		}

	}

770
	/* Keep threads busy until the end of the data */
771
772
773
774
775
776
777
778
	do {

		int i;

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

			char line[1024];
			struct process_args *pargs;
779
780
781
782
			int done;

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

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

787
			/* Has the thread finished yet? */
788
			pargs = worker_args[i];
789
790
791
792
			pthread_mutex_lock(&pargs->control_mutex);
			done = pargs->done;
			pthread_mutex_unlock(&pargs->control_mutex);
			if ( !done ) continue;
793

794
795
			/* Results will be processed after checking if
			 * there are any more images to process. */
796

797
			/* Get next filename */
798
			rval = fgets(line, 1023, fh);
799
800
801
			/* In this case, the result of the last file
			 * file will be processed when the thread is
			 * joined. */
802
			if ( rval == NULL ) break;
803
804

			/* Record the result */
Thomas White's avatar
Thomas White committed
805
806
			n_indexable += pargs->indexable;
			n_sane += pargs->sane;
807

808
			chomp(line);
809
			snprintf(pargs->filename, 1023, "%s%s", prefix, line);
810
811

			n_images++;
812
813
814
815
816
817
818

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

819
		}
Thomas White's avatar
Thomas White committed
820

Thomas White's avatar
Thomas White committed
821
822
	} while ( rval != NULL );

823
	/* Join threads */
824
	for ( i=0; i<nthreads; i++ ) {
Thomas White's avatar
Thomas White committed
825

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

828
829
		/* Tell the thread to exit */
		struct process_args *pargs = worker_args[i];
Thomas White's avatar
Thomas White committed
830
		pthread_mutex_lock(&pargs->control_mutex);
831
		pargs->finish = 1;
Thomas White's avatar
Thomas White committed
832
		pthread_mutex_unlock(&pargs->control_mutex);
Thomas White's avatar
Thomas White committed
833

834
835
		/* Wait for it to join */
		pthread_join(workers[i], NULL);
Thomas White's avatar
Thomas White committed
836
837
		worker_active[i] = 0;

Thomas White's avatar
Thomas White committed
838
839
		n_indexable += pargs->indexable;
		n_sane += pargs->sane;
Thomas White's avatar
Thomas White committed
840

Thomas White's avatar
Thomas White committed
841
	free:
842
		if ( worker_args[i]->filename != NULL ) {
843
844
			free(worker_args[i]->filename);
		}
845
		free(worker_args[i]);
Thomas White's avatar
Thomas White committed
846

847
848
	}

849
850
	cleanup_indexing(ipriv);

851
	free(prefix);
Thomas White's avatar
Thomas White committed
852
853
	free(det->panels);
	free(det);
854
	cell_free(cell);
Thomas White's avatar
Thomas White committed
855
856
	fclose(fh);

Thomas White's avatar
Thomas White committed
857
858
	STATUS("There were %i images.  %i could be indexed, of which %i"
	       " looked sane.\n", n_images, n_indexable, n_sane);
Thomas White's avatar
Thomas White committed
859

Thomas White's avatar
Thomas White committed
860
861
862
863
	if ( gctx != NULL ) {
		cleanup_gpu(gctx);
	}

Thomas White's avatar
Thomas White committed
864
865
	return 0;
}