render_hkl.c 11.5 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
23
24
/*
 * 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>
#include <sys/types.h>
#include <sys/wait.h>
Thomas White's avatar
Thomas White committed
25
26
#include <cairo.h>
#include <cairo-pdf.h>
27
28
29
30
31

#include "utils.h"
#include "reflections.h"


32
33
34
#define MAX_PROC (256)


35
36
37
38
static void show_help(const char *s)
{
	printf("Syntax: %s [options] <file.hkl>\n\n", s);
	printf(
39
"Render intensity lists using POV-ray.\n"
40
"\n"
41
42
43
44
45
46
"  -h, --help       Display this help message.\n"
"      --povray     Render a 3D animation using POV-ray.\n"
"      --zone-axis  Render a 2D zone axis pattern.\n"
"  -j <n>           Run <n> instances of POV-ray in parallel.\n"
"  -p, --pdb=<file> PDB file from which to get the unit cell.\n"
);
47
48
49
}


Thomas White's avatar
Thomas White committed
50
static void render_za(UnitCell *cell, double *ref, unsigned int *c)
51
{
Thomas White's avatar
Thomas White committed
52
53
54
55
56
	cairo_surface_t *surface;
	cairo_t *dctx;
	double max_u, max_v, max_res, max_intensity, scale;
	double sep_u, sep_v, max_r;
	double as, bs, theta;
57
58
59
	double asx, asy, asz;
	double bsx, bsy, bsz;
	double csx, csy, csz;
Thomas White's avatar
Thomas White committed
60
61
	signed int h, k;
	float wh, ht;
62

Thomas White's avatar
Thomas White committed
63
64
	wh = 1024;
	ht = 1024;
65

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

Thomas White's avatar
Thomas White committed
68
69
70
71
72
	if ( cairo_surface_status(surface) != CAIRO_STATUS_SUCCESS ) {
		fprintf(stderr, "Couldn't create Cairo surface\n");
		cairo_surface_destroy(surface);
		return;
	}
73

Thomas White's avatar
Thomas White committed
74
	dctx = cairo_create(surface);
75

Thomas White's avatar
Thomas White committed
76
77
78
79
	/* 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);
80

Thomas White's avatar
Thomas White committed
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
	max_u = 0.0;  max_v = 0.0;  max_intensity = 0.0;
	max_res = 0.0;

	/* Work out reciprocal lattice spacings and angles for this cut */
	cell_get_reciprocal(cell, &asx, &asy, &asz,
	                          &bsx, &bsy, &bsz,
	                          &csx, &csy, &csz);
	theta = angle_between(asx, asy, asz, bsx, bsy, bsz);
	as = modulus(asx, asy, asz) / 1e9;
	bs = modulus(bsx, bsy, bsz) / 1e9;
	STATUS("theta=%f\n", rad2deg(theta));

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

		double u, v, intensity, res;
		int ct;

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

		intensity = lookup_intensity(ref, h, k, 0) / (float)ct;

		res = resolution(cell, h, k, 0);
		if ( res > max_res ) max_res = res;

		if ( intensity != 0 ) {
			u =  (double)h*as*sin(theta);
			v =  (double)h*as*cos(theta) + k*bs;
			if ( fabs(u) > fabs(max_u) ) max_u = fabs(u);
			if ( fabs(v) > fabs(max_v) ) max_v = fabs(v);
			if ( fabs(intensity) > fabs(max_intensity) )
						max_intensity = fabs(intensity);
114
115
116
		}

	}
Thomas White's avatar
Thomas White committed
117
	}
118

Thomas White's avatar
Thomas White committed
119
120
121
122
123
124
125
126
	max_res /= 1e9;
	max_u /= 0.5;
	max_v /= 0.5;
	printf("Maximum resolution is %f nm^-1\n", max_res);

	if ( max_intensity <= 0.0 ) {
		max_r = 4.0;
		goto out;
127
128
	}

Thomas White's avatar
Thomas White committed
129
130
131
132
133
	/* Choose whichever scaling factor gives the smallest value */
	scale = ((double)wh-50.0) / (2*max_u);
	if ( ((double)ht-50.0) / (2*max_v) < scale ) {
		scale = ((double)ht-50.0) / (2*max_v);
	}
134

