Mercurial > hg > octave-image

changeset 158:b1ba33b88de2

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Changes to handle MxNx3 images of different numeric classes: * imread.m: Return the appropriate numeric class. Colour images are of size MxNx3, gray images MxN. Uses ImageMagick's "identify" to find the image type (b&w/grayscale/colour). * jpgread.cc: Error handler modified to prevent segfaults. Colour images are returned as MxNx3 matrices, gray images as MxN. Conversion from colour to grey-level is not done anymore (rather use im2gray). * pngread.cc: No segfault on loading invalid file types. Handle different depths correctly. Split into pngread, pngwrite and pngcanvas. Return the appropriate numeric classes and matrix sizes. * imread.m: Prefer jpgread and pngread to ImageMagick.
author sjvdw
date Sat, 10 Sep 2005 16:07:47 +0000
parents 4aaa6584062c
children 3b34641f1dc9
files Makefile imread.m jpgread.cc pngcanvas.h pngread.cc pngwrite.cc
diffstat 6 files changed, 638 insertions(+), 735 deletions(-) [+]
line wrap: on
line diff
--- a/Makefile
+++ b/Makefile
@@ -22,6 +22,6 @@
 	$(MKOCTFILE) $< -lpng
 
 pngwrite.oct:
-	$(LN_S) pngread.oct $@
+	${MKOCTFILE) $< -lpng
 
 clean: ; -$(RM) *.o octave-core core *.oct *~
--- a/imread.m
+++ b/imread.m
@@ -5,392 +5,141 @@
 ## the Free Software Foundation; either version 2, or (at your option)
 ## any later version. USE THIS SOFTWARE AT YOUR OWN RISK.
 
-#IMREAD: read images into octave from various file formats
-#
-# Note: this requires the ImageMagick "convert" utility
-#       get this from www.imagemagick.org if required
-#       additional documentation of options is available from the
-#       convert man page
-#
-# BASIC USAGE:
-# img    = imread( fname )
-#                 - img is a greyscale (0-255) of image in fname
-# [im,map]=imread( fname )
-#                 - map is a matrix of [r,g,b], 0-1 triples
-#                 - img is a matrix on indeces into map
-# [r,g,b]= imread( fname )
-#                 - r,g,b are red,green,blue (0-255) compondents
-#
-# Formats for image fname
-#   1. simple guess from extention ie "fig.jpg", "blah.gif"
-#   2. specify explicitly             "jpg:fig.jpg", "gif:blah.gif"
-#   3. specify subimage for multi-image format "tiff:file.tif[3]"
-#   4. raw images (row major format) specify geometry
-#                                      "raw:img[256x180]"
-#
-# IMREAD OPTIONS:
-# imread will support most of the options for convert.1
-#
-# img    = imread( fname , options)
-# [r,g,b]= imread( fname , options)
-#
-# where options is a string matrix (or list) of options
-#
-# example:   options= ["-rotate 25";
-#                      "-crop 200x200+150+150";
-#                      "-sample 200%" ];
-#   will rotate, crop, and then expand the image.
-#   note that the order of operations is important
-#
-# The following options are supported
-#  -antialias           remove pixel-aliasing
-#  -blur geometry       blur the image
-#  -border geometry     surround image with a border of color
-#  -bordercolor color   border color
-#  -colors value        maximum number of colors in the image
-#  -colorspace type     alternate image colorspace
-#  -contrast            enhance or reduce the image contrast
-#  -crop geometry       preferred size and location of the cropped image
-#  -despeckle           reduce the speckles within an image
-#  -dither              apply Floyd/Steinberg error diffusion to image
-#  -draw string         annotate the image with a graphic primitive
-#  -edge radius         apply a filter to detect edges in the image
-#  -emboss radius       emboss an image
-#  -enhance             apply a digital filter to enhance a noisy image
-#  -equalize            perform histogram equalization to an image
-#  -filter type         use this filter when resizing an image
-#  -flip                flip image in the vertical direction
-#  -flop                flop image in the horizontal direction
-#  -font name           font for rendering text
-#  -frame geometry      surround image with an ornamental border
-#  -fuzz distance       colors within this distance are considered equal
-#  -gamma value         level of gamma correction
-#  -geometry geometry   perferred size or location of the image
-#  -gaussian geometry   gaussian blur an image
-#  -gravity type        vertical and horizontal text placement
-#  -implode amount      implode image pixels about the center
-#  -intent type         Absolute, Perceptual, Relative, or Saturation
-#  -interlace type      None, Line, Plane, or Partition
-#  -level value         adjust the level of image contrast
-#  -map filename        transform image colors to match this set of colors
-#  -median radius       apply a median filter to the image
-#  -modulate value      vary the brightness, saturation, and hue
-#  -monochrome          transform image to black and white
-#  -mosaic              create an mosaic from an image sequence
-#  -negate              replace every pixel with its complementary color
-#  -noise radius        add or reduce noise in an image
-#  -normalize           transform image to span the full range of colors
-#  -paint radius        simulate an oil painting
-#  -raise value         lighten/darken image edges to create a 3-D effect
-#  -region geometry     apply options to a portion of the image
-#  -roll geometry       roll an image vertically or horizontally
-#  -rotate degrees      apply Paeth rotation to the image
-#  -sample geometry     scale image with pixel sampling
-#  -scale geometry      resize image
-#  -scene value         image scene number
-#  -segment values      segment an image
-#  -seed value          pseudo-random number generator seed value
-#  -shade degrees       shade the image using a distant light source
-#  -sharpen geometry    sharpen the image
-#  -shave geometry      shave pixels from the image edges
-#  -shear geometry      slide one edge of the image along the X or Y axis
-#  -size geometry       width and height of image
-#  -solarize threshold  negate all pixels above the threshold level
-#  -spread amount       displace image pixels by a random amount
-#  -swirl degrees       swirl image pixels about the center
-#  -texture filename    name of texture to tile onto the image background
-#  -threshold value     threshold the image
-#  -transparent color   make this color transparent within the image
-#  -treedepth value     depth of the color tree
-#  -type type           image type
-#  -unsharp geometry    sharpen the image
-#  -wave geometry       alter an image along a sine wave
+## usage: I = imread(fname)
+##
+## Read images from various file formats.
+##
+##   The size and numeric class of the output depends on the
+##   format of the image.  A colour image is returned as an
+##   MxNx3 matrix.  Grey-level and black-and-white images are
+##   of size MxN.
+##     The colour depth of the image determines the numeric
+##   class of the output: 'uint8' or 'uint16' for grey
+##   and colour, and 'logical' for black and white.
+##
+##   Note: For image formats other than jpeg and png, the
+##         ImageMagick "convert" and "identify" utilities
+##         are needed. ImageMagick can be found at
+##
+##         www.imagemagick.org
+##
+
+## Author: Andy Adler
+##
+## Modified: Stefan van der Walt <stefan@sun.ac.za>
+## Date: 24 January 2005
+
+function varargout = imread(filename, options)
+    if (nargin != 1)
+	usage("I = imread(filename)")
+    endif
+
+    if !isstr(filename)
+	error("imread: filename must be a string")
+    endif
 
