indexamajig.c 21.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;
Thomas White's avatar
Thomas White committed
64
	struct detector *det;
65
	IndexingMethod indm;
66
	IndexingPrivate *ipriv;
67
68
	const double *intensities;
	struct gpu_context *gctx;
69
70
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;
	int hit;
	int peaks_sane;
77
78
};

79

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


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


158
static struct image *get_simage(struct image *template, int alternate)
159
{
160
	struct image *image;
161
	struct panel panels[2];
162

163
164
	image = malloc(sizeof(*image));

165
	/* Simulate a diffraction pattern */
166
167
168
	image->twotheta = NULL;
	image->data = NULL;
	image->det = template->det;
169
170
171
	image->flags = NULL;
	image->f0_available = 0;
	image->f0 = 1.0;
172
173

	/* View head-on (unit cell is tilted) */
174
175
176
177
	image->orientation.w = 1.0;
	image->orientation.x = 0.0;
	image->orientation.y = 0.0;
	image->orientation.z = 0.0;
178
179
180

	/* Detector geometry for the simulation
	 * - not necessarily the same as the original. */
181
182
	image->width = 1024;
	image->height = 1024;
Thomas White's avatar
Thomas White committed
183
	image->det->n_panels = 2;
184

185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
	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;
203
		panels[1].cy = 525.0;
204
205
206
		panels[1].clen = 56.7e-2;  /* 56.7 cm */
		panels[1].res = 13333.3;   /* 75 microns/pixel */

Thomas White's avatar
Thomas White committed
207
		image->det->panels = panels;
208

209
210
211
	} else {

		/* Copy pointer to old geometry */
Thomas White's avatar
Thomas White committed
212
		image->det->panels = template->det->panels;
213
214

	}
215

216
	image->lambda = ph_en_to_lambda(eV_to_J(1.8e3));
217
	image->features = template->features;
218
	image->filename = template->filename;
219
	image->indexed_cell = template->indexed_cell;
Thomas White's avatar
Thomas White committed
220
	image->f0 = template->f0;
221

Thomas White's avatar
Thomas White committed
222
223
224
225
	/* Prevent muppetry */
	image->hits = NULL;
	image->n_hits = 0;

226
227
	return image;
}
228
229


230
static void simulate_and_write(struct image *simage, struct gpu_context **gctx,
Thomas White's avatar
Thomas White committed
231
                               const double *intensities, UnitCell *cell)
232
{
233
	/* Set up GPU if necessary */
Thomas White's avatar
Thomas White committed
234
	if ( (gctx != NULL) && (*gctx == NULL) ) {
Thomas White's avatar
Thomas White committed
235
		*gctx = setup_gpu(0, simage, intensities);
236
237
	}

Thomas White's avatar
Thomas White committed
238
	if ( (gctx != NULL) && (*gctx != NULL) ) {
239
		get_diffraction_gpu(*gctx, simage, 24, 24, 40, cell);
240
	} else {
241
		get_diffraction(simage, 24, 24, 40,
Thomas White's avatar
Thomas White committed
242
		                intensities, NULL, cell, 0,
243
		                GRADIENT_MOSAIC);
244
	}
245
	record_image(simage, 0);
246

247
	hdf5_write("simulated.h5", simage->data, simage->width, simage->height,
248
249
250
251
		   H5T_NATIVE_FLOAT);
}


