partialator.c 12.5 KB
Newer Older
1
/*
2
 * partialator.c
3
 *
4
 * Scaling and post refinement for coherent nanocrystallography
5
 *
Thomas White's avatar
Thomas White committed
6
 * (c) 2006-2011 Thomas White <taw@physics.org>
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 <assert.h>
Thomas White's avatar
Tidy up    
Thomas White committed
24
#include <pthread.h>
25

Thomas White's avatar
Thomas White committed
26
#include "utils.h"
27
#include "hdf5-file.h"
Thomas White's avatar
Thomas White committed
28
#include "symmetry.h"
Thomas White's avatar
Thomas White committed
29
30
#include "reflections.h"
#include "stream.h"
31
#include "geometry.h"
32
#include "peaks.h"
Thomas White's avatar
Thomas White committed
33
#include "thread-pool.h"
Thomas White's avatar
Thomas White committed
34
#include "beam-parameters.h"
35
#include "post-refinement.h"
36
#include "hrs-scaling.h"
Thomas White's avatar
Thomas White committed
37
#include "reflist.h"
38
39


40
41
/* Maximum number of iterations of NLSq to do for each image per macrocycle. */
#define MAX_CYCLES (100)
42

43

44
45
46
47
static void show_help(const char *s)
{
	printf("Syntax: %s [options]\n\n", s);
	printf(
48
"Scaling and post refinement for coherent nanocrystallography.\n"
49
50
51
"\n"
"  -h, --help                 Display this help message.\n"
"\n"
Thomas White's avatar
Thomas White committed
52
53
"  -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
54
"  -o, --output=<filename>    Output filename.  Default: facetron.hkl.\n"
Thomas White's avatar
Thomas White committed
55
"  -g. --geometry=<file>      Get detector geometry from file.\n"
56
57
58
59
"  -b, --beam=<file>          Get beam parameters from file, which provides\n"
"                              initial values for parameters, and nominal\n"
"                              wavelengths if no per-shot value is found in \n"
"                              an HDF5 file.\n"
Thomas White's avatar
Thomas White committed
60
61
62
"  -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"
63
"  -y, --symmetry=<sym>       Merge according to symmetry <sym>.\n"
64
"  -n, --iterations=<n>       Run <n> cycles of scaling and post-refinement.\n"
65
"\n"
Thomas White's avatar
Thomas White committed
66
67
68
69
"  -j <n>                     Run <n> analyses in parallel.\n");
}


Thomas White's avatar
Thomas White committed
70
struct refine_args
Thomas White's avatar
Thomas White committed
71
{
Thomas White's avatar
Thomas White committed
72
73
74
75
	const char *sym;
	ReflItemList *obs;
	double *i_full;
	struct image *image;
76
	FILE *graph;
77
	FILE *pgraph;
Thomas White's avatar
Thomas White committed
78
79
80
81
82
83
84
};


static void refine_image(int mytask, void *tasks)
{
	struct refine_args *all_args = tasks;
	struct refine_args *pargs = &all_args[mytask];
85
86
87
	struct image *image = pargs->image;
	double nominal_photon_energy = pargs->image->beam->photon_energy;
	struct hdfile *hdfile;
Thomas White's avatar
Thomas White committed
88
	int i;
89
	double dev, last_dev;
Thomas White's avatar
Thomas White committed
90
	RefList *reflections;
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107

	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, nominal_photon_energy) ) {
		ERROR("Couldn't read '%s'\n", image->filename);
		hdfile_close(hdfile);
		return;
	}

Thomas White's avatar
Thomas White committed
108
109
110
111
112
113
	double a, b, c, al, be, ga;
	cell_get_parameters(image->indexed_cell, &a, &b, &c, &al, &be, &ga);
	STATUS("Initial cell: %5.2f %5.2f %5.2f nm %5.2f %5.2f %5.2f deg\n",
	       a/1.0e-9, b/1.0e-9, c/1.0e-9,
	       rad2deg(al), rad2deg(be), rad2deg(ga));

