render_hkl.c 11.7 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
		max_r = 4.0;
187
		STATUS("Couldn't find max value.\n");
Thomas White's avatar
Thomas White committed
188
		goto out;
189
190
	}

Thomas White's avatar
Thomas White committed
191
	/* Choose whichever scaling factor gives the smallest value */
192
193
194
	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;
195

Thomas White's avatar
Thomas White committed
196
197
	sep_u = as * scale * cos(theta);
	sep_v = bs * scale;
198
199
	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
200
201
202
203

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

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

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

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

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

			signed int he, ke, le;
235
			float r, g, b;
236
237
238
239
			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
240

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

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

		}
Thomas White's avatar
Thomas White committed
250
251

	}
252
253
	}

Thomas White's avatar
Thomas White committed
254
255
256
257
258
259
260
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);

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
290
	/* 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
291
292
293
	cairo_surface_finish(surface);
	cairo_destroy(dctx);
}
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
332

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;
}
333
#endif
Thomas White's avatar
Thomas White committed
334
335
336
337
338
339
340
341
342
343


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
344
	int config_sqrt = 0;
345
	int config_colkey = 0;
Thomas White's avatar
Thomas White committed
346
	unsigned int nproc = 1;
347
	char *pdb = NULL;
348
	int r = 0;
349
	double boost = 1.0;
350
	char *sym = NULL;
351
352
	char *weighting = NULL;
	int wght;
353
354
	int colscale;
	char *cscale = NULL;
Thomas White's avatar
Thomas White committed
355
	unsigned int *cts;
Thomas White's avatar
Thomas White committed
356
357
358
359
360
361

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

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

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

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

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

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

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

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

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

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

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

	}

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

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

422
423
424
425
426
427
428
429
430
431
	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;
432
433
	} else if ( strcmp(weighting, "counts") == 0 ) {
		wght = WGHT_COUNTS;
434
435
	} else if ( strcmp(weighting, "rawcts") == 0 ) {
		wght = WGHT_RAWCOUNTS;
436
437
438
439
	} else if ( strcmp(weighting, "rawcount") == 0 ) {
		wght = WGHT_RAWCOUNTS;
	} else if ( strcmp(weighting, "rawcounts") == 0 ) {
		wght = WGHT_RAWCOUNTS;
440
441
442
443
	} else {
		ERROR("Unrecognised weighting '%s'\n", weighting);
		return 1;
	}
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
	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);
463

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

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

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

	if ( config_povray ) {
485
		r = povray_render_animation(cell, ref, cts, nproc);
Thomas White's avatar
Thomas White committed
486
	} else if ( config_zoneaxis ) {
487
#ifdef HAVE_CAIRO
488
		render_za(cell, ref, cts, boost, sym, wght, colscale);
489
490
491
492
493
#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
494
495
496
497
	} else {
		ERROR("Try again with either --povray or --zone-axis.\n");
	}

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

501
	return r;
502
}