facetron.c 6.86 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
/*
 * facetron.c
 *
 * Profile fitting for coherent nanocrystallography
 *
 * (c) 2006-2010 Thomas White <taw@physics.org>
 *
 * 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>
23
#include <errno.h>
24
25
26

#include "image.h"
#include "cell.h"
27
28
#include "hdf5-file.h"
#include "detector.h"
29
#include "geometry.h"
30
31
32
33
34
35
36
37
38
39


static void show_help(const char *s)
{
	printf("Syntax: %s [options]\n\n", s);
	printf(
"Profile fitting for coherent nanocrystallography.\n"
"\n"
"  -h, --help                 Display this help message.\n"
"\n"
40
41
42
"  -i, --input=<filename>     Specify the name of the input stream.\n"
"                              Can be '-' for stdin.\n"
"  -g. --geometry=<file>   Get detector geometry from file.\n"
43
44
45
46
);
}


Thomas White's avatar
Thomas White committed
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
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
132
133
134
135
136
137
138
139
140
141
142
143
144
static UnitCell *read_orientation_matrix(FILE *fh)
{
	float u, v, w;
	struct rvec as, bs, cs;
	UnitCell *cell;
	char line[1024];

	if ( fgets(line, 1023, fh) == NULL ) return NULL;
	if ( sscanf(line, "astar = %f %f %f", &u, &v, &w) != 3 ) {
		ERROR("Couldn't read a-star\n");
		return NULL;
	}
	as.u = u*1e9;  as.v = v*1e9;  as.w = w*1e9;
	if ( fgets(line, 1023, fh) == NULL ) return NULL;
	if ( sscanf(line, "bstar = %f %f %f", &u, &v, &w) != 3 ) {
		ERROR("Couldn't read b-star\n");
		return NULL;
	}
	bs.u = u*1e9;  bs.v = v*1e9;  bs.w = w*1e9;
	if ( fgets(line, 1023, fh) == NULL ) return NULL;
	if ( sscanf(line, "cstar = %f %f %f", &u, &v, &w) != 3 ) {
		ERROR("Couldn't read c-star\n");
		return NULL;
	}
	cs.u = u*1e9;  cs.v = v*1e9;  cs.w = w*1e9;
	cell = cell_new_from_axes(as, bs, cs);

	return cell;
}


static int find_chunk(FILE *fh, UnitCell **cell, char **filename)
{
	char line[1024];
	char *rval = NULL;

	do {

		rval = fgets(line, 1023, fh);
		if ( rval == NULL ) continue;

		chomp(line);

		if ( strncmp(line, "Reflections from indexing", 25) != 0 ) {
			continue;
		}

		*filename = strdup(line+29);

		/* Skip two lines (while checking for errors) */
		rval = fgets(line, 1023, fh);
		if ( rval == NULL ) continue;
		rval = fgets(line, 1023, fh);
		if ( rval == NULL ) continue;

		*cell = read_orientation_matrix(fh);
		if ( *cell == NULL ) {
			STATUS("Got filename but no cell for %s\n", *filename);
			continue;
		}

		return 0;

	} while ( rval != NULL );

	return 1;
}


static void get_trial_cell(UnitCell *out, UnitCell *in, int i, double step)
{
	double asx, asy, asz;
	double bsx, bsy, bsz;
	double csx, csy, csz;

	cell_get_reciprocal(in, &asx, &asy, &asz, &bsx, &bsy, &bsz,
	                    &csx, &csy, &csz);

	switch ( i ) {
	case  0 : asx += step; break;
	case  1 : asx -= step; break;
	case  2 : asy += step; break;
	case  3 : asy -= step; break;
	case  4 : asz += step; break;
	case  5 : asz -= step; break;
	case  6 : bsx += step; break;
	case  7 : bsx -= step; break;
	case  8 : bsy += step; break;
	case  9 : bsy -= step; break;
	case 10 : bsz += step; break;
	case 11 : bsz -= step; break;
	case 12 : csx += step; break;
	case 13 : csx -= step; break;
	case 14 : csy += step; break;
	case 15 : csy -= step; break;
	case 16 : csz += step; break;
	case 17 : csz -= step; break;
	default : break;
145
146
	}

Thomas White's avatar
Thomas White committed
147
	cell_set_reciprocal(out, asx, asy, asz, bsx, bsy, bsz, csx, csy, csz);
148
149
150
}


Thomas White's avatar
Thomas White committed
151
152
static void try_refine(struct image *image, UnitCell *cell,
                       double da, double dw, double step)
