render_hkl.c 5.22 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
	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 */
81
	if ( cell_get_reciprocal(cell, &asx, &asy, &asz,
Thomas White's avatar
Thomas White committed
82
	                          &bsx, &bsy, &bsz,
83
84
85
86
	                          &csx, &csy, &csz) ) {
		ERROR("Couldn't get reciprocal parameters\n");
		return;
	}
Thomas White's avatar
Thomas White committed
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
	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);
113
114
115
		}

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

Thomas White's avatar
Thomas White committed
118
119
120
121
122
123
124
125
	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;
126
127
	}

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

Thomas White's avatar
Thomas White committed
134
135
136
137
138
139
140
141
142
143
144
145
146
147
	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;
Thomas White's avatar
Thomas White committed
148
		val = 3.0*intensity/max_intensity;
Thomas White's avatar
Thomas White committed
149
150
151
152
153
154
155
156
157
158
159

		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);

	}
160
161
	}

Thomas White's avatar
Thomas White committed
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
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;
184
	char *pdb = NULL;
185
	int r = 0;
Thomas White's avatar
Thomas White committed
186
187
188
189
190
191

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

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

		switch (c) {
Thomas White's avatar
Thomas White committed
200
		case 'h' :
Thomas White's avatar
Thomas White committed
201
202
203
			show_help(argv[0]);
			return 0;

Thomas White's avatar
Thomas White committed
204
		case 'j' :
Thomas White's avatar
Thomas White committed
205
206
207
			nproc = atoi(optarg);
			break;

208
209
210
211
		case 'p' :
			pdb = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
212
		case 0 :
Thomas White's avatar
Thomas White committed
213
214
			break;

Thomas White's avatar
Thomas White committed
215
		default :
Thomas White's avatar
Thomas White committed
216
217
218
219
220
			return 1;
		}

	}

221
222
223
224
	if ( pdb == NULL ) {
		pdb = strdup("molecule.pdb");
	}

Thomas White's avatar
Thomas White committed
225
226
	infile = argv[optind];

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

247
	free(pdb);
248

249
	return r;
250
}