render_hkl.c 11.6 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
/*
 * render_hkl.c
 *
 * Draw pretty renderings of reflection lists
 *
 * (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>
23
#ifdef HAVE_CAIRO
Thomas White's avatar
Thomas White committed
24
25
#include <cairo.h>
#include <cairo-pdf.h>
26
#endif
27
28
29

#include "utils.h"
#include "reflections.h"
30
#include "povray.h"
31
#include "symmetry.h"
32
#include "render.h"
33

34
35
36
37
enum {
	WGHT_I,
	WGHT_SQRTI,
	WGHT_COUNTS,
38
	WGHT_RAWCOUNTS,
39
};
40

41
42
43
44
static void show_help(const char *s)
{
	printf("Syntax: %s [options] <file.hkl>\n\n", s);
	printf(
45
"Render intensity lists in various ways.\n"
46
"\n"
47
"      --povray            Render a 3D animation using POV-ray.\n"
48
#ifdef HAVE_CAIRO
49
"      --zone-axis         Render a 2D zone axis pattern.\n"
50
#endif
51
"\n"
52
53
54
"      --boost=<val>       Squash colour scale by <val>.\n"
"  -p, --pdb=<file>        PDB file from which to get the unit cell.\n"
"  -y, --symmetry=<sym>    Expand reflections according to point group <sym>.\n"
55
56
57
58
59
60
61
"\n"
"  -c, --colscale=<scale>  Use the given colour scale.  Choose from:\n"
"                           mono    : Greyscale, black is zero.\n"
"                           invmono : Greyscale, white is zero.\n"
"                           colour  : Colours scale:\n"
"                                     black-blue-pink-red-orange-yellow-white\n"
"\n"
62
63
"  -w  --weighting=<wght>  Colour/shade the reciprocal lattice points\n"
"                           according to:\n"
64
65
66
67
68
69
"                            I      : the intensity of the reflection.\n"
"                            sqrtI  : the square root of the intensity.\n"
"                            count  : the number of hits for the reflection.\n"
"                                     (after correcting for 'epsilon')\n"
"                            rawcts : the raw number of hits for the\n"
"                                     reflection (no 'epsilon' correction).\n"
70
71
72
73
"\n"
"      --colour-key        Draw (only) the key for the current colour scale.\n"
"  -j <n>                  Run <n> instances of POV-ray in parallel.\n"
"  -h, --help              Display this help message.\n"
74
);
75
76
77
}


78
#ifdef HAVE_CAIRO
79
static void render_za(UnitCell *cell, double *ref, unsigned int *c,
80
                      double boost, const char *sym, int wght, int colscale)
81
{
Thomas White's avatar
Thomas White committed
82
83
	cairo_surface_t *surface;
	cairo_t *dctx;
84
	double max_u, max_v, max_res, max_val;
85
	double scale_u, scale_v, scale;
Thomas White's avatar
Thomas White committed
86
	double sep_u, sep_v, max_r;
87
88
	double u, v;
	signed int max_h, max_k;
Thomas White's avatar
Thomas White committed
89
	double as, bs, theta;
90
91
92
	double asx, asy, asz;
	double bsx, bsy, bsz;
	double csx, csy, csz;
Thomas White's avatar
Thomas White committed
93
94
	signed int h, k;
	float wh, ht;
95

Thomas White's avatar
Thomas White committed
96
97
	wh = 1024;
	ht = 1024;
98

Thomas White's avatar
Thomas White committed
99
	surface = cairo_pdf_surface_create("za.pdf", wh, ht);
100

Thomas White's avatar
Thomas White committed
101
102
103
104
105
	if ( cairo_surface_status(surface) != CAIRO_STATUS_SUCCESS ) {
		fprintf(stderr, "Couldn't create Cairo surface\n");
		cairo_surface_destroy(surface);
		return;
	}
106

Thomas White's avatar
Thomas White committed
107
	dctx = cairo_create(surface);
108

Thomas White's avatar
Thomas White committed
109
110
111
112
	/* Black background */
	cairo_rectangle(dctx, 0.0, 0.0, wh, ht);
	cairo_set_source_rgb(dctx, 0.0, 0.0, 0.0);
	cairo_fill(dctx);
