indexamajig.c 12.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
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>
Thomas White's avatar
Thomas White committed
25
26
27

#include "utils.h"
#include "hdf5-file.h"
Thomas White's avatar
Thomas White committed
28
#include "index.h"
29
#include "peaks.h"
30
#include "diffraction.h"
Thomas White's avatar
Thomas White committed
31
#include "diffraction-gpu.h"
32
#include "detector.h"
33
#include "sfac.h"
Thomas White's avatar
Thomas White committed
34
#include "filters.h"
35
#include "reflections.h"
Thomas White's avatar
Thomas White committed
36
37
38
39
40
41
42
43
44
45
46
47


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
48
49
50
"      --indexing=<method> Use 'method' for indexing.  Choose from:\n"
"                           none     : no indexing\n"
"                           dirax    : invoke DirAx\n"
51
52
"\n"
"      --verbose           Be verbose about indexing.\n"
53
54
"      --gpu               Use the GPU to speed up the simulation.\n"
"\n"
55
"      --near-bragg        Output a list of reflection intensities to stdout.\n"
56
57
58
59
60
"                           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"
61
62
"      --simulate          Simulate the diffraction pattern using the indexed\n"
"                           unit cell.\n"
63
64
65
66
67
68
69
"      --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"
70
71
72
73
74
75
"\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"
76
77
"                           The intensities in this list are from the\n"
"                           centroid/integration procedure.\n"
78
79
80
"      --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
81
82
83
"     --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"
84
" -x, --prefix=<p>         Prefix filenames from input file with 'p'.\n"
Thomas White's avatar
Thomas White committed
85
);
Thomas White's avatar
Thomas White committed
86
87
88
}


89
static struct image *get_simage(struct image *template, int alternate)
90
{
91
	struct image *image;
92
	struct panel panels[2];
93

94
95
	image = malloc(sizeof(*image));

96
	/* Simulate a diffraction pattern */
97
98
99
	image->twotheta = NULL;
	image->data = NULL;
	image->det = template->det;
100
101

	/* View head-on (unit cell is tilted) */
102
103
104
105
	image->orientation.w = 1.0;
	image->orientation.x = 0.0;
	image->orientation.y = 0.0;
	image->orientation.z = 0.0;
106
107
108

	/* Detector geometry for the simulation
	 * - not necessarily the same as the original. */
109
110
111
	image->width = 1024;
	image->height = 1024;
	image->det.n_panels = 2;
112

113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
	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;
136

137
138
139
140
141
142
	} else {

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

	}
143

144
	image->lambda = ph_en_to_lambda(eV_to_J(1.8e3));
145
	image->features = template->features;
146
147
148

	return image;
}
149
150


151
static void simulate_and_write(struct image *simage, struct gpu_context **gctx,
152
153
                               double *intensities, unsigned int *counts,
                               UnitCell *cell)
154
{
155
	/* Set up GPU if necessary */
Thomas White's avatar
Thomas White committed
156
	if ( (gctx != NULL) && (*gctx == NULL) ) {
157
		*gctx = setup_gpu(0, simage, intensities, counts);
158
159
	}

Thomas White's avatar
Thomas White committed
160
	if ( (gctx != NULL) && (*gctx != NULL) ) {
161
		get_diffraction_gpu(*gctx, simage, 24, 24, 40, cell);
162
	} else {
163
164
		get_diffraction(simage, 24, 24, 40,
		                intensities, counts, cell, 0);
165
	}
166
	record_image(simage, 0);
167

168
	hdf5_write("simulated.h5", simage->data, simage->width, simage->height,
169
170
171
172
		   H5T_NATIVE_FLOAT);
}


173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
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
static int process_pattern(const 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, double *intensities,
		          unsigned int *counts, struct gpu_context *gctx)
{
	struct hdfile *hdfile;
	struct image image;
	struct image *simage;
	float *data_for_measurement;
	size_t data_size;

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

	#include "geometry-lcls.tmp"

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

	hdfile = hdfile_open(filename);
	if ( hdfile == NULL ) {
		return 0;
	} else if ( hdfile_set_first_image(hdfile, "/") ) {
		ERROR("Couldn't select path\n");
		return 0;
	}

	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);
	H5close();

	if ( image.indexed_cell == NULL ) return 0;
	return 1;
}


