indexamajig.c 11.1 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"
Thomas White's avatar
Thomas White committed
83
);
Thomas White's avatar
Thomas White committed
84
85
86
}


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

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

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

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

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

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

135
136
137
138
139
140
	} else {

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

	}
141

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

	return image;
}
147
148


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

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

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


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

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

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

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

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

Thomas White's avatar
Thomas White committed
235
236
237
238
239
		case 'z' : {
			indm_str = strdup(optarg);
			break;
		}

240
241
242
243
244
		case 'q' : {
			intfile = strdup(optarg);
			break;
		}

Thomas White's avatar
Thomas White committed
245
246
247
248
249
		case 'p' : {
			pdb = strdup(optarg);
			break;
		}

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

275
	if ( intfile != NULL ) {
276
		intensities = read_reflections(intfile, counts);
277
278
	} else {
		intensities = NULL;
279
		counts = NULL;
280
281
	}

Thomas White's avatar
Thomas White committed
282
283
284
285
	if ( pdb == NULL ) {
		pdb = strdup("molecule.pdb");
	}

Thomas White's avatar
Thomas White committed
286
287
	if ( indm_str == NULL ) {
		STATUS("You didn't specify an indexing method, so I won't"
288
289
290
		       " 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
291
292
293
294
295
296
297
298
299
		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;
	}
300
	free(indm_str);
Thomas White's avatar
Thomas White committed
301

Thomas White's avatar
Thomas White committed
302
303
	cell = load_cell_from_pdb(pdb);
	free(pdb);
304
305
	if ( cell == NULL ) {
		ERROR("Couldn't read unit cell (from molecule.pdb)\n");
306
307
308
		return 1;
	}

Thomas White's avatar
Thomas White committed
309
310
311
312
313
314
	n_images = 0;
	n_hits = 0;
	do {

		char line[1024];
		struct hdfile *hdfile;
315
		struct image *simage;
316
317
		float *data_for_measurement;
		size_t data_size;
Thomas White's avatar
Thomas White committed
318
319

		rval = fgets(line, 1023, fh);
Thomas White's avatar
Thomas White committed
320
		if ( rval == NULL ) continue;
Thomas White's avatar
Thomas White committed
321
322
		chomp(line);

Thomas White's avatar
Thomas White committed
323
		image.features = NULL;
Thomas White's avatar
Thomas White committed
324
		image.data = NULL;
325
		image.indexed_cell = NULL;
Thomas White's avatar
Thomas White committed
326

Thomas White's avatar
Thomas White committed
327
		#include "geometry-lcls.tmp"
328

Thomas White's avatar
Thomas White committed
329
330
		STATUS("Processing '%s'\n", line);

331
332
		n_images++;

Thomas White's avatar
Thomas White committed
333
334
		hdfile = hdfile_open(line);
		if ( hdfile == NULL ) {
335
			continue;
Thomas White's avatar
Thomas White committed
336
337
		} else if ( hdfile_set_first_image(hdfile, "/") ) {
			ERROR("Couldn't select path\n");
338
			continue;
Thomas White's avatar
Thomas White committed
339
340
341
342
		}

		hdf5_read(hdfile, &image);

343
344
		if ( config_cmfilter ) {
			filter_cm(&image);
Thomas White's avatar
Thomas White committed
345
346
		}

347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
		/* 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;
			}
			}

		}

368
		/* Perform 'fine' peak search */
369
		search_peaks(&image);
Thomas White's avatar
Thomas White committed
370

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

373
		if ( config_dumpfound ) dump_peaks(&image);
374

375
376
377
378
379
		/* 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
380

381
382
		/* Calculate orientation matrix (by magic) */
		if ( config_writedrx || (indm != INDEXING_NONE) ) {
383
			index_pattern(&image, cell, indm, config_nomatch,
384
385
			              config_verbose);
		}
Thomas White's avatar
Thomas White committed
386

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

390
		n_hits++;
391

392
		simage = get_simage(&image, config_alternate);
393

394
395
		/* Measure intensities if requested */
		if ( config_nearbragg ) {
396
			/* Use original data (temporarily) */
397
			simage->data = data_for_measurement;
398
399
			output_intensities(simage, image.indexed_cell);
			simage->data = NULL;
400
		}
401

402
		/* Simulate if requested */
Thomas White's avatar
Thomas White committed
403
404
		if ( config_simulate ) {
			if ( config_gpu ) {
405
				simulate_and_write(simage, &gctx, intensities,
406
				                   counts, cell);
Thomas White's avatar
Thomas White committed
407
			} else {
408
				simulate_and_write(simage, NULL, intensities,
409
				                   counts, cell);
Thomas White's avatar
Thomas White committed
410
411
			}
		}
Thomas White's avatar
Thomas White committed
412

413
414
415
416
		/* Finished with alternate image */
		if ( simage->twotheta != NULL ) free(simage->twotheta);
		if ( simage->data != NULL ) free(simage->data);
		free(simage);
417

418
419
420
		/* Only free cell if found */
		free(image.indexed_cell);

Thomas White's avatar
Thomas White committed
421
done:
Thomas White's avatar
Thomas White committed
422
		free(image.data);
423
		free(image.det.panels);
424
		image_feature_list_free(image.features);
425
		free(data_for_measurement);
Thomas White's avatar
Thomas White committed
426
		hdfile_close(hdfile);
427
		H5close();
Thomas White's avatar
Thomas White committed
428

Thomas White's avatar
Thomas White committed
429
430
431
432
433
434
435
	} 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
436
437
438
439
	if ( gctx != NULL ) {
		cleanup_gpu(gctx);
	}

Thomas White's avatar
Thomas White committed
440
441
	return 0;
}