On 11-02-08 10:03 PM, Simon Urbanek wrote:
> Ben,
>
> did you actually look at the result of your function with useRaster=TRUE ? ;)
> [Hint: don't use an image that is symmetric]
>
> Apart from that nice bug there are more issues as well, try
> image(matrix(1:4,2),col=1:3)
> The underlying issue is that as.raster() is not quite what you would hope.
> Unfortunately I'm not aware of an easy fix (that doesn't involve going
back to RGB decomposition).
>
> In general, I think it's a nice option, but I don't think you'll get away
> with only a few lines...
>
> Cheers,
> Simon
>
>
>
> On Feb 8, 2011, at 8:49 PM, Ben Bolker wrote:
>
>>
>> Has anyone yet tried incorporating rasterImage into the base image()
>> function? It seems to make a *huge* difference, with
>> a very small number of added/changed lines of code. (Of course I have
>> barely tested it at all.)
>> Is there any reason this *shouldn't* go into the next release?
>>
>>> source("image.R")
>>> z <- matrix(runif(1e6),nrow=1000)
>>> image(z)
>>> image(z,useRaster=TRUE)
>>
>> (Patch against SVN 54284 attached; people can contact me if it doesn't
>> go through and they want it.)
>>
>> Ben Bolker
>>
>> <image_diff.txt>______________________________________________
>> [email protected] mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-devel
Trying again. Rotated counterclockwise within R (although this *could*
be coded in C if speed were important?)
Some brute-force testing suggests it is *slightly* slower for small
images (7 vs 8 seconds for 1000 reps) and still much faster (and
produces much smaller images, which don't suffer from antialiasing junk
in my PDF viewer) for large images.
Index: R/image.R
===================================================================
--- R/image.R (revision 54297)
+++ R/image.R (working copy)
@@ -24,7 +24,8 @@
ylim = range(y),
col = heat.colors(12), add = FALSE,
xaxs = "i", yaxs = "i", xlab, ylab,
- breaks, oldstyle=FALSE, ...)
+ breaks, oldstyle=FALSE,
+ useRaster=FALSE, ...)
{
if (missing(z)) {
if (!missing(x)) {
@@ -97,5 +98,21 @@
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)")
- .Internal(image(as.double(x), as.double(y), as.integer(zi), col))
+ if (useRaster) {
+ if (is.numeric(col)) {
+ p <- palette()
+ pl <- length(p)
+ col <- p[((col-1) %% pl)+1]
+ }
+ zc <- col[zi+1]
+ dim(zc) <- dim(z)
+ ##Rotate a square matrix 90 degrees COUNTERclockwise.
+ ## inspired by
+ ##
<http://tuoq.blogspot.com/2008/03/rotate-matrix-90-degrees-clockwise.html>
+ zc <- t(zc)[ncol(zc):1,]
+ rasterImage(as.raster(zc),
+ xlim[1],ylim[1],xlim[2],ylim[2],interpolate=FALSE)
+ } else {
+ .Internal(image(as.double(x), as.double(y), as.integer(zi), col))
+ }
}
require(grDevices) # for colours
## read modified image.R file, from wherever, if necessary
## source("~/R/r-devel/R/src/library/graphics/R/image.R")
tmpf <- function(ifun=image) {
ifun(matrix(1:25,5),col=1:25)
x <- y <- seq(-4*pi, 4*pi, len=27)
r <- sqrt(outer(x^2, y^2, "+"))
ifun(z = z <- cos(r^2)*exp(-r/6), col=gray((0:32)/32))
ifun(z, axes = FALSE, main = "Math can be beautiful ...",
xlab = expression(cos(r^2) * e^{-r/6}))
contour(z, add = TRUE, drawlabels = FALSE)
## Volcano data visualized as matrix. Need to transpose and flip
## matrix horizontally.
ifun(t(volcano)[ncol(volcano):1,])
## A prettier display of the volcano
x <- 10*(1:nrow(volcano))
y <- 10*(1:ncol(volcano))
ifun(x, y, volcano, col = terrain.colors(100), axes = FALSE)
contour(x, y, volcano, levels = seq(90, 200, by = 5),
add = TRUE, col = "peru")
axis(1, at = seq(100, 800, by = 100))
axis(2, at = seq(100, 600, by = 100))
box()
title(main = "Maunga Whau Volcano", font.main = 4)
}
pdf("image1.pdf")
tmpf()
dev.off()
pdf("image2.pdf")
tmpf(ifun=function(...) image(...,useRaster=TRUE))
dev.off()
st1 <- system.time(replicate(1000,image(matrix(1:25,5),col=1:25)))
st2 <- system.time(replicate(1000,image(matrix(1:25,5),col=1:25,
useRaster=TRUE)))
## 21.3 (standard) vs 22.9 (raster) seconds
st3 <- system.time(replicate(1000,image(matrix(1:10000,100))))
st4 <- system.time(replicate(1000,image(matrix(1:10000,100),
useRaster=TRUE)))
## 149 (standard) vs 31 (raster) seconds
# File src/library/graphics/R/image.R
# Part of the R package, http://www.R-project.org
#
# 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 General Public License for more details.
#
# A copy of the GNU General Public License is available at
# http://www.r-project.org/Licenses/
image <- function(x, ...) UseMethod("image")
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)) { # midpoints
dx <- 0.5*diff(x)
x <- c(x[1L] - dx[1L], x[-length(x)]+dx,
x[length(x)]+dx[length(x)-1])
}
if (length(y) > 1 && length(y) == ncol(z)) { # midpoints
dy <- 0.5*diff(y)
y <- c(y[1L] - dy[1L], 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[1L] == 0) c(-1, 1)
else zlim[1L] + c(-.4, .4)*abs(zlim[1L])
z <- (z - zlim[1L])/diff(zlim)
zi <- if (oldstyle) floor((nc - 1) * z + 0.5)
else floor((nc - 1e-5) * z + 1e-7)
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, ...)
## need plot set up before we do this
if (length(x) <= 1) x <- par("usr")[1L: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) {
if (is.numeric(col)) {
p <- palette()
pl <- length(p)
col <- p[((col-1) %% pl)+1]
}
zc <- col[zi+1]
dim(zc) <- dim(z)
##Rotate a square matrix 90 degrees COUNTERclockwise.
## inspired by
##
<http://tuoq.blogspot.com/2008/03/rotate-matrix-90-degrees-clockwise.html>
zc <- t(zc)[ncol(zc):1,]
rasterImage(as.raster(zc),
xlim[1],ylim[1],xlim[2],ylim[2],interpolate=FALSE)
} else {
.Internal(image(as.double(x), as.double(y), as.integer(zi), col))
}
}
______________________________________________
[email protected] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