113

114
	max_u = 0.0;  max_v = 0.0;  max_val = 0.0;
115
	max_res = 0.0;  max_h = 0;  max_k = 0;
Thomas White's avatar
Thomas White committed
116
117

	/* Work out reciprocal lattice spacings and angles for this cut */
118
	if ( cell_get_reciprocal(cell, &asx, &asy, &asz,
Thomas White's avatar
Thomas White committed
119
	                          &bsx, &bsy, &bsz,
120
121
122
123
	                          &csx, &csy, &csz) ) {
		ERROR("Couldn't get reciprocal parameters\n");
		return;
	}
Thomas White's avatar
Thomas White committed
124
125
126
127
128
129
130
	theta = angle_between(asx, asy, asz, bsx, bsy, bsz);
	as = modulus(asx, asy, asz) / 1e9;
	bs = modulus(bsx, bsy, bsz) / 1e9;

	for ( h=-INDMAX; h<INDMAX; h++ ) {
	for ( k=-INDMAX; k<INDMAX; k++ ) {

131
		double u, v, val, res;
Thomas White's avatar
Thomas White committed
132
		int ct;
133
		int nequiv, p;
Thomas White's avatar
Thomas White committed
134
135
136
137

		ct = lookup_count(c, h, k, 0);
		if ( ct < 1 ) continue;

138
139
		switch ( wght ) {
		case WGHT_I :
Thomas White's avatar
Thomas White committed
140
			val = lookup_intensity(ref, h, k, 0);
141
142
			break;
		case WGHT_SQRTI :
Thomas White's avatar
Thomas White committed
143
			val = lookup_intensity(ref, h, k, 0);
144
145
146
			val = (val>0.0) ? sqrt(val) : 0.0;
			break;
		case WGHT_COUNTS :
147
148
149
150
			val = (float)ct;
			val /= (float)num_equivs(h, k, 0, sym);
			break;
		case WGHT_RAWCOUNTS :
151
152
			val = (float)ct;
			break;
Thomas White's avatar
Thomas White committed
153
154
155
		default :
			ERROR("Invalid weighting.\n");
			abort();
Thomas White's avatar
Thomas White committed
156
		}
157

158
159
160
161
162
163
164
165
166
167
168
169
		nequiv = num_equivs(h, k, 0, sym);
		for ( p=0; p<nequiv; p++ ) {

			signed int he, ke, le;
			get_equiv(h, k, 0, &he, &ke, &le, sym, p);

			u =  (double)he*as*sin(theta);
			v =  (double)he*as*cos(theta) + ke*bs;
			if ( fabs(u) > fabs(max_u) ) max_u = fabs(u);
			if ( fabs(v) > fabs(max_v) ) max_v = fabs(v);
			if ( fabs(val) > fabs(max_val) ) {
				max_val = fabs(val);
170
			}
171
172
173
174
175
			if ( fabs(he) > max_h ) max_h = fabs(he);
			if ( fabs(ke) > max_k ) max_k = fabs(ke);
			res = resolution(cell, he, ke, 0);
			if ( res > max_res ) max_res = res;

176
177
178
		}

	}
Thomas White's avatar
Thomas White committed
179
	}
180

Thomas White's avatar
Thomas White committed
181
	max_res /= 1e9;
182
183
	printf("Maximum resolution is 1/d = %5.3f nm^-1, d = %5.3f nm\n",
	       max_res*2.0, 1.0/(max_res*2.0));
Thomas White's avatar
Thomas White committed
184

185
	if ( max_val <= 0.0 ) {
Thomas White's avatar
Thomas White committed
186
187
		max_r = 4.0;
		goto out;
188
189
	}

Thomas White's avatar
Thomas White committed
190
	/* Choose whichever scaling factor gives the smallest value */
