indexamajig.c 5.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
29
#include "intensities.h"
#include "ewald.h"
30
#include "peaks.h"
31
32
#include "diffraction.h"
#include "detector.h"
33
#include "sfac.h"
Thomas White's avatar
Thomas White committed
34
35
36
37
38
39
40
41
42
43
44
45


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
46
47
48
49
50
51
"      --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
52
"      --dump-peaks        Write the results of the peak search to stdout.\n"
53
54
55
"      --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
56
57
58
59
60
61
62
63
64
65
66
67
"\n");
}


int main(int argc, char *argv[])
{
	int c;
	char *filename = NULL;
	FILE *fh;
	char *rval;
	int n_images;
	int n_hits;
68
	int config_noindex = 0;
Thomas White's avatar
Thomas White committed
69
	int config_dumpfound = 0;
70
	int config_nearbragg = 0;
Thomas White's avatar
Thomas White committed
71
	int config_writedrx = 0;
72
	int config_simulate = 0;
Thomas White's avatar
Thomas White committed
73
74
	IndexingMethod indm;
	char *indm_str = NULL;
Thomas White's avatar
Thomas White committed
75
76
77
78
79

	/* Long options */
	const struct option longopts[] = {
		{"help",               0, NULL,               'h'},
		{"input",              1, NULL,               'i'},
80
		{"no-index",           0, &config_noindex,     1},
Thomas White's avatar
Thomas White committed
81
		{"dump-peaks",         0, &config_dumpfound,   1},
82
		{"near-bragg",         0, &config_nearbragg,   1},
Thomas White's avatar
Thomas White committed
83
84
		{"write-drx",          0, &config_writedrx,    1},
		{"indexing",           1, NULL,               'z'},
85
		{"simulate",           0, &config_simulate,    1},
Thomas White's avatar
Thomas White committed
86
87
88
89
		{0, 0, NULL, 0}
	};

	/* Short options */
90
	while ((c = getopt_long(argc, argv, "hi:w", longopts, NULL)) != -1) {
Thomas White's avatar
Thomas White committed
91
92
93
94
95
96
97
98
99
100
101
102

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

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

Thomas White's avatar
Thomas White committed
103
104
105
106
107
		case 'z' : {
			indm_str = strdup(optarg);
			break;
		}

Thomas White's avatar
Thomas White committed
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
		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
133
134
135
136
137
138
139
140
141
142
143
144
145
146
	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
147
148
149
150
151
152
153
154
155
	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
156
		if ( rval == NULL ) continue;
Thomas White's avatar
Thomas White committed
157
158
		chomp(line);

Thomas White's avatar
Thomas White committed
159
		image.features = NULL;
160
		image.molecule = load_molecule();
Thomas White's avatar
Thomas White committed
161
		image.data = NULL;
162
		image.indexed_cell = NULL;
Thomas White's avatar
Thomas White committed
163

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

166
167
		n_images++;

Thomas White's avatar
Thomas White committed
168
169
		hdfile = hdfile_open(line);
		if ( hdfile == NULL ) {
170
			continue;
Thomas White's avatar
Thomas White committed
171
172
		} else if ( hdfile_set_first_image(hdfile, "/") ) {
			ERROR("Couldn't select path\n");
173
			continue;
Thomas White's avatar
Thomas White committed
174
175
176
177
		}

		hdf5_read(hdfile, &image);

178
		/* Perform 'fine' peak search */
179
		search_peaks(&image);
Thomas White's avatar
Thomas White committed
180

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

183
184
			if ( config_dumpfound ) dump_peaks(&image);

Thomas White's avatar
Thomas White committed
185
186
187
188
189
			/* 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
190

191
			/* Calculate orientation matrix (by magic) */
Thomas White's avatar
Thomas White committed
192
193
194
195
196
			if ( config_writedrx || (indm != INDEXING_NONE) ) {
				index_pattern(&image, indm);
			}

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

199
200
			n_hits++;

Thomas White's avatar
Thomas White committed
201
202
			/* Simulation or intensity measurements both require
			 * Ewald sphere vectors */
203
			if ( config_nearbragg || config_simulate ) {
204

Thomas White's avatar
Thomas White committed
205
206
207
208
209
210
				/* Simulate a diffraction pattern */
				image.sfacs = NULL;
				image.qvecs = NULL;
				image.twotheta = NULL;
				image.hdr = NULL;

211
212
213
214
215
216
217
				/* 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;
				get_ewald(&image);

218
219
			}

Thomas White's avatar
Thomas White committed
220
			/* Measure intensities if requested */
221
			if ( config_nearbragg ) {
222
				output_intensities(&image, image.indexed_cell);
223
			}
Thomas White's avatar
Thomas White committed
224

Thomas White's avatar
Thomas White committed
225
			/* Simulate pattern if requested */
226
227
			if ( config_simulate ) {

228
229
				image.data = NULL;

230
231
232
233
234
235
236
237
238
239
240
241
				get_diffraction(&image, 8, 8, 8);
				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,
				           image.width, image.height);

			}

Thomas White's avatar
Thomas White committed
242
243
		}

Thomas White's avatar
Thomas White committed
244
done:
Thomas White's avatar
Thomas White committed
245
		free(image.data);
246
		image_feature_list_free(image.features);
Thomas White's avatar
Thomas White committed
247
		hdfile_close(hdfile);
248
		H5close();
Thomas White's avatar
Thomas White committed
249

Thomas White's avatar
Thomas White committed
250
251
252
253
254
255
256
257
258
	} while ( rval != NULL );

	fclose(fh);

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

	return 0;
}