facetron.c 16.4 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>
Thomas White's avatar
Tidy up    
Thomas White committed
24
#include <pthread.h>
25
26
27
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_linalg.h>
28

Thomas White's avatar
Thomas White committed
29
#include "utils.h"
30
#include "hdf5-file.h"
Thomas White's avatar
Thomas White committed
31
#include "symmetry.h"
Thomas White's avatar
Thomas White committed
32
33
#include "reflections.h"
#include "stream.h"
34
#include "geometry.h"
35
#include "peaks.h"
Thomas White's avatar
Thomas White committed
36
#include "thread-pool.h"
Thomas White's avatar
Thomas White committed
37
#include "beam-parameters.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
/* Refineable parameters */
enum {
	REF_SCALE,
Thomas White's avatar
Thomas White committed
46
	REF_DIV,
47
48
49
50
	NUM_PARAMS
};


51
52
53
54
static void show_help(const char *s)
{
	printf("Syntax: %s [options]\n\n", s);
	printf(
Thomas White's avatar
Thomas White committed
55
"Post-refinement and profile fitting for coherent nanocrystallography.\n"
56
57
58
"\n"
"  -h, --help                 Display this help message.\n"
"\n"
Thomas White's avatar
Thomas White committed
59
60
"  -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
61
"  -o, --output=<filename>    Output filename.  Default: facetron.hkl.\n"
Thomas White's avatar
Thomas White committed
62
"  -g. --geometry=<file>      Get detector geometry from file.\n"
Thomas White's avatar
Thomas White committed
63
64
65
"  -b, --beam=<file>          Get beam parameters from file (provides initial\n"
"                              values for parameters, and nominal wavelengths\n"
"                              if no per-shot value is found in an HDF5 file.\n"
Thomas White's avatar
Thomas White committed
66
67
68
"  -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"
69
70
71
"  -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
72
73
74
75
"  -j <n>                     Run <n> analyses in parallel.\n");
}


Thomas White's avatar
Thomas White committed
76
struct refine_args
Thomas White's avatar
Thomas White committed
77
{
Thomas White's avatar
Thomas White committed
78
79
80
81
	const char *sym;
	ReflItemList *obs;
	double *i_full;
	struct image *image;
82
	FILE *graph;
Thomas White's avatar
Thomas White committed
83
84
85
};


86
87
88
89
/* Return the gradient of parameter 'k' given the current status of 'image'. */
static double gradient(struct image *image, int k,
                       struct cpeak spot, double I_partial)
{
Thomas White's avatar
Thomas White committed
90
91
92
93
94
	double ds;
	double nom, den;

	ds = 2.0 * resolution(image->indexed_cell, spot.h, spot.k, spot.l);

95
96
97
98
99
	switch ( k ) {

	case REF_SCALE :
		return I_partial;

Thomas White's avatar
Thomas White committed
100
101
102
103
104
	case REF_DIV :
		nom = sqrt(2.0) * ds * sin(image->div);
		den = sqrt(1.0 - cos(image->div));
		return nom/den;

105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
	}

	ERROR("No gradient defined for parameter %i\n", k);
	abort();
}


/* Apply the given shift to the 'k'th parameter of 'image'. */
static void apply_shift(struct image *image, int k, double shift)
{
	switch ( k ) {

	case REF_SCALE :
		image->osf += shift;
		break;

Thomas White's avatar
Thomas White committed
121
122
123
124
125
	case REF_DIV :
		STATUS("Shifting div by %e\n", shift);
		image->div += shift;
		break;

126
	default :
Thomas White's avatar
Thomas White committed
127
		ERROR("No shift defined for parameter %i\n", k);
128
129
130
131
132
133
134
		abort();

	}
}


static double mean_partial_dev(struct image *image, struct cpeak *spots, int n,
135
                               const char *sym, double *i_full, FILE *graph)
