indexamajig.c 10.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
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"
80
81
"     --intensities=<file>  Specify file containing reflection intensities\n"
"                            to use.\n"
Thomas White's avatar
Thomas White committed
82
);
Thomas White's avatar
Thomas White committed
83
84
85
}


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

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

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

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

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

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

134
135
136
137
138
139
	} else {

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

	}
140

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

	return image;
}
146
147


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

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

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


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

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

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

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

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

Thomas White's avatar
Thomas White committed
232
233
234
235
236
		case 'z' : {
			indm_str = strdup(optarg);
			break;
		}

237
238
239
240
241
		case 'q' : {
			intfile = strdup(optarg);
			break;
		}

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

267
	if ( intfile != NULL ) {
268
		intensities = read_reflections(intfile, counts);
269
270
	} else {
		intensities = NULL;
271
		counts = NULL;
272
273
	}

Thomas White's avatar
Thomas White committed
274
275
	if ( indm_str == NULL ) {
		STATUS("You didn't specify an indexing method, so I won't"
276
277
278
		       " 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
279
280
281
282
283
284
285
286
287
		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;
	}
288
	free(indm_str);
Thomas White's avatar
Thomas White committed
289

290
291
292
	cell = load_cell_from_pdb("molecule.pdb");
	if ( cell == NULL ) {
		ERROR("Couldn't read unit cell (from molecule.pdb)\n");
293
294
295
		return 1;
	}

Thomas White's avatar
Thomas White committed
296
297
298
299
300
301
	n_images = 0;
	n_hits = 0;
	do {

		char line[1024];
		struct hdfile *hdfile;
302
		struct image *simage;
303
304
		float *data_for_measurement;
		size_t data_size;
Thomas White's avatar
Thomas White committed
305
306

		rval = fgets(line, 1023, fh);
Thomas White's avatar
Thomas White committed
307
		if ( rval == NULL ) continue;
Thomas White's avatar
Thomas White committed
308
309
		chomp(line);

Thomas White's avatar
Thomas White committed
310
		image.features = NULL;
Thomas White's avatar
Thomas White committed
311
		image.data = NULL;
312
		image.indexed_cell = NULL;
Thomas White's avatar
Thomas White committed
313

Thomas White's avatar
Thomas White committed
314
		#include "geometry-lcls.tmp"
315

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

318
319
		n_images++;

Thomas White's avatar
Thomas White committed
320
321
		hdfile = hdfile_open(line);
		if ( hdfile == NULL ) {
322
			continue;
Thomas White's avatar
Thomas White committed
323
324
		} else if ( hdfile_set_first_image(hdfile, "/") ) {
			ERROR("Couldn't select path\n");
325
			continue;
Thomas White's avatar
Thomas White committed
326
327
328
329
		}

		hdf5_read(hdfile, &image);

330
331
		if ( config_cmfilter ) {
			filter_cm(&image);
Thomas White's avatar
Thomas White committed
332
333
		}

334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
		/* 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;
			}
			}

		}

355
		/* Perform 'fine' peak search */
356
		search_peaks(&image);
Thomas White's avatar
Thomas White committed
357

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

360
		if ( config_dumpfound ) dump_peaks(&image);
361

362
363
364
365
366
		/* 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
367

368
369
		/* Calculate orientation matrix (by magic) */
		if ( config_writedrx || (indm != INDEXING_NONE) ) {
370
			index_pattern(&image, cell, indm, config_nomatch,
371
372
			              config_verbose);
		}
Thomas White's avatar
Thomas White committed
373

374
375
		/* No cell at this point?  Then we're done. */
		if ( image.indexed_cell == NULL ) goto done;
376

377
		n_hits++;
378

379
		simage = get_simage(&image, config_alternate);
380

381
382
		/* Measure intensities if requested */
		if ( config_nearbragg ) {
383
			/* Use original data (temporarily) */
384
			simage->data = data_for_measurement;
385
386
			output_intensities(simage, image.indexed_cell);
			simage->data = NULL;
387
		}
388

389
		/* Simulate if requested */
Thomas White's avatar
Thomas White committed
390
391
		if ( config_simulate ) {
			if ( config_gpu ) {
392
				simulate_and_write(simage, &gctx, intensities,
393
				                   counts, cell);
Thomas White's avatar
Thomas White committed
394
			} else {
395
				simulate_and_write(simage, NULL, intensities,
396
				                   counts, cell);
Thomas White's avatar
Thomas White committed
397
398
			}
		}
Thomas White's avatar
Thomas White committed
399

400
401
402
403
		/* Finished with alternate image */
		if ( simage->twotheta != NULL ) free(simage->twotheta);
		if ( simage->data != NULL ) free(simage->data);
		free(simage);
404

405
406
407
		/* Only free cell if found */
		free(image.indexed_cell);

Thomas White's avatar
Thomas White committed
408
done:
Thomas White's avatar
Thomas White committed
409
		free(image.data);
410
		free(image.det.panels);
411
		image_feature_list_free(image.features);
412
		free(data_for_measurement);
Thomas White's avatar
Thomas White committed
413
		hdfile_close(hdfile);
414
		H5close();
Thomas White's avatar
Thomas White committed
415

Thomas White's avatar
Thomas White committed
416
417
418
419
420
421
422
	} 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
423
424
425
426
	if ( gctx != NULL ) {
		cleanup_gpu(gctx);
	}

Thomas White's avatar
Thomas White committed
427
428
	return 0;
}