Thomas White's avatar
Thomas White committed
294
295
296
int main(int argc, char *argv[])
{
	int c;
Thomas White's avatar
Thomas White committed
297
	struct gpu_context *gctx = NULL;
Thomas White's avatar
Thomas White committed
298
299
300
301
302
	char *filename = NULL;
	FILE *fh;
	char *rval;
	int n_images;
	int n_hits;
303
	int config_noindex = 0;
Thomas White's avatar
Thomas White committed
304
	int config_dumpfound = 0;
305
	int config_nearbragg = 0;
Thomas White's avatar
Thomas White committed
306
	int config_writedrx = 0;
307
	int config_simulate = 0;
308
309
	int config_cmfilter = 0;
	int config_noisefilter = 0;
310
	int config_nomatch = 0;
Thomas White's avatar
Thomas White committed
311
	int config_gpu = 0;
312
	int config_verbose = 0;
313
	int config_alternate = 0;
Thomas White's avatar
Thomas White committed
314
315
	IndexingMethod indm;
	char *indm_str = NULL;
316
317
318
	UnitCell *cell;
	double *intensities = NULL;
	char *intfile = NULL;
319
	unsigned int *counts = NULL;
Thomas White's avatar
Thomas White committed
320
	char *pdb = NULL;
321
	char *prefix = NULL;
Thomas White's avatar
Thomas White committed
322
323
324
325
326

	/* Long options */
	const struct option longopts[] = {
		{"help",               0, NULL,               'h'},
		{"input",              1, NULL,               'i'},
Thomas White's avatar
Thomas White committed
327
		{"gpu",                0, &config_gpu,         1},
328
		{"no-index",           0, &config_noindex,     1},
Thomas White's avatar
Thomas White committed
329
		{"dump-peaks",         0, &config_dumpfound,   1},
330
		{"near-bragg",         0, &config_nearbragg,   1},
Thomas White's avatar
Thomas White committed
331
332
		{"write-drx",          0, &config_writedrx,    1},
		{"indexing",           1, NULL,               'z'},
333
		{"simulate",           0, &config_simulate,    1},
334
335
		{"filter-cm",          0, &config_cmfilter,    1},
		{"filter-noise",       0, &config_noisefilter, 1},
336
		{"no-match",           0, &config_nomatch,     1},
337
		{"verbose",            0, &config_verbose,     1},
338
		{"alternate",          0, &config_alternate,   1},
339
		{"intensities",        1, NULL,               'q'},
Thomas White's avatar
Thomas White committed
340
		{"pdb",                1, NULL,               'p'},
341
		{"prefix",             1, NULL,               'x'},
Thomas White's avatar
Thomas White committed
342
343
344
345
		{0, 0, NULL, 0}
	};

	/* Short options */
Thomas White's avatar
Thomas White committed
346
	while ((c = getopt_long(argc, argv, "hi:wp:", longopts, NULL)) != -1) {
Thomas White's avatar
Thomas White committed
347
348
349
350
351
352
353
354
355
356
357
358

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

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

Thomas White's avatar
Thomas White committed
359
360
361
362
363
		case 'z' : {
			indm_str = strdup(optarg);
			break;
		}

364
365
366
367
368
		case 'q' : {
			intfile = strdup(optarg);
			break;
		}

Thomas White's avatar
Thomas White committed
369
370
371
372
373
		case 'p' : {
			pdb = strdup(optarg);
			break;
		}

374
375
376
377
378
		case 'x' : {
			prefix = strdup(optarg);
			break;
		}

Thomas White's avatar
Thomas White committed
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
		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;
	}

404
	if ( intfile != NULL ) {
405
		intensities = read_reflections(intfile, counts);
406
407
	} else {
		intensities = NULL;
408
		counts = NULL;
409
410
	}

Thomas White's avatar
Thomas White committed
411
412
413
414
	if ( pdb == NULL ) {
		pdb = strdup("molecule.pdb");
	}

415
416
417
418
	if ( prefix == NULL ) {
		prefix = "";
	}

Thomas White's avatar
Thomas White committed
419
420
	if ( indm_str == NULL ) {
		STATUS("You didn't specify an indexing method, so I won't"
421
422
423
		       " 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
424
425
426
427
428
429
430
431
432
		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;
	}
433
	free(indm_str);
Thomas White's avatar
Thomas White committed
434

Thomas White's avatar
Thomas White committed
435
436
	cell = load_cell_from_pdb(pdb);
	free(pdb);
437
438
	if ( cell == NULL ) {
		ERROR("Couldn't read unit cell (from molecule.pdb)\n");
439
440
441
		return 1;
	}

Thomas White's avatar
Thomas White committed
442
443
	gsl_set_error_handler_off();

Thomas White's avatar
Thomas White committed
444
445
446
447
448
	n_images = 0;
	n_hits = 0;
	do {

		char line[1024];
449
		char prefixed[1024];
Thomas White's avatar
Thomas White committed
450
451

		rval = fgets(line, 1023, fh);
Thomas White's avatar
Thomas White committed
452
		if ( rval == NULL ) continue;
Thomas White's avatar
Thomas White committed
453
		chomp(line);
454
455
		snprintf(prefixed, 1023, "%s%s", prefix, line);

456
457
		n_images++;

458
459
460
461
462
463
464
		n_hits += process_pattern(line, cell, config_cmfilter,
		                          config_noisefilter,
		                          config_writedrx, config_dumpfound,
		                          config_verbose, config_alternate,
		                          config_nearbragg, config_gpu,
		                          config_simulate, config_nomatch, indm,
		                          intensities, counts, gctx);
Thomas White's avatar
Thomas White committed
465

Thomas White's avatar
Thomas White committed
466
467
468
469
470
471
472
	} while ( rval != NULL );

	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
473
474
475
476
	if ( gctx != NULL ) {
		cleanup_gpu(gctx);
	}

Thomas White's avatar
Thomas White committed
477
478
	return 0;
}