diff --git a/examples/mesh/README.rst b/examples/mesh/README.rst
new file mode 100644
index 000000000..99e569fed
--- /dev/null
+++ b/examples/mesh/README.rst
@@ -0,0 +1,2 @@
+Mesh Examples
+=============
diff --git a/examples/mesh/image_surface.py b/examples/mesh/image_surface.py
index c089f0f90..fce3c4958 100644
--- a/examples/mesh/image_surface.py
+++ b/examples/mesh/image_surface.py
@@ -10,19 +10,18 @@
 
 import imageio.v3 as iio
 import fastplotlib as fpl
-import numpy as np
 import scipy.ndimage
 
 im = iio.imread("imageio:astronaut.png")
 
-figure = fpl.Figure(size=(700, 560), cameras='3d', controller_types='orbit')
+figure = fpl.Figure(size=(700, 560), cameras="3d", controller_types="orbit")
 
 
 # Create the height map from the image
 z = im.mean(axis=2)
 z = scipy.ndimage.gaussian_filter(z, 5)  # 2nd arg is sigma
 
-mesh = figure[0, 0].add_surface(z, colors="magenta", cmap=im)
+mesh = figure[0, 0].add_surface(z, cmap=im)
 mesh.world_object.local.scale_y = -1
 
 
diff --git a/examples/mesh/mesh.py b/examples/mesh/mesh.py
index b22ae61db..4c8de088d 100644
--- a/examples/mesh/mesh.py
+++ b/examples/mesh/mesh.py
@@ -9,11 +9,10 @@
 # sphinx_gallery_pygfx_docs = 'screenshot'
 
 import fastplotlib as fpl
-import numpy as np
 import pygfx as gfx
 
 
-figure = fpl.Figure(size=(700, 560), cameras='3d', controller_types='orbit')
+figure = fpl.Figure(size=(700, 560), cameras="3d", controller_types="orbit")
 
 
 # Load geometry using Pygfx's geometry util
@@ -21,9 +20,9 @@
 positions = geo.positions.data
 indices = geo.indices.data
 
-mesh = figure[0, 0].add_mesh(positions, indices, colors="magenta")
-
+mesh = fpl.MeshGraphic(positions, indices, colors="magenta")
 
+figure[0, 0].add_graphic(mesh)
 figure[0, 0].axes.grids.xy.visible = True
 figure[0, 0].camera.show_object(mesh.world_object, (1, 1, -1), up=(0, 0, 1))
 
