render_hkl: Multiple improvements

6 files changed, 160 insertions, 43 deletions

diff --git a/src/cubeit.c b/src/cubeit.c
index 33442ac3..51f4a57c 100644
--- a/src/cubeit.c
+++ b/src/cubeit.c
@@ -361,8 +361,8 @@ static void write_slice(const char *filename, double *vals, int z,
 	for ( y=0; y<h; y++ ) {
 	for ( x=0; x<w; x++ ) {
 
-		float r, g, b;
-		float val;
+		double r, g, b;
+		double val;
 
 		val = vals[xs*ys*z + xs*y + x];
 
diff --git a/src/povray.c b/src/povray.c
index 8d0ec048..660f1c59 100644
--- a/src/povray.c
+++ b/src/povray.c
@@ -30,7 +30,8 @@
 
 
 int povray_render_animation(UnitCell *cell, RefList *list, unsigned int nproc,
-                            const char *sym, int wght, double boost)
+                            const char *sym, int wght, double boost,
+                            double scale_top)
 {
 	FILE *fh;
 	double asx, asy, asz;
@@ -201,6 +202,12 @@ int povray_render_animation(UnitCell *cell, RefList *list, unsigned int nproc,
 	}
 	max /= boost;
 
+	/* Use manual scale top if specified */
+	if ( scale_top > 0.0 ) {
+		max = scale_top;
+	}
+
+
 	for ( refl = first_refl(list, &iter);
 	      refl != NULL;
 	      refl = next_refl(refl, iter) ) {
diff --git a/src/povray.h b/src/povray.h
index 0b0b3446..9dbf109d 100644
--- a/src/povray.h
+++ b/src/povray.h
@@ -21,6 +21,6 @@
 
 extern int povray_render_animation(UnitCell *cell, RefList *list,
                                    unsigned int nproc, const char *sym,
-                                   int wght, double boost);
+                                   int wght, double boost, double scale_top);
 
 #endif	/* POVRAY_H */
diff --git a/src/render.c b/src/render.c
index 63357083..46f93566 100644
--- a/src/render.c
+++ b/src/render.c
@@ -38,11 +38,12 @@
 #include "utils.h"
 
 
-static void render_rgb(float val, float max, float *rp, float *gp, float *bp)
+static void render_rgb(double val, double max,
+                       double *rp, double *gp, double *bp)
 {
 	int s;
-	float p;
-	float r, g, b;
+	double p;
+	double r, g, b;
 
 	s = val / (max/6);
 	p = fmod(val, max/6.0);
@@ -94,9 +95,23 @@ static void render_rgb(float val, float max, float *rp, float *gp, float *bp)
 }
 
 
-static void render_mono(float val, float max, float *rp, float *gp, float *bp)
+static void render_ratio(double val, double max,
+                         double *rp, double *gp, double *bp)
 {
-	float p;
+	if ( val <= 1.0 ) {
+		render_rgb(val, 2.0, rp, gp, bp);
+	} else {
+		/* Your homework is to simplify this expression */
+		val = ((val-1.0)/(max-1.0)) * (max/2.0) + max/2.0;
+		render_rgb(val, max, rp, gp, bp);
+	}
+}
+
+
+static void render_mono(double val, double max,
+                        double *rp, double *gp, double *bp)
+{
+	double p;
 	p = val / max;
 	if ( val < 0.0 ) p = 0.0;
 	if ( val > max ) p = 1.0;
@@ -106,10 +121,10 @@ static void render_mono(float val, float max, float *rp, float *gp, float *bp)
 }
 
 
-static void render_invmono(float val, float max,
-                           float *rp, float *gp, float *bp)
+static void render_invmono(double val, double max,
+                           double *rp, double *gp, double *bp)
 {
-	float p;
+	double p;
 	p = val / max;
 	p = 1.0 - p;
 	if ( val < 0.0 ) p = 1.0;
@@ -120,8 +135,8 @@ static void render_invmono(float val, float max,
 }
 
 
-void render_scale(float val, float max, int scale,
-                  float *rp, float *gp, float *bp)
+void render_scale(double val, double max, int scale,
+                  double *rp, double *gp, double *bp)
 {
 	switch ( scale ) {
 	case SCALE_COLOUR :
@@ -133,6 +148,9 @@ void render_scale(float val, float max, int scale,
 	case SCALE_INVMONO :
 		render_invmono(val, max, rp, gp, bp);
 		break;
+	case SCALE_RATIO :
+		render_ratio(val, max, rp, gp, bp);
+		break;
 	}
 }
 
@@ -232,8 +250,8 @@ static GdkPixbuf *render_panel(struct image *image,
 	for ( y=0; y<h; y++ ) {
 	for ( x=0; x<w; x++ ) {
 
-		float val;
-		float r, g, b;
+		double val;
+		double r, g, b;
 
 		val = hdr[x+w*y];
 		render_scale(val, max, scale, &r, &g, &b);
@@ -302,7 +320,7 @@ GdkPixbuf *render_get_colour_scale(size_t w, size_t h, int scale)
 
 	for ( y=0; y<h; y++ ) {
 
-		float r, g, b;
+		double r, g, b;
 		int val;
 
 		val = y;
@@ -369,7 +387,7 @@ int render_tiff_int16(struct image *image, const char *filename, double boost)
 	TIFF *th;
 	int16_t *line;
 	int x, y;
-	float max;
+	double max;
 
 	th = TIFFOpen(filename, "w");
 	if ( th == NULL ) return 1;
@@ -389,7 +407,7 @@ int render_tiff_int16(struct image *image, const char *filename, double boost)
 	max = 0.0;
 	for ( y=0; y<image->height; y++ ) {
 	for ( x=0;x<image->width; x++ ) {
-		float val;
+		double val;
 		val = image->data[x+image->height*y];
 		if ( val > max ) max = val;
 	}
@@ -399,10 +417,10 @@ int render_tiff_int16(struct image *image, const char *filename, double boost)
 	for ( y=0; y<image->height; y++ ) {
 		for ( x=0;x<image->width; x++ ) {
 
-			float val;
+			double val;
 
 			val = image->data[x+(image->height-1-y)*image->width];
-			val *= ((float)boost/max);
+			val *= ((double)boost/max);
 
 			/* Clamp to 16-bit range,
 			 * and work round inability of most readers to deal
diff --git a/src/render.h b/src/render.h
index e8017de8..fbfbbdb2 100644
--- a/src/render.h
+++ b/src/render.h
@@ -25,12 +25,13 @@
 enum {
 	SCALE_COLOUR,
 	SCALE_MONO,
-	SCALE_INVMONO
+	SCALE_INVMONO,
+	SCALE_RATIO
 };
 
 /* Colour scale lookup */
-extern void render_scale(float val, float max, int scale,
-                         float *rp, float *gp, float *bp);
+extern void render_scale(double val, double max, int scale,
+                         double *rp, double *gp, double *bp);
 
 
 #ifdef HAVE_GTK
@@ -38,8 +39,8 @@ extern void render_scale(float val, float max, int scale,
 #include <gdk-pixbuf/gdk-pixbuf.h>
 
 extern GdkPixbuf **render_panels(struct image *image,
-                                  int binning, int scale, double boost,
-                                  int *n_pixbufs);
+                                 int binning, int scale, double boost,
+                                 int *n_pixbufs);
 
 extern GdkPixbuf *render_get_colour_scale(size_t w, size_t h, int scale);
 
diff --git a/src/render_hkl.c b/src/render_hkl.c
index 44ec23fe..8ac73a72 100644
--- a/src/render_hkl.c
+++ b/src/render_hkl.c
@@ -143,7 +143,7 @@ static void draw_circles(signed int xh, signed int xk, signed int xl,
 
 		if ( dctx != NULL ) {
 
-			float r, g, b;
+			double r, g, b;
 
 			cairo_arc(dctx, ((double)cx)+u*scale,
 			                ((double)cy)+v*scale,
@@ -238,7 +238,7 @@ static void render_za(UnitCell *cell, RefList *list,
                       double boost, const char *sym, int wght, int colscale,
                       signed int xh, signed int xk, signed int xl,
                       signed int yh, signed int yk, signed int yl,
-                      const char *outfile)
+                      const char *outfile, double scale_top)
 {
 	cairo_surface_t *surface;
 	cairo_t *dctx;
@@ -259,6 +259,7 @@ static void render_za(UnitCell *cell, RefList *list,
 	cairo_text_extents_t size;
 	double cx, cy;
 	const double border = 200.0;
+	int png;
 
 	/* Vector product to determine the zone axis. */
 	zh = xk*yl - xl*yk;
@@ -302,6 +303,11 @@ static void render_za(UnitCell *cell, RefList *list,
 		return;
 	}
 
+	/* Use manual scale top if specified */
+	if ( scale_top > 0.0 ) {
+		max_val = scale_top;
+	}
+
 	/* Choose whichever scaling factor gives the smallest value */
 	scale_u = ((double)wh-border) / (2.0*max_u);
 	scale_v = ((double)ht-border) / (2.0*max_v);
@@ -318,7 +324,15 @@ static void render_za(UnitCell *cell, RefList *list,
 	}
 
 	if ( outfile == NULL ) outfile = "za.pdf";
-	surface = cairo_pdf_surface_create(outfile, wh, ht);
+
+	if ( strcmp(outfile+strlen(outfile)-4, ".png") == 0 ) {
+		png = 1;
+		surface = cairo_image_surface_create(CAIRO_FORMAT_ARGB32,
+		                                     wh, ht);
+	} else {
+		png = 0;
+		surface = cairo_pdf_surface_create(outfile, wh, ht);
+	}
 
 	if ( cairo_surface_status(surface) != CAIRO_STATUS_SUCCESS ) {
 		ERROR("Couldn't create Cairo surface\n");
@@ -390,20 +404,34 @@ static void render_za(UnitCell *cell, RefList *list,
 	render_overlined_indices(dctx, yh, yk, yl);
 	cairo_fill(dctx);
 
+	if ( png ) {
+		int r = cairo_surface_write_to_png(surface, outfile);
+		if ( r != CAIRO_STATUS_SUCCESS ) {
+			ERROR("Failed to write PNG to '%s'\n", outfile);
+		}
+	}
+
 	cairo_surface_finish(surface);
 	cairo_destroy(dctx);
 }
 
 
-static int render_key(int colscale)
+static int render_key(int colscale, double scale_top)
 {
 	cairo_surface_t *surface;
 	cairo_t *dctx;
-	float wh, ht;
-	float y;
+	double top, wh, ht, y;
+	double slice;
 
-	wh = 128;
-	ht = 1024;
+	wh = 128.0;
+	ht = 1024.0;
+	slice = 1.0;
+
+	if ( scale_top > 0.0 ) {
+		top = scale_top;
+	} else {
+		top = 1.0;
+	}
 
 	surface = cairo_pdf_surface_create("key.pdf", wh, ht);
 
@@ -415,19 +443,61 @@ static int render_key(int colscale)
 
 	dctx = cairo_create(surface);
 
-	for ( y=0; y<ht; y++ ) {
+	for ( y=0.0; y<ht; y+=slice ) {
 
-		float r, g, b;
+		double r, g, b;
+		double val;
+		double v = y;
 
-		cairo_rectangle(dctx, 0.0, y, wh, y+1.0);
+		cairo_rectangle(dctx, 0.0, ht-y, wh/2.0, slice);
+
+		if ( colscale == SCALE_RATIO ) {
+			if ( v < ht/2.0 ) {
+				val = v/(ht/2.0);
+			} else {
+				val = (((v-ht/2.0)/(ht/2.0))*(top-1.0))+1.0;
+			}
+		} else {
+			val = v/ht;
+		}
 
-		render_scale(ht-y, ht, colscale, &r, &g, &b);
+		render_scale(val, top, colscale, &r, &g, &b);
 		cairo_set_source_rgb(dctx, r, g, b);
 
+		cairo_stroke_preserve(dctx);
 		cairo_fill(dctx);
 
 	}
 
+	if ( colscale == SCALE_RATIO ) {
+
+		cairo_text_extents_t size;
+		char tmp[32];
+
+		cairo_rectangle(dctx, 0.0, ht/2.0-2.0, wh/2.0, 4.0);
+		cairo_set_source_rgb(dctx, 0.0, 0.0, 0.0);
+		cairo_stroke_preserve(dctx);
+		cairo_fill(dctx);
+
+		cairo_set_font_size(dctx, 20.0);
+		cairo_text_extents(dctx, "1.0", &size);
+		cairo_move_to(dctx, wh/2.0+5.0, ht/2.0+size.height/2.0);
+		cairo_show_text(dctx, "1.0");
+
+		cairo_set_font_size(dctx, 20.0);
+		cairo_text_extents(dctx, "0.0", &size);
+		cairo_move_to(dctx, wh/2.0+5.0, ht-5.0);
+		cairo_show_text(dctx, "0.0");
+
+		cairo_set_font_size(dctx, 20.0);
+		snprintf(tmp, 31, "%.1f", top);
+		cairo_text_extents(dctx, tmp, &size);
+		cairo_move_to(dctx, wh/2.0+5.0, size.height+5.0);
+		cairo_show_text(dctx, tmp);
+
+	}
+
+
 	cairo_surface_finish(surface);
 	cairo_destroy(dctx);
 
@@ -487,6 +557,8 @@ int main(int argc, char *argv[])
 	char *down = NULL;
 	char *right = NULL;
 	char *outfile = NULL;
+	double scale_top = -1.0;
+	char *endptr;
 
 	/* Long options */
 	const struct option longopts[] = {
@@ -503,6 +575,7 @@ int main(int argc, char *argv[])
 		{"right",              1, NULL,               'r'},
 		{"counts",             0, &config_sqrt,        1},
 		{"colour-key",         0, &config_colkey,      1},
+		{"scale-top",          1, NULL,                2},
 		{0, 0, NULL, 0}
 	};
 
@@ -551,6 +624,21 @@ int main(int argc, char *argv[])
 			outfile = strdup(optarg);
 			break;
 
+		case 2 :
+			errno = 0;
+			scale_top = strtod(optarg, &endptr);
+			if ( !( (optarg[0] != '\0') && (endptr[0] == '\0') )
+			   || (errno != 0) )
+			{
+				ERROR("Invalid scale top('%s')\n", optarg);
+				return 1;
+			}
+			if ( scale_top < 0.0 ) {
+				ERROR("Scale top must be positive.\n");
+				return 1;
+			}
+			break;
+
 		case 0 :
 			break;
 
@@ -560,8 +648,9 @@ int main(int argc, char *argv[])
 
 	}
 
-	if ( pdb == NULL ) {
-		pdb = strdup("molecule.pdb");
+	if ( (pdb == NULL) && !config_colkey ) {
+		ERROR("You must specify the PDB containing the unit cell.\n");
+		return 1;
 	}
 
 	if ( sym == NULL ) {
@@ -604,6 +693,8 @@ int main(int argc, char *argv[])
 		colscale = SCALE_COLOUR;
 	} else if ( strcmp(cscale, "color") == 0 ) {
 		colscale = SCALE_COLOUR;
+	} else if ( strcmp(cscale, "ratio") == 0 ) {
+		colscale = SCALE_RATIO;
 	} else {
 		ERROR("Unrecognised colour scale '%s'\n", cscale);
 		return 1;
@@ -611,7 +702,7 @@ int main(int argc, char *argv[])
 	free(cscale);
 
 	if ( config_colkey ) {
-		return render_key(colscale);
+		return render_key(colscale, scale_top);
 	}
 
 	if ( config_zoneaxis ) {
@@ -659,10 +750,10 @@ int main(int argc, char *argv[])
 
 	if ( config_povray ) {
 		r = povray_render_animation(cell, list,
-		                            nproc, sym, wght, boost);
+		                            nproc, sym, wght, boost, scale_top);
 	} else if ( config_zoneaxis ) {
 		render_za(cell, list, boost, sym, wght, colscale,
-		          rh, rk, rl, dh, dk, dl, outfile);
+		          rh, rk, rl, dh, dk, dl, outfile, scale_top);
 	} else {
 		ERROR("Try again with either --povray or --zone-axis.\n");
 	}