indexamajig.c 10.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"
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"
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
);
Thomas White's avatar
Thomas White committed
81
82
83
}


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

89
90
	image = malloc(sizeof(*image));

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

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

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

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

132
133
134
135
136
137
	} else {

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

	}
138

139
140
	image->lambda = ph_en_to_lambda(eV_to_J(1.8e3));

141
142
143
	image->molecule = copy_molecule(template->molecule);
	free(image->molecule->cell);
	image->molecule->cell = cell_new_from_cell(template->indexed_cell);
144
145
146

	return image;
}
147
148


149
150
151
static void simulate_and_write(struct image *simage,
                               struct gpu_context **gctx)
{
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, simage->molecule);
155
156
	}

Thomas White's avatar
Thomas White committed
157
	if ( (gctx != NULL) && (*gctx != NULL) ) {
158
		get_diffraction_gpu(*gctx, simage, 24, 24, 40);
159
	} else {
160
		get_diffraction(simage, 24, 24, 40, 0, 0);
161
	}
162
	if ( simage->molecule == NULL ) {
163
164
165
		ERROR("Couldn't open molecule.pdb\n");
		return;
	}
166
	record_image(simage, 0);
167

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


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

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

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

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

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

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

Thomas White's avatar
Thomas White committed
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
		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
260
261
	if ( indm_str == NULL ) {
		STATUS("You didn't specify an indexing method, so I won't"
262
263
264
		       " 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
265
266
267
268
269
270
271
272
273
		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;
	}
274
	free(indm_str);
Thomas White's avatar
Thomas White committed
275

276
277
278
279
280
281
282
	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
283
284
285
286
287
288
	n_images = 0;
	n_hits = 0;
	do {

		char line[1024];
		struct hdfile *hdfile;
289
		struct image *simage;
290
291
		float *data_for_measurement;
		size_t data_size;
Thomas White's avatar
Thomas White committed
292
293

		rval = fgets(line, 1023, fh);
Thomas White's avatar
Thomas White committed
294
		if ( rval == NULL ) continue;
Thomas White's avatar
Thomas White committed
295
296
		chomp(line);

Thomas White's avatar
Thomas White committed
297
		image.features = NULL;
Thomas White's avatar
Thomas White committed
298
		image.data = NULL;
299
		image.indexed_cell = NULL;
Thomas White's avatar
Thomas White committed
300

Thomas White's avatar
Thomas White committed
301
		#include "geometry-lcls.tmp"
302

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

305
306
		n_images++;

Thomas White's avatar
Thomas White committed
307
308
		hdfile = hdfile_open(line);
		if ( hdfile == NULL ) {
309
			continue;
Thomas White's avatar
Thomas White committed
310
311
		} else if ( hdfile_set_first_image(hdfile, "/") ) {
			ERROR("Couldn't select path\n");
312
			continue;
Thomas White's avatar
Thomas White committed
313
314
315
316
		}

		hdf5_read(hdfile, &image);

317
318
		if ( config_cmfilter ) {
			filter_cm(&image);
Thomas White's avatar
Thomas White committed
319
320
		}

321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
		/* 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;
			}
			}

		}

342
		/* Perform 'fine' peak search */
343
		search_peaks(&image);
Thomas White's avatar
Thomas White committed
344

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

347
		if ( config_dumpfound ) dump_peaks(&image);
348

349
350
351
352
353
		/* 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
354

355
356
357
358
359
		/* 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
360

361
362
		/* No cell at this point?  Then we're done. */
		if ( image.indexed_cell == NULL ) goto done;
363

364
		n_hits++;
365

366
		simage = get_simage(&image, config_alternate);
367

368
369
		/* Measure intensities if requested */
		if ( config_nearbragg ) {
370
			/* Use original data (temporarily) */
371
			simage->data = data_for_measurement;
372
373
			output_intensities(simage, image.indexed_cell);
			simage->data = NULL;
374
		}
375

376
		/* Simulate if requested */
Thomas White's avatar
Thomas White committed
377
378
		if ( config_simulate ) {
			if ( config_gpu ) {
379
				simulate_and_write(simage, &gctx);
Thomas White's avatar
Thomas White committed
380
			} else {
381
				simulate_and_write(simage, NULL);
Thomas White's avatar
Thomas White committed
382
383
			}
		}
Thomas White's avatar
Thomas White committed
384

385
386
387
		/* Finished with alternate image */
		if ( simage->twotheta != NULL ) free(simage->twotheta);
		if ( simage->data != NULL ) free(simage->data);
388
		free_molecule(simage->molecule);
389
		free(simage);
390

391
392
393
		/* Only free cell if found */
		free(image.indexed_cell);

Thomas White's avatar
Thomas White committed
394
done:
Thomas White's avatar
Thomas White committed
395
		free(image.data);
396
		free(image.det.panels);
397
		image_feature_list_free(image.features);
398
		free(data_for_measurement);
Thomas White's avatar
Thomas White committed
399
		hdfile_close(hdfile);
400
		H5close();
Thomas White's avatar
Thomas White committed
401

Thomas White's avatar
Thomas White committed
402
403
404
	} while ( rval != NULL );

	fclose(fh);
405
	free_molecule(image.molecule);
Thomas White's avatar
Thomas White committed
406
407
408
409

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

Thomas White's avatar
Thomas White committed
410
411
412
413
	if ( gctx != NULL ) {
		cleanup_gpu(gctx);
	}

Thomas White's avatar
Thomas White committed
414
415
	return 0;
}