render_hkl.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
/*
 * 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
80
static void render_za(UnitCell *cell, ReflItemList *items,
                      double *ref, unsigned int *c,
81
                      double boost, const char *sym, int wght, int colscale)
82
{
Thomas White's avatar
Thomas White committed
83
84
	cairo_surface_t *surface;
	cairo_t *dctx;
85
	double max_u, max_v, max_res, max_val;
86
	double scale_u, scale_v, scale;
Thomas White's avatar
Thomas White committed
87
	double sep_u, sep_v, max_r;
88
	double u, v;
89
	signed int max_ux, max_uy;
Thomas White's avatar
Thomas White committed
90
	double as, bs, theta;
91
92
93
	double asx, asy, asz;
	double bsx, bsy, bsz;
	double csx, csy, csz;
Thomas White's avatar
Thomas White committed
94
	float wh, ht;
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
	signed int xh, xk, xl;
	signed int yh, yk, yl;
	signed int zh, zk, zl;
	double xx, xy, xz;
	double yx, yy, yz;
	signed int xi, yi;
	char tmp[256];
	cairo_text_extents_t size;
	double cx, cy;
	const double border = 200.0;

	xh = 1;  xk = -1;  xl = 0;
	yh = 0;  yk = 0;  yl = 1;

	zh = xk*yl - xl*yk;
	zk = - xh*yl + xl*yh;
	zl = xh*yk - xk*yh;

	STATUS("Zone axis is %i %i %i\n", zh, zk, zl);
114

Thomas White's avatar
Thomas White committed
115
116
	wh = 1024;
	ht = 1024;
117
118
	cx = 512.0;
	cy = 500.0;
119

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

Thomas White's avatar
Thomas White committed
122
123
124
125
126
	if ( cairo_surface_status(surface) != CAIRO_STATUS_SUCCESS ) {
		fprintf(stderr, "Couldn't create Cairo surface\n");
		cairo_surface_destroy(surface);
		return;
	}
127

Thomas White's avatar
Thomas White committed
128
	dctx = cairo_create(surface);
129

Thomas White's avatar
Thomas White committed
130
131
132
133
	/* 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);
134

135
	max_u = 0.0;  max_v = 0.0;  max_val = 0.0;
136
	max_res = 0.0;  max_ux = 0;  max_uy = 0;
Thomas White's avatar
Thomas White committed
137
138

	/* Work out reciprocal lattice spacings and angles for this cut */
139
	if ( cell_get_reciprocal(cell, &asx, &asy, &asz,
Thomas White's avatar
Thomas White committed
140
	                          &bsx, &bsy, &bsz,
141
142
143
144
	                          &csx, &csy, &csz) ) {
		ERROR("Couldn't get reciprocal parameters\n");
		return;
	}
