New release 1.7.0

6 files changed, 969 insertions, 222 deletions

diff --git a/render/Jamfile b/render/Jamfile
index fe975ea..ca199ab 100644
--- a/render/Jamfile
+++ b/render/Jamfile
@@ -18,12 +18,13 @@ InstallBin  $(DESTDIR)$(PREFIX)/bin : $(Executables) ;
 #InstallFile $(DESTDIR)$(PREFIX)/h : $(Headers) ;
 #InstallLib  $(DESTDIR)$(PREFIX)/lib : $(Libraries) ;
 
-HDRS = ../h ../numlib $(TIFFINC) ;
+HDRS = ../h ../numlib $(TIFFINC) $(PNGINC) ;
 
 # 2D Rendering library
 Library librender : render.c thscreen.c ;
 
-Main timage : timage.c : : : : : librender ../numlib/libnum $(TIFFLIB) $(JPEGLIB) ;
+Main timage : timage.c : : : : : librender ../numlib/libnum
+	          $(TIFFLIB) $(JPEGLIB) $(PNGLIB) $(ZLIB) ;
 
 
 if $(BUILD_JUNK) {
diff --git a/render/render.c b/render/render.c
index cbeb78d..de34ebd 100644
--- a/render/render.c
+++ b/render/render.c
@@ -14,18 +14,31 @@
  * see the License.txt file for licencing details.
  */
 
+/*
+ * TTBD: Should make this much more self contained in how
+ * it deals with errors - return an error code & string,
+ * and clean up resourcse.
+ */
+
 #undef DEBUG
+#undef CCTEST_PATTERN		/* Create ccast upsampling test pattern if */
+							/* "ARGYLL_CCAST_TEST_PATTERN" env variable is set */
+#undef TEST_SCREENING		/* For testing by making screen visible */
+
+#define OVERLAP 0.1			/* Stocastic screening overlap between levels */
 
 #define verbo stdout
 
 #include <stdio.h>
 #include <math.h>
 #include <stdlib.h>
+#include <fcntl.h>
 #include "copyright.h"
 #include "aconfig.h"
 #include "sort.h"
 #include "numlib.h"
 #include "tiffio.h"
+#include "png.h"
 #include "render.h"
 #include "thscreen.h"
 
@@ -85,6 +98,131 @@ static void cvt_Lab_to_CIELAB16(double *out, double *in) {
 }
 
 /* ------------------------------------------------------------- */
+/* PNG memory write support */
+typedef struct {
+	unsigned char *buf;
+	png_size_t len;			/* Current length of the buffer */
+	png_size_t off;			/* Current offset of the buffer */
+} png_mem_info;
+
+static void mem_write_data(png_structp png_ptr, png_bytep data, png_size_t length) {
+	png_mem_info *s = (png_mem_info *)png_get_io_ptr(png_ptr);
+
+	if ((s->off + length) > s->len) {			/* Need more space */
+		png_size_t more = (s->off + length) - s->len;
+
+		if (more < (1024 * 80))
+			more = 1024 * 50 - 32;		/* Increase 50K at a time */
+		s->len += more;
+
+		if ((s->buf = realloc(s->buf, s->len)) == NULL) {
+			png_error(png_ptr, "malloc failed in mem_write_data");
+		}
+	}
+	memcpy(s->buf + s->off, data, length);
+	s->off += length;
+}
+
+static void mem_flush_data(png_structp png_ptr) {
+	return;
+}
+
+/* ------------------------------------------------------------- */
+#ifdef CCTEST_PATTERN
+
+#define SG 9		// Spacing is 9 pixels
+//#define SG 64		// 10 lines per half screen
+
+static void test_value(render2d *s, tdata_t *outbuf, int xx, int yy) {
+	int x = xx, y = yy;
+	int i, j, v;
+	int oval[3];	/* 8 bit value */
+
+	if (x < s->pw/2) {
+		if (y < s->ph/2) {
+			/* Generate vertical white stripes */
+			/* Stripes are 5 pixels apart in groups of 2 of the same level, */
+			/* declining by 1 each group */
+			if ((x % SG) == 0) {
+				i = x / (2 * SG);
+				v = 255 - i;
+				if (v < 0)
+					v = 0;
+				oval[0] = oval[1] = oval[2] = v;
+			} else {
+				oval[0] = oval[1] = oval[2] = 0;
+			}
+		} else {
+			y -= s->ph/2;
+			/* Generate white dots */
+			if ((x % SG) == 0 && (y % SG) == 0) {
+				i = x / (2 * SG);
+				j = y / (2 * SG);
+
+				i += j * s->pw/2/(2 * SG);
+				oval[0] = oval[1] = oval[2] = 0;
+
+				for (j = 5; j >= 0; j--) {
+					oval[0] = (oval[0] << 1) | ((i >> (j * 3 + 0)) & 1);
+					oval[1] = (oval[1] << 1) | ((i >> (j * 3 + 1)) & 1);
+					oval[2] = (oval[2] << 1) | ((i >> (j * 3 + 2)) & 1);
+				}
+				oval[0] = 255 - oval[0];
+				oval[1] = 255 - oval[1];
+				oval[2] = 255 - oval[2];
+			} else {
+				oval[0] = oval[1] = oval[2] = 0;
+			}
+		}
+	} else {
+		x -= s->pw/2;
+		if (y < s->ph/2) {
+			/* Generate horizontal white stripes */
+			/* Stripes are 5 pixels apart in groups of 2 of the same level, */
+			/* declining by 1 each group */
+			if ((y % SG) == 0) {
+				j = y / (2 * SG);
+				v = 255 - j;
+				oval[0] = oval[1] = oval[2] = v;
+			} else {
+				oval[0] = oval[1] = oval[2] = 0;
+			}
+		} else {
+			y -= s->ph/2;
+			/* Generate black dots */
+			if ((x % SG) == 0 && (y % SG) == 0) {
+				i = x / (2 * SG);
+				j = y / (2 * SG);
+				i += j * s->pw/2/(2 * SG);
+				oval[0] = oval[1] = oval[2] = 0;
+
+				for (j = 5; j >= 0; j--) {
+					oval[0] = (oval[0] << 1) | ((i >> (j * 3 + 0)) & 1);
+					oval[1] = (oval[1] << 1) | ((i >> (j * 3 + 1)) & 1);
+					oval[2] = (oval[2] << 1) | ((i >> (j * 3 + 2)) & 1);
+				}
+			} else {
+				oval[0] = oval[1] = oval[2] = 255;
+			}
+		}
+	}
+
+	if (s->dpth == bpc8_2d) {
+		unsigned char *p = ((unsigned char *)outbuf) + xx * s->ncc;
+		p[0] = oval[0];
+		p[1] = oval[1];
+		p[2] = oval[2];
+	} else {
+		unsigned short *p = ((unsigned short *)outbuf) + xx * s->ncc;
+		p[0] = oval[0] * 256;
+		p[1] = oval[1] * 256;
+		p[2] = oval[2] * 256;
+	}
+}
+#undef SG
+#endif
+
+/* ------------------------------------------------------------- */
 /* Main class implementation */
 
 /* Free ourselves and all primitives */
@@ -156,9 +294,16 @@ static int colordiff(render2d *s, color2d c1, color2d c2) {
 #define MIXPOW 1.3
 #define OSAMLS 16
 
-/* Render and write to a TIFF file */
+/* Render and write to a TIFF or PNG file or memory buffer */
 /* Return NZ on error */
-static int render2d_write(render2d *s, char *filename, int comprn) {
+static int render2d_write(
+	render2d *s,
+	char *filename,			/* Name of file for file output */
+	int comprn,				/* nz to use compression */
+	unsigned char **obuf,	/* pointer to returned buffer for mem output. Free after use */
+	size_t *olen,			/* pointer to returned length of data in buffer */
+	rend_format fmt			/* Output format, tiff/png, file/memory */
+) {
 	TIFF *wh = NULL;
 	uint16 samplesperpixel = 0, bitspersample = 0;
 	uint16 extrasamples = 0;	/* Extra "alpha" samples */
@@ -166,9 +311,19 @@ static int render2d_write(render2d *s, char *filename, int comprn) {
 	uint16 photometric = 0;
 	uint16 inkset = 0xffff;
 	char *inknames = NULL;
-	tdata_t *outbuf;
-	unsigned char *tempbuf = NULL;		/* 16 bit buffer for dithering */
+	tdata_t *outbuf = NULL;
+
+	FILE *png_fp = NULL;
+	png_mem_info png_minfo = { NULL, 0, 0 };
+	png_structp png_ptr = NULL;
+	png_infop png_info = NULL;
+	png_uint_32 png_width = 0, png_height = 0;
+	int png_bit_depth = 0, png_color_type = 0;
+	int png_samplesperpixel = 0;
+
+	unsigned char *dithbuf16 = NULL;	/* 16 bit buffer for dithering */
 	thscreens *screen = NULL;			/* dithering object */
+	int foundfg;						/* Found a forground object in this line */
 	prim2d *th, **pthp;
 	prim2d **xlist, **ylist;	/* X, Y sorted start lists */
 	int xli, yli;				/* Indexes into X, Y list */
@@ -181,96 +336,232 @@ static int render2d_write(render2d *s, char *filename, int comprn) {
 	double rx0, rx1, ry0, ry1;	/* Box being processed, newest sample is rx1, ry1 */
 	int x, y;					/* Pixel x & y index */
 
+#ifdef CCTEST_PATTERN		// For testing by making screen visible
+#pragma message("######### render.c TEST_PATTERN defined ! ##")
+	int do_test_pattern = 0;
+
+	if (getenv("ARGYLL_CCAST_TEST_PATTERN") != NULL) {
+		verbose(0, "Substituting ChromeCast test pattern\n");
+		do_test_pattern = 1;
+	} else { 
+		static int verbed = 0;
+		if (!verbed) {
+			verbose(0, "Set ARGYLL_CCAST_TEST_PATTERN to enable test pattern\n");
+			verbed = 1;
+		}
+	}
+#endif
+
 	if ((so = new_sobol(2)) == NULL)
 		return 1;
 
-	switch (s->csp) {
-		case w_2d:			/* Video style grey */
-			samplesperpixel = 1;
-			photometric = PHOTOMETRIC_MINISBLACK;
-			break;
-		case k_2d:			/* Printing style grey */
-			samplesperpixel = 1;
-			photometric = PHOTOMETRIC_MINISWHITE;
-			break;
-		case lab_2d:		/* TIFF CIE L*a*b* */
-			samplesperpixel = 3;
-			photometric = PHOTOMETRIC_CIELAB;
-			break;
-		case rgb_2d:		/* RGB */
-			samplesperpixel = 3;
-			photometric = PHOTOMETRIC_RGB;
-			break;
-		case cmyk_2d:		/* CMYK */
-			samplesperpixel = 4;
-			photometric = PHOTOMETRIC_SEPARATED;
-			inkset = INKSET_CMYK;
-			inknames = "cyan\000magenta\000yellow\000\000";
-			break;
-		case ncol_2d:		/* N color */
-			samplesperpixel = s->ncc;
-			extrasamples = 0;
-			photometric = PHOTOMETRIC_SEPARATED;
-			inkset = 0;			// ~~99 should fix this
-			inknames = NULL;	// ~~99 should fix this
-			break;
-		case ncol_a_2d:		/* N color with extras in alpha */
-			samplesperpixel = s->ncc;
-			extrasamples = 0;
-			if (samplesperpixel > 4) {
-				extrasamples = samplesperpixel - 4;	/* Call samples > 4 "alpha" samples */
-				for (j = 0; j < extrasamples; j++)
-					extrainfo[j] = EXTRASAMPLE_UNASSALPHA;
+	if (fmt == tiff_file) {
+		switch (s->csp) {
+			case w_2d:			/* Video style grey */
+				samplesperpixel = 1;
+				photometric = PHOTOMETRIC_MINISBLACK;
+				break;
+			case k_2d:			/* Printing style grey */
+				samplesperpixel = 1;
+				photometric = PHOTOMETRIC_MINISWHITE;
+				break;
+			case lab_2d:		/* TIFF CIE L*a*b* */
+				samplesperpixel = 3;
+				photometric = PHOTOMETRIC_CIELAB;
+				break;
+			case rgb_2d:		/* RGB */
+				samplesperpixel = 3;
+				photometric = PHOTOMETRIC_RGB;
+				break;
+			case cmyk_2d:		/* CMYK */
+				samplesperpixel = 4;
+				photometric = PHOTOMETRIC_SEPARATED;
+				inkset = INKSET_CMYK;
+				inknames = "cyan\000magenta\000yellow\000\000";
+				break;
+			case ncol_2d:		/* N color */
+				samplesperpixel = s->ncc;
+				extrasamples = 0;
+				photometric = PHOTOMETRIC_SEPARATED;
+				inkset = 0;			// ~~99 should fix this
+				inknames = NULL;	// ~~99 should fix this
+				break;
+			case ncol_a_2d:		/* N color with extras in alpha */
+				samplesperpixel = s->ncc;
+				extrasamples = 0;
+				if (samplesperpixel > 4) {
+					extrasamples = samplesperpixel - 4;	/* Call samples > 4 "alpha" samples */
+					for (j = 0; j < extrasamples; j++)
+						extrainfo[j] = EXTRASAMPLE_UNASSALPHA;
+				}
+				photometric = PHOTOMETRIC_SEPARATED;
+				inkset = 0;			// ~~99 should fix this
+				inknames = NULL;	// ~~99 should fix this
+				break;
+			default:
+				error("render2d: Illegal colorspace for TIFF file '%s'",filename);
+		}
+		if (samplesperpixel != s->ncc)
+			error("render2d: mismatched number of color components");
+
+		switch (s->dpth) {
+			case bpc8_2d:		/* 8 bits per component */
+				bitspersample = 8;
+				break;
+			case bpc16_2d:		/* 16 bits per component */
+				bitspersample = 16;
+				break;
+			default:
+				error("render2d: Illegal bits per component for TIFF file '%s'",filename);
+		}
+
+		if ((wh = TIFFOpen(filename, "w")) == NULL)
+			error("render2d: Can\'t create TIFF file '%s'!",filename);
+		
+		TIFFSetField(wh, TIFFTAG_IMAGEWIDTH,  s->pw);
+		TIFFSetField(wh, TIFFTAG_IMAGELENGTH, s->ph);
+		TIFFSetField(wh, TIFFTAG_ORIENTATION, ORIENTATION_TOPLEFT);
+		TIFFSetField(wh, TIFFTAG_SAMPLESPERPIXEL, samplesperpixel);
+		TIFFSetField(wh, TIFFTAG_BITSPERSAMPLE, bitspersample);
+		TIFFSetField(wh, TIFFTAG_PLANARCONFIG, PLANARCONFIG_CONTIG);
+		TIFFSetField(wh, TIFFTAG_PHOTOMETRIC, photometric);
+		if (extrasamples > 0)
+			TIFFSetField(wh, TIFFTAG_EXTRASAMPLES, extrasamples, extrainfo);
+
+		if (inknames != NULL) {
+			int inlen = zzstrlen(inknames);
+			TIFFSetField(wh, TIFFTAG_INKSET, inkset);
+			TIFFSetField(wh, TIFFTAG_INKNAMES, inlen, inknames);
+		}
+		TIFFSetField(wh, TIFFTAG_COMPRESSION, COMPRESSION_NONE);
+		TIFFSetField(wh, TIFFTAG_RESOLUTIONUNIT, RESUNIT_CENTIMETER);
+		TIFFSetField(wh, TIFFTAG_XRESOLUTION, 10.0 * s->hres);	/* Cvt. to pixels/cm */
+		TIFFSetField(wh, TIFFTAG_YRESOLUTION, 10.0 * s->vres);
+		TIFFSetField(wh, TIFFTAG_XPOSITION, 0.1 * s->lm);		/* Cvt. to cm */
+		TIFFSetField(wh, TIFFTAG_YPOSITION, 0.1 * s->tm);
+		if (comprn) {
+			TIFFSetField(wh, TIFFTAG_COMPRESSION, COMPRESSION_LZW);
+		}
+		TIFFSetField(wh, TIFFTAG_IMAGEDESCRIPTION, "Test chart created with Argyll");
+
+		/* Allocate one TIFF line buffer */
+		outbuf = _TIFFmalloc(TIFFScanlineSize(wh));
+
+	} else if (fmt == png_file
+	        || fmt == png_mem) {
+		char *nmode = "w";
+
+#if !defined(O_CREAT) && !defined(_O_CREAT)
+# error "Need to #include fcntl.h!"
+#endif
+#if defined(O_BINARY) || defined(_O_BINARY)
+		nmode = "wb";
+#endif
+		png_width = s->pw;
+		png_height = s->ph;
+
+		switch (s->dpth) {
+			case bpc8_2d:		/* 8 bits per component */
+				png_bit_depth = 8;
+				break;
+			case bpc16_2d:		/* 16 bits per component */
+				png_bit_depth = 16;
+				break;
+			default:
+				error("render2d: Illegal bits per component for PNG file '%s'",filename);
+		}
+
+		switch (s->csp) {
+			case w_2d:			/* Video style grey */
+				png_color_type = PNG_COLOR_TYPE_GRAY;
+				png_samplesperpixel = 1;
+				break;
+			case rgb_2d:		/* RGB */
+				png_color_type = PNG_COLOR_TYPE_RGB;
+				png_samplesperpixel = 3;
+				break;
+			default:
+				error("render2d: Illegal colorspace for PNG file '%s'",filename);
+		}
+
+		if (fmt == png_file) {
+			if ((png_fp = fopen(filename, nmode)) == NULL)
+				error("render2d: Can\'t create PNG file '%s'!",filename);
+		}
+
+		if ((png_ptr = png_create_write_struct(PNG_LIBPNG_VER_STRING,
+		                                    NULL, NULL, NULL)) == NULL)
+			error("render2d: png_create_write_struct failed");
+
+		if (fmt == png_file) {
+			png_init_io(png_ptr, png_fp);
+		}
+
+		if ((png_info = png_create_info_struct(png_ptr)) == NULL) {
+			png_destroy_write_struct(&png_ptr, &png_info);
+			error("render2d: png_create_info_struct failed");
+		}
+
+	   if (setjmp(png_jmpbuf(png_ptr))) {
+			a1loge(g_log, 1, "%s -> %s: render2d libpng write error\n", filename);
+			png_destroy_info_struct(png_ptr, &png_info);
+			png_destroy_write_struct(&png_ptr, &png_info);
+			if (png_fp != NULL)
+				fclose(png_fp);
+			else {
+				free(png_minfo.buf);
 			}
-			photometric = PHOTOMETRIC_SEPARATED;
-			inkset = 0;			// ~~99 should fix this
-			inknames = NULL;	// ~~99 should fix this
-			break;
-		default:
-			error("render2d: Illegal colorspace for file '%s'",filename);
-	}
-	if (samplesperpixel != s->ncc)
-		error("render2d: mismatched number of color components");
+			return (1);
+		}
 
-	switch (s->dpth) {
-		case bpc8_2d:		/* 8 bits per component */
-			bitspersample = 8;
-			break;
-		case bpc16_2d:		/* 16 bits per component */
-			bitspersample = 16;
-			break;
-		default:
-			error("render2d: Illegal bits per component for file '%s'",filename);
-	}
+		if (fmt == png_mem) {
+			png_set_write_fn(png_ptr, &png_minfo, mem_write_data, mem_flush_data);  
+		}
 
-	if ((wh = TIFFOpen(filename, "w")) == NULL)
-		error("render2d: Can\'t create TIFF file '%s'!",filename);
-	
-	TIFFSetField(wh, TIFFTAG_IMAGEWIDTH,  s->pw);
-	TIFFSetField(wh, TIFFTAG_IMAGELENGTH, s->ph);
-	TIFFSetField(wh, TIFFTAG_ORIENTATION, ORIENTATION_TOPLEFT);
-	TIFFSetField(wh, TIFFTAG_SAMPLESPERPIXEL, samplesperpixel);
-	TIFFSetField(wh, TIFFTAG_BITSPERSAMPLE, bitspersample);
-	TIFFSetField(wh, TIFFTAG_PLANARCONFIG, PLANARCONFIG_CONTIG);
-	TIFFSetField(wh, TIFFTAG_PHOTOMETRIC, photometric);
-	if (extrasamples > 0)
-		TIFFSetField(wh, TIFFTAG_EXTRASAMPLES, extrasamples, extrainfo);
-
-	if (inknames != NULL) {
-		int inlen = zzstrlen(inknames);
-		TIFFSetField(wh, TIFFTAG_INKSET, inkset);
-		TIFFSetField(wh, TIFFTAG_INKNAMES, inlen, inknames);
-	}
-	TIFFSetField(wh, TIFFTAG_COMPRESSION, COMPRESSION_NONE);
-	TIFFSetField(wh, TIFFTAG_RESOLUTIONUNIT, RESUNIT_CENTIMETER);
-	TIFFSetField(wh, TIFFTAG_XRESOLUTION, 10.0 * s->hres);	/* Cvt. to pixels/cm */
-	TIFFSetField(wh, TIFFTAG_YRESOLUTION, 10.0 * s->vres);
-	TIFFSetField(wh, TIFFTAG_XPOSITION, 0.1 * s->lm);		/* Cvt. to cm */
-	TIFFSetField(wh, TIFFTAG_YPOSITION, 0.1 * s->tm);
-	if (comprn) {
-		TIFFSetField(wh, TIFFTAG_COMPRESSION, COMPRESSION_LZW);
+		png_set_IHDR(png_ptr, png_info, png_width, png_height, png_bit_depth,
+			png_color_type, PNG_INTERLACE_NONE, PNG_COMPRESSION_TYPE_BASE, PNG_FILTER_TYPE_BASE);
+
+
+		/* pix/mm to pix/meter */
+		png_set_pHYs(png_ptr, png_info, (png_uint_32)(1000.0 * s->hres + 0.5),
+		                                (png_uint_32)(1000.0 * s->vres + 0.5),
+		                                PNG_RESOLUTION_METER);
+		
+		/* mm to um */
+		png_set_oFFs(png_ptr, png_info, (png_uint_32)(1000.0 * 0.1 * s->lm),  
+		            (png_uint_32)(1000.0 * s->tm), PNG_OFFSET_MICROMETER);
+
+		{
+			png_text txt;
+			txt.compression = PNG_TEXT_COMPRESSION_NONE;
+			txt.key = "Description";
+			txt.text = "Test chart created with Argyll";
+			txt.text_length = strlen(txt.text);
+#ifdef PNG_iTXt_SUPPORTED
+			txt.itxt_length = 0;
+			txt.lang = NULL;
+			txt.lang_key = NULL;
+#endif
+			png_set_text(png_ptr, png_info, &txt, 1);
+		}
+
+		/* Write the header */
+		png_write_info(png_ptr, png_info);
+
+		/* PNG expects network (BE) 16 bit values */
+		if (png_bit_depth == 16) {
+			png_uint_16 tt = 0x0001;
+			if (*((unsigned char *)&tt) == 0x01) 	/* Little endian */
+				png_set_swap(png_ptr);
+		}
+
+		/* Allocate one PNG line buffer */
+		if ((outbuf = malloc((png_bit_depth >> 3) * png_samplesperpixel * s->pw) ) == NULL)
+			error("malloc of PNG line buffer failed");
+
+	} else {
+		error("render2d: Illegal output format %d",fmt);
 	}
-	TIFFSetField(wh, TIFFTAG_IMAGEDESCRIPTION, "Test chart created with Argyll");
 
 	/* Allocate pixel value storage for aliasing detection */
 	if ((_pixv0 = malloc(sizeof(color2d) * (s->pw+2))) == NULL)
@@ -280,23 +571,24 @@ static int render2d_write(render2d *s, char *filename, int comprn) {
 		return 1;
 	pixv1 = _pixv1+1;
 
-	/* Allocate one TIFF line buffer */
-	outbuf = _TIFFmalloc(TIFFScanlineSize(wh));
-
 	if (s->dpth == bpc8_2d && s->dither) {
-#ifdef NEVER		// For testing by making screen visible
-# define LEVELS 16
+#ifdef TEST_SCREENING		// For testing by making screen visible
+#pragma message("######### render TEST_SCREENING defined! ##")
+# define LEVELS 4
 		int i, olevs[LEVELS];
-		for (i = 0; i < LEVELS; i++)
+		for (i = 0; i < LEVELS; i++) {
 			olevs[i] = (int)(i/(LEVELS-1.0) * 255.0 + 0.5);
+		}
 		if ((screen = new_thscreens(0, s->ncc, 1.0, 79, scie_16, 8, LEVELS, olevs,
-		                           scoo_l, 0.1, NULL, NULL)) == NULL)
+		                           scoo_l, OVERLAP, s->pw, NULL, NULL,
+		                           s->dither == 2 ? 1 : 0, s->quant, s->qcntx)) == NULL)
 #else
 		if ((screen = new_thscreens(0, s->ncc, 1.0, 79, scie_16, 8, 256, NULL,
-		                           scoo_l, 0.1, NULL, NULL)) == NULL)
+		                           scoo_l, OVERLAP, s->pw, NULL, NULL,
+		                           s->dither == 2 ? 1 : 0, s->quant, s->qcntx)) == NULL)
 #endif
 			return 1;
-		if ((tempbuf = malloc(s->pw * s->ncc * 2)) == NULL)
+		if ((dithbuf16 = malloc(s->pw * s->ncc * 2)) == NULL)
 			return 1;
 	}
 
@@ -323,12 +615,14 @@ static int render2d_write(render2d *s, char *filename, int comprn) {
 	yli = 0;
 	s->yl = NULL;
 
-	/* Render each line and write it. */
+	/* Render each line and write it, in raster order. */
 	/* We sample +- half a pixel around the pixel we want. */
 	/* We make the active element list encompass this region, */
 	/* so that we can super sample it for anti-aliasing. */
 	for (y = -1; y < s->ph; y++) {
+		foundfg = 0;
 
+		/* Convert to coordinate order */
 		ry0 = (((s->ph-1) - y) - 0.5) / s->vres;
 		ry1 = (((s->ph-1) - y) + 0.5) / s->vres;
 
@@ -389,6 +683,8 @@ static int render2d_write(render2d *s, char *filename, int comprn) {
 					*pthp = th->xl;
 					th = th->xl;
 				} else {
+//printf("x %d y %d, rx1 %f, ry0 %f\n",x,y,rx1, ry0);
+
 					if (th->rend(th, rv, rx1, ry0) && th->ix > pixv1[x][PRIX2D]) {
 						/* Overwrite the current color */
 						/* (This is where we should handle depth and opacity */
@@ -400,7 +696,7 @@ static int render2d_write(render2d *s, char *filename, int comprn) {
 					th = th->xl;
 				}
 			}
-			/* Check if anti-aliasing is neded for previous lines previous pixel */
+			/* Check if anti-aliasing is needed for previous lines previous pixel */
 			if (y >= 0 && x >= 0) {
 				color2d cc;
 
@@ -409,15 +705,17 @@ static int render2d_write(render2d *s, char *filename, int comprn) {
 				cc[PRIX2D] = pixv1[x][PRIX2D];
 
 				/* See if anti aliasing is needed */
-				if ((pixv0[x+0][PRIX2D] != cc[PRIX2D] && colordiff(s, pixv0[x+0], cc))
-				 || (pixv0[x-1][PRIX2D] != cc[PRIX2D] && colordiff(s, pixv0[x-1], cc))
-				 || (pixv1[x-1][PRIX2D] != cc[PRIX2D] && colordiff(s, pixv1[x-1], cc))) {
+				if (!s->noavg
+				 && ((pixv0[x+0][PRIX2D] != cc[PRIX2D] && colordiff(s, pixv0[x+0], cc))
+				  || (pixv0[x-1][PRIX2D] != cc[PRIX2D] && colordiff(s, pixv0[x-1], cc))
+				  || (pixv1[x-1][PRIX2D] != cc[PRIX2D] && colordiff(s, pixv1[x-1], cc)))) {
 					double nn = 0;
 
 					so->reset(so);
 
 					for (j = 0; j < s->ncc; j++)
 						cc[j] = 0.0;
+					cc[PRIX2D] = -1;
 
 					/* Compute the sample value by re-sampling the region */
 					/* around the pixel. */
@@ -447,6 +745,8 @@ static int render2d_write(render2d *s, char *filename, int comprn) {
 						}
 						for (j = 0; j < s->ncc; j++)
 							cc[j] += pow(ccc[j], MIXPOW);
+						if (ccc[PRIX2D] > cc[PRIX2D])
+							cc[PRIX2D] = ccc[PRIX2D];	/* Note if not BG */
 					}
 					for (j = 0; j < s->ncc; j++)
 						cc[j] = pow(cc[j]/nn, 1.0/MIXPOW);
@@ -456,6 +756,12 @@ static int render2d_write(render2d *s, char *filename, int comprn) {
 					cc[1] = 0.0;
 					cc[2] = 1.0;
 #endif
+				} else if (s->noavg) {
+					/* Compute output value directly from primitive */
+
+					for (j = 0; j < s->ncc; j++)
+						cc[j] = cc[j];
+
 				} else {
 
 					/* Compute output value as mean of surrounding samples */
@@ -465,14 +771,24 @@ static int render2d_write(render2d *s, char *filename, int comprn) {
 						      + pixv0[x][j]
 						      + pixv1[x-1][j];
 						cc[j] = cc[j] * 0.25;
-							
 					}
+					/* Note if not BG */
+					if (pixv0[x-1][PRIX2D] > cc[PRIX2D])
+						cc[PRIX2D] = pixv0[x-1][PRIX2D];
+					if (pixv0[x][PRIX2D] > cc[PRIX2D])
+						cc[PRIX2D] = pixv0[x][PRIX2D];
+					if (pixv1[x-1][PRIX2D] > cc[PRIX2D])
+						cc[PRIX2D] = pixv1[x-1][PRIX2D];
 				}
+				if (cc[PRIX2D] != -1)		/* Line is no longer background */
+					foundfg = 1;
 
 				/* Translate from render value to output pixel value */
 				if (s->dpth == bpc8_2d) {
+					/* if dithering and dithering all or found FG in line */
 					if (s->dither) {
-						unsigned short *p = ((unsigned short *)tempbuf) + x * s->ncc;
+						unsigned short *p = ((unsigned short *)dithbuf16) + x * s->ncc;
+
 						if (s->csp == lab_2d) {
 							cvt_Lab_to_CIELAB16(cc, cc);
 							for (j = 0; j < s->ncc; j++)
@@ -507,12 +823,81 @@ static int render2d_write(render2d *s, char *filename, int comprn) {
 		}
 
 		if (y >= 0) {
-			if (s->dpth == bpc8_2d && s->dither)
-				screen->screen(screen, s->pw, 1, 0, y, tempbuf, s->pw * s->ncc * 2,
-				                       (unsigned char *)outbuf, s->pw * s->ncc);
+			/* if dithering and dithering all or found FG in line */
+			if (s->dpth == bpc8_2d && s->dither) {
+				// If we need to screen this line
+				if (!s->dithfgo || foundfg) {
+					/* If we are dithering only the foreground colors, */
+					/* Subsitute the quantized un-dithered color for any */
+					/* pixels soley from the background */
+					if (s->dithfgo) {
+						int st, ed;
+						unsigned short *ip = ((unsigned short *)dithbuf16);
+						unsigned char *op = ((unsigned char *)outbuf);
+
+						/* Copy pixels up to first non-BG */
+						for (x = 0; x < s->pw; x++, ip += s->ncc, op += s->ncc) {
+							if (pixv1[x][PRIX2D] == -1) {
+								for (j = 0; j < s->ncc; j++)
+									op[j] = (ip[j] * 255 + 128)/65535;
+							} else {
+								st = x;
+								break;
+							}
+						}
+						if (x < s->pw) {	/* If there are FG pixels */
+
+							/* Copy down to first non-BG */
+							ip = ((unsigned short *)dithbuf16) + (s->pw-1) * s->ncc;
+							op = ((unsigned char *)outbuf) + (s->pw-1) * s->ncc;
+							for (x = s->pw-1; x > st; x--, ip -= s->ncc, op -= s->ncc) {
+								if (pixv1[x][PRIX2D] == -1) {
+									for (j = 0; j < s->ncc; j++)
+										op[j] = (ip[j] * 255 + 128)/65535;
+								} else {
+									ed = x;
+									break;
+								}
+							}
+							/* Screen just the FG pixels */
+							ip = ((unsigned short *)dithbuf16) + st * s->ncc;
+							op = ((unsigned char *)outbuf) + st * s->ncc;
+							screen->screen(screen, ed-st+1, 1, st, y,
+						                       op, s->pw * s->ncc,
+						                       (unsigned char*)ip, s->pw * s->ncc);
+						}
+
+					/* Dither/screen the whole lot */
+					} else {
+						screen->screen(screen, s->pw, 1, 0, y,
+						                       (unsigned char *)outbuf, s->pw * s->ncc,
+						                       dithbuf16, s->pw * s->ncc);
+					}
+				// Don't need to screen this line - quantize from 16 bit
+				} else {
+					unsigned short *ip = ((unsigned short *)dithbuf16);
+					unsigned char *op = ((unsigned char *)outbuf);
+					for (x = 0; x < s->pw; x++, ip += s->ncc, op += s->ncc) {
+						for (j = 0; j < s->ncc; j++)
+							op[j] = (ip[j] * 255 + 128)/65535;
+					}
+				}
+			}
 
-			if (TIFFWriteScanline(wh, outbuf, y, 0) < 0)
-				error ("Failed to write TIFF file '%s' line %d",filename,y);
+#ifdef CCTEST_PATTERN		// Substitute the testing pattern
+			if (do_test_pattern) {
+				for (x = 0; x < s->pw; x++)
+					test_value(s, outbuf, x, y);
+			}
+#endif
+			if (fmt == tiff_file) {
+				if (TIFFWriteScanline(wh, outbuf, y, 0) < 0)
+					error ("Failed to write TIFF file '%s' line %d",filename,y);
+			} else if (fmt == png_file
+			        || fmt == png_mem) {
+				png_bytep pixdata = (png_bytep)outbuf;
+				png_write_rows(png_ptr, &pixdata, 1);
+			}
 		}
 
 		/* Shuffle the pointers */
@@ -529,12 +914,29 @@ static int render2d_write(render2d *s, char *filename, int comprn) {
 	free(_pixv0);
 	free(_pixv1);
 
-	if (tempbuf != NULL)
-		free(tempbuf);
+	if (dithbuf16 != NULL)
+		free(dithbuf16);
 	if (screen != NULL)
 		screen->del(screen);
-	_TIFFfree(outbuf);
-	TIFFClose(wh);		/* Close Output file */
+
+	if (fmt == tiff_file) {
+		_TIFFfree(outbuf);
+		TIFFClose(wh);		/* Close Output file */
+
+	} else if (fmt == png_file
+	        || fmt == png_mem) {
+
+		free(outbuf);
+		png_write_end(png_ptr, NULL);
+//		png_destroy_info_struct(png_ptr, &png_info);
+		png_destroy_write_struct(&png_ptr, &png_info);
+		if (fmt == png_file) {
+			fclose(png_fp);
+		} else if (fmt == png_mem) {
+			*obuf = png_minfo.buf;
+			*olen = png_minfo.off;
+		}
+	}
 
 	so->del(so);
 
@@ -551,7 +953,9 @@ double vres,	/* horizontal resolution in pixels/mm */
 colort2d csp,	/* Color type */
 int nd,			/* Number of channels if c = ncol */
 depth2d dpth,	/* Pixel depth */
-int dither		/* Dither flag */
+int dither,		/* Dither flag, 1 = ordered, 2 = error diffusion, | 0x8000 to dither FG only  */
+void (*quant)(void *qcntx, double *out, double *in), /* optional quantization func. for edith */
+void *qcntx
 ) {
 	render2d *s;
 
@@ -577,7 +981,11 @@ int dither		/* Dither flag */
 	s->vres = vres;
 	s->csp = csp;
 	s->dpth = dpth;
-	s->dither = dither;
+	s->dither  = 0x0fff & dither;
+	s->noavg   = 0x4000 & dither;
+	s->dithfgo = 0x8000 & dither;
+	s->quant = quant;
+	s->qcntx = qcntx;
 
 	s->del = render2d_del;
 	s->set_defc = render2d_set_defc;
@@ -637,13 +1045,31 @@ static int rect2d_rend(prim2d *ss, color2d rv, double x, double y) {
 	 || x < s->rx0 || x > s->rx1)
 		return 0;
 
-	for (j = 0; j < s->ncc; j++)
-		rv[j] = s->c[j];
+	if (s->dpat == NULL) {
+		for (j = 0; j < s->ncc; j++)
+			rv[j] = s->c[j];
+
+	/* We have a dither pattern */
+	} else {
+		int xi = ((int)floor(x)) % s->dp_w;
+		int yi = ((int)floor(y)) % s->dp_h;
+		double *val = (*s->dpat)[xi][yi];
+		for (j = 0; j < s->ncc; j++)
+			rv[j] = val[j];
+	}
+
 	rv[PRIX2D] = s->ix;
 
 	return 1;
 }
 
+static void rect2d_del(prim2d *ss) {
+	rect2d *s = (rect2d *)ss;
+	if (s->dpat != NULL)
+		free(s->dpat);
+	prim2d_del(ss);
+}
+
 prim2d *new_rect2d(
 render2d *ss,
 double x,
@@ -664,7 +1090,7 @@ color2d c
 	y -= ss->bm;
 
 	s->ncc = ss->ncc;
-	s->del = prim2d_del; 
+	s->del = rect2d_del; 
 	s->rend = rect2d_rend; 
 
 	/* Set bounding box */
@@ -685,6 +1111,13 @@ color2d c
 	return (prim2d *)s;
 }
 
+/* Allocate pat using malloc(sizeof(double) * MXPATSIZE * MXPATSIZE * TOTC2D) */
+void set_rect2d_dpat(rect2d *s, double (*pat)[MXPATSIZE][MXPATSIZE][TOTC2D], int w, int h) {
+	s->dpat = pat;
+	s->dp_w = w;
+	s->dp_h = h;
+}
+
 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 /* Vertex shaded rectangle */
 
@@ -1039,7 +1472,8 @@ color2d c
 
 
 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
-/* Primitive Macros. */
+/* Primitive Macros. These shapes are composed of */
+/* underlying primitives */
 
 /* add a dashed line */
 void add_dashed_line2d(
diff --git a/render/render.h b/render/render.h
index e7bcdee..c060ebc 100644
--- a/render/render.h
+++ b/render/render.h
@@ -10,7 +10,7 @@
  * Author:  Graeme W. Gill
  * Date:    28/12/2005
  *
- * Copyright 2005, 2008, 2012 Graeme W. Gill
+ * Copyright 2005, 2008, 2012, 2014 Graeme W. Gill
  * All rights reserved.
  *
  * This material is licenced under the GNU AFFERO GENERAL PUBLIC LICENSE Version 3 :-
@@ -28,7 +28,9 @@
 
 #define MXCH2D  8			/* Maximum color channels */
 #define TOTC2D  (MXCH2D+1)	/* Maximum total components */
-#define PRIX2D  (MXCH2D)	/* Index of primitive */
+#define PRIX2D  (MXCH2D)	/* Index of primitive kept with color value */
+
+#define MXPATSIZE 4			/* Maximum color pattern size */
 
 /* Color type */
 /* Shouldn't this be an xcolorants mask ? */
@@ -55,7 +57,7 @@ typedef enum {
 	rowman_s = 0,	/* Rownman, single stroke */
 	rowman_d = 1,	/* Rownman, double stroke */
 	rowman_t = 2,	/* Rownman, triple stroke */
-	timesr = 3,		/* Times Roman */
+	timesr   = 3,	/* Times Roman */
 	timesr_b = 4,	/* Times Roman, Bold */
 	futura_l = 5,	/* Futura, Light */
 	futura_m = 6	/* Futura, Medium */
@@ -87,11 +89,15 @@ struct _prim2d {
 struct _rect2d {
 	PRIM_STRUCT
 	double rx0, ry0, rx1, ry1;	/* Rectangle verticies */
-	color2d c;			/* Color of rectangle */
+	color2d c;					/* Color of rectangle (if dpat == NULL) */
+	double (*dpat)[MXPATSIZE][MXPATSIZE][TOTC2D];
+	int dp_w, dp_h;
 }; typedef struct _rect2d rect2d;
 
 prim2d *new_rect2d(struct _render2d *s, double x, double y, double w, double h, color2d c);
 
+void set_rect2d_dpat(struct _rect2d *s, double (*pat)[MXPATSIZE][MXPATSIZE][TOTC2D], int w, int h);
+
 /* ------------------------------------ */
 /* Vertex shaded rectange */
 struct _rectvs2d {
@@ -145,7 +151,7 @@ color2d c);
 /* Add a text character at the given location using lines */
 void add_char2d(
 struct _render2d *s,
-double *xinc,		/* Add increment to next character */
+double *xinc,		/* Return increment in position for next character */
 double *yinc,
 font2d fo,			/* Font to use */
 char ch,			/* Character code to be printed */
@@ -158,7 +164,7 @@ color2d c			/* Color of text */
 /* Add a string from the given location using lines. */
 void add_string2d(
 struct _render2d *s,
-double *xinc,		/* Add increment to next character */
+double *xinc,		/* Return increment in position for next character */
 double *yinc,
 font2d fo,			/* Font to use */
 char *string,		/* Character code to be printed */
@@ -171,7 +177,7 @@ color2d c			/* Color of text */
 /* Return the total width of the string without adding it */
 void meas_string2d(
 struct _render2d *s,
-double *xinc,		/* Add increment to next character */
+double *xinc,		/* Return increment in position for next character */
 double *yinc,
 font2d fo,			/* Font to use */
 char *string,		/* Character code to be printed */
@@ -180,6 +186,15 @@ int or				/* Orintation, 0 = right, 1 = down, 2 = left, 3 = up */
 );
 
 /* ------------------------------------ */
+
+/* Type of output to save to. */
+typedef enum {
+	tiff_file,		/* Write a TIFF format file */
+	png_file,		/* Write a PNG format file */
+	png_mem			/* Write a PNG image to a memory buffer */
+} rend_format;
+
+/* ------------------------------------ */
 /* Render object */
 
 struct _render2d {
@@ -194,7 +209,11 @@ struct _render2d {
 	colort2d csp;			/* Color space */
 	int      ncc;			/* Number of color components */
 	depth2d  dpth;			/* Depth of the components */
-	int    dither;     		/* Dither flag */
+	int    dither;     		/* Dither flag, 1 = ordered, 2 = error diffusion */
+	int    noavg;     		/* Don't anti-alias or average 4 pixels together */
+	int    dithfgo;			/* Dither F.G. only flag */
+	void (*quant)(void *qcntx, double *out, double *in); /* optional quantization func. for edith */
+	void *qcntx;
 
 	color2d defc;			/* Default color value */
 
@@ -218,8 +237,10 @@ struct _render2d {
 
 	void (*add)(struct _render2d *s, prim2d *p);		/* Add a primitive */
 
-	int (*write)(struct _render2d *s, char *filename, int comprn);
-														/* Render and write to a TIFF file */
+	int (*write)(struct _render2d *s, char *filename, int comprn,
+		unsigned char **obuf, size_t *olen,
+	    rend_format fmt);
+												/* Render and write to a TIFF or PNG file */
 }; typedef struct _render2d render2d;
 
 /* Constructor */
@@ -232,7 +253,10 @@ render2d *new_render2d(
 	colort2d csp,	/* Color type */
 	int nd,			/* Number of channels if c = ncol */
 	depth2d dpth,	/* Pixel depth */
-	int dither		/* Dither flag */
+	int dither,		/* Dither flag, 1 = ordered, 2 = error diffusion, | 0x8000 to dither FG only */
+					/* | 0x4000 don't anti-alias by averaging pixels together. */
+	void (*quant)(void *qcntx, double *out, double *in), /* optional quantization func. for edith */
+	void *qcntx
 );
 
 #endif /* RENDER2D_H */
diff --git a/render/thscreen.c b/render/thscreen.c
index c0a5548..8523e38 100644
--- a/render/thscreen.c
+++ b/render/thscreen.c
@@ -2,7 +2,7 @@
 /*
  * render2d
  *
- * Threshold screen pixel processing object.
+ * Threshold or Error diffusion screen pixel processing object.
  * (Simplified from DPS code)
  *
  * Author:  Graeme W. Gill
@@ -32,6 +32,8 @@
 /* Configuration: */
 #undef DEBUG
 
+#undef CHECK_EXPECTED_ED_LEVELS		/* Output expected quantized levels for checkking */
+
 /* ----------------------------------------------------------- */
 
 #ifdef DEBUG
@@ -46,29 +48,190 @@
 
 #include "screens.h"	/* Pre-generated screen patterns */
 
-/* Screen a single color plane */
-void screen_thscreens(			/* Pointer to dither function */
+/* Threshold screen lines of multiplane pixels */
+void screen_thscreens(
 	thscreens *t,			/* Screening object pointer */
 	int width, int height,	/* Width and height to screen in pixels */
 	int xoff, int yoff,		/* Offset into screening pattern */
-	unsigned char *in,		/* Input pixel buffer */
-	unsigned long ipitch,	/* Increment between input lines */
 	unsigned char *out,		/* Output pixel buffer */
-	unsigned long opitch	/* Increment between output lines */
+	unsigned long opitch,	/* Increment between output lines in components */
+	unsigned char *in,		/* Input pixel buffer */
+	unsigned long ipitch	/* Increment between input lines in components */
 ) {
 	int i;
 	for (i = 0; i < t->np; i++)
-		t->sc[i]->screen(t->sc[i], width, height, xoff, yoff, in + 2 * i, t->np, ipitch,
-		                                                      out + i, t->np, opitch);
+		t->sc[i]->screen(t->sc[i], width, height, xoff, yoff,
+		                           out + i, t->np, opitch,
+		                           in + 2 * i, t->np, ipitch);
+}
+
+/* Error diffusion screen lines of multiplane pixels */
+void screen_edscreens(
+	thscreens *t,			/* Screening object pointer */
+	int width, int height,	/* Width and height to screen in pixels */
+	int xoff, int yoff,		/* Offset into screening pattern, [xoff + width < mxwidth] */
+	unsigned char *out,		/* Output pixel buffer */
+	unsigned long opitch,	/* Increment between output lines in components */
+	unsigned char *_in,		/* Input pixel buffer */
+	unsigned long ipitch	/* Increment between input lines in components */
+) {
+	unsigned short *in = (unsigned short *)_in;	/* Pointer to input pixel sized values */
+	unsigned short *ein = in + height * ipitch;	/* Vertical end pixel marker */
+	unsigned short *ein1;				/* Horizontal end pixel markers */
+	int xo, yo;				/* Threshold screen offset */
+	int x, j;
+
+	/* Limit width to mxwidth */
+	if ((xoff + width) > t->mxwidth) {
+		width = t->mxwidth - xoff;
+		if (width < 0)
+			return;
+	}
+
+	/* If not sequential, clear error buffer */
+	if (yoff != (t->lastyoff+1)) {
+		for (x = -1; x <= t->mxwidth; x++) {
+			for (j = 0; j < t->np; j++)
+				t->ebuf[j][x] = 0.0;
+		}
+	}
+
+	/* Clear "next to right" error */
+	for (j = 0; j < t->np; j++) {
+		t->ebuf[j][-2] = 0.0;
+	}
+
+	t->lastyoff = yoff;
+
+	/* For each line: */
+	for (; in < ein; in += ipitch, ein1 += ipitch, out += opitch, yoff++) {
+		unsigned short *ip;	/* Horizontal input pointer */
+		unsigned char *op;	/* Horizontal output pointer */
+		int xinc, pinc;
+
+		/* Do in serpentine order */
+		if (yoff & 1) {
+			xinc = -1;
+			x = xoff + width-1;			/* x is index into error buffer */
+			pinc = -t->np;
+			ein1 = in + pinc;
+			ip = in + t->np * (width-1);
+			op = out + t->np * (width-1);
+		} else {
+			xinc = 1;
+			x = xoff;
+			pinc = t->np;
+			ein1 = in + t->np * width;
+			ip = in;
+			op = out;
+		}
+
+		/* For each pixel */
+		for (; ip != ein1; ip += pinc, op += pinc, x += xinc) {
+			double ov[THMXCH2D], tv[THMXCH2D], ev[THMXCH2D];
+
+			/* For each plane */
+			for (j = 0; j < t->np; j++) {
+				tv[j] = t->luts[j][ip[j]] / 65535.0;		/* 0.0 - 1.0 value */
+
+				/* Value + accumulated error */
+				ov[j] = tv[j] = tv[j] + t->ebuf[j][x];
+
+				/* Limit */
+				if (ov[j] > 1.0)
+					ov[j] = 1.0;
+				else if (ov[j] < 0.0)
+					ov[j] = 0.0;
+
+				/* Output encode */
+				op[j] = t->oevalues[(int)(ov[j] * (t->oelev-1.0) + 0.5)];
+			}
+
+#ifdef CHECK_EXPECTED_ED_LEVELS
+#pragma message("######### render/thscreen.c CHECK_EXPECTED_ED_LEVELS defined ! ##")
+			// Put expected values in output to check levels
+			t->quant(t->qcntx, ev, ov);
+			for (j = 0; j < t->np; j++)
+				op[j] = t->oevalues[(int)(ev[j] * (t->oelev-1.0) + 0.5)];
+#endif
+
+			/* Quantize to values that it actually will be */
+			if (t->quant != NULL)
+				t->quant(t->qcntx, ov, ov);
+			else {
+				for (j = 0; j < t->np; j++)
+					ov[j] = floor(ov[j] * (t->oelev-1) + 0.5)/(t->oelev-1.0);
+			}
+
+			/* Compute the error to the target */
+			for (j = 0; j < t->np; j++) {
+
+				/* Error to target */
+				ev[j] = tv[j] - ov[j];
+			}
+
+			/* Distribute the error */
+			for (j = 0; j < t->np; j++) {
+#ifdef NEVER
+				/* Classic error diffusion */
+				t->ebuf[j][x-xinc] += 0.1875 * ev[j];				/* Lower left */
+				t->ebuf[j][x] = t->ebuf[j][-2] + 0.3125 * ev[j];	/* Lower */
+				t->ebuf[j][-2] = 0.0625 * ev[j];					/* Lower right */
+				t->ebuf[j][x+xinc] += 0.4375 * ev[j];				/* Right */
+#else
+				/* Using random placement error distribution */
+				double rav;
+				int ii;
+				t->so->next(t->so, &rav);		/* For some order */
+				rav *= 4.0;
+				rav += d_rand(0.0, 2.5);		/* For some randomness */
+				ii = (int)(rav);
+				if (ii > 3)
+				 	ii -= 4;
+				t->ebuf[j][x] = t->ebuf[j][-2];
+				t->ebuf[j][-2] = 0.0;
+				switch (ii) {
+					case 0:
+						t->ebuf[j][x-xinc] += ev[j];		/* Lower left */
+						break;
+					case 1:
+						t->ebuf[j][x] += ev[j];				/* Lower */
+						break;
+					case 2:
+						t->ebuf[j][-2] += ev[j];			/* Lower right */
+						break;
+					case 3:
+						t->ebuf[j][x+xinc] += ev[j];		/* Right */
+						break;
+				}
+#endif
+			}
+		}
+	}
 }
 
 /* Delete a thscreens */
 void del_thscreens(thscreens *t) {
 	int i;
 
-	for (i = 0; i < t->np; i++)
-		t->sc[i]->del(t->sc[i]);
-	free(t->sc);
+	if (t->sc != NULL) {
+		for (i = 0; i < t->np; i++) {
+			if (t->sc[i] != NULL)
+				t->sc[i]->del(t->sc[i]);
+		}
+		free(t->sc);
+	}
+	if (t->ebuf != NULL) {
+		free_fmatrix(t->ebuf, 0, t->np-1, -2, t->mxwidth);
+	}
+
+	if (t->luts != NULL) {
+		free_imatrix(t->luts, 0, t->np-1, 0, 65535);
+	}
+
+	if (t->so != NULL)
+		t->so->del(t->so);
+
 	free(t);
 }
 
@@ -77,7 +240,7 @@ thscreens *new_thscreens(
 	int exact,				/* Return only exact matches */
 	int nplanes,			/* Number of planes to screen */
 	double asp,				/* Target aspect ratio (== dpiX/dpiY) */
-	int size,				/* Target size */
+	int size,				/* Target screen size */
 	sc_iencoding ie,		/* Input encoding - must be scie_16 */
 	int oebpc,				/* Output encoding bits per component - must be 8 */
 	int oelev,				/* Output encoding levels. Must be <= 2 ^ oebpc */
@@ -85,8 +248,12 @@ thscreens *new_thscreens(
 							/* Must be oelev entries. Default is 0 .. oelev-1 */
 	sc_oorder oo,			/* Output bit ordering */
 	double overlap,			/* Overlap between levels, 0 - 1.0 */
+	int mxwidth,			/* max width in pixels of raster to be screened */ 
 	void   **cntx,			/* List of contexts for lookup table callback */
-	double (**lutfunc)(void *cntx, double in)	/* List of callback function, NULL if none */
+	double (**lutfunc)(void *cntx, double in),	/* List of callback functions, NULL if none */
+	int edif,				/* nz if using error diffusion */
+	void (*quant)(void *qcntx, double *out, double *in), /* optional quantization func. for edif */
+	void *qcntx
 ) {
 	thscreens *t;
 	int i, bi = -1;
@@ -115,18 +282,40 @@ thscreens *new_thscreens(
 	}
 
 	t->np = nplanes; 	/* Number of planes */
+	t->edif = edif;		/* Error diffusion */
+	t->quant = quant;	/* Optional quantization function */
+	t->qcntx = qcntx;
 
-	DBG(("thscreens no planes = %d\n",t->np));
+	t->mxwidth = mxwidth;
+	t->lastyoff = -1;
 
-	t->screen = screen_thscreens;
-	t->del = del_thscreens;
+	/* Allocate and initialise a next line error buffer. */
+	/* we allow 2 extra locations for pixels to the left and right of the current one: */
+	/* [-1] for the one to the below left when we are at x = 0, */
+	/* [-2] for the one below right, before we use [x] */
+	if (t->edif)
+		t->ebuf = fmatrixz(0, t->np-1, -2, t->mxwidth);
 
-	if ((t->sc = malloc(sizeof(thscreen *) * t->np)) == NULL) {
-		free(t);
-		DBG(("thscreens: malloc of thscreens->sc[] failed\n"));
-		return NULL;
+	t->oebpc = oebpc;
+	t->oelev = oelev;
+	if (oevalues != NULL) {
+		for (i = 0; i < t->oelev; i++) {
+			if (oevalues[i] >= (1 << t->oebpc)) {
+				DBG(("new_thscreens() oevalues[%d] value %d can't fit in %d bits\n",i,oevalues[i],t->oebpc));
+				free(t);
+				return NULL;
+			}
+			t->oevalues[i] = oevalues[i];
+		}
+	} else {
+		for (i = 0; i < t->oelev; i++)
+			t->oevalues[i] = i;
 	}
 
+	DBG(("thscreens no planes = %d\n",t->np));
+
+	t->del = del_thscreens;
+
 	DBG(("thscreens: searching amongst %d screens, exact = %d\n",NO_SCREENS,exact));
 
 	DBG(("thscreens: looking for non-exact match\n"));
@@ -163,57 +352,96 @@ thscreens *new_thscreens(
 	if (bi < 0)			/* Strange */
 		return NULL;
 
-	/* Create each screening object from one defined screen. */
-	/* Use the 0'th plane screen */
-	/* Stagger the screens with a round of 9 offset */
-	for (i = 0; i < t->np; i++) {
-		int xoff = ((i % 3) * screens[bi].width)/3;
-		int yoff = (((i/3) % 3) * screens[bi].height)/3;
-		void   *cx = NULL;
-		double (*lf)(void *cntx, double in) = 0;
-		if (cntx != NULL)
-			cx = cntx[i];
-		if (lutfunc != NULL)
-			lf = lutfunc[i];
-
-		DBG(("thscreens: creating plane %d/%d thscreen, offset %d %d\n",i,t->np,xoff,yoff));
-		if ((t->sc[i] = new_thscreen(screens[bi].width, screens[bi].height, xoff, yoff,
-		                          screens[bi].asp, swap, screens[bi].list[0],
-		                          ie, oebpc, oelev, oevalues, oo, overlap,
-		                          cx, lf)) == NULL) {
-			for (--i; i >= 0; i--)
-				t->sc[i]->del(t->sc[i]);
-			free(t->sc);
+	if (t->edif) {
+		int j;
+		int npix;
+
+		t->screen = screen_edscreens;
+
+		t->luts = imatrix(0, t->np-1, 0, 65535);
+
+		/* Create a suitable LUT from the given function */
+		/* Input is either 8 or 16 bits, output is always 16 bits */
+		for (j = 0; j < t->np; j++) {
+			for (i = 0; i < 65536; i++) {
+				if (lutfunc != NULL && lutfunc[j] != NULL) {
+					double v = i/65535.0;
+					v = lutfunc[j](cntx[j], v);
+					t->luts[j][i] = (int)(v * 65535.0 + 0.5);
+				} else
+					t->luts[j][i] = i;
+			}
+		}
+
+		if ((t->so = new_sobol(1)) == NULL) {
+			DBG(("thscreens: new_sobol() failed\n"));
+			return NULL;
+		}
+
+
+	} else {
+
+		t->screen = screen_thscreens;
+
+		if ((t->sc = malloc(sizeof(thscreen *) * t->np)) == NULL) {
 			free(t);
-			DBG(("thscreens: new_thscreen() failed\n"));
+			DBG(("thscreens: malloc of thscreens->sc[] failed\n"));
 			return NULL;
 		}
+
+		/* Create each screening object from one defined screen. */
+		/* Use the 0'th plane screen */
+		/* Stagger the screens with a round of 9 offset */
+		for (i = 0; i < t->np; i++) {
+			int xoff = ((i % 3) * screens[bi].width)/3;
+			int yoff = (((i/3) % 3) * screens[bi].height)/3;
+			void   *cx = NULL;
+			double (*lf)(void *cntx, double in) = 0;
+			if (cntx != NULL)
+				cx = cntx[i];
+			if (lutfunc != NULL)
+				lf = lutfunc[i];
+	
+			DBG(("thscreens: creating plane %d/%d thscreen, offset %d %d\n",i,t->np,xoff,yoff));
+			if ((t->sc[i] = new_thscreen(screens[bi].width, screens[bi].height, xoff, yoff,
+			                          screens[bi].asp, swap, screens[bi].list[0],
+			                          ie, oebpc, oelev, oevalues, oo, overlap,
+			                          cx, lf)) == NULL) {
+				for (--i; i >= 0; i--)
+					t->sc[i]->del(t->sc[i]);
+				free(t->sc);
+				free(t);
+				DBG(("thscreens: new_thscreen() failed\n"));
+				return NULL;
+			}
+		}
 	}
 	DBG(("thscreens: returning nonexact match\n"));
+
 	return t;
 }
 
 /* ----------------------------------------------------------- */
-/* The kernel screening routin */
+/* The kernel stocastic screening routine */
 
 void thscreen16_8(
 	struct _thscreen *t,	/* Screening object pointer */
 	int width, int height,	/* Width and height to screen in pixels */
 	int xoff, int yoff,		/* Offset into screening pattern (must be +ve) */
-	unsigned char *_in,		/* Input pixel buffer */
-	unsigned long ipinc,	/* Increment between input pixels */
-	unsigned long ipitch,	/* Increment between input lines */
 	unsigned char *out,		/* Output pixel buffer */
-	unsigned long opinc,	/* Increment between output pixels */
-	unsigned long opitch	/* Increment between output lines */
+	unsigned long opinc,	/* Increment between output pixels in components */
+	unsigned long opitch,	/* Increment between output lines in components */
+	unsigned char *_in,		/* Input pixel buffer */
+	unsigned long ipinc,	/* Increment between input pixels in components */
+	unsigned long ipitch	/* Increment between input lines in components */
 ) {
 	unsigned short *in = (unsigned short *)_in;	/* Pointer to input pixel sized values */
 	int *lut = t->lut;			/* Copy of 8 or 16 -> 16 bit lookup table */
+	unsigned short *ein = in + height * ipitch;	/* Vertical end pixel marker */
+	unsigned short *ein1;				/* Horizontal end pixel markers */
 	unsigned char **oth, **eth; /* Current lines start, origin and end in screening table. */
 	int thtsize;				/* Overall size of threshold table */
 	unsigned char **eeth;		/* Very end of threshold table */
-	unsigned short *ein = in + height * ipitch;	/* Vertical end pixel marker */
-	unsigned short *ein1;				/* Horizontal end pixel markers */
 
 	{
 		unsigned char **sth;				/* Start point of line intable */
diff --git a/render/thscreen.h b/render/thscreen.h
index d3ac9b1..9a504bc 100644
--- a/render/thscreen.h
+++ b/render/thscreen.h
@@ -5,20 +5,22 @@
 /*
  * render2d
  *
- * Threshold screen pixel processing object.
+ * Threshold or Error diffusion screen pixel processing object.
  * (Simplified from DPS code)
  *
  * Author:  Graeme W. Gill
  * Date:    11/7/2005
  * Version: 1.00
  *
- * Copyright 2005, 2012 Graeme W. Gill
+ * Copyright 2005, 2012, 2014 Graeme W. Gill
  * All rights reserved.
  * This material is licenced under the GNU AFFERO GENERAL PUBLIC LICENSE Version 3 :-
  * see the License.txt file for licencing details.
  *
  */
 
+#define THMXCH2D  8           /* Maximum color channels */
+
 /* Light Separation in screening flag */
 typedef enum {
 	scls_false = 0,		/* Don't do light ink separation during screening. */
@@ -45,15 +47,31 @@ struct _thscreens {
 	int np;						/* Number of planes */
 	struct _thscreen **sc;		/* List of screens */
 
-	/* Screen a single color plane */
+	int oebpc;					/* Output encoding bits per component, 1,2,4 or 8 */
+	int oelev;					/* Output encoding levels. Must be <= 2 ^ oebpc */
+	int oevalues[256];			/* Output encoding values for each level */
+
+	int edif;					/* nz if using error diffusion */
+	int **luts;					/* Lookup tables */
+	int mxwidth;				/* max width in pixels of raster to be screened */ 
+	int lastyoff;				/* Last y offset */
+	float **ebuf;				/* Error buffer for each plane */
+								/* ebuf[][-1] is used for next pixel error */
+	
+	void (*quant)(void *qcntx, double *out, double *in); /* optional quantization func. for edif */
+	void *qcntx;					/* Context for quant */
+
+	sobol *so;					/* Random number generator for error diffusion */
+
+	/* Screen pixel values */
 	void (* screen)(			/* Pointer to dither function */
 		struct _thscreens *t,	/* Screening object pointer */
 		int width, int height,	/* Width and height to screen in pixels */
 		int xoff, int yoff,		/* Offset into screening pattern */
-		unsigned char *in,		/* Input pixel buffer */
-		unsigned long ipitch,	/* Increment between input lines */
 		unsigned char *out,		/* Output pixel buffer */
-		unsigned long opitch);	/* Increment between output lines */
+		unsigned long opitch,	/* Increment between output lines in components */
+		unsigned char *in,		/* Input pixel buffer */
+		unsigned long ipitch);	/* Increment between input lines in components */
 
 	void (* del)(				/* Destructor */
 		struct _thscreens *t);	/* Screening objects pointer */
@@ -68,7 +86,7 @@ thscreens *new_thscreens(
 	int exact,				/* Return only exact matches */
 	int nplanes,			/* Number of planes to screen */
 	double asp,				/* Target aspect ratio (== dpiX/dpiY) */
-	int size,				/* Target size */
+	int size,				/* Target screen size */
 	sc_iencoding ie,		/* Input encoding - must be scie_16 */
 	int oebpc,				/* Output encoding bits per component - must be 8 */
 	int oelev,				/* Output encoding levels. Must be <= 2 ^ oebpc */
@@ -76,8 +94,12 @@ thscreens *new_thscreens(
 							/* Must be oelev entries. Default is 0 .. oelev-1 */
 	sc_oorder oo,			/* Output bit ordering */
 	double overlap,			/* Overlap between levels, 0 - 1.0 */
+	int mxwidth,			/* max width in pixels of raster to be screened */ 
 	void   **cntx,			/* List of contexts for lookup table callback */
-	double (**lutfunc)(void *cntx, double in)	/* List of callback function, NULL if none */
+	double (**lutfunc)(void *cntx, double in),	/* List of callback function, NULL if none */
+	int edif,				/* nz if using error diffusion */
+	void (*quant)(void *qcntx, double *out, double *in), /* optional quantization func. for edif */
+	void *qcntx
 );
 
 /* ---------------------------- */
@@ -121,12 +143,12 @@ struct _thscreen {
 		struct _thscreen *t,	/* Screening object pointer */
 		int width, int height,	/* Width and height to screen in pixels */
 		int xoff, int yoff,		/* Offset into screening pattern */
-		unsigned char *in,		/* Input pixel buffer */
-		unsigned long ipinc,	/* Increment between input pixels */
-		unsigned long ipitch,	/* Increment between input lines */
 		unsigned char *out,		/* Output pixel buffer */
-		unsigned long opinc,	/* Increment between output pixels */
-		unsigned long opitch);	/* Increment between output lines */
+		unsigned long opinc,	/* Increment between output pixels in components */
+		unsigned long opitch,	/* Increment between output lines in components */
+		unsigned char *in,		/* Input pixel buffer */
+		unsigned long ipinc,	/* Increment between input pixels in components */
+		unsigned long ipitch);	/* Increment between input lines in components */
 
 	void (* del)(				/* Destructor */
 		struct _thscreen *t);	/* Screening object pointer */
@@ -136,8 +158,8 @@ struct _thscreen {
 /* Create a new thscreen object */
 /* Return NULL on error */
 thscreen *new_thscreen(
-	int width,					/* width in pixels */
-	int height,					/* Height in pixels */
+	int width,					/* width in pixels of screen */
+	int height,					/* Height in pixels of screen */
 	int xoff, int yoff,			/* Pattern offsets into width & height (must be +ve) */
 	double asp,					/* Aspect ratio (== dpiX/dpiY) */
 	int swap,					/* Swap X & Y to invert aspect ratio */
diff --git a/render/timage.c b/render/timage.c
index c4404f7..2b610d1 100644
--- a/render/timage.c
+++ b/render/timage.c
@@ -32,39 +32,43 @@
 #include <stdlib.h>
 #include <string.h>
 #include <math.h>
+#include <fcntl.h>
 #include "copyright.h"
 #include "aconfig.h"
 #include "numsup.h"
 #include "render.h"
 
 #define DEF_DPI 200
-#define DITHER 0			/* Test 8 bit didthering */
+#define DITHER 0			/* 1 for test 8 bit dithering, 2 for test error diffusion */
+							/* 0x8001 for dithering FG only, 0x8002 for err. diff. FG only */
+#undef PNG_MEM				/* Test PNG save to memory */
 
 void
 usage(void) {
 	fprintf(stderr,"Create test images, default hex RGB surface and wedge, Version %s\n",ARGYLL_VERSION_STR);
 	fprintf(stderr,"Author: Graeme W. Gill, licensed under the AGPL Version 3\n");
-	fprintf(stderr,"usage: timage [-options] outfile.tif\n");
-//	fprintf(stderr," -v             Verbose\n");
-	fprintf(stderr," -t             Generate rectangular gamut boundary test chart\n");
-	fprintf(stderr," -p steps       Generate a colorspace step chart with L* steps^2\n");
-	fprintf(stderr," -r res         Resolution in DPI (default %d)\n",DEF_DPI);
-	fprintf(stderr," -s             Smooth blend\n");
-	fprintf(stderr," -x             16 bit output\n");
-	fprintf(stderr," -4             CMYK output\n");
-	fprintf(stderr," -g prop        Percentage towards grey (default 0%%)\n");
-//	fprintf(stderr," -D	            Debug primitives plot */
-	fprintf(stderr," outfile.tif    Profile to check against\n");
+	fprintf(stderr,"usage: timage [-options] outfile.[tif|png]\n");
+//	fprintf(stderr," -v                Verbose\n");
+	fprintf(stderr," -t                Generate rectangular gamut boundary test chart\n");
+	fprintf(stderr," -p steps          Generate a colorspace step chart with L* steps^2\n");
+	fprintf(stderr," -r res            Resolution in DPI (default %d)\n",DEF_DPI);
+	fprintf(stderr," -s                Smooth blend\n");
+	fprintf(stderr," -x                16 bit output\n");
+	fprintf(stderr," -4                CMYK output\n");
+	fprintf(stderr," -g prop           Percentage towards grey (default 0%%)\n");
+	fprintf(stderr," -P                Save as PNG file (deffault TIFF)\n");
+//	fprintf(stderr," -D	               Debug primitives plot */
+	fprintf(stderr," outfile.[tif|png] Output TIFF or PNG file\n");
 	exit(1);
 	}
 
-int main(int argc, char *argv[])
-{
+int main(int argc, char *argv[]) {
 	int fa,nfa;				/* current argument we're looking at */
 	int verb = 0;
 	int rchart = 0;			/* Rectangular chart */
 	int schart = 0;			/* Step chart with steps^2 */
 	int smooth = 0;			/* Use smooth blending */
+	rend_format fmt = tiff_file;	/* Output filr format */
 	int debugchart = 0;		/* Debug chart */
 	double res = DEF_DPI;
 	depth2d depth = bpc8_2d;
@@ -104,20 +108,20 @@ int main(int argc, char *argv[])
 			if (argv[fa][1] == '?')
 				usage();
 
-			else if (argv[fa][1] == 'v' || argv[fa][1] == 'V')
+			else if (argv[fa][1] == 'v')
 				verb = 1;
 
 			/* Rectangular chart */
-			else if (argv[fa][1] == 't' || argv[fa][1] == 'T') {
+			else if (argv[fa][1] == 't') {
 				rchart = 1;
 				schart = 0;
 
 			/* Smooth blending */
-			} else if (argv[fa][1] == 's' || argv[fa][1] == 'S')
+			} else if (argv[fa][1] == 's')
 				smooth = 1;
 
 			/* 16 bit depth */
-			else if (argv[fa][1] == 'x' || argv[fa][1] == 'X')
+			else if (argv[fa][1] == 'x')
 				depth = bpc16_2d;
 
 			/* cmyk */
@@ -125,7 +129,7 @@ int main(int argc, char *argv[])
 				cmyk = 1;
 
 			/* step chart */
-			else if (argv[fa][1] == 'p' || argv[fa][1] == 'P') {
+			else if (argv[fa][1] == 'p') {
 				fa = nfa;
 				if (na == NULL) usage();
 				schart = atoi(na);
@@ -133,13 +137,18 @@ int main(int argc, char *argv[])
 				rchart = 0;
 			}
 
+			/* PNG file */
+			else if (argv[fa][1] == 'P') {
+				fmt = png_file;
+			}
+
 			/* debug chart */
 			else if (argv[fa][1] == 'D') {
 				debugchart = 1;
 			}
 
 			/* resolution */
-			else if (argv[fa][1] == 'r' || argv[fa][1] == 'R') {
+			else if (argv[fa][1] == 'r') {
 				fa = nfa;
 				if (na == NULL) usage();
 				res = atof(na);
@@ -147,7 +156,7 @@ int main(int argc, char *argv[])
 			}
 
 			/* grey blend */
-			else if (argv[fa][1] == 'g' || argv[fa][1] == 'G') {
+			else if (argv[fa][1] == 'g') {
 				fa = nfa;
 				if (na == NULL) usage();
 				gbf = 1.0 - 0.01 * atof(na);
@@ -175,7 +184,7 @@ int main(int argc, char *argv[])
 		if (cmyk)
 			error("CMYK not supported for test chart");
 
-		if ((r = new_render2d(w, h, NULL, res, res, rgb_2d, 0, depth, DITHER)) == NULL) {
+		if ((r = new_render2d(w, h, NULL, res, res, rgb_2d, 0, depth, DITHER, NULL, NULL)) == NULL) {
 			error("new_render2d() failed");
 		}
 	
@@ -300,7 +309,7 @@ int main(int argc, char *argv[])
 		h = (1.0 + 2.0 * bb) * hh;
 		w = (4.0 * bb + 0.25 + 2.0 * r3o2) * hh;
 	
-		if ((r = new_render2d(w, h, NULL, res, res, cmyk ? cmyk_2d : rgb_2d, 0, depth, DITHER)) == NULL) {
+		if ((r = new_render2d(w, h, NULL, res, res, cmyk ? cmyk_2d : rgb_2d, 0, depth, DITHER, NULL, NULL)) == NULL) {
 			error("new_render2d() failed");
 		}
 	
@@ -568,7 +577,7 @@ int main(int argc, char *argv[])
 		h = (1.0 + 2.0 * bb) * hh;
 		w = (2.0 * bb + 0.20 * 7.0) * hh;
 	
-		if ((r = new_render2d(w, h, NULL, res, res, rgb_2d, 0, depth, DITHER)) == NULL) {
+		if ((r = new_render2d(w, h, NULL, res, res, rgb_2d, 0, depth, DITHER, NULL, NULL)) == NULL) {
 			error("new_render2d() failed");
 		}
 	
@@ -638,7 +647,7 @@ int main(int argc, char *argv[])
 		bs = (bb * hh)/(schart + 1.0);
 		ss = hh * (1.0 - bb)/schart;
 	
-		if ((r = new_render2d(w, h, NULL, res, res, lab_2d, 0, depth, DITHER)) == NULL) {
+		if ((r = new_render2d(w, h, NULL, res, res, lab_2d, 0, depth, DITHER, NULL, NULL)) == NULL) {
 			error("new_render2d() failed");
 		}
 	
@@ -673,7 +682,36 @@ int main(int argc, char *argv[])
 		}
 	}
 
-	r->write(r, outname,1);
+#ifdef PNG_MEM
+	{
+		char *nmode = "w";
+		FILE *fp;
+		unsigned char *buf;
+		size_t len, wlen;
+
+#if !defined(O_CREAT) && !defined(_O_CREAT)
+# error "Need to #include fcntl.h!"
+#endif
+#if defined(O_BINARY) || defined(_O_BINARY)
+		nmode = "wb";
+#endif
+
+		if (r->write(r, "MemoryBuf", 1, &buf, &len, png_mem))
+			error("render->write failed");
+
+		if ((fp = fopen(outname, nmode)) == NULL)
+			error("render2d: open '%s' for writing",outname);
+
+		if (len != (wlen = fwrite(buf, 1, len, fp)))
+			error("render2d: writing %u bytes to '%s' failed (wrote %u)",len,outname,wlen);
+
+		if (fclose(fp))
+			error("render2d: failed to close after writing",outname);
+	}
+#else
+	if (r->write(r, outname, 1, NULL, NULL, fmt))
+		error("render->write failed");
+#endif
 	r->del(r);
 
 	return 0;