252
static struct process_result process_image(struct process_args *pargs)
253
254
255
256
257
{
	struct hdfile *hdfile;
	struct image image;
	struct image *simage;
	float *data_for_measurement;
258
	struct process_result result;
259
	size_t data_size;
260
261
262
263
264
265
266
267
268
269
270
271
	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
272
	int config_polar = pargs->config_polar;
273
274
275
	IndexingMethod indm = pargs->indm;
	const double *intensities = pargs->intensities;
	struct gpu_context *gctx = pargs->gctx;
276
277
278
279

	image.features = NULL;
	image.data = NULL;
	image.indexed_cell = NULL;
Thomas White's avatar
Thomas White committed
280
	image.id = pargs->id;
281
	image.filename = filename;
Thomas White's avatar
Thomas White committed
282
283
	image.hits = NULL;
	image.n_hits = 0;
Thomas White's avatar
Thomas White committed
284
	image.det = pargs->det;
285

286
287
288
289
290
291
	/* 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;

292
	STATUS("Processing '%s'\n", image.filename);
293

294
295
	result.peaks_sane = 0;
	result.hit = 0;
296

297
298
	hdfile = hdfile_open(filename);
	if ( hdfile == NULL ) {
299
		return result;
300
301
	} else if ( hdfile_set_first_image(hdfile, "/") ) {
		ERROR("Couldn't select path\n");
302
		return result;
303
304
	}

305
	hdf5_read(hdfile, &image, pargs->config_satcorr);
306
307
308
309
310
311
312
313
314
315
316
317
318

	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 {
319
		memcpy(data_for_measurement, image.data, data_size);
320
321
322
323
	}

	/* Perform 'fine' peak search */
	search_peaks(&image);
Thomas White's avatar
Thomas White committed
324
325
326
327
328
329

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

330
331
	if ( image_feature_count(image.features) < 5 ) goto done;

332
	if ( config_dumpfound ) dump_peaks(&image, pargs->output_mutex);
333
334
335
336
337
338
339
340
341
342

	/* 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,
343
		              config_verbose, pargs->ipriv);
344
345
346
347
348
	}

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

Thomas White's avatar
Thomas White committed
349
	/* Sanity check */
350
351
	if ( pargs->config_sanity
	  && !peak_sanity_check(&image, image.indexed_cell) ) {
Thomas White's avatar
Thomas White committed
352
353
354
		STATUS("Failed peak sanity check.\n");
		goto done;
	} else {
355
		result.peaks_sane = 1;
Thomas White's avatar
Thomas White committed
356
357
	}

358
359
	/* Measure intensities if requested */
	if ( config_nearbragg ) {
360
		output_intensities(&image, image.indexed_cell,
Thomas White's avatar
Thomas White committed
361
		                   pargs->output_mutex, config_polar,
362
		                   pargs->config_sa);
363
364
	}

365
366
	simage = get_simage(&image, config_alternate);

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

	/* 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 */
	free(image.indexed_cell);

done:
	free(image.data);
Thomas White's avatar
Thomas White committed
390
	free(image.flags);
391
	image_feature_list_free(image.features);
Thomas White's avatar
Thomas White committed
392
	free(image.hits);
393
394
	hdfile_close(hdfile);

395
	if ( image.indexed_cell == NULL ) {
396
		result.hit = 0;
397
	} else {
398
		result.hit = 1;
399
400
	}
	return result;
401
402
403
}


404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
static void *worker_thread(void *pargsv)
{
	struct process_args *pargs = pargsv;
	int finish;

	do {

		struct process_result result;
		int wakeup;

		result = process_image(pargs);

		pthread_mutex_lock(&pargs->control_mutex);
		pargs->hit = result.hit;
		pargs->peaks_sane = result.peaks_sane;
		pargs->done = 1;
		pargs->start = 0;
		pthread_mutex_unlock(&pargs->control_mutex);

		/* Go to sleep until told to exit or process next image */
		do {

			pthread_mutex_lock(&pargs->control_mutex);
			/* Either of these can result in the thread waking up */
			wakeup = pargs->start || pargs->finish;
			finish = pargs->finish;
			pthread_mutex_unlock(&pargs->control_mutex);
			usleep(20000);

		} while ( !wakeup );

	} while ( !pargs->finish );

	return NULL;
}


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

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

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

		switch (c) {
Thomas White's avatar
Thomas White committed
518
		case 'h' :
Thomas White's avatar
Thomas White committed
519
520
521
			show_help(argv[0]);
			return 0;

Thomas White's avatar
Thomas White committed
522
		case 'i' :
Thomas White's avatar
Thomas White committed
523
524
525
			filename = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
526
		case 'z' :
Thomas White's avatar
Thomas White committed
527
528
529
			indm_str = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
530
		case 'q' :
531
532
533
			intfile = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
534
		case 'p' :
Thomas White's avatar
Thomas White committed
535
536
537
			pdb = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
538
		case 'x' :
539
540
541
			prefix = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
542
		case 'j' :
543
544
545
			nthreads = atoi(optarg);
			break;

Thomas White's avatar
Thomas White committed
546
		case 'g' :
Thomas White's avatar
Thomas White committed
547
548
549
			geometry = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
550
		case 0 :
Thomas White's avatar
Thomas White committed
551
552
			break;

Thomas White's avatar
Thomas White committed
553
		default :
Thomas White's avatar
Thomas White committed
554
555
556
557
558
559
560
561
562
563
564
565
566
567
			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
568
		ERROR("Failed to open input file '%s'\n", filename);
Thomas White's avatar
Thomas White committed
569
570
		return 1;
	}