Thomas White's avatar
Thomas White committed
114
115
	/* FIXME: Don't do this each time */
	reflections = find_intersections(image, image->indexed_cell, 0);
116
117
118
119
	dev = +INFINITY;
	i = 0;
	do {
		last_dev = dev;
Thomas White's avatar
Thomas White committed
120
		dev = pr_iterate(image, pargs->i_full, pargs->sym, reflections);
121
		STATUS("Iteration %2i: mean dev = %5.2f\n", i, dev);
122
		i++;
123
	} while ( (fabs(last_dev - dev) > 1.0) && (i < MAX_CYCLES) );
Thomas White's avatar
Thomas White committed
124
	mean_partial_dev(image, reflections, pargs->sym,
125
	                 pargs->i_full, pargs->graph);
126
127
128
	if ( pargs->pgraph ) {
		fprintf(pargs->pgraph, "%5i %5.2f\n", mytask, dev);
	}
129
130
131
132

	free(image->data);
	if ( image->flags != NULL ) free(image->flags);
	hdfile_close(hdfile);
Thomas White's avatar
Thomas White committed
133
	reflist_free(reflections);
134
135
136
137

	/* Muppet proofing */
	image->data = NULL;
	image->flags = NULL;
138
139
140
}


141
142
143
144
static void refine_all(struct image *images, int n_total_patterns,
                       struct detector *det, const char *sym,
                       ReflItemList *obs, double *i_full, int nthreads,
                       FILE *graph, FILE *pgraph)
145
{
146
147
	struct refine_args *tasks;
	int i;
Thomas White's avatar
Thomas White committed
148

149
150
151
152
153
154
155
156
157
158
159
	tasks = malloc(n_total_patterns * sizeof(struct refine_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].image = &images[i];
		tasks[i].graph = graph;
		tasks[i].pgraph = pgraph;

	}
Thomas White's avatar
Thomas White committed
160

161
162
163
164
165
166
167
	run_thread_range(n_total_patterns, nthreads, "Refining",
	                 refine_image, tasks);

	free(tasks);
}


Thomas White's avatar
Thomas White committed
168
static void uniquify(Reflection *refl, const char *sym)
169
{
Thomas White's avatar
Thomas White committed
170
	signed int h, k, l;
171
172
	signed int ha, ka, la;

Thomas White's avatar
Thomas White committed
173
174
175
	get_indices(refl, &h, &k, &l);
	get_asymm(h, k, l, &ha, &ka, &la, sym);
	set_indices(refl, h, k, l);
176
177
178
}


Thomas White's avatar
Thomas White committed
179
/* FIXME: Get rid of this */
180
181
static void integrate_image(struct image *image, ReflItemList *obs,
                            const char *sym)
Thomas White's avatar
Thomas White committed
182
{
Thomas White's avatar
Thomas White committed
183
184
	RefList *reflections;
	Reflection *refl;
185
	struct hdfile *hdfile;
186
	double nominal_photon_energy = image->beam->photon_energy;
187

188
189
190
191
192
193
194
195
196
197
	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;
	}

198
	if ( hdf5_read(hdfile, image, 0, nominal_photon_energy) ) {
199
200
201
202
203
		ERROR("Couldn't read '%s'\n", image->filename);
		hdfile_close(hdfile);
		return;
	}

Thomas White's avatar
Thomas White committed
204
	/* Figure out which spots should appear in this pattern */
Thomas White's avatar
Thomas White committed
205
	reflections = find_intersections(image, image->indexed_cell, 0);
206
207

	/* For each reflection, estimate the partiality */
