add support for writing surface meshes in VTK format

NAMESPACE

@@ -58,6 +58,7 @@ export(write.fs.morph)
 export(write.fs.patch)
 export(write.fs.surface)
 export(write.fs.surface.asc)
+export(write.fs.surface.vtk)
 export(write.fs.weight)
 importFrom(grDevices,rgb)
 importFrom(pkgfilecache,ensure_files_available)

R/write_fs_surface.R

@@ -129,3 +129,56 @@ write.fs.surface.asc <- function(filepath, vertex_coords, faces) {
 
   return(invisible('tris'));
 }
+
+
+#' @title Write mesh to file in VTK ASCII format
+#'
+#' @param filepath string. Full path to the output surface file, should end with '.vtk', but that is not enforced.
+#'
+#' @param vertex_coords n x 3 matrix of doubles. Each row defined the x,y,z coords for a vertex.
+#'
+#' @param faces n x 3 matrix of integers. Each row defined the 3 vertex indices that make up the face. WARNING: Vertex indices should be given in R-style, i.e., the index of the first vertex is 1. However, they will be written in FreeSurfer style, i.e., all indices will have 1 substracted, so that the index of the first vertex will be zero.
+#'
+#' @return string the format that was written. One of "tris" or "quads". Currently only triangular meshes are supported, so always 'tris'.
+#'
+#' @family mesh functions
+#'
+#' @examples
+#' \donttest{
+#'     # Read a surface from a file:
+#'     surface_file = system.file("extdata", "lh.tinysurface",
+#'      package = "freesurferformats", mustWork = TRUE);
+#'     mesh = read.fs.surface(surface_file);
+#'
+#'     # Now save it:
+#'     write.fs.surface.vtk(tempfile(fileext=".vtk"), mesh$vertices, mesh$faces);
+#' }
+#'
+#' @export
+write.fs.surface.vtk <- function(filepath, vertex_coords, faces) {
+
+  fh = file(filepath, "w");
+
+  num_verts = nrow(vertex_coords);
+  num_faces = nrow(faces);
+
+  # write header
+  writeLines(c("# vtk DataFile Version 1.0", "fsbrain output", "ASCII", "DATASET POLYDATA", sprintf("POINTS %d float", num_verts)), fh);
+  close(fh);
+
+
+  # Append the vertex data
+  write.table(vertex_coords, file = filepath, append = TRUE, quote = FALSE, sep = " ", row.names = FALSE, col.names = FALSE);
+
+  fh = file(filepath, "a");
+  writeLines(c(sprintf("POLYGONS %d %d", num_faces, num_faces * 4L)), fh);
+  close(fh);
+
+  # Append the face data
+
+  faces = faces - 1L;       # from R to 0-based indices
+  faces = cbind(3L, faces); # in VTK format, each face line starts with the number of vertices
+  write.table(faces, file = filepath, append = TRUE, quote = FALSE, sep = " ", row.names = FALSE, col.names = FALSE);
+
+  return(invisible('tris'));
+}