indexamajig.c 19.4 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
 *
 * Part of CrystFEL - crystallography with a FEL
 *
 */


#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

17
#define _GNU_SOURCE 1
Thomas White's avatar
Thomas White committed
18
19
20
21
22
23
#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
24
#include <hdf5.h>
Thomas White's avatar
Thomas White committed
25
#include <gsl/gsl_errno.h>
26
#include <pthread.h>
Thomas White's avatar
Thomas White committed
27
#include <sys/time.h>
Thomas White's avatar
Thomas White committed
28
29
30

#include "utils.h"
#include "hdf5-file.h"
Thomas White's avatar
Thomas White committed
31
#include "index.h"
32
#include "peaks.h"
33
#include "diffraction.h"
Thomas White's avatar
Thomas White committed
34
#include "diffraction-gpu.h"
35
#include "detector.h"
36
#include "sfac.h"
Thomas White's avatar
Thomas White committed
37
#include "filters.h"
38
#include "reflections.h"
Thomas White's avatar
Thomas White committed
39
40


41
42
43
44
45
#define MAX_THREADS (96)

struct process_args
{
	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
66
67
68
69
70
71
72
73
	IndexingMethod indm;
	const double *intensities;
	const unsigned int *counts;
	struct gpu_context *gctx;
};

struct process_result
{
	int hit;
Thomas White's avatar
Thomas White committed
74
	int peaks_sane;
75
76
77
};


Thomas White's avatar
Thomas White committed
78
79
80
81
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"
"  -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
88
"\n"
Thomas White's avatar
Thomas White committed
89
90
91
"      --indexing=<method> Use 'method' for indexing.  Choose from:\n"
"                           none     : no indexing\n"
"                           dirax    : invoke DirAx\n"
Thomas White's avatar
Thomas White committed
92
"  -g. --geometry=<file>   Get detector geometry from file.\n"
Thomas White's avatar
Thomas White committed
93
94
"\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"
95
"      --near-bragg        Output a list of reflection intensities to stdout.\n"
Thomas White's avatar
Thomas White committed
96
97
98
99
100
101
102
103
104
105
"                           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"
106
"      --simulate          Simulate the diffraction pattern using the indexed\n"
Thomas White's avatar
Thomas White committed
107
108
109
110
"                           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"
111
112
113
114
115
116
117
"      --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"
118
119
120
121
122
"      --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"
123
124
"                           The intensities in this list are from the\n"
"                           centroid/integration procedure.\n"
125
126
127
"      --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"
128
"      --unpolarized       Don't correct for the polarisation of the X-rays.\n"
129
130
"      --check-sanity      Check that indexed locations approximately correspond\n"
"                           with detected peaks.\n"
131
132
"      --sat-corr          Correct values of saturated peaks using a table\n"
"                           included in the HDF5 file.\n"
133
134
"      --no-sa             Don't correct for the differing solid angles of\n"
"                           the pixels.\n"
Thomas White's avatar
Thomas White committed
135
136
137
138
139
"\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
140
141
142
"     --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"
143
" -x, --prefix=<p>         Prefix filenames from input file with 'p'.\n"
Thomas White's avatar
Thomas White committed
144
);
Thomas White's avatar
Thomas White committed
145
146
147
}


148
static struct image *get_simage(struct image *template, int alternate)
149
{
150
	struct image *image;
151
	struct panel panels[2];
152

153
154
	image = malloc(sizeof(*image));

155
	/* Simulate a diffraction pattern */
156
157
158
	image->twotheta = NULL;
	image->data = NULL;
	image->det = template->det;
159
160

	/* View head-on (unit cell is tilted) */
161
162
163
164
	image->orientation.w = 1.0;
	image->orientation.x = 0.0;
	image->orientation.y = 0.0;
	image->orientation.z = 0.0;
165
166
167

	/* Detector geometry for the simulation
	 * - not necessarily the same as the original. */
168
169
	image->width = 1024;
	image->height = 1024;
Thomas White's avatar
Thomas White committed
170
	image->det->n_panels = 2;
171

172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
	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;
190
		panels[1].cy = 525.0;
191
192
193
		panels[1].clen = 56.7e-2;  /* 56.7 cm */
		panels[1].res = 13333.3;   /* 75 microns/pixel */

Thomas White's avatar
Thomas White committed
194
		image->det->panels = panels;
195

196
197
198
	} else {

		/* Copy pointer to old geometry */
Thomas White's avatar
Thomas White committed
199
		image->det->panels = template->det->panels;
200
201

	}
202

