42

In R, it seems none of the currently available options (e.g. image) allow for fast real-time display of 2D raster graphics. I was interested for example to make a real-time interactive Mandelbrot viewer, where I was hoping to display 1920x1080 raw hex color matrices as raster images to achieve ca. 5-10 fps (calculating the Mandelbrot images themselves now achieves ca. 20-30 fps at moderate zooms). Or even better I would like to be able to pass an intensity pixel map to the graphics card & have the intensity to colour mapping also done using OpenGL acceleration. Using image() with option useRaster=TRUE, plot.raster or even grid.raster() doesn't cut it as displaying the raster image takes ca. 1/4 of a second (more than actually calculating the frame), so I am on the lookout for a more performant option, ideally using SDL or OpenGL acceleration.

I noticed that one should be able to call SDL, SDL_image and GL/OpenGL functions from R using the rdyncall foreign function interface package, which should have much better performance.

Although this package is archived on CRAN, an updated version is available here (and we are trying to get it back on CRAN)

library(remotes)
remotes::install_github("hongyuanjia/rdyncall")
library(rdyncall)

The SDL, SDL_image and SDL_mixer DLLs (version 1.2) (on Windows) have to be installed first from https://libsdl.org/release/, https://www.libsdl.org/projects/SDL_image/release/ and https://www.libsdl.org/projects/SDL_mixer/release/)(the 64 bit DLLs are to be put underR-4.2.1/bin/x64`) or using

source("https://raw.githubusercontent.com/Jean-Romain/lidRviewer/master/sdl.R")

On Ubuntu they can be installed using

sudo apt-get install libsdl1.2-dev libsdl-image1.2-dev libsdl-mixer1.2

Some demos of how to make SDL & OpenGL calls are available at https://dyncall.org/demos/soulsalicious/index.html (1980s computer-game style starfield, with music included).

