render_hkl: Also do 2D

1 files changed, 195 insertions, 47 deletions

diff --git a/src/render_hkl.c b/src/render_hkl.c
index c4555242..e18ae90d 100644
--- a/src/render_hkl.c
+++ b/src/render_hkl.c
@@ -22,6 +22,8 @@
 #include <getopt.h>
 #include <sys/types.h>
 #include <sys/wait.h>
+#include <cairo.h>
+#include <cairo-pdf.h>
 
 #include "utils.h"
 #include "reflections.h"
@@ -37,76 +39,152 @@ static void show_help(const char *s)
 "Render intensity lists using POV-ray.\n"
 "\n"
 "  -h, --help      Display this help message.\n"
+"      --povray    Render a 3D animation using POV-ray.\n"
 "  -j <n>          Run <n> instances of POV-ray in parallel.\n"
+"      --zone-axis Render a 2D zone axis pattern.\n"
 "\n");
 }
 
 
-int main(int argc, char *argv[])
+static void render_za(UnitCell *cell, double *ref, unsigned int *c)
 {
-	int c;
-	UnitCell *cell;
-	char *infile;
-	double *ref1;
-	signed int h, k, l;
-	unsigned int *c1;
-	FILE *fh;
+	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;
 	double asx, asy, asz;
 	double bsx, bsy, bsz;
 	double csx, csy, csz;
-	pid_t pids[MAX_PROC];
-	float max;
-	int nproc = 1;
-	int i;
+	signed int h, k;
+	float wh, ht;
 
-	/* Long options */
-	const struct option longopts[] = {
-		{"help",               0, NULL,               'h'},
-		{0, 0, NULL, 0}
-	};
+	wh = 1024;
+	ht = 1024;
 
-	/* Short options */
-	while ((c = getopt_long(argc, argv, "hj:", longopts, NULL)) != -1) {
+	surface = cairo_pdf_surface_create("za.pdf", wh, ht);
 
-		switch (c) {
-		case 'h' : {
-			show_help(argv[0]);
-			return 0;
-		}
+	if ( cairo_surface_status(surface) != CAIRO_STATUS_SUCCESS ) {
+		fprintf(stderr, "Couldn't create Cairo surface\n");
+		cairo_surface_destroy(surface);
+		return;
+	}
 
-		case 'j' : {
-			nproc = atoi(optarg);
-			break;
-		}
+	dctx = cairo_create(surface);
 
-		case 0 : {
-			break;
-		}
+	/* 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);
 
-		default : {
-			return 1;
-		}
+	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);
 		}
 
 	}
+	}
 
-	if ( (nproc > MAX_PROC) || (nproc < 1) ) {
-		ERROR("Number of processes is invalid.\n");
-		return 1;
+	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;
 	}
 
-	infile = argv[optind];
+	/* 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);
+	}
 
-	cell = load_cell_from_pdb("molecule.pdb");
-	c1 = new_list_count();
-	ref1 = read_reflections(infile, c1, NULL);
-	if ( ref1 == NULL ) {
-		ERROR("Couldn't open file '%s'\n", infile);
-		return 1;
+	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);
+
+	}
 	}
 
+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;
+
 	fh = fopen("render.pov", "w");
-	fprintf(fh, "/* POV-Ray scene written by Synth3D */\n\n");
+	fprintf(fh, "/* POV-Ray scene written by CrystFEL */\n\n");
 	fprintf(fh, "#include \"colors.inc\"\n");
 	fprintf(fh, "#include \"textures.inc\"\n\n");
 	fprintf(fh, "global_settings {\n");
@@ -232,9 +310,9 @@ int main(int argc, char *argv[])
 		int s;
 		float val, p, r, g, b, trans;
 
-		if ( !lookup_count(c1, h, k, l) ) continue;
+		if ( !lookup_count(c, h, k, l) ) continue;
 
-		val = lookup_intensity(ref1, h, k, l);
+		val = lookup_intensity(ref, h, k, l);
 
 		val = max-val;
 
@@ -344,6 +422,76 @@ int main(int argc, char *argv[])
 		int r;
 		waitpid(pids[i], &r, 0);
 	}
+}
+
+
+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;
+
+	/* Long options */
+	const struct option longopts[] = {
+		{"help",               0, NULL,               'h'},
+		{"povray",             0, &config_povray,      1},
+		{"zone-axis",          0, &config_zoneaxis,    1},
+		{0, 0, NULL, 0}
+	};
+
+	/* Short options */
+	while ((c = getopt_long(argc, argv, "hj:", longopts, NULL)) != -1) {
+
+		switch (c) {
+		case 'h' : {
+			show_help(argv[0]);
+			return 0;
+		}
+
+		case 'j' : {
+			nproc = atoi(optarg);
+			break;
+		}
+
+		case 0 : {
+			break;
+		}
+
+		default : {
+			return 1;
+		}
+		}
+
+	}
+
+	if ( (nproc > MAX_PROC) || (nproc < 1) ) {
+		ERROR("Number of processes is invalid.\n");
+		return 1;
+	}
+
+	infile = argv[optind];
+
+	cell = load_cell_from_pdb("molecule.pdb");
+	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");
+	}
+
 
 	return 0;
 }