145
146
147
148
149
150
151
152
153
154
155
156
	xx = xh*asx + xk*bsx + xl*csx;
	xy = xh*asy + xk*bsy + xl*csy;
	xz = xh*asz + xk*bsz + xl*csz;
	yx = yh*asx + yk*bsx + yl*csx;
	yy = yh*asy + yk*bsy + yl*csy;
	yz = yh*asz + yk*bsz + yl*csz;
	theta = angle_between(xx, xy, xz, yx, yy, yz);
	as = modulus(xx, xy, xz) / 1e9;
	bs = modulus(yx, yy, yz) / 1e9;

	for ( xi=-INDMAX; xi<INDMAX; xi++ ) {
	for ( yi=-INDMAX; yi<INDMAX; yi++ ) {
Thomas White's avatar
Thomas White committed
157

158
		double u, v, val, res;
Thomas White's avatar
Thomas White committed
159
		int ct;
160
		int nequiv, p;
161
162
163
164
165
		signed int h, k, l;

		h = xi*xh + yi*yh;
		k = xi*xk + yi*yk;
		l = xi*xl + yi*yl;
Thomas White's avatar
Thomas White committed
166

167
168
169
170
171
172
173
		/* Do we have this reflection? */
		if ( find_unique_equiv(items, h, k, l, sym, &h, &k, &l) == 0 ) {
			continue;
		}

		/* Reflection in the zone? */
		if ( zh*h + zk*k + zl*l != 0 ) continue;
Thomas White's avatar
Thomas White committed
174

175
176
		switch ( wght ) {
		case WGHT_I :
177
			val = lookup_intensity(ref, h, k, l);
178
179
			break;
		case WGHT_SQRTI :
180
			val = lookup_intensity(ref, h, k, l);
181
182
183
			val = (val>0.0) ? sqrt(val) : 0.0;
			break;
		case WGHT_COUNTS :
184
			val = (float)ct;
185
			val /= (float)num_equivs(h, k, l, sym);
186
187
			break;
		case WGHT_RAWCOUNTS :
188
189
			val = (float)ct;
			break;
Thomas White's avatar
Thomas White committed
190
191
192
		default :
			ERROR("Invalid weighting.\n");
			abort();
Thomas White's avatar
Thomas White committed
193
		}
194

195
		nequiv = num_equivs(h, k, l, sym);
196
197
198
		for ( p=0; p<nequiv; p++ ) {

			signed int he, ke, le;
199
200
201
202
203
204
205
206
207
			signed int ux, uy;

			get_equiv(h, k, l, &he, &ke, &le, sym, p);

			ux = he*xh + ke*xk + le*xl;
			uy = he*yh + ke*yk + le*yl;

			u = (double)ux*as*sin(theta);
			v = (double)ux*as*cos(theta) + uy*bs;
208
209
210
211
212

			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);
213
			}
214
215
216
			if ( fabs(ux) > max_ux ) max_ux = fabs(ux);
			if ( fabs(uy) > max_uy ) max_uy = fabs(uy);
			res = resolution(cell, he, ke, le);
217
218
			if ( res > max_res ) max_res = res;

219
220
221
		}

	}
Thomas White's avatar
Thomas White committed
222
	}
223

Thomas White's avatar
Thomas White committed
224
	max_res /= 1e9;
225
226
	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
227

228
	if ( max_val <= 0.0 ) {
Thomas White's avatar
Thomas White committed
229
		max_r = 4.0;
230
		STATUS("Couldn't find max value.\n");
Thomas White's avatar
Thomas White committed
231
		goto out;
232
233
	}

Thomas White's avatar
Thomas White committed
234
	/* Choose whichever scaling factor gives the smallest value */
235
236
	scale_u = ((double)wh-border) / (2.0*max_u);
	scale_v = ((double)ht-border) / (2.0*max_v);
237
	scale = (scale_u < scale_v) ? scale_u : scale_v;
238

239
240
241
	sep_u = (double)scale*as*sin(theta);
	sep_v = (double)scale*as*cos(theta) + scale*bs;
	/* We are interested in the smaller of the two separations */
242
	max_r = (sep_u < sep_v) ? sep_u : sep_v;
243
	max_r /= 2.0;  /* Max radius if half the separation */
244
	max_r -= 1.0;  /* Add a tiny separation between circles */
245
246
247
	if ( max_r < 1.0 ) {
		ERROR("Circle radius is probably too small (%f).\n", max_r);
	}
Thomas White's avatar
Thomas White committed
248

