indexamajig.c 4.65 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
"      --no-index          Do everything else (including fine peak search and\n"
"                           writing 'xfel.drx' if DirAx is being used), but\n"
"                           don't actually index.\n"
"      --dirax             Use DirAx for indexing.\n"
Thomas White's avatar
Thomas White committed
50
"      --dump-peaks        Write the results of the peak search to stdout.\n"
51
52
53
"      --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
54
55
56
57
58
59
60
61
62
63
64
65
"\n");
}


int main(int argc, char *argv[])
{
	int c;
	char *filename = NULL;
	FILE *fh;
	char *rval;
	int n_images;
	int n_hits;
66
	int config_noindex = 0;
Thomas White's avatar
Thomas White committed
67
	int config_dumpfound = 0;
Thomas White's avatar
Thomas White committed
68
	int config_dirax = 0;
69
70
	int config_nearbragg = 0;
	int config_simulate = 0;
Thomas White's avatar
Thomas White committed
71
72
73
74
75

	/* Long options */
	const struct option longopts[] = {
		{"help",               0, NULL,               'h'},
		{"input",              1, NULL,               'i'},
76
		{"no-index",           0, &config_noindex,     1},
Thomas White's avatar
Thomas White committed
77
		{"dump-peaks",         0, &config_dumpfound,   1},
78
		{"near-bragg",         0, &config_nearbragg,   1},
Thomas White's avatar
Thomas White committed
79
		{"dirax",              0, &config_dirax,       1},
80
		{"simulate",           0, &config_simulate,    1},
Thomas White's avatar
Thomas White committed
81
82
83
84
		{0, 0, NULL, 0}
	};

	/* Short options */
85
	while ((c = getopt_long(argc, argv, "hi:w", longopts, NULL)) != -1) {
Thomas White's avatar
Thomas White committed
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131

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

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

		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;
	}

	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
132
		if ( rval == NULL ) continue;
Thomas White's avatar
Thomas White committed
133
134
		chomp(line);

Thomas White's avatar
Thomas White committed
135
		image.features = NULL;
136
		image.molecule = load_molecule();
Thomas White's avatar
Thomas White committed
137
		image.data = NULL;
138
		image.indexed_cell = NULL;
Thomas White's avatar
Thomas White committed
139

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

142
143
		n_images++;

Thomas White's avatar
Thomas White committed
144
145
		hdfile = hdfile_open(line);
		if ( hdfile == NULL ) {
146
			continue;
Thomas White's avatar
Thomas White committed
147
148
		} else if ( hdfile_set_first_image(hdfile, "/") ) {
			ERROR("Couldn't select path\n");
149
			continue;
Thomas White's avatar
Thomas White committed
150
151
152
153
		}

		hdf5_read(hdfile, &image);

154
		/* Perform 'fine' peak search */
155
		search_peaks(&image);
Thomas White's avatar
Thomas White committed
156

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

159
			n_hits++;
Thomas White's avatar
Thomas White committed
160

161
162
			if ( config_dumpfound ) dump_peaks(&image);

163
164
			/* Not indexing?  Then there's nothing left to do. */
			if ( config_noindex ) goto done;
Thomas White's avatar
Thomas White committed
165

166
			/* Calculate orientation matrix (by magic) */
167
168
			index_pattern(&image, config_noindex, config_dirax);
			if ( image.indexed_cell == NULL ) goto done;
169
170

			if ( config_nearbragg || config_simulate ) {
171

Thomas White's avatar
Thomas White committed
172
173
174
175
176
177
				/* Simulate a diffraction pattern */
				image.sfacs = NULL;
				image.qvecs = NULL;
				image.twotheta = NULL;
				image.hdr = NULL;

178
179
180
181
182
183
184
				/* 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);

185
186
187
			}

			if ( config_nearbragg ) {
188
				/* Read h,k,l,I */
189
				output_intensities(&image, image.indexed_cell);
190
			}
Thomas White's avatar
Thomas White committed
191

192
193
			if ( config_simulate ) {

194
195
				image.data = NULL;

196
197
198
199
200
201
202
203
204
205
206
207
				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
208
209
		}

Thomas White's avatar
Thomas White committed
210
done:
Thomas White's avatar
Thomas White committed
211
		free(image.data);
212
		image_feature_list_free(image.features);
Thomas White's avatar
Thomas White committed
213
		hdfile_close(hdfile);
214
		H5close();
Thomas White's avatar
Thomas White committed
215

Thomas White's avatar
Thomas White committed
216
217
218
219
220
221
222
223
224
	} while ( rval != NULL );

	fclose(fh);

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

	return 0;
}