indexamajig.c 8.86 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"
Thomas White's avatar
Thomas White committed
28
#include "intensities.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"
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"
Thomas White's avatar
Thomas White committed
52
"      --write-drx         Write 'xfel.drx' for visualisation of reciprocal\n"
Thomas White's avatar
Thomas White committed
53
"                           space.  Implied by any indexing method other than\n"
Thomas White's avatar
Thomas White committed
54
55
"                           'none'.  Beware: the units in this file are\n"
"                           reciprocal Angstroms.\n"
Thomas White's avatar
Thomas White committed
56
"      --dump-peaks        Write the results of the peak search to stdout.\n"
57
58
59
"      --near-bragg        Output a list of reflection intensities to stdout.\n"
"      --simulate          Simulate the diffraction pattern using the indexed\n"
"                           unit cell.\n"
Thomas White's avatar
Thomas White committed
60
"      --gpu               Use the GPU to speed up the simulation.\n"
Thomas White's avatar
Thomas White committed
61
62
"      --clean-image       Perform common-mode noise subtraction and\n"
"                           background removal on images before proceeding.\n"
63
64
65
"      --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
66
);
Thomas White's avatar
Thomas White committed
67
68
69
}


70
static struct image *get_simage(struct image *template, int alternate)
71
{
72
	struct image *image;
73
	struct panel panels[2];
74

75
76
	image = malloc(sizeof(*image));

77
	/* Simulate a diffraction pattern */
78
79
80
	image->twotheta = NULL;
	image->data = NULL;
	image->det = template->det;
81
82

	/* View head-on (unit cell is tilted) */
83
84
85
86
	image->orientation.w = 1.0;
	image->orientation.x = 0.0;
	image->orientation.y = 0.0;
	image->orientation.z = 0.0;
87
88
89

	/* Detector geometry for the simulation
	 * - not necessarily the same as the original. */
90
91
92
	image->width = 1024;
	image->height = 1024;
	image->det.n_panels = 2;
93

94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
	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;
117

118
119
120
121
122
123
	} else {

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

	}
124

125
126
	image->lambda = ph_en_to_lambda(eV_to_J(1.8e3));

127
128
129
	image->molecule = copy_molecule(template->molecule);
	free(image->molecule->cell);
	image->molecule->cell = cell_new_from_cell(template->indexed_cell);
130
131
132

	return image;
}
133
134


135
136
137
static void simulate_and_write(struct image *simage,
                               struct gpu_context **gctx)
{
138
	/* Set up GPU if necessary */
Thomas White's avatar
Thomas White committed
139
	if ( (gctx != NULL) && (*gctx == NULL) ) {
140
		*gctx = setup_gpu(0, simage, simage->molecule);
141
142
	}

Thomas White's avatar
Thomas White committed
143
	if ( (gctx != NULL) && (*gctx != NULL) ) {
144
		get_diffraction_gpu(*gctx, simage, 24, 24, 40);
145
	} else {
146
		get_diffraction(simage, 24, 24, 40, 0, 0);
147
	}
148
	if ( simage->molecule == NULL ) {
149
150
151
		ERROR("Couldn't open molecule.pdb\n");
		return;
	}
152
	record_image(simage, 0);
153

154
	hdf5_write("simulated.h5", simage->data, simage->width, simage->height,
155
156
157
158
		   H5T_NATIVE_FLOAT);
}