203
	image->lambda = ph_en_to_lambda(eV_to_J(1.8e3));
204
	image->features = template->features;
205
	image->filename = template->filename;
206
	image->indexed_cell = template->indexed_cell;
Thomas White's avatar
Thomas White committed
207
	image->f0 = template->f0;
208

Thomas White's avatar
Thomas White committed
209
210
211
212
	/* Prevent muppetry */
	image->hits = NULL;
	image->n_hits = 0;

213
214
	return image;
}
215
216


217
static void simulate_and_write(struct image *simage, struct gpu_context **gctx,
218
219
                               const double *intensities,
                               const unsigned int *counts, UnitCell *cell)
220
{
221
	/* Set up GPU if necessary */
Thomas White's avatar
Thomas White committed
222
	if ( (gctx != NULL) && (*gctx == NULL) ) {
223
		*gctx = setup_gpu(0, simage, intensities, counts);
224
225
	}

Thomas White's avatar
Thomas White committed
226
	if ( (gctx != NULL) && (*gctx != NULL) ) {
227
		get_diffraction_gpu(*gctx, simage, 24, 24, 40, cell);
228
	} else {
229
		get_diffraction(simage, 24, 24, 40,
230
231
		                intensities, counts, NULL, cell, 0,
		                GRADIENT_MOSAIC);
232
	}
233
	record_image(simage, 0);
234

235
	hdf5_write("simulated.h5", simage->data, simage->width, simage->height,
236
237
238
239
		   H5T_NATIVE_FLOAT);
}