Thomas White's avatar
Thomas White committed
208
209
210
	for ( refl = first_refl(reflections);
	      refl != NULL;
	      refl = next_refl(refl) ) {
211
212
213
214

		signed int h, k, l;
		float i_partial;
		float xc, yc;
Thomas White's avatar
Thomas White committed
215
		double x, y;
216

Thomas White's avatar
Thomas White committed
217
218
		uniquify(refl, sym);
		get_indices(refl, &h, &k, &l);
219
220
221

		/* Don't attempt to use spots with very small
		 * partialities, since it won't be accurate. */
Thomas White's avatar
Thomas White committed
222
		if ( get_partiality(refl) < 0.1 ) continue;
223

Thomas White's avatar
Thomas White committed
224
225
226
		/* Actual measurement of this reflection from this pattern? */
		get_detector_pos(refl, &x, &y);
		if ( integrate_peak(image, x, y,
Thomas White's avatar
Thomas White committed
227
		                    &xc, &yc, &i_partial, NULL, NULL, 1, 0) ) {
Thomas White's avatar
Thomas White committed
228
			delete_refl(refl);
Thomas White's avatar
Thomas White committed
229
230
			continue;
		}
231

Thomas White's avatar
Thomas White committed
232
		set_int(refl, i_partial);
Thomas White's avatar
Thomas White committed
233
234
235

		if ( !find_item(obs, h, k, l) ) add_item(obs, h, k, l);

236
	}
Thomas White's avatar
Thomas White committed
237
	image->reflections = reflections;
238

239
240
	free(image->data);
	if ( image->flags != NULL ) free(image->flags);
241
	hdfile_close(hdfile);
Thomas White's avatar
Thomas White committed
242
243
244
245

	/* Muppet proofing */
	image->data = NULL;
	image->flags = NULL;
Thomas White's avatar
Thomas White committed
246
247
248
}


249
250
251
252
253
254
255
/* Decide which reflections can be scaled */
static void select_scalable_reflections(struct image *images, int n)
{
	int m;

	for ( m=0; m<n; m++ ) {

Thomas White's avatar
Thomas White committed
256
		Reflection *refl;
257

Thomas White's avatar
Thomas White committed
258
259
260
		for ( refl = first_refl(images[m].reflections);
		      refl != NULL;
		      refl = next_refl(refl) ) {
261
262

			int scalable = 1;
Thomas White's avatar
Thomas White committed
263
			double v;
264

Thomas White's avatar
Thomas White committed
265
266
267
			if ( get_partiality(refl) < 0.1 ) scalable = 0;
			v = fabs(get_intensity(refl));
			if ( v < 0.1 ) scalable = 0;
268

Thomas White's avatar
Thomas White committed
269
			set_scalable(refl, scalable);
270
271
272
273
274
275
276

		}

	}
}


