indexamajig.c 15.8 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
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
#define MAX_THREADS (96)

struct process_args
{
	char *filename;
	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;
	IndexingMethod indm;
	const double *intensities;
	const unsigned int *counts;
	struct gpu_context *gctx;
};

struct process_result
{
	int hit;
	struct process_args *pargs;
};


Thomas White's avatar
Thomas White committed
70
71
72
73
74
75
76
77
78
79
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
80
81
82
"      --indexing=<method> Use 'method' for indexing.  Choose from:\n"
"                           none     : no indexing\n"
"                           dirax    : invoke DirAx\n"
83
84
"\n"
"      --verbose           Be verbose about indexing.\n"
85
"      --gpu               Use the GPU to speed up the simulation.\n"
86
"  -j <n>                  Run <n> analyses in parallel.\n"
87
"\n"
88
"      --near-bragg        Output a list of reflection intensities to stdout.\n"
89
90
91
92
93
"                           The intensities in this list are the sum of\n"
"                           the values in a 7x7 square centered on the pixel\n"
"                           closest to the Bragg condition.  Only pixels with\n"
"                           fractional indices within 0.1 of the Bragg\n"
"                           condition will be counted.\n"
94
95
"      --simulate          Simulate the diffraction pattern using the indexed\n"
"                           unit cell.\n"
96
97
98
99
100
101
102
"      --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"
103
104
105
106
107
108
"\n"
"      --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"
109
110
"                           The intensities in this list are from the\n"
"                           centroid/integration procedure.\n"
111
112
113
"      --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"
Thomas White's avatar
Thomas White committed
114
115
116
"     --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"
117
" -x, --prefix=<p>         Prefix filenames from input file with 'p'.\n"
Thomas White's avatar
Thomas White committed
118
);
Thomas White's avatar
Thomas White committed
119
120
121
}


122
static struct image *get_simage(struct image *template, int alternate)
123
{
124
	struct image *image;
125
	struct panel panels[2];
126

127
128
	image = malloc(sizeof(*image));

129
	/* Simulate a diffraction pattern */
130
131
132
	image->twotheta = NULL;
	image->data = NULL;
	image->det = template->det;
133
134

	/* View head-on (unit cell is tilted) */
135
136
137
138
	image->orientation.w = 1.0;
	image->orientation.x = 0.0;
	image->orientation.y = 0.0;
	image->orientation.z = 0.0;
139
140
141

	/* Detector geometry for the simulation
	 * - not necessarily the same as the original. */
142
143
144
	image->width = 1024;
	image->height = 1024;
	image->det.n_panels = 2;
145

146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
	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;
		panels[1].cy = 772.1;
		panels[1].clen = 56.7e-2;  /* 56.7 cm */
		panels[1].res = 13333.3;   /* 75 microns/pixel */

		image->det.panels = panels;
169

170
171
172
173
174
175
	} else {

		/* Copy pointer to old geometry */
		image->det.panels = template->det.panels;

	}
176

177
	image->lambda = ph_en_to_lambda(eV_to_J(1.8e3));
178
	image->features = template->features;
179
180
181

	return image;
}
182
183


184
static void simulate_and_write(struct image *simage, struct gpu_context **gctx,
185
186
                               const double *intensities,
                               const unsigned int *counts, UnitCell *cell)
187
{
188
	/* Set up GPU if necessary */
Thomas White's avatar
Thomas White committed
189
	if ( (gctx != NULL) && (*gctx == NULL) ) {
190
		*gctx = setup_gpu(0, simage, intensities, counts);
191
192
	}

Thomas White's avatar
Thomas White committed
193
	if ( (gctx != NULL) && (*gctx != NULL) ) {
194
		get_diffraction_gpu(*gctx, simage, 24, 24, 40, cell);
195
	} else {
196
197
		get_diffraction(simage, 24, 24, 40,
		                intensities, counts, cell, 0);
198
	}
199
	record_image(simage, 0);
200

201
	hdf5_write("simulated.h5", simage->data, simage->width, simage->height,
202
203
204
205
		   H5T_NATIVE_FLOAT);
}


