v4.4+ supports images, but accepts 3-typles (x,y,value) while essentially ignoring the (x,y) part. The documentation indicates that "It is assumed that in the viewing plane the image data forms an equidistant sampling grid in the viewing plane along two, not necessarily orthogonal, directions.". If making such a drastic assumption, it doesn't make sense to accept the full (x,y) column values at all -- for a 2k x 2k image with a 3-tuple, the pixel shape information is oversampled by a factor of 2 million.
Clearly this needs attention for the future. I came across this particular bug because, while working on the new PDL::Graphics::Gnuplot interface module to PDL (http://pdl.perl.org), I wanted to adapt the FITS WCS display system that is currently used with PGPLOT. WCS is the "World Coordinate System", and it follows a particular standard for communicating how scientific image coordinates (e.g. "arcseconds" in a telescope image plane or "meters" on a map) map to pixel coordinates. PGPLOT uses the parallelogram assumption as does Gnuplot, but allows the user to send a transform matrix rather than laboriously enumerating the (ignored) position of every pixel. Contrariwise, allowing fully general grid layout (subject to the usual caveats that stupid grids will probably yield stupid results) would allow fully general nonlinear transformations of images without pre-render resampling, which would be a big win.