indexamajig.c 24.2 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
#define MAX_THREADS (96)

42
43
44
45
46
47
48

enum {
	PEAK_ZAEF,
	PEAK_HDF5,
};


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

	/* 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
84
85
	int indexable;
	int sane;
86
87
88
89

	/* Output stream */
	pthread_mutex_t *output_mutex;  /* Protects the output stream */
	FILE *ofh;
90
91
92
};


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


182
static struct image *get_simage(struct image *template, int alternate)
183
{
184
	struct image *image;
185
	struct panel panels[2];
186

187
188
	image = malloc(sizeof(*image));

189
	/* Simulate a diffraction pattern */
190
191
192
	image->twotheta = NULL;
	image->data = NULL;
	image->det = template->det;
193
194
195
	image->flags = NULL;
	image->f0_available = 0;
	image->f0 = 1.0;
196
197

	/* View head-on (unit cell is tilted) */
198
199
200
201
	image->orientation.w = 1.0;
	image->orientation.x = 0.0;
	image->orientation.y = 0.0;
	image->orientation.z = 0.0;
202
203
204

	/* Detector geometry for the simulation
	 * - not necessarily the same as the original. */
205
206
	image->width = 1024;
	image->height = 1024;
Thomas White's avatar
Thomas White committed
207
	image->det->n_panels = 2;
208

209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
	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;
227
		panels[1].cy = 525.0;
228
229
230
		panels[1].clen = 56.7e-2;  /* 56.7 cm */
		panels[1].res = 13333.3;   /* 75 microns/pixel */

Thomas White's avatar
Thomas White committed
231
		image->det->panels = panels;
232

233
234
235
	} else {

		/* Copy pointer to old geometry */
Thomas White's avatar
Thomas White committed
236
		image->det->panels = template->det->panels;
237
238

	}
239

240
	image->lambda = ph_en_to_lambda(eV_to_J(1.8e3));
241
	image->features = template->features;
242
	image->filename = template->filename;
243
	image->indexed_cell = template->indexed_cell;
Thomas White's avatar
Thomas White committed
244
	image->f0 = template->f0;
245

Thomas White's avatar
Thomas White committed
246
247
248
249
	/* Prevent muppetry */
	image->hits = NULL;
	image->n_hits = 0;

250
251
	return image;
}
252
253


254
static void simulate_and_write(struct image *simage, struct gpu_context **gctx,
Thomas White's avatar
Thomas White committed
255
                               const double *intensities, UnitCell *cell)
256
{
257
	/* Set up GPU if necessary */
Thomas White's avatar
Thomas White committed
258
	if ( (gctx != NULL) && (*gctx == NULL) ) {
Thomas White's avatar
Thomas White committed
259
		*gctx = setup_gpu(0, simage, intensities);
260
261
	}

Thomas White's avatar
Thomas White committed
262
	if ( (gctx != NULL) && (*gctx != NULL) ) {
263
		get_diffraction_gpu(*gctx, simage, 24, 24, 40, cell);
264
	} else {
265
		get_diffraction(simage, 24, 24, 40,
Thomas White's avatar
Thomas White committed
266
		                intensities, NULL, cell, 0,
267
		                GRADIENT_MOSAIC);
268
	}
269
	record_image(simage, 0);
270

271
	hdf5_write("simulated.h5", simage->data, simage->width, simage->height,
272
273
274
275
		   H5T_NATIVE_FLOAT);
}