Thomas White's avatar
Thomas White committed
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
	sep_u = as * scale * cos(theta);
	sep_v = bs * scale;
	max_r = ((sep_u < sep_v)?sep_u:sep_v);

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

		double u, v, intensity, val;
		int ct;

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

		intensity = lookup_intensity(ref, h, k, 0) / (float)ct;
		val = 10.0*intensity/max_intensity;

		u = (double)h*as*sin(theta);
		v = (double)h*as*cos(theta) + k*bs;

		cairo_arc(dctx, ((double)wh/2)+u*scale*2,
				((double)ht/2)+v*scale*2, max_r, 0, 2*M_PI);

		cairo_set_source_rgb(dctx, val, val, val);
		cairo_fill(dctx);

	}
161
162
	}

Thomas White's avatar
Thomas White committed
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
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);

	cairo_surface_finish(surface);
	cairo_destroy(dctx);
}


static void povray_render_animation(UnitCell *cell, double *ref,
                                    unsigned int *c, unsigned int nproc)
{
	FILE *fh;
	double asx, asy, asz;
	double bsx, bsy, bsz;
	double csx, csy, csz;
	pid_t pids[MAX_PROC];
	float max;
	int i;
	signed int h, k, l;

187
	fh = fopen("render.pov", "w");
Thomas White's avatar
Thomas White committed
188
	fprintf(fh, "/* POV-Ray scene written by CrystFEL */\n\n");
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
	fprintf(fh, "#include \"colors.inc\"\n");
	fprintf(fh, "#include \"textures.inc\"\n\n");
	fprintf(fh, "global_settings {\n");
	fprintf(fh, "	assumed_gamma 1.0\n");
	fprintf(fh, "	ambient_light 5.0\n");
	fprintf(fh, "}\n\n");

	/* First quarter */
	fprintf(fh, "#if ( (clock >= 0) & (clock <= 124) )\n");
	fprintf(fh, "camera { location <0.0, -3.0, 0.0>"
	            " sky z direction 1.1*y\n"
	            " right -x*(image_width/image_height)\n"
	            " look_at <0.0, 0.0, 0.0> }\n\n");
	fprintf(fh, "#end\n");

	/* Second quarter */
205
	fprintf(fh, "#if ( (clock >= 125) & (clock <= 249) )\n");
206
207
	fprintf(fh, "camera { location <0.0,"
	            " -(2.0+cos(radians((clock-125)*(180/125)))), 0.0>"
208
209
210
211
212
213
214
	            " sky z direction 1.1*y\n"
	            " right -x*(image_width/image_height)\n"
	            " look_at <0.0, 0.0, 0.0> }\n\n");
	fprintf(fh, "#end\n");

	/* Third quarter */
	fprintf(fh, "#if ( (clock >= 250) & (clock <= 374) )\n");
Thomas White's avatar
Thomas White committed
215
	fprintf(fh, "camera { location <0.0, -1.0, 0.0>"
216
217
218
219
220
221
222
	            " sky z direction 1.1*y\n"
	            " right -x*(image_width/image_height)\n"
	            " look_at <0.0, 0.0, 0.0> }\n\n");
	fprintf(fh, "#end\n");

	/* Fourth quarter */
	fprintf(fh, "#if ( (clock >= 375) & (clock <= 500) )\n");
223
224
	fprintf(fh, "camera { location <0.0,"
	            " -(2.0+cos(radians((clock-375)*(180/125)+180))), 0.0>"
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
	            " sky z direction 1.1*y\n"
	            " right -x*(image_width/image_height)\n"
	            " look_at <0.0, 0.0, 0.0> }\n\n");
	fprintf(fh, "#end\n");

	fprintf(fh, "light_source { <-3.0 -3.0 3.0> White }\n");
	fprintf(fh, "light_source { <+3.0 -3.0 3.0> White }\n");
	fprintf(fh, "light_source { <0.0, -3.0, 0.0> 2*White }\n");
	fprintf(fh, "plane {z,-2.0 pigment { rgb <0.0, 0.0, 0.1> } }\n");
	fprintf(fh, "plane {-z,-2.0 pigment { rgb <0.0, 0.0, 0.05> } }\n\n");

	cell_get_reciprocal(cell, &asx, &asy, &asz,
	                          &bsx, &bsy, &bsz,
	                          &csx, &csy, &csz);

240
241
242
243
	fprintf(fh, "#declare WCA = (720/19);\n");
	fprintf(fh, "#declare WCL = (360/8.5);\n");
	fprintf(fh, "#declare TA = (4.875);\n");
	fprintf(fh, "#declare TB = (1.125);\n");
244
245
246
247
248
249
250
251

	fprintf(fh, "#declare TRANS = \n");
	fprintf(fh, "transform {\n");

	/* First half */

	/* Acceleration */
	fprintf(fh, "#if ( clock <= 24 )\n"
252
	"rotate <0, 0, 0.5*WCA*(clock/25)*(clock/25)>\n"
253
254
255
256
	"#end\n"

	/* Cruise */
	"#if ( (clock >= 25) & (clock <= 224) )\n"
257
	"rotate <0, 0, (WCA/2)+WCA*((clock-25)/25)>\n"
258
259
260
261
262
263
	"#end\n"

	/* Overlap */

	/* Deceleration */
	"#if ( (clock >= 225) & (clock <= 274) )\n"
264
265
	"rotate <0, 0, 360-WCA + WCA*((clock-225)/25) "
	" - 0.5*(WCA/2)*((clock-225)/25)*((clock-225)/25) >\n"
266
267
268
269
	"#end\n"

	/* Acceleration */
	"#if ( (clock >= 225) & (clock <= 274) )\n"
270
	"rotate <0.5*(WCL/2)*((clock-225)/25)*((clock-225)/25), 0, 0>\n"
271
272
273
274
275
	"#end\n"

	/* Second half */

	/* Cruise */
276
277
	"#if ( (clock >= 275) & (clock <= 396) )\n"
	"rotate <WCL + WCL*((clock-275)/25), 0, 0>\n"
278
279
	"#end\n"

280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
	/* Deceleration to pause */
	"#if ( (clock >= 397) & (clock <= 421) )\n"
	"rotate <(1+TA)*WCL+ WCL*((clock-397)/25) "
	" - 0.5*WCL*((clock-397)/25)*((clock-397)/25), 0, 0 >\n"
	"#end\n"

	/* Acceleration after pause */
	"#if ( (clock >= 422) & (clock <= 446) )\n"
	"rotate <(1.5+TA)*WCL"
	" + 0.5*WCL*((clock-422)/25)*((clock-422)/25), 0, 0>\n"
	"#end\n"

	/* Final Cruise */
	"#if ( (clock >= 447) & (clock <= 474) )\n"
	"rotate <(2+TA)*WCL + WCL*((clock-447)/25), 0, 0>\n"
	"#end\n"

	/* Final Deceleration */
298
	"#if ( (clock >= 475) & (clock <= 499) )\n"
299
300
	"rotate <(2+TA+TB)*WCL + WCL*((clock-475)/25) "
	" - 0.5*WCL*((clock-475)/25)*((clock-475)/25), 0, 0 >\n"
301
302
303
304
305
306
307
308
309
310
311
312
313
	"#end\n");

	fprintf(fh, "}\n");

	max = 0.5e6;
	for ( h=-INDMAX; h<INDMAX; h++ ) {
	for ( k=-INDMAX; k<INDMAX; k++ ) {
	for ( l=-INDMAX; l<INDMAX; l++ ) {

		float radius, x, y, z;
		int s;
		float val, p, r, g, b, trans;

Thomas White's avatar
Thomas White committed
314
		if ( !lookup_count(c, h, k, l) ) continue;
315

Thomas White's avatar
Thomas White committed
316
		val = lookup_intensity(ref, h, k, l);
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333

		val = max-val;

		s = val / (max/6);
		p = fmod(val, max/6);
		p /= (max/6);

		r = 0;	g = 0;	b = 0;

		if ( (val < 0.0) ) {
			s = 0;
			p = 1.0;
		}
		if ( (val > max) ) {
			s = 6;
		}
		switch ( s ) {
Thomas White's avatar
Thomas White committed
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
		case 0 :   /* Black to blue */
			r = 0.0;  g = 0.0;  b = p;
			break;
		case 1 :   /* Blue to green */
			r = 0.0;  g = p;  b = 1.0-p;
			break;
		case 2 :   /* Green to red */
			r =p;  g = 1.0-p;  b = 0.0;
			break;
		case 3 :   /* Red to Orange */
			r = 1.0;  g = 0.5*p;  b = 0.0;
			break;
		case 4 :   /* Orange to Yellow */
			r = 1.0;  g = 0.5 + 0.5*p;  b = 0.0;
			break;
		case 5 :   /* Yellow to White */
			r = 1.0;  g = 1.0;  b = 1.0*p;
			break;
		case 6 :   /* Pixel has hit the maximum value */
			r = 1.0;  g = 1.0;  b = 1.0;
			break;
355
356
357
358
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
		}

		val = max-val;

		if ( val <= 0.0 ) continue;
		radius = 0.1 * sqrt(sqrt(val))/1e2;
		radius -= 0.005;
		if ( radius > 0.03 ) radius = 0.03;
		if ( radius <= 0.0 ) continue;
		trans = (0.03-radius)/0.03;
		radius += 0.002;

		x = asx*h + bsx*k + csx*l;
		y = asy*h + bsy*k + csy*l;
		z = asz*h + bsz*k + csz*l;

		fprintf(fh, "sphere { <%.5f, %.5f, %.5f>, %.5f "
		            "texture{pigment{color rgb <%f, %f, %f>"
		            " transmit %f} "
		            "finish { reflection 0.1 } } \n"
			    "transform { TRANS }\n"
			    "}\n",
			x/1e9, y/1e9, z/1e9, radius, r, g, b, trans);

	}
	}
	}

	fprintf(fh, "\n");
	fclose(fh);

386
387
388
389
390
391
392
393
394
395
	for ( i=0; i<nproc; i++ ) {

		pids[i] = fork();
		if ( !( (pids[i] != 0) && (pids[i] != -1) ) ) {
			if ( pids[i] == -1 ) {
				ERROR("fork() failed.\n");
			} else {

				char minf[256];
				char maxf[256];
Thomas White's avatar
Thomas White committed
396
				float nf, xf, nsec;
397

Thomas White's avatar
Thomas White committed
398
399
400
				nsec = 500.0 / (float)nproc;
				nf = nsec * (float)i;
				xf = (nsec * (float)i + nsec) - 1.0;
401

Thomas White's avatar
Thomas White committed
402
403
				snprintf(minf, 255, "+SF%i", (int)nf);
				snprintf(maxf, 255, "+EF%i", (int)xf);
404
405
406
407
408

				/* Forked successfully, child process */
				execlp("povray", "", "+W1024", "+H768",
				       "+Irender.pov", "+Orender.png",
				       "+KFI0", "+KFF499", "+KI0", "+KF499",
409
				       minf, maxf, "-D", NULL);
410
411
412
413
414
415

			}
		} /* else start the next one */
	}

	for ( i=0; i<nproc; i++ ) {
416
		int r;
417
		waitpid(pids[i], &r, 0);
418
	}
Thomas White's avatar
Thomas White committed
419
420
421
422
423
424
425
426
427
428
429
430
431
}