136
137
{
	int h;
138
	double delta_I = 0.0;
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162

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

		signed int hind, kind, lind;
		signed int ha, ka, la;
		double I_full;
		float I_partial;
		float xc, yc;

		hind = spots[h].h;
		kind = spots[h].k;
		lind = spots[h].l;

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

		/* Actual measurement of this reflection from this
		 * pattern? */
		/* FIXME: Coordinates aren't whole numbers */
		if ( integrate_peak(image, spots[h].x, spots[h].y,
		                    &xc, &yc, &I_partial, NULL, NULL, 1, 1) ) {
			continue;
		}
163
		I_partial *= image->osf;
164
165
166

		get_asymm(hind, kind, lind, &ha, &ka, &la, sym);
		I_full = lookup_intensity(i_full, ha, ka, la);
167
		delta_I += fabs(I_partial - spots[h].p * I_full);
168

169
		if ( graph != NULL ) {
Thomas White's avatar
Thomas White committed
170
171
172
			fprintf(graph, "%3i %3i %3i %5.2f (at %5.2f,%5.2f)"
			               " %5.2f %5.2f\n",
			       hind, kind, lind, I_partial/spots[h].p, xc, yc,
173
174
175
			       spots[h].p, I_partial / I_full);
		}

176
177
178
179
180
181
182
	}

	return delta_I / (double)n;
}


static double iterate(struct image *image, double *i_full, const char *sym,
183
                      struct cpeak **pspots, int *n)
184
185
186
187
{
	gsl_matrix *M;
	gsl_vector *v;
	gsl_vector *shifts;
188
	int h, param;
189
	struct cpeak *spots = *pspots;
190

191
192
	M = gsl_matrix_calloc(NUM_PARAMS, NUM_PARAMS);
	v = gsl_vector_calloc(NUM_PARAMS);
193
194

	/* Construct the equations, one per reflection in this image */
195
	for ( h=0; h<*n; h++ ) {
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218

		signed int hind, kind, lind;
		signed int ha, ka, la;
		double I_full, delta_I;
		float I_partial;
		float xc, yc;
		int k;

		hind = spots[h].h;
		kind = spots[h].k;
		lind = spots[h].l;

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

		/* Actual measurement of this reflection from this
		 * pattern? */
		/* FIXME: Coordinates aren't whole numbers */
		if ( integrate_peak(image, spots[h].x, spots[h].y,
		                    &xc, &yc, &I_partial, NULL, NULL, 1, 1) ) {
			continue;
		}
219
		I_partial *= image->osf;
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253

		get_asymm(hind, kind, lind, &ha, &ka, &la, sym);
		I_full = lookup_intensity(i_full, ha, ka, la);
		delta_I = I_partial - spots[h].p * I_full;

		for ( k=0; k<NUM_PARAMS; k++ ) {

			int g;
			double v_c;

			for ( g=0; g<NUM_PARAMS; g++ ) {

				double M_curr, M_c;

				M_curr = gsl_matrix_get(M, g, k);

				M_c = gradient(image, g, spots[h], I_partial)
				    * gradient(image, k, spots[h], I_partial);
				M_c *= pow(I_full, 2.0);

				gsl_matrix_set(M, g, k, M_curr + M_c);

			}

			v_c = delta_I * I_full * gradient(image, k, spots[h],
			                                  I_partial);
			gsl_vector_set(v, k, v_c);

		}

	}

	shifts = gsl_vector_alloc(NUM_PARAMS);
	gsl_linalg_HH_solve(M, v, shifts);
254
255
256
	for ( param=0; param<NUM_PARAMS; param++ ) {
		double shift = gsl_vector_get(shifts, param);
		apply_shift(image, param, shift);
257
258
	}

259
260
261
262
	gsl_matrix_free(M);
	gsl_vector_free(v);
	gsl_vector_free(shifts);

263
	free(spots);
264
265
266
	spots = find_intersections(image, image->indexed_cell, n, 0);
	*pspots = spots;
	return mean_partial_dev(image, spots, *n, sym, i_full, NULL);
267
268
269
}