191
192
193
	scale_u = ((double)(wh/2.0)-50.0) / max_u;
	scale_v = ((double)(ht/2.0)-50.0) / max_v;
	scale = (scale_u < scale_v) ? scale_u : scale_v;
194

Thomas White's avatar
Thomas White committed
195
196
	sep_u = as * scale * cos(theta);
	sep_v = bs * scale;
197
198
	max_r = (sep_u < sep_v) ? sep_u : sep_v;
	max_r -= 1.0;  /* Add a tiny separation between circles */
Thomas White's avatar
Thomas White committed
199
200
201
202

	for ( h=-INDMAX; h<INDMAX; h++ ) {
	for ( k=-INDMAX; k<INDMAX; k++ ) {

203
		double u, v, val;
Thomas White's avatar
Thomas White committed
204
		int ct;
205
		int nequiv, p;
Thomas White's avatar
Thomas White committed
206
207

		ct = lookup_count(c, h, k, 0);
208
		if ( ct < 1 ) continue;  /* Must have at least one count */
Thomas White's avatar
Thomas White committed
209

210
211
		switch ( wght ) {
		case WGHT_I :
212
			val = lookup_intensity(ref, h, k, 0);
213
214
			break;
		case WGHT_SQRTI :
215
			val = lookup_intensity(ref, h, k, 0);
216
217
218
			val = (val>0.0) ? sqrt(val) : 0.0;
			break;
		case WGHT_COUNTS :
219
220
221
222
			val = (float)ct;
			val /= (float)num_equivs(h, k, 0, sym);
			break;
		case WGHT_RAWCOUNTS :
223
224
			val = (float)ct;
			break;
Thomas White's avatar
Thomas White committed
225
226
227
		default :
			ERROR("Invalid weighting.\n");
			abort();
Thomas White's avatar
Thomas White committed
228
		}
229

230
231
232
233
		nequiv = num_equivs(h, k, 0, sym);
		for ( p=0; p<nequiv; p++ ) {

			signed int he, ke, le;
234
			float r, g, b;
235
236
237
238
			get_equiv(h, k, 0, &he, &ke, &le, sym, p);

			u = (double)he*as*sin(theta);
			v = (double)he*as*cos(theta) + ke*bs;
Thomas White's avatar
Thomas White committed
239

240
241
			cairo_arc(dctx, ((double)wh/2)+u*scale,
					((double)ht/2)+v*scale, max_r,
242
					0, 2*M_PI);
Thomas White's avatar
Thomas White committed
243

244
245
			render_scale(val, max_val/boost, colscale, &r, &g, &b);
			cairo_set_source_rgb(dctx, r, g, b);
246
247
248
			cairo_fill(dctx);

		}
Thomas White's avatar
Thomas White committed
249
250

	}
251
252
	}

Thomas White's avatar
Thomas White committed
253
254
255
256
257
258
259
out:
	/* Centre marker */
	cairo_arc(dctx, (double)wh/2,
			(double)ht/2, max_r, 0, 2*M_PI);
	cairo_set_source_rgb(dctx, 1.0, 0.0, 0.0);
	cairo_fill(dctx);

260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
	/* Draw indexing lines */
	cairo_set_line_width(dctx, 4.0);
	cairo_move_to(dctx, (double)wh/2, (double)ht/2);
	u = (double)max_h*as*sin(theta);
	v = (double)max_h*as*cos(theta) + 0*bs;
	cairo_line_to(dctx, ((double)wh/2)+u*scale,
	                    ((double)ht/2)+v*scale);
	cairo_set_source_rgb(dctx, 0.0, 1.0, 0.0);
	cairo_stroke(dctx);

	cairo_move_to(dctx,((double)wh/2)+u*scale-40.0,
	                   ((double)ht/2)+v*scale+40.0);
	                   cairo_set_font_size(dctx, 40.0);
	cairo_show_text(dctx, "h");
	cairo_fill(dctx);

	cairo_move_to(dctx, (double)wh/2, (double)ht/2);
	u = 0.0;
	v = (double)max_k*bs;
	cairo_line_to(dctx, ((double)wh/2)+u*scale,
	                    ((double)ht/2)+v*scale);
	cairo_set_source_rgb(dctx, 0.0, 1.0, 0.0);
	cairo_stroke(dctx);

	cairo_move_to(dctx,((double)wh/2)+u*scale-40.0,
	                   ((double)ht/2)+v*scale-40.0);
	                   cairo_set_font_size(dctx, 40.0);
	cairo_show_text(dctx, "k");
	cairo_fill(dctx);