Thomas White's avatar
Thomas White committed
571
	free(filename);
Thomas White's avatar
Thomas White committed
572

573
	if ( intfile != NULL ) {
Thomas White's avatar
Thomas White committed
574
575
576
		ReflItemList *items;
		items = read_reflections(intfile, intensities, NULL, NULL);
		delete_items(items);
577
578
579
580
	} else {
		intensities = NULL;
	}

Thomas White's avatar
Thomas White committed
581
582
583
584
	if ( pdb == NULL ) {
		pdb = strdup("molecule.pdb");
	}

585
	if ( prefix == NULL ) {
Thomas White's avatar
Thomas White committed
586
		prefix = strdup("");
587
588
	}

589
	if ( (nthreads == 0) || (nthreads > MAX_THREADS) ) {
590
591
592
593
		ERROR("Invalid number of threads.\n");
		return 1;
	}

Thomas White's avatar
Thomas White committed
594
595
	if ( indm_str == NULL ) {
		STATUS("You didn't specify an indexing method, so I won't"
596
597
598
		       " 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
599
600
601
602
603
		indm = INDEXING_NONE;
	} else if ( strcmp(indm_str, "none") == 0 ) {
		indm = INDEXING_NONE;
	} else if ( strcmp(indm_str, "dirax") == 0) {
		indm = INDEXING_DIRAX;
604
605
	} else if ( strcmp(indm_str, "template") == 0) {
		indm = INDEXING_TEMPLATE;
Thomas White's avatar
Thomas White committed
606
607
608
609
	} else {
		ERROR("Unrecognised indexing method '%s'\n", indm_str);
		return 1;
	}
610
	free(indm_str);
Thomas White's avatar
Thomas White committed
611

Thomas White's avatar
Thomas White committed
612
613
614
615
616
617
618
619
620
621
622
623
	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);

624
	if ( (!config_nomatch) || (indm == INDEXING_TEMPLATE) ) {
625
626
		cell = load_cell_from_pdb(pdb);
		if ( cell == NULL ) {
627
628
629
			ERROR("Couldn't read unit cell (from %s)\n", pdb);
			return 1;
		}
630
631
632
	} else {
		STATUS("No cell needed because --no-match was used.\n");
		cell = NULL;
633
	}
634
	free(pdb);
635

636
637
638
639
640
641
642
643
	/* 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);
644
	ipriv = prepare_indexing(indm, cell, prepare_filename, det);
645
646
647
648
649
	if ( ipriv == NULL ) {
		ERROR("Failed to prepare indexing.\n");
		return 1;
	}

Thomas White's avatar
Thomas White committed
650
	gsl_set_error_handler_off();
Thomas White's avatar
Thomas White committed
651
652
	n_images = 0;
	n_hits = 0;
Thomas White's avatar
Thomas White committed
653
	n_sane = 0;
654

655
656
657
658
659
660
	for ( i=0; i<nthreads; i++ ) {
		worker_args[i] = malloc(sizeof(struct process_args));
		worker_args[i]->filename = malloc(1024);
		worker_active[i] = 0;
	}

661
	/* Start threads off */