249
250
	for ( xi=-INDMAX; xi<INDMAX; xi++ ) {
	for ( yi=-INDMAX; yi<INDMAX; yi++ ) {
Thomas White's avatar
Thomas White committed
251

252
		double u, v, val;
Thomas White's avatar
Thomas White committed
253
		int ct;
254
		int nequiv, p;
255
		signed int h, k, l;
Thomas White's avatar
Thomas White committed
256

257
258
259
260
261
262
263
264
265
266
267
268
		h = xi*xh + yi*yh;
		k = xi*xk + yi*yk;
		l = xi*xl + yi*yl;

		/* Do we have this reflection? */
		if ( find_unique_equiv(items, h, k, l, sym, &h, &k, &l) == 0 ) {
			continue;
		}


		/* Reflection in the zone? */
		if ( zh*h + zk*k + zl*l != 0 ) continue;
Thomas White's avatar
Thomas White committed
269

270
271
		switch ( wght ) {
		case WGHT_I :
272
			val = lookup_intensity(ref, h, k, l);
273
274
			break;
		case WGHT_SQRTI :
275
			val = lookup_intensity(ref, h, k, l);
276
277
278
			val = (val>0.0) ? sqrt(val) : 0.0;
			break;
		case WGHT_COUNTS :
279
			val = (float)ct;
280
			val /= (float)num_equivs(h, k, l, sym);
281
282
			break;
		case WGHT_RAWCOUNTS :
283
284
			val = (float)ct;
			break;
Thomas White's avatar
Thomas White committed
285
286
287
		default :
			ERROR("Invalid weighting.\n");
			abort();
Thomas White's avatar
Thomas White committed
288
		}
289

290
		nequiv = num_equivs(h, k, l, sym);
291
292
293
		for ( p=0; p<nequiv; p++ ) {

			signed int he, ke, le;
294
			signed int ux, uy;
295
			float r, g, b;
296
297
298
299
			get_equiv(h, k, l, &he, &ke, &le, sym, p);

			ux = he*xh + ke*xk + le*xl;
			uy = he*yh + ke*yk + le*yl;
300

301
302
			u = (double)ux*as*sin(theta);
			v = (double)ux*as*cos(theta) + uy*bs;
Thomas White's avatar
Thomas White committed
303

304
305
306
			cairo_arc(dctx, ((double)cx)+u*scale,
			                ((double)cy)+v*scale,
			                max_r, 0, 2*M_PI);
Thomas White's avatar
Thomas White committed
307

308
309
			render_scale(val, max_val/boost, colscale, &r, &g, &b);
			cairo_set_source_rgb(dctx, r, g, b);
310
311
312
			cairo_fill(dctx);

		}
Thomas White's avatar
Thomas White committed
313
314

	}
315
316
	}

Thomas White's avatar
Thomas White committed
317
318
out:
	/* Centre marker */
319
320
	cairo_arc(dctx, (double)cx,
			(double)cy, max_r, 0, 2*M_PI);
Thomas White's avatar
Thomas White committed
321
322
323
	cairo_set_source_rgb(dctx, 1.0, 0.0, 0.0);
	cairo_fill(dctx);

324
325
	/* Draw indexing lines */
	cairo_set_line_width(dctx, 4.0);
326
327
328
329
	cairo_move_to(dctx, (double)cx, (double)cy);
	u = (double)max_ux*as*sin(theta);
	v = (double)max_ux*as*cos(theta);
	cairo_line_to(dctx, cx+u*scale, cy+v*scale);
330
331
332
	cairo_set_source_rgb(dctx, 0.0, 1.0, 0.0);
	cairo_stroke(dctx);

333
334
335
336
337
338
	cairo_set_font_size(dctx, 40.0);
	snprintf(tmp, 255, "%i%i%i", xh, xk, xl);
	cairo_text_extents(dctx, tmp, &size);

	cairo_move_to(dctx, cx+u*scale + 20.0, cy+v*scale + size.height/2.0);
	cairo_show_text(dctx, tmp);
339
340
	cairo_fill(dctx);

341
342
343
344
	snprintf(tmp, 255, "%i%i%i", yh, yk, yl);
	cairo_text_extents(dctx, tmp, &size);

	cairo_move_to(dctx, (double)cx, (double)cy);
345
	u = 0.0;
346
347
	v = (double)max_uy*bs;
	cairo_line_to(dctx, cx+u*scale, cy+v*scale);
348
349
350
	cairo_set_source_rgb(dctx, 0.0, 1.0, 0.0);
	cairo_stroke(dctx);

351
352
353
	cairo_move_to(dctx, cx+u*scale - size.width/2.0,
	                    cy+v*scale + size.height + 20.0);
	cairo_show_text(dctx, tmp);
354
355
	cairo_fill(dctx);

Thomas White's avatar
Thomas White committed
356
357
358
	cairo_surface_finish(surface);
	cairo_destroy(dctx);
}
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397

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;
}
398
#endif
Thomas White's avatar
Thomas White committed
399
400
401
402
403
404
405
406
407
408


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
409
	int config_sqrt = 0;
410
	int config_colkey = 0;
Thomas White's avatar
Thomas White committed
411
	unsigned int nproc = 1;
412
	char *pdb = NULL;
413
	int r = 0;
414
	double boost = 1.0;
415
	char *sym = NULL;
416
417
	char *weighting = NULL;
	int wght;
418
419
	int colscale;
	char *cscale = NULL;
Thomas White's avatar
Thomas White committed
420
	unsigned int *cts;