Thomas White's avatar
Thomas White committed
290
291
292
	cairo_surface_finish(surface);
	cairo_destroy(dctx);
}
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331

static int render_key(int colscale)
{
	cairo_surface_t *surface;
	cairo_t *dctx;
	float wh, ht;
	float y;

	wh = 128;
	ht = 1024;

	surface = cairo_pdf_surface_create("key.pdf", wh, ht);

	if ( cairo_surface_status(surface) != CAIRO_STATUS_SUCCESS ) {
		fprintf(stderr, "Couldn't create Cairo surface\n");
		cairo_surface_destroy(surface);
		return 1;
	}

	dctx = cairo_create(surface);

	for ( y=0; y<ht; y++ ) {

		float r, g, b;

		cairo_rectangle(dctx, 0.0, y, wh, y+1.0);

		render_scale(ht-y, ht, colscale, &r, &g, &b);
		cairo_set_source_rgb(dctx, r, g, b);

		cairo_fill(dctx);

	}

	cairo_surface_finish(surface);
	cairo_destroy(dctx);

	return 0;
}
332
#endif
Thomas White's avatar
Thomas White committed
333
334
335
336
337
338
339
340
341
342


int main(int argc, char *argv[])
{
	int c;
	UnitCell *cell;
	char *infile;
	double *ref;
	int config_povray = 0;
	int config_zoneaxis = 0;
Thomas White's avatar
Thomas White committed
343
	int config_sqrt = 0;
344
	int config_colkey = 0;
Thomas White's avatar
Thomas White committed
345
	unsigned int nproc = 1;
346
	char *pdb = NULL;
347
	int r = 0;
348
	double boost = 1.0;
349
	char *sym = NULL;
350
351
	char *weighting = NULL;
	int wght;
352
353
	int colscale;
	char *cscale = NULL;
Thomas White's avatar
Thomas White committed
354
	unsigned int *cts;
Thomas White's avatar
Thomas White committed
355
356
357
358
359
360

	/* Long options */
	const struct option longopts[] = {
		{"help",               0, NULL,               'h'},
		{"povray",             0, &config_povray,      1},
		{"zone-axis",          0, &config_zoneaxis,    1},
361
		{"pdb",                1, NULL,               'p'},
362
		{"boost",              1, NULL,               'b'},
363
		{"symmetry",           1, NULL,               'y'},
364
		{"weighting",          1, NULL,               'w'},
365
		{"colscale",           1, NULL,               'c'},
366
		{"counts",             0, &config_sqrt,        1},
367
		{"colour-key",         0, &config_colkey,      1},
Thomas White's avatar
Thomas White committed
368
369
370
371
		{0, 0, NULL, 0}
	};

	/* Short options */
372
373
	while ((c = getopt_long(argc, argv, "hj:p:w:c:y:",
	                        longopts, NULL)) != -1) {
Thomas White's avatar
Thomas White committed
374
375

		switch (c) {
Thomas White's avatar
Thomas White committed
376
		case 'h' :
Thomas White's avatar
Thomas White committed
377
378
379
			show_help(argv[0]);
			return 0;

Thomas White's avatar
Thomas White committed
380
		case 'j' :
Thomas White's avatar
Thomas White committed
381
382
383
			nproc = atoi(optarg);
			break;

384
385
386
387
		case 'p' :
			pdb = strdup(optarg);
			break;

388
389
390
391
		case 'b' :
			boost = atof(optarg);
			break;

392
393
394
395
		case 'y' :
			sym = strdup(optarg);
			break;

396
397
398
399
		case 'w' :
			weighting = strdup(optarg);
			break;

400
401
402
403
		case 'c' :
			cscale = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
404
		case 0 :
Thomas White's avatar
Thomas White committed
405
406
			break;

Thomas White's avatar
Thomas White committed
407
		default :
Thomas White's avatar
Thomas White committed
408
409
410
411
412
			return 1;
		}

	}

413
414
415
416
	if ( pdb == NULL ) {
		pdb = strdup("molecule.pdb");
	}

417
418
419
420
	if ( sym == NULL ) {
		sym = strdup("1");
	}

421
422
423
424
425
426
427
428
429
430
	if ( weighting == NULL ) {
		weighting = strdup("I");
	}

	if ( strcmp(weighting, "I") == 0 ) {
		wght = WGHT_I;
	} else if ( strcmp(weighting, "sqrtI") == 0 ) {
		wght = WGHT_SQRTI;
	} else if ( strcmp(weighting, "count") == 0 ) {
		wght = WGHT_COUNTS;
431
432
	} else if ( strcmp(weighting, "counts") == 0 ) {
		wght = WGHT_COUNTS;
433
434
	} else if ( strcmp(weighting, "rawcts") == 0 ) {
		wght = WGHT_RAWCOUNTS;
435
436
437
438
	} else if ( strcmp(weighting, "rawcount") == 0 ) {
		wght = WGHT_RAWCOUNTS;
	} else if ( strcmp(weighting, "rawcounts") == 0 ) {
		wght = WGHT_RAWCOUNTS;
439
440
441
442
	} else {
		ERROR("Unrecognised weighting '%s'\n", weighting);
		return 1;
	}
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
	free(weighting);

	if ( cscale == NULL ) {
		cscale = strdup("mono");
	}

	if ( strcmp(cscale, "mono") == 0 ) {
		colscale = SCALE_MONO;
	} else if ( strcmp(cscale, "invmono") == 0 ) {
		colscale = SCALE_INVMONO;
	} else if ( strcmp(cscale, "colour") == 0 ) {
		colscale = SCALE_COLOUR;
	} else if ( strcmp(cscale, "color") == 0 ) {
		colscale = SCALE_COLOUR;
	} else {
		ERROR("Unrecognised colour scale '%s'\n", cscale);
		return 1;
	}
	free(cscale);
462

463
464
465
466
	if ( config_colkey ) {
		return render_key(colscale);
	}

Thomas White's avatar
Thomas White committed
467
468
	infile = argv[optind];

469
	cell = load_cell_from_pdb(pdb);
470
471
472
473
	if ( cell == NULL ) {
		ERROR("Couldn't load unit cell from %s\n", pdb);
		return 1;
	}
Thomas White's avatar
Thomas White committed
474
	ref = new_list_intensity();
Thomas White's avatar
Thomas White committed
475
	cts = new_list_count();
Thomas White's avatar
Thomas White committed
476
477
	ReflItemList *items = read_reflections(infile, ref, NULL, cts);
	delete_items(items);
Thomas White's avatar
Thomas White committed
478
479
480
481
482
483
	if ( ref == NULL ) {
		ERROR("Couldn't open file '%s'\n", infile);
		return 1;
	}

	if ( config_povray ) {
484
		r = povray_render_animation(cell, ref, cts, nproc);
Thomas White's avatar
Thomas White committed
485
	} else if ( config_zoneaxis ) {
486
#ifdef HAVE_CAIRO
487
		render_za(cell, ref, cts, boost, sym, wght, colscale);
488
489
490
491
492
#else
		ERROR("This version of CrystFEL was compiled without Cairo");
		ERROR(" support, which is required to plot a zone axis");
		ERROR(" pattern.  Sorry!\n");
#endif
Thomas White's avatar
Thomas White committed
493
494
495
496
	} else {
		ERROR("Try again with either --povray or --zone-axis.\n");
	}

497
	free(pdb);
498
	free(sym);
499

500
	return r;
501
}