facetron.c 9.76 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>
Thomas White's avatar
Thomas White committed
23
#include <assert.h>
24

Thomas White's avatar
Thomas White committed
25
#include "utils.h"
26
#include "hdf5-file.h"
Thomas White's avatar
Thomas White committed
27
#include "symmetry.h"
Thomas White's avatar
Thomas White committed
28
29
#include "reflections.h"
#include "stream.h"
30
#include "geometry.h"
31
#include "peaks.h"
Thomas White's avatar
Thomas White committed
32
#include "thread-pool.h"
33
34
35
36
37
38


static void show_help(const char *s)
{
	printf("Syntax: %s [options]\n\n", s);
	printf(
Thomas White's avatar
Thomas White committed
39
"Post-refinement and profile fitting for coherent nanocrystallography.\n"
40
41
42
"\n"
"  -h, --help                 Display this help message.\n"
"\n"
Thomas White's avatar
Thomas White committed
43
44
"  -i, --input=<filename>     Specify the name of the input 'stream'.\n"
"                              (must be a file, not e.g. stdin)\n"
Thomas White's avatar
Thomas White committed
45
"  -o, --output=<filename>    Output filename.  Default: facetron.hkl.\n"
Thomas White's avatar
Thomas White committed
46
47
48
49
"  -g. --geometry=<file>      Get detector geometry from file.\n"
"  -x, --prefix=<p>           Prefix filenames from input file with <p>.\n"
"      --basename             Remove the directory parts of the filenames.\n"
"      --no-check-prefix      Don't attempt to correct the --prefix.\n"
50
51
52
"  -y, --symmetry=<sym>       Merge according to symmetry <sym>.\n"
"  -n, --iterations=<n>       Run <n> cycles of post-refinement.\n"
"\n"
Thomas White's avatar
Thomas White committed
53
54
55
56
"  -j <n>                     Run <n> analyses in parallel.\n");
}


Thomas White's avatar
Thomas White committed
57
struct refine_args
Thomas White's avatar
Thomas White committed
58
{
Thomas White's avatar
Thomas White committed
59
60
61
62
63
64
65
66
67
68
69
	const char *sym;
	ReflItemList *obs;
	double *i_full;
	struct image *image;
};


static void refine_image(int mytask, void *tasks)
{
	struct refine_args *all_args = tasks;
	struct refine_args *pargs = &all_args[mytask];
70
71
72
73
	/* Do, er, something. */
}


Thomas White's avatar
Thomas White committed
74
struct integrate_args
75
{
Thomas White's avatar
Thomas White committed
76
77
78
79
80
81
82
83
84
85
86
87
88
	const char *sym;
	ReflItemList *obs;
	double *i_full;
	unsigned int *cts;
	pthread_mutex_t *list_lock;
	struct image *image;
};