206
static void *process_image(void *pargsv)
207
{
208
	struct process_args *pargs = pargsv;
209
210
211
212
213
	struct hdfile *hdfile;
	struct image image;
	struct image *simage;
	float *data_for_measurement;
	size_t data_size;
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
	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;
	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;
231
232
233
234
235
236
237

	image.features = NULL;
	image.data = NULL;
	image.indexed_cell = NULL;

	STATUS("Processing '%s'\n", filename);

238
239
240
241
	result = malloc(sizeof(*result));
	if ( result == NULL ) return NULL;
	result->pargs = pargs;

242
243
	hdfile = hdfile_open(filename);
	if ( hdfile == NULL ) {
244
245
		result->hit = 0;
		return result;
246
247
	} else if ( hdfile_set_first_image(hdfile, "/") ) {
		ERROR("Couldn't select path\n");
248
249
		result->hit = 0;
		return result;
250
251
	}

252
253
	#include "geometry-lcls.tmp"

254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
	hdf5_read(hdfile, &image);

	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 {

		int x, y;

		for ( x=0; x<image.width; x++ ) {
		for ( y=0; y<image.height; y++ ) {
			float val;
			val = image.data[x+image.width*y];
			data_for_measurement[x+image.width*y] = val;
		}
		}

	}

	/* Perform 'fine' peak search */
	search_peaks(&image);
	if ( image_feature_count(image.features) < 5 ) goto done;

	if ( config_dumpfound ) dump_peaks(&image);

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

	simage = get_simage(&image, config_alternate);

	/* Measure intensities if requested */
	if ( config_nearbragg ) {
		/* Use original data (temporarily) */
		simage->data = data_for_measurement;
		output_intensities(simage, image.indexed_cell);
		simage->data = NULL;
	}

	/* Simulate if requested */
	if ( config_simulate ) {
		if ( config_gpu ) {
			simulate_and_write(simage, &gctx, intensities,
			                   counts, cell);
		} else {
			simulate_and_write(simage, NULL, intensities,
			                   counts, cell);
		}
	}

	/* 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);
	free(image.det.panels);
	image_feature_list_free(image.features);
	free(data_for_measurement);
	hdfile_close(hdfile);

338
339
340
341
342
343
	if ( image.indexed_cell == NULL ) {
		result->hit = 0;
	} else {
		result->hit = 1;
	}
	return result;
344
345
346
}


Thomas White's avatar
Thomas White committed
347
348
349
int main(int argc, char *argv[])
{
	int c;
Thomas White's avatar
Thomas White committed
350
	struct gpu_context *gctx = NULL;
Thomas White's avatar
Thomas White committed
351
352
	char *filename = NULL;
	FILE *fh;
353
	char *rval = NULL;
Thomas White's avatar
Thomas White committed
354
355
	int n_images;
	int n_hits;
356
	int config_noindex = 0;
Thomas White's avatar
Thomas White committed
357
	int config_dumpfound = 0;
358
	int config_nearbragg = 0;
Thomas White's avatar
Thomas White committed
359
	int config_writedrx = 0;
360
	int config_simulate = 0;
361
362
	int config_cmfilter = 0;
	int config_noisefilter = 0;
363
	int config_nomatch = 0;
Thomas White's avatar
Thomas White committed
364
	int config_gpu = 0;
365
	int config_verbose = 0;
366
	int config_alternate = 0;
Thomas White's avatar
Thomas White committed
367
368
	IndexingMethod indm;
	char *indm_str = NULL;
369
370
371
	UnitCell *cell;
	double *intensities = NULL;
	char *intfile = NULL;
372
	unsigned int *counts = NULL;
Thomas White's avatar
Thomas White committed
373
	char *pdb = NULL;
374
	char *prefix = NULL;
375
376
377
	int nthreads = 1;
	pthread_t workers[MAX_THREADS];
	struct process_args *worker_args[MAX_THREADS];
Thomas White's avatar
Thomas White committed
378
	int worker_active[MAX_THREADS];
379
	int i;
Thomas White's avatar
Thomas White committed
380
381
382
383
384

	/* Long options */
	const struct option longopts[] = {
		{"help",               0, NULL,               'h'},
		{"input",              1, NULL,               'i'},
Thomas White's avatar
Thomas White committed
385
		{"gpu",                0, &config_gpu,         1},
386
		{"no-index",           0, &config_noindex,     1},
Thomas White's avatar
Thomas White committed
387
		{"dump-peaks",         0, &config_dumpfound,   1},
388
		{"near-bragg",         0, &config_nearbragg,   1},
Thomas White's avatar
Thomas White committed
389
390
		{"write-drx",          0, &config_writedrx,    1},
		{"indexing",           1, NULL,               'z'},
391
		{"simulate",           0, &config_simulate,    1},
392
393
		{"filter-cm",          0, &config_cmfilter,    1},
		{"filter-noise",       0, &config_noisefilter, 1},
394
		{"no-match",           0, &config_nomatch,     1},
395
		{"verbose",            0, &config_verbose,     1},
396
		{"alternate",          0, &config_alternate,   1},
397
		{"intensities",        1, NULL,               'q'},
Thomas White's avatar
Thomas White committed
398
		{"pdb",                1, NULL,               'p'},
399
		{"prefix",             1, NULL,               'x'},
Thomas White's avatar
Thomas White committed
400
401
402
403
		{0, 0, NULL, 0}
	};

	/* Short options */
404
	while ((c = getopt_long(argc, argv, "hi:wp:j:", longopts, NULL)) != -1) {
Thomas White's avatar
Thomas White committed
405
406
407
408
409
410
411
412
413
414
415
416

		switch (c) {
		case 'h' : {
			show_help(argv[0]);
			return 0;
		}

		case 'i' : {
			filename = strdup(optarg);
			break;
		}

Thomas White's avatar
Thomas White committed
417
418
419
420
421
		case 'z' : {
			indm_str = strdup(optarg);
			break;
		}

422
423
424
425
426
		case 'q' : {
			intfile = strdup(optarg);
			break;
		}

Thomas White's avatar
Thomas White committed
427
428
429
430
431
		case 'p' : {
			pdb = strdup(optarg);
			break;
		}

432
433
434
435
436
		case 'x' : {
			prefix = strdup(optarg);
			break;
		}

437
438
439
440
441
		case 'j' : {
			nthreads = atoi(optarg);
			break;
		}

Thomas White's avatar
Thomas White committed
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
		case 0 : {
			break;
		}

		default : {
			return 1;
		}
		}

	}

	if ( filename == NULL ) {
		filename = strdup("-");
	}
	if ( strcmp(filename, "-") == 0 ) {
		fh = stdin;
	} else {
		fh = fopen(filename, "r");
	}
	free(filename);
	if ( fh == NULL ) {
		ERROR("Failed to open input file\n");
		return 1;
	}

467
	if ( intfile != NULL ) {
468
		intensities = read_reflections(intfile, counts);
469
470
	} else {
		intensities = NULL;
471
		counts = NULL;
472
473
	}

Thomas White's avatar
Thomas White committed
474
475
476
477
	if ( pdb == NULL ) {
		pdb = strdup("molecule.pdb");
	}

478
	if ( prefix == NULL ) {
Thomas White's avatar
Thomas White committed
479
		prefix = strdup("");
480
481
	}

482
483
484
485
486
	if ( nthreads == 0 ) {
		ERROR("Invalid number of threads.\n");
		return 1;
	}

Thomas White's avatar
Thomas White committed
487
488
	if ( indm_str == NULL ) {
		STATUS("You didn't specify an indexing method, so I won't"
489
490
491
		       " 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
492
493
494
495
496
497
498
499
500
		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;
	}
501
	free(indm_str);
Thomas White's avatar
Thomas White committed
502

Thomas White's avatar
Thomas White committed
503
504
	cell = load_cell_from_pdb(pdb);
	free(pdb);
505
506
	if ( cell == NULL ) {
		ERROR("Couldn't read unit cell (from molecule.pdb)\n");
507
508
509
		return 1;
	}

Thomas White's avatar
Thomas White committed
510
	gsl_set_error_handler_off();
Thomas White's avatar
Thomas White committed
511
512
	n_images = 0;
	n_hits = 0;
513
514
515

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

		char line[1024];
518
519
520
		char *prefixed;
		struct process_args *pargs;
		int r;
Thomas White's avatar
Thomas White committed
521

Thomas White's avatar
Thomas White committed
522
523
		worker_active[i] = 0;

Thomas White's avatar
Thomas White committed
524
		rval = fgets(line, 1023, fh);
Thomas White's avatar
Thomas White committed
525
		if ( rval == NULL ) continue;
Thomas White's avatar
Thomas White committed
526
		chomp(line);
527
		prefixed = malloc(1024);
528
529
		snprintf(prefixed, 1023, "%s%s", prefix, line);

530
531
		n_images++;

532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
		pargs = malloc(sizeof(*pargs));
		pargs->filename = prefixed;
		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;
		pargs->cell = cell;
		pargs->indm = indm;
		pargs->intensities = intensities;
		pargs->counts = counts;
		pargs->gctx = gctx;
		worker_args[i] = pargs;

Thomas White's avatar
Thomas White committed
551
		worker_active[i] = 1;
552
553
		r = pthread_create(&workers[i], NULL, process_image, pargs);
		if ( r != 0 ) {
Thomas White's avatar
Thomas White committed
554
			worker_active[i] = 0;
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
			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];
			char *prefixed;
			int r;
			struct process_result *result = NULL;
			struct timespec t;
Thomas White's avatar
Thomas White committed
572
			struct timeval tv;
573
574
			struct process_args *pargs;

Thomas White's avatar
Thomas White committed
575
576
577
578
579
			if ( !worker_active[i] ) continue;

			gettimeofday(&tv, NULL);
			t.tv_sec = tv.tv_sec;
			t.tv_nsec = tv.tv_usec * 1000 + 20000;
580
581
582
583
584

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

Thomas White's avatar
Thomas White committed
585
586
			worker_active[i] = 0;

587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
			if ( result != NULL ) {
				n_hits += result->hit;
				free(result);
			}

			rval = fgets(line, 1023, fh);
			if ( rval == NULL ) break;
			chomp(line);
			prefixed = malloc(1024);
			snprintf(prefixed, 1023, "%s%s", prefix, line);

			pargs = worker_args[i];
			free(pargs->filename);
			pargs->filename = prefixed;
			/* Other arguments unchanged */

Thomas White's avatar
Thomas White committed
603
			worker_active[i] = 1;
604
605
606
			r = pthread_create(&workers[i], NULL, process_image,
			                   pargs);
			if ( r != 0 ) {
Thomas White's avatar
Thomas White committed
607
				worker_active[i] = 0;
608
609
610
611
612
				ERROR("Couldn't start thread %i\n", i);
			}

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

Thomas White's avatar
Thomas White committed
614
615
	} while ( rval != NULL );

Thomas White's avatar
Thomas White committed
616
	/* Catch all remaining threads */
617
	for ( i=0; i<nthreads; i++ ) {
Thomas White's avatar
Thomas White committed
618
619
620
621
622
623
624
625
626
627
628
629
630
631

		struct process_result *result = NULL;

		if ( !worker_active[i] ) continue;

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

		worker_active[i] = 0;

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

632
633
		free(worker_args[i]->filename);
		free(worker_args[i]);
Thomas White's avatar
Thomas White committed
634

635
636
637
638
	}

	free(prefix);
	free(cell);
Thomas White's avatar
Thomas White committed
639
640
641
642
643
	fclose(fh);

	STATUS("There were %i images.\n", n_images);
	STATUS("%i hits were found.\n", n_hits);

Thomas White's avatar
Thomas White committed
644
645
646
647
	if ( gctx != NULL ) {
		cleanup_gpu(gctx);
	}

Thomas White's avatar
Thomas White committed
648
649
	return 0;
}