Thomas White's avatar
Thomas White committed
276
static void process_image(struct process_args *pargs)
277
278
279
280
281
282
{
	struct hdfile *hdfile;
	struct image image;
	struct image *simage;
	float *data_for_measurement;
	size_t data_size;
283
284
285
286
287
288
289
290
291
292
293
294
	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
295
	int config_polar = pargs->config_polar;
296
297
298
	IndexingMethod indm = pargs->indm;
	const double *intensities = pargs->intensities;
	struct gpu_context *gctx = pargs->gctx;
299
300
301
302

	image.features = NULL;
	image.data = NULL;
	image.indexed_cell = NULL;
Thomas White's avatar
Thomas White committed
303
	image.id = pargs->id;
304
	image.filename = filename;
Thomas White's avatar
Thomas White committed
305
306
	image.hits = NULL;
	image.n_hits = 0;
Thomas White's avatar
Thomas White committed
307
	image.det = pargs->det;
308

309
310
311
312
313
314
	/* 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;

315
	STATUS("Processing '%s'\n", image.filename);
316

Thomas White's avatar
Thomas White committed
317
318
	pargs->sane = 0;
	pargs->indexable = 0;
319

320
321
	hdfile = hdfile_open(filename);
	if ( hdfile == NULL ) {
Thomas White's avatar
Thomas White committed
322
		return;
323
324
	} else if ( hdfile_set_first_image(hdfile, "/") ) {
		ERROR("Couldn't select path\n");
Thomas White's avatar
Thomas White committed
325
		return;
326
327
	}

328
	hdf5_read(hdfile, &image, pargs->config_satcorr);
329
330
331
332
333
334
335
336
337
338
339
340
341

	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 {
342
		memcpy(data_for_measurement, image.data, data_size);
343
344
	}

345
346
347
348
349
350
351
352
353
354
355
356
357
	switch ( pargs->peaks )
	{
	case PEAK_HDF5 :
		/* Get peaks from HDF5 */
		if ( get_peaks(&image, hdfile) ) {
			ERROR("Failed to get peaks from HDF5 file.\n");
			return;
		}
		break;
	case PEAK_ZAEF :
		search_peaks(&image, pargs->threshold);
		break;
	}
Thomas White's avatar
Thomas White committed
358
359
360
361
362
363

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

364
365
366
	if ( config_dumpfound ) {
		dump_peaks(&image, pargs->ofh, pargs->output_mutex);
	}
367
368
369
370
371
372
373
374
375
376

	/* 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,
377
		              config_verbose, pargs->ipriv);
378
379
380
381
	}

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

Thomas White's avatar
Thomas White committed
384
	/* Sanity check */
385
	if ( pargs->config_sanity
386
	  && !peak_sanity_check(&image, image.indexed_cell, 0, 0.1) ) {
Thomas White's avatar
Thomas White committed
387
388
389
		STATUS("Failed peak sanity check.\n");
		goto done;
	} else {
Thomas White's avatar
Thomas White committed
390
		pargs->sane = 1;
Thomas White's avatar
Thomas White committed
391
392
	}

393
394
	/* Measure intensities if requested */
	if ( config_nearbragg ) {
395
		output_intensities(&image, image.indexed_cell,
Thomas White's avatar
Thomas White committed
396
		                   pargs->output_mutex, config_polar,
397
		                   pargs->config_sa, pargs->config_closer,
398
		                   pargs->ofh, 0, 0.1);
399
400
	}

401
402
	simage = get_simage(&image, config_alternate);

403
404
405
	/* Simulate if requested */
	if ( config_simulate ) {
		if ( config_gpu ) {
406
			pthread_mutex_lock(pargs->gpu_mutex);
407
			simulate_and_write(simage, &gctx, intensities,
Thomas White's avatar
Thomas White committed
408
			                   image.indexed_cell);
409
			pthread_mutex_unlock(pargs->gpu_mutex);
410
411
		} else {
			simulate_and_write(simage, NULL, intensities,
Thomas White's avatar
Thomas White committed
412
			                   image.indexed_cell);
413
414
415
416
417
418
419
420
421
		}
	}

	/* 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 */
422
	cell_free(image.indexed_cell);
423
424
425

done:
	free(image.data);
Thomas White's avatar
Thomas White committed
426
	free(image.flags);
427
	image_feature_list_free(image.features);
