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


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

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

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

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

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

145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
	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;
168

169
170
171
172
173
174
	} else {

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

	}
175

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

	return image;
}
181
182


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

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

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


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

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

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

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

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

251
252
	#include "geometry-lcls.tmp"

253
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
	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);

330
331
332
	/* Free detector panel records */
	free(image.det.panels);

333
334
335
336
337
338
339
done:
	free(image.data);
	free(image.det.panels);
	image_feature_list_free(image.features);
	free(data_for_measurement);
	hdfile_close(hdfile);

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


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

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

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

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

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

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

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

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

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

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

Thomas White's avatar
Thomas White committed
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
		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;
	}

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

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

479
480
481
482
	if ( prefix == NULL ) {
		prefix = "";
	}

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

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

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

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

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

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

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

529
530
		n_images++;

531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
		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;

		r = pthread_create(&workers[i], NULL, process_image, pargs);
		if ( r != 0 ) {
			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;
			struct process_args *pargs;

			t.tv_sec = 0;
			t.tv_nsec = 20000; /* 20 ms */

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

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

			r = pthread_create(&workers[i], NULL, process_image,
			                   pargs);
			if ( r != 0 ) {
				ERROR("Couldn't start thread %i\n", i);
			}

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

Thomas White's avatar
Thomas White committed
603
604
	} while ( rval != NULL );

605
606
607
608
609
610
611
	for ( i=0; i<nthreads; i++ ) {
		free(worker_args[i]->filename);
		free(worker_args[i]);
	}

	free(prefix);
	free(cell);
Thomas White's avatar
Thomas White committed
612
613
614
615
616
	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
617
618
619
620
	if ( gctx != NULL ) {
		cleanup_gpu(gctx);
	}

Thomas White's avatar
Thomas White committed
621
622
	return 0;
}