Am I correct that with this package one should be able to display a 2D image raster using SDL & opengl acceleration? If so, has anyone any thoughts how to do this (I'm asking because I have no experience in using either SDL or OpenGL)?

To open a 1920 x 1080 SDL window I think I have to use (gathered from some OpenGL examples and windowed.R script in https://dyncall.org/demos/soulsalicious/soulsalicious.tar.gz, fullscreen also possible, see fullscreen.R)

init <- function()
{
  require(rdyncall)
  dynport(SDL)
  SDL_Init(SDL_INIT_VIDEO)
  dynport(GL)
  dynport(GLU)
  dynport(SDL_image)
  
  SDL_GL_SetAttribute(SDL_GL_RED_SIZE,8)
  SDL_GL_SetAttribute(SDL_GL_GREEN_SIZE,8)
  SDL_GL_SetAttribute(SDL_GL_BLUE_SIZE,8)
  SDL_GL_SetAttribute(SDL_GL_DOUBLEBUFFER,1)
  x_res <- 1920
  y_res <- 1080
  win <- SDL_SetVideoMode(x_res, y_res, 32, 
                          SDL_HWSURFACE + SDL_OPENGL + SDL_DOUBLEBUF)
  SDL_WM_SetCaption("SDL surface",NULL)
  glEnable(GL_TEXTURE_2D)
  # Set the projection matrix for the image
  glMatrixMode(GL_PROJECTION)
  glLoadIdentity()
  x_min=1
  x_max=x_res
  y_min=1
  y_max=y_res
  glOrtho(x_min, x_max, y_min, y_max, -1, 1)
  # Set the modelview matrix for the image
  glMatrixMode(GL_MODELVIEW)
  glLoadIdentity()
}
init()

I gather I should then set up some pixel transfer map using something like

glPixelMapfv(GL_PIXEL_MAP_I_TO_R, nb_colors, map_colors)
glPixelMapfv(GL_PIXEL_MAP_I_TO_G, nb_colors, map_colors)
glPixelMapfv(GL_PIXEL_MAP_I_TO_B, nb_colors, map_colors)

then create a buffer for the pixel data using another pixels <- glPixelMapfv call & draw the pixel data to the screen using glDrawPixels and swap the back and front buffers to display the image using SDL_GL_SwapBuffers(win) and then wait for the user to close the window & then clean up using SDL_Quit() etc. Trouble is I have no OpenGL or SDL experience, so would anybody know how to carry out these last few steps? (I am using SDL version 1.2 here)

Some timings of non-OpenGL options which are too slow for my application:

# some example data & desired colour mapping of [0-1] ranged data matrix
library(RColorBrewer)
ncol=1080
cols=colorRampPalette(RColorBrewer::brewer.pal(11, "RdYlBu"))(ncol)
colfun=colorRamp(RColorBrewer::brewer.pal(11, "RdYlBu"))
col = rgb(colfun(seq(0,1, length.out = ncol)), max = 255)
mat=matrix(seq(1:1080)/1080,nrow=1920,ncol=1080,byrow=TRUE)
mat2rast = function(mat, col) {
  idx = findInterval(mat, seq(0, 1, length.out = length(col)))
  colors = col[idx]
  rastmat = t(matrix(colors, ncol = ncol(mat), nrow = nrow(mat), byrow = TRUE))
  class(rastmat) = "raster"
  return(rastmat)
}
system.time(mat2rast(mat, col)) # 0.24s

# plot.raster method - one of the best?
par(mar=c(0, 0, 0, 0))
system.time(plot(mat2rast(mat, col), asp=NA)) # 0.26s

# grid graphics - tie with plot.raster?
library(grid)
system.time(grid.raster(mat2rast(mat, col),interpolate=FALSE)) # 0.28s

# base R image()
par(mar=c(0, 0, 0, 0))
system.time(image(mat,axes=FALSE,useRaster=TRUE,col=cols)) # 0.74s # note Y is flipped to compared to 2 options above - but not so important as I can fill matrix the way I want

# ggplot2 - just for the record...
df=expand.grid(y=1:1080,x=1:1920)
df$z=seq(1,1080)/1080
library(ggplot2)
system.time({q <- qplot(data=df,x=x,y=y,fill=z,geom="raster") + 
                scale_x_continuous(expand = c(0,0)) + 
                scale_y_continuous(expand = c(0,0)) +
                scale_fill_gradientn(colours = cols) + 
                theme_void() + theme(legend.position="none"); print(q)}) # 11s 
Tom Wenseleers
  • 7,535
  • 7
  • 63
  • 103
  • Try quadmesh package for conversion of raster to efficient rgl form. – mdsumner Jan 06 '18 at 03:57
  • Well that doesn't solve my problem I reckon - I don't want to get a 3D height map of my matrix, but just a 2D colour-mapped levelplot/heatmap/raster. – Tom Wenseleers Jan 06 '18 at 04:06
  • 1
    You can set z = 0. I'm not suggesting it solves your problem! I don't know the solution and am keen to find one too. Using raster or GDAL to provide level-of-detail specific to the current window is about the best I've seen, but it's hard to string the pieces together – mdsumner Jan 06 '18 at 04:34
  • btw, I don't think you are driving magick here, there's no conversion to that format (surely plot doesn't do that implicitly??) – mdsumner Jan 06 '18 at 04:38
  • This invokes magick in the ways I understand. I don't otherwise know how to convert array/matrix data to magick, though this is only one of many file-based ways. https://gist.github.com/mdsumner/7751139dccfb63d648c8e60da158529b – mdsumner Jan 06 '18 at 05:01
  • Ha sorry I see now - is there no way to get raster data into magick without going to file first, because that's going to be really inefficient? – Tom Wenseleers Jan 06 '18 at 12:00
  • It's the first thing I tried, but still not something I've been able to put enough thought into to chase down :| – mdsumner Jan 06 '18 at 13:13
  • Ah, no, this works - it's a browser Viewer though, not the native device image_read(mat2rast(mat2, col)) – mdsumner Jan 06 '18 at 13:26
  • Ha thanks for clearing that up! – Tom Wenseleers Jan 06 '18 at 13:34
  • @TomWenseleers It's been almost an year since your last edit to your question, meanwhile did you manage to do this? I am also interested... – Ramiro Magno Nov 18 '18 at 14:46
  • 1
    No I had been kind of hoping for an answer here, but unfortunately I didn't get one so far... – Tom Wenseleers Nov 19 '18 at 08:45

2 Answers2

1

According to the RGL package introduction, it is:

a visualization device system for R, using OpenGL as the rendering backend. An rgl device at its core is a real-time 3D engine written in C++. It provides an interactive viewpoint navigation facility (mouse + wheel support) and an R programming interface.

As RGL is a real time 3D engine, I expect that using RGL for 2D will give you a fast display.

Please note that this is an old project so I am not sure that it fit your requirement.

You can take a look on this paper and see some result images in this gallery

tripleee
  • 175,061
  • 34
  • 275
  • 318
A. STEFANI
  • 6,707
  • 1
  • 23
  • 48
  • I tried the RGL library, but it's only for 3D OpenGL, not 2D, which I would need here... So I'm afraid this will not work (I tried that option before)... – Tom Wenseleers Apr 17 '18 at 12:26
  • That rdyncall package I think would be the way to go, but I'm afraid I would need some help to use it, as I'm not an OpenGL expert... – Tom Wenseleers Apr 17 '18 at 12:27
  • Do you try to plot data in 3D on a (XY) plan/grid and fix the camera (position and direction) on the Z axe toward the XY grid ? – A. STEFANI Apr 17 '18 at 12:38
  • 1
    Yes that's what I tried, but my data is not 3D so this is not efficient - there are very efficient OpenGL methods for pure 2D graphics which should be used here... – Tom Wenseleers Apr 17 '18 at 12:54
  • " there are very efficient OpenGL methods for pure 2D graphics which should be used here" no there are not. OpenGL does not have this distinction. Drawing is _always_ 2D, because framebuffers are 2D. The 3D transformation stuff is in the shaders and as such totally optional. There are a few bits and pieces in the render pipeline which are meant for 3D data - e.g. the clipping and one could argue also the z-test, but their existence doesn't make 2D rendering less efficient. – derhass Apr 17 '18 at 19:39
  • 1
    @derhass Well at least the small subset of OpenGL that's used in the RGL package wasn't sufficient for what I wanted to do, but feel free to prove me wrong... Using SDL and OpenGL directly using rdyncall on the other hand should work (looking at the demo of the random field linked above), but I haven't figured out how to get it to work for my specific problem... – Tom Wenseleers Apr 18 '18 at 08:05
-1

Five years after first having posted this question finally found the answer:

A faster option than image to display a raster image in R is to use the nativeRaster format of the nara package plus grid.raster. Faster still though is to use an SDL+OpenGl solution which can be called via the rdyncall foreign function interface - that solution is ca. 10x faster still than by using nativeRaster in combination with grid.raster. Theoretically, this should allow for 200+ fps in 640 x 480 resolution, fullscreen or in a window. I am using an OpenGL Texture mapping approach - I haven't tried using an SDL_BlitSurface approach instead - this will give different timings still. And doing the intensity to colour mapping using SDL or OpenGL presumably would also be possible, though I couldn't immediately figure out how to do that (as I now do this beforehand, this would still be a big bottleneck). If anyone knows how to achieve that in pure OpenGL (maybe using glPixelMapfv as mentioned in the question above), please let me know.

First we install the requird packages & calculate a nice example Mandelbrot fractal image :

# 0. load required packages ####
library(remotes)
remotes::install_github("hongyuanjia/rdyncall")
library(rdyncall)
# install SDL libraries (SDL, SDL_image & SDL_mixer, version 1.2) 
# from https://libsdl.org/release/, 
# https://www.libsdl.org/projects/SDL_image/release/ and 
# https://www.libsdl.org/projects/SDL_mixer/release/)
# 64 bit DLLs are to be put under R-4.2.3/bin/x64
# you can use
# source("https://raw.githubusercontent.com/Jean-Romain/lidRviewer/master/sdl.R") 
# on Ubuntu install using sudo apt-get install libsdl1.2-dev libsdl-image1.2-dev libsdl-mixer1.2

remotes::install_github("coolbutuseless/nara")
library(nara)
remotes::install_github("tomwenseleers/mandelExplorer") # for nice example images
library(mandelExplorer)
library(microbenchmark)

# 1. create nice 640x480 image - here Mandelbrot fractal ####
xlims=c(-0.74877,-0.74872)
ylims=c(0.065053,0.065103)
x_res=640L
y_res=480L
nb_iter=as.double(nrofiterations(xlims))
m <- matrix(mandelRcpp2(as.double(xlims[[1]]), as.double(xlims[[2]]), # openmp+SIMD version
                        as.double(ylims[[1]]), as.double(ylims[[2]]), 
                        x_res, 
                        y_res, 
                        nb_iter), 
            nrow=as.integer(x_res))
m[m==0] <- nb_iter
m <- equalizeman(m, nb_iter, rng = c(0, 0.95), levels = 1E4)^(1/8) # equalize colours & apply gamma correction
dim(m) # x_res x y_res, grayscale matrix normalized between 0 and 1
range(m) # 0 1

A. Timings of display of raster image using default image & dbcairo graphics window :

# A. display raster using image - not terribly fast (not fast enough for real-time rendering at good framerate) ####
x11(type = 'dbcairo', antialias = 'none', width = 6*x_res/y_res, height = 6) # Setup a fast graphics device that can render quickly
par(mar=c(0, 0, 0, 0))
microbenchmark(image(m, col=palettes[[2]], asp=y_res/x_res, axes=FALSE, useRaster=TRUE), unit='s') 
# 0.08s, this includes colour mapping+raster display

B. Timings of second solution to display raster images using the nativeRaster functionality of the nara package & grid.raster : 2.3x faster than image

# B. display raster user nara package, using nara::nativeRaster & grid.raster - 2.3x faster than image ####

microbenchmark({ col = palettes[[2]]
              idx = findInterval(m, seq(0, 1, length.out = length(col)))
              colors = col[idx] # matrix of hex colour values
              natrast = nara::raster_to_nr(matrix(t(colors), ncol = x_res,  # nativeRaster 
                                           nrow = y_res, byrow = TRUE))
              grid::grid.raster(natrast, interpolate = FALSE) }, unit='s') 
# 0.035s, 2.3x faster than image(), including the colour mapping + raster display

# timings just for colour mapping part & converting to nativeRaster
microbenchmark({ col = palettes[[2]]
idx = findInterval(m, seq(0, 1, length.out = length(col)))
colors = col[idx] # matrix of hex colour values
natrast = nara::raster_to_nr(matrix(t(colors), ncol = x_res,  # nativeRaster 
                                    nrow = y_res, byrow = TRUE)) }, unit='s')
# 0.027s for colour mapping part & converting to nativeRaster

# timings just for display of nativeRaster
microbenchmark({ grid::grid.raster(natrast, interpolate = FALSE); 
                 dev.flush() }, unit='s')
# 0.004s just for displaying nativeRaster

C. Timings of display of raster images using third solution via SDL & OpenGL using rdyncall package: to display nativeRaster 10x faster than grid.graphics :

# C. display nativeRaster or rgb array using SDL & OpenGL and rdyncall: to display nativeRaster 10x faster than grid.graphics ####

dynport(SDL)
dynport(GL)
dynport(GLU)

# initialize SDL and create a window
if(SDL_Init(SDL_INIT_VIDEO) != 0) {
  stop("Failed to initialize SDL")
}

# set OpenGL attributes
SDL_GL_SetAttribute(SDL_GL_RED_SIZE, 8)
SDL_GL_SetAttribute(SDL_GL_GREEN_SIZE, 8)
SDL_GL_SetAttribute(SDL_GL_BLUE_SIZE, 8)
SDL_GL_SetAttribute(SDL_GL_ALPHA_SIZE, 8)
SDL_GL_SetAttribute(SDL_GL_DOUBLEBUFFER, 1)

# set up screen: graphics window or run full screen
screen <- SDL_SetVideoMode(x_res, y_res, 0, SDL_OPENGL+SDL_DOUBLEBUF) # for windowed mode
# screen <- SDL_SetVideoMode(x_res, y_res, 0, SDL_OPENGL+SDL_DOUBLEBUF+SDL_FULLSCREEN) # for fullscreen mode

if(is.null(screen)) {
  stop("Failed to set video mode")
}

# initialize projection matrix
glMatrixMode(GL_PROJECTION)
glLoadIdentity()
glOrtho(0, x_res, y_res, 0, -1, 1)

# initialize modelview matrix
glMatrixMode(GL_MODELVIEW)
glLoadIdentity()

# generate a texture
texId <- as.integer(0)
glGenTextures(1, texId)
glBindTexture(GL_TEXTURE_2D, texId)

# set texture parameters
glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST)
glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST)

# convert colour vector to matrix
colorsMatrix <- matrix(colors,  
                       nrow = y_res,
                       ncol = x_res, 
                       byrow = TRUE)

# create image matrix (created as nativeRaster to make sure it is a contiguous block of memory)  
imageMatrix <- nara::raster_to_nr(colorsMatrix)

# timings just for nativeRaster display

microbenchmark({ 

# upload the image data to the texture
glTexImage2D(GL_TEXTURE_2D, 0, 4, 
             x_res, y_res, 0, GL_RGBA, GL_UNSIGNED_BYTE, imageMatrix)

# clear the window
glClear(GL_COLOR_BUFFER_BIT)

# enable textures
glEnable(GL_TEXTURE_2D)

# draw a quad with the texture
glBegin(GL_QUADS)
glTexCoord2f(0, 1); glVertex2f(0, 0)
glTexCoord2f(1, 1); glVertex2f(x_res, 0)
glTexCoord2f(1, 0); glVertex2f(x_res, y_res)
glTexCoord2f(0, 0); glVertex2f(0, y_res)
glEnd()

# disable textures
glDisable(GL_TEXTURE_2D)

# swap buffers to display the result
SDL_GL_SwapBuffers() }, 
unit='s') 
# 0.0004s - 10x faster than grid.raster

# close window
SDL_Quit()

Output images are shown here.

(Full disclosure: the OpenGL+SDL+rdyncall code was originally based on ChatGPT4 output; the benchmarks and the rest of the answer were manually added; the code was manually modified though to prevent this answer from being taken down by moderators)

Tom Wenseleers
  • 7,535
  • 7
  • 63
  • 103