facetron.c 13.1 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
24
25
#include <pthread.h>
#include <sys/time.h>
#include <assert.h>
26

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


#define MAX_THREADS (256)

struct process_args
{
40
	struct image *image;
Thomas White's avatar
Thomas White committed
41
42
43
44
45
46
47

	/* Thread control */
	pthread_mutex_t control_mutex;  /* Protects the scary stuff below */
	int start;
	int finish;
	int done;

48
49
50
	/* Analysis routine */
	void (*func)(struct process_args *);

51
	/* Analysis parameters */
Thomas White's avatar
Thomas White committed
52
	const char *sym;
53
	pthread_mutex_t *list_lock;  /* Protects 'obs', 'i_full' and 'cts' */
54
55
	ReflItemList *obs;
	double *i_full;
56
	unsigned int *cts;
Thomas White's avatar
Thomas White committed
57
};
58
59
60
61
62
63


static void show_help(const char *s)
{
	printf("Syntax: %s [options]\n\n", s);
	printf(
Thomas White's avatar
Thomas White committed
64
"Post-refinement and profile fitting for coherent nanocrystallography.\n"
65
66
67
"\n"
"  -h, --help                 Display this help message.\n"
"\n"
Thomas White's avatar
Thomas White committed
68
69
"  -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
70
"  -o, --output=<filename>    Output filename.  Default: facetron.hkl.\n"
Thomas White's avatar
Thomas White committed
71
72
73
74
"  -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"
75
76
77
"  -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
78
79
80
81
"  -j <n>                     Run <n> analyses in parallel.\n");
}


82
static void refine_image(struct process_args *pargs)
Thomas White's avatar
Thomas White committed
83
{
84
85
86
87
88
89
90
91
92
93
	/* Do, er, something. */
}


static void integrate_image(struct process_args *pargs)
{
	struct reflhit *spots;
	int j, n;
	struct hdfile *hdfile;
	struct image *image = pargs->image;
94

95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
	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
111
	/* Figure out which spots should appear in this pattern */
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
149
	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);

	}

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

	/* Muppet proofing */
	image->data = NULL;
	image->flags = NULL;
Thomas White's avatar
Thomas White committed
158
159
160
161
162
163
164
165
166
167
168
169
}


static void *worker_thread(void *pargsv)
{
	struct process_args *pargs = pargsv;
	int finish;

	do {

		int wakeup;

Thomas White's avatar
Thomas White committed
170
171
172
173
174
		/* Acknowledge start */
		pthread_mutex_lock(&pargs->control_mutex);
		pargs->start = 0;
		pthread_mutex_unlock(&pargs->control_mutex);

175
		pargs->func(pargs);
Thomas White's avatar
Thomas White committed
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195

		pthread_mutex_lock(&pargs->control_mutex);
		pargs->done = 1;
		pthread_mutex_unlock(&pargs->control_mutex);

		/* Go to sleep until told to exit or process next image */
		do {

			pthread_mutex_lock(&pargs->control_mutex);
			/* Either of these can result in the thread waking up */
			wakeup = pargs->start || pargs->finish;
			finish = pargs->finish;
			pthread_mutex_unlock(&pargs->control_mutex);
			usleep(20000);

		} while ( !wakeup );

	} while ( !pargs->finish );

	return NULL;
196
197
198
}


199
200
201
202
203
static void munch_threads(struct image *images, int n_total_patterns,
                          struct detector *det, const char *sym,
                          ReflItemList *obs, double *i_full, unsigned int *cts,
                          int nthreads, void (*func)(struct process_args *),
                          const char *text)
