indexamajig.c 6.89 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
71
72
}


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

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

	/* Short options */
109
	while ((c = getopt_long(argc, argv, "hi:w", longopts, NULL)) != -1) {
Thomas White's avatar
Thomas White committed
110
111
112
113
114
115
116
117
118
119
120
121

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

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

Thomas White's avatar
Thomas White committed
122
123
124
125
126
		case 'z' : {
			indm_str = strdup(optarg);
			break;
		}

Thomas White's avatar
Thomas White committed
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
		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
152
153
154
155
156
157
158
159
160
161
162
163
164
165
	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
166
167
168
169
170
171
172
173
174
	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
175
		if ( rval == NULL ) continue;
Thomas White's avatar
Thomas White committed
176
177
		chomp(line);

Thomas White's avatar
Thomas White committed
178
		image.features = NULL;
179
		image.molecule = load_molecule();
Thomas White's avatar
Thomas White committed
180
		image.data = NULL;
181
		image.indexed_cell = NULL;
Thomas White's avatar
Thomas White committed
182

Thomas White's avatar
Thomas White committed
183
		#include "geometry-lcls.tmp"
184

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

187
188
		n_images++;

Thomas White's avatar
Thomas White committed
189
190
		hdfile = hdfile_open(line);
		if ( hdfile == NULL ) {
191
			continue;
Thomas White's avatar
Thomas White committed
192
193
		} else if ( hdfile_set_first_image(hdfile, "/") ) {
			ERROR("Couldn't select path\n");
194
			continue;
Thomas White's avatar
Thomas White committed
195
196
197
198
		}

		hdf5_read(hdfile, &image);

Thomas White's avatar
Thomas White committed
199
200
201
202
		if ( config_clean ) {
			clean_image(&image);
		}

203
		/* Perform 'fine' peak search */
204
		search_peaks(&image);
Thomas White's avatar
Thomas White committed
205

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

208
209
			if ( config_dumpfound ) dump_peaks(&image);

Thomas White's avatar
Thomas White committed
210
211
212
213
214
			/* 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
215

216
			/* Calculate orientation matrix (by magic) */
Thomas White's avatar
Thomas White committed
217
			if ( config_writedrx || (indm != INDEXING_NONE) ) {
218
219
				index_pattern(&image, indm, config_nomatch,
				              config_verbose);
Thomas White's avatar
Thomas White committed
220
221
222
			}

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

225
226
			n_hits++;

227
228
229
230
231
			/* Measure intensities if requested */
			if ( config_nearbragg ) {
				output_intensities(&image, image.indexed_cell);
			}

Thomas White's avatar
Thomas White committed
232
233
			/* Simulation or intensity measurements both require
			 * Ewald sphere vectors */
234
			if ( config_simulate ) {
235

Thomas White's avatar
Thomas White committed
236
237
				/* Simulate a diffraction pattern */
				image.twotheta = NULL;
238
				image.data = NULL;
Thomas White's avatar
Thomas White committed
239

240
241
242
243
244
245
				/* 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;

246
247
				image.molecule->cell = image.indexed_cell;

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

				hdf5_write("simulated.h5", image.data,
265
				           image.width, image.height,
266
				           H5T_NATIVE_FLOAT);
267
268
269

			}

Thomas White's avatar
Thomas White committed
270
271
		}

Thomas White's avatar
Thomas White committed
272
done:
Thomas White's avatar
Thomas White committed
273
		free(image.data);
274
		image_feature_list_free(image.features);
Thomas White's avatar
Thomas White committed
275
		hdfile_close(hdfile);
276
		H5close();
Thomas White's avatar
Thomas White committed
277

Thomas White's avatar
Thomas White committed
278
279
280
281
282
283
284
	} 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
285
286
287
288
	if ( gctx != NULL ) {
		cleanup_gpu(gctx);
	}

Thomas White's avatar
Thomas White committed
289
290
	return 0;
}