Thomas White's avatar
Thomas White committed
421
422
423
424
425
426

	/* Long options */
	const struct option longopts[] = {
		{"help",               0, NULL,               'h'},
		{"povray",             0, &config_povray,      1},
		{"zone-axis",          0, &config_zoneaxis,    1},
427
		{"pdb",                1, NULL,               'p'},
428
		{"boost",              1, NULL,               'b'},
429
		{"symmetry",           1, NULL,               'y'},
430
		{"weighting",          1, NULL,               'w'},
431
		{"colscale",           1, NULL,               'c'},
432
		{"counts",             0, &config_sqrt,        1},
433
		{"colour-key",         0, &config_colkey,      1},
Thomas White's avatar
Thomas White committed
434
435
436
437
		{0, 0, NULL, 0}
	};

	/* Short options */
438
439
	while ((c = getopt_long(argc, argv, "hj:p:w:c:y:",
	                        longopts, NULL)) != -1) {
Thomas White's avatar
Thomas White committed
440
441

		switch (c) {
Thomas White's avatar
Thomas White committed
442
		case 'h' :
Thomas White's avatar
Thomas White committed
443
444
445
			show_help(argv[0]);
			return 0;

Thomas White's avatar
Thomas White committed
446
		case 'j' :
Thomas White's avatar
Thomas White committed
447
448
449
			nproc = atoi(optarg);
			break;

450
451
452
453
		case 'p' :
			pdb = strdup(optarg);
			break;

454
455
456
457
		case 'b' :
			boost = atof(optarg);
			break;

458
459
460
461
		case 'y' :
			sym = strdup(optarg);
			break;

462
463
464
465
		case 'w' :
			weighting = strdup(optarg);
			break;

466
467
468
469
		case 'c' :
			cscale = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
470
		case 0 :
Thomas White's avatar
Thomas White committed
471
472
			break;

Thomas White's avatar
Thomas White committed
473
		default :
Thomas White's avatar
Thomas White committed
474
475
476
477
478
			return 1;
		}

	}

479
480
481
482
	if ( pdb == NULL ) {
		pdb = strdup("molecule.pdb");
	}

483
484
485
486
	if ( sym == NULL ) {
		sym = strdup("1");
	}

487
488
489
490
491
492
493
494
495
496
	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;
497
498
	} else if ( strcmp(weighting, "counts") == 0 ) {
		wght = WGHT_COUNTS;
499
500
	} else if ( strcmp(weighting, "rawcts") == 0 ) {
		wght = WGHT_RAWCOUNTS;
501
502
503
504
	} else if ( strcmp(weighting, "rawcount") == 0 ) {
		wght = WGHT_RAWCOUNTS;
	} else if ( strcmp(weighting, "rawcounts") == 0 ) {
		wght = WGHT_RAWCOUNTS;
505
506
507
508
	} else {
		ERROR("Unrecognised weighting '%s'\n", weighting);
		return 1;
	}
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
	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);
528

529
530
531
532
	if ( config_colkey ) {
		return render_key(colscale);
	}

Thomas White's avatar
Thomas White committed
533
534
	infile = argv[optind];

535
	cell = load_cell_from_pdb(pdb);
536
537
538
539
	if ( cell == NULL ) {
		ERROR("Couldn't load unit cell from %s\n", pdb);
		return 1;
	}
Thomas White's avatar
Thomas White committed
540
	ref = new_list_intensity();
Thomas White's avatar
Thomas White committed
541
	cts = new_list_count();
Thomas White's avatar
Thomas White committed
542
	ReflItemList *items = read_reflections(infile, ref, NULL, cts);
Thomas White's avatar
Thomas White committed
543
544
545
546
547
548
	if ( ref == NULL ) {
		ERROR("Couldn't open file '%s'\n", infile);
		return 1;
	}

	if ( config_povray ) {
549
		r = povray_render_animation(cell, ref, cts, nproc);
Thomas White's avatar
Thomas White committed
550
	} else if ( config_zoneaxis ) {
551
#ifdef HAVE_CAIRO
552
		render_za(cell, items, ref, cts, boost, sym, wght, colscale);
553
554
555
556
557
#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
558
559
560
561
	} else {
		ERROR("Try again with either --povray or --zone-axis.\n");
	}

562
	free(pdb);
563
	free(sym);
564
	delete_items(items);
565

566
	return r;
567
}