render_hkl.c 5.16 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>
Thomas White's avatar
Thomas White committed
23
24
#include <cairo.h>
#include <cairo-pdf.h>
25
26
27

#include "utils.h"
#include "reflections.h"
28
#include "povray.h"
29
30


31
32
33
34
static void show_help(const char *s)
{
	printf("Syntax: %s [options] <file.hkl>\n\n", s);
	printf(
35
"Render intensity lists in various ways.\n"
36
"\n"
37
38
39
40
41
42
"  -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"
);
43
44
45
}


Thomas White's avatar
Thomas White committed
46
static void render_za(UnitCell *cell, double *ref, unsigned int *c)
47
{
Thomas White's avatar
Thomas White committed
48
49
50
51
52
	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;
53
54
55
	double asx, asy, asz;
	double bsx, bsy, bsz;
	double csx, csy, csz;
Thomas White's avatar
Thomas White committed
56
57
	signed int h, k;
	float wh, ht;
58

Thomas White's avatar
Thomas White committed
59
60
	wh = 1024;
	ht = 1024;
61

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

Thomas White's avatar
Thomas White committed
64
65
66
67
68
	if ( cairo_surface_status(surface) != CAIRO_STATUS_SUCCESS ) {
		fprintf(stderr, "Couldn't create Cairo surface\n");
		cairo_surface_destroy(surface);
		return;
	}
69

Thomas White's avatar
Thomas White committed
70
	dctx = cairo_create(surface);
71

Thomas White's avatar
Thomas White committed
72
73
74
75
	/* 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);
76

Thomas White's avatar
Thomas White committed
77
78
79
80
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
	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);
110
111
112
		}

	}
Thomas White's avatar
Thomas White committed
113
	}
114

Thomas White's avatar
Thomas White committed
115
116
117
118
119
120
121
122
	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;
123
124
	}

Thomas White's avatar
Thomas White committed
125
126
127
128
129
	/* 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);
	}
130

Thomas White's avatar
Thomas White committed
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
	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);

	}
157
158
	}

Thomas White's avatar
Thomas White committed
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
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);
}


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;
181
	char *pdb = NULL;
182
	int r = 0;
Thomas White's avatar
Thomas White committed
183
184
185
186
187
188

	/* Long options */
	const struct option longopts[] = {
		{"help",               0, NULL,               'h'},
		{"povray",             0, &config_povray,      1},
		{"zone-axis",          0, &config_zoneaxis,    1},
189
		{"pdb",                1, NULL,               'p'},
Thomas White's avatar
Thomas White committed
190
191
192
193
		{0, 0, NULL, 0}
	};

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

		switch (c) {
Thomas White's avatar
Thomas White committed
197
		case 'h' :
Thomas White's avatar
Thomas White committed
198
199
200
			show_help(argv[0]);
			return 0;

Thomas White's avatar
Thomas White committed
201
		case 'j' :
Thomas White's avatar
Thomas White committed
202
203
204
			nproc = atoi(optarg);
			break;

205
206
207
208
		case 'p' :
			pdb = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
209
		case 0 :
Thomas White's avatar
Thomas White committed
210
211
			break;

Thomas White's avatar
Thomas White committed
212
		default :
Thomas White's avatar
Thomas White committed
213
214
215
216
217
			return 1;
		}

	}

218
219
220
221
	if ( pdb == NULL ) {
		pdb = strdup("molecule.pdb");
	}

Thomas White's avatar
Thomas White committed
222
223
	infile = argv[optind];

224
	cell = load_cell_from_pdb(pdb);
225
226
227
228
	if ( cell == NULL ) {
		ERROR("Couldn't load unit cell from %s\n", pdb);
		return 1;
	}
Thomas White's avatar
Thomas White committed
229
230
231
232
233
234
235
236
	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 ) {
237
		r = povray_render_animation(cell, ref, cts, nproc);
Thomas White's avatar
Thomas White committed
238
239
240
241
242
243
	} else if ( config_zoneaxis ) {
		render_za(cell, ref, cts);
	} else {
		ERROR("Try again with either --povray or --zone-axis.\n");
	}

244
	free(pdb);
245

246
	return r;
247
}