Thomas White's avatar
Thomas White committed
204
{
Thomas White's avatar
Thomas White committed
205
206
	pthread_t workers[MAX_THREADS];
	struct process_args *worker_args[MAX_THREADS];
207
	pthread_mutex_t list_lock = PTHREAD_MUTEX_INITIALIZER;
Thomas White's avatar
Thomas White committed
208
209
	int worker_active[MAX_THREADS];
	int i;
Thomas White's avatar
Thomas White committed
210
	int n_done = 0;
211
	int n_started = 0;
212

213
	/* Initialise worker arguments with the unchanging data */
Thomas White's avatar
Thomas White committed
214
	for ( i=0; i<nthreads; i++ ) {
Thomas White's avatar
Thomas White committed
215

Thomas White's avatar
Thomas White committed
216
217
218
219
		worker_args[i] = malloc(sizeof(struct process_args));
		worker_active[i] = 0;
		pthread_mutex_init(&worker_args[i]->control_mutex, NULL);
		worker_args[i]->sym = sym;
220
221
222
223
224
		worker_args[i]->obs = obs;
		worker_args[i]->i_full = i_full;
		worker_args[i]->cts = cts;
		worker_args[i]->list_lock = &list_lock;
		worker_args[i]->func = func;
225

Thomas White's avatar
Thomas White committed
226
	}
227

Thomas White's avatar
Thomas White committed
228
229
	/* Start threads off */
	for ( i=0; i<nthreads; i++ ) {
Thomas White's avatar
Thomas White committed
230

Thomas White's avatar
Thomas White committed
231
232
		struct process_args *pargs;
		int r;
233

234
235
		if ( n_started == n_total_patterns ) break;

Thomas White's avatar
Thomas White committed
236
		pargs = worker_args[i];
237
		pargs->image = &images[n_started++];
Thomas White's avatar
Thomas White committed
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276

		pthread_mutex_lock(&pargs->control_mutex);
		pargs->done = 0;
		pargs->start = 1;
		pargs->finish = 0;
		pthread_mutex_unlock(&pargs->control_mutex);

		worker_active[i] = 1;
		r = pthread_create(&workers[i], NULL, worker_thread, pargs);
		if ( r != 0 ) {
			worker_active[i] = 0;
			ERROR("Couldn't start thread %i\n", i);
		}

	}

	/* Keep threads busy until the end of the data */
	do {

		int i;

		for ( i=0; i<nthreads; i++ ) {

			struct process_args *pargs;
			int done;

			/* Spend time working, not managing threads */
			usleep(100000);

			/* Are we using this thread record at all? */
			if ( !worker_active[i] ) continue;

			/* Has the thread finished yet? */
			pargs = worker_args[i];
			pthread_mutex_lock(&pargs->control_mutex);
			done = pargs->done;
			pthread_mutex_unlock(&pargs->control_mutex);
			if ( !done ) continue;

277
			/* Reset "done" flag */
Thomas White's avatar
Thomas White committed
278
			pthread_mutex_lock(&pargs->control_mutex);
279
			pargs->done = 0;
Thomas White's avatar
Thomas White committed
280
			pthread_mutex_unlock(&pargs->control_mutex);
281

Thomas White's avatar
Thomas White committed
282
			n_done++;
283
			progress_bar(n_done, n_total_patterns, text);
Thomas White's avatar
Thomas White committed
284

285
286
			/* If there are no more patterns, "done" will remain
			 * zero, so the last pattern will not be re-counted. */
287
			if ( n_started == n_total_patterns ) break;
Thomas White's avatar
Thomas White committed
288

289
290
			/* Start work on the next pattern */
			pargs->image = &images[n_started++];
Thomas White's avatar
Thomas White committed
291
292
293
294
295
296
			pthread_mutex_lock(&pargs->control_mutex);
			pargs->start = 1;
			pthread_mutex_unlock(&pargs->control_mutex);

		}

297
	} while ( n_started < n_total_patterns );
Thomas White's avatar
Thomas White committed
298
299
300
301

	/* Join threads */
	for ( i=0; i<nthreads; i++ ) {

302
		if ( !worker_active[i] ) continue;
Thomas White's avatar
Thomas White committed
303
304
305
306
307
308
309
310
311
312

		/* Tell the thread to exit */
		struct process_args *pargs = worker_args[i];
		pthread_mutex_lock(&pargs->control_mutex);
		pargs->finish = 1;
		pthread_mutex_unlock(&pargs->control_mutex);

		/* Wait for it to join */
		pthread_join(workers[i], NULL);

313
314
		if ( pargs->done ) {
			n_done++;
315
			progress_bar(n_done, n_total_patterns, text);
316
		} /* else this thread was not busy */
317
318

	}
319
320
321
322

	for ( i=0; i<nthreads; i++ ) {
		free(worker_args[i]);
	}
323
324
325
}


326
327
328
static void refine_all(struct image *images, int n_total_patterns,
                       struct detector *det, const char *sym,
                       ReflItemList *obs, double *i_full, int nthreads)
329
{
330
331
	munch_threads(images, n_total_patterns, det, sym, obs, i_full, NULL,
	              nthreads, refine_image, "Refining");
332
333
334
335
336
}


static void estimate_full(struct image *images, int n_total_patterns,
                          struct detector *det, const char *sym,
337
                          ReflItemList *obs, double *i_full, int nthreads)
338
339
{
	int i;
340
	unsigned int *cts;
341

342
343
	cts = new_list_count();
	clear_items(obs);
344

345
346
	munch_threads(images, n_total_patterns, det, sym, obs, i_full, cts,
	              nthreads, integrate_image, "Integrating");
347

348
349
	/* Divide the totals to get the means */
	for ( i=0; i<num_items(obs); i++ ) {
350

351
352
		struct refl_item *it;
		double total;
353

354
355
356
357
		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);
358

359
	}
360
361

	free(cts);
Thomas White's avatar
Thomas White committed
362
}
363

Thomas White's avatar
Thomas White committed
364