240
static void *process_image(void *pargsv)
241
{
242
	struct process_args *pargs = pargsv;
243
244
245
246
247
	struct hdfile *hdfile;
	struct image image;
	struct image *simage;
	float *data_for_measurement;
	size_t data_size;
248
249
250
251
252
253
254
255
256
257
258
259
	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
260
	int config_polar = pargs->config_polar;
261
262
263
264
265
	IndexingMethod indm = pargs->indm;
	const double *intensities = pargs->intensities;
	const unsigned int *counts = pargs->counts;
	struct gpu_context *gctx = pargs->gctx;
	struct process_result *result;
266
267
268
269

	image.features = NULL;
	image.data = NULL;
	image.indexed_cell = NULL;
Thomas White's avatar
Thomas White committed
270
	image.id = pargs->id;
271
	image.filename = filename;
Thomas White's avatar
Thomas White committed
272
273
	image.hits = NULL;
	image.n_hits = 0;
Thomas White's avatar
Thomas White committed
274
	image.det = pargs->det;
275

276
277
278
279
280
281
	/* 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;

282
	STATUS("Processing '%s'\n", image.filename);
283

284
285
	result = malloc(sizeof(*result));
	if ( result == NULL ) return NULL;
286
287
	result->peaks_sane = 0;
	result->hit = 0;
288

289
290
	hdfile = hdfile_open(filename);
	if ( hdfile == NULL ) {
291
		return result;
292
293
	} else if ( hdfile_set_first_image(hdfile, "/") ) {
		ERROR("Couldn't select path\n");
294
		return result;
295
296
	}

297
	hdf5_read(hdfile, &image, pargs->config_satcorr);
298
299
300
301
302
303
304
305
306
307
308
309
310

	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 {
311
		memcpy(data_for_measurement, image.data, data_size);
312
313
314
315
	}

	/* Perform 'fine' peak search */
	search_peaks(&image);
Thomas White's avatar
Thomas White committed
316
317
318
319
320
321

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

322
323
	if ( image_feature_count(image.features) < 5 ) goto done;

324
	if ( config_dumpfound ) dump_peaks(&image, pargs->output_mutex);
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340

	/* 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,
		              config_verbose);
	}

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

Thomas White's avatar
Thomas White committed
341
	/* Sanity check */
342
343
	if ( pargs->config_sanity
	  && !peak_sanity_check(&image, image.indexed_cell) ) {
Thomas White's avatar
Thomas White committed
344
345
346
347
348
349
		STATUS("Failed peak sanity check.\n");
		goto done;
	} else {
		result->peaks_sane = 1;
	}

350
351
	/* Measure intensities if requested */
	if ( config_nearbragg ) {
352
		output_intensities(&image, image.indexed_cell,
Thomas White's avatar
Thomas White committed
353
		                   pargs->output_mutex, config_polar,
354
		                   pargs->config_sa);
355
356
	}

357
358
	simage = get_simage(&image, config_alternate);

359
360
361
	/* Simulate if requested */
	if ( config_simulate ) {
		if ( config_gpu ) {
362
			pthread_mutex_lock(pargs->gpu_mutex);
363
			simulate_and_write(simage, &gctx, intensities,
364
			                   counts, image.indexed_cell);
365
			pthread_mutex_unlock(pargs->gpu_mutex);
366
367
		} else {
			simulate_and_write(simage, NULL, intensities,
368
			                   counts, image.indexed_cell);
369
370
371
372
373
374
375
376
377
378
379
380
381
		}
	}

	/* 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
382
	free(image.flags);
383
	image_feature_list_free(image.features);
Thomas White's avatar
Thomas White committed
384
	free(image.hits);
385
386
	hdfile_close(hdfile);

387
388
389
390
391
392
	if ( image.indexed_cell == NULL ) {
		result->hit = 0;
	} else {
		result->hit = 1;
	}
	return result;
393
394
395
}


Thomas White's avatar
Thomas White committed
396
397
398
int main(int argc, char *argv[])
{
	int c;
Thomas White's avatar
Thomas White committed
399
	struct gpu_context *gctx = NULL;
Thomas White's avatar
Thomas White committed
400
401
	char *filename = NULL;
	FILE *fh;
402
	char *rval = NULL;
Thomas White's avatar
Thomas White committed
403
404
	int n_images;
	int n_hits;
Thomas White's avatar
Thomas White committed
405
	int n_sane;
406
	int config_noindex = 0;
Thomas White's avatar
Thomas White committed
407
	int config_dumpfound = 0;
408
	int config_nearbragg = 0;
Thomas White's avatar
Thomas White committed
409
	int config_writedrx = 0;
410
	int config_simulate = 0;
411
412
	int config_cmfilter = 0;
	int config_noisefilter = 0;
413
	int config_nomatch = 0;
Thomas White's avatar
Thomas White committed
414
	int config_gpu = 0;
415
	int config_verbose = 0;
416
	int config_alternate = 0;
Thomas White's avatar
Thomas White committed
417
	int config_polar = 1;
418
	int config_sanity = 0;
419
	int config_satcorr = 0;
420
	int config_sa = 1;
Thomas White's avatar
Thomas White committed
421
422
	struct detector *det;
	char *geometry = NULL;
Thomas White's avatar
Thomas White committed
423
424
	IndexingMethod indm;
	char *indm_str = NULL;
425
426
427
	UnitCell *cell;
	double *intensities = NULL;
	char *intfile = NULL;
428
	unsigned int *counts;
Thomas White's avatar
Thomas White committed
429
	char *pdb = NULL;
430
	char *prefix = NULL;
431
432
433
	int nthreads = 1;
	pthread_t workers[MAX_THREADS];
	struct process_args *worker_args[MAX_THREADS];
Thomas White's avatar
Thomas White committed
434
	int worker_active[MAX_THREADS];
435
	int i;
436
	pthread_mutex_t output_mutex = PTHREAD_MUTEX_INITIALIZER;
437
	pthread_mutex_t gpu_mutex = PTHREAD_MUTEX_INITIALIZER;
Thomas White's avatar
Thomas White committed
438
439
440
441
442

	/* Long options */
	const struct option longopts[] = {
		{"help",               0, NULL,               'h'},
		{"input",              1, NULL,               'i'},
Thomas White's avatar
Thomas White committed
443
		{"gpu",                0, &config_gpu,         1},
444
		{"no-index",           0, &config_noindex,     1},
Thomas White's avatar
Thomas White committed
445
		{"dump-peaks",         0, &config_dumpfound,   1},
446
		{"near-bragg",         0, &config_nearbragg,   1},
Thomas White's avatar
Thomas White committed
447
448
		{"write-drx",          0, &config_writedrx,    1},
		{"indexing",           1, NULL,               'z'},
Thomas White's avatar
Thomas White committed
449
		{"geometry",           1, NULL,               'g'},
450
		{"simulate",           0, &config_simulate,    1},
451
452
		{"filter-cm",          0, &config_cmfilter,    1},
		{"filter-noise",       0, &config_noisefilter, 1},
453
		{"no-match",           0, &config_nomatch,     1},
454
		{"verbose",            0, &config_verbose,     1},
455
		{"alternate",          0, &config_alternate,   1},
456
		{"intensities",        1, NULL,               'q'},
Thomas White's avatar
Thomas White committed
457
		{"pdb",                1, NULL,               'p'},
458
		{"prefix",             1, NULL,               'x'},
Thomas White's avatar
Thomas White committed
459
		{"unpolarized",        0, &config_polar,       0},
460
		{"check-sanity",       0, &config_sanity,      1},
461
		{"sat-corr",           0, &config_satcorr,     1},
462
		{"no-sa",              0, &config_sa,          0},
Thomas White's avatar
Thomas White committed
463
464
465
466
		{0, 0, NULL, 0}
	};

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

		switch (c) {
Thomas White's avatar
Thomas White committed
471
		case 'h' :
Thomas White's avatar
Thomas White committed
472
473
474
			show_help(argv[0]);
			return 0;

Thomas White's avatar
Thomas White committed
475
		case 'i' :
Thomas White's avatar
Thomas White committed
476
477
478
			filename = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
479
		case 'z' :
Thomas White's avatar
Thomas White committed
480
481
482
			indm_str = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
483
		case 'q' :
484
485
486
			intfile = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
487
		case 'p' :
Thomas White's avatar
Thomas White committed
488
489
490
			pdb = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
491
		case 'x' :
492
493
494
			prefix = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
495
		case 'j' :
496
497
498
			nthreads = atoi(optarg);
			break;

Thomas White's avatar
Thomas White committed
499
		case 'g' :
Thomas White's avatar
Thomas White committed
500
501
502
			geometry = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
503
		case 0 :
Thomas White's avatar
Thomas White committed
504
505
			break;

Thomas White's avatar
Thomas White committed
506
		default :
Thomas White's avatar
Thomas White committed
507
508
509
510
511
512
513
514
515
516
517
518
519
520
			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
521
		ERROR("Failed to open input file '%s'\n", filename);
Thomas White's avatar
Thomas White committed
522
523
		return 1;
	}
Thomas White's avatar
Thomas White committed
524
	free(filename);
Thomas White's avatar
Thomas White committed
525

526
	if ( intfile != NULL ) {
527
		counts = new_list_count();
528
		intensities = read_reflections(intfile, counts, NULL);
529
530
	} else {
		intensities = NULL;
531
		counts = NULL;
532
533
	}

Thomas White's avatar
Thomas White committed
534
535
536
537
	if ( pdb == NULL ) {
		pdb = strdup("molecule.pdb");
	}

538
	if ( prefix == NULL ) {
Thomas White's avatar
Thomas White committed
539
		prefix = strdup("");
540
541
	}

542
	if ( (nthreads == 0) || (nthreads > MAX_THREADS) ) {
543
544
545
546
		ERROR("Invalid number of threads.\n");
		return 1;
	}

Thomas White's avatar
Thomas White committed
547
548
	if ( indm_str == NULL ) {
		STATUS("You didn't specify an indexing method, so I won't"
549
550
551
		       " 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
552
553
554
555
556
557
558
559
560
		indm = INDEXING_NONE;
	} else if ( strcmp(indm_str, "none") == 0 ) {
		indm = INDEXING_NONE;
	} else if ( strcmp(indm_str, "dirax") == 0) {
		indm = INDEXING_DIRAX;
	} else {
		ERROR("Unrecognised indexing method '%s'\n", indm_str);
		return 1;
	}
561
	free(indm_str);
Thomas White's avatar
Thomas White committed
562

Thomas White's avatar
Thomas White committed
563
564
565
566
567
568
569
570
571
572
573
574
	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);

Thomas White's avatar
Thomas White committed
575
	cell = load_cell_from_pdb(pdb);
576
	if ( cell == NULL ) {
Thomas White's avatar
Thomas White committed
577
		ERROR("Couldn't read unit cell (from %s)\n", pdb);
578
579
		return 1;
	}
580
	free(pdb);
581

Thomas White's avatar
Thomas White committed
582
	gsl_set_error_handler_off();
Thomas White's avatar
Thomas White committed
583
584
	n_images = 0;
	n_hits = 0;
Thomas White's avatar
Thomas White committed
585
	n_sane = 0;
586

587
588
589
590
591
592
	for ( i=0; i<nthreads; i++ ) {
		worker_args[i] = malloc(sizeof(struct process_args));
		worker_args[i]->filename = malloc(1024);
		worker_active[i] = 0;
	}

593
594
	/* Initially, fire off the full number of threads */
	for ( i=0; i<nthreads; i++ ) {
Thomas White's avatar
Thomas White committed
595
596

		char line[1024];
597
598
		struct process_args *pargs;
		int r;
Thomas White's avatar
Thomas White committed
599

600
		pargs = worker_args[i];
Thomas White's avatar
Thomas White committed
601

Thomas White's avatar
Thomas White committed
602
		rval = fgets(line, 1023, fh);
Thomas White's avatar
Thomas White committed
603
		if ( rval == NULL ) continue;
Thomas White's avatar
Thomas White committed
604
		chomp(line);
605
		snprintf(pargs->filename, 1023, "%s%s", prefix, line);
606

607
608
		n_images++;

609
		pargs->output_mutex = &output_mutex;
610
		pargs->gpu_mutex = &gpu_mutex;
611
612
613
614
615
616
617
618
619
620
		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
621
		pargs->config_polar = config_polar;
622
		pargs->config_sanity = config_sanity;
623
		pargs->config_satcorr = config_satcorr;
624
		pargs->config_sa = config_sa;
625
		pargs->cell = cell;
Thomas White's avatar
Thomas White committed
626
		pargs->det = det;
627
628
629
630
		pargs->indm = indm;
		pargs->intensities = intensities;
		pargs->counts = counts;
		pargs->gctx = gctx;
Thomas White's avatar
Thomas White committed
631
		pargs->id = i;
632

Thomas White's avatar
Thomas White committed
633
		worker_active[i] = 1;
634
635
		r = pthread_create(&workers[i], NULL, process_image, pargs);
		if ( r != 0 ) {
Thomas White's avatar
Thomas White committed
636
			worker_active[i] = 0;
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
			ERROR("Couldn't start thread %i\n", i);
		}

	}

	/* Start new threads as old ones finish */
	do {

		int i;

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

			char line[1024];
			int r;
			struct process_result *result = NULL;
			struct timespec t;
Thomas White's avatar
Thomas White committed
653
			struct timeval tv;
654
655
			struct process_args *pargs;

Thomas White's avatar
Thomas White committed
656
657
			if ( !worker_active[i] ) continue;

658
659
			pargs = worker_args[i];

Thomas White's avatar
Thomas White committed
660
661
662
			gettimeofday(&tv, NULL);
			t.tv_sec = tv.tv_sec;
			t.tv_nsec = tv.tv_usec * 1000 + 20000;
663
664
665
666
667

			r = pthread_timedjoin_np(workers[i], (void *)&result,
			                         &t);
			if ( r != 0 ) continue; /* Not ready yet */

Thomas White's avatar
Thomas White committed
668
669
			worker_active[i] = 0;

670
671
			if ( result != NULL ) {
				n_hits += result->hit;
Thomas White's avatar
Thomas White committed
672
				n_sane += result->peaks_sane;
673
674
675
676
677
678
				free(result);
			}

			rval = fgets(line, 1023, fh);
			if ( rval == NULL ) break;
			chomp(line);
679
			snprintf(pargs->filename, 1023, "%s%s", prefix, line);
680

Thomas White's avatar
Thomas White committed
681
			worker_active[i] = 1;
682
683
684
			r = pthread_create(&workers[i], NULL, process_image,
			                   pargs);
			if ( r != 0 ) {
Thomas White's avatar
Thomas White committed
685
				worker_active[i] = 0;
686
687
688
689
690
				ERROR("Couldn't start thread %i\n", i);
			}

			n_images++;
		}
Thomas White's avatar
Thomas White committed
691

Thomas White's avatar
Thomas White committed
692
693
	} while ( rval != NULL );