static void integrate_image(int mytask, void *tasks)
{
	struct integrate_args *all_args = tasks;
	struct integrate_args *pargs = &all_args[mytask];
89
90
91
92
	struct reflhit *spots;
	int j, n;
	struct hdfile *hdfile;
	struct image *image = pargs->image;
93

94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
	hdfile = hdfile_open(image->filename);
	if ( hdfile == NULL ) {
		ERROR("Couldn't open '%s'\n", image->filename);
		return;
	} else if ( hdfile_set_image(hdfile, "/data/data0") ) {
		ERROR("Couldn't select path\n");
		hdfile_close(hdfile);
		return;
	}

	if ( hdf5_read(hdfile, pargs->image, 0) ) {
		ERROR("Couldn't read '%s'\n", image->filename);
		hdfile_close(hdfile);
		return;
	}

Thomas White's avatar
Thomas White committed
110
	/* Figure out which spots should appear in this pattern */
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
145
146
147
148
	spots = find_intersections(image, image->indexed_cell,
	                           image->div, image->bw, &n, 0);

	/* For each reflection, estimate the partiality */
	for ( j=0; j<n; j++ ) {

		signed int h, k, l;
		float i_partial;
		double p;
		float xc, yc;

		h = spots[j].h;
		k = spots[j].k;
		l = spots[j].l;

		/* Calculated partiality of this spot in this pattern */
		p = partiality(image, h, k, l);

		/* Don't attempt to use spots with very small
		 * partialities, since it won't be accurate. */
		if ( p < 0.1 ) continue;

		/* Actual measurement of this reflection from this
		 * pattern? */
		/* FIXME: Coordinates aren't whole numbers */
		if ( integrate_peak(image, spots[j].x, spots[j].y,
		                    &xc, &yc, &i_partial, 1, 1) ) continue;

		pthread_mutex_lock(pargs->list_lock);
		integrate_intensity(pargs->i_full, h, k, l, i_partial);
		integrate_count(pargs->cts, h, k, l, 1);
		if ( !find_item(pargs->obs, h, k, l) ) {
			add_item(pargs->obs, h, k, l);
		}
		pthread_mutex_unlock(pargs->list_lock);

	}

149
150
	free(image->data);
	if ( image->flags != NULL ) free(image->flags);
151
	hdfile_close(hdfile);
Thomas White's avatar
Thomas White committed
152
153
154
155
156
	free(spots);

	/* Muppet proofing */
	image->data = NULL;
	image->flags = NULL;
Thomas White's avatar
Thomas White committed
157
158
159
}


Thomas White's avatar
Thomas White committed
160
161
162
static void refine_all(struct image *images, int n_total_patterns,
                       struct detector *det, const char *sym,
                       ReflItemList *obs, double *i_full, int nthreads)
Thomas White's avatar
Thomas White committed
163
{
Thomas White's avatar
Thomas White committed
164
	struct refine_args *tasks;
Thomas White's avatar
Thomas White committed
165
	int i;
166

Thomas White's avatar
Thomas White committed
167
168
	tasks = malloc(n_total_patterns * sizeof(struct refine_args));
	for ( i=0; i<n_total_patterns; i++ ) {
Thomas White's avatar
Thomas White committed
169

Thomas White's avatar
Thomas White committed
170
171
172
173
		tasks[i].sym = sym;
		tasks[i].obs = obs;
		tasks[i].i_full = i_full;
		tasks[i].image = &images[i];
Thomas White's avatar
Thomas White committed
174
175
176

	}

177
178
	run_thread_range(n_total_patterns, nthreads, "Refining",
	                 refine_image, tasks);
Thomas White's avatar
Thomas White committed
179

Thomas White's avatar
Thomas White committed
180
	free(tasks);
181
182
183
184
185
}


static void estimate_full(struct image *images, int n_total_patterns,
                          struct detector *det, const char *sym,
186
                          ReflItemList *obs, double *i_full, int nthreads)
187
188
{
	int i;
189
	unsigned int *cts;
Thomas White's avatar
Thomas White committed
190
191
	struct integrate_args *tasks;
	pthread_mutex_t list_lock = PTHREAD_MUTEX_INITIALIZER;
192

193
194
	cts = new_list_count();
	clear_items(obs);
195

Thomas White's avatar
Thomas White committed
196
197
198
199
200
201
202
203
204
205
206
207
	tasks = malloc(n_total_patterns * sizeof(struct integrate_args));
	for ( i=0; i<n_total_patterns; i++ ) {

		tasks[i].sym = sym;
		tasks[i].obs = obs;
		tasks[i].i_full = i_full;
		tasks[i].cts = cts;
		tasks[i].list_lock = &list_lock;
		tasks[i].image = &images[i];

	}

208
209
	run_thread_range(n_total_patterns, nthreads, "Integrating",
	                 integrate_image, tasks);
Thomas White's avatar
Thomas White committed
210
211

	free(tasks);
212

213
214
	/* Divide the totals to get the means */
	for ( i=0; i<num_items(obs); i++ ) {
215

216
217
		struct refl_item *it;
		double total;
218

219
220
221
222
		it = get_item(obs, i);
		total = lookup_intensity(i_full, it->h, it->k, it->l);
		total /= lookup_count(cts, it->h, it->k, it->l);
		set_intensity(i_full, it->h, it->k, it->l, total);
223

224
	}
225
226

	free(cts);
Thomas White's avatar
Thomas White committed
227
}
228

