process_hkl.c 13.7 KB
Newer Older
1
2
3
4
5
/*
 * process_hkl.c
 *
 * Assemble and process FEL Bragg intensities
 *
6
7
 * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY,
 *                  a research centre of the Helmholtz Association.
8
 * Copyright © 2012 Lorenzo Galli
9
10
11
 *
 * Authors:
 *   2009-2012 Thomas White <taw@physics.org>
12
13
 *   2011      Andrew Martin <andrew.martin@desy.de>
 *   2012      Lorenzo Galli <lorenzo.galli@desy.de>
14
 *
Thomas White's avatar
Thomas White committed
15
16
17
18
19
20
21
22
23
24
25
26
27
28
 * This file is part of CrystFEL.
 *
 * CrystFEL is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * CrystFEL is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with CrystFEL.  If not, see <http://www.gnu.org/licenses/>.
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
 *
 */


#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>

#include "utils.h"
#include "statistics.h"
46
#include "reflist-utils.h"
47
#include "symmetry.h"
48
#include "stream.h"
49
50
#include "reflist.h"
#include "image.h"
51
52
53
54
55


static void show_help(const char *s)
{
	printf("Syntax: %s [options]\n\n", s);
Thomas White's avatar
Thomas White committed
56
57
58
	printf(
"Assemble and process FEL Bragg intensities.\n"
"\n"
59
60
61
"  -h, --help                Display this help message.\n"
"  -i, --input=<filename>    Specify input filename (\"-\" for stdin).\n"
"  -o, --output=<filename>   Specify output filename for merged intensities\n"
62
"                             Default: processed.hkl).\n"
Thomas White's avatar
Thomas White committed
63
"  -y, --symmetry=<sym>      Merge according to point group <sym>.\n"
Thomas White's avatar
Thomas White committed
64
"\n"
65
"      --start-after=<n>     Skip n patterns at the start of the stream.\n"
Thomas White's avatar
Thomas White committed
66
"      --stop-after=<n>      Stop after processing n patterns.\n"
67
68
"  -g, --histogram=<h,k,l>   Calculate the histogram of measurements for this\n"
"                             reflection.\n"
69
"  -z, --hist-parameters     Set the range for the histogram and the number of\n"
Thomas White's avatar
Thomas White committed
70
"          =<min,max,nbins>   bins. \n"
71
"\n"
72
73
"      --scale               Scale each pattern for best fit with the current\n"
"                             model.\n"
74
"      --reference=<file>    Compare against intensities from <file> when\n"
Thomas White's avatar
Thomas White committed
75
"                             scaling. \n"
76
"      --no-polarisation     Disable polarisation correction.\n"
77
78
"      --min-measurements=<n> Require at least <n> measurements before a\n"
"                             reflection appears in the output.  Default: 2\n"
79
);
80
81
82
}


Thomas White's avatar
Thomas White committed
83
84
static void plot_histogram(double *vals, int n, float hist_min, float hist_max,
                           int nbins)
85
86
87
88
89
{
	int i;
	double max = -INFINITY;
	double min = +INFINITY;
	double step;
90
	int histo[nbins];
91
92
93
94
95
96
97
98
	FILE *fh;

	fh = fopen("histogram.dat", "w");
	if ( fh == NULL ) {
		ERROR("Couldn't open 'histogram.dat'\n");
		return;
	}

99
100
101
102
103
104
105
106
	if ( hist_min == hist_max ) {
		for ( i=0; i<n; i++ ) {
			if ( vals[i] > max ) max = vals[i];
			if ( vals[i] < min ) min = vals[i];
		}
	} else {
		min = hist_min;
		max = hist_max;
107
	}
Thomas White's avatar
Thomas White committed
108
	STATUS("min max nbins: %f %f %i\n", min, max, nbins);
109
	min--;  max++;
110

111
	for ( i=0; i<nbins; i++ ) {
112
113
114
		histo[i] = 0;
	}

115
	step = (max-min)/nbins;
116
117
118

	for ( i=0; i<n; i++ ) {
		int bin;
119
120
121
122
		if ( (vals[i] > min) && (vals[i] < max) ) {
			bin = (vals[i]-min)/step;
			histo[bin]++;
		}
123
124
	}

125
	for ( i=0; i<nbins; i++ ) {
126
127
128
129
130
131
132
		fprintf(fh, "%f %i\n", min+step*i, histo[i]);
	}

	fclose(fh);
}