-function [out1,out2,out3]= imread(filename, options );
-
-try save_empty_list_elements_ok= empty_list_elements_ok;
-catch save_empty_list_elements_ok= 0;
-end
-try save_warn_empty_list_elements= warn_empty_list_elements;
-catch save_warn_empty_list_elements= 0;
-end
-unwind_protect
-empty_list_elements_ok= 1;
-warn_empty_list_elements= 0;
-
-# some older versions of octave didn't seem handle piped output correctly
-usepipe=1;
+    fn = file_in_path(IMAGEPATH, filename);
+    if isempty(fn)
+	error("imread: cannot find %s", filename);
+    endif
 
-if (nargin == 0)
-  usage (["img    =  imread (filename,options)\n", ...
-          "[r,g,b]=  imread (filename,options)\n", ...
-          "[img,map]=imread (filename,options)" ]);
-elseif (! ischar (filename))
-  error ("imread: expecting filename as a string");
-endif
-
-if any(filename==':') || any(filename=='[')
-   # we've using special imagemagick characters, so we don't
-   # do any octave path processing
-   fname= filename;
-else
-   fname= file_in_path(IMAGEPATH,filename);
-   if isempty(fname)
-      error(['imread: file ' filename ' not found']);
-   end
-end
+    ## divert jpegs and pngs to "jpgread" and "pngread"
+    [ig, ig, ext] = fileparts(fn);
+    ext = upper(ext);
+    if ( file_in_loadpath("jpgread.oct") &&
+	(strcmp(ext, ".JPG") || strcmp(ext, ".JPEG")) )
+	varargout{1} = jpgread(fn);
+	break
+    endif
+    if ( file_in_loadpath("pngread.oct") && (strcmp(ext, ".PNG")) )
+	varargout{1} = pngread(fn);
+	break
+    endif	
 
-# Put together the options string
-# TODO: add some error checking here
-if nargin==1
-   option_str="";
-else
-   option_str="";
-   if     ( is_list( options ) )
-      for i= 1:length(options)
-         option_str=[option_str," ", nth(options,i) ];
-      end
-   elseif ( ischar( options ) )
-      for i= 1:size(options,1)
-         option_str=[option_str," ", options(i,:) ];
-      end
-   else
-     error ("imread: expecting options as a string");
-   end
-end
-#
-# decode the nargout to output what the user wants
-#
-if     nargout==3;
-    wantgreyscale= 0;
-    wantmap= 0;
-    outputtype="ppm";
-elseif nargout==2;
-    wantgreyscale= 0; wantmap= 1;
-    outputtype="ppm";
-else
-    wantgreyscale= 1; wantmap= 0;
-    outputtype="pgm";
-end
+    [ident, sys] = system(sprintf("identify -verbose %s | grep -e Depth -e Type",
+				  fn));
+    if (sys != 0)
+	error("imread: error running ImageMagick's 'identify' on %s", fn)
+    endif
+    depth = re_grab("Depth: ([[:digit:]]{1,2})", ident);
+    imtype = re_grab("Type: ([[:alpha:]]*)", ident);
 
-if usepipe
-   pname= sprintf("convert %s '%s' %s:- 2>/dev/null",
-                  option_str, fname, outputtype);
-   fid= popen(pname ,'r');
-else
-   tnam= tmpnam();
-   cmd= sprintf("convert %s '%s' '%s:%s' 2>/dev/null ",
-                  option_str, fname, outputtype, tnam);
-   system(cmd);
-   fid= fopen(tnam,"rb");
-end
+    depth = str2num(depth);
+    if isempty(depth) || (pow2(nextpow2(depth)) != depth)
+	error("imread: invalid image depth %s", depth)
+    endif
 
-
-#
-# can we open the pipe?
-# if not 1) The file format is wrong and the conver program has bailed out
-#        2) The apropriate converter program hasn't been installed
-#
-   if fid<0;
-      if ~usepipe
-         unlink(tnam);
-      end
-      error('could not read image data. Is ImageMagick installed?');
-   end
+    if !(strcmp(imtype, "Bilevel") || strcmp(imtype, "Grayscale") ||
+	 strcmp(imtype, "TrueColor"))
+	error("imread: unknown image type '%s'", imtype);
+    endif
 
-# get file type
-   line= fgetl( fid );
-   if     strcmp(line, 'P1');   bpp=1; spp=1; bindata=0;
-   elseif strcmp(line, 'P4');   bpp=1; spp=1; bindata=1;
-   elseif strcmp(line, 'P2');   bpp=8; spp=1; bindata=0;
-   elseif strcmp(line, 'P5');   bpp=8; spp=1; bindata=1;
-   elseif strcmp(line, 'P3');   bpp=8; spp=3; bindata=0;
-   elseif strcmp(line, 'P6');   bpp=8; spp=3; bindata=1;
-   else
-      if usepipe
-         pclose(fid);
-      else
-         fclose(fid); unlink(tnam);
-      end
-      error(['Image format error for ',fname,':line=', char(line)]);
-   end
+    switch (imtype)
+	case("Bilevel")
+	    fmt = "pgm";
+	case("Grayscale")
+	    fmt = "pgm";
+	case("TrueColor")
+	    fmt = "ppm";
+    endswitch
+    
+    ## Why are pipes so slow?
+    ##    cmd = sprintf("convert -flatten -strip %s %s:-", fn, fmt);
 
-# ignore comments
-   line= fgetl( fid );
-   while length(line)==0 || line(1)=='#'
-      line= fgetl( fid );
-   end
-
-# get width, height
-   [wid, hig]= sscanf( line, '%d %d', 'C' );
-
-# get max component value
-   if bpp==8
-     max_component= sscanf( fgetl(fid), '%d' );
-     if (max_component > 255)
-       # The PGM standard supports only values below 65536
-       bpp = 16;
-     endif
-   end
-
-   if bindata
-      if(bpp == 16)
-         # PGM format has MSB first, i.e. big endian
-         data = fread(fid, "uint16", 0, "ieee-be");
-      else
-         data = fread(fid);
-      endif
-      numdata= size(data,1);
+    tmpf = tmpnam();
+    cmd = sprintf("convert -flatten -strip '%s' '%s:%s' 2>/dev/null",
+		  fn, fmt, tmpf);
 