Thomas White's avatar
Thomas White committed
270
271
272
273
static void refine_image(int mytask, void *tasks)
{
	struct refine_args *all_args = tasks;
	struct refine_args *pargs = &all_args[mytask];
274
275
276
277
	struct image *image = pargs->image;
	double nominal_photon_energy = pargs->image->beam->photon_energy;
	struct hdfile *hdfile;
	struct cpeak *spots;
278
	int n, i;
279
	double dev, last_dev;
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297

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

	spots = find_intersections(image, image->indexed_cell, &n, 0);
298
299
300
301
	dev = +INFINITY;
	i = 0;
	do {
		last_dev = dev;
302
		dev = iterate(image, pargs->i_full, pargs->sym, &spots, &n);
303
		STATUS("Iteration %2i: mean dev = %5.2f\n", i, dev);
304
		i++;
305
	} while ( (fabs(last_dev - dev) > 1.0) && (i < MAX_CYCLES) );
306
307
	mean_partial_dev(image, spots, n, pargs->sym,
	                 pargs->i_full, pargs->graph);
308
309
310
311

	free(image->data);
	if ( image->flags != NULL ) free(image->flags);
	hdfile_close(hdfile);
312
	free(spots);
313
314
315
316

	/* Muppet proofing */
	image->data = NULL;
	image->flags = NULL;
317
318
319
}


Thomas White's avatar
Thomas White committed
320
struct integrate_args
321
{
Thomas White's avatar
Thomas White committed
322
323
324
325
326
327
328
329
330
331
332
333
334
	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];
Thomas White's avatar
Thomas White committed
335
	struct cpeak *spots;
336
337
338
	int j, n;
	struct hdfile *hdfile;
	struct image *image = pargs->image;
Thomas White's avatar
Thomas White committed
339
	double nominal_photon_energy = pargs->image->beam->photon_energy;
340

341
342
343
344
345
346
347
348
349
350
	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;
	}

Thomas White's avatar
Thomas White committed
351
	if ( hdf5_read(hdfile, pargs->image, 0, nominal_photon_energy) ) {
352
353
354
355
356
		ERROR("Couldn't read '%s'\n", image->filename);
		hdfile_close(hdfile);
		return;
	}

Thomas White's avatar
Thomas White committed
357
	/* Figure out which spots should appear in this pattern */
358
	spots = find_intersections(image, image->indexed_cell, &n, 0);
359
360
361
362
363

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

		signed int h, k, l;
364
		signed int ha, ka, la;
365
366
		float i_partial;
		float xc, yc;
Thomas White's avatar
Thomas White committed
367
		float i_full_est;
368
369
370
371
372
373
374

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

		/* Don't attempt to use spots with very small
		 * partialities, since it won't be accurate. */
Thomas White's avatar
Thomas White committed
375
		if ( spots[j].p < 0.1 ) continue;
376
377
378
379
380

		/* Actual measurement of this reflection from this
		 * pattern? */
		/* FIXME: Coordinates aren't whole numbers */
		if ( integrate_peak(image, spots[j].x, spots[j].y,
Thomas White's avatar
Thomas White committed
381
382
383
		                    &xc, &yc, &i_partial, NULL, NULL, 1, 1) ) {
			continue;
		}
384
		i_partial *= image->osf;
385

Thomas White's avatar
Thomas White committed
386
		i_full_est = i_partial / spots[j].p;
Thomas White's avatar
Thomas White committed
387

388
389
		get_asymm(h, k, l, &ha, &ka, &la, pargs->sym);

390
		pthread_mutex_lock(pargs->list_lock);
391
392
393
394
		integrate_intensity(pargs->i_full, ha, ka, la, i_full_est);
		integrate_count(pargs->cts, ha, ka, la, 1);
		if ( !find_item(pargs->obs, ha, ka, la) ) {
			add_item(pargs->obs, ha, ka, la);
395
396
397
398
399
		}
		pthread_mutex_unlock(pargs->list_lock);

	}

400
401
	free(image->data);
	if ( image->flags != NULL ) free(image->flags);
402
	hdfile_close(hdfile);
Thomas White's avatar
Thomas White committed
403
404
405
406
407
	free(spots);

	/* Muppet proofing */
	image->data = NULL;
	image->flags = NULL;
Thomas White's avatar
Thomas White committed
408
409
410
}


Thomas White's avatar
Thomas White committed
411
412
static void refine_all(struct image *images, int n_total_patterns,
                       struct detector *det, const char *sym,
413
414
                       ReflItemList *obs, double *i_full, int nthreads,
                       FILE *graph)