diff --git a/examples/mesh/polygon_animation.py b/examples/mesh/polygon_animation.py
new file mode 100644
index 000000000..6d4bc7bf0
--- /dev/null
+++ b/examples/mesh/polygon_animation.py
@@ -0,0 +1,76 @@
+"""
+Polygon animation
+=================
+
+Polygon animation example that changes the polygon data. Random points are generated by sampling from a
+2D gaussian and a polygon is updated to visualize a convex hull for the sampled points.
+
+"""
+
+# test_example = false
+# sphinx_gallery_pygfx_docs = 'animate 8s'
+
+import numpy as np
+from scipy.spatial import ConvexHull
+import fastplotlib as fpl
+
+
+def points_to_hull(points) -> np.ndarray:
+    hull = ConvexHull(points, qhull_options="Qs")
+    return points[hull.vertices]
+
+
+figure = fpl.Figure(size=(700, 560))
+
+
+cov = np.array([[1, 0], [0, 1]])
+
+# sample points from a 2d gaussian
+samples1 = np.random.multivariate_normal((0, 0), cov, size=20)
+samples2 = np.random.multivariate_normal((5, 0), cov, size=50)
+
+# add the convex hull as a polygon
+polygon1 = figure[0, 0].add_polygon(
+    points_to_hull(samples1), colors="cyan", alpha=0.7, alpha_mode="blend"
+)
+# add the sampled points
+scatter1 = figure[0, 0].add_scatter(
+    samples1, sizes=8, colors="blue", alpha=0.7, alpha_mode="blend"
+)
+
+# add the second gaussian and convex hull polygon
+polygon2 = figure[0, 0].add_polygon(
+    points_to_hull(samples2), colors="magenta", alpha=0.7, alpha_mode="blend"
+)
+scatter2 = figure[0, 0].add_scatter(
+    samples2, sizes=8, colors="r", alpha=0.7, alpha_mode="blend"
+)
+
+
+def animate():
+    # set new scatter data
+    scatter1.data[:, :-1] += np.random.normal(0, 0.05, size=samples1.size).reshape(
+        samples1.shape
+    )
+    # set convex hull with new polygon vertices
+    polygon1.data = points_to_hull(scatter1.data[:, :-1])
+
+    # set the other scatter and polygon
+    scatter2.data[:, :-1] += np.random.normal(0, 0.05, size=samples2.size).reshape(
+        samples2.shape
+    )
+    polygon2.data = points_to_hull(scatter2.data[:, :-1])
+
+
+figure.show()
+figure[0, 0].camera.width = 10
+figure[0, 0].camera.height = 10
+
+figure.add_animations(animate)
+
+
+# NOTE: fpl.loop.run() should not be used for interactive sessions
+# See the "JupyterLab and IPython" section in the user guide
+if __name__ == "__main__":
+    print(__doc__)
+    fpl.loop.run()
diff --git a/examples/mesh/polygons.py b/examples/mesh/polygons.py
new file mode 100644
index 000000000..616c2e0fb
--- /dev/null
+++ b/examples/mesh/polygons.py
@@ -0,0 +1,61 @@
+"""
+Polygons
+========
+
+An example with polygons.
+
+"""
+
+# test_example = True
+# sphinx_gallery_pygfx_docs = 'screenshot'
+
+import fastplotlib as fpl
+import numpy as np
+from cmap import Colormap
+
+figure = fpl.Figure(size=(700, 560))
+
+
+def make_circle(center, radius: float, n_points: int = 75) -> np.ndarray:
+    theta = np.linspace(0, 2 * np.pi, n_points, endpoint=False)
+    xs = radius * np.sin(theta)
+    ys = radius * np.cos(theta)
+
+    return np.column_stack([xs, ys]) + np.asarray(center)[None]
+
+
+# define vertices for some polygons
+circle_data = make_circle(center=(0, 0), radius=5)
+octogon_data = make_circle(center=(15, 0), radius=7, n_points=8)
+rectangle_data = np.array([[10, 10], [20, 10], [20, 15], [10, 15]])
+triangle_data = np.array(
+    [
+        [-5, 8],
+        [5, 8],
+        [0, 15],
+        [-5, 8],
+    ]
+)
+
+# add polygons
+figure[0, 0].add_polygon(circle_data, name="circle")
+figure[0, 0].add_polygon(
+    octogon_data,
+    colors=Colormap("jet").lut(8), # set vertex colors from jet cmap
+    name="octogon"
+)
+figure[0, 0].add_polygon(
+    rectangle_data,
+    colors=["r", "r", "cyan", "y"], # manually specify vertex colors
+    name="rectangle"
+)
+figure[0, 0].add_polygon(triangle_data, colors="m")
+
+figure.show()
+
+
+# NOTE: fpl.loop.run() should not be used for interactive sessions
+# See the "JupyterLab and IPython" section in the user guide
+if __name__ == "__main__":
+    print(__doc__)
+    fpl.loop.run()
diff --git a/examples/mesh/surface_earth.py b/examples/mesh/surface_earth.py
new file mode 100644
index 000000000..c2e137bc8
--- /dev/null
+++ b/examples/mesh/surface_earth.py
@@ -0,0 +1,95 @@
+"""
+Earth sphere animation
+======================
+
+Example showing how to create a sphere with an image of the Earth and rotate it around its 23.44° axis of rotation
+with respect to the ecliptic (the xz plane in the visualization).
+
+"""
+
+# test_example = false
+# sphinx_gallery_pygfx_docs = 'animate 8s'
+
+import fastplotlib as fpl
+import numpy as np
+import imageio.v3 as iio
+import pylinalg as la
+
+
+figure = fpl.Figure(size=(700, 560), cameras="3d", controller_types="orbit")
+
+# create a sphere from spherical coordinates
+# see this for reference: https://mathworld.wolfram.com/SphericalCoordinates.html
+# phi and theta are swapped in this example w.r.t. the wolfram alpha description
+radius = 10
+nx = 101
+phi = np.linspace(0, np.pi * 2, num=nx, dtype=np.float32)
+ny = 51
+theta = np.linspace(0, np.pi, num=ny, dtype=np.float32)
+
+phi_grid, theta_grid = np.meshgrid(phi, theta)
+
+# convert to cartesian coordinates
+theta_grid_sin = np.sin(theta_grid)
+x = radius * np.cos(phi_grid) * theta_grid_sin * -1
+y = radius * np.cos(theta_grid)
+z = radius * np.sin(phi_grid) * theta_grid_sin
+
+# get texture coords to map the image onto the mesh positions
+u = phi_grid / (np.pi * 2)
+v = 1 - (theta_grid / np.pi)
+texcoords = np.dstack([u, v]).reshape(-1, 2)
+
+# get an image of the earth from nasa
+image = iio.imread(
+    "https://svs.gsfc.nasa.gov/vis/a000000/a003600/a003615/flat_earth_Largest_still.0330.jpg"
+)
+# images coordinate systems are typically inverted in y, so flip the image
+image = np.ascontiguousarray(np.flipud(image))
+
+# create a sphere
+sphere = figure[0, 0].add_surface(
+    np.dstack([x, y, z]),
+    mode="phong",
+    colors="magenta",
+    cmap=image,
+    mapcoords=texcoords,
+)
+
+# display xz plane as a grid
+figure[0, 0].axes.grids.xz.visible = True
+figure.show()
+
+# view from top right angle
+figure[0, 0].camera.show_object(sphere.world_object, (-0.5, -0.25, -1), up=(0, 1, 0))
+figure[0, 0].camera.zoom = 1.25
+
+# create quaternion for 23.44 degrees axial tilt
+axial_tilt = la.quat_from_euler((np.radians(23.44), 0), order="XY")
+
+# a line to indicate the axial tilt
+figure[0, 0].add_line(
+    np.array([[0, -20, 0], [0, 20, 0]]), rotation=axial_tilt, colors="magenta"
+)
+
+rot = 1
+
+
+def rotate():
+    # rotate by 1 degree
+    global rot
+    rot += 1
+    rot_quat = la.quat_from_euler((0, np.radians(rot)), order="XY")
+
+    # apply rotation w.r.t. axial tilt
+    sphere.rotation = la.quat_mul(axial_tilt, rot_quat)
+
+
+figure[0, 0].add_animations(rotate)
+
+
+# NOTE: fpl.loop.run() should not be used for interactive sessions
+# See the "JupyterLab and IPython" section in the user guide
+if __name__ == "__main__":
+    print(__doc__)
+    fpl.loop.run()
diff --git a/examples/mesh/surface_ellipsoid.py b/examples/mesh/surface_ellipsoid.py
new file mode 100644
index 000000000..6d7cdae7b
--- /dev/null
+++ b/examples/mesh/surface_ellipsoid.py
@@ -0,0 +1,55 @@
+"""
+Ellipsoid surface
+=================
+
+Simple example of a sphere surface mesh with a colormap indicating z values.
+
+"""
+
+# test_example = false
+# sphinx_gallery_pygfx_docs = 'screenshot'
+
+import fastplotlib as fpl
+import numpy as np
+
+figure = fpl.Figure(size=(700, 560), cameras="3d", controller_types="orbit")
+
+# create an ellipsoid from spherical coordinates
+# see this for reference: https://mathworld.wolfram.com/SphericalCoordinates.html
+# phi and theta are swapped in this example w.r.t. the wolfram alpha description
+radius = 10
+
+nx = 101
+phi = np.linspace(0, np.pi * 2, num=nx, dtype=np.float32)
+ny = 51
+theta = np.linspace(0, np.pi, num=ny, dtype=np.float32)
+
+phi_grid, theta_grid = np.meshgrid(phi, theta)
+
+# convert to cartesian coordinates
+theta_grid_sin = np.sin(theta_grid)
+x = radius * np.cos(phi_grid) * theta_grid_sin * -1
+y = radius * np.cos(theta_grid)
+
+# elongate along z axis
+z = radius * 2 * np.sin(phi_grid) * theta_grid_sin
+
+sphere = figure[0, 0].add_surface(
+    np.dstack([x, y, z]),
+    mode="phong",
+    cmap="bwr",  # by default, providing a colormap name will map the colors to z values
+)
+
+# display xz plane as a grid
+figure[0, 0].axes.grids.xy.visible = True
+figure.show()
+
+# view from top right angle
+figure[0, 0].camera.show_object(sphere.world_object, (1, 1, -1), up=(0, 0, 1))
+
+
+# NOTE: fpl.loop.run() should not be used for interactive sessions
+# See the "JupyterLab and IPython" section in the user guide
+if __name__ == "__main__":
+    print(__doc__)
+    fpl.loop.run()
diff --git a/examples/mesh/surface_gaussian.py b/examples/mesh/surface_gaussian.py
new file mode 100644
index 000000000..6a9fb0f1d
--- /dev/null
+++ b/examples/mesh/surface_gaussian.py
@@ -0,0 +1,45 @@
+"""
+Gaussian kernel as a surface
+============================
+
+Example showing a gaussian kernel as a surface mesh
+"""
+
+# test_example = true
+# sphinx_gallery_pygfx_docs = 'screenshot'
+
+import fastplotlib as fpl
+import numpy as np
+
+
+figure = fpl.Figure(size=(700, 560), cameras="3d", controller_types="orbit")
+
+
+def gaus2d(x=0, y=0, mx=0, my=0, sx=1, sy=1):
+    return (
+        1.0
+        / (2.0 * np.pi * sx * sy)
+        * np.exp(
+            -((x - mx) ** 2.0 / (2.0 * sx**2.0) + (y - my) ** 2.0 / (2.0 * sy**2.0))
+        )
+    )
+
+
+r = np.linspace(0, 10, num=200)
+x, y = np.meshgrid(r, r)
+z = gaus2d(x, y, mx=5, my=5, sx=1, sy=1) * 50
+
+mesh = figure[0, 0].add_surface(
+    np.dstack([x, y, z]), mode="phong", cmap="jet"
+)
+
+# figure[0, 0].axes.grids.xy.visible = True
+figure[0, 0].camera.show_object(mesh.world_object, (-2, 2, -2), up=(0, 0, 1))
+figure.show()
+
+
+# NOTE: fpl.loop.run() should not be used for interactive sessions
+# See the "JupyterLab and IPython" section in the user guide
+if __name__ == "__main__":
+    print(__doc__)
+    fpl.loop.run()
diff --git a/examples/mesh/surface.py b/examples/mesh/surface_height.py
similarity index 70%
rename from examples/mesh/surface.py
rename to examples/mesh/surface_height.py
index 12d8d9b1f..1e1db7ffe 100644
--- a/examples/mesh/surface.py
+++ b/examples/mesh/surface_height.py
@@ -13,20 +13,18 @@
 import pygfx as gfx
 
 
-figure = fpl.Figure(size=(700, 560), cameras='3d', controller_types='orbit')
+figure = fpl.Figure(size=(700, 560), cameras="3d", controller_types="orbit")
 
 
 t = np.linspace(0, 6, 100).astype(np.float32)
 x = np.sin(t)
-y = np.cos(t*2)
+y = np.cos(t * 2)
 z = (x.reshape(1, -1) * x.reshape(-1, 1)) * 50  # 100x100
 
-
-mesh = figure[0, 0].add_surface(z, colors="magenta", cmap='jet')
-
+surface = figure[0, 0].add_surface(z, cmap="bwr")
 
 # figure[0, 0].axes.grids.xy.visible = True
-figure[0, 0].camera.show_object(mesh.world_object, (-2, 2, -3), up=(0, 0, 1))
+figure[0, 0].camera.show_object(surface.world_object, (-2, 2, -3), up=(0, 0, 1))
 figure.show()
 
 
