How to display MeshRegion without vertices

Question

I'm writing a library for handling very large and complicated (boundary) mesh regions on the order of ~150,000 vertices, ~450,000 edges, and ~300,000 triangles. The regions themselves are the cortical surfaces of brains (derived from MRI experiments), and the convention for such surfaces is to plot only the colored triangles with no edges or vertices.

I've noticed that, while it's pretty straightforward to plot the region as a Graphics3D form in which each triangle is rendered as something like Polygon[vertexCoordinates, VertexColors -> colorData], these rather large graphics3d objects don't rotate or render at all quickly. In fact, they're painfully slow when compared to similar routines in other languages, including Matlab.

My suspicion was that a MeshRegion or BoundaryMeshRegion would include some under-the-hood optimization and render/rotate more quickly. In some tests, this seems to be true; in fact, the MeshRegion graphics render and rotate almost instantaneously by comparison with the Graphics3D objects.

My question is this: can I display the MeshRegions without displaying the vertices and lines? I can edit the color of the triangles easily enough via the MeshCellStyle option, but I can't find any way to entirely eliminate the points and lines. Even MeshCellStyle -> {{0,All} -> Opacity[0], {1,All} -> Opacity[0]} does not entirely clear them from the display, and they continue to obscure the visualization of the triangles. Weirdly, however, this seems to work for smaller/simpler meshes.

For example; the following produces two (very slightly) different-looking meshes:

GraphicsRow[
  Map[
    MeshRegion[
      {{0, 0, 0}, {1, 0, 0}, {0, 1, 0}, {0, 0, 1}},
      Polygon[{{1, 2, 3}, {1, 2, 4}, {1, 3, 4}, {2, 3, 4}}],
      MeshCellStyle -> Append[#, {2, All} -> Gray],
      Lighting -> "Neutral"] &,
    {{}, {{0, All} -> Opacity[0], {1, All} -> Opacity[0]}}],
  ImageSize -> 10*72]

However, this code produces a mesh that still contain a lot of blue...

MeshRegion[
  pialMeshLH,
  MeshCellStyle -> {
    {0, All} -> Opacity[0],
    {1, All} -> Opacity[0],
    {2, All} -> Gray}]

Even more confusingly, the following code renderes something incomprehensibly gray and blue:

MeshRegion[
  pialMeshLH,
  MeshCellStyle -> {
    {0, All} -> Gray,
    {1, All} -> Gray,
    {2, All} -> Gray}]

Ignoring scale/zoom, here's approximately how I would like this rendering to look (this rendered via Graphics3D):

The brain shown here is not a real person's brain, but is the pial surface of the left hemisphere of the fsaverage brain in the (open source) FreeSurfer software suite. I've created an approximately 9MB binary file that contains just the vertex and face data for this mesh and uploaded it to filedropper (you can download it here: http://www.filedropper.com/fsmesh); the following function will read it into Mathematica, given the filename:

ReadFSMesh[filename_String] := Block[
  {$ByteOrdering = 1},
  With[
    {dat = Import[filename, "Real64"]},
    With[
      {n = Round@dat[[1]], m = Round@dat[[2]]},
      MeshRegion[
        Partition[dat[[3 ;; (3*n + 2)]], 3],
        Polygon[Round@Partition[dat[[(3*n + 3) ;; All]], 3]]]]]];

Ideally, I'd like to know if it's possible to either (A) speed up the rendering and rotation of complex Graphics3D objects or (B) remove the edges and vertices displayed in the MeshRegion so as to speed up their display.

Thanks in advance.

Answer 1

I suppose MeshRegion and Graphics3D have similar methods of the visualization. However, Graphics3D is more convenient and GraphicsComplex can increase the performance:

mesh = ReadFSMesh@"fsmesh.bin";
brain = GraphicsComplex[MeshCoordinates@#, MeshCells[#, 2]] &@mesh;

Graphics3D[{EdgeForm[], Gray, brain}, Lighting -> "Neutral", Boxed -> False]

Also you can increase the quality with VertexNormals. Since you don't know them you can restore them like in my answer here. I assume your surface consists of triangles only with a proper (clockwise or counterclockwise) order of vertices. Vertex normals are calculated as an average normal of neighborhood triangles.

normals[p_List, cells_] := Module[{t0, t1, t2, P, ux, uy, uz, vx, vy, vz},
  {t0, t1, t2} = Transpose@cells[[All, 1]];
  P = SparseArray[Join @@ (Transpose@{#, Range@Length@#} & /@ {t0, t1, t2}) -> 
     ConstantArray[1, 3 Length@t0]];
  {ux, uy, uz} = Transpose[p[[t1]] - p[[t0]]];
  {vx, vy, vz} = Transpose[p[[t2]] - p[[t0]]];
  P.Transpose@{uy vz - uz vy, uz vx - ux vz, ux vy - uy vx}]

brain2 = GraphicsComplex[MeshCoordinates@#, MeshCells[#, 2], 
     VertexNormals -> normals[MeshCoordinates@#, MeshCells[#, 2]]] &@mesh;

Graphics3D[{EdgeForm[], Gray, brain2}, Lighting -> "Neutral", Boxed -> False]