read-write support for neuroglancer meshes

R/neuron-io.R

@@ -68,17 +68,21 @@ read.neuron<-function(f, format=NULL, class=c("neuron", "ngraph"), ...){
     writeBin(httr::content(filecontents,type = 'raw'), con = f)
   } else url=NULL
   #if(!file.exists(f)) stop("Unable to read file: ",f)
-  if(is.null(format))
-    format=tolower(sub(".*\\.([^.]+$)","\\1",basename(f)))
-  if(format=="rds")
+  # use ext as backup if format is missing
+  ext=if(is.null(format)) 
+    tolower(sub(".*\\.([^.]+$)","\\1",basename(f))) else ""
+  if(!is.null(format) && format %in% c("ply", "obj", "ngmesh")) 
+    format=paste0("neuron.", format)
+  if(isTRUE(format=="rds") || ext=='rds')
     n=readRDS(f)
-  else if(format=="rda"){
-    objname=load(f,envir=environment())
+  else if(isTRUE(format=="rda") || ext=='rda'){
+    objname=load(f, envir=environment())
     if(length(objname)>1) stop("More than 1 object in file:",f)
     n=get(objname,envir=environment())
   } else {
     class=match.arg(class, choices = c("neuron", "ngraph"))
-    ffs=getformatreader(f, class = class)
+    # note that
+    ffs=getformatreader(f, format=format, class = class)
     if(is.null(ffs)) {
       # as a special test, check to see if this looks like an swc file
       # we don't do this by default since is.swc is a little slow
@@ -181,6 +185,8 @@ read.neurons<-function(paths, pattern=NULL, neuronnames=NULL, format=NULL,
       on.exit(unlink(paths))
       download.file(url, destfile = paths)
     }
+    # the neurons inside the zip file will not have format zip
+    if(isTRUE(format=='zip')) format=NULL
     neurons_dir <- file.path(tempfile(pattern = "user_neurons"))
     on.exit(unlink(neurons_dir, recursive=TRUE), add = TRUE)
     unzip(paths, exdir=neurons_dir)
@@ -423,7 +429,19 @@ registerformat<-function(format=NULL,ext=format,read=NULL,write=NULL,magic=NULL,
 #' write.neuron(Cell07PNs[[1]], swc)
 #' stopifnot(isTRUE(getformatreader(swc)$format=='swc'))
 #' unlink(swc)
-getformatreader<-function(file, class=NULL){
+getformatreader<-function(file, format=NULL, class=NULL) {
+  if(!is.null(format)) {
+    # the format should exactly specify what we need
+    ff=fileformats(format = format, class=class, read = TRUE, rval = 'all')
+    if(length(ff)>1) {
+      if(is.null(class))
+        stop("Please specificy a class to fully identify the reader for format:", format)
+      warning("there are multiple readers for format:", format)
+      ff=ff[length(ff)]
+    }
+    return(ff[[1]])
+  }
+
   formatsforclass<-fileformats(class=class, read = TRUE)
   if(!length(formatsforclass)) return(NULL)
 
@@ -639,10 +657,10 @@ is.swc<-function(f, TrustSuffix=TRUE) {
 #' @param dir Path to directory (this will replace \code{dirname(file)} if
 #'   specified)
 #' @param format Unique abbreviation of one of the registered file formats for
-#'   neurons including 'swc', 'hxlineset', 'hxskel' (skeletons) and 'ply', 'obj'
-#'   (neuron meshes). If no format can be extracted from the filename or the
-#'   \code{ext} parameter, then it defaults to 'swc' for skeletons and 'ply' for
-#'   meshes.
+#'   neurons including 'swc', 'hxlineset', 'hxskel' (skeletons) and 'ply',
+#'   'obj', 'ngmesh' (neuron meshes). If no format can be extracted from the
+#'   filename or the \code{ext} parameter, then it defaults to 'swc' for
+#'   skeletons and 'ply' for meshes.
 #' @param ext Will replace the default extension for the filetype and should
 #'   include the period e.g. \code{ext='.amiramesh'} or \code{ext='_reg.swc'}.
 #'   The special value of ext=NA will prevent the extension from being changed
@@ -677,6 +695,13 @@ is.swc<-function(f, TrustSuffix=TRUE) {
 #' write.neuron(MBL.surf, file = 'MBL.surf.ply')
 #' # ... or Wavefront obj format
 #' write.neuron(MBL.surf, file = 'MBL.surf.obj')
+#' # ... or Neuroglancer precomputed (legacy) mesh format
+#' write.neuron(MBL.surf, file = 'MBL.surf.ngmesh')
+#' # ... same if you really don't want to use (or add) the default extension
+#' write.neuron(MBL.surf, file = 'MBL', format='ngmesh', ext=NA)
+#' # but you'll need to specify the extension
+#' read.neuron('MBL', format='ngmesh')
+#' 
 #' # specify the format directly if not evident from file suffix
 #' # not recommended though as will probably just cause trouble when reading
 #' write.neuron(MBL.surf, file = 'MBL.surf', format='obj')
@@ -707,8 +732,7 @@ write.neuron<-function(n, file=NULL, dir=NULL, format=NULL, ext=NULL,
   if(!is.null(format)) {
     # TODO - one day it should be possible to have one file format associated 
     # with different R classes
-    if(format=='obj') format='neuron.obj'
-    else if(format=='ply') format='neuron.ply'
+    if(format%in% c('obj', 'ply', 'ngmesh')) format=paste0('neuron.', format)
   }
   fw=try(getformatwriter(format=format, file=file, ext=ext, class='neuron'), silent = T)
   if(inherits(fw, 'try-error')) {

R/neuron-mesh.R

@@ -56,5 +56,67 @@ write.neuron.mesh <- function(x, file, format=c("ply", "obj"), ...) {
   else if(format=="obj")
     Rvcg::vcgObjWrite(x, filename=file, writeNormals=F, ...)
   else stop("Unknown format")
+}
+
+read.neuron.ngmesh <- function(x, format=c('mesh3d', "raw"), ...) {
+  format=match.arg(format)
+  bytes <- if(inherits(x, 'response')) {
+    httr::content(x, as = 'raw')
+  } else if(is.list(x) && !is.null(x$content) && is.raw(x$content)) {
+    x$content
+  } else if(is.character(x)) {
+    if(grepl("^http[s]{0,1}://", x)) {
+    res=httr::GET(x, ...)
+    httr::stop_for_status(res)
+    httr::content(res, as = 'raw')
+    } else readBin(x, what=raw(), n=file.size(x))
+  } else stop("Invalid input. I can accept an httr or curl response, a file or url!")
+
+  decode_neuroglancer_mesh(bytes, format = format)
+}
+
+decode_neuroglancer_mesh <- function(bytes, format=c('mesh3d', "raw")) {
+  format=match.arg(format)
+  con=rawConnection(bytes)
+  on.exit(close(con))
+
+  nverts=readBin(con, what = 'int', size = 4, n=1)
+  verts=readBin(con, what='numeric', n=nverts*3, size=4)
+  nidxs=length(bytes)/4-1L-length(verts)
+  idx=readBin(con, what='int', n=nidxs, size=4)
+
+  if(format=='raw') {
+    structure(list(
+      v = matrix(verts, ncol = 3, byrow = T),
+      i = matrix(idx, ncol = 3, byrow = T)
+    ),
+    class = 'ngmesh')
+  } else{
+    rgl::tmesh3d(
+      matrix(verts, nrow=3, byrow = F),
+      matrix(idx+1L, nrow=3, byrow = F),
+      homogeneous = F)
+  }
+}
+
+write.neuron.ngmesh <- function(x, file, ...) {
+
+  if(!inherits(x, 'mesh3d')) {
+    x=tryCatch(as.mesh3d(x), error=function(e) stop("Unable to convert x to mesh3d object! Only neuron meshes can be written in ",format," format!"))
+  }
 
+  nverts=as.integer(ncol(x$vb))
+  stopifnot(isTRUE(nverts>0))
+  verts <- c(x$vb[1:3, , drop = FALSE])
+  mode(verts)='double'
+  idxs=c(x$it-1L)
+  mode(idxs)='integer'
+  if(is.character(file)) {
+    file=file(file, open = 'wb')
+    on.exit(close(file))
+  }
+  writeBin(nverts, file, size = 4, endian='little')
+  writeBin(verts, file, size=4, endian='little')
+  writeBin(idxs, file, size=4, endian='little')
 }
+

R/zzz.R

@@ -30,7 +30,9 @@ update_igraph <- FALSE
                  read=read.neuron.mesh, write=write.neuron.obj, class='neuron')
   registerformat('neuron.ply', ext='.ply', magic=is.ply, magiclen = 3,
                  read=read.neuron.mesh, write=write.neuron.ply, class='neuron')
-
+  registerformat('neuron.ngmesh', ext='.ngmesh', 
+                 read=read.neuron.ngmesh, write=write.neuron.ngmesh, class='neuron')
+
   # image formats
   registerformat('nrrd', ext=c('.nrrd','.nhdr'), read=read.im3d.nrrd, 
                  write=write.nrrd, magic=is.nrrd,

tests/testthat/test-neuron-io.R

@@ -201,14 +201,14 @@ test_that("we can read hxskel format neurons",{
 })
 
 test_that("we can read multiple neurons from a zip archive", {
-  files_to_zip <- c("testdata/neuron/testneuron_am3d.am", "testdata/neuron/testneuron_lineset.am")
+  files_to_zip <- test_path(c("testdata/neuron/testneuron_am3d.am", "testdata/neuron/testneuron_lineset.am"))
   # swallow extraneous warning
   expect_warning(neurons <- read.neurons(files_to_zip,
                                          neuronnames = function(f) tools::file_path_sans_ext(basename(f))),
                  regexp = "specifies radius")
   zip_file <- paste0(tempfile(), ".zip")
   on.exit(unlink(zip_file, recursive=TRUE))
-  zip(zip_file, files_to_zip)
+  zip(zip_file, files_to_zip, flags = "-r9Xq")
   expect_warning(zip_neurons <- read.neurons(zip_file, format="zip",
                                              neuronnames = function(f) tools::file_path_sans_ext(basename(f))),
                  regexp = "specifies radius")

tests/testthat/test-neuron-mesh.R

@@ -1,12 +1,30 @@
-skip_if_not_installed('Rvcg')
-
 bl=neuronlist(icosahedron3d(), tetrahedron3d())
 names(bl)=c("a","b")
-td=tempfile()
-setup(dir.create(td))
-teardown(unlink(td, recursive = TRUE))
+
+test_that("read/write works with ngmesh files", {
+  td3=tempfile()
+  setup(dir.create(td3))
+  teardown(unlink(td3, recursive = TRUE))
+
+  expect_silent(ff1 <- write.neurons(bl, dir=td3, format='ngmesh'))
+  expect_is(bl2 <- read.neurons(td3, format = 'ngmesh'), 'neuronlist')
+  expect_equal(summary(bl), summary(bl2))
+  expect_error(write.neurons(Cell07PNs[1:3], format = 'ngmesh'))
+
+  expect_silent(write.neuron(MBL.surf, file = file.path(td3, 'MBL'), ext=NA, format = 'ngmesh'))
+  expect_equal(tolerance = 1e-6,
+    read.neuron(file.path(td3, 'MBL'), format='ngmesh'),
+    as.mesh3d(MBL.surf))
+})
+
 
 test_that("read/write works", {
+  skip_if_not_installed('Rvcg')
+
+  td=tempfile()
+  setup(dir.create(td))
+  teardown(unlink(td, recursive = TRUE))
+
   expect_silent(ff1 <- write.neurons(bl, dir=td, format='ply'))
   md5.1=tools::md5sum(ff1)
   expect_warning(ff2 <- write.neurons(bl, dir=td, Force=TRUE), regexp = 'ply')
@@ -28,13 +46,14 @@ test_that("read/write works", {
   expect_equal(read.neuron(file.path(td, "MBL.surf.ply")), as.mesh3d(MBL.surf), tolerance = 1e-6)
 })
 
-skip_if_not_installed('readobj')
-
-td2=tempfile()
-setup(dir.create(td2))
-teardown(unlink(td2, recursive = TRUE))
 
 test_that("read/write works with obj files", {
+  skip_if_not_installed('readobj')
+
+  td2=tempfile()
+  setup(dir.create(td2))
+  teardown(unlink(td2, recursive = TRUE))
+
   expect_silent(ff1 <- write.neurons(bl, dir=td2, format='obj'))
   expect_is(bl2 <- read.neurons(td2), 'neuronlist')
   expect_equal(summary(bl), summary(bl2))