133
134
static double scale_intensities(RefList *reference, RefList *new,
                              const SymOpList *sym)
135
{
136
137
138
	double s;
	double top = 0.0;
	double bot = 0.0;
139
140
	Reflection *refl;
	RefListIterator *iter;
141

142
143
	for ( refl = first_refl(new, &iter);
	      refl != NULL;
144
145
	      refl = next_refl(refl, iter) )
	{
146

147
148
		double i1, i2;
		signed int hu, ku, lu;
149
		signed int h, k, l;
150
		Reflection *reference_version;
151

152
		get_indices(refl, &h, &k, &l);
153
		get_asymm(sym, h, k, l, &hu, &ku, &lu);
Thomas White's avatar
Thomas White committed
154

155
156
		reference_version = find_refl(reference, hu, ku, lu);
		if ( reference_version == NULL ) continue;
157

158
159
		i1 = get_intensity(reference_version);
		i2 = get_intensity(refl);
160

161
162
163
		/* Calculate LSQ estimate of scaling factor */
		top += i1 * i2;
		bot += i2 * i2;
164

165
	}
166

167
	s = top / bot;
168

169
	return s;
170
171
172
}


173
174
175
static int merge_pattern(RefList *model, struct image *new, RefList *reference,
                         const SymOpList *sym,
                         double *hist_vals, signed int hist_h,
176
177
                         signed int hist_k, signed int hist_l, int *hist_n,
                         int config_nopolar)
178
{
179
180
	Reflection *refl;
	RefListIterator *iter;
181
182
183
184
	double scale;
	double asx, asy, asz;
	double bsx, bsy, bsz;
	double csx, csy, csz;
185

186
187
188
189
190
191
	if ( reference != NULL ) {
		scale = scale_intensities(reference, new->reflections, sym);
	} else {
		scale = 1.0;
	}
	if ( isnan(scale) ) return 1;
192

193
194
195
196
197
198
199
200
201
202
203
	cell_get_reciprocal(new->indexed_cell, &asx, &asy, &asz,
	                                       &bsx, &bsy, &bsz,
	                                       &csx, &csy, &csz);

	for ( refl = first_refl(new->reflections, &iter);
	      refl != NULL;
	      refl = next_refl(refl, iter) )
	{
		double intensity;
		double xl, yl, zl;
		double pol, pa, pb, phi, tt, ool;
204
		signed int h, k, l;
205
206
207
		int cur_redundancy;
		double cur_intensity, cur_sumsq;
		Reflection *model_version;
208

209
		get_indices(refl, &h, &k, &l);
210

211
212
		/* Put into the asymmetric unit for the target group */
		get_asymm(sym, h, k, l, &h, &k, &l);
213

214
215
216
217
		model_version = find_refl(model, h, k, l);
		if ( model_version == NULL ) {
			model_version = add_refl(model, h, k, l);
		}
218

219
		intensity = scale * get_intensity(refl);
220

221
		if ( !config_nopolar ) {
222

223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
			/* Polarisation correction assuming 100% polarisation
			 * along the x direction */
			xl = h*asx + k*bsx + l*csx;
			yl = h*asy + k*bsy + l*csy;
			zl = h*asz + k*bsz + l*csz;

			ool = 1.0 / new->lambda;
			tt = angle_between(0.0, 0.0, 1.0,  xl, yl, zl+ool);
			phi = atan2(yl, xl);
			pa = pow(sin(phi)*sin(tt), 2.0);
			pb = pow(cos(tt), 2.0);
			pol = 1.0 - 2.0*(1.0-pa) + (1.0+pb);
			intensity /= pol;

		}
238

239
240
		cur_intensity = get_intensity(model_version);
		set_intensity(model_version, cur_intensity + intensity);
241

242
243
		cur_redundancy = get_redundancy(model_version);
		set_redundancy(model_version, cur_redundancy+1);
244

245
246
		cur_sumsq = get_temp1(model_version);
		set_temp1(model_version, cur_sumsq + pow(intensity, 2.0));
Thomas White's avatar
Thomas White committed
247

248
		if ( hist_vals != NULL ) {
Thomas White's avatar
Thomas White committed
249

250
251
252
253
			if ( (h==hist_h) && (k==hist_k) && (l==hist_l) ) {
				hist_vals[*hist_n] = intensity;
				*hist_n += 1;
			}
Thomas White's avatar
Thomas White committed
254

255
		}
Thomas White's avatar
Thomas White committed
256

257
	}
258

259
	return 0;
260
261
262
}