Thomas White's avatar
Thomas White committed
159
160
161
int main(int argc, char *argv[])
{
	int c;
Thomas White's avatar
Thomas White committed
162
	struct gpu_context *gctx = NULL;
Thomas White's avatar
Thomas White committed
163
164
165
166
167
	char *filename = NULL;
	FILE *fh;
	char *rval;
	int n_images;
	int n_hits;
168
	int config_noindex = 0;
Thomas White's avatar
Thomas White committed
169
	int config_dumpfound = 0;
170
	int config_nearbragg = 0;
Thomas White's avatar
Thomas White committed
171
	int config_writedrx = 0;
172
	int config_simulate = 0;
Thomas White's avatar
Thomas White committed
173
	int config_clean = 0;
174
	int config_nomatch = 0;
Thomas White's avatar
Thomas White committed
175
	int config_gpu = 0;
176
	int config_verbose = 0;
177
	int config_alternate = 0;
Thomas White's avatar
Thomas White committed
178
179
	IndexingMethod indm;
	char *indm_str = NULL;
180
	struct image image;
Thomas White's avatar
Thomas White committed
181
182
183
184
185

	/* Long options */
	const struct option longopts[] = {
		{"help",               0, NULL,               'h'},
		{"input",              1, NULL,               'i'},
Thomas White's avatar
Thomas White committed
186
		{"gpu",                0, &config_gpu,         1},
187
		{"no-index",           0, &config_noindex,     1},
Thomas White's avatar
Thomas White committed
188
		{"dump-peaks",         0, &config_dumpfound,   1},
189
		{"near-bragg",         0, &config_nearbragg,   1},
Thomas White's avatar
Thomas White committed
190
191
		{"write-drx",          0, &config_writedrx,    1},
		{"indexing",           1, NULL,               'z'},
192
		{"simulate",           0, &config_simulate,    1},
Thomas White's avatar
Thomas White committed
193
		{"clean-image",        0, &config_clean,       1},
194
		{"no-match",           0, &config_nomatch,     1},
195
		{"verbose",            0, &config_verbose,     1},
196
		{"alternate",          0, &config_alternate,   1},
Thomas White's avatar
Thomas White committed
197
198
199
200
		{0, 0, NULL, 0}
	};

	/* Short options */
201
	while ((c = getopt_long(argc, argv, "hi:w", longopts, NULL)) != -1) {
Thomas White's avatar
Thomas White committed
202
203
204
205
206
207
208
209
210
211
212
213

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

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

Thomas White's avatar
Thomas White committed
214
215
216
217
218
		case 'z' : {
			indm_str = strdup(optarg);
			break;
		}

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

Thomas White's avatar
Thomas White committed
244
245
	if ( indm_str == NULL ) {
		STATUS("You didn't specify an indexing method, so I won't"
246
247
248
		       " 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
249
250
251
252
253
254
255
256
257
		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;
	}
258
	free(indm_str);
Thomas White's avatar
Thomas White committed
259

260
261
262
263
264
265
266
	image.molecule = load_molecule();
	if ( image.molecule == NULL ) {
		ERROR("You must provide molecule.pdb in the working "
		      "directory.\n");
		return 1;
	}

Thomas White's avatar
Thomas White committed
267
268
269
270
271
272
	n_images = 0;
	n_hits = 0;
	do {

		char line[1024];
		struct hdfile *hdfile;
273
		struct image *simage;
Thomas White's avatar
Thomas White committed
274
275

		rval = fgets(line, 1023, fh);
Thomas White's avatar
Thomas White committed
276
		if ( rval == NULL ) continue;
Thomas White's avatar
Thomas White committed
277
278
		chomp(line);

Thomas White's avatar
Thomas White committed
279
		image.features = NULL;
Thomas White's avatar
Thomas White committed
280
		image.data = NULL;
281
		image.indexed_cell = NULL;
Thomas White's avatar
Thomas White committed
282

Thomas White's avatar
Thomas White committed
283
		#include "geometry-lcls.tmp"
284

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

287
288
		n_images++;

Thomas White's avatar
Thomas White committed
289
290
		hdfile = hdfile_open(line);
		if ( hdfile == NULL ) {
291
			continue;
Thomas White's avatar
Thomas White committed
292
293
		} else if ( hdfile_set_first_image(hdfile, "/") ) {
			ERROR("Couldn't select path\n");
294
			continue;
Thomas White's avatar
Thomas White committed
295
296
297
298
		}

		hdf5_read(hdfile, &image);

Thomas White's avatar
Thomas White committed
299
300
301
302
		if ( config_clean ) {
			clean_image(&image);
		}

303
		/* Perform 'fine' peak search */
304
		search_peaks(&image);
Thomas White's avatar
Thomas White committed
305

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

308
		if ( config_dumpfound ) dump_peaks(&image);
309

310
311
312
313
314
		/* 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
315

316
317
318
319
320
		/* Calculate orientation matrix (by magic) */
		if ( config_writedrx || (indm != INDEXING_NONE) ) {
			index_pattern(&image, indm, config_nomatch,
			              config_verbose);
		}
Thomas White's avatar
Thomas White committed
321

322
323
		/* No cell at this point?  Then we're done. */
		if ( image.indexed_cell == NULL ) goto done;
324

325
		n_hits++;
326

327
		simage = get_simage(&image, config_alternate);
328

329
330
		/* Measure intensities if requested */
		if ( config_nearbragg ) {
331
332
333
334
			/* Use original data (temporarily) */
			simage->data = image.data;
			output_intensities(simage, image.indexed_cell);
			simage->data = NULL;
335
		}
336

337
		/* Simulate if requested */
Thomas White's avatar
Thomas White committed
338
339
		if ( config_simulate ) {
			if ( config_gpu ) {
340
				simulate_and_write(simage, &gctx);
Thomas White's avatar
Thomas White committed
341
			} else {
342
				simulate_and_write(simage, NULL);
Thomas White's avatar
Thomas White committed
343
344
			}
		}
Thomas White's avatar
Thomas White committed
345

346
347
348
		/* Finished with alternate image */
		if ( simage->twotheta != NULL ) free(simage->twotheta);
		if ( simage->data != NULL ) free(simage->data);
349
		free_molecule(simage->molecule);
350
		free(simage);
351

352
353
354
		/* Only free cell if found */
		free(image.indexed_cell);

Thomas White's avatar
Thomas White committed
355
done:
Thomas White's avatar
Thomas White committed
356
		free(image.data);
357
		free(image.det.panels);
358
		image_feature_list_free(image.features);
Thomas White's avatar
Thomas White committed
359
		hdfile_close(hdfile);
360
		H5close();
Thomas White's avatar
Thomas White committed
361

Thomas White's avatar
Thomas White committed
362
363
364
	} while ( rval != NULL );

	fclose(fh);
365
	free_molecule(image.molecule);
Thomas White's avatar
Thomas White committed
366
367
368
369

	STATUS("There were %i images.\n", n_images);
	STATUS("%i hits were found.\n", n_hits);

Thomas White's avatar
Thomas White committed
370
371
372
373
	if ( gctx != NULL ) {
		cleanup_gpu(gctx);
	}

Thomas White's avatar
Thomas White committed
374
375
	return 0;
}