indexamajig.c 6.63 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
50
"      --indexing=<method> Use 'method' for indexing.  Choose from:\n"
"                           none     : no indexing\n"
"                           dirax    : invoke DirAx\n"
"      --write-drx         Write 'xfel.drx' for visualisation of reciprocal\n"
Thomas White's avatar
Thomas White committed
51
"                           space.  Implied by any indexing method other than\n"
Thomas White's avatar
Thomas White committed
52
"                           'none'.\n"
Thomas White's avatar
Thomas White committed
53
"      --dump-peaks        Write the results of the peak search to stdout.\n"
54
55
56
"      --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
57
"      --gpu               Use the GPU to speed up the simulation.\n"
Thomas White's avatar
Thomas White committed
58
59
"      --clean-image       Perform common-mode noise subtraction and\n"
"                           background removal on images before proceeding.\n"
60
61
62
"      --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
63
);
Thomas White's avatar
Thomas White committed
64
65
66
67
68
69
}


int main(int argc, char *argv[])
{
	int c;
Thomas White's avatar
Thomas White committed
70
	struct gpu_context *gctx = NULL;
Thomas White's avatar
Thomas White committed
71
72
73
74
75
	char *filename = NULL;
	FILE *fh;
	char *rval;
	int n_images;
	int n_hits;
76
	int config_noindex = 0;
Thomas White's avatar
Thomas White committed
77
	int config_dumpfound = 0;
78
	int config_nearbragg = 0;
Thomas White's avatar
Thomas White committed
79
	int config_writedrx = 0;
80
	int config_simulate = 0;
Thomas White's avatar
Thomas White committed
81
	int config_clean = 0;
82
	int config_nomatch = 0;
Thomas White's avatar
Thomas White committed
83
	int config_gpu = 0;
Thomas White's avatar
Thomas White committed
84
85
	IndexingMethod indm;
	char *indm_str = NULL;
Thomas White's avatar
Thomas White committed
86
87
88
89
90

	/* Long options */
	const struct option longopts[] = {
		{"help",               0, NULL,               'h'},
		{"input",              1, NULL,               'i'},
Thomas White's avatar
Thomas White committed
91
		{"gpu",                0, &config_gpu,         1},
92
		{"no-index",           0, &config_noindex,     1},
Thomas White's avatar
Thomas White committed
93
		{"dump-peaks",         0, &config_dumpfound,   1},
94
		{"near-bragg",         0, &config_nearbragg,   1},
Thomas White's avatar
Thomas White committed
95
96
		{"write-drx",          0, &config_writedrx,    1},
		{"indexing",           1, NULL,               'z'},
97
		{"simulate",           0, &config_simulate,    1},
Thomas White's avatar
Thomas White committed
98
		{"clean-image",        0, &config_clean,       1},
99
		{"no-match",           0, &config_nomatch,     1},
Thomas White's avatar
Thomas White committed
100
101
102
103
		{0, 0, NULL, 0}
	};

	/* Short options */
104
	while ((c = getopt_long(argc, argv, "hi:w", longopts, NULL)) != -1) {
Thomas White's avatar
Thomas White committed
105
106
107
108
109
110
111
112
113
114
115
116

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

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

Thomas White's avatar
Thomas White committed
117
118
119
120
121
		case 'z' : {
			indm_str = strdup(optarg);
			break;
		}

Thomas White's avatar
Thomas White committed
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
		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
147
148
149
150
151
152
153
154
155
156
157
158
159
160
	if ( indm_str == NULL ) {
		STATUS("You didn't specify an indexing method, so I won't"
		       " try to index anything.  If that isn't what you\n"
		       " wanted, re-run with --indexing=<method>.\n");
		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;
	}

Thomas White's avatar
Thomas White committed
161
162
163
164
165
166
167
168
169
	n_images = 0;
	n_hits = 0;
	do {

		char line[1024];
		struct hdfile *hdfile;
		struct image image;

		rval = fgets(line, 1023, fh);
Thomas White's avatar
Thomas White committed
170
		if ( rval == NULL ) continue;
Thomas White's avatar
Thomas White committed
171
172
		chomp(line);

Thomas White's avatar
Thomas White committed
173
		image.features = NULL;
174
		image.molecule = load_molecule();
Thomas White's avatar
Thomas White committed
175
		image.data = NULL;
176
		image.indexed_cell = NULL;
Thomas White's avatar
Thomas White committed
177

Thomas White's avatar
Thomas White committed
178
		#include "geometry-lcls.tmp"
179

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

182
183
		n_images++;

Thomas White's avatar
Thomas White committed
184
185
		hdfile = hdfile_open(line);
		if ( hdfile == NULL ) {
186
			continue;
Thomas White's avatar
Thomas White committed
187
188
		} else if ( hdfile_set_first_image(hdfile, "/") ) {
			ERROR("Couldn't select path\n");
189
			continue;
Thomas White's avatar
Thomas White committed
190
191
192
193
		}

		hdf5_read(hdfile, &image);

Thomas White's avatar
Thomas White committed
194
195
196
197
		if ( config_clean ) {
			clean_image(&image);
		}

198
		/* Perform 'fine' peak search */
199
		search_peaks(&image);
Thomas White's avatar
Thomas White committed
200

201
		if ( image_feature_count(image.features) > 5 ) {
Thomas White's avatar
Thomas White committed
202

203
204
			if ( config_dumpfound ) dump_peaks(&image);

Thomas White's avatar
Thomas White committed
205
206
207
208
209
			/* 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
210

211
			/* Calculate orientation matrix (by magic) */
Thomas White's avatar
Thomas White committed
212
			if ( config_writedrx || (indm != INDEXING_NONE) ) {
213
				index_pattern(&image, indm, config_nomatch);
Thomas White's avatar
Thomas White committed
214
215
216
			}

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

219
220
			n_hits++;

221
222
223
224
225
			/* Measure intensities if requested */
			if ( config_nearbragg ) {
				output_intensities(&image, image.indexed_cell);
			}

Thomas White's avatar
Thomas White committed
226
227
			/* Simulation or intensity measurements both require
			 * Ewald sphere vectors */
228
			if ( config_simulate ) {
229

Thomas White's avatar
Thomas White committed
230
231
				/* Simulate a diffraction pattern */
				image.twotheta = NULL;
232
				image.data = NULL;
Thomas White's avatar
Thomas White committed
233

234
235
236
237
238
239
				/* View head-on (unit cell is tilted) */
				image.orientation.w = 1.0;
				image.orientation.x = 0.0;
				image.orientation.y = 0.0;
				image.orientation.z = 0.0;

240
241
				image.molecule->cell = image.indexed_cell;

Thomas White's avatar
Thomas White committed
242
243
244
				if ( config_gpu ) {
					if ( gctx == NULL ) {
						gctx = setup_gpu(0, &image,
245
246
						                image.molecule,
						                24, 24, 40);
Thomas White's avatar
Thomas White committed
247
					}
248
					get_diffraction_gpu(gctx, &image);
Thomas White's avatar
Thomas White committed
249
				} else {
250
					get_diffraction(&image, 8, 8, 8, 0, 0);
Thomas White's avatar
Thomas White committed
251
				}
252
253
254
255
				if ( image.molecule == NULL ) {
					ERROR("Couldn't open molecule.pdb\n");
					return 1;
				}
256
				record_image(&image, 0);
257
258

				hdf5_write("simulated.h5", image.data,
259
				           image.width, image.height,
260
				           H5T_NATIVE_FLOAT);
261
262
263

			}

Thomas White's avatar
Thomas White committed
264
265
		}

Thomas White's avatar
Thomas White committed
266
done:
Thomas White's avatar
Thomas White committed
267
		free(image.data);
268
		image_feature_list_free(image.features);
Thomas White's avatar
Thomas White committed
269
		hdfile_close(hdfile);
270
		H5close();
Thomas White's avatar
Thomas White committed
271

Thomas White's avatar
Thomas White committed
272
273
274
275
276
277
278
	} 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
279
280
281
282
	if ( gctx != NULL ) {
		cleanup_gpu(gctx);
	}

Thomas White's avatar
Thomas White committed
283
284
	return 0;
}