Thomas White's avatar
Thomas White committed
428
	free(image.hits);
429
430
431
432
	hdfile_close(hdfile);
}


433
434
435
436
437
438
439
440
441
static void *worker_thread(void *pargsv)
{
	struct process_args *pargs = pargsv;
	int finish;

	do {

		int wakeup;

Thomas White's avatar
Thomas White committed
442
		process_image(pargs);
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466

		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
467
468
469
int main(int argc, char *argv[])
{
	int c;
Thomas White's avatar
Thomas White committed
470
	struct gpu_context *gctx = NULL;
Thomas White's avatar
Thomas White committed
471
	char *filename = NULL;
472
	char *outfile = NULL;
Thomas White's avatar
Thomas White committed
473
	FILE *fh;
474
	FILE *ofh;
475
	char *rval = NULL;
Thomas White's avatar
Thomas White committed
476
	int n_images;
Thomas White's avatar
Thomas White committed
477
	int n_indexable;
Thomas White's avatar
Thomas White committed
478
	int n_sane;
479
	int config_noindex = 0;
Thomas White's avatar
Thomas White committed
480
	int config_dumpfound = 0;
481
	int config_nearbragg = 0;
Thomas White's avatar
Thomas White committed
482
	int config_writedrx = 0;
483
	int config_simulate = 0;
484
485
	int config_cmfilter = 0;
	int config_noisefilter = 0;
486
	int config_nomatch = 0;
Thomas White's avatar
Thomas White committed
487
	int config_gpu = 0;
488
	int config_verbose = 0;
489
	int config_alternate = 0;
Thomas White's avatar
Thomas White committed
490
	int config_polar = 1;
491
	int config_sanity = 0;
492
	int config_satcorr = 1;
493
	int config_sa = 1;
494
	int config_checkprefix = 1;
495
	int config_closer = 1;
496
	float threshold = 800.0;
Thomas White's avatar
Thomas White committed
497
498
	struct detector *det;
	char *geometry = NULL;
Thomas White's avatar
Thomas White committed
499
500
	IndexingMethod indm;
	char *indm_str = NULL;
501
502
503
	UnitCell *cell;
	double *intensities = NULL;
	char *intfile = NULL;
Thomas White's avatar
Thomas White committed
504
	char *pdb = NULL;
505
	char *prefix = NULL;
506
507
	char *speaks = NULL;
	int peaks;
508
509
510
	int nthreads = 1;
	pthread_t workers[MAX_THREADS];
	struct process_args *worker_args[MAX_THREADS];
Thomas White's avatar
Thomas White committed
511
	int worker_active[MAX_THREADS];
512
	int i;
513
	pthread_mutex_t output_mutex = PTHREAD_MUTEX_INITIALIZER;
514
	pthread_mutex_t gpu_mutex = PTHREAD_MUTEX_INITIALIZER;
515
516
517
	char prepare_line[1024];
	char prepare_filename[1024];
	IndexingPrivate *ipriv;
Thomas White's avatar
Thomas White committed
518
519
520
521
522

	/* Long options */
	const struct option longopts[] = {
		{"help",               0, NULL,               'h'},
		{"input",              1, NULL,               'i'},
523
		{"output",             1, NULL,               'o'},
Thomas White's avatar
Thomas White committed
524
		{"gpu",                0, &config_gpu,         1},
525
		{"no-index",           0, &config_noindex,     1},
Thomas White's avatar
Thomas White committed
526
		{"dump-peaks",         0, &config_dumpfound,   1},
527
		{"peaks",              1, NULL,                2},
528
		{"near-bragg",         0, &config_nearbragg,   1},
Thomas White's avatar
Thomas White committed
529
530
		{"write-drx",          0, &config_writedrx,    1},
		{"indexing",           1, NULL,               'z'},
Thomas White's avatar
Thomas White committed
531
		{"geometry",           1, NULL,               'g'},
532
		{"simulate",           0, &config_simulate,    1},
533
534
		{"filter-cm",          0, &config_cmfilter,    1},
		{"filter-noise",       0, &config_noisefilter, 1},
535
		{"no-match",           0, &config_nomatch,     1},
536
		{"verbose",            0, &config_verbose,     1},
537
		{"alternate",          0, &config_alternate,   1},
538
		{"intensities",        1, NULL,               'q'},
Thomas White's avatar
Thomas White committed
539
		{"pdb",                1, NULL,               'p'},
540
		{"prefix",             1, NULL,               'x'},
Thomas White's avatar
Thomas White committed
541
		{"unpolarized",        0, &config_polar,       0},
542
		{"check-sanity",       0, &config_sanity,      1},
543
544
		{"no-sat-corr",        0, &config_satcorr,     0},
		{"sat-corr",           0, &config_satcorr,     1}, /* Compat */
545
		{"no-sa",              0, &config_sa,          0},
546
		{"threshold",          1, NULL,               't'},
547
		{"no-check-prefix",    0, &config_checkprefix, 0},
548
		{"no-closer-peak",     0, &config_closer,      0},
Thomas White's avatar
Thomas White committed
549
550
551
552
		{0, 0, NULL, 0}
	};

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

		switch (c) {
Thomas White's avatar
Thomas White committed
557
		case 'h' :
Thomas White's avatar
Thomas White committed
558
559
560
			show_help(argv[0]);
			return 0;

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

565
566
567
568
		case 'o' :
			outfile = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
569
		case 'z' :
Thomas White's avatar
Thomas White committed
570
571
572
			indm_str = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
573
		case 'q' :
574
575
576
			intfile = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
577
		case 'p' :
Thomas White's avatar
Thomas White committed
578
579
580
			pdb = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
581
		case 'x' :
582
583
584
			prefix = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
585
		case 'j' :
586
587
588
			nthreads = atoi(optarg);
			break;

Thomas White's avatar
Thomas White committed
589
		case 'g' :
Thomas White's avatar
Thomas White committed
590
591
592
			geometry = strdup(optarg);
			break;

593
594
595
596
		case 't' :
			threshold = strtof(optarg, NULL);
			break;

597
598
599
600
		case 2 :
			speaks = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
601
		case 0 :
Thomas White's avatar
Thomas White committed
602
603
			break;

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

624
625
626
627
628
629
630
631
632
633
634
635
636
637
	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);

638
639
640
641
642
643
644
645
646
647
648
649
650
651
	if ( speaks == NULL ) {
		speaks = strdup("zaef");
		STATUS("You didn't specify a peak detection method.\n");
		STATUS("I'm using 'zaef' for you.\n");
	}
	if ( strcmp(speaks, "zaef") == 0 ) {
		peaks = PEAK_ZAEF;
	} else if ( strcmp(speaks, "hdf5") == 0 ) {
		peaks = PEAK_HDF5;
	} else {
		ERROR("Unrecognised peak detection method '%s'\n", speaks);
		return 1;
	}

652
	if ( intfile != NULL ) {
Thomas White's avatar
Thomas White committed
653
		ReflItemList *items;
654
655
		items = read_reflections(intfile, intensities,
		                         NULL, NULL, NULL);
Thomas White's avatar
Thomas White committed
656
		delete_items(items);
657
658
659
660
	} else {
		intensities = NULL;
	}

Thomas White's avatar
Thomas White committed
661
662
663
664
	if ( pdb == NULL ) {
		pdb = strdup("molecule.pdb");
	}

665
	if ( prefix == NULL ) {
Thomas White's avatar
Thomas White committed
666
		prefix = strdup("");
667
	} else {
668
669
670
		if ( config_checkprefix ) {
			prefix = check_prefix(prefix);
		}
671
672
	}

673
	if ( (nthreads == 0) || (nthreads > MAX_THREADS) ) {
674
675
676
677
		ERROR("Invalid number of threads.\n");
		return 1;
	}

Thomas White's avatar
Thomas White committed
678
679
	if ( indm_str == NULL ) {
		STATUS("You didn't specify an indexing method, so I won't"
680
681
682
		       " 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
683
684
685
686
687
		indm = INDEXING_NONE;
	} else if ( strcmp(indm_str, "none") == 0 ) {
		indm = INDEXING_NONE;
	} else if ( strcmp(indm_str, "dirax") == 0) {
		indm = INDEXING_DIRAX;
688
689
	} else if ( strcmp(indm_str, "template") == 0) {
		indm = INDEXING_TEMPLATE;
Thomas White's avatar
Thomas White committed
690
691
692
693
	} else {
		ERROR("Unrecognised indexing method '%s'\n", indm_str);
		return 1;
	}
694
	free(indm_str);
Thomas White's avatar
Thomas White committed
695

Thomas White's avatar
Thomas White committed
696
697
698
699
700
701
702
703
704
705
706
707
	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);

708
	if ( (!config_nomatch) || (indm == INDEXING_TEMPLATE) ) {
709
710
		cell = load_cell_from_pdb(pdb);
		if ( cell == NULL ) {
711
712
713
			ERROR("Couldn't read unit cell (from %s)\n", pdb);
			return 1;
		}
714
715
716
	} else {
		STATUS("No cell needed because --no-match was used.\n");
		cell = NULL;
717
	}
718
	free(pdb);
719

720
721
722
723
724
	/* Start by writing the entire command line to stdout */
	fprintf(ofh, "Command line:");
	for ( i=0; i<argc; i++ ) {
		fprintf(ofh, " %s", argv[i]);
	}
725
726
	fprintf(ofh, "\n");
	fflush(ofh);
727

728
729
730
731
732
733
734
735
	/* 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);
736
	ipriv = prepare_indexing(indm, cell, prepare_filename, det);
737
738
739
740
741
	if ( ipriv == NULL ) {
		ERROR("Failed to prepare indexing.\n");
		return 1;
	}

Thomas White's avatar
Thomas White committed
742
	gsl_set_error_handler_off();
Thomas White's avatar
Thomas White committed
743
	n_images = 0;
Thomas White's avatar
Thomas White committed
744
	n_indexable = 0;
Thomas White's avatar
Thomas White committed
745
	n_sane = 0;
746

747
748
749
	for ( i=0; i<nthreads; i++ ) {
		worker_args[i] = malloc(sizeof(struct process_args));
		worker_args[i]->filename = malloc(1024);
750
		worker_args[i]->ofh = ofh;
751
		worker_args[i]->peaks = peaks;
752
753
754
		worker_active[i] = 0;
	}

755
	/* Start threads off */
756
	for ( i=0; i<nthreads; i++ ) {
Thomas White's avatar
Thomas White committed
757
758

		char line[1024];
759
760
		struct process_args *pargs;
		int r;
Thomas White's avatar
Thomas White committed
761

762
		pargs = worker_args[i];
Thomas White's avatar
Thomas White committed
763

764
765
766
767
768
769
770
		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
771
		chomp(line);
772
		snprintf(pargs->filename, 1023, "%s%s", prefix, line);
773

774
775
		n_images++;

776
		pargs->output_mutex = &output_mutex;
777
		pargs->gpu_mutex = &gpu_mutex;
778
		pthread_mutex_init(&pargs->control_mutex, NULL);
779
780
781
782
783
784
785
786
787
788
		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
789
		pargs->config_polar = config_polar;
790
		pargs->config_sanity = config_sanity;
791
		pargs->config_satcorr = config_satcorr;
792
		pargs->config_sa = config_sa;
793
		pargs->config_closer = config_closer;
794
		pargs->cell = cell;
Thomas White's avatar
Thomas White committed
795
		pargs->det = det;
796
		pargs->ipriv = ipriv;
797
798
799
		pargs->indm = indm;
		pargs->intensities = intensities;
		pargs->gctx = gctx;
800
		pargs->threshold = threshold;
Thomas White's avatar
Thomas White committed
801
		pargs->id = i;
802
803
804
805
806
		pthread_mutex_lock(&pargs->control_mutex);
		pargs->done = 0;
		pargs->start = 1;
		pargs->finish = 0;
		pthread_mutex_unlock(&pargs->control_mutex);
807

Thomas White's avatar
Thomas White committed
808
		worker_active[i] = 1;
809
		r = pthread_create(&workers[i], NULL, worker_thread, pargs);
810
		if ( r != 0 ) {
Thomas White's avatar
Thomas White committed
811
			worker_active[i] = 0;
812
813
814
815
816
			ERROR("Couldn't start thread %i\n", i);
		}

	}

817
	/* Keep threads busy until the end of the data */
818
819
820
821
822
823
824
825
	do {

		int i;

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

			char line[1024];
			struct process_args *pargs;
826
827
828
829
			int done;

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

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

834
			/* Has the thread finished yet? */
835
			pargs = worker_args[i];
836
837
838
839
			pthread_mutex_lock(&pargs->control_mutex);
			done = pargs->done;
			pthread_mutex_unlock(&pargs->control_mutex);
			if ( !done ) continue;
840

841
842
			/* Results will be processed after checking if
			 * there are any more images to process. */
843

844
			/* Get next filename */
845
			rval = fgets(line, 1023, fh);
846
847
848
			/* In this case, the result of the last file
			 * file will be processed when the thread is
			 * joined. */
849
			if ( rval == NULL ) break;
850
851

			/* Record the result */
Thomas White's avatar
Thomas White committed
852
853
			n_indexable += pargs->indexable;
			n_sane += pargs->sane;
854

855
			chomp(line);
856
			snprintf(pargs->filename, 1023, "%s%s", prefix, line);
857
858

			n_images++;
859
860
861
862
863
864
865

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

866
		}
Thomas White's avatar
Thomas White committed
867

Thomas White's avatar
Thomas White committed
868
869
	} while ( rval != NULL );

870
	/* Join threads */
871
	for ( i=0; i<nthreads; i++ ) {
Thomas White's avatar
Thomas White committed
872

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

875
876
		/* Tell the thread to exit */
		struct process_args *pargs = worker_args[i];
Thomas White's avatar
Thomas White committed
877
		pthread_mutex_lock(&pargs->control_mutex);
878
		pargs->finish = 1;
Thomas White's avatar
Thomas White committed
879
		pthread_mutex_unlock(&pargs->control_mutex);
Thomas White's avatar
Thomas White committed
880

881
882
		/* Wait for it to join */
		pthread_join(workers[i], NULL);
Thomas White's avatar
Thomas White committed
883
884
		worker_active[i] = 0;

Thomas White's avatar
Thomas White committed
885
886
		n_indexable += pargs->indexable;
		n_sane += pargs->sane;
Thomas White's avatar
Thomas White committed
887

Thomas White's avatar
Thomas White committed
888
	free:
889
		if ( worker_args[i]->filename != NULL ) {
890
891
			free(worker_args[i]->filename);
		}
892
		free(worker_args[i]);
Thomas White's avatar
Thomas White committed
893

894
895
	}

896
897
	cleanup_indexing(ipriv);

898
	free(prefix);
Thomas White's avatar
Thomas White committed
899
900
	free(det->panels);
	free(det);
901
	cell_free(cell);
Thomas White's avatar
Thomas White committed
902
903
	fclose(fh);

Thomas White's avatar
Thomas White committed
904
905
	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
906

Thomas White's avatar
Thomas White committed
907
908
909
910
	if ( gctx != NULL ) {
		cleanup_gpu(gctx);
	}

Thomas White's avatar
Thomas White committed
911
912
	return 0;
}