Implement ContourPlot3D

mathics/builtin/drawing/plot_plot3d.py

@@ -113,6 +113,8 @@ def eval(
             if isinstance(self, ParametricPlot3D) and len(plot_options.ranges) == 1:
                 # ParametricPlot3D with one independent variable generating a curve
                 default_plot_points = (1000,)
+            elif isinstance(self, ContourPlot3D):
+                default_plot_points = (50, 50, 50)
             elif plot.use_vectorized_plot:
                 default_plot_points = (200, 200)
             else:
@@ -237,6 +239,40 @@ class ContourPlot(_Plot3D):
     graphics_class = Graphics
 
 
+class ContourPlot3D(_Plot3D):
+    """
+    <url>:Isosurface: https://en.wikipedia.org/wiki/Isosurface</url> (
+    <url>:WMA link: https://reference.wolfram.com/language/ref/ContourPlot3D.html</url>)
+    <dl>
+      <dt>'ContourPlot3D'[$f(x,y,z)$, {$x$, $x_{min}$, $x_{max}$}, {$y$, $y_{min}$, $y_{max}$, {$y$, $y_{min}$, $y_{max}$}]
+      <dd>creates a three-dimensional contour plot of $f(x,y,z)$ over the specified region on $x$, $y$, and $z$.
+
+      See <url>:Drawing Option and Option Values:
+    /doc/reference-of-built-in-symbols/graphics-and-drawing/drawing-options-and-option-values
+    </url> for a list of Plot options.
+
+    </dl>
+
+    >> ContourPlot3D[x ^ 2 + y ^ 2 - z ^ 2, {x, -1, 1}, {y, -1, 1}, {z, -1, 1}]
+     = ContourPlot3D[x ^ 2 + y ^ 2 - z ^ 2, {x, -1, 1}, {y, -1, 1}, {z, -1, 1}]
+
+    Multiple isosurfaces (3d contours) of a second degree equation form conical suraces, hyperboloids in this case.
+    """
+
+    requires = ["skimage"]
+    summary_text = "creates a 3d contour plot"
+    expected_args = 4
+    options = _Plot3D.options3d | {
+        "Contours": "Automatic",
+        "BoxRatios": "{1,1,1}",
+        "Mesh": "None",
+    }
+    # TODO: other options?
+
+    many_functions = False
+    graphics_class = Graphics3D
+
+
 class DensityPlot(_Plot3D):
     """
     <url>:heat map:https://en.wikipedia.org/wiki/Heat_map</url>(<url>:WMA link: https://reference.wolfram.com/language/ref/DensityPlot.html</url>)

mathics/eval/drawing/colors.py

@@ -181,9 +181,11 @@ def gradient_palette(
 ]
 
 
-def palette_color_directive(palette, i):
+def palette_color_directive(palette, i, opacity=None):
     """returns a directive in a form suitable for GraphicsGenerator.add_directives"""
     """ for setting the color of an entire entity such as a line or surface """
     rgb = palette[i % len(palette)]
     rgb = [c / 255.0 for c in rgb]
+    if opacity is not None:
+        rgb.append(opacity)
     return [SymbolRGBColor, *rgb]

mathics/eval/drawing/plot3d.py

@@ -515,6 +515,13 @@ def eval_ContourPlot(
     return None
 
 
+def eval_ContourPlot3D(
+    plot_options,
+    evaluation: Evaluation,
+):
+    return None
+
+
 def eval_ParametricPlot3D(
     plot_options,
     evaluation: Evaluation,

mathics/eval/drawing/plot3d_vectorized.py

@@ -14,6 +14,7 @@
 from mathics.core.systemsymbols import (
     SymbolAbsoluteThickness,
     SymbolEqual,
+    SymbolFull,
     SymbolNone,
     SymbolRGBColor,
     SymbolSubtract,
@@ -222,6 +223,33 @@ def emit(graphics, i, xyzs, _, quads):
     return make_surfaces(plot_options, evaluation, dim=2, is_complex=False, emit=emit)
 
 
+def equations_to_contours(plot_options):
+    """
+    For contour plots, convert functions of the form f==g to f-g,
+    and adjust contours ot [0] and disable background
+    """
+    for i, f in enumerate(plot_options.functions):
+        if hasattr(f, "head") and f.head == SymbolEqual:
+            f = Expression(SymbolSubtract, *f.elements[0:2])
+            plot_options.functions[i] = f
+            plot_options.contours = [0]
+            plot_options.background = False
+
+
+def choose_contour_levels(plot_options, vmin, vmax, default):
+    levels = plot_options.contours
+    if isinstance(levels, str):
+        # TODO: need to pick "nice" number so levels have few digits
+        # an odd number ensures there is a contour at 0 if range is balanced
+        levels = default
+    if isinstance(levels, int):
+        # computed contour levels have equal distance between them,
+        # and half that between first/last contours and vmin/vmax
+        dv = (vmax - vmin) / levels
+        levels = vmin + np.arange(levels) * dv + dv / 2
+    return levels
+
+
 @Timer("eval_ContourPlot")
 def eval_ContourPlot(
     plot_options,
@@ -230,21 +258,14 @@ def eval_ContourPlot(
     import skimage.measure
 
     # whether to show a background similar to density plot, except quantized
-    background = len(plot_options.functions) == 1
-    contour_levels = plot_options.contours
+    plot_options.background = len(plot_options.functions) == 1
 
-    # convert fs of the form a==b to a-b, inplicit contour level 0
-    plot_options.functions = list(plot_options.functions)  # so we can modify it
-    for i, f in enumerate(plot_options.functions):
-        if hasattr(f, "head") and f.head == SymbolEqual:
-            f = Expression(SymbolSubtract, *f.elements[0:2])
-            plot_options.functions[i] = f
-            contour_levels = [0]
-            background = False
+    # adjust functions, contours, and background for functions of the form f==g
+    equations_to_contours(plot_options)
 
     def emit(graphics, i, xyzs, _, quads):
         # set line color
-        if background:
+        if plot_options.background:
             # showing a background, so just black lines
             color_directive = [SymbolRGBColor, 0, 0, 0]
         else:
@@ -259,17 +280,8 @@ def emit(graphics, i, xyzs, _, quads):
         zs = xyzs[:, 2]  # this is a linear list matching with quads
 
         # process contour_levels
-        levels = contour_levels
         zmin, zmax = np.min(zs), np.max(zs)
-        if isinstance(levels, str):
-            # TODO: need to pick "nice" number so levels have few digits
-            # an odd number ensures there is a contour at 0 if range is balanced
-            levels = 9
-        if isinstance(levels, int):
-            # computed contour levels have equal distance between them,
-            # and half that between first/last contours and zmin/zmax
-            dz = (zmax - zmin) / levels
-            levels = zmin + np.arange(levels) * dz + dz / 2
+        levels = choose_contour_levels(plot_options, zmin, zmax, default=9)
 
         # one contour line per contour level
         for level in levels:
@@ -285,7 +297,7 @@ def emit(graphics, i, xyzs, _, quads):
                 graphics.add_linexyzs(segment)
 
         # background is solid colors between contour lines
-        if background:
+        if plot_options.background:
             with Timer("contour background"):
                 # use masks and fancy indexing to assign (lo+hi)/2 to all zs between lo and hi
                 zs[zs <= levels[0]] = zmin
@@ -301,6 +313,81 @@ def emit(graphics, i, xyzs, _, quads):
     return make_surfaces(plot_options, evaluation, dim=2, is_complex=False, emit=emit)
 
 
+@Timer("eval_ContourPlot3D")
+def eval_ContourPlot3D(
+    plot_options,
+    evaluation: Evaluation,
+):
+    import skimage.measure
+
+    # adjust functions, contours, and background for functions of the form f==g
+    equations_to_contours(plot_options)
+
+    # pull out options
+    _, xmin, xmax = plot_options.ranges[0]
+    _, ymin, ymax = plot_options.ranges[1]
+    _, zmin, zmax = plot_options.ranges[2]
+    nx, ny, nz = plot_options.plot_points
+    names = [strip_context(str(r[0])) for r in plot_options.ranges]
+
+    # compile the functions
+    with Timer("compile"):
+        compiled_functions = compile_exprs(evaluation, plot_options.functions, names)
+
+    # construct (nx,ny,nz) grid
+    xs = np.linspace(xmin, xmax, nx)
+    ys = np.linspace(ymin, ymax, ny)
+    zs = np.linspace(zmin, zmax, nz)
+    xs, ys, zs = np.meshgrid(xs, ys, zs)
+
+    # just one function
+    function = compiled_functions[0]
+
+    # compute function values fs over the grid
+    with Timer("compute fs"):
+        fs = function(**{n: v for n, v in zip(names, [xs, ys, zs])})
+
+    # process contour_levels
+    fmin, fmax = np.min(fs), np.max(fs)
+    levels = choose_contour_levels(plot_options, fmin, fmax, default=7)
+
+    # find contour for each level and emit it
+    graphics = GraphicsGenerator(dim=3)
+    for i, level in enumerate(levels):
+        color_directive = palette_color_directive(palette3, i)
+        graphics.add_directives(color_directive)
+
+        # find contour for this level
+        with Timer("3d contours"):
+            verts, faces, normals, values = skimage.measure.marching_cubes(fs, level)
+            verts[:, (0, 1)] = verts[:, (1, 0)]  # skimage bug?
+
+        # marching_cubes gives back coordinates relative to grid unit, so rescale to x, y, z
+        offset = [xmin, ymin, zmin]
+        scale = [
+            (xmax - xmin) / (nx - 1),
+            (ymax - ymin) / (ny - 1),
+            (zmax - zmin) / (nz - 1),
+        ]
+        verts = np.array(offset) + verts * np.array(scale)
+
+        # WL is 1-based
+        faces += 1
+
+        # emit faces as GraphicsComplex
+        graphics.add_complex(verts, lines=None, polys=faces, colors=None)
+
+        # emit mesh as GraphicsComplex
+        # TODO: this should share vertices with previous GraphicsComplex
+        if plot_options.mesh is SymbolFull:
+            # TODO: each segment emitted twice - is there reasonable way to fix?
+            lines = np.array([faces[:, [0, 1]], faces[:, [1, 2]], faces[:, [2, 0]]])
+            graphics.add_directives([SymbolRGBColor, 0, 0, 0])
+            graphics.add_complex(verts, lines=lines, polys=None, colors=None)
+
+    return graphics
+
+
 @Timer("complex colors")
 def complex_colors(zs, s=None):
     # hue depends on phase

test/builtin/drawing/test_plot_detail.py

@@ -370,15 +370,16 @@ def test_yaml(parms):
     one_test(**parms)
 
 
-def do_test_all(fns):
+def do_test_all(fns, names=None):
     # several of these tests failed on pyodide due to apparent differences
     # in numpy (and/or the blas library backing it) between pyodide and other platforms
     # including numerical instability, different data types (integer vs real)
     # the tests above seem so far to be ok on pyodide, but generally they are
     # simpler than these doc_tests
     if not pyodide:
         for parms in all_yaml_tests_generator(fns):
-            one_test(**parms)
+            if not names or parms["name"] in names:
+                one_test(**parms)
 
 
 if __name__ == "__main__":
@@ -391,7 +392,14 @@ def do_test_all(fns):
     UPDATE_MODE = args.update
 
     try:
-        do_test_all(args.files)
+        if args.files:
+            for fn in args.files:
+                split = fn.split(":")
+                fn = split[0]
+                names = split[1].split(",") if len(split) > 1 else None
+                do_test_all([fn], names)
+        else:
+            do_test_all(YAML_TESTS)
     except AssertionError as oops:
         print(oops)
         print("FAIL")