Thomas White's avatar
Thomas White committed
694
	/* Catch all remaining threads */
695
	for ( i=0; i<nthreads; i++ ) {
Thomas White's avatar
Thomas White committed
696
697
698

		struct process_result *result = NULL;

Thomas White's avatar
Thomas White committed
699
		if ( !worker_active[i] ) goto free;
Thomas White's avatar
Thomas White committed
700
701
702
703
704
705
706
707
708
709

		pthread_join(workers[i], (void *)&result);

		worker_active[i] = 0;

		if ( result != NULL ) {
			n_hits += result->hit;
			free(result);
		}

Thomas White's avatar
Thomas White committed
710
	free:
711
		if ( worker_args[i]->filename != NULL ) {
712
713
			free(worker_args[i]->filename);
		}
714
		free(worker_args[i]);
Thomas White's avatar
Thomas White committed
715

716
717
718
	}

	free(prefix);
Thomas White's avatar
Thomas White committed
719
720
	free(det->panels);
	free(det);
721
	free(cell);
Thomas White's avatar
Thomas White committed
722
723
724
	fclose(fh);

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

Thomas White's avatar
Thomas White committed
727
728
729
730
	if ( gctx != NULL ) {
		cleanup_gpu(gctx);
	}

Thomas White's avatar
Thomas White committed
731
732
	return 0;
}