indexamajig.c 11.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
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
25
26

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


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


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

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

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

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

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

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

136
137
138
139
140
141
	} else {

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

	}
142

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

	return image;
}
148
149


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

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

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


Thomas White's avatar
Thomas White committed
172
173
174
int main(int argc, char *argv[])
{
	int c;
Thomas White's avatar
Thomas White committed
175
	struct gpu_context *gctx = NULL;
Thomas White's avatar
Thomas White committed
176
177
178
179
180
	char *filename = NULL;
	FILE *fh;
	char *rval;
	int n_images;
	int n_hits;
181
	int config_noindex = 0;
Thomas White's avatar
Thomas White committed
182
	int config_dumpfound = 0;
183
	int config_nearbragg = 0;
Thomas White's avatar
Thomas White committed
184
	int config_writedrx = 0;
185
	int config_simulate = 0;
186
187
	int config_cmfilter = 0;
	int config_noisefilter = 0;
188
	int config_nomatch = 0;
Thomas White's avatar
Thomas White committed
189
	int config_gpu = 0;
190
	int config_verbose = 0;
191
	int config_alternate = 0;
Thomas White's avatar
Thomas White committed
192
193
	IndexingMethod indm;
	char *indm_str = NULL;
194
	struct image image;
195
196
197
	UnitCell *cell;
	double *intensities = NULL;
	char *intfile = NULL;
198
	unsigned int *counts = NULL;
Thomas White's avatar
Thomas White committed
199
	char *pdb = NULL;
200
	char *prefix = NULL;
Thomas White's avatar
Thomas White committed
201
202
203
204
205

	/* Long options */
	const struct option longopts[] = {
		{"help",               0, NULL,               'h'},
		{"input",              1, NULL,               'i'},
Thomas White's avatar
Thomas White committed
206
		{"gpu",                0, &config_gpu,         1},
207
		{"no-index",           0, &config_noindex,     1},
Thomas White's avatar
Thomas White committed
208
		{"dump-peaks",         0, &config_dumpfound,   1},
209
		{"near-bragg",         0, &config_nearbragg,   1},
Thomas White's avatar
Thomas White committed
210
211
		{"write-drx",          0, &config_writedrx,    1},
		{"indexing",           1, NULL,               'z'},
212
		{"simulate",           0, &config_simulate,    1},
213
214
		{"filter-cm",          0, &config_cmfilter,    1},
		{"filter-noise",       0, &config_noisefilter, 1},
215
		{"no-match",           0, &config_nomatch,     1},
216
		{"verbose",            0, &config_verbose,     1},
217
		{"alternate",          0, &config_alternate,   1},
218
		{"intensities",        1, NULL,               'q'},
Thomas White's avatar
Thomas White committed
219
		{"pdb",                1, NULL,               'p'},
220
		{"prefix",             1, NULL,               'x'},
Thomas White's avatar
Thomas White committed
221
222
223
224
		{0, 0, NULL, 0}
	};

	/* Short options */
Thomas White's avatar
Thomas White committed
225
	while ((c = getopt_long(argc, argv, "hi:wp:", longopts, NULL)) != -1) {
Thomas White's avatar
Thomas White committed
226
227
228
229
230
231
232
233
234
235
236
237

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

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

Thomas White's avatar
Thomas White committed
238
239
240
241
242
		case 'z' : {
			indm_str = strdup(optarg);
			break;
		}

243
244
245
246
247
		case 'q' : {
			intfile = strdup(optarg);
			break;
		}

Thomas White's avatar
Thomas White committed
248
249
250
251
252
		case 'p' : {
			pdb = strdup(optarg);
			break;
		}

253
254
255
256
257
		case 'x' : {
			prefix = strdup(optarg);
			break;
		}

Thomas White's avatar
Thomas White committed
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
		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;
	}

283
	if ( intfile != NULL ) {
284
		intensities = read_reflections(intfile, counts);
285
286
	} else {
		intensities = NULL;
287
		counts = NULL;
288
289
	}

Thomas White's avatar
Thomas White committed
290
291
292
293
	if ( pdb == NULL ) {
		pdb = strdup("molecule.pdb");
	}

294
295
296
297
	if ( prefix == NULL ) {
		prefix = "";
	}

Thomas White's avatar
Thomas White committed
298
299
	if ( indm_str == NULL ) {
		STATUS("You didn't specify an indexing method, so I won't"
300
301
302
		       " 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
303
304
305
306
307
308
309
310
311
		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;
	}
312
	free(indm_str);
Thomas White's avatar
Thomas White committed
313

Thomas White's avatar
Thomas White committed
314
315
	cell = load_cell_from_pdb(pdb);
	free(pdb);
316
317
	if ( cell == NULL ) {
		ERROR("Couldn't read unit cell (from molecule.pdb)\n");
318
319
320
		return 1;
	}

Thomas White's avatar
Thomas White committed
321
322
323
324
325
326
	n_images = 0;
	n_hits = 0;
	do {

		char line[1024];
		struct hdfile *hdfile;
327
		struct image *simage;
328
329
		float *data_for_measurement;
		size_t data_size;
330
		char prefixed[1024];
Thomas White's avatar
Thomas White committed
331
332

		rval = fgets(line, 1023, fh);
Thomas White's avatar
Thomas White committed
333
		if ( rval == NULL ) continue;
Thomas White's avatar
Thomas White committed
334
335
		chomp(line);

Thomas White's avatar
Thomas White committed
336
		image.features = NULL;
Thomas White's avatar
Thomas White committed
337
		image.data = NULL;
338
		image.indexed_cell = NULL;
Thomas White's avatar
Thomas White committed
339

Thomas White's avatar
Thomas White committed
340
		#include "geometry-lcls.tmp"
341

342
343
344
		snprintf(prefixed, 1023, "%s%s", prefix, line);

		STATUS("Processing '%s'\n", prefixed);
Thomas White's avatar
Thomas White committed
345

346
347
		n_images++;

348
		hdfile = hdfile_open(prefixed);
Thomas White's avatar
Thomas White committed
349
		if ( hdfile == NULL ) {
350
			continue;
Thomas White's avatar
Thomas White committed
351
352
		} else if ( hdfile_set_first_image(hdfile, "/") ) {
			ERROR("Couldn't select path\n");
353
			continue;
Thomas White's avatar
Thomas White committed
354
355
356
357
		}

		hdf5_read(hdfile, &image);

358
359
		if ( config_cmfilter ) {
			filter_cm(&image);
Thomas White's avatar
Thomas White committed
360
361
		}

362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
		/* 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;
			}
			}

		}

383
		/* Perform 'fine' peak search */
384
		search_peaks(&image);
Thomas White's avatar
Thomas White committed
385

386
		if ( image_feature_count(image.features) < 5 ) goto done;
Thomas White's avatar
Thomas White committed
387

388
		if ( config_dumpfound ) dump_peaks(&image);
389

390
391
392
393
394
		/* Not indexing nor writing xfel.drx?
		 * Then there's nothing left to do. */
		if ( (!config_writedrx) && (indm == INDEXING_NONE) ) {
			goto done;
		}
Thomas White's avatar
Thomas White committed
395

396
397
		/* Calculate orientation matrix (by magic) */
		if ( config_writedrx || (indm != INDEXING_NONE) ) {
398
			index_pattern(&image, cell, indm, config_nomatch,
399
400
			              config_verbose);
		}
Thomas White's avatar
Thomas White committed
401

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

405
		n_hits++;
406

407
		simage = get_simage(&image, config_alternate);
408

409
410
		/* Measure intensities if requested */
		if ( config_nearbragg ) {
411
			/* Use original data (temporarily) */
412
			simage->data = data_for_measurement;
413
414
			output_intensities(simage, image.indexed_cell);
			simage->data = NULL;
415
		}
416

417
		/* Simulate if requested */
Thomas White's avatar
Thomas White committed
418
419
		if ( config_simulate ) {
			if ( config_gpu ) {
420
				simulate_and_write(simage, &gctx, intensities,
421
				                   counts, cell);
Thomas White's avatar
Thomas White committed
422
			} else {
423
				simulate_and_write(simage, NULL, intensities,
424
				                   counts, cell);
Thomas White's avatar
Thomas White committed
425
426
			}
		}
Thomas White's avatar
Thomas White committed
427

428
429
430
431
		/* Finished with alternate image */
		if ( simage->twotheta != NULL ) free(simage->twotheta);
		if ( simage->data != NULL ) free(simage->data);
		free(simage);
432

433
434
435
		/* Only free cell if found */
		free(image.indexed_cell);

Thomas White's avatar
Thomas White committed
436
done:
Thomas White's avatar
Thomas White committed
437
		free(image.data);
438
		free(image.det.panels);
439
		image_feature_list_free(image.features);
440
		free(data_for_measurement);
Thomas White's avatar
Thomas White committed
441
		hdfile_close(hdfile);
442
		H5close();
Thomas White's avatar
Thomas White committed
443

Thomas White's avatar
Thomas White committed
444
445
446
447
448
449
450
	} 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
451
452
453
454
	if ( gctx != NULL ) {
		cleanup_gpu(gctx);
	}

Thomas White's avatar
Thomas White committed
455
456
	return 0;
}