spherical volume rendering refactored

.gitignore

@@ -1,3 +1,5 @@
+# temp files
+*.swp
 # Byte-compiled / optimized / DLL files
 __pycache__/
 *.py[cod]
@@ -8,6 +10,8 @@ __pycache__/
 
 # Cythonized files
 yt_idv/utilities.c
+yt_idv/coordinate_utilities.c
+
 
 # imgui files
 imgui.ini

README.md

@@ -29,6 +29,7 @@ loaded in yt.  It is written to provide both scripting and interactive access.
     * Block and grid outlines
 * Support for sub-selections of data via the yt data selection interface
 * Integration with the [ipywidgets](https://ipywidgets.readthedocs.org/) ``Image`` widget.
+* Direct volume rendering of block-AMR data in spherical coordinates.
 
 ## Examples
 

docs/coordinate_systems.rst

@@ -0,0 +1,55 @@
+.. highlight:: shell
+
+================================
+Non-Cartesian Coordinate Systems
+================================
+
+Initial support for volume rendering of data defined in 3D non-cartesian coordinate systems
+was added in yt_idv 0.5.0 for block-AMR data in spherical coordinates. While not all
+rendering methods and annotations are supported for spherical coordinates, yt_idv can
+directly calculate maximum intensity projections, integrative projections and projection with custom
+transfer functions without any pre-interpolation or re-gridding (with some caveats).
+
+PUT A NICE SCREEN SHOT HERE.
+
+---------------------------------------------------
+Volume Rendering of Data in a Spherical Coordinates
+---------------------------------------------------
+
+Overview
+
+#. pre-calculation of cartesian bounding boxes of the amr blocks, accelerated with cython.
+#. load block data as textures in normalized native coordinates (just like cartesian)
+#. pass down cartesian boudning boxes and spherical bounding boxes of blocks as vertex attributes
+#. standard slab test with the cartesian bounding boxes for ray-element intersections
+#. step along ray, calculate spherical coordinates (in the shader!), evalulate texture (discarding points outside the actual speherical volume element).
+
+
+Mention: Pre-processor switches in order to re-use carteisan shaders efficiently.
+
+The main limitation in this approach is in stepping along the ray: intersection with the
+ cartesian bounding box of a spherical element does not gaurentee intersection with the enclosed
+spherical elemnt. So a fairly large number of sample points along the ray must be used, increasing
+the computational cost to achieve the same fidelity as rendering data that is natively in
+cartesian coordinates.
+
+----------------------------
+Notes on further development
+----------------------------
+
+Further contributions are welcome for adding support for the remaining 3d non-cartesian coordinate systems
+that yt supports that are not yet supported here (3d cylindrical, 3d geographic) as well as for adding
+support for non-cartesian coordinate systems in additional yt_idv components.
+
+****************************************
+Supporting additional coordinate systems
+****************************************
+
+Describe coordinate_utilities.pyx (and .pxd), how to add a new coordinate system.
+
+***************************************
+Supporting additional rendering methods
+***************************************
+
+Adding additional coordiante system conversions to shaders (with pre-processor
+directives).

docs/index.rst

@@ -9,6 +9,7 @@ Welcome to yt_idv's documentation!
    installation
    usage
    scene
+   coordinate_systems
    examples
    modules
    contributing

examples/amr_spherical_volume_rendering.py

@@ -0,0 +1,127 @@
+import argparse
+
+import numpy as np
+import yt
+
+import yt_idv
+
+# yt reminder: phi is the polar angle (0 to 2pi)
+# theta is the angle from north (0 to pi)
+# coord ordering here will be r, phi, theta
+bbox_options = {
+    "partial": np.array([[0.5, 1.0], [0.0, np.pi / 3], [np.pi / 4, np.pi / 2]]),
+    "whole": np.array([[0.0, 1.0], [0.0, 2 * np.pi], [0, np.pi]]),
+    "shell": np.array([[0.7, 1.0], [0.0, 2 * np.pi], [0, np.pi]]),
+    "north_hemi": np.array([[0.1, 1.0], [0.0, 2 * np.pi], [0, 0.5 * np.pi]]),
+    "north_shell": np.array([[0.8, 1.0], [0.0, 2 * np.pi], [0, 0.5 * np.pi]]),
+    "south_hemi": np.array([[0.1, 1.0], [0.0, 2 * np.pi], [0.5 * np.pi, np.pi]]),
+    "ew_hemi": np.array([[0.1, 1.0], [0.0, np.pi], [0.0, np.pi]]),
+    "quadrant_shell": np.array([[0.6, 1.0], [0.0, np.pi / 2], [0.0, np.pi / 2]]),
+}
+
+if __name__ == "__main__":
+
+    parser = argparse.ArgumentParser(
+        prog="amr_spherical_volume_rendering",
+        description="Loads an example spherical dataset in yt_idv",
+    )
+
+    msg = f"The geometry subset to generate: one of {list(bbox_options.keys())}"
+    parser.add_argument("-d", "--domain", default="partial", help=msg)
+    msg = (
+        "The field to plot. Provide a comma-separated string with field_type,field "
+        "e.g., to plot the field tuple ('index', 'phi'): \n "
+        "    $ python amr_spherical_volume_rendering.py -f index,x "
+        "\nIf a single string is provided, a field type of gas is assumed."
+    )
+    parser.add_argument("-f", "--field", default="index,phi", help=msg)
+    parser.add_argument(
+        "-np", "--nprocs", default=64, help="number of grids to decompose domain to"
+    )
+    parser.add_argument(
+        "-sz", "--size", default=256, help="dimensions, will be (size, size size)"
+    )
+
+    args = parser.parse_args()
+
+    sz = (int(args.size),) * 3
+    fake_data = {"density": np.random.random(sz)}
+
+    field = str(args.field).split(",")
+    if len(field) == 1:
+        field = ("gas", str(field).strip())
+    elif len(field) == 2:
+        field = (field[0].strip(), field[1].strip())
+    else:
+        raise RuntimeError(
+            "Unexpected field formatting. Provide a single string"
+            " to provide just a field (will assume field type "
+            " of 'gas', or a comma separated string to provide a "
+            "field type and a field"
+        )
+
+    if args.domain not in bbox_options:
+        raise RuntimeError(
+            f"domain must be one of {list(bbox_options.keys())}, found {args.domain}"
+        )
+    bbox = bbox_options[args.domain]
+
+    nprocs = int(args.nprocs)
+
+    ds = yt.load_uniform_grid(
+        fake_data,
+        sz,
+        bbox=bbox,
+        nprocs=nprocs,
+        geometry="spherical",
+        axis_order=("r", "phi", "theta"),
+        length_unit=1,
+    )
+
+    phi_c = ds.quan(ds.domain_center[ds.coordinates.axis_id["phi"]].d, "")
+    theta_c = ds.quan(ds.domain_center[ds.coordinates.axis_id["theta"]].d, "")
+    rmax = ds.domain_right_edge[ds.coordinates.axis_id["r"]]
+    phi_f = ds.quan(15.0 * np.pi / 180.0, "")
+    theta_f = ds.quan(15.0 * np.pi / 180.0, "")
+    min_val = ds.quan(0.1, "")
+
+    def _tube(field, data):
+        phi = data["index", "phi"]
+        theta = data["index", "theta"]
+        tube = (1 - min_val) * np.exp(-(((theta - theta_c) / theta_f) ** 2))
+        tube = tube * np.exp(-(((phi - phi_c) / phi_f) ** 2))
+        return tube + min_val
+
+    ds.add_field(
+        name=("stream", "tube"),
+        function=_tube,
+        sampling_type="local",
+    )
+
+    def _r_rev(field, data):
+        r = data["index", "r"]
+        return rmax - r
+
+    ds.add_field(
+        name=("stream", "r_rev"),
+        function=_r_rev,
+        sampling_type="local",
+    )
+
+    if field not in ds.field_list + ds.derived_field_list:
+        spaces = " " * 8
+        fld_list_str = f"\n{spaces}".join(str(fld) for fld in ds.field_list)
+        drv_fld_list_str = f"\n{spaces}".join(str(fld) for fld in ds.derived_field_list)
+        raise RuntimeError(
+            f"field {field} not in field_list or derived_field_list:\n"
+            f"\n    ds.field_list:\n{spaces}{fld_list_str}"
+            f"\n    ds.derived_field_list:\n{spaces}{drv_fld_list_str}"
+        )
+
+    rc = yt_idv.render_context(height=800, width=800, gui=True)
+    sg = rc.add_scene(ds, field, no_ghost=True)
+    rc.scene.components[0].sample_factor = 20.0
+    rc.scene.components[0].cmap_log = False
+    rc.scene.components[0]._reset_cmap_bounds()
+
+    rc.run()

yt_idv/coordinate_utilities.pxd

@@ -0,0 +1,57 @@
+cimport numpy as np
+
+
+cdef inline np.float64_t fmax(np.float64_t f0, np.float64_t f1) noexcept nogil:
+    if f0 > f1: return f0
+    return f1
+
+cdef inline np.float64_t fmin(np.float64_t f0, np.float64_t f1) noexcept nogil:
+    if f0 < f1: return f0
+    return f1
+
+cdef (np.float64_t, np.float64_t, np.float64_t) _spherical_to_cartesian(np.float64_t r,
+                           np.float64_t theta,
+                           np.float64_t phi) noexcept nogil
+
+
+cdef (np.float64_t, np.float64_t, np.float64_t) _cartesian_to_spherical(np.float64_t x,
+                           np.float64_t y,
+                           np.float64_t z) noexcept nogil
+
+cdef class MixedCoordBBox:
+    cdef int get_cartesian_bbox(self,
+                                np.float64_t pos0,
+                                np.float64_t pos1,
+                                np.float64_t pos2,
+                                np.float64_t dpos0,
+                                np.float64_t dpos1,
+                                np.float64_t dpos2,
+                                np.float64_t xyz_i[3],
+                                np.float64_t dxyz_i[3]
+                                ) noexcept nogil
+
+
+cdef class SphericalMixedCoordBBox(MixedCoordBBox):
+    cdef int get_cartesian_bbox(
+                        self,
+                        np.float64_t pos0,
+                        np.float64_t pos1,
+                        np.float64_t pos2,
+                        np.float64_t dpos0,
+                        np.float64_t dpos1,
+                        np.float64_t dpos2,
+                        np.float64_t[3] xyz_i,
+                        np.float64_t[3] dxyz_i
+                        ) noexcept nogil
+
+    cdef void _get_cartesian_bbox(
+                    self,
+                    np.float64_t pos0,
+                    np.float64_t pos1,
+                    np.float64_t pos2,
+                    np.float64_t dpos0,
+                    np.float64_t dpos1,
+                    np.float64_t dpos2,
+                    np.float64_t[3] xyz_i,
+                    np.float64_t[3] dxyz_i
+                    ) noexcept nogil