seung-lab/zmesh

Marching Cubes & Mesh Simplification on multi-label 3D images.

zmesh: Multi-Label Marching Cubes & Mesh Simplification

from zmesh import Mesher

labels = ... # some dense volumetric labeled image
mesher = Mesher( (4,4,40) ) # anisotropy of image

# initial marching cubes pass
# close controls whether meshes touching
# the image boundary are left open or closed
mesher.mesh(labels, close=False) 

meshes = []
for obj_id in mesher.ids():
  meshes.append(
    mesher.get(
      obj_id, 
      normals=False, # whether to calculate normals or not

      # tries to reduce triangles by this factor
      # 0 disables simplification
      reduction_factor=100, 

      # Max tolerable error in physical distance
      # note: if max_error is not set, the max error
      # will be set equivalent to one voxel along the 
      # smallest dimension.
      max_error=8,
      # whether meshes should be centered in the voxel
      # on (0,0,0) [False] or (0.5,0.5,0.5) [True]
      voxel_centered=False, 
    )
  )
  mesher.erase(obj_id) # delete high res mesh

mesher.clear() # clear memory retained by mesher

mesh = meshes[0]
mesh = mesher.simplify(
  mesh, 
  # same as reduction_factor in get
  reduction_factor=100, 
  # same as max_error in get
  max_error=40, 
  compute_normals=False, # whether to also compute face normals
) # apply simplifier to a pre-existing mesh

# compute normals without simplifying
mesh = mesher.compute_normals(mesh) 

mesh.vertices
mesh.faces 
mesh.normals
mesh.triangles() # compute triangles from vertices and faces

# Extremely common obj format
with open('iconic_doge.obj', 'wb') as f:
  f.write(mesh.to_obj())

# Common binary format
with open('iconic_doge.ply', 'wb') as f:
  f.write(mesh.to_ply())

# Neuroglancer Precomputed format
with open('10001001:0', 'wb') as f:
  f.write(mesh.to_precomputed())

Note: As of the latest version, mesher.get_mesh has been deprecated in favor of mesher.get which fixes a long standing bug where you needed to transpose your data in order to get a mesh in the correct orientation.

Installation

If binaries are not available for your system, ensure you have a C++ compiler installed.

pip install zmesh

Performance Tuning & Notes

• The mesher will consume about double memory in 64 bit mode if the size of the object exceeds <1023, 1023, 511> on the x, y, or z axes. This is due to a limitation of the 32-bit format.
• The mesher is ambidextrous, it can handle C or Fortran order arrays.
• The maximum vertex range supported .simplify after converting to voxel space is 2²⁰ (appx. 1M) due to the packed 64-bit vertex format.

Related Projects

• zi_lib - zmesh makes heavy use of Aleks' C++ library.
• Igneous - Visualization of connectomics data using cloud computing.

Credits

Thanks to Aleks Zlateski for creating and sharing this beautiful mesher.

Later changes by Will Silversmith, Nico Kemnitz, and Jingpeng Wu.

References

1. W. Lorensen and H. Cline. "Marching Cubes: A High Resolution 3D Surface Construction Algorithm". pp 163-169. Computer Graphics, Volume 21, Number 4, July 1987. (link)
2. M. Garland and P. Heckbert. "Surface simplification using quadric error metrics". SIGGRAPH '97: Proceedings of the 24th annual conference on Computer graphics and interactive techniques. Pages 209–216. August 1997. doi: 10.1145/258734.258849 (link)
3. H. Hoppe. "New Quadric Metric for Simplifying Meshes with Appearance Attributes". IEEE Visualization 1999 Conference. pp. 59-66. doi: 10.1109/VISUAL.1999.809869 (link)

zmesh's Issues

Simplify Method for Existing Mesh

Not too difficult, but requires writing a conversion for vertex/faces to triangles.

Problems simplifying meshes

Hi Will,

I somewhat randomly stumbled on your library in search of something to simplify a flywire mesh. The example code seems pretty straight forward yet I don't get the expected results:

>>> from cloudvolume import CloudVolume
>>> vol = CloudVolume('graphene://[...]/fly_v26')
>>> mesh = vol.mesh.get(720575940610453042)[720575940610453042]
>>> # Very large mesh -> hard to work with
>>> mesh
Mesh(vertices<990112>, faces<1887184>, normals<0>, segid=720575940610453042, encoding_type=<draco>)

>>> from zmesh import Mesher
>>> # Assuming isotropy since mesh has coordinates in nm
>>> mesher = Mesher( (1,1,1) )
>>> # Tried ramping up max error
>>> simp = mesher.simplify(mesh, 100, max_error=4e6, compute_normals=False)
>>> # Lost some vertices but no faces
>>> simp
Mesh(vertices<938272>, faces<1887184>, normals<None>)

Any pointers?

Thanks a bunch! :)

Add mesh data equality tests

Lots of people are depending on zmesh now, so we should make sure updates don't change the output unintentionally.

Remove zi_lib Dependency on Boost