-      if bpp==1
-         data= rem( floor( (data*ones(1,8)) ./ ...
-                 (ones(length(data),1)*[128 64 32 16 8 4 2 1]) ), 2)';
-      end
-   else
-      numdata= wid*hig*spp;
-      data= zeros( numdata,1 );
-      dptr= 1;
-
-         line= fgetl( fid );
-      while !feof( fid)
-         rdata= sscanf( line ,' %d');
-         nptr= dptr + size(rdata,1);
-         data( dptr:nptr-1 ) = rdata;
-         dptr= nptr;
-         line= fgetl( fid );
-      end # feof
-   end #if bindata
+    [ignored, sys] = system(cmd);
+    
+    if (sys != 0)
+	error("imread: error running ImageMagick's 'convert'");
+	unlink(tmpf);
+    endif
+    
+    try
+	##fid = popen(cmd, "r");
+	fid = fopen(tmpf, "rb");
+    catch
+	unlink(tmpf);
+	error("imread: could not open temporary file %s", tmpf)
+    end_try_catch
 
-   if usepipe
-      pclose(fid);
-   else
-      fclose(fid); unlink(tnam);
-   end
-
-   if spp==1
-      greyimg= reshape( data(:), wid, hig )';
-   elseif spp==3
-      redimg= reshape( data(1:spp:numdata), wid, hig )';
-      grnimg= reshape( data(2:spp:numdata), wid, hig )';
-      bluimg= reshape( data(3:spp:numdata), wid, hig )';
-   else
-      error(sprintf("imread: don't know how to handle pnm with spp=%d",spp));
-   end
-
-
-#   This section outputs the image in the desired output format
-# if the user requested the colour map, the we regenerate it from
-# the image.
-#
-#   Of course, 1) This may result in huge colourmaps
-#              2) The colourmap will not be in the same order as
-#                   in the original file
-
-if wantgreyscale
+    fgetl(fid); # P5 or P6 (pgm or ppm)
+    [width, height] = sscanf(fgetl(fid), "%d %d", "C");
+    fgetl(fid); # ignore max components
 