diff --git a/examples/mesh/surface_ripple.py b/examples/mesh/surface_ripple.py
new file mode 100644
index 000000000..ac556bd1b
--- /dev/null
+++ b/examples/mesh/surface_ripple.py
@@ -0,0 +1,62 @@
+"""
+Surface animation
+=================
+
+Example of a surface ripple animation by setting the z-height data on every render.
+
+"""
+
+# test_example = false
+# sphinx_gallery_pygfx_docs = 'animate 6s'
+
+import fastplotlib as fpl
+import numpy as np
+
+
+figure = fpl.Figure(size=(700, 560), cameras="3d", controller_types="orbit")
+
+
+def create_ripple(shape=(100, 100), phase=0.0, freq=np.pi / 4, ampl=1.0):
+    m, n = shape
+    y, x = np.ogrid[-m / 2 : m / 2, -n / 2 : n / 2]
+    r = np.sqrt(x**2 + y**2)
+    z = (ampl * np.sin(freq * r + phase)) / np.sqrt(r + 1)
+
+    return z * 8
+
+
+z = create_ripple()
+
+# set the clim vmax
+max_z = create_ripple(phase=(np.pi / 4) - (np.pi / 2)).max()
+
+surface = figure[0, 0].add_surface(
+    z, mode="basic", cmap="viridis", clim=(-max_z, max_z)
+)
+
+figure[0, 0].camera.show_object(surface.world_object, (-1, 3, -1), up=(0, 0, 1))
+figure.show()
+
+figure[0, 0].camera.zoom = 1.15
+
+phase = 0.0
+
+
+def animate():
+    global phase
+
+    z = create_ripple(phase=phase)
+
+    surface.data = z
+
+    phase -= 0.1
+
+
+figure[0, 0].add_animations(animate)
+
+
+# NOTE: fpl.loop.run() should not be used for interactive sessions
+# See the "JupyterLab and IPython" section in the user guide
+if __name__ == "__main__":
+    print(__doc__)
+    fpl.loop.run()
diff --git a/examples/mesh/surface_sphere_ripple.py b/examples/mesh/surface_sphere_ripple.py
new file mode 100644
index 000000000..6caa03465
--- /dev/null
+++ b/examples/mesh/surface_sphere_ripple.py
@@ -0,0 +1,81 @@
+"""
+Sphere ripple animation
+=======================
+
+Example of a sphere with a ripple effect by setting the data on every render.
+
+"""
+
+# test_example = false
+# sphinx_gallery_pygfx_docs = 'animate 6s'
+
+import fastplotlib as fpl
+import numpy as np
+
+figure = fpl.Figure(size=(700, 560), cameras="3d", controller_types="orbit")
+
+# create an ellipsoid from spherical coordinates
+# see this for reference: https://mathworld.wolfram.com/SphericalCoordinates.html
+# phi and theta are swapped in this example w.r.t. the wolfram alpha description
+radius = 10
+nx = 250
+phi = np.linspace(0, np.pi * 2, num=nx, dtype=np.float32)
+ny = 250
+theta = np.linspace(0, np.pi, num=ny, dtype=np.float32)
+
+phi_grid, theta_grid = np.meshgrid(phi, theta)
+
+# convert to cartesian coordinates
+theta_grid_sin = np.sin(theta_grid)
+x = radius * np.cos(phi_grid) * theta_grid_sin * -1
+y = radius * np.cos(theta_grid)
+
+ripple_amplitude = 1.0
+ripple_frequency = 20.0
+ripple = ripple_amplitude * np.sin(ripple_frequency * theta_grid)
+
+z_ref = radius * np.sin(phi_grid) * theta_grid_sin
+z = z_ref * (1 + ripple / radius)
+
+sphere = figure[0, 0].add_surface(
+    np.dstack([x, y, z]),
+    mode="phong",
+    colors="red",
+    cmap="jet",
+)
+
+# display xz plane as a grid
+figure[0, 0].axes.grids.xy.visible = True
+figure.show()
+
+figure[0, 0].camera.show_object(sphere.world_object, (10, 1, -1), up=(0, 0, 1))
+figure[0, 0].camera.zoom = 1.3
+
+
+start = 0
+
+
+def animate():
+    global start
+    theta = np.linspace(start, start + np.pi, num=ny, dtype=np.float32)
+    _, theta_grid = np.meshgrid(phi, theta)
+    ripple = ripple_amplitude * np.sin(ripple_frequency * theta_grid)
+
+    z = z_ref * (1 + ripple / radius)
+
+    sphere.data = np.dstack([x, y, z])
+
+    start += 0.005
+
+    if start > np.pi * 2:
+        start = 0
+
+
+figure[0, 0].add_animations(animate)
+
+
+# NOTE: fpl.loop.run() should not be used for interactive sessions
+# See the "JupyterLab and IPython" section in the user guide
+if __name__ == "__main__":
+    print(__doc__)
+    fpl.loop.run()
diff --git a/examples/mesh/surface_terrain.py b/examples/mesh/surface_terrain.py
new file mode 100644
index 000000000..394e65216
--- /dev/null
+++ b/examples/mesh/surface_terrain.py
@@ -0,0 +1,36 @@
+"""
+Elevation map of the earth
+==========================
+
+Surface graphic showing elevation map of the earth
+"""
+
+# test_example = false
+# sphinx_gallery_pygfx_docs = 'code'
+
+import imageio.v3 as iio
+import fastplotlib as fpl
+import numpy as np
+
+# grayscale image of the earth where the pixel value indicates elevation
+elevation = iio.imread("https://neo.gsfc.nasa.gov/archive/bluemarble/bmng/topography/srtm_ramp2.world.5400x2700.jpg").astype(np.float32)
+elevation /= 2
+
+figure = fpl.Figure(size=(700, 560), cameras="3d", controller_types="orbit")
+
+mesh = figure[0, 0].add_surface(elevation, cmap="terrain")
+mesh.world_object.local.scale_y = -1
+
+
+figure[0, 0].axes.grids.xy.visible = True
+figure[0, 0].camera.show_object(mesh.world_object, (-4, 2, -1), up=(0, 0, 1))
+figure.show()
+
+figure[0, 0].camera.zoom = 2.5
+
+
+# NOTE: fpl.loop.run() should not be used for interactive sessions
+# See the "JupyterLab and IPython" section in the user guide
+if __name__ == "__main__":
+    print(__doc__)
+    fpl.loop.run()
diff --git a/examples/tests/testutils.py b/examples/tests/testutils.py
index 7b70defdb..18ad8ed41 100644
--- a/examples/tests/testutils.py
+++ b/examples/tests/testutils.py
@@ -25,6 +25,7 @@
     "line/*.py",
     "line_collection/*.py",
     "vectors/*.py"
+    "mesh/*.py",
     "gridplot/*.py",
     "window_layouts/*.py",
     "events/*.py",
diff --git a/fastplotlib/graphics/__init__.py b/fastplotlib/graphics/__init__.py
index 04fa95001..3d01e4a35 100644
--- a/fastplotlib/graphics/__init__.py
+++ b/fastplotlib/graphics/__init__.py
@@ -4,7 +4,7 @@
 from .image import ImageGraphic
 from .image_volume import ImageVolumeGraphic
 from ._vectors import VectorsGraphic
-from .mesh import MeshGraphic
+from .mesh import MeshGraphic, SurfaceGraphic, PolygonGraphic
 from .text import TextGraphic
 from .line_collection import LineCollection, LineStack
 
@@ -17,6 +17,8 @@
     "ImageVolumeGraphic",
     "VectorsGraphic",
     "MeshGraphic",
+    "SurfaceGraphic",
+    "PolygonGraphic",
     "TextGraphic",
     "LineCollection",
     "LineStack",