662
	for ( i=0; i<nthreads; i++ ) {
Thomas White's avatar
Thomas White committed
663
664

		char line[1024];
665
666
		struct process_args *pargs;
		int r;
Thomas White's avatar
Thomas White committed
667

668
		pargs = worker_args[i];
Thomas White's avatar
Thomas White committed
669

670
671
672
673
674
675
676
		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
677
		chomp(line);
678
		snprintf(pargs->filename, 1023, "%s%s", prefix, line);
679

680
681
		n_images++;

682
		pargs->output_mutex = &output_mutex;
683
		pargs->gpu_mutex = &gpu_mutex;
684
		pthread_mutex_init(&pargs->control_mutex, NULL);
685
686
687
688
689
690
691
692
693
694
		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
695
		pargs->config_polar = config_polar;
696
		pargs->config_sanity = config_sanity;
697
		pargs->config_satcorr = config_satcorr;
698
		pargs->config_sa = config_sa;
699
		pargs->cell = cell;
Thomas White's avatar
Thomas White committed
700
		pargs->det = det;
701
		pargs->ipriv = ipriv;
702
703
704
		pargs->indm = indm;
		pargs->intensities = intensities;
		pargs->gctx = gctx;
Thomas White's avatar
Thomas White committed
705
		pargs->id = i;
706
707
708
709
710
		pthread_mutex_lock(&pargs->control_mutex);
		pargs->done = 0;
		pargs->start = 1;
		pargs->finish = 0;
		pthread_mutex_unlock(&pargs->control_mutex);
711

Thomas White's avatar
Thomas White committed
712
		worker_active[i] = 1;
713
		r = pthread_create(&workers[i], NULL, worker_thread, pargs);
714
		if ( r != 0 ) {
Thomas White's avatar
Thomas White committed
715
			worker_active[i] = 0;
716
717
718
719
720
			ERROR("Couldn't start thread %i\n", i);
		}

	}

721
	/* Keep threads busy until the end of the data */
722
723
724
725
726
727
728
729
	do {

		int i;

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

			char line[1024];
			struct process_args *pargs;
730
731
732
733
			int done;

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

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

738
			/* Has the thread finished yet? */
739
			pargs = worker_args[i];
740
741
742
743
			pthread_mutex_lock(&pargs->control_mutex);
			done = pargs->done;
			pthread_mutex_unlock(&pargs->control_mutex);
			if ( !done ) continue;
744

745
746
			/* Results will be processed after checking if
			 * there are any more images to process. */
747

748
			/* Get next filename */
749
			rval = fgets(line, 1023, fh);
750
751
752
			/* In this case, the result of the last file
			 * file will be processed when the thread is
			 * joined. */
753
			if ( rval == NULL ) break;
754
755
756
757
758

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

759
			chomp(line);
760
			snprintf(pargs->filename, 1023, "%s%s", prefix, line);
761
762

			n_images++;
763
764
765
766
767
768
769

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

770
		}
Thomas White's avatar
Thomas White committed
771

Thomas White's avatar
Thomas White committed
772
773
	} while ( rval != NULL );

774
	/* Join threads */
775
	for ( i=0; i<nthreads; i++ ) {
Thomas White's avatar
Thomas White committed
776

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

779
780
		/* Tell the thread to exit */
		struct process_args *pargs = worker_args[i];
Thomas White's avatar
Thomas White committed
781
		pthread_mutex_lock(&pargs->control_mutex);
782
		pargs->finish = 1;
Thomas White's avatar
Thomas White committed
783
		pthread_mutex_unlock(&pargs->control_mutex);
Thomas White's avatar
Thomas White committed
784

785
786
		/* Wait for it to join */
		pthread_join(workers[i], NULL);
Thomas White's avatar
Thomas White committed
787
788
		worker_active[i] = 0;

789
790
		n_hits += pargs->hit;
		n_sane += pargs->peaks_sane;
Thomas White's avatar
Thomas White committed
791

Thomas White's avatar
Thomas White committed
792
	free:
793
		if ( worker_args[i]->filename != NULL ) {
794
795
			free(worker_args[i]->filename);
		}
796
		free(worker_args[i]);
Thomas White's avatar
Thomas White committed
797

798
799
800
	}

	free(prefix);
Thomas White's avatar
Thomas White committed
801
802
	free(det->panels);
	free(det);
803
	free(cell);
Thomas White's avatar
Thomas White committed
804
805
806
	fclose(fh);

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

Thomas White's avatar
Thomas White committed
809
810
811
812
	if ( gctx != NULL ) {
		cleanup_gpu(gctx);
	}

Thomas White's avatar
Thomas White committed
813
814
	return 0;
}