Thomas White's avatar
Thomas White committed
365
366
367
368
int main(int argc, char *argv[])
{
	int c;
	char *infile = NULL;
Thomas White's avatar
Thomas White committed
369
	char *outfile = NULL;
Thomas White's avatar
Thomas White committed
370
371
	char *geomfile = NULL;
	char *prefix = NULL;
Thomas White's avatar
Thomas White committed
372
373
	char *sym = NULL;
	FILE *fh;
Thomas White's avatar
Thomas White committed
374
375
376
377
	int nthreads = 1;
	int config_basename = 0;
	int config_checkprefix = 1;
	struct detector *det;
Thomas White's avatar
Thomas White committed
378
379
380
381
	double *i_full;
	ReflItemList *obs;
	int i;
	int n_total_patterns;
382
383
	struct image *images;
	int n_iter = 10;
Thomas White's avatar
Thomas White committed
384
385
386
387
388

	/* Long options */
	const struct option longopts[] = {
		{"help",               0, NULL,               'h'},
		{"input",              1, NULL,               'i'},
Thomas White's avatar
Thomas White committed
389
		{"output",             1, NULL,               'o'},
Thomas White's avatar
Thomas White committed
390
391
392
393
		{"geometry",           1, NULL,               'g'},
		{"prefix",             1, NULL,               'x'},
		{"basename",           0, &config_basename,    1},
		{"no-check-prefix",    0, &config_checkprefix, 0},
394
		{"symmetry",           1, NULL,               'y'},
395
		{"iterations",         1, NULL,               'n'},
Thomas White's avatar
Thomas White committed
396
397
398
399
		{0, 0, NULL, 0}
	};

	/* Short options */
Thomas White's avatar
Thomas White committed
400
	while ((c = getopt_long(argc, argv, "hi:g:x:j:y:o:",
401
402
	                        longopts, NULL)) != -1)
	{
Thomas White's avatar
Thomas White committed
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424

		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;

425
426
427
428
		case 'y' :
			sym = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
429
430
431
432
		case 'o' :
			outfile = strdup(optarg);
			break;

433
434
435
436
		case 'n' :
			n_iter = atoi(optarg);
			break;

Thomas White's avatar
Thomas White committed
437
438
439
440
441
442
443
444
445
		case 0 :
			break;

		default :
			return 1;
		}

	}

Thomas White's avatar
Thomas White committed
446
	/* Sanitise input filename and open */
Thomas White's avatar
Thomas White committed
447
448
449
450
451
452
453
454
455
456
457
458
459
460
	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
461
462
463
464
465
466
	/* Sanitise output filename */
	if ( outfile == NULL ) {
		outfile = strdup("facetron.hkl");
	}

	/* Sanitise prefix */
Thomas White's avatar
Thomas White committed
467
468
469
470
471
472
473
474
	if ( prefix == NULL ) {
		prefix = strdup("");
	} else {
		if ( config_checkprefix ) {
			prefix = check_prefix(prefix);
		}
	}

Thomas White's avatar
Thomas White committed
475
	/* Get detector geometry */
Thomas White's avatar
Thomas White committed
476
477
478
479
480
481
482
	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
483
484
485
486
487
488
489
	/* 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);

490
491
492
493
494
495
496
497
498
499
500
501
	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;
502
		char *fnamereal;
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518

		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;
		}
519
		fnamereal = malloc(1024);
520
521
522
523
524
		snprintf(fnamereal, 1023, "%s%s", prefix, filename);

		images[i].filename = fnamereal;
		images[i].div = 0.5e-3;
		images[i].bw = 0.001;
525
526
527
528
529
		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;
530

Thomas White's avatar
Thomas White committed
531
532
533
534
		/* Muppet proofing */
		images[i].data = NULL;
		images[i].flags = NULL;

535
536
537
538
539
540
541
542
543
		free(filename);

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

	}
	fclose(fh);
	free(prefix);

	/* Make initial estimates */
544
545
	estimate_full(images, n_total_patterns, det, sym, obs, i_full,
	              nthreads);
546

Thomas White's avatar
Thomas White committed
547
	/* Iterate */
548
	for ( i=0; i<n_iter; i++ ) {
Thomas White's avatar
Thomas White committed
549

550
		STATUS("Post refinement iteration %i of %i\n", i+1, n_iter);
Thomas White's avatar
Thomas White committed
551
552

		/* Refine the geometry of all patterns to get the best fit */
553
554
		refine_all(images, n_total_patterns, det, sym, obs, i_full,
		           nthreads);
Thomas White's avatar
Thomas White committed
555
556

		/* Re-estimate all the full intensities */
557
558
		estimate_full(images, n_total_patterns, det, sym, obs, i_full,
		              nthreads);
Thomas White's avatar
Thomas White committed
559
560
561
562
563

	}

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

Thomas White's avatar
Thomas White committed
565
566
567
	/* Clean up */
	free(i_full);
	delete_items(obs);
Thomas White's avatar
Thomas White committed
568
	free(sym);
Thomas White's avatar
Thomas White committed
569
	free(outfile);
570
571
572
573
574
575
576
	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);
577
578
579

	return 0;
}