From lachmann at eva.mpg.de  Wed Aug  3 23:51:27 2011
From: lachmann at eva.mpg.de (Michael Lachmann)
Date: Wed, 3 Aug 2011 23:51:27 +0200
Subject: [Rd] one way to solve bad looking density plots in postscript
In-Reply-To: <CANLFJPp1N0YWq5u3H91p+z=BLyZKvyRF9obKWHzL4710OyJniw@mail.gmail.com>
References: <9318F8C2-1470-4AFB-A9E1-EC9CBF940150@eva.mpg.de>
	<64058E1C-3122-4D67-9894-897EEFE4F9E2@eva.mpg.de>
	<CANLFJPp1N0YWq5u3H91p+z=BLyZKvyRF9obKWHzL4710OyJniw@mail.gmail.com>
Message-ID: <80E0BF0D-20C0-4934-A775-4B0649DF3822@eva.mpg.de>

Sorry, I must have made a mistake before.
In my R, the attached replacement for image.default DOES put the axes last, and does look better. 

I also tested my code for the change in the eps header on some linux machines (I'm using a mac), and there the grid effects were less pronounced, and my code did not help...

anyway, here is the plot.default function that does seem to put the box after the image is drawn.


Michael


The only change is adding box(...) at the end, and doing the first call to plot with bty="n".
--
image.default = function (x = seq(0, 1, length.out = nrow(z)), y = seq(0, 1, 
    length.out = ncol(z)), z, zlim = range(z[is.finite(z)]), 
    xlim = range(x), ylim = range(y), col = heat.colors(12), 
    add = FALSE, xaxs = "i", yaxs = "i", xlab, ylab, breaks, 
    oldstyle = FALSE, useRaster = FALSE, ...) 
{
    if (missing(z)) {
        if (!missing(x)) {
            if (is.list(x)) {
                z <- x$z
                y <- x$y
                x <- x$x
            }
            else {
                if (is.null(dim(x))) 
                  stop("argument must be matrix-like")
                z <- x
                x <- seq.int(0, 1, length.out = nrow(z))
            }
            if (missing(xlab)) 
                xlab <- ""
            if (missing(ylab)) 
                ylab <- ""
        }
        else stop("no 'z' matrix specified")
    }
    else if (is.list(x)) {
        xn <- deparse(substitute(x))
        if (missing(xlab)) 
            xlab <- paste(xn, "x", sep = "$")
        if (missing(ylab)) 
            ylab <- paste(xn, "y", sep = "$")
        y <- x$y
        x <- x$x
    }
    else {
        if (missing(xlab)) 
            xlab <- if (missing(x)) 
                ""
            else deparse(substitute(x))
        if (missing(ylab)) 
            ylab <- if (missing(y)) 
                ""
            else deparse(substitute(y))
    }
    if (any(!is.finite(x)) || any(!is.finite(y))) 
        stop("'x' and 'y' values must be finite and non-missing")
    if (any(diff(x) <= 0) || any(diff(y) <= 0)) 
        stop("increasing 'x' and 'y' values expected")
    if (!is.matrix(z)) 
        stop("'z' must be a matrix")
    if (length(x) > 1 && length(x) == nrow(z)) {
        dx <- 0.5 * diff(x)
        x <- c(x[1] - dx[1], x[-length(x)] + dx, x[length(x)] + 
            dx[length(x) - 1])
    }
    if (length(y) > 1 && length(y) == ncol(z)) {
        dy <- 0.5 * diff(y)
        y <- c(y[1] - dy[1], y[-length(y)] + dy, y[length(y)] + 
            dy[length(y) - 1])
    }
    if (missing(breaks)) {
        nc <- length(col)
        if (!missing(zlim) && (any(!is.finite(zlim)) || diff(zlim) < 
            0)) 
            stop("invalid z limits")
        if (diff(zlim) == 0) 
            zlim <- if (zlim[1] == 0) 
                c(-1, 1)
            else zlim[1] + c(-0.4, 0.4) * abs(zlim[1])
        z <- (z - zlim[1])/diff(zlim)
        zi <- if (oldstyle) 
            floor((nc - 1) * z + 0.5)
        else floor((nc - 1e-05) * z + 1e-07)
        zi[zi < 0 | zi >= nc] <- NA
    }
    else {
        if (length(breaks) != length(col) + 1) 
            stop("must have one more break than colour")
        if (any(!is.finite(breaks))) 
            stop("breaks must all be finite")
        zi <- .C("bincode", as.double(z), length(z), as.double(breaks), 
            length(breaks), code = integer(length(z)), (TRUE), 
            (TRUE), nok = TRUE, NAOK = TRUE, DUP = FALSE, PACKAGE = "base")$code - 
            1
    }
    if (!add) 
        plot(NA, NA, xlim = xlim, ylim = ylim, type = "n", xaxs = xaxs, 
            yaxs = yaxs, xlab = xlab, ylab = ylab,bty="n", ...)
    if (length(x) <= 1) 
        x <- par("usr")[1:2]
    if (length(y) <= 1) 
        y <- par("usr")[3:4]
    if (length(x) != nrow(z) + 1 || length(y) != ncol(z) + 1) 
        stop("dimensions of z are not length(x)(-1) times length(y)(-1)")
    if (useRaster) {
        dx <- diff(x)
        dy <- diff(y)
        if ((length(dx) && !isTRUE(all.equal(dx, rep(dx[1], length(dx))))) || 
            (length(dy) && !isTRUE(all.equal(dy, rep(dy[1], length(dy)))))) 
            stop("useRaster=TRUE can only be used with a regular grid")
        if (!is.character(col)) {
            p <- palette()
            pl <- length(p)
            col <- as.integer(col)
            col[col < 1] <- NA
            col <- p[((col - 1)%%pl) + 1]
        }
        zc <- col[zi + 1]
        dim(zc) <- dim(z)
        zc <- t(zc)[ncol(zc):1, ]
        rasterImage(as.raster(zc), min(x), min(y), max(x), max(y), 
            interpolate = FALSE)
    }
    else .Internal(image(as.double(x), as.double(y), as.integer(zi), 
        col))
    box(...)
}

--

On 3 Aug 2011, at 11:14PM, baptiste auguie wrote:

> Hi,
> 
> I was going to suggest panel.last to redraw a box around the image,
> but for some reason it does not seem to come after the image neither
> (this is perhaps to be expected from the note in ?plot.default).
> 
> plot(1,1, panel.last={box(lwd=50, col="#0000FF")})
> 
> image(volcano, panel.last={box(lwd=50, col="#0000FF")})
> 
> sessionInfo()
> R version 2.13.1 (2011-07-08)
> Platform: x86_64-apple-darwin9.8.0/x86_64 (64-bit)
> 
> locale:
> [1] en_NZ.UTF-8/en_NZ.UTF-8/C/C/en_NZ.UTF-8/en_NZ.UTF-8
> 
> attached base packages:
> [1] stats     graphics  grDevices utils     datasets  grid
> methods   base
> 
> other attached packages:
> [1] ggplot2_0.8.9 proto_0.3-9.2 reshape_0.8.4 plyr_1.6
> 
> loaded via a namespace (and not attached):
> [1] tools_2.13.1
> 
> Best,
> 
> baptiste
> 
> On 4 August 2011 00:23, Michael Lachmann <lachmann at eva.mpg.de> wrote:
>> Some more digging.
>> 
>> 1.
>> The following code will fix density plots for me:
>> ---
>> .ps.prolog=grDevices:::.ps.prolog
>> i=grep("/p2",.ps.prolog)
>> .ps.prolog[i] = "/p2 { bg gsave fill grestore 0.001 setlinewidth stroke
>> newpath } def"
>> ---
>> 
>> 2.
>> It seems that it doesn't matter where in image.default() the box is drawn, before or after
>> .Internal(image(as.double(x), as.double(y), as.integer(zi),
>>        col))
>> is called, the resulting eps has the box in front.
>> 
>> Michael
>> 
>> On 3 Aug 2011, at 9:25AM, Michael Lachmann wrote:
>> 
>>> When R generates density plots and these are exported to postscript(
>>> a=matrix(1:100,10,10);image(a,col=rainbow(100);dev.copy2eps(file="image.eps")
>>> )
>>> The result often looks bad when rendered on screen. The help page states that this is because programs use anti-aliasing. That seems to be true - turning off anti-aliasing for gs (-dGraphicsAlphaBits=1) of in OSX's preview makes the plots look really smooth, but makes everything else look bad (personal opinion..). The plots do look ok in acrobat reader. I think it would be much better if R corrected this problem - even if it is not totally R's fault.
>>> It seems that using the option useRaster=T in image() solved this problem, but creates other problems for OSX's Preview (it seems that OSX's preview first anti-aliases the raster, and then scales it... creating a mess).
>>> 
>>> Density plots produced by gnuplot do not seem to have this problem:
>>> ---
>>> set pm3d map
>>> set pm3d at b
>>> set ticslevel 0.8
>>> set isosample 40,40
>>> set output "gtest.eps"
>>> set term postscript eps color
>>> splot [-3:3] [-3:3] x*x*exp(-x*x)*y*y*exp(-y*y)
>>> --
>>> But I haven't figured out why that is. Maybe someone who understands more about postscript can. Maybe it is something about the order that the rectangles are rendered? I did notice that rectangles are plotted with slightly different sizes - 50, 51, 50, 51 and so on. Is that it?
>>> 
>>> After a lot of experimentation, I found that a small change in the eps file can correct the output.
>>> If in the eps file produced above, you change the line
>>> /p2  { gsave bg fill grestore newpath } def
>>> with
>>> /p2  { bg gsave fill grestore stroke newpath } def
>>> and add
>>> 0.0001 setlinewidth
>>> two lines before the next p2 - i.e. before we start plotting the rectangles of the image:
>>> change:
>>> /bg { 1 0 0 setrgb } def
>>> 59.04 73.44 41.47 37.15 r p2
>>> to:
>>> 0.0001 setlinewidth
>>> /bg { 1 0 0 setrgb } def
>>> 59.04 73.44 41.47 37.15 r p2
>>> 
>>> What this does is plot the outline of each rectangle in addition to filling it.
>>> The page at:
>>> http://pages.uoregon.edu/noeckel/MathematicaGraphics.html#ExportGraphics
>>> Claims that Mathematica also has/had this problem, and how it can be solved there (sadly, I don't know enough Mathematica to understand that solution)
>>> 
>>> The output of the image would look even better if the axes were drawn AFTER the density rectangles, not before. That would cause the rectangles not to overwrite part of the lines of the axes. But that is probably a change in the image() routine, not in the postscript driver....
>>> 
>>> Thanks for listening,
>>> 
>>> Michael Lachmann
>>> 
>>> ______________________________________________
>>> R-devel at r-project.org mailing list
>>> https://stat.ethz.ch/mailman/listinfo/r-devel
>>> 
>>> 
>>> 
>> 
>> ______________________________________________
>> R-devel at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-devel
>> 
> 
> 
> 



[Rd] one way to solve bad looking density plots in postscript
Michael Lachmann