# R matrices and image

code
Published

2014-04-17

## Matrix and image

In R, matrices are ordered row-wise:

``(m <- matrix(1:12, nrow = 3, ncol = 4))``
``````     [,1] [,2] [,3] [,4]
[1,]    1    4    7   10
[2,]    2    5    8   11
[3,]    3    6    9   12``````

The image() function presents this as the transpose of what we see printed.

``````m[] <- 0
m[2, 1] <- -10
m[3, 2] <- 30
m``````
``````     [,1] [,2] [,3] [,4]
[1,]    0    0    0    0
[2,]  -10    0    0    0
[3,]    0   30    0    0``````
``t(m[, ncol(m):1])``
``````     [,1] [,2] [,3]
[1,]    0    0    0
[2,]    0    0    0
[3,]    0    0   30
[4,]    0  -10    0``````

… Notice that image interprets the z matrix as a table of f(x[i], y[j]) values, so that the x axis corresponds to row number and the y axis to column number, with column 1 at the bottom, i.e. a 90 degree counter-clockwise rotation of the conventional printed layout of a matrix. …

``image(m)`` ## Data placement with image()

This is fairly obvious, each cell is painted as a discrete block with cell centres evenly spaced between 0 and 1.

``````m <- matrix(1:12, 3)
image(m)`````` We didn’t give it any coordinates to position the image, so it made some up.

``````image(m, main = "input coordinates are cell centres")
xx <- seq.int(0, 1, length.out = nrow(m))
yy <- seq.int(0, 1, length.out = ncol(m))
abline(h = yy, v = xx, lty = 2)`````` This lends itself to a convenient data structure.

``````dat <- list(x = xx, y = yy, z = m)
image(dat)
text(expand.grid(xx, yy), lab = as.vector(m))`````` ``## points(expand.grid(xx, yy))``

The function image() has some hidden tricks.

``````xcorner <- seq.int(0, 1, length.out = nrow(m) + 1L)
ycorner <- seq.int(0, 1, length.out = ncol(m) + 1L)
print(xcorner)``````
`` 0.0000000 0.3333333 0.6666667 1.0000000``
``print(ycorner)``
`` 0.00 0.25 0.50 0.75 1.00``
``##  0.00 0.25 0.50 0.75 1.00``
``````image(xcorner, ycorner, m, main = "input coordinates are cell corners")
abline(h = ycorner, v = xcorner)`````` We can even use non-regular coordinates.

``````ycorner <- 1.5^seq_along(ycorner)
image(xcorner, ycorner, m)
abline(h = ycorner, v = xcorner)`````` ## Under the hood

``print(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 = hcl.colors(12, "YlOrRd",
rev = TRUE), add = FALSE, xaxs = "i", yaxs = "i", xlab,
ylab, breaks, oldstyle = FALSE, useRaster, ...)
{
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 <- deparse1(substitute(x))
if (missing(xlab))
xlab <- paste0(xn, "\$x")
if (missing(ylab))
ylab <- paste0(xn, "\$y")
y <- x\$y
x <- x\$x
}
else {
if (missing(xlab))
xlab <- if (missing(x))
""
else deparse1(substitute(x))
if (missing(ylab))
ylab <- if (missing(y))
""
else deparse1(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 (!typeof(z) %in% c("logical", "integer", "double"))
stop("'z' must be numeric or logical")
if (length(x) > 1 && length(x) == nrow(z)) {
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)) {
dy <- 0.5 * diff(y)
y <- c(y[1L] - dy[1L], y[-length(y)] + dy, y[length(y)] +
dy[length(y) - 1L])
}
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(-0.4, 0.4) * abs(zlim[1L])
z <- (z - zlim[1L])/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")
if (is.unsorted(breaks)) {
warning("unsorted 'breaks' will be sorted before use")
breaks <- sort(breaks)
}
zi <- .bincode(z, breaks, TRUE, TRUE) - 1L
}
plot(xlim, ylim, xlim = xlim, ylim = ylim, type = "n",
xaxs = xaxs, yaxs = yaxs, xlab = xlab, ylab = ylab,
...)
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)")
check_irregular <- function(x, y) {
dx <- diff(x)
dy <- diff(y)
(length(dx) && !isTRUE(all.equal(dx, rep(dx, length(dx))))) ||
(length(dy) && !isTRUE(all.equal(dy, rep(dy, length(dy)))))
}
if (missing(useRaster)) {
useRaster <- getOption("preferRaster", FALSE)
if (useRaster && check_irregular(x, y))
useRaster <- FALSE
if (useRaster) {
useRaster <- FALSE
ras <- dev.capabilities("rasterImage")\$rasterImage
if (identical(ras, "yes"))
useRaster <- TRUE
if (identical(ras, "non-missing"))
useRaster <- all(!is.na(zi))
}
}
if (useRaster) {
if (check_irregular(x, y))
stop(gettextf("%s can only be used with a regular grid",
sQuote("useRaster = TRUE")), domain = NA)
if (!is.character(col)) {
col <- as.integer(col)
if (any(!is.na(col) & col < 0L))
stop("integer colors must be non-negative")
col[col < 1L] <- NA_integer_
p <- palette()
col <- p[((col - 1L)%%length(p)) + 1L]
}
zc <- col[zi + 1L]
dim(zc) <- dim(z)
zc <- t(zc)[ncol(zc):1L, , drop = FALSE]
rasterImage(as.raster(zc), min(x), min(y), max(x), max(y),
interpolate = FALSE)
}
else .External.graphics(C_image, x, y, zi, col)
invisible()
}
<bytecode: 0x55b5b27527a8>
<environment: namespace:graphics>``````

This is like looping with rect()

``````op <- par(mfrow = c(1, 2))
## life is hard
cols <- topo.colors(25)
scale <- round((m - min(m))/diff(range(m)) * (length(cols) - 1) + 1)
plot(NA, type = "n", xlim = range(xcorner), ylim = range(ycorner), asp = 1)
for (i in seq_along(xcorner[-1L])) {
for (j in seq_along(ycorner[-1L])) {
rect(xleft = xcorner[i], ybottom = ycorner[j], xright = xcorner[i +
1L], ytop = ycorner[j + 1L], col = cols[scale[i, j]], angle = 45 *
(i + j)%%2, density = 20, lwd = 2)
}

}

## life is good
image(list(x = xcorner, y = ycorner, z = m), col = topo.colors(25), asp = 1)`````` ``par(op)``

## “Raster graphics” (not the raster package)

Relatively recently native image-graphics support was added to R.

Old style

``````m <- matrix(1:12, nrow = 3)
xcorner <- seq.int(0, 1, length.out = nrow(m) + 1L)
ycorner <- seq.int(0, 1, length.out = ncol(m) + 1L)
image(xcorner, ycorner, m, col = topo.colors(25))`````` 