153
{
Thomas White's avatar
Thomas White committed
154
155
156
157
158
159
160
161
162
163
164
	struct reflhit *hits;
	int np;
	double itot;
	UnitCell *trial_cell;
	int did_something;

	trial_cell = cell_new();

	hits = find_intersections(image, cell, da, dw, &np, 0);
	itot = integrate_all(image, hits, np);
	STATUS("%f\n", itot);
165
166
167

	do {

Thomas White's avatar
Thomas White committed
168
		int i;
169

Thomas White's avatar
Thomas White committed
170
		did_something = 0;
171

Thomas White's avatar
Thomas White committed
172
		for ( i=0; i<18; i++ ) {
173

Thomas White's avatar
Thomas White committed
174
175
			struct reflhit *trial_hits;
			double trial_itot;
176

Thomas White's avatar
Thomas White committed
177
			get_trial_cell(trial_cell, cell, i, step);
178

Thomas White's avatar
Thomas White committed
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
			trial_hits = find_intersections(image, trial_cell,
				                        da, dw, &np, 0);
			trial_itot = integrate_all(image, hits, np);

			if ( trial_itot > itot ) {

				double asx, asy, asz;
				double bsx, bsy, bsz;
				double csx, csy, csz;

				cell_get_reciprocal(trial_cell, &asx, &asy,
				                    &asz, &bsx, &bsy, &bsz,
				                    &csx, &csy, &csz);
				cell_set_reciprocal(cell, asx, asy, asz, bsx,
				                    bsy, bsz, csx, csy, csz);

				itot = trial_itot;
				free(hits);
				hits = trial_hits;

				did_something = 1;

			} else {

				free(trial_hits);

			}
206
207
		}

Thomas White's avatar
Thomas White committed
208
	} while ( did_something );
209

Thomas White's avatar
Thomas White committed
210
211
	free(trial_cell);
}
212

Thomas White's avatar
Thomas White committed
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231

/* Predict peaks */
static void pre_refine(struct image *image, UnitCell *cell,
                       double *da, double *dw)
{
	double asx, asy, asz;
	double bsx, bsy, bsz;
	double csx, csy, csz;
	double istep, step;

	/* Start by changing parameters by 1% */
	cell_get_reciprocal(cell, &asx, &asy,
	                    &asz, &bsx, &bsy, &bsz,
	                    &csx, &csy, &csz);
	istep = (asx+asy+asz+bsx+bsy+bsz+csx+csy+csz)/900.0;

	for ( step=istep; step>istep/100.0; step-=istep/100.0 ) {
		try_refine(image, cell, *da, *dw, step);
	}
232
233
234
}


235
236
237
238
int main(int argc, char *argv[])
{
	int c;
	char *infile = NULL;
239
	FILE *fh;
240
	UnitCell *cell;
241
242
243
	char *filename;
	struct detector *det;
	char *geometry = NULL;
244
245
246
247
248

	/* Long options */
	const struct option longopts[] = {
		{"help",               0, NULL,               'h'},
		{"input",              1, NULL,               'i'},
249
		{"geometry",           1, NULL,               'g'},
250
251
252
253
		{0, 0, NULL, 0}
	};

	/* Short options */
254
	while ((c = getopt_long(argc, argv, "hi:g:", longopts, NULL)) != -1) {
255
256
257
258
259
260
261
262
263
264

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

		case 'i' :
			infile = strdup(optarg);
			break;

265
266
		case 'g' :
			geometry = strdup(optarg);
267
268
269
270
271
272
273
274
275
276
277
			break;

		case 0 :
			break;

		default :
			return 1;
		}

	}

278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
	if ( infile == NULL ) infile = strdup("-");
	if ( strcmp(infile, "-") == 0 ) {
		fh = stdin;
	} else {
		fh = fopen(infile, "r");
	}
	free(infile);
	if ( fh == NULL ) {
		ERROR("Couldn't open input stream '%s'\n", infile);
		return ENOENT;
	}

	det = get_detector_geometry(geometry);
	if ( det == NULL ) {
		ERROR("Failed to read detector geometry from '%s'\n", geometry);
293
294
		return 1;
	}
295
296
297
298
299
300
301
	free(geometry);

	/* Loop over all successfully indexed patterns */
	while ( find_chunk(fh, &cell, &filename) == 0 ) {

		struct image image;
		struct hdfile *hdfile;
Thomas White's avatar
Thomas White committed
302
303
		double da, dw;
		int np;
304
305
306
307
308
309
310
311
312
313
314
315
316
317

		STATUS("Integrating intensities from '%s'\n", filename);

		image.det = det;

		hdfile = hdfile_open(filename);
		if ( hdfile == NULL ) {
			return ENOENT;
		} else if ( hdfile_set_image(hdfile, "/data/data0") ) {
			ERROR("Couldn't select path\n");
			return ENOENT;
		}

		hdf5_read(hdfile, &image, 0);
Thomas White's avatar
Thomas White committed
318
319
320
321
322
323
324

		da = 5.0e-3;     /* Initial divergence */
		dw = 3.0/100.0;  /* Initial bandwidth */

		pre_refine(&image, cell, &da, &dw);

		find_intersections(&image, cell, da, dw, &np, 1);
325
326

		hdfile_close(hdfile);
327
		cell_free(cell);
328
329
330
331
332
		free(filename);
		free(image.data);
		free(image.flags);

	}
333

334
335
336
	free(det->panels);
	free(det);
	fclose(fh);
337
338
339

	return 0;
}