Thomas White's avatar
Thomas White committed
229

Thomas White's avatar
Thomas White committed
230
231
232
233
int main(int argc, char *argv[])
{
	int c;
	char *infile = NULL;
Thomas White's avatar
Thomas White committed
234
	char *outfile = NULL;
Thomas White's avatar
Thomas White committed
235
236
	char *geomfile = NULL;
	char *prefix = NULL;
Thomas White's avatar
Thomas White committed
237
238
	char *sym = NULL;
	FILE *fh;
Thomas White's avatar
Thomas White committed
239
240
241
242
	int nthreads = 1;
	int config_basename = 0;
	int config_checkprefix = 1;
	struct detector *det;
Thomas White's avatar
Thomas White committed
243
244
245
246
	double *i_full;
	ReflItemList *obs;
	int i;
	int n_total_patterns;
247
248
	struct image *images;
	int n_iter = 10;
Thomas White's avatar
Thomas White committed
249
250
251
252
253

	/* Long options */
	const struct option longopts[] = {
		{"help",               0, NULL,               'h'},
		{"input",              1, NULL,               'i'},
Thomas White's avatar
Thomas White committed
254
		{"output",             1, NULL,               'o'},
Thomas White's avatar
Thomas White committed
255
256
257
258
		{"geometry",           1, NULL,               'g'},
		{"prefix",             1, NULL,               'x'},
		{"basename",           0, &config_basename,    1},
		{"no-check-prefix",    0, &config_checkprefix, 0},
259
		{"symmetry",           1, NULL,               'y'},
260
		{"iterations",         1, NULL,               'n'},
Thomas White's avatar
Thomas White committed
261
262
263
264
		{0, 0, NULL, 0}
	};

	/* Short options */
Thomas White's avatar
Thomas White committed
265
	while ((c = getopt_long(argc, argv, "hi:g:x:j:y:o:",
266
267
	                        longopts, NULL)) != -1)
	{
Thomas White's avatar
Thomas White committed
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289

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

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

		case 'g' :
			geomfile = strdup(optarg);
			break;

		case 'x' :
			prefix = strdup(optarg);
			break;

		case 'j' :
			nthreads = atoi(optarg);
			break;

290
291
292
293
		case 'y' :
			sym = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
294
295
296
297
		case 'o' :
			outfile = strdup(optarg);
			break;

298
299
300
301
		case 'n' :
			n_iter = atoi(optarg);
			break;

Thomas White's avatar
Thomas White committed
302
303
304
305
306
307
308
309
310
		case 0 :
			break;

		default :
			return 1;
		}

	}

Thomas White's avatar
Thomas White committed
311
	/* Sanitise input filename and open */
Thomas White's avatar
Thomas White committed
312
313
314
315
316
317
318
319
320
321
322
323
324
325
	if ( infile == NULL ) {
		infile = strdup("-");
	}
	if ( strcmp(infile, "-") == 0 ) {
		fh = stdin;
	} else {
		fh = fopen(infile, "r");
	}
	if ( fh == NULL ) {
		ERROR("Failed to open input file '%s'\n", infile);
		return 1;
	}
	free(infile);

Thomas White's avatar
Thomas White committed
326
327
328
329
330
331
	/* Sanitise output filename */
	if ( outfile == NULL ) {
		outfile = strdup("facetron.hkl");
	}

	/* Sanitise prefix */
Thomas White's avatar
Thomas White committed
332
333
334
335
336
337
338
339
	if ( prefix == NULL ) {
		prefix = strdup("");
	} else {
		if ( config_checkprefix ) {
			prefix = check_prefix(prefix);
		}
	}