Thomas White's avatar
Thomas White committed
277
278
279
280
int main(int argc, char *argv[])
{
	int c;
	char *infile = NULL;
Thomas White's avatar
Thomas White committed
281
	char *outfile = NULL;
Thomas White's avatar
Thomas White committed
282
283
	char *geomfile = NULL;
	char *prefix = NULL;
Thomas White's avatar
Thomas White committed
284
285
	char *sym = NULL;
	FILE *fh;
Thomas White's avatar
Thomas White committed
286
287
288
289
	int nthreads = 1;
	int config_basename = 0;
	int config_checkprefix = 1;
	struct detector *det;
Thomas White's avatar
Thomas White committed
290
	unsigned int *cts;
Thomas White's avatar
Thomas White committed
291
292
293
	ReflItemList *obs;
	int i;
	int n_total_patterns;
294
295
	struct image *images;
	int n_iter = 10;
Thomas White's avatar
Thomas White committed
296
	struct beam_params *beam = NULL;
Thomas White's avatar
Thomas White committed
297
	double *I_full;
Thomas White's avatar
Thomas White committed
298
299
300
301
302

	/* Long options */
	const struct option longopts[] = {
		{"help",               0, NULL,               'h'},
		{"input",              1, NULL,               'i'},
Thomas White's avatar
Thomas White committed
303
		{"output",             1, NULL,               'o'},
Thomas White's avatar
Thomas White committed
304
		{"geometry",           1, NULL,               'g'},
Thomas White's avatar
Thomas White committed
305
		{"beam",               1, NULL,               'b'},
Thomas White's avatar
Thomas White committed
306
307
308
		{"prefix",             1, NULL,               'x'},
		{"basename",           0, &config_basename,    1},
		{"no-check-prefix",    0, &config_checkprefix, 0},
309
		{"symmetry",           1, NULL,               'y'},
310
		{"iterations",         1, NULL,               'n'},
Thomas White's avatar
Thomas White committed
311
312
313
314
		{0, 0, NULL, 0}
	};

	/* Short options */
Thomas White's avatar
Thomas White committed
315
	while ((c = getopt_long(argc, argv, "hi:g:x:j:y:o:b:",
316
317
	                        longopts, NULL)) != -1)
	{
Thomas White's avatar
Thomas White committed
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339

		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;

340
341
342
343
		case 'y' :
			sym = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
344
345
346
347
		case 'o' :
			outfile = strdup(optarg);
			break;

348
349
350
351
		case 'n' :
			n_iter = atoi(optarg);
			break;

Thomas White's avatar
Thomas White committed
352
353
354
355
356
357
358
359
360
		case 'b' :
			beam = get_beam_parameters(optarg);
			if ( beam == NULL ) {
				ERROR("Failed to load beam parameters"
				      " from '%s'\n", optarg);
				return 1;
			}
			break;

Thomas White's avatar
Thomas White committed
361
362
363
364
365
366
367
368
369
		case 0 :
			break;

		default :
			return 1;
		}

	}

Thomas White's avatar
Thomas White committed
370
	/* Sanitise input filename and open */
Thomas White's avatar
Thomas White committed
371
372
373
374
375
376
377
378
379
380
381
382
383
384
	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
385
386
387
388
389
390
	/* Sanitise output filename */
	if ( outfile == NULL ) {
		outfile = strdup("facetron.hkl");
	}

	/* Sanitise prefix */
Thomas White's avatar
Thomas White committed
391
392
393
394
395
396
397
398
	if ( prefix == NULL ) {
		prefix = strdup("");
	} else {
		if ( config_checkprefix ) {
			prefix = check_prefix(prefix);
		}
	}

Thomas White's avatar
Thomas White committed
399
400
	if ( sym == NULL ) sym = strdup("1");

Thomas White's avatar
Thomas White committed
401
	/* Get detector geometry */
Thomas White's avatar
Thomas White committed
402
403
404
405
406
407
408
	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
409
410
411
412
413
	if ( beam == NULL ) {
		ERROR("You must provide a beam parameters file.\n");
		return 1;
	}

Thomas White's avatar
Thomas White committed
414
415
416
	n_total_patterns = count_patterns(fh);
	STATUS("There are %i patterns to process\n", n_total_patterns);

417
418
419
420
421
422
423
424
	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);
Thomas White's avatar
Thomas White committed
425
	obs = new_items();
426
427
428
429
	for ( i=0; i<n_total_patterns; i++ ) {

		UnitCell *cell;
		char *filename;
430
		char *fnamereal;
431
432
433
434
435
436
437
438
439
440
441
442

		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;
443
			tmp = safe_basename(filename);
444
445
446
			free(filename);
			filename = tmp;
		}
447
		fnamereal = malloc(1024);
448
		snprintf(fnamereal, 1023, "%s%s", prefix, filename);
449
		free(filename);
450
451

		images[i].filename = fnamereal;
Thomas White's avatar
Thomas White committed
452
453
		images[i].div = beam->divergence;
		images[i].bw = beam->bandwidth;
454
		images[i].det = det;
Thomas White's avatar
Thomas White committed
455
		images[i].beam = beam;
456
		images[i].osf = 1.0;
457
		images[i].profile_radius = 0.005e9;
458