I think most of the time Boost is needed it was to supplement future additions to the standard library (when ZI_NO_BOOST is specified, tr1 is used, but tr1 is also too old). Most of this can be done via the standard library now since c++11.

Simplification is independent of reduction_factor

Regardless of the reduction factor, mesher.simplify seems to simplify what seems like as much as possible:

import zmesh
import numpy as np

for order in ['C','F']:
    labels = np.zeros((11, 17, 19), dtype=np.uint8, order=order)
    labels[1:-1, 1:-1, 1:-1] = 1

    mesher = zmesh.Mesher((4, 4, 40))
    mesher.mesh(labels)

    original_mesh = mesher.get_mesh(1, normals=False)
    mesh = mesher.simplify(
        original_mesh, reduction_factor=2, max_error=40, compute_normals=True
    )
    print(f"Actual reduction factor for order {order}: {len(original_mesh.faces)/len(mesh.faces)}")

Actual reduction factor for order C: 24.636363636363637
Actual reduction factor for order F: 20.452830188679247

The actual reduction seems to be slightly affected by whether it is order C or F.

Increase 32-bit X Representation Space to 1023 Voxels

For some reason, the last bit of the X representation in 32-bit mode is overflowing in marching_cubes.hpp. This prevents us from using all 11 bits, so we are limited to 511 voxels in the X dimension. There's a few ways this could happen:

• Two vertices interact using addition and it triggers a carry / overflow.
• Mysterious conversion from uint32 to int32.

How to combine surfaces from chunks (or how to deal with a real large region)

I am trying to generate the surface mesh of a very larger object (the bounding box is around 8000x8000x1400 voxels), which requires a very large memory. I am wondering if there is a way to combine the meshes generated from each chunck? (like 800x800x140 voxels) Thanks

User Friendly Mesh Output

Igneous does the following to make user friendly meshes. It would be good to integrate some of this logic. The output of zmesh.get_mesh is fixed point premultiplied by 2, which is not obvious...

  def _create_mesh(self, obj_id):
    mesh = self._mesher.get_mesh(
      obj_id,
      simplification_factor=self.options['simplification_factor'],
      max_simplification_error=self.options['max_simplification_error']
    )

    self._mesher.erase(obj_id)

    vertices = self._update_vertices(
      np.array(mesh['points'], dtype=np.float32)
    )
    vertex_index_format = [
      np.uint32(len(vertices) / 3), # Number of vertices (3 coordinates)
      vertices,
      np.array(mesh['faces'], dtype=np.uint32)
    ]
    return b''.join([array.tobytes() for array in vertex_index_format])

  def _update_vertices(self, points):
    # zi_lib meshing multiplies vertices by 2.0 to avoid working with floats,
    # but we need to recover the exact position for display
    # Note: points are already multiplied by resolution, but missing the offset
    points /= 2.0
    resolution = self._volume.resolution
    xmin, ymin, zmin = self._bounds.minpt - self.options['low_padding']
    points[0::3] = points[0::3] + xmin * resolution.x
    points[1::3] = points[1::3] + ymin * resolution.y
    points[2::3] = points[2::3] + zmin * resolution.z
    return points

Build manylinux Docker with Boost

Why do we need to transpose?

In Igneous, we were writing (for a Fortran array):

self._mesher.mesh(data[..., 0].T)

Why is that? How can we avoid that?

how to export step from obj?

how to export step from obj?

Is it possible to export IGES,STEP using zmesh package ?

Is it possible to export IGES,STEP using zmesh package ?

Reading Out Meshes is Slow

Marching cubes is super fast, but very dumbly, just reading the meshes out is slow. We should fix this.

Questions on obj file format conversion and nrrd file format support

Hi Will,

More related to generating mesh in precomputed file format workflow Questions on understanding the workflow for generating mesh "precomputed" data #406, does zmesh support converting a collection of obj file format to one or more precomputed files? Basically looking to do the other way around of the following code. Also, is nrrd file format supported to generate mesh?

# Extremely common obj format
with open('iconic_doge.obj', 'wb') as f:
  f,write(mesh.to_obj())

# Common binary format
with open('iconic_doge.ply', 'wb') as f:
  f,write(mesh.to_ply())

# Neuroglancer Precomputed format
with open('10001001:0', 'wb') as f:
  f.write(mesh.to_precomputed())

Thank you,
-m

Create Official First Release

Transparently Handle C vs Fortran Order

OSX wheels for Python 3.7

Hi there! I was wondering if there was any chance for wheels for OSX + Python 3.7? :)

Neuroglancer precomputed format

Hello,

how should one use that precomputed mesh for Neuroglancer? When starting Neuroglancer up, it throws an error:

"Error retrieving chunk 13955040: Error: Error parsing "fragments" property: Expected array, received: ".../meshes/13955040"."

'meshes' is the directory of the mesh data and we are trying to get the mesh of id 13955040.

As per Neuroglancer documentation, the ":0" should be a JSON metadata file, and not the mesh bytes file. Did zmesh bypass this in any way?

Any response is much appreciated.

Best regards,
Andrei Mancu