indexamajig.c 22.5 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

	/* 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
77
78
	int indexable;
	int sane;
79
80
81
};


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


165
static struct image *get_simage(struct image *template, int alternate)
166
{
167
	struct image *image;
168
	struct panel panels[2];
169

170
171
	image = malloc(sizeof(*image));

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

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

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

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

Thomas White's avatar
Thomas White committed
214
		image->det->panels = panels;
215

216
217
218
	} else {

		/* Copy pointer to old geometry */
Thomas White's avatar
Thomas White committed
219
		image->det->panels = template->det->panels;
220
221

	}
222

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

Thomas White's avatar
Thomas White committed
229
230
231
232
	/* Prevent muppetry */
	image->hits = NULL;
	image->n_hits = 0;

233
234
	return image;
}
235
236


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

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

254
	hdf5_write("simulated.h5", simage->data, simage->width, simage->height,
255
256
257
258
		   H5T_NATIVE_FLOAT);
}


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

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

292
293
294
295
296
297
	/* 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;

298
	STATUS("Processing '%s'\n", image.filename);
299

Thomas White's avatar
Thomas White committed
300
301
	pargs->sane = 0;
	pargs->indexable = 0;
302

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

311
	hdf5_read(hdfile, &image, pargs->config_satcorr);
312
313
314
315
316
317
318
319
320
321
322
323
324

	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 {
325
		memcpy(data_for_measurement, image.data, data_size);
326
327
328
	}

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

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

336
	if ( config_dumpfound ) dump_peaks(&image, pargs->output_mutex);
337
338
339
340
341
342
343
344
345
346

	/* 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,
347
		              config_verbose, pargs->ipriv);
348
349
350
351
	}

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

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

363
364
	/* Measure intensities if requested */
	if ( config_nearbragg ) {
365
		output_intensities(&image, image.indexed_cell,
Thomas White's avatar
Thomas White committed
366
		                   pargs->output_mutex, config_polar,
367
368
		                   pargs->config_sa, pargs->config_closer,
		                   0, 0.1);
369
370
	}

371
372
	simage = get_simage(&image, config_alternate);

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

	/* 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 */
392
	cell_free(image.indexed_cell);
393
394
395

done:
	free(image.data);
Thomas White's avatar
Thomas White committed
396
	free(image.flags);
397
	image_feature_list_free(image.features);
Thomas White's avatar
Thomas White committed
398
	free(image.hits);
399
400
401
402
	hdfile_close(hdfile);
}