int main(int argc, char *argv[])
{
	int c;
	UnitCell *cell;
	char *infile;
	double *ref;
	unsigned int *cts;
	int config_povray = 0;
	int config_zoneaxis = 0;
	unsigned int nproc = 1;
432
	char *pdb = NULL;
Thomas White's avatar
Thomas White committed
433
434
435
436
437
438

	/* Long options */
	const struct option longopts[] = {
		{"help",               0, NULL,               'h'},
		{"povray",             0, &config_povray,      1},
		{"zone-axis",          0, &config_zoneaxis,    1},
439
		{"pdb",                1, NULL,               'p'},
Thomas White's avatar
Thomas White committed
440
441
442
443
		{0, 0, NULL, 0}
	};

	/* Short options */
444
	while ((c = getopt_long(argc, argv, "hj:p:", longopts, NULL)) != -1) {
Thomas White's avatar
Thomas White committed
445
446

		switch (c) {
Thomas White's avatar
Thomas White committed
447
		case 'h' :
Thomas White's avatar
Thomas White committed
448
449
450
			show_help(argv[0]);
			return 0;

Thomas White's avatar
Thomas White committed
451
		case 'j' :
Thomas White's avatar
Thomas White committed
452
453
454
			nproc = atoi(optarg);
			break;

455
456
457
458
		case 'p' :
			pdb = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
459
		case 0 :
Thomas White's avatar
Thomas White committed
460
461
			break;

Thomas White's avatar
Thomas White committed
462
		default :
Thomas White's avatar
Thomas White committed
463
464
465
466
467
			return 1;
		}

	}

468
469
470
471
	if ( pdb == NULL ) {
		pdb = strdup("molecule.pdb");
	}

Thomas White's avatar
Thomas White committed
472
473
474
475
476
477
478
	if ( (nproc > MAX_PROC) || (nproc < 1) ) {
		ERROR("Number of processes is invalid.\n");
		return 1;
	}

	infile = argv[optind];

479
	cell = load_cell_from_pdb(pdb);
Thomas White's avatar
Thomas White committed
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
	cts = new_list_count();
	ref = read_reflections(infile, cts, NULL);
	if ( ref == NULL ) {
		ERROR("Couldn't open file '%s'\n", infile);
		return 1;
	}

	if ( config_povray ) {
		povray_render_animation(cell, ref, cts, nproc);
	} else if ( config_zoneaxis ) {
		render_za(cell, ref, cts);
	} else {
		ERROR("Try again with either --povray or --zone-axis.\n");
	}

495
	free(pdb);
496
497
498

	return 0;
}