Thomas White's avatar
Thomas White committed
415
{
Thomas White's avatar
Thomas White committed
416
	struct refine_args *tasks;
Thomas White's avatar
Thomas White committed
417
	int i;
418

Thomas White's avatar
Thomas White committed
419
420
	tasks = malloc(n_total_patterns * sizeof(struct refine_args));
	for ( i=0; i<n_total_patterns; i++ ) {
Thomas White's avatar
Thomas White committed
421

Thomas White's avatar
Thomas White committed
422
423
424
425
		tasks[i].sym = sym;
		tasks[i].obs = obs;
		tasks[i].i_full = i_full;
		tasks[i].image = &images[i];
426
		tasks[i].graph = graph;
Thomas White's avatar
Thomas White committed
427
428
429

	}

430
431
	run_thread_range(n_total_patterns, nthreads, "Refining",
	                 refine_image, tasks);
Thomas White's avatar
Thomas White committed
432

Thomas White's avatar
Thomas White committed
433
	free(tasks);
434
435
436
437
438
}


static void estimate_full(struct image *images, int n_total_patterns,
                          struct detector *det, const char *sym,
Thomas White's avatar
Thomas White committed
439
440
                          ReflItemList *obs, double *i_full, unsigned int *cts,
                          int nthreads)
441
442
{
	int i;
Thomas White's avatar
Thomas White committed
443
444
	struct integrate_args *tasks;
	pthread_mutex_t list_lock = PTHREAD_MUTEX_INITIALIZER;
445

446
	clear_items(obs);
447

Thomas White's avatar
Thomas White committed
448
449
450
451
452
453
454
455
456
457
458
459
	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];

	}

460
461
	run_thread_range(n_total_patterns, nthreads, "Integrating",
	                 integrate_image, tasks);
Thomas White's avatar
Thomas White committed
462
463

	free(tasks);
464

465
466
	/* Divide the totals to get the means */
	for ( i=0; i<num_items(obs); i++ ) {
467

468
469
		struct refl_item *it;
		double total;
470

471
472
473
474
		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);
475

476
	}
Thomas White's avatar
Thomas White committed
477
}
478

Thomas White's avatar
Thomas White committed
479

Thomas White's avatar
Thomas White committed
480
481
482
483
int main(int argc, char *argv[])
{
	int c;
	char *infile = NULL;
Thomas White's avatar
Thomas White committed
484
	char *outfile = NULL;
Thomas White's avatar
Thomas White committed
485
486
	char *geomfile = NULL;
	char *prefix = NULL;
Thomas White's avatar
Thomas White committed
487
488
	char *sym = NULL;
	FILE *fh;
Thomas White's avatar
Thomas White committed
489
490
491
492
	int nthreads = 1;
	int config_basename = 0;
	int config_checkprefix = 1;
	struct detector *det;
Thomas White's avatar
Thomas White committed
493
	double *i_full;
Thomas White's avatar
Thomas White committed
494
	unsigned int *cts;
Thomas White's avatar
Thomas White committed
495
496
497
	ReflItemList *obs;
	int i;
	int n_total_patterns;
498
499
	struct image *images;
	int n_iter = 10;
Thomas White's avatar
Thomas White committed
500
	struct beam_params *beam = NULL;
Thomas White's avatar
Thomas White committed
501
502
503
504
505

	/* Long options */
	const struct option longopts[] = {
		{"help",               0, NULL,               'h'},
		{"input",              1, NULL,               'i'},
Thomas White's avatar
Thomas White committed
506
		{"output",             1, NULL,               'o'},
Thomas White's avatar
Thomas White committed
507
		{"geometry",           1, NULL,               'g'},
Thomas White's avatar
Thomas White committed
508
		{"beam",               1, NULL,               'b'},
Thomas White's avatar
Thomas White committed
509
510
511
		{"prefix",             1, NULL,               'x'},
		{"basename",           0, &config_basename,    1},
		{"no-check-prefix",    0, &config_checkprefix, 0},
512
		{"symmetry",           1, NULL,               'y'},
513
		{"iterations",         1, NULL,               'n'},
Thomas White's avatar
Thomas White committed
514
515
516
517
		{0, 0, NULL, 0}
	};

	/* Short options */
Thomas White's avatar
Thomas White committed
518
	while ((c = getopt_long(argc, argv, "hi:g:x:j:y:o:b:",
519
520
	                        longopts, NULL)) != -1)
	{
Thomas White's avatar
Thomas White committed
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542

		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;

543
544
545
546
		case 'y' :
			sym = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
547
548
549
550
		case 'o' :
			outfile = strdup(optarg);
			break;

551
552
553
554
		case 'n' :
			n_iter = atoi(optarg);
			break;

Thomas White's avatar
Thomas White committed
555
556
557
558
559
560
561
562
563
		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
564
565
566
567
568
569
570
571
572
		case 0 :
			break;

		default :
			return 1;
		}

	}