Thomas White's avatar
Thomas White committed
340
	/* Get detector geometry */
Thomas White's avatar
Thomas White committed
341
342
343
344
345
346
347
	det = get_detector_geometry(geomfile);
	if ( det == NULL ) {
		ERROR("Failed to read detector geometry from '%s'\n", geomfile);
		return 1;
	}
	free(geomfile);

Thomas White's avatar
Thomas White committed
348
349
350
351
352
353
354
	/* Prepare for iteration */
	i_full = new_list_intensity();
	obs = new_items();

	n_total_patterns = count_patterns(fh);
	STATUS("There are %i patterns to process\n", n_total_patterns);

355
356
357
358
359
360
361
362
363
364
365
366
	images = malloc(n_total_patterns * sizeof(struct image));
	if ( images == NULL ) {
		ERROR("Couldn't allocate memory for images.\n");
		return 1;
	}

	/* Fill in what we know about the images so far */
	rewind(fh);
	for ( i=0; i<n_total_patterns; i++ ) {

		UnitCell *cell;
		char *filename;
367
		char *fnamereal;
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383

		if ( find_chunk(fh, &cell, &filename) == 1 ) {
			ERROR("Couldn't get all of the filenames and cells"
			      " from the input stream.\n");
			return 1;
		}

		images[i].indexed_cell = cell;

		/* Mangle the filename now */
		if ( config_basename ) {
			char *tmp;
			tmp = strdup(basename(filename));
			free(filename);
			filename = tmp;
		}
384
		fnamereal = malloc(1024);
385
386
387
388
389
		snprintf(fnamereal, 1023, "%s%s", prefix, filename);

		images[i].filename = fnamereal;
		images[i].div = 0.5e-3;
		images[i].bw = 0.001;
390
391
392
393
394
		images[i].orientation.w = 1.0;
		images[i].orientation.x = 0.0;
		images[i].orientation.y = 0.0;
		images[i].orientation.z = 0.0;
		images[i].det = det;
395

Thomas White's avatar
Thomas White committed
396
397
398
399
		/* Muppet proofing */
		images[i].data = NULL;
		images[i].flags = NULL;

400
401
402
403
404
405
406
407
408
		free(filename);

		progress_bar(i, n_total_patterns-1, "Loading pattern data");

	}
	fclose(fh);
	free(prefix);

	/* Make initial estimates */
409
410
	estimate_full(images, n_total_patterns, det, sym, obs, i_full,
	              nthreads);
411

Thomas White's avatar
Thomas White committed
412
	/* Iterate */
413
	for ( i=0; i<n_iter; i++ ) {
Thomas White's avatar
Thomas White committed
414

415
		STATUS("Post refinement iteration %i of %i\n", i+1, n_iter);
Thomas White's avatar
Thomas White committed
416
417

		/* Refine the geometry of all patterns to get the best fit */
418
419
		refine_all(images, n_total_patterns, det, sym, obs, i_full,
		           nthreads);
Thomas White's avatar
Thomas White committed
420
421

		/* Re-estimate all the full intensities */
422
423
		estimate_full(images, n_total_patterns, det, sym, obs, i_full,
		              nthreads);
Thomas White's avatar
Thomas White committed
424
425
426
427
428

	}

	/* Output results */
	write_reflections(outfile, obs, i_full, NULL, NULL, NULL);
Thomas White's avatar
Thomas White committed
429

Thomas White's avatar
Thomas White committed
430
431
432
	/* Clean up */
	free(i_full);
	delete_items(obs);
Thomas White's avatar
Thomas White committed
433
	free(sym);
Thomas White's avatar
Thomas White committed
434
	free(outfile);
435
436
437
438
439
440
441
	free(det->panels);
	free(det);
	for ( i=0; i<n_total_patterns; i++ ) {
		cell_free(images[i].indexed_cell);
		free(images[i].filename);
	}
	free(images);
442
443
444

	return 0;
}