403
404
405
406
407
408
409
410
411
static void *worker_thread(void *pargsv)
{
	struct process_args *pargs = pargsv;
	int finish;

	do {

		int wakeup;

Thomas White's avatar
Thomas White committed
412
		process_image(pargs);
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436

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

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

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

		switch (c) {
Thomas White's avatar
Thomas White committed
521
		case 'h' :
Thomas White's avatar
Thomas White committed
522
523
524
			show_help(argv[0]);
			return 0;

Thomas White's avatar
Thomas White committed
525
		case 'i' :
Thomas White's avatar
Thomas White committed
526
527
528
			filename = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
529
		case 'z' :
Thomas White's avatar
Thomas White committed
530
531
532
			indm_str = strdup(optarg);
			break;

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

Thomas White's avatar
Thomas White committed
537
		case 'p' :
Thomas White's avatar
Thomas White committed
538
539
540
			pdb = strdup(optarg);
			break;

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

Thomas White's avatar
Thomas White committed
545
		case 'j' :
546
547
548
			nthreads = atoi(optarg);
			break;

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

553
554
555
556
		case 't' :
			threshold = strtof(optarg, NULL);
			break;

Thomas White's avatar
Thomas White committed
557
		case 0 :
Thomas White's avatar
Thomas White committed
558
559
			break;

Thomas White's avatar
Thomas White committed
560
		default :
Thomas White's avatar
Thomas White committed
561
562
563
564
565
			return 1;
		}

	}

566
567
568
569
570
571
572
	/* 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
573
574
575
576
577
578
579
580
581
	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
582
		ERROR("Failed to open input file '%s'\n", filename);
Thomas White's avatar
Thomas White committed
583
584
		return 1;
	}
Thomas White's avatar
Thomas White committed
585
	free(filename);
Thomas White's avatar
Thomas White committed
586

587
	if ( intfile != NULL ) {
Thomas White's avatar
Thomas White committed
588
		ReflItemList *items;
589
590
		items = read_reflections(intfile, intensities,
		                         NULL, NULL, NULL);
Thomas White's avatar
Thomas White committed
591
		delete_items(items);
592
593
594
595
	} else {
		intensities = NULL;
	}

Thomas White's avatar
Thomas White committed
596
597
598
599
	if ( pdb == NULL ) {
		pdb = strdup("molecule.pdb");
	}

600
	if ( prefix == NULL ) {
Thomas White's avatar
Thomas White committed
601
		prefix = strdup("");
602
	} else {
603
604
605
		if ( config_checkprefix ) {
			prefix = check_prefix(prefix);
		}
606
607
	}

608
	if ( (nthreads == 0) || (nthreads > MAX_THREADS) ) {
609
610
611
612
		ERROR("Invalid number of threads.\n");
		return 1;
	}

Thomas White's avatar
Thomas White committed
613
614
	if ( indm_str == NULL ) {
		STATUS("You didn't specify an indexing method, so I won't"
615
616
617
		       " 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
618
619
620
621
622
		indm = INDEXING_NONE;
	} else if ( strcmp(indm_str, "none") == 0 ) {
		indm = INDEXING_NONE;
	} else if ( strcmp(indm_str, "dirax") == 0) {
		indm = INDEXING_DIRAX;
623
624
	} else if ( strcmp(indm_str, "template") == 0) {
		indm = INDEXING_TEMPLATE;
Thomas White's avatar
Thomas White committed
625
626
627
628
	} else {
		ERROR("Unrecognised indexing method '%s'\n", indm_str);
		return 1;
	}
629
	free(indm_str);
Thomas White's avatar
Thomas White committed
630

Thomas White's avatar
Thomas White committed
631
632
633
634
635
636
637
638
639
640
641
642
	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);

643
	if ( (!config_nomatch) || (indm == INDEXING_TEMPLATE) ) {
644
645
		cell = load_cell_from_pdb(pdb);
		if ( cell == NULL ) {
646
647
648
			ERROR("Couldn't read unit cell (from %s)\n", pdb);
			return 1;
		}
649
650
651
	} else {
		STATUS("No cell needed because --no-match was used.\n");
		cell = NULL;
652
	}
653
	free(pdb);
654

655
656
657
658
659
660
661
662
	/* 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);
663
	ipriv = prepare_indexing(indm, cell, prepare_filename, det);
664
665
666
667
668
	if ( ipriv == NULL ) {
		ERROR("Failed to prepare indexing.\n");
		return 1;
	}

Thomas White's avatar
Thomas White committed
669
	gsl_set_error_handler_off();
Thomas White's avatar
Thomas White committed
670
	n_images = 0;
Thomas White's avatar
Thomas White committed
671
	n_indexable = 0;
Thomas White's avatar
Thomas White committed
672
	n_sane = 0;
673

674
675
676
677
678
679
	for ( i=0; i<nthreads; i++ ) {
		worker_args[i] = malloc(sizeof(struct process_args));
		worker_args[i]->filename = malloc(1024);
		worker_active[i] = 0;
	}

680
	/* Start threads off */
681
	for ( i=0; i<nthreads; i++ ) {
Thomas White's avatar
Thomas White committed
682
683

		char line[1024];
684
685
		struct process_args *pargs;
		int r;
Thomas White's avatar
Thomas White committed
686

687
		pargs = worker_args[i];
Thomas White's avatar
Thomas White committed
688

689
690
691
692
693
694
695
		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
696
		chomp(line);
697
		snprintf(pargs->filename, 1023, "%s%s", prefix, line);
698

699
700
		n_images++;

701
		pargs->output_mutex = &output_mutex;
702
		pargs->gpu_mutex = &gpu_mutex;
703
		pthread_mutex_init(&pargs->control_mutex, NULL);
704
705
706
707
708
709
710
711
712
713
		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
714
		pargs->config_polar = config_polar;
715
		pargs->config_sanity = config_sanity;
716
		pargs->config_satcorr = config_satcorr;
717
		pargs->config_sa = config_sa;
718
		pargs->config_closer = config_closer;
719
		pargs->cell = cell;
Thomas White's avatar
Thomas White committed
720
		pargs->det = det;
721
		pargs->ipriv = ipriv;
722
723
724
		pargs->indm = indm;
		pargs->intensities = intensities;
		pargs->gctx = gctx;
725
		pargs->threshold = threshold;
Thomas White's avatar
Thomas White committed
726
		pargs->id = i;
727
728
729
730
731
		pthread_mutex_lock(&pargs->control_mutex);
		pargs->done = 0;
		pargs->start = 1;
		pargs->finish = 0;
		pthread_mutex_unlock(&pargs->control_mutex);
732

Thomas White's avatar
Thomas White committed
733
		worker_active[i] = 1;
734
		r = pthread_create(&workers[i], NULL, worker_thread, pargs);
735
		if ( r != 0 ) {
Thomas White's avatar
Thomas White committed
736
			worker_active[i] = 0;
737
738
739
740
741
			ERROR("Couldn't start thread %i\n", i);
		}

	}

742
	/* Keep threads busy until the end of the data */
743
744
745
746
747
748
749
750
	do {

		int i;

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

			char line[1024];
			struct process_args *pargs;
751
752
753
754
			int done;

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

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

759
			/* Has the thread finished yet? */
760
			pargs = worker_args[i];
761
762
763
764
			pthread_mutex_lock(&pargs->control_mutex);
			done = pargs->done;
			pthread_mutex_unlock(&pargs->control_mutex);
			if ( !done ) continue;
765

766
767
			/* Results will be processed after checking if
			 * there are any more images to process. */
768

769
			/* Get next filename */
770
			rval = fgets(line, 1023, fh);
771
772
773
			/* In this case, the result of the last file
			 * file will be processed when the thread is
			 * joined. */
774
			if ( rval == NULL ) break;
775
776

			/* Record the result */
Thomas White's avatar
Thomas White committed
777
778
			n_indexable += pargs->indexable;
			n_sane += pargs->sane;
779

780
			chomp(line);
781
			snprintf(pargs->filename, 1023, "%s%s", prefix, line);
782
783

			n_images++;
784
785
786
787
788
789
790

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

791
		}
Thomas White's avatar
Thomas White committed
792

Thomas White's avatar
Thomas White committed
793
794
	} while ( rval != NULL );