263
static void merge_all(FILE *fh, RefList *model, RefList *reference,
264
                      int config_startafter, int config_stopafter,
Thomas White's avatar
Thomas White committed
265
                      const SymOpList *sym,
266
                      int n_total_patterns,
267
268
                      double *hist_vals, signed int hist_h,
                      signed int hist_k, signed int hist_l,
269
                      int *hist_i, int config_nopolar, int min_measurements)
270
{
271
	int rval;
Thomas White's avatar
Thomas White committed
272
	int n_patterns = 0;
273
274
275
	int n_used = 0;
	Reflection *refl;
	RefListIterator *iter;
276

277
278
279
	if ( skip_some_files(fh, config_startafter) ) {
		ERROR("Failed to skip first %i files.\n", config_startafter);
		return;
280
281
	}

282
283
	do {

284
		struct image image;
285

286
287
		image.det = NULL;

288
289
		/* Get data from next chunk */
		rval = read_chunk(fh, &image);
Thomas White's avatar
Thomas White committed
290
		if ( rval ) break;
291

292
293
294
		n_patterns++;

		if ( (image.reflections != NULL) && (image.indexed_cell) ) {
295

296
			int r;
297

298
299
			r = merge_pattern(model, &image, reference, sym,
			                  hist_vals, hist_h, hist_k, hist_l,
300
			                  hist_i, config_nopolar);
301

302
			if ( r == 0 ) n_used++;
303
304
305

		}

Thomas White's avatar
Thomas White committed
306
		free(image.filename);
307
308
309
		reflist_free(image.reflections);
		image_feature_list_free(image.features);
		cell_free(image.indexed_cell);
310

311
312
		progress_bar(n_patterns, n_total_patterns-config_startafter,
		             "Merging");
313

314
315
316
317
		if ( config_stopafter ) {
			if ( n_patterns == config_stopafter ) break;
		}

318
	} while ( rval == 0 );
319

320
321
322
323
324
325
326
327
	for ( refl = first_refl(model, &iter);
	      refl != NULL;
	      refl = next_refl(refl, iter) )
	{
		double intensity, sumsq, esd;
		int red;

		red = get_redundancy(refl);
328
		if ( red < min_measurements ) {
329
330
			set_redundancy(refl, 0);
			continue;
Thomas White's avatar
Thomas White committed
331
		}
332

333
334
		intensity = get_intensity(refl) / red;
		set_intensity(refl, intensity);
Thomas White's avatar
Thomas White committed
335

336
337
338
		sumsq = get_temp1(refl) / red;
		esd = sqrt(sumsq - pow(intensity, 2.0)) / sqrt(red);
		set_esd_intensity(refl, esd);
339

340
	}
Thomas White's avatar
Thomas White committed
341

342
	STATUS("%i of the patterns could be used.\n", n_used);
343
344
345
}