Thomas White's avatar
Thomas White committed
459
460
461
462
		/* Muppet proofing */
		images[i].data = NULL;
		images[i].flags = NULL;

Thomas White's avatar
Thomas White committed
463
464
		/* Get reflections from this image.
		 * FIXME: Use the ones from the stream */
465
		integrate_image(&images[i], obs, sym);
466

467
468
469
470
471
472
		progress_bar(i, n_total_patterns-1, "Loading pattern data");

	}
	fclose(fh);
	free(prefix);

Thomas White's avatar
Thomas White committed
473
474
	cts = new_list_count();

475
	/* Make initial estimates */
476
	STATUS("Performing initial scaling.\n");
477
	select_scalable_reflections(images, n_total_patterns);
Thomas White's avatar
Thomas White committed
478
	I_full = scale_intensities(images, n_total_patterns, sym, obs);
479

Thomas White's avatar
Thomas White committed
480
	/* Iterate */
481
	for ( i=0; i<n_iter; i++ ) {
Thomas White's avatar
Thomas White committed
482

483
484
		FILE *fhg;
		FILE *fhp;
485
486
		char filename[1024];

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

489
490
491
492
493
494
495
496
497
498
		snprintf(filename, 1023, "p-iteration-%i.dat", i+1);
		fhg = fopen(filename, "w");
		if ( fhg == NULL ) {
			ERROR("Failed to open '%s'\n", filename);
			/* Nothing will be written later */
		}

		snprintf(filename, 1023, "g-iteration-%i.dat", i+1);
		fhp = fopen(filename, "w");
		if ( fhp == NULL ) {
499
500
501
502
			ERROR("Failed to open '%s'\n", filename);
			/* Nothing will be written later */
		}

Thomas White's avatar
Thomas White committed
503
		/* Refine the geometry of all patterns to get the best fit */
Thomas White's avatar
Thomas White committed
504
		refine_all(images, n_total_patterns, det, sym, obs, I_full,
505
		           nthreads, fhg, fhp);
Thomas White's avatar
Thomas White committed
506
507

		/* Re-estimate all the full intensities */
Thomas White's avatar
Thomas White committed
508
		free(I_full);
509
		select_scalable_reflections(images, n_total_patterns);
Thomas White's avatar
Thomas White committed
510
		I_full = scale_intensities(images, n_total_patterns, sym, obs);
Thomas White's avatar
Thomas White committed
511

512
513
		fclose(fhg);
		fclose(fhp);
Thomas White's avatar
Thomas White committed
514
515
516
517
518

	}

	STATUS("Final scale factors:\n");
	for ( i=0; i<n_total_patterns; i++ ) {
Thomas White's avatar
Thomas White committed
519
		STATUS("%4i : %5.2f\n", i, images[i].osf);
Thomas White's avatar
Thomas White committed
520
521
522
	}

	/* Output results */
Thomas White's avatar
Thomas White committed
523
	write_reflections(outfile, obs, I_full, NULL, NULL, cts, NULL);
Thomas White's avatar
Thomas White committed
524

Thomas White's avatar
Thomas White committed
525
	/* Clean up */
Thomas White's avatar
Thomas White committed
526
	for ( i=0; i<n_total_patterns; i++ ) {
Thomas White's avatar
Thomas White committed
527
		reflist_free(images[i].reflections);
Thomas White's avatar
Thomas White committed
528
529
	}
	free(I_full);
Thomas White's avatar
Thomas White committed
530
	delete_items(obs);
Thomas White's avatar
Thomas White committed
531
	free(sym);
Thomas White's avatar
Thomas White committed
532
	free(outfile);
533
534
	free(det->panels);
	free(det);
535
	free(beam);
Thomas White's avatar
Thomas White committed
536
	free(cts);
537
538
539
540
541
	for ( i=0; i<n_total_patterns; i++ ) {
		cell_free(images[i].indexed_cell);
		free(images[i].filename);
	}
	free(images);
542
543
544

	return 0;
}