-   if exist('greyimg')
-      out1= greyimg;
-   elseif exist('idximg')
-      greymap= mean(map')';
-      out1= reshape( greymap( idximg ) , size(idximg,1), size(idximg,2) );
-   else
-      out1= ( redimg+grnimg+bluimg ) / 3 ;
-   end
-
-elseif wantmap
-
-   if exist('idximg')
-      out1= idximg;
-      out2= map;
-   elseif exist('greyimg')
-      [simg, sidx] = sort( greyimg(:) );
-      [jnkval, sidx] = sort( sidx );
-
-      dimg= [1; diff(simg)>0 ];
-      cimg= cumsum( dimg );
-      out1= reshape( cimg( sidx ) , hig, wid);
-      out2= ( simg(find( dimg ))*[1,1,1] - 1)/256;
-   else
-#
-# Generate a colourmap for RGB images by packing RGB into
-# an integer and assigning the unique values to a colormap.
-#
-      c_range= 256; # use 65536 for 48 bit color
-      [simg, sidx] = sort( round(redimg(:)) + ...
-                  c_range *round(grnimg(:)) + ...
-                c_range^2 *round(bluimg(:)) );
-      [jnkval, sidx] = sort( sidx );
-
-      dimg= [1; diff(simg)>0 ];
-      cimg= cumsum( dimg );
-      out1= reshape( cimg( sidx ) , hig, wid );
-      tmpv= simg(find( dimg ));
-      out2= [ rem(tmpv,c_range), ...
-              rem(floor(tmpv/c_range), c_range), ...
-              floor(tmpv/c_range^2) ]/c_range;
-   end
-
-
-else
+    if(depth == 16)
+	## PGM format has MSB first, i.e. big endian
+	data = fread(fid, "uint16", 0, "ieee-be");
+    else
+        data = fread(fid);
+    endif
+    
+    fclose(fid);
+    unlink(tmpnam);
 
-   if exist('greyimg')
-      out1= greyimg;
-      out2= greyimg;
-      out3= greyimg;
-   else
-      out1= redimg;
-      out2= grnimg;
-      out3= bluimg;
-   end
-
-end
-
-unwind_protect_cleanup
-empty_list_elements_ok= save_empty_list_elements_ok;
-warn_empty_list_elements= save_warn_empty_list_elements;
-end_unwind_protect
+    varargout = {};
+    switch (imtype)
+ 	case ("Bilevel")
+ 	    varargout{1} = logical(reshape(data, height, width)');
+ 	case ("Grayscale") 
+ 	    varargout{1} = uint8(reshape(data, height, width)');
+ 	case ("TrueColor")
+ 	    varargout{1} = cat(3, reshape(data(1:3:end), height, width)',
+ 			       reshape(data(2:3:end), height, width)',
+ 			       reshape(data(3:3:end), height, width)');
+ 	    eval(sprintf("varargout{1} = uint%d(varargout{1});", depth));
+		
+    endswitch
+endfunction
 
-#
-# $Log$
-# Revision 1.14  2005/09/08 03:12:29  pkienzle
-# [for Bill Denney] setstr -> char
-#
-# Revision 1.13  2005/09/08 02:00:17  pkienzle
-# [for Bill Denney] isstr -> ischar
-#
-# Revision 1.12  2005/07/21 16:03:01  aadler
-# Improve error messages and use pipes
-#
-# Revision 1.11  2005/05/25 03:43:40  pkienzle
-# Author/Copyright consistency
-#
-# Revision 1.10  2005/04/25 01:05:28  aadler
-# added GPL copyrights
-#
-# Revision 1.9  2004/02/18 14:54:10  pkienzle
-# Colour values from ppm are 0-255 not 1-256
-#
-# Revision 1.8  2003/11/14 16:55:10  tpikonen
-# Fix endianess bug in 16-bit reads.
-#
-# Revision 1.6  2003/09/12 14:22:42  adb014
-# Changes to allow use with latest CVS of octave (do_fortran_indexing, etc)
-#
-# Revision 1.5  2003/07/25 19:11:41  pkienzle
-# Make sure all files names referenced in system calls are wrapped in quotes
-# to protect against spaces in the path.
-#
-# Revision 1.4  2002/11/27 08:40:11  pkienzle
-# author/license updates
-#
-# Revision 1.3  2002/03/19 18:14:13  aadler
-# unfortunately using popen seems to create problems, mostly
-# on win32, but also on linux, so we need to move to a tmpfile approach
-#
-# Revision 1.2  2002/03/17 05:26:14  aadler
-# now accept filenames with spaces
-#
-# Revision 1.1  2002/03/11 01:56:47  aadler
-# general image read/write functionality using imagemagick utilities
-#
-#
+function value = re_grab(re, str)
+    idx = regexp(re, str);
+    if !isempty(idx)
+	idx = idx(2,:);
+	value = substr(str, idx(1), diff(idx)+1);	
+    else
+	value = "";
+    endif    
+endfunction
--- a/jpgread.cc
+++ b/jpgread.cc
@@ -9,32 +9,39 @@
  // $Id$                                     
  //
 
-#include <octave/oct.h>
-#include <iostream>
+/* Modified: Stefan van der Walt <stefan@sun.ac.za>
+ * Date: 27 January 2004
+ * - Manual error handler to prevent segfaults in Octave.
+ * - Use uint8NDArray for output.
+ */
 
 /*
- * Simple jpeg reading MEX-file. 
- *
  * Compilation:
  * First, try
  *   mkoctfile jpgread.cc -ljpeg
  *
- * If this doesn't work, then do
- *
- * Calls the jpeg library which is part of 
+ * If this doesn't work, install the jpeg library which is part of
  * "The Independent JPEG Group's JPEG software" collection.
  *
- * The jpeg library came from,
+ * The jpeg library came from
  *
  * ftp://ftp.uu.net/graphics/jpeg/jpegsrc.v6.tar.gz
+ *
+ * Extract and build the library:
  * tar xvfz jpegsrc.v6.tar.gz
  * cd jpeg-6b
  * ./configure
  * make
  * make test
- * mkoctfile jpgread.cc -I<jpeg-6b include> -L<jpeg-6b lib> -ljpeg
+ *
+ * Compile this file using:
+ * mkoctfile jpgread.cc -I<jpeg-6b include dir> -L<jpeg-6b lib dir> -ljpeg
  */
 
+#include <octave/oct.h>
+#include <iostream>
+#include <csetjmp>
+
 #ifdef __cplusplus
 extern "C" {
 #endif
@@ -45,137 +52,130 @@
 } //extern "C"
 #endif
 
-DEFUN_DLD (jpgread, args, nargout ,
-"JPGREAD Read a JPEG file from disk. \n\
-       [R,G,B] = jpgread('filename') reads the specified file \n\
-       and returns the Red, Green, and Blue intensity matrices.\n\
-\n\
-       [X] = jpgread('filename') reads the specified file \n\
-       and returns the average intensity matrix.")
-{ 
-   octave_value_list retval;
-   int nargin  = args.length();
+struct oct_error_mgr {
+    struct jpeg_error_mgr pub;    /* "public" fields */
+    jmp_buf setjmp_buffer;        /* for return to caller */
+};
 
-   FILE * infile;
+typedef struct oct_error_mgr * oct_error_ptr;
 
-   JSAMPARRAY buffer;
-   long row_stride;
-   struct jpeg_decompress_struct cinfo;
-   struct jpeg_error_mgr jerr;
-
-//
-// We bail out if the input parameters are bad
-//
-   if (nargin != 1 || !args(0).is_string() ) {
-      print_usage ("jpgread");
-      return retval;
-   }
-
-   if (nargout != 1 && nargout != 3 ) {
-      print_usage ("jpgread");
-      return retval;
-   }
+METHODDEF(void)
+oct_error_exit (j_common_ptr cinfo)
+{
+    /* cinfo->err really points to an oct_error_mgr struct, so coerce pointer */
+    oct_error_ptr octerr = (oct_error_ptr) cinfo->err;
+    
+    /* Format error message and send to interpreter */
+    char errmsg[JMSG_LENGTH_MAX];
+    (octerr->pub.format_message)(cinfo, errmsg);
+    error("jpgread: %s", errmsg);
+    
+    /* Return control to the setjmp point */
+    longjmp(octerr->setjmp_buffer, 1);
+}
 
-//
-// Open jpg file
-//
-   std::string filename = args(0).string_value();
-   if ((infile = fopen(filename.c_str(), "rb")) == NULL) {
-      error("Couldn't open file");
-      return retval;
-   }
-
-//
-// Initialize the jpeg library
-//
-   cinfo.err = jpeg_std_error(&jerr);
-   jpeg_create_decompress(&cinfo);
-
-//
-// Read the jpg header to get info about size and color depth
-//
-   jpeg_stdio_src(&cinfo, infile);
-   jpeg_read_header(&cinfo, TRUE);
-   jpeg_start_decompress(&cinfo);
-   /*
-   if (cinfo.output_components == 1) { // Grayscale 
-      jpeg_destroy_decompress(&cinfo);
-      error("Grayscale jpegs not supported");
-   }
-   */
-  
-//
-// Allocate buffer for one scan line
-//
-
-   row_stride = cinfo.output_width * cinfo.output_components;
-   buffer = (*cinfo.mem->alloc_sarray)
-                 ((j_common_ptr) &cinfo, JPOOL_IMAGE, row_stride, 1);
+DEFUN_DLD (jpgread, args, nargout ,
+"usage: I = jpgread('filename')\n\
+\n\
+  Read a JPEG file from disk.\n\
+\n\
+  For a grey-level image, the output is an MxN matrix. For a\n\
+  colour image, three such matrices are returned (MxNx3),\n\
+  representing the red, green and blue components. The output\n\
+  is of class 'uint8'.\n\
+\n\
+  See also: imread, im2double, im2gray, im2rgb.")
+{ 
+    octave_value_list retval;
+    int nargin  = args.length();
+    
+    FILE * infile;
+    
+    JSAMPARRAY buffer;
+    long row_stride;
+    struct jpeg_decompress_struct cinfo;
+    struct oct_error_mgr jerr;
 
-//
-// Create 3 matrices, One each for the Red, Green, and Blue component of the image.
-// or create 1 matrix for the avg intensity
-//
-// Now, loop through each of the scanlines. For each, copy the image
-// data from the buffer, convert to double.
-//
-
-   if (nargout == 3) {
-      Matrix red   ( cinfo.output_height , cinfo.output_width );
-      Matrix green ( cinfo.output_height , cinfo.output_width );
-      Matrix blue  ( cinfo.output_height , cinfo.output_width );
-
-      for (unsigned long j=0; cinfo.output_scanline < cinfo.output_height; j++) {
-         jpeg_read_scanlines(&cinfo, buffer,1);
-   
-         if (cinfo.output_components == 1) { // Grayscale 
-            for (unsigned long i=0; i<cinfo.output_width; i++) {
-               red(j,i)   =
-               green(j,i) =
-               blue(j,i)  = (double) buffer[0][i*3+2];
-            } // for i
-	 } else { // not Grayscale
-            for (unsigned long i=0; i<cinfo.output_width; i++) {
-               red(j,i)   = (double) buffer[0][i*3+0];
-               green(j,i) = (double) buffer[0][i*3+1];
-               blue(j,i)  = (double) buffer[0][i*3+2];
-            } // for i
-	 } // not grayscale
-      } //for j
-
-      retval(0)= red;
-      retval(1)= green;
-      retval(2)= blue;
-
-   } // if nargout==3
-   else {
-      Matrix avg   ( cinfo.output_height , cinfo.output_width );
-
-      for (unsigned long j=0; cinfo.output_scanline < cinfo.output_height; j++) {
-         jpeg_read_scanlines(&cinfo, buffer,1);
-   
-         if (cinfo.output_components == 1) { // Grayscale 
-            for (unsigned long i=0; i<cinfo.output_width; i++) {
-               avg(j,i)   = (double) buffer[0][i];
-            } // for i
-	 } else { // not Grayscale
-            for (unsigned long i=0; i<cinfo.output_width; i++) {
-               avg(j,i)   = ((double) buffer[0][i*3+0] +
-                                         buffer[0][i*3+1] +
-                                         buffer[0][i*3+2] ) / 3;
-            } // for i
-	} // not Grayscale
-      } //for j
-
-      retval(0)= avg;
-   } // if nargout
-
-//
-// Clean up
-//
-   jpeg_finish_decompress(&cinfo);
-   jpeg_destroy_decompress(&cinfo);
-   fclose(infile);
-   
-   return retval;
-};
+    //
+    // We bail out if the input parameters are bad
+    //
+    if ((nargin != 1) || !args(0).is_string() || (nargout != 1)) {
+	print_usage ("jpgread");
+	return retval;
+    }    
+    
+    //
+    // Open jpg file
+    //
+    std::string filename = args(0).string_value();
+    if ((infile = fopen(filename.c_str(), "rb")) == NULL) {
+	error("jpgread: couldn't open file %s", filename.c_str());
+	return retval;
+    }
+    
+    //
+    // Initialize the jpeg library
+    //
+    cinfo.err = jpeg_std_error(&jerr.pub);
+    jerr.pub.error_exit = oct_error_exit;
+    if (setjmp(jerr.setjmp_buffer)) {
+	/* If we get here, the JPEG code has signaled an error.
+	 * We need to clean up the JPEG object, close the input file, and return.
+	 */
+	jpeg_destroy_decompress(&cinfo);
+	fclose(infile);
+	return retval;
+    }
+    
+    jpeg_create_decompress(&cinfo);
+    
+    //
+    // Read the jpg header to get info about size and color depth
+    //
+    jpeg_stdio_src(&cinfo, infile);
+    jpeg_read_header(&cinfo, TRUE);
+    jpeg_start_decompress(&cinfo);
+    
+    //
+    // Allocate buffer for one scan line
+    //
+    row_stride = cinfo.output_width * cinfo.output_components;
+    buffer = (*cinfo.mem->alloc_sarray)
+	((j_common_ptr) &cinfo, JPOOL_IMAGE, row_stride, 1);
+    
+    //
+    // Create an NDArray for the output.  Loop through each of the
+    // scanlines and copy the image data from the buffer.
+    //
+    
+    dim_vector dim = dim_vector();
+    dim.resize(3);
+    dim(0) = cinfo.output_height;
+    dim(1) = cinfo.output_width;
+    dim(2) = cinfo.output_components;
+    uint8NDArray out = uint8NDArray(dim, 0);
+    
+    Array<int> coord = Array<int> (3);
+    for (unsigned long j=0; cinfo.output_scanline < cinfo.output_height; j++) {
+	jpeg_read_scanlines(&cinfo, buffer, 1);
+	
+	coord(0) = j;
+	for (unsigned long i=0; i<cinfo.output_width; i++) {
+	    coord(1) = i;
+	    for (unsigned int c = 0; c < cinfo.output_components; c++) {
+		coord(2) = c;
+		out(coord) = buffer[0][i*cinfo.output_components+c];
+	    }
+	}
+    }
+    retval.append(out.squeeze());
+    
+    //
+    // Clean up
+    //
+    jpeg_finish_decompress(&cinfo);
+    jpeg_destroy_decompress(&cinfo);
+    fclose(infile);
+    
+    return retval;
+}
new file mode 100644
--- /dev/null
+++ b/pngcanvas.h
@@ -0,0 +1,76 @@
+/*
+ *  readpng.h
+ *
+ *  Copyright (C) 2003 Nadav Rotem <nadav256@hotmail.com>
+ *
+ *  This program is free software; you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation; either version 2 of the License, or
+ *  (at your option) any later version.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU Library General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program; if not, write to the Free Software
+ *  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
+ */
+
+/*
+ * Modified: Stefan van der Walt <stefan@sun.ac.za>
+ * Date: 28 January 2005
+ * - Fix bugs, restructure
+ */
+
+typedef struct
+{
+  int width;
+  int height;
+  int bit_depth;
+  int color_type;
+  unsigned char **row_pointers;
+} canvas;
+
+//////////////Libcanvas///////////
+canvas *new_canvas(int width, int height, int stride)
+{
+  // Default stride if none given
+  if (stride==0) stride=width*4;
+
+  // Clean allocation of canvas structure
+  canvas *can=new(canvas);
+  unsigned char *image_data = new unsigned char[stride*height];
+  unsigned char **row_pointers = new unsigned char *[height];
+  if (can == NULL || image_data == NULL || row_pointers == NULL) 
+    {
+      if (can == NULL) delete can;
+      if (image_data == NULL) delete[] image_data;
+      if (row_pointers == NULL) delete[] row_pointers;
+      return NULL;
+    }
+
+  // Fill in canvas structure
+  can->width=width;
+  can->height=height;
+  can->bit_depth=8;
+  can->color_type=PNG_COLOR_TYPE_RGB_ALPHA;
+  can->row_pointers = row_pointers;
+  for (int i=0; i < height; i++) row_pointers[i] = image_data + i*stride;
+
+  return can;
+}
+
+void delete_canvas(canvas *can)
+{
+  
+  if (can!=NULL)
+    {
+      delete[] can->row_pointers[0]; 
+      delete[] can->row_pointers;
+      delete can;
+    }
+  return;
+}
+
--- a/pngread.cc
+++ b/pngread.cc
@@ -1,5 +1,5 @@
 /*
- *  readpng.c
+ *  pngread.cc
  *
  *  Copyright (C) 2003 Nadav Rotem <nadav256@hotmail.com>
  *
@@ -18,10 +18,9 @@
  *  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
  */
 
-
 /*
 
-Loads PNG files to octave using libpng;
+Load PNG files to octave using libpng;
 
        PNG  (Portable  Network  Graphics) is an extensible file format for the
        lossless, portable, well-compressed storage of raster images. PNG  pro-
@@ -32,23 +31,27 @@
 
 */
 
-#include <png.h>
+/*
+ * Modified: Stefan van der Walt <stefan@sun.ac.za>
+ * Date: 28 January 2005
+ * - Fix bugs, restructure
+ */
+
 #include <octave/oct.h>
 
-typedef struct
-{
-  int width;
-  int height;
-  int bit_depth;
-  int color_type;
-  unsigned char **row_pointers;
-} canvas;
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+#include "png.h"
+
+#ifdef __cplusplus
+} //extern "C"
+#endif
+
+#include "pngcanvas.h"
 
 canvas *load_canvas(char *filename);
-void save_canvas(canvas *can,char *filename);
-canvas *new_canvas(int x, int y, int stride=0);
-void delete_canvas(canvas *can);
-
 
 DEFUN_DLD (pngread, args, nargout ,"\
 pngread reads a png file from disk.\n\
@@ -61,7 +64,7 @@
   //
   // We bail out if the input parameters are bad
   //
-  if (nargin != 1 || !args(0).is_string() || nargout < 3) {
+  if (nargin != 1 || !args(0).is_string()) {
     print_usage ("pngread");
     return retval;
   }
@@ -70,145 +73,60 @@
   // Load png file
   //
   canvas *pic=load_canvas((char *)args(0).string_value().c_str());
-  
-  // octave_stdout << "Canvas [" << pic->width << "x" << pic->height << "]\n";
- 
-  if (nargout >= 3) {
-    Matrix red   ( pic->height , pic->width );
-    Matrix green ( pic->height , pic->width );
-    Matrix blue  ( pic->height , pic->width );
-    Matrix alpha ( pic->height , pic->width );
-    
-    for (unsigned long j=0; j<pic->height; j++) 
-      for (unsigned long i=0; i<pic->width; i++) 
-	{
-	  red(j,i) = pic->row_pointers[j][i*4+0];
-	  green(j,i) = pic->row_pointers[j][i*4+1]           ;
-	  blue(j,i) = pic->row_pointers[j][i*4+2];
-	  alpha(j,i) = pic->row_pointers[j][i*4+3];
-	}
-    retval(3)= alpha;
-    retval(2)= blue;
-    retval(1)= green;
-    retval(0)= red;    
-  }
+  if (!pic) return retval;
+
+  dim_vector dim = dim_vector();
+  dim.resize(3);
+  dim(0) = pic->height;
+  dim(1) = pic->width;
+  dim(2) = 3;
 
-  return retval;
-}
-
+  if ( (pic->color_type == PNG_COLOR_TYPE_GRAY) ||
+       (pic->color_type == PNG_COLOR_TYPE_GRAY_ALPHA) ||
+       ((pic->color_type == PNG_COLOR_TYPE_PALETTE) && (pic->bit_depth == 1)) )
+      dim(2) = 1;
 
-DEFUN_DLD (pngwrite, args, ,"\
-pngwrite writes a png file to the disk.\n\
-    pngwrite('filename',R,G,B,A) writes the specified file\n\
-    using the Red, Green, Blue and Alpha matrices.\n\
-    \n\
-    Data must be [0 255] or the high bytes will be lost.")
-{
-   octave_value_list retval;
-   int nargin  = args.length();
-   
-   //
-   // We bail out if the input parameters are bad
-   //
-   if (nargin < 4 || !args(0).is_string() ) {
-     print_usage ("pngwrite");
-     return retval;
-   }
+  if (pic->bit_depth > 1 && pic->bit_depth < 8)
+      pic->bit_depth = 8;
 
-   Matrix red  = args(1).matrix_value();
-   Matrix green= args(2).matrix_value();
-   Matrix blue = args(3).matrix_value();
-   Matrix alpha= args(4).matrix_value();
-   
-   long image_width  = args(1).columns();
-   long image_height = args(1).rows();
-   
-   if ( args(2).columns() != image_width  ||
-	args(3).columns() != image_width  ||
-	args(4).columns() != image_width  ||
-	args(2).rows()    != image_height ||
-	args(3).rows()    != image_height ||
-	args(4).rows()    != image_height )  
-     {
-       error("pngwrite R,G,B,A matrix sizes aren't the same");
-       return retval;
-     }
-   
-   canvas *pic=new_canvas(image_width,image_height);
-   if (pic == NULL) 
-     {
-       error("pngwrite out of memory");
-       return retval;
-     }
-
-   // octave_stdout << "Canvas [" << pic->width << "x" << pic->height << "]\n";
-
-   for(int i=0; i<pic->width; i++) 
-     for(int j=0; j<pic->height; j++) 
-       {
-	 pic->row_pointers[j][i*4+0]=(unsigned char)(red(j,i));
-	 pic->row_pointers[j][i*4+1]=(unsigned char)(green(j,i));
-	 pic->row_pointers[j][i*4+2]=(unsigned char)(blue(j,i));
-	 pic->row_pointers[j][i*4+3]=(unsigned char)(alpha(j,i));
-       }
-   
-   save_canvas(pic,(char *)args(0).string_value().c_str());
-   delete_canvas(pic);
+  NDArray out = NDArray(dim, 0);
+  
+  Array<int> coord = Array<int> (3);
+  Matrix alpha ( pic->height , pic->width );
+  
+  for (unsigned long j=0; j < pic->height; j++) {
+      coord(0) = j;
+      for (unsigned long i=0; i < pic->width; i++) {
+	  coord(1) = i;
 
-   return retval;
-}
-
-//////////////Libcanvas///////////
-canvas *new_canvas(int width, int height, int stride)
-{
-  // Default stride if none given
-  if (stride==0) stride=width*4;
-
-  // Clean allocation of canvas structure
-  canvas *can=new(canvas);
-  unsigned char *image_data = new unsigned char[stride*height];
-  unsigned char **row_pointers = new unsigned char *[height];
-  if (can == NULL || image_data == NULL || row_pointers == NULL) 
-    {
-      if (can == NULL) delete can;
-      if (image_data == NULL) delete[] image_data;
-      if (row_pointers == NULL) delete[] row_pointers;
-      return NULL;
-    }
+	  for (int c = 0; c < dim(2); c++) {
+	      coord(2) = c;
+	      out(coord) = pic->row_pointers[j][i*3+c];
+	  }
+	  alpha(j,i) = pic->row_pointers[j][i*3+3];
+      }
+  }
+  out = out.squeeze();
 
-  // Fill in canvas structure
-  can->width=width;
-  can->height=height;
-  can->bit_depth=8;
-  can->color_type=PNG_COLOR_TYPE_RGB_ALPHA;
-  can->row_pointers = row_pointers;
-  for (int i=0; i < height; i++) row_pointers[i] = image_data + i*stride;
-
-  return can;
-}
-
+  switch (pic->bit_depth) {
+  case 1: retval.append((boolNDArray)out); break;
+  case 8: retval.append((uint8NDArray)out); break;
+  case 16: retval.append((uint16NDArray)out); break;
+  default: retval.append(out);
+  }
+  retval.append(alpha);
 
-
-void delete_canvas(canvas *can)
-{
-  
-  if (can!=NULL)
-    {
-      delete[] can->row_pointers[0]; 
-      delete[] can->row_pointers;
-      delete can;
-    }
-  return;
+  delete_canvas(pic);
+  return retval;
 }
 
 canvas *load_canvas(char *filename)
 {
-  FILE *infile;
   png_structp png_ptr;
   png_infop info_ptr;
   
-  infile = fopen(filename,"r");
-  if (NULL==infile) {
+  FILE *infile = fopen(filename,"r");
+  if (!infile) {
     error("pngread could not open file %s", filename); 
     return NULL;
   }
@@ -221,9 +139,9 @@
     return NULL;
   }
   
-  png_ptr=png_create_read_struct(PNG_LIBPNG_VER_STRING,NULL,NULL,NULL);
+  png_ptr = png_create_read_struct(PNG_LIBPNG_VER_STRING,NULL,NULL,NULL);
   if (!png_ptr) {
-    error("pngread out of mem"); 
+    error("pngread out of memory"); 
     fclose(infile);
     return NULL;
   }
@@ -235,96 +153,77 @@
     fclose(infile);
     return NULL;
   }
-  
-  // XXX FIXME XXX do we need to check for failure in any of these steps?
-  png_init_io(png_ptr,infile);
-  png_set_sig_bytes(png_ptr,8);
-  png_read_info(png_ptr,info_ptr);
+
+  /* Set error handling */
+  if (setjmp(png_jmpbuf(png_ptr))) {
+      error("pngread: libpng exited abnormally");
+      png_destroy_read_struct(&png_ptr, &info_ptr, png_infopp_NULL);
+      fclose(infile);
+      return NULL;
+  }
+
+  png_init_io(png_ptr, infile);
+  png_set_sig_bytes(png_ptr, 8);
+  png_read_info(png_ptr, info_ptr);
   
   png_uint_32 width,height;
-  int color_type,bit_depth;
+  int color_type, bit_depth;
   png_get_IHDR(png_ptr,info_ptr,&width,&height,
-	       &color_type,&bit_depth,NULL,NULL,NULL);
-  
-  // Set to RGB - ALPHA - 8bit
-  if (color_type == PNG_COLOR_TYPE_PALETTE)
-    png_set_palette_to_rgb(png_ptr);
+	       &bit_depth,&color_type,NULL,NULL,NULL);
   
-  if (color_type == PNG_COLOR_TYPE_GRAY && bit_depth < 8)
-    png_set_gray_1_2_4_to_8(png_ptr);
+  /* Transform grayscale images with depths < 8-bit to 8-bit, change
+   * paletted images to RGB, and add a full alpha channel if there is
+   * transparency information in a tRNS chunk.
+   */
+  if (color_type == PNG_COLOR_TYPE_PALETTE) {
+      png_set_palette_to_rgb(png_ptr);
+  }
+  if (color_type == PNG_COLOR_TYPE_GRAY && bit_depth < 8) {
+      png_set_gray_1_2_4_to_8(png_ptr);
+  }
+  if (png_get_valid(png_ptr, info_ptr, PNG_INFO_tRNS)) { //add alpha
+      png_set_tRNS_to_alpha(png_ptr);
+  }
   
-  if (png_get_valid(png_ptr, info_ptr,PNG_INFO_tRNS)) //add alpha
-    png_set_tRNS_to_alpha(png_ptr);
-  
-  if (bit_depth == 16)
-    png_set_strip_16(png_ptr);
-  
+  // Assume black background
+  if (color_type & PNG_COLOR_MASK_ALPHA) {
+      png_set_strip_alpha(png_ptr);
+  }
+
+  // Always transform image to RGB
   if (color_type == PNG_COLOR_TYPE_GRAY 
       || color_type == PNG_COLOR_TYPE_GRAY_ALPHA) 
-    png_set_gray_to_rgb(png_ptr);
+      png_set_gray_to_rgb(png_ptr);  
+
+  if (bit_depth < 8) {
+      png_set_packing(png_ptr);
+  }
   
-  if (color_type ==2)
-    png_set_filler(png_ptr,0xff, PNG_FILLER_AFTER);
-  
+  // For now, use 8-bit only
+  if (bit_depth == 16) {
+      png_set_strip_16(png_ptr);
+  }
+
   png_read_update_info(png_ptr,info_ptr);
 
+
   // Read the data from the file
-  int stride=png_get_rowbytes(png_ptr,info_ptr);
-  canvas *can = new_canvas(width,height,stride);
-  if (can != NULL) {
-    png_read_image(png_ptr,can->row_pointers);
+  int stride = png_get_rowbytes(png_ptr, info_ptr);
+  canvas *can = new_canvas(width, height, stride);
+  
+  if (can) {
+    png_read_image(png_ptr, can->row_pointers);
   } else {
-    error("pngread out of mem");
+    error("pngread out of memory");
   }
 
   png_read_end(png_ptr,NULL);
   png_destroy_read_struct(&png_ptr, &info_ptr, (png_infopp)NULL);
 
   fclose(infile);
-  
+
+  // Set color type and depth. Used to determine octave output arguments.
+  can->color_type = color_type;  
+  can->bit_depth = bit_depth;
   return can;
 }
-
-void save_canvas(canvas *can,char *filename)
-{
-  FILE            *fp;
-  png_structp     png_ptr;
-  png_infop       info_ptr;
-  
-  fp = fopen(filename, "wb");
-  if (fp == NULL) {
-    error("pngwrite could not open %s", filename);
-    return;
-  }
-
-  // XXX FIXME XXX do we need to check for failure in any of these steps?
-  png_ptr = png_create_write_struct(PNG_LIBPNG_VER_STRING, NULL, NULL, NULL);
-  info_ptr = png_create_info_struct(png_ptr);             
-  png_init_io(png_ptr, fp);
-  png_set_compression_level(png_ptr, 3);
-  
-  png_set_IHDR(png_ptr, info_ptr, can->width, can->height,         
-	       can->bit_depth, can->color_type, PNG_INTERLACE_NONE,
-	       PNG_COMPRESSION_TYPE_DEFAULT, PNG_FILTER_TYPE_DEFAULT);
-  
-  
-  png_set_gAMA(png_ptr, info_ptr, 0.7);
-  
-  time_t          gmt;
-  png_time        mod_time;
-  png_text        text_ptr[2];
-  time(&gmt);
-  png_convert_from_time_t(&mod_time, gmt);
-  png_set_tIME(png_ptr, info_ptr, &mod_time);
-  text_ptr[0].key = "Created by";
-  text_ptr[0].text = "Octave";
-  text_ptr[0].compression = PNG_TEXT_COMPRESSION_NONE;
-  
-  png_set_text(png_ptr, info_ptr, text_ptr, 1);
-  
-  png_write_info(png_ptr, info_ptr);
-  png_write_image(png_ptr, can->row_pointers);
-  png_write_end(png_ptr, info_ptr);                      
-  png_destroy_write_struct(&png_ptr, &info_ptr);          
-  fclose(fp); 
-}
new file mode 100644
--- /dev/null
+++ b/pngwrite.cc
@@ -0,0 +1,179 @@
+/*
+ *  pngwrite.cc
+ *
+ *  Copyright (C) 2003 Nadav Rotem <nadav256@hotmail.com>
+ *
+ *  This program is free software; you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation; either version 2 of the License, or
+ *  (at your option) any later version.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU Library General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program; if not, write to the Free Software
+ *  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
+ */
+
+/*
+
+Write PNG files to disk from octave
+
+       PNG  (Portable  Network  Graphics) is an extensible file format for the
+       lossless, portable, well-compressed storage of raster images. PNG  pro-
+       vides  a patent-free replacement for GIF and can also replace many com-
+       mon uses of TIFF. Indexed-color, grayscale, and  truecolor  images  are
+       supported,  plus  an optional alpha channel. Sample depths range from 1
+       to 16 bits.
+
+*/
+
+/*
+ * Modified: Stefan van der Walt <stefan@sun.ac.za>
+ * Date: 28 January 2005
+ * - Fix bugs, restructure
+ */
+
+#include <octave/oct.h>
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+#include "png.h"
+
+#ifdef __cplusplus
+} //extern "C"
+#endif
+
+#include "pngcanvas.h"
+
+void save_canvas(canvas *can, char *filename);
+
+DEFUN_DLD (pngwrite, args, ,"\
+pngwrite writes a png file to the disk.\n\
+    pngwrite('filename',R,G,B,A) writes the specified file\n\
+    using the Red, Green, Blue and Alpha matrices.\n\
+    \n\
+    Data must be [0 255] or the high bytes will be lost.")
+{
+   octave_value_list retval;
+   int nargin  = args.length();
+   
+   //
+   // We bail out if the input parameters are bad
+   //
+   if (nargin < 5 || !args(0).is_string() ) {
+     print_usage ("pngwrite");
+     return retval;
+   }
+
+   Matrix red  = args(1).matrix_value();
+   Matrix green= args(2).matrix_value();
+   Matrix blue = args(3).matrix_value();
+   Matrix alpha= args(4).matrix_value();
+   
+   long image_width  = args(1).columns();
+   long image_height = args(1).rows();
+   
+   if ( args(2).columns() != image_width  ||
+	args(3).columns() != image_width  ||
+	args(4).columns() != image_width  ||
+	args(2).rows()    != image_height ||
+	args(3).rows()    != image_height ||
+	args(4).rows()    != image_height )  
+     {
+       error("pngwrite R,G,B,A matrix sizes aren't the same");
+       return retval;
+     }
+   
+   octave_stdout << "we reached point A" << std::endl;
+   canvas *pic=new_canvas(image_width, image_height, image_width*4);
+   if (!pic) {
+       error("pngwrite out of memory");
+       return retval;
+   }
+
+   octave_stdout << "we reached point B" << std::endl;
+   for(int i=0; i < pic->width; i++) {
+       for(int j=0; j < pic->height; j++) {
+	   pic->row_pointers[j][i*4+0]=(unsigned char)(red(j,i));
+	   pic->row_pointers[j][i*4+1]=(unsigned char)(green(j,i));
+	   pic->row_pointers[j][i*4+2]=(unsigned char)(blue(j,i));
+	   pic->row_pointers[j][i*4+3]=(unsigned char)(alpha(j,i));
+       }
+   }
+   
+   octave_stdout << "we reached point C" << std::endl;
+
+   save_canvas(pic,(char *)args(0).string_value().c_str());
+   delete_canvas(pic);
+
+   return retval;
+}
+
+void save_canvas(canvas *can,char *filename)
+{
+  FILE            *fp;
+  png_structp     png_ptr;
+  png_infop       info_ptr;
+  
+  fp = fopen(filename, "wb");
+  if (fp == NULL) {
+    error("pngwrite could not open %s", filename);
+    return;
+  }
+
+  png_ptr = png_create_write_struct(PNG_LIBPNG_VER_STRING, NULL, NULL, NULL);
+  if (!png_ptr) {
+      fclose(fp);
+      error("pngwrite: cannot create write structure");
+      return;
+  }
+
+  info_ptr = png_create_info_struct(png_ptr);
+  if (!info_ptr) {
+      fclose(fp);
+      error("pngwrite: cannot not create image structure");
+      png_destroy_write_struct(&png_ptr, png_infopp_NULL);
+      return;
+  }
+
+  if (setjmp(png_jmpbuf(png_ptr))) {
+      fclose(fp);
+      png_destroy_write_struct(&png_ptr, &info_ptr);
+      error("pngread: libpng exited abnormally");
+      return;
+  }
+
+  png_init_io(png_ptr, fp);
+  png_set_compression_level(png_ptr, 3);
+  
+  png_set_IHDR(png_ptr, info_ptr, can->width, can->height,         
+	       can->bit_depth, can->color_type, PNG_INTERLACE_NONE,
+	       PNG_COMPRESSION_TYPE_DEFAULT, PNG_FILTER_TYPE_DEFAULT);
+  
+  
+  png_set_gAMA(png_ptr, info_ptr, 0.7);
+  
+  time_t          gmt;
+  png_time        mod_time;
+  png_text        text_ptr[2];
+  time(&gmt);
+  png_convert_from_time_t(&mod_time, gmt);
+  png_set_tIME(png_ptr, info_ptr, &mod_time);
+  text_ptr[0].key = "Created by";
+  text_ptr[0].text = "Octave";
+  text_ptr[0].compression = PNG_TEXT_COMPRESSION_NONE;
+  
+  png_set_text(png_ptr, info_ptr, text_ptr, 1);
+  
+  png_write_info(png_ptr, info_ptr);
+  png_write_image(png_ptr, can->row_pointers);
+  png_write_end(png_ptr, info_ptr);                      
+  png_destroy_write_struct(&png_ptr, &info_ptr);          
+  fclose(fp); 
+}