795
	/* Join threads */
796
	for ( i=0; i<nthreads; i++ ) {
Thomas White's avatar
Thomas White committed
797

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

800
801
		/* Tell the thread to exit */
		struct process_args *pargs = worker_args[i];
Thomas White's avatar
Thomas White committed
802
		pthread_mutex_lock(&pargs->control_mutex);
803
		pargs->finish = 1;
Thomas White's avatar
Thomas White committed
804
		pthread_mutex_unlock(&pargs->control_mutex);
Thomas White's avatar
Thomas White committed
805

806
807
		/* Wait for it to join */
		pthread_join(workers[i], NULL);
Thomas White's avatar
Thomas White committed
808
809
		worker_active[i] = 0;

Thomas White's avatar
Thomas White committed
810
811
		n_indexable += pargs->indexable;
		n_sane += pargs->sane;
Thomas White's avatar
Thomas White committed
812

Thomas White's avatar
Thomas White committed
813
	free:
814
		if ( worker_args[i]->filename != NULL ) {
815
816
			free(worker_args[i]->filename);
		}
817
		free(worker_args[i]);
Thomas White's avatar
Thomas White committed
818

819
820
	}

821
822
	cleanup_indexing(ipriv);

823
	free(prefix);
Thomas White's avatar
Thomas White committed
824
825
	free(det->panels);
	free(det);
826
	cell_free(cell);
Thomas White's avatar
Thomas White committed
827
828
	fclose(fh);

Thomas White's avatar
Thomas White committed
829
830
	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
831

Thomas White's avatar
Thomas White committed
832
833
834
835
	if ( gctx != NULL ) {
		cleanup_gpu(gctx);
	}

Thomas White's avatar
Thomas White committed
836
837
	return 0;
}