diff --git a/fastplotlib/graphics/features/__init__.py b/fastplotlib/graphics/features/__init__.py
index f745f10c8..cb0b3ab2e 100644
--- a/fastplotlib/graphics/features/__init__.py
+++ b/fastplotlib/graphics/features/__init__.py
@@ -1,10 +1,20 @@
-from ._positions_graphics import (
+from ._positions import (
     VertexColors,
     UniformColor,
     SizeSpace,
     VertexPositions,
     VertexCmap,
 )
+from ._mesh import (
+    MeshVertexPositions,
+    MeshIndices,
+    MeshCmap,
+    SurfaceData,
+    PolygonData,
+    resolve_cmap_mesh,
+    surface_data_to_mesh,
+    triangulate_polygon,
+)
 from ._line import Thickness
 from ._scatter import (
     VertexMarkers,
@@ -71,6 +81,10 @@
     "SizeSpace",
     "VertexPositions",
     "VertexCmap",
+    "MeshVertexPositions",
+    "MeshIndices",
+    "MeshCmap",
+    "SurfaceData",
     "Thickness",
     "VertexMarkers",
     "UniformMarker",
diff --git a/fastplotlib/graphics/features/_base.py b/fastplotlib/graphics/features/_base.py
index 5dec9f1e5..779310476 100644
--- a/fastplotlib/graphics/features/_base.py
+++ b/fastplotlib/graphics/features/_base.py
@@ -289,6 +289,12 @@ def _update_range(
             # the first dimension corresponding to n_datapoints
             key: int | np.ndarray[int | bool] | slice = key[0]
 
+        if isinstance(key, slice):
+            if key == slice(None):
+                # directly update full, don't need to figure out chunks
+                self.buffer.update_full()
+                return
+
         offset, size = self._parse_offset_size(key, upper_bound)
         self.buffer.update_range(offset=offset, size=size)
 
diff --git a/fastplotlib/graphics/features/_mesh.py b/fastplotlib/graphics/features/_mesh.py
new file mode 100644
index 000000000..2ec2af3af
--- /dev/null
+++ b/fastplotlib/graphics/features/_mesh.py
@@ -0,0 +1,297 @@
+from typing import Any, Sequence
+
+import numpy as np
+import pygfx
+
+from ._base import (
+    GraphicFeature,
+    GraphicFeatureEvent,
+    to_gpu_supported_dtype,
+    block_reentrance,
+)
+
+from ._positions import VertexPositions
+from ...utils.functions import get_cmap
+from ...utils.triangulation import triangulate
+
+
+def resolve_cmap_mesh(cmap) -> pygfx.TextureMap | None:
+    """Turn a user-provided in a pygfx.TextureMap, supporting 1D, 2D and 3D data."""
+
+    if cmap is None:
+        pygfx_cmap = None
+    elif isinstance(cmap, pygfx.TextureMap):
+        pygfx_cmap = cmap
+    elif isinstance(cmap, pygfx.Texture):
+        pygfx_cmap = pygfx.TextureMap(cmap)
+    elif isinstance(cmap, (str, dict)):
+        pygfx_cmap = pygfx.cm.create_colormap(get_cmap(cmap))
+    else:
+        map = np.asarray(cmap)
+        if map.ndim == 2:  # 1D plus color
+            pygfx_cmap = pygfx.cm.create_colormap(cmap)
+        else:
+            tex = pygfx.Texture(map, dim=map.ndim - 1)
+            pygfx_cmap = pygfx.TextureMap(tex)
+
+    return pygfx_cmap
+
+
+class MeshVertexPositions(VertexPositions):
+    """Manages mesh vertex positions, same as VertexPosition but data must be of shape [n, 3]"""
+
+    def _fix_data(self, data):
+        if data.ndim != 2 or data.shape[1] != 3:
+            raise ValueError(
+                f"mesh vertex positions must be of shape: [n_vertices, 3], you passed an array of shape: {data.shape}"
+            )
+
+        return to_gpu_supported_dtype(data)
+
+
+class MeshIndices(VertexPositions):
+    event_info_spec = [
+        {
+            "dict key": "key",
+            "type": "slice, index (int) or numpy-like fancy index",
+            "description": "key at which vertex indices were indexed/sliced",
+        },
+        {
+            "dict key": "value",
+            "type": "int | float | array-like",
+            "description": "new data values for indices that were changed",
+        },
+    ]
+
+    def __init__(
+        self, data: Any, isolated_buffer: bool = True, property_name: str = "indices"
+    ):
+        """
+        Manages the vertex indices buffer shown in the graphic.
+        Supports fancy indexing if the data array also supports it.
+        """
+
+        data = self._fix_data(data)
+        super().__init__(
+            data, isolated_buffer=isolated_buffer, property_name=property_name
+        )
+
+    def _fix_data(self, data):
+        if data.shape == (3,):
+            pass
+        elif data.ndim != 2 or data.shape[1] not in (3, 4):
+            raise ValueError(
+                f"indices must be of shape: [n_vertices, 3] or [n_vertices, 4], "
+                f"you passed an array of shape: {data.shape}"
+            )
+
+        return data.astype("i4")
+
+
+class MeshCmap(GraphicFeature):
+    event_info_spec = [
+        {
+            "dict key": "value",
+            "type": "str | dict | pygfx.TextureMap | pygfx.Texture | np.ndarray",
+            "description": "new cmap",
+        },
+    ]
+
+    def __init__(
+        self,
+        value: str | dict | pygfx.TextureMap | pygfx.Texture | np.ndarray | None,
+        property_name: str = "cmap",
+    ):
+        """Manages a mesh colormap"""
+
+        self._value = value
+        super().__init__(property_name=property_name)
+
+    @property
+    def value(
+        self,
+    ) -> str | dict | pygfx.TextureMap | pygfx.Texture | np.ndarray | None:
+        return self._value
+
+    @block_reentrance
+    def set_value(
+        self,
+        graphic,
+        value: str | dict | pygfx.TextureMap | pygfx.Texture | np.ndarray | None,
+    ):
+        graphic.world_object.material.map = resolve_cmap_mesh(value)
+        self._value = value
+
+        event = GraphicFeatureEvent(type=self._property_name, info={"value": value})
+        self._call_event_handlers(event)
+
+
+def surface_data_to_mesh(data: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
+    """
+    surface data to mesh positions and indices
+
+    expects data that is of shape: [m, n, 3] or [m, n]
+    """
+
+    data = np.asarray(data)
+
+    if data.ndim == 2:
+        # "image" of z values passed
+        # [m, n] -> [n_vertices, 3]
+        y = (
+            np.arange(data.shape[0])
+            .reshape(data.shape[0], 1)
+            .repeat(data.shape[1], axis=1)
+        )
+        x = (
+            np.arange(data.shape[1])
+            .reshape(1, data.shape[1])
+            .repeat(data.shape[0], axis=0)
+        )
+        positions = np.column_stack((x.ravel(), y.ravel(), data.ravel()))
+    else:
+        if data.ndim != 3:
+            raise ValueError(
+                f"expect data that is of shape: [m, n, 3], [m, n]\n"
+                f"you passed: {data.shape}"
+            )
+        if data.shape[2] != 3:
+            raise ValueError(
+                f"expect data that is of shape: [m, n, 3], [m, n]\n"
+                f"you passed: {data.shape}"
+            )
+
+        # [m, n, 3] -> [n_vertices, 3]
+        positions = data.reshape(-1, 3)
+
+    # Create faces
+    w = data.shape[1]
+    i = np.arange(data.shape[0] - 1)
+    j = np.arange(w - 1)
+
+    j, i = np.meshgrid(j, i, indexing="ij")
+    start = j.ravel() + w * i.ravel()
+
+    indices = np.column_stack([start, start + 1, start + w + 1, start + w])
+
+    return positions, indices
+
+
+class SurfaceData(GraphicFeature):
+    event_info_spec = [
+        {
+            "dict key": "value",
+            "type": "np.ndarray",
+            "description": "new surface data",
+        },
+    ]
+
+    def __init__(self, value: np.ndarray | Sequence, property_name: str = "data"):
+        self._value = np.asarray(value, dtype=np.float32)
+        super().__init__(property_name=property_name)
+
+    @property
+    def value(self) -> np.ndarray:
+        return self._value
+
+    @block_reentrance
+    def set_value(self, graphic, value: np.ndarray):
+        positions, indices = surface_data_to_mesh(value)
+
+        graphic.positions = positions
+        graphic.indices = indices
+
+        # if cmap is a 1D texture we need to set the texcoords again using new z values
+        if graphic.world_object.material.map is not None:
+            if graphic.world_object.material.map.texture.dim == 1:
+                mapcoords = positions[:, 2]
+
+                if graphic.clim is None:
+                    clim = mapcoords.min(), mapcoords.max()
+                else:
+                    clim = graphic.clim
+                mapcoords = (mapcoords - clim[0]) / (clim[1] - clim[0])
+                graphic.mapcoords = mapcoords
+
+        self._value = value
+
+        event = GraphicFeatureEvent(type=self._property_name, info={"value": value})
+        self._call_event_handlers(event)
+
+
+def triangulate_polygon(data: np.ndarray | Sequence):
+    """vertices of shape [n_vertices , 2] -> positions, indices"""
+    data = np.asarray(data, dtype=np.float32)
+
+    err_msg = (
+        f"polygon vertex data must be of shape [n_vertices, 2], you passed: {data}"
+    )
+
+    if data.ndim != 2:
+        raise ValueError(err_msg)
+    if data.shape[1] != 2:
+        raise ValueError(err_msg)
+
+    if len(data) >= 3:
+        indices = triangulate(data)
+    else:
+        indices = np.arange((0, 3), np.int32)
+
+    data = np.column_stack([data, np.zeros(data.shape[0], dtype=np.float32)])
+
+    return data, indices
+
+
+class PolygonData(GraphicFeature):
+    event_info_spec = [
+        {
+            "dict key": "value",
+            "type": "np.ndarray",
+            "description": "new polygon vertex data",
+        },
+    ]
+
+    def __init__(self, value: np.ndarray, property_name: str = "data"):
+        self._value = np.asarray(value, dtype=np.float32)
+        super().__init__(property_name=property_name)
+
+    @property
+    def value(self) -> np.ndarray:
+        return self._value
+
+    @block_reentrance
+    def set_value(self, graphic, value: np.ndarray | Sequence):
+        value = np.asarray(value, dtype=np.float32)
+
+        positions, indices = triangulate_polygon(value)
+
+        geometry = graphic.world_object.geometry
+
+        # Need larger buffer?
+        if len(positions) > geometry.positions.nitems:
+            arr = np.zeros((geometry.positions.nitems * 2, 3), np.float32)
+            geometry.positions = pygfx.Buffer(arr)
+        if len(indices) > geometry.indices.nitems:
+            arr = np.zeros((geometry.indices.nitems * 2, 3), np.int32)
+            geometry.indices = pygfx.Buffer(arr)
+
+        geometry.positions.data[: len(positions)] = positions
+        geometry.positions.data[len(positions) :] = (
+            positions[-1] if len(positions) else (0, 0, 0)
+        )
+        geometry.positions.draw_range = 0, len(positions)
+        geometry.positions.update_full()
+
+        geometry.indices.data[: len(indices)] = indices
+        geometry.indices.data[len(indices) :] = 0
+        geometry.indices.draw_range = 0, len(indices)
+        geometry.indices.update_full()
+
+        # send event
+        if len(self._event_handlers) < 1:
+            return
+
+        event = GraphicFeatureEvent(self._property_name, {"value": self.value})
+
+        # calls any events
+        self._call_event_handlers(event)
diff --git a/fastplotlib/graphics/features/_positions_graphics.py b/fastplotlib/graphics/features/_positions.py
similarity index 99%
rename from fastplotlib/graphics/features/_positions_graphics.py
rename to fastplotlib/graphics/features/_positions.py
index ae57e77d7..295d22417 100644
--- a/fastplotlib/graphics/features/_positions_graphics.py
+++ b/fastplotlib/graphics/features/_positions.py
@@ -245,8 +245,6 @@ def __init__(
         )
 
     def _fix_data(self, data):
-        # data = to_gpu_supported_dtype(data)
-
         if data.ndim == 1:
             # if user provides a 1D array, assume these are y-values
             data = np.column_stack([np.arange(data.size, dtype=data.dtype), data])
diff --git a/fastplotlib/graphics/features/utils.py b/fastplotlib/graphics/features/utils.py
index 408610e1e..aa4022052 100644
--- a/fastplotlib/graphics/features/utils.py
+++ b/fastplotlib/graphics/features/utils.py
@@ -34,8 +34,9 @@ def parse_colors(
         elif colors.ndim == 2:
             if not (colors.shape[1] in (3, 4) and colors.shape[0] == n_colors):
                 raise ValueError(
-                    "Valid array color arguments must be a single RGBA array or a stack of "
-                    "RGB or RGBA arrays for each datapoint in the shape [n_datapoints, 3] or [n_datapoints, 4]"
+                    f"Valid array color arguments must be a single RGBA array or a stack of "
+                    f"RGB or RGBA arrays for each datapoint in the shape [n_datapoints, 3] or [n_datapoints, 4].\n"
+                    f"n_datapoints is: {n_colors}, you passed a colors array of shape: {colors.shape}"
                 )
             data = colors
         else:
diff --git a/fastplotlib/graphics/mesh.py b/fastplotlib/graphics/mesh.py
index a44038371..94441db98 100644
--- a/fastplotlib/graphics/mesh.py
+++ b/fastplotlib/graphics/mesh.py
@@ -1,63 +1,43 @@
-from typing import Sequence, Any
+from typing import Sequence, Any, Literal
 
 import numpy as np
 
 import pygfx
 
 from ._positions_base import Graphic
-from .selectors import (
-    LinearRegionSelector,
-    LinearSelector,
-    RectangleSelector,
-    PolygonSelector,
-)
 from .features import (
-    BufferManager,
-    VertexPositions,
+    MeshVertexPositions,
+    MeshIndices,
+    MeshCmap,
+    SurfaceData,
+    surface_data_to_mesh,
     VertexColors,
     UniformColor,
-    VertexCmap,
+    resolve_cmap_mesh,
+    VolumeSlicePlane,
+    PolygonData,
+    triangulate_polygon,
 )
-from ..utils.functions import get_cmap
-from ..utils import quick_min_max
-
-
-def resolve_cmap(cmap):
-    """Turn a user-provided in a pygfx.TextureMap, supporting 1D, 2D and 3D data."""
-    if cmap is None:
-        pygfx_cmap = None
-    elif isinstance(cmap, pygfx.TextureMap):
-        pygfx_cmap = cmap
-    elif isinstance(cmap, pygfx.Texture):
-        pygfx_cmap = pygfx.TextureMap(cmap)
-    elif isinstance(cmap, (str, dict)):
-        pygfx_cmap = pygfx.cm.create_colormap(get_cmap(cmap))
-    else:
-        map = np.asarray(cmap)
-        if map.ndim == 2:  # 1D plus color
-            pygfx_cmap = pygfx.cm.create_colormap(cmap)
-        else:
-            tex = pygfx.Texture(map, dim=map.ndim - 1)
-            pygfx_cmap = pygfx.TextureMap(tex)
-
-    return pygfx_cmap
 
 
 class MeshGraphic(Graphic):
     _features = {
-        "positions": VertexPositions,
-        "indices": BufferManager,
-        "mapcoords": (BufferManager, None),
+        "positions": MeshVertexPositions,
+        "indices": MeshIndices,
         "colors": (VertexColors, UniformColor),
+        "cmap": MeshCmap,
     }
 
     def __init__(
         self,
         positions: Any,
         indices: Any,
+        mode: Literal["basic", "phong", "slice"] = "phong",
+        plane: tuple[float, float, float, float] = (0., 0., 1., 0.),
         colors: str | np.ndarray | Sequence = "w",
         mapcoords: Any = None,
-        cmap: str = None,
+        cmap: str | dict | pygfx.Texture | pygfx.TextureMap | np.ndarray = None,
+        clim: tuple[float, float] = None,
         isolated_buffer: bool = True,
         **kwargs,
     ):
@@ -73,6 +53,15 @@ def __init__(
             The indices into the positions that make up the triangles. Each 3
             subsequent indices form a triangle.
 
+        mode: one of "basic", "phong", "slice", default "phong"
+            * basic: illuminate mesh with only ambient lighting
+            * phong: phong lighting model, good for most use cases, see https://en.wikipedia.org/wiki/Phong_shading
+            * slice: display a slice of the mesh at the specified ``plane``
+
+        plane: (float, float, float, float), default (0., 0., 1., 0.)
+            Slice mesh at this plane. Sets (a, b, c, d) in the equation the defines a plane: ax + by + cz + d = 0.
+            Used only if `mode` = "slice". The plane is defined in world space.
+
         colors: str, array, or iterable, default "w"
             A uniform color, or the per-position colors.
 
@@ -86,6 +75,13 @@ def __init__(
             "colors". For supported colormaps see the ``cmap`` library
             catalogue: https://cmap-docs.readthedocs.io/en/stable/catalog/
             Both 1D and 2D colormaps are supported, though the mapcoords has to match the dimensionality.
+            An image can also be used, this is basically a 2D colormap.
+
+        isolated_buffer: bool, default True
+            If True, initialize a buffer with the same shape as the input data and then
+            set the data, useful if the data arrays are ready-only such as memmaps.
+            If False, the input array is itself used as the buffer - useful if the
+            array is large. In almost all cases this should be ``True``.
 
         **kwargs
             passed to :class:`.Graphic`
@@ -94,32 +90,38 @@ def __init__(
 
         super().__init__(**kwargs)
 
-        if isinstance(positions, VertexPositions):
+        if isinstance(positions, MeshVertexPositions):
             self._positions = positions
         else:
-            self._positions = VertexPositions(
+            self._positions = MeshVertexPositions(
                 positions, isolated_buffer=isolated_buffer, property_name="positions"
             )
 
-        if isinstance(positions, BufferManager):
+        if isinstance(positions, MeshIndices):
             self._indices = indices
         else:
-            self._indices = BufferManager(
+            self._indices = MeshIndices(
                 indices, isolated_buffer=isolated_buffer, property_name="indices"
             )
 
-        if mapcoords is None:
-            self._mapcoords = None
-        elif isinstance(mapcoords, BufferManager):
-            self._mapcoords = mapcoords
+        self._cmap = MeshCmap(cmap)
+
+        # Apply contrast limits. Would be nice if Pygfx mesh material had clim too! But
+        # for now we apply it as a pre-processing step.
+        if clim is None and mapcoords is not None:
+            clim = mapcoords.min(), mapcoords.max()
+
+        if mapcoords is not None:
+            mapcoords = (mapcoords - clim[0]) / (clim[1] - clim[0])
+            self._mapcoords = pygfx.Buffer(np.asarray(mapcoords, dtype=np.float32))
         else:
-            self._mapcoords = mapcoords = BufferManager(
-                mapcoords, isolated_buffer=isolated_buffer, property_name="mapcoords"
-            )
+            self._mapcoords = None
+
+        self._clim = clim
 
         uniform_color = "w"
         per_vertex_colors = False
-        pygfx_cmap = resolve_cmap(cmap)
+
         if cmap is None:
             if colors is None:
                 uniform_color = "w"
@@ -130,7 +132,6 @@ def __init__(
             elif isinstance(colors, VertexColors):
                 per_vertex_colors = True
                 self._colors = colors
-                self._colors._shared += 1
             else:
                 per_vertex_colors = True
                 self._colors = VertexColors(
@@ -140,19 +141,35 @@ def __init__(
         geometry = pygfx.Geometry(
             positions=self._positions.buffer, indices=self._indices._buffer
         )
-        material = pygfx.MeshPhongMaterial(
+
+        valid_modes = ["basic", "phong", "slice"]
+        if mode not in valid_modes:
+            raise ValueError(f"mode must be one of: {valid_modes}\nYou passed: {mode}")
+        self._mode = mode
+
+        material_cls = getattr(pygfx, f"Mesh{mode.capitalize()}Material")
+
+        if mode == "slice":
+            self._plane = VolumeSlicePlane(plane)
+            add_kwargs = {"plane": self._plane.value}
+        else:
+            # for basic and phong, maybe later we can add more of the properties
+            add_kwargs = {}
+
+        material = material_cls(
             color_mode="uniform",
             color=uniform_color,
             pick_write=True,
+            **add_kwargs,
         )
 
         # Set all the data
         if per_vertex_colors:
             geometry.colors = self._colors.buffer
-        if mapcoords is not None:
-            geometry.texcoords = self._mapcoords.buffer
-        if pygfx_cmap is not None:
-            material.map = pygfx_cmap
+        if self._mapcoords is not None:
+            geometry.texcoords = self._mapcoords
+        if cmap is not None:
+            material.map = resolve_cmap_mesh(cmap)
 
         # Decide on color mode
         # uniform = None  #: Use the uniform color (usually ``material.color``).
@@ -160,7 +177,7 @@ def __init__(
         # face = None  #: Use the per-face color specified in the geometry  (usually  ``geometry.colors``).
         # vertex_map = None  #: Use per-vertex texture coords (``geometry.texcoords``), and sample these in ``material.map``.
         # face_map = None  #: Use per-face texture coords (``geometry.texcoords``), and sample these in ``material.map``.
-        if mapcoords is not None and pygfx_cmap is not None:
+        if mapcoords is not None and cmap is not None:
             material.color_mode = "vertex_map"
         elif per_vertex_colors:
             material.color_mode = "vertex"
@@ -170,3 +187,287 @@ def __init__(
         world_object: pygfx.Mesh = pygfx.Mesh(geometry=geometry, material=material)
 
         self._set_world_object(world_object)
+
+    @property
+    def mode(self) -> Literal["basic", "phong", "slice"]:
+        """get mesh rendering mode"""
+        return self._mode
+
+    @property
+    def positions(self) -> MeshVertexPositions:
+        """Get or set the vertex positions"""
+        return self._positions
+
+    @positions.setter
+    def positions(self, new_positions):
+        self._positions[:] = new_positions
+
+    @property
+    def indices(self) -> MeshIndices:
+        """Get or set the vertex indices"""
+        return self._indices
+
+    @indices.setter
+    def indices(self, mew_indices):
+        self._indices[:] = mew_indices
+
+    @property
+    def mapcoords(self) -> np.ndarray | None:
+        """get or set the mapcoords"""
+        if self._mapcoords is not None:
+            return self._mapcoords.data
+
+    @mapcoords.setter
+    def mapcoords(self, new_mapcoords: np.ndarray | None):
+        if new_mapcoords is None:
+            self.world_object.geometry.texcoords = None
+            self._mapcoords = None
+            return
+
+        if new_mapcoords.shape == self._mapcoords.data.shape:
+            self._mapcoords.data[:] = new_mapcoords
+            self._mapcoords.update_full()
+        else:
+            # allocate new buffer
+            self._mapcoords = pygfx.Buffer(np.asarray(new_mapcoords, dtype=np.float32))
+            self.world_object.geometry.texcoords = self._mapcoords
+
+    @property
+    def clim(self) -> tuple[float, float] | None:
+        """get or set the colormap limits"""
+        return self._clim
+
+    @clim.setter
+    def clim(self, new_clim: tuple[float, float]):
+        if len(new_clim) != 2:
+            raise ValueError("clim must be a: tuple[float, float]")
+
+        self._clim = tuple(new_clim)
+
+        self.mapcoords = (self.mapcoords - self.clim[0]) / (self.clim[1] - self.clim[0])
+
+    @property
+    def colors(self) -> VertexColors | pygfx.Color:
+        """Get or set the colors"""
+        if isinstance(self._colors, VertexColors):
+            return self._colors
+
+        elif isinstance(self._colors, UniformColor):
+            return self._colors.value
+
+    @colors.setter
+    def colors(self, value: str | np.ndarray | Sequence[float] | Sequence[str]):
+        if isinstance(self._colors, VertexColors):
+            self._colors[:] = value
+
+        elif isinstance(self._colors, UniformColor):
+            self._colors.set_value(self, value)
+
+    @property
+    def cmap(self) -> str | dict | pygfx.Texture | pygfx.TextureMap | np.ndarray | None:
+        """get or set the cmap"""
+        if self._cmap is not None:
+            return self._cmap.value
+
+    @cmap.setter
+    def cmap(
+        self,
+        new_cmap: str | dict | pygfx.Texture | pygfx.TextureMap | np.ndarray | None,
+    ):
+        self._cmap.set_value(self, new_cmap)
+
+    @property
+    def plane(self) -> tuple[float, float, float, float] | None:
+        """Get or set the current slice plane. Valid only for ``"slice"`` render mode."""
+        if self.mode != "slice":
+            return
+
+        return self._plane.value
+
+    @plane.setter
+    def plane(self, value: tuple[float, float, float, float]):
+        if self.mode != "slice":
+            raise TypeError("`plane` property is only valid for `slice` render mode.")
+
+        self._plane.set_value(self, value)
+
+
+class SurfaceGraphic(MeshGraphic):
+    _features = {
+        "data": SurfaceData,
+        "colors": (VertexColors, UniformColor),
+        "cmap": MeshCmap,
+    }
+
+    def __init__(
+        self,
+        data: np.ndarray,
+        mode: Literal["basic", "phong", "slice"] = "phong",
+        colors: str | np.ndarray | Sequence = "w",
+        mapcoords: Any = None,
+        cmap: str | dict | pygfx.Texture | pygfx.TextureMap | np.ndarray = None,
+        clim: tuple[float, float] | None = None,
+        **kwargs,
+    ):
+        """
+        Create a Surface mesh Graphic
+
+        Parameters
+        ----------
+        data: array-like
+            A height-map (an image where the values indicate height, i.e. z values).
+            Can also be a [m, n, 3] to explicitly specify the x and y values in addition to the z values.
+            [m, n, 3] is a dstack of (x, y, z) values that form a grid on the xy plane.
+
+        mode: one of "basic", "phong", "slice", default "phong"
+            * basic: illuminate mesh with only ambient lighting
+            * phong: phong lighting model, good for most use cases, see https://en.wikipedia.org/wiki/Phong_shading
+
+        colors: str, array, or iterable, default "w"
+            A uniform color, or the per-position colors.
+
+        mapcoords: array-like
+            The per-position coordinates to which to apply the colormap (a.k.a. texcoords).
+            These can e.g. be some domain-specific value (mapped to [0..1] using ``clim``).
+            If not given, they will be the depth (z-coordinate) of the surface.
+
+        cmap: str, optional
+            Apply a colormap to the mesh, this overrides any argument passed to
+            "colors". For supported colormaps see the ``cmap`` library
+            catalogue: https://cmap-docs.readthedocs.io/en/stable/catalog/
+            Both 1D and 2D colormaps are supported, though the mapcoords has to match the dimensionality.
+
+        clim: tuple[float, float]
+            The colormap limits. If the mapcoords has values between e.g. 5 and 90, you want to set the clim
+            to e.g. (5, 90) or (0, 100) to determine how the values map onto the colormap.
+
+        **kwargs
+             passed to :class:`.Graphic`
+
+        """
+
+        self._data = SurfaceData(data)
+
+        positions, indices = surface_data_to_mesh(data)
+
+        cmap_tex_view = resolve_cmap_mesh(cmap)
+        if (cmap_tex_view is not None) and (mapcoords is None):
+            if cmap_tex_view.texture.dim == 1:  # 1d
+                mapcoords = positions[:, 2]
+
+            elif cmap_tex_view.texture.dim == 2:
+                mapcoords = np.column_stack((positions[:, 0], positions[:, 1])).astype(
+                    np.float32
+                )
+
+        super().__init__(
+            positions,
+            indices,
+            mode=mode,
+            colors=colors,
+            mapcoords=mapcoords,
+            cmap=cmap,
+            clim=clim,
+            **kwargs,
+        )
+
+    @property
+    def data(self) -> np.ndarray:
+        """get or set the surface data"""
+        return self._data.value
+
+    @data.setter
+    def data(self, new_data: np.ndarray):
+        self._data.set_value(self, new_data)
+
+
+class PolygonGraphic(MeshGraphic):
+    _features = {
+        "data": SurfaceData,
+        "colors": (VertexColors, UniformColor),
+        "cmap": MeshCmap,
+    }
+
+    def __init__(
+        self,
+        data: np.ndarray,
+        mode: Literal["basic", "phong"] = "basic",
+        colors: str | np.ndarray | Sequence = "w",
+        mapcoords: Any = None,
+        cmap: str | dict | pygfx.Texture | pygfx.TextureMap | np.ndarray = None,
+        clim: tuple[float, float] | None = None,
+        **kwargs,
+    ):
+        """
+        Create a polygon mesh graphic.
+
+        The data are always in the 'xy' plane. Set a rotation to display the polygon in another plane or in 3D space.
+
+        Parameters
+        ----------
+        data: array-like
+            The polygon vertices, must be of shape: [n_vertices, 2]
+
+        mode: one of "basic", "phong", "slice", default "phong"
+            * basic: illuminate mesh with only ambient lighting
+            * phong: phong lighting model, good for most use cases, see https://en.wikipedia.org/wiki/Phong_shading
+
+        colors: str, array, or iterable, default "w"
+            A uniform color, or the per-position colors.
+
+        mapcoords: array-like
+            The per-position coordinates to which to apply the colormap (a.k.a. texcoords).
+            These can e.g. be some domain-specific value (mapped to [0..1] using ``clim``).
+            If not given, they will be the depth (z-coordinate) of the surface.
+
+        cmap: str, optional
+            Apply a colormap to the mesh, this overrides any argument passed to
+            "colors". For supported colormaps see the ``cmap`` library
+            catalogue: https://cmap-docs.readthedocs.io/en/stable/catalog/
+            Both 1D and 2D colormaps are supported, though the mapcoords has to match the dimensionality.
+
+        clim: tuple[float, float]
+            The colormap limits. If the mapcoords has values between e.g. 5 and 90, you want to set the clim
+            to e.g. (5, 90) or (0, 100) to determine how the values map onto the colormap.
+
+        **kwargs
+             passed to :class:`.Graphic`
+        """
+
+        positions, indices = triangulate_polygon(data)
+
+        self._data = PolygonData(positions)
+
+        super().__init__(
+            positions,
+            indices,
+            mode=mode,
+            colors=colors,
+            mapcoords=mapcoords,
+            cmap=cmap,
+            clim=clim,
+            **kwargs,
+        )
+
+    @property
+    def data(self) -> np.ndarray:
+        """get or set the polygon vertex data"""
+        return self._data.value
+
+    @data.setter
+    def data(self, new_data: np.ndarray | Sequence):
+        self._data.set_value(self, new_data)
+
+    @property
+    def clim(self) -> tuple[float, float] | None:
+        """get or set the colormap limits"""
+        return self._clim
+
+    @clim.setter
+    def clim(self, new_clim: tuple[float, float]):
+        if len(new_clim) != 2:
+            raise ValueError("clim must be a: tuple[float, float]")
+
+        self._clim = tuple(new_clim)
+
+        self.mapcoords = (self.mapcoords - self.clim[0]) / (self.clim[1] - self.clim[0])
diff --git a/fastplotlib/layouts/_graphic_methods_mixin.py b/fastplotlib/layouts/_graphic_methods_mixin.py
index 403cd0a22..ae227a3f0 100644
--- a/fastplotlib/layouts/_graphic_methods_mixin.py
+++ b/fastplotlib/layouts/_graphic_methods_mixin.py
@@ -8,7 +8,6 @@
 
 from ..graphics import *
 from ..graphics._base import Graphic
-from ..graphics.mesh import resolve_cmap as resolve_mesh_cmap
 
 
 class GraphicMethodsMixin:
@@ -437,19 +436,26 @@ def add_mesh(
         self,
         positions: Any,
         indices: Any,
-        mapcoords: Any = None,
+        mode: Literal["basic", "phong", "slice"] = "phong",
+        plane: tuple[float, float, float, float] = (0, 0, 1, 0),
         colors: Union[str, numpy.ndarray, Sequence] = "w",
-        cmap: str = None,
+        mapcoords: Any = None,
+        cmap: (
+            str
+            | dict
+            | pygfx.resources._texture.Texture
+            | pygfx.resources._texturemap.TextureMap
+            | numpy.ndarray
+        ) = None,
         isolated_buffer: bool = True,
         **kwargs,
     ) -> MeshGraphic:
         """
 
-        Create a mesh Graphic
+        Create a mesh Graphic.
 
         Parameters
         ----------
-
         positions: array-like
             The 3D positions of the vertices.
 
@@ -457,6 +463,15 @@ def add_mesh(
             The indices into the positions that make up the triangles. Each 3
             subsequent indices form a triangle.
 
+        mode: one of "basic", "phong", "slice", default "phong"
+            * basic: illuminate mesh with only ambient lighting
+            * phong: phong lighting model, good for most use cases, see https://en.wikipedia.org/wiki/Phong_shading
+            * slice: display a slice of the mesh at the specified ``plane``
+
+        plane: (float, float, float, float), default (0, 0, -1, 0)
+            Slice mesh at this plane. Sets (a, b, c, d) in the equation the defines a plane: ax + by + cz + d = 0.
+            Used only if `mode` = "slice". The plane is defined in world space.
+
         colors: str, array, or iterable, default "w"
             A uniform color, or the per-position colors.
 
@@ -470,6 +485,13 @@ def add_mesh(
             "colors". For supported colormaps see the ``cmap`` library
             catalogue: https://cmap-docs.readthedocs.io/en/stable/catalog/
             Both 1D and 2D colormaps are supported, though the mapcoords has to match the dimensionality.
+            An image can also be used, this is basically a 2D colormap.
+
+        isolated_buffer: bool, default True
+            If True, initialize a buffer with the same shape as the input data and then
+            set the data, useful if the data arrays are ready-only such as memmaps.
+            If False, the input array is itself used as the buffer - useful if the
+            array is large. In almost all cases this should be ``True``.
 
         **kwargs
             passed to :class:`.Graphic`
@@ -480,6 +502,8 @@ def add_mesh(
             MeshGraphic,
             positions,
             indices,
+            mode,
+            plane,
             colors,
             mapcoords,
             cmap,
@@ -487,25 +511,36 @@ def add_mesh(
             **kwargs,
         )
 
-    def add_surface(
+    def add_polygon(
         self,
-        data: Any,
+        data: numpy.ndarray,
+        mode: Literal["basic", "phong"] = "basic",
         colors: Union[str, numpy.ndarray, Sequence] = "w",
         mapcoords: Any = None,
-        cmap: str = None,
+        cmap: (
+            str
+            | dict
+            | pygfx.resources._texture.Texture
+            | pygfx.resources._texturemap.TextureMap
+            | numpy.ndarray
+        ) = None,
         clim: tuple[float, float] | None = None,
-        isolated_buffer: bool = True,
         **kwargs,
-    ) -> MeshGraphic:
+    ) -> PolygonGraphic:
         """
 
-        Create a mesh Graphic
+        Create a polygon mesh graphic.
+
+        The data are always in the 'xy' plane. Set a rotation to display the polygon in another plane or in 3D space.
 
         Parameters
         ----------
         data: array-like
-            A height-map (an image where the values indicate height, i.e. z values).
-            Can also be a 3-tuple to explicitly specify the x and y values in addition to the z values.
+            The polygon vertices, must be of shape: [n_vertices, 2]
+
+        mode: one of "basic", "phong", "slice", default "phong"
+            * basic: illuminate mesh with only ambient lighting
+            * phong: phong lighting model, good for most use cases, see https://en.wikipedia.org/wiki/Phong_shading
 
         colors: str, array, or iterable, default "w"
             A uniform color, or the per-position colors.
@@ -528,111 +563,9 @@ def add_surface(
         **kwargs
              passed to :class:`.Graphic`
 
-
         """
-
-        def check_z(z):
-            if z.ndim != 2:
-                raise ValueError("Z must be a 2D array.")
-
-        # In VisVis (https://github.com/almarklein/visvis/blob/main/visvis/functions/surf.py)
-        # the tuple can be 2-element or 4-element, to pass a color per z-value.
-        # I disabled that by commenting related logic. Maybe it can be enabled someday.
-
-        if isinstance(data, (tuple, list)):
-            if len(data) == 1:
-                z = numpy.asanyarray(data[0])
-                check_z(z)
-                y = numpy.arange(z.shape[0])
-                x = numpy.arange(z.shape[1])
-                c = None
-            # elif len(data) == 2:
-            #     z = numpy.asanyarray(data[0])
-            #     c = numpy.asanyarray(data[1])
-            #     check_z(z)
-            #     y = numpy.arange(z.shape[0])
-            #     x = numpy.arange(z.shape[1])
-            elif len(args) == 3:
-                x = numpy.asanyarray(data[0])
-                y = numpy.asanyarray(data[1])
-                z = numpy.asanyarray(data[2])
-                check_z(z)
-                c = None
-            # elif len(args) == 4:
-            #     x = numpy.asanyarray(data[0])
-            #     y = numpy.asanyarray(data[1])
-            #     z = numpy.asanyarray(data[2])
-            #     c = numpy.asanyarray(data[3])
-            #     check_z(z)
-            else:
-                raise ValueError(
-                    "Surface tuple has invalid number of elements (need 1-4)."
-                )
-        else:
-            z = numpy.asanyarray(data)
-            check_z(z)
-            y = numpy.arange(z.shape[0])
-            x = numpy.arange(z.shape[1])
-            c = None
-
-        # Set y vertices
-        if y.shape == (z.shape[0],):
-            y = y.reshape(z.shape[0], 1).repeat(z.shape[1], axis=1)
-        elif y.shape != z.shape:
-            raise ValueError(
-                "Y must have same shape as Z, or be 1D with length of rows of Z."
-            )
-
-        # Set x vertices
-        if x.shape == (z.shape[1],):
-            x = x.reshape(1, z.shape[1]).repeat(z.shape[0], axis=0)
-        elif x.shape != z.shape:
-            raise ValueError(
-                "X must have same shape as Z, or be 1D with length of columns of Z."
-            )
-
-        # Set vertices
-        positions = numpy.column_stack((x.ravel(), y.ravel(), z.ravel()))
-
-        # Create texcoords
-        cmap = resolve_mesh_cmap(cmap)
-        if mapcoords is None:
-            if cmap.texture.dim == 1:  # 1d
-                mapcoords = z.ravel()
-            elif cmap.texture.dim == 2:
-                mapcoords = numpy.column_stack((x.ravel(), y.ravel())).astype(
-                    numpy.float32
-                )
-            # Apply contrast limits. Would be nice if Pygfx mesh material had clim too! But
-            # for now we apply it as a pre-processing step.
-            if clim is None and mapcoords is not None:
-                clim = mapcoords.min(), mapcoords.max()
-            mapcoords = (mapcoords - clim[0]) / (clim[1] - clim[0])
-        else:
-            raise ValueError("C must have same shape as Z, or be 3D array.")
-
-        # Create faces
-        w = z.shape[1]
-        i = numpy.arange(z.shape[0] - 1)
-        indices = numpy.row_stack(
-            [
-                numpy.column_stack(
-                    (j + w * i, j + 1 + w * i, j + 1 + w * (i + 1), j + w * (i + 1))
-                )
-                for j in range(w - 1)
-            ]
-        )
-        indices = indices.astype("i4")
-
         return self._create_graphic(
-            MeshGraphic,
-            positions,
-            indices,
-            colors,
-            mapcoords,
-            cmap,
-            isolated_buffer,
-            **kwargs,
+            PolygonGraphic, data, mode, colors, mapcoords, cmap, clim, **kwargs
         )
 
     def add_scatter(
@@ -789,6 +722,63 @@ def add_scatter(
             **kwargs,
         )
 
+    def add_surface(
+        self,
+        data: numpy.ndarray,
+        mode: Literal["basic", "phong", "slice"] = "phong",
+        colors: Union[str, numpy.ndarray, Sequence] = "w",
+        mapcoords: Any = None,
+        cmap: (
+            str
+            | dict
+            | pygfx.resources._texture.Texture
+            | pygfx.resources._texturemap.TextureMap
+            | numpy.ndarray
+        ) = None,
+        clim: tuple[float, float] | None = None,
+        **kwargs,
+    ) -> SurfaceGraphic:
+        """
+
+        Create a Surface mesh Graphic
+
+        Parameters
+        ----------
+        data: array-like
+            A height-map (an image where the values indicate height, i.e. z values).
+            Can also be a [m, n, 3] to explicitly specify the x and y values in addition to the z values.
+
+        mode: one of "basic", "phong", "slice", default "phong"
+            * basic: illuminate mesh with only ambient lighting
+            * phong: phong lighting model, good for most use cases, see https://en.wikipedia.org/wiki/Phong_shading
+
+        colors: str, array, or iterable, default "w"
+            A uniform color, or the per-position colors.
+
+        mapcoords: array-like
+            The per-position coordinates to which to apply the colormap (a.k.a. texcoords).
+            These can e.g. be some domain-specific value (mapped to [0..1] using ``clim``).
+            If not given, they will be the depth (z-coordinate) of the surface.
+
+        cmap: str, optional
+            Apply a colormap to the mesh, this overrides any argument passed to
+            "colors". For supported colormaps see the ``cmap`` library
+            catalogue: https://cmap-docs.readthedocs.io/en/stable/catalog/
+            Both 1D and 2D colormaps are supported, though the mapcoords has to match the dimensionality.
+
+        clim: tuple[float, float]
+            The colormap limits. If the mapcoords has values between e.g. 5 and 90, you want to set the clim
+            to e.g. (5, 90) or (0, 100) to determine how the values map onto the colormap.
+
+        **kwargs
+             passed to :class:`.Graphic`
+
+
+        """
+        return self._create_graphic(
+            SurfaceGraphic, data, mode, colors, mapcoords, cmap, clim, **kwargs
+        )
+
     def add_text(
         self,
         text: str,