346
347
348
349
int main(int argc, char *argv[])
{
	int c;
	char *filename = NULL;
350
	char *output = NULL;
351
	FILE *fh;
352
	RefList *model;
Thomas White's avatar
Thomas White committed
353
	int config_maxonly = 0;
354
	int config_startafter = 0;
Thomas White's avatar
Thomas White committed
355
	int config_stopafter = 0;
356
	int config_sum = 0;
357
	int config_scale = 0;
358
	unsigned int n_total_patterns;
Thomas White's avatar
Thomas White committed
359
360
	char *sym_str = NULL;
	SymOpList *sym;
361
362
	char *histo = NULL;
	signed int hist_h, hist_k, hist_l;
363
364
	signed int hist_nbins=50;
	float hist_min=0.0, hist_max=0.0;
365
366
	double *hist_vals = NULL;
	int hist_i;
367
	int space_for_hist = 0;
368
	char *histo_params = NULL;
369
	int config_nopolar = 0;
370
371
	char *rval;
	int min_measurements = 2;
372
373
374
375
376

	/* Long options */
	const struct option longopts[] = {
		{"help",               0, NULL,               'h'},
		{"input",              1, NULL,               'i'},
Thomas White's avatar
Thomas White committed
377
		{"output",             1, NULL,               'o'},
Thomas White's avatar
Thomas White committed
378
		{"max-only",           0, &config_maxonly,     1},
Thomas White's avatar
Thomas White committed
379
		{"output-every",       1, NULL,               'e'},
Thomas White's avatar
Thomas White committed
380
		{"stop-after",         1, NULL,               's'},
381
		{"start-after",        1, NULL,               'f'},
382
		{"sum",                0, &config_sum,         1},
383
		{"scale",              0, &config_scale,       1},
384
385
		{"no-polarisation",    0, &config_nopolar,     1},
		{"no-polarization",    0, &config_nopolar,     1},
386
		{"symmetry",           1, NULL,               'y'},
387
		{"histogram",          1, NULL,               'g'},
388
		{"hist-parameters",    1, NULL,               'z'},
389
		{"min-measurements",   1, NULL,                2},
390
391
392
393
		{0, 0, NULL, 0}
	};

	/* Short options */
394
	while ((c = getopt_long(argc, argv, "hi:e:o:y:g:f:b:z:",
395
	                        longopts, NULL)) != -1) {
396
397

		switch (c) {
398
399

			case 'h' :
400
401
402
			show_help(argv[0]);
			return 0;

403
			case 'i' :
404
405
406
			filename = strdup(optarg);
			break;

407
			case 'o' :
408
409
410
			output = strdup(optarg);
			break;

411
			case 's' :
Thomas White's avatar
Thomas White committed
412
413
414
			config_stopafter = atoi(optarg);
			break;

415
			case 'f' :
416
417
418
			config_startafter = atoi(optarg);
			break;

419
			case 'y' :
Thomas White's avatar
Thomas White committed
420
			sym_str = strdup(optarg);
421
422
			break;

423
			case 'g' :
424
425
426
			histo = strdup(optarg);
			break;

427
			case 'z' :
428
429
430
			histo_params = strdup(optarg);
			break;

431
432
433
434
435
436
437
438
439
			case 2 :
			errno = 0;
			min_measurements = strtod(optarg, &rval);
			if ( *rval != '\0' ) {
				ERROR("Invalid value for --min-measurements.\n");
				return 1;
			}
			break;

440
441
442
			case '?' :
			break;

443
444
445
446
447
			case 0 :
			break;

			default :
			ERROR("Unhandled option '%c'\n", c);
448
449
450
451
452
453
454
455
456
457
458
			break;

		}

	}

	if ( filename == NULL ) {
		ERROR("Please specify filename using the -i option\n");
		return 1;
	}

459
460
461
462
	if ( output == NULL ) {
		output = strdup("processed.hkl");
	}

Thomas White's avatar
Thomas White committed
463
464
465
	if ( sym_str == NULL ) sym_str = strdup("1");
	sym = get_pointgroup(sym_str);
	free(sym_str);
466

467
468
469
470
471
472
473
474
475
476
477
	/* Open the data stream */
	if ( strcmp(filename, "-") == 0 ) {
		fh = stdin;
	} else {
		fh = fopen(filename, "r");
	}
	free(filename);
	if ( fh == NULL ) {
		ERROR("Failed to open input file\n");
		return 1;
	}
Thomas White's avatar
Tidy-up    
Thomas White committed
478

479
480
	/* Count the number of patterns in the file */
	n_total_patterns = count_patterns(fh);
Thomas White's avatar
Thomas White committed
481
482
483
484
	if ( n_total_patterns == 0 ) {
		ERROR("No patterns to process.\n");
		return 1;
	}
485
486
	STATUS("There are %i patterns to process\n", n_total_patterns);
	rewind(fh);
487

488
489
	model = reflist_new();

490
	if ( histo != NULL ) {
Thomas White's avatar
Thomas White committed
491

492
		int r;
Thomas White's avatar
Thomas White committed
493

494
495
496
497
498
		r = sscanf(histo, "%i,%i,%i", &hist_h, &hist_k, &hist_l);
		if ( r != 3 ) {
			ERROR("Invalid indices for '--histogram'\n");
			return 1;
		}
Thomas White's avatar
Thomas White committed
499

Thomas White's avatar
Thomas White committed
500
		space_for_hist = n_total_patterns * num_equivs(sym, NULL);
501
502
503
		hist_vals = malloc(space_for_hist * sizeof(double));
		free(histo);
		STATUS("Histogramming %i %i %i -> ", hist_h, hist_k, hist_l);
Thomas White's avatar
Thomas White committed
504

505
		/* Put into the asymmetric cell for the target group */
Thomas White's avatar
Thomas White committed
506
507
		get_asymm(sym, hist_h, hist_k, hist_l,
		          &hist_h, &hist_k, &hist_l);
508
		STATUS("%i %i %i\n", hist_h, hist_k, hist_l);
Thomas White's avatar
Thomas White committed
509

510
511
	}

512
	if ( histo_params != NULL ) {
Thomas White's avatar
Thomas White committed
513

514
		int rr;
Thomas White's avatar
Thomas White committed
515
516
517

		rr = sscanf(histo_params, "%f,%f,%i", &hist_min, &hist_max,
		                                      &hist_nbins);
518
519
520
521
522
		if ( rr != 3 ) {
			ERROR("Invalid parameters for '--hist-parameters'\n");
			return 1;
		}
		free(histo_params);
Thomas White's avatar
Thomas White committed
523
524
525
		if ( hist_max <= hist_min ) {
			ERROR("Invalid range for '--hist-parameters'. "
			      "Make sure that 'max' is greater than 'min'.\n");
526
527
			return 1;
		}
Thomas White's avatar
Thomas White committed
528

529
530
	}

531
	hist_i = 0;
532
	merge_all(fh, model, NULL, config_startafter, config_stopafter,
Thomas White's avatar
Thomas White committed
533
	          sym, n_total_patterns, hist_vals, hist_h, hist_k, hist_l,
534
	          &hist_i, config_nopolar, min_measurements);
Thomas White's avatar
Thomas White committed
535
536
537
538
	if ( ferror(fh) ) {
		ERROR("Stream read error.\n");
		return 1;
	}
539
	rewind(fh);
540

541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
	if ( config_scale ) {

		RefList *reference;

		STATUS("Extra pass for scaling...\n");

		reference = copy_reflist(model);

		reflist_free(model);
		model = reflist_new();

		rewind(fh);

		merge_all(fh, model, reference,
			  config_startafter, config_stopafter, sym,
			  n_total_patterns,
557
			  hist_vals, hist_h, hist_k, hist_l, &hist_i,
558
			  config_nopolar, min_measurements);
559
560
561
562
563
564
565
566

		if ( ferror(fh) ) {
			ERROR("Stream read error.\n");
			return 1;
		}

		reflist_free(reference);

Thomas White's avatar
Thomas White committed
567
	}
568

569
570
571
572
573
574
575
576
577
578
579
	if ( space_for_hist && (hist_i >= space_for_hist) ) {
		ERROR("Histogram array was too small!\n");
	}

	if ( hist_vals != NULL ) {
		STATUS("%i %i %i was seen %i times.\n", hist_h, hist_k, hist_l,
		                                        hist_i);
		plot_histogram(hist_vals, hist_i, hist_min, hist_max,
		               hist_nbins);
	}

580
	write_reflist(output, model);
581
582

	fclose(fh);
583

584
	free(sym);
585
	reflist_free(model);
Thomas White's avatar
Thomas White committed
586
	free(output);
587

588
589
	return 0;
}