Thomas White's avatar
Thomas White committed
573
	/* Sanitise input filename and open */
Thomas White's avatar
Thomas White committed
574
575
576
577
578
579
580
581
582
583
584
585
586
587
	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
588
589
590
591
592
593
	/* Sanitise output filename */
	if ( outfile == NULL ) {
		outfile = strdup("facetron.hkl");
	}

	/* Sanitise prefix */
Thomas White's avatar
Thomas White committed
594
595
596
597
598
599
600
601
	if ( prefix == NULL ) {
		prefix = strdup("");
	} else {
		if ( config_checkprefix ) {
			prefix = check_prefix(prefix);
		}
	}

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

Thomas White's avatar
Thomas White committed
604
	/* Get detector geometry */
Thomas White's avatar
Thomas White committed
605
606
607
608
609
610
611
	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
612
613
614
615
616
	if ( beam == NULL ) {
		ERROR("You must provide a beam parameters file.\n");
		return 1;
	}

Thomas White's avatar
Thomas White committed
617
618
619
620
621
622
623
	/* 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);

624
625
626
627
628
629
630
631
632
633
634
635
	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;
636
		char *fnamereal;
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652

		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;
		}
653
		fnamereal = malloc(1024);
654
655
656
		snprintf(fnamereal, 1023, "%s%s", prefix, filename);

		images[i].filename = fnamereal;
Thomas White's avatar
Thomas White committed
657
658
		images[i].div = beam->divergence;
		images[i].bw = beam->bandwidth;
659
		images[i].det = det;
Thomas White's avatar
Thomas White committed
660
		images[i].beam = beam;
661
		images[i].osf = 1.0;
662

Thomas White's avatar
Thomas White committed
663
664
665
666
		/* Muppet proofing */
		images[i].data = NULL;
		images[i].flags = NULL;

667
668
669
670
671
672
673
674
		free(filename);

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

	}
	fclose(fh);
	free(prefix);

Thomas White's avatar
Thomas White committed
675
676
	cts = new_list_count();

677
	/* Make initial estimates */
Thomas White's avatar
Thomas White committed
678
	estimate_full(images, n_total_patterns, det, sym, obs, i_full, cts,
679
	              nthreads);
680

Thomas White's avatar
Thomas White committed
681
	/* Iterate */
682
	for ( i=0; i<n_iter; i++ ) {
Thomas White's avatar
Thomas White committed
683

684
685
686
		FILE *fh;
		char filename[1024];

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

689
690
691
692
693
694
695
		snprintf(filename, 1023, "iteration-%i.dat", i+1);
		fh = fopen(filename, "w");
		if ( fh == NULL ) {
			ERROR("Failed to open '%s'\n", filename);
			/* Nothing will be written later */
		}

Thomas White's avatar
Thomas White committed
696
		/* Refine the geometry of all patterns to get the best fit */
697
		refine_all(images, n_total_patterns, det, sym, obs, i_full,
698
		           nthreads, fh);
Thomas White's avatar
Thomas White committed
699
700

		/* Re-estimate all the full intensities */
701
		estimate_full(images, n_total_patterns, det, sym, obs, i_full,
Thomas White's avatar
Thomas White committed
702
		              cts, nthreads);
Thomas White's avatar
Thomas White committed
703

704
		fclose(fh);
Thomas White's avatar
Thomas White committed
705
706
707
	}

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

Thomas White's avatar
Thomas White committed
710
711
712
	/* Clean up */
	free(i_full);
	delete_items(obs);
Thomas White's avatar
Thomas White committed
713
	free(sym);
Thomas White's avatar
Thomas White committed
714
	free(outfile);
715
716
	free(det->panels);
	free(det);
717
	free(beam);
Thomas White's avatar
Thomas White committed
718
	free(cts);
719
720
721
722
723
	for ( i=0; i<n_total_patterns; i++ ) {
		cell_free(images[i].indexed_cell);
		free(images[i].filename);
	}
	free(images);
724
725
726

	return 0;
}