indexamajig.c 6.66 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
51
52
"      --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"
"                           space.  Implied by any indexing method other than"
"                           '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
58
"      --clean-image       Perform common-mode noise subtraction and\n"
"                           background removal on images before proceeding.\n"
59
60
61
"      --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
62
63
64
65
66
67
68
"\n");
}


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

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

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

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

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

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

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

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

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

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

181
182
		n_images++;

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

		hdf5_read(hdfile, &image);

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

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

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

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

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

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

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

218
219
			n_hits++;

Thomas White's avatar
Thomas White committed
220
221
			/* Simulation or intensity measurements both require
			 * Ewald sphere vectors */
222
			if ( config_nearbragg || config_simulate ) {
223

Thomas White's avatar
Thomas White committed
224
225
226
227
228
				/* Simulate a diffraction pattern */
				image.sfacs = NULL;
				image.twotheta = NULL;
				image.hdr = NULL;

229
230
231
232
233
234
				/* 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;

235
236
			}

Thomas White's avatar
Thomas White committed
237
			/* Measure intensities if requested */
238
			if ( config_nearbragg ) {
239
				output_intensities(&image, image.indexed_cell);
240
			}
Thomas White's avatar
Thomas White committed
241

Thomas White's avatar
Thomas White committed
242
			/* Simulate pattern if requested */
243
244
			if ( config_simulate ) {

245
246
				image.data = NULL;

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

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

			}

Thomas White's avatar
Thomas White committed
269
270
		}

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

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

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