very messy... but improvements????

examples/amr_spherical_volume_rendering.py

@@ -12,8 +12,8 @@
 # coord ordering here will be r, phi, theta
 
 bbox_options = {
-    "partial": np.array([[0.5, 1.0], [0.0, np.pi / 4], [np.pi / 4, np.pi / 2]]),
-    "whole": np.array([[0.1, 1.0], [0.0, 2 * np.pi], [0, np.pi]]),
+    "partial": np.array([[0.5, 1.0], [0.0, np.pi / 3], [np.pi / 4, np.pi / 2]]),
+    "whole": np.array([[0., 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]]),
     "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]]),
@@ -35,16 +35,17 @@
         fake_data,
         sz,
         bbox=bbox,
-        nprocs=4096,
+        nprocs=64,
         geometry=("spherical", ("r", "phi", "theta")),
         length_unit="m",
     )
 
     rc = yt_idv.render_context(height=800, width=800, gui=True)
     dd = ds.all_data()
     dd.max_level = 1
-    sg = rc.add_scene(ds, ("index", "r"), no_ghost=True)
+    # sg = rc.add_scene(ds, ("index", "r"), no_ghost=True)
     # sg = rc.add_scene(ds, ("index", "theta"), no_ghost=True)
-    # sg = rc.add_scene(ds, ("index", "phi"), no_ghost=True)
+    sg = rc.add_scene(ds, ("index", "phi"), no_ghost=True)
+    # sg = rc.add_scene(ds, ("stream", "density"), no_ghost=True)
     sg.camera.focus = [0.0, 0.0, 0.0]
     rc.run()

yt_idv/scene_data/block_collection.py

@@ -119,6 +119,11 @@ def _set_geometry_attributes(self, le, re):
         if self._yt_geom_str == "cartesian":
             return
         elif self._yt_geom_str == "spherical":
+            from yt_idv.coordinate_utilities  import (
+                SphericalMixedCoordBBox,
+                cartesian_bboxes,
+            )
+
             axis_id = self.data_source.ds.coordinates.axis_id
             phi_plane_le = phi_normal_planes(le, axis_id, cast_type="f4")
             phi_plane_re = phi_normal_planes(re, axis_id, cast_type="f4")
@@ -128,6 +133,50 @@ def _set_geometry_attributes(self, le, re):
             self.vertex_array.attributes.append(
                 VertexAttribute(name="phi_plane_re", data=phi_plane_re)
             )
+            # cartesian bbox: very ugly, just testing...
+            print("calculating cartesian bounding boxes.")
+            widths = re - le
+            centers = (le + re) / 2.0
+            bbox_handler = SphericalMixedCoordBBox()
+            r = centers[:, axis_id['r']].astype('float64')
+            theta = centers[:, axis_id['theta']].astype('float64')
+            phi = centers[:, axis_id['phi']].astype('float64')
+            dr = widths[:, axis_id['r']].astype('float64')
+            dtheta = widths[:, axis_id['theta']].astype('float64')
+            dphi = widths[:, axis_id['phi']].astype('float64')
+            x = np.full(r.shape, np.nan, dtype="float64")
+            y = np.full(r.shape, np.nan, dtype="float64")
+            z = np.full(r.shape, np.nan, dtype="float64")
+            dx = np.full(r.shape, np.nan, dtype="float64")
+            dy = np.full(r.shape, np.nan, dtype="float64")
+            dz = np.full(r.shape, np.nan, dtype="float64")
+            cartesian_bboxes(bbox_handler, r, theta, phi, dr, dtheta, dphi, x, y, z, dx, dy, dz)
+            le_cart = np.column_stack([x - dx/2., y - dy/2., z - dz/2.])
+            re_cart = np.column_stack([x + dx / 2., y + dy / 2., z + dz / 2.])
+
+            # normalize to viewport in (0, 1)
+            domain_le = []
+            domain_re = []
+            domain_wid = []
+            print("scaling cart bounds")
+            for idim in range(3):
+                domain_le.append(le_cart[:, idim].min())
+                domain_re.append(re_cart[:, idim].max())
+                domain_wid.append(domain_re[idim] - domain_le[idim])
+                le_cart[:, idim] = (le_cart[:, idim] - domain_le[idim]) / domain_wid[idim]
+                re_cart[:, idim] = (re_cart[:, idim] - domain_le[idim]) / domain_wid[idim]
+
+            le_cart = np.asarray(le_cart)
+            re_cart = np.asarray(re_cart)
+
+
+            self.vertex_array.attributes.append(
+                VertexAttribute(name="le_cart", data=le_cart.astype('f4'))
+            )
+            self.vertex_array.attributes.append(
+                VertexAttribute(name="re_cart", data=re_cart.astype('f4'))
+            )
+
         else:
             raise NotImplementedError(
                 f"{self.name} does not implement {self._yt_geom_str} geometries."

yt_idv/shaders/grid_expand.geom.glsl

@@ -6,12 +6,16 @@ layout ( triangle_strip, max_vertices = 14 ) out;
 flat in vec3 vdx[];
 flat in vec3 vleft_edge[];
 flat in vec3 vright_edge[];
+flat in vec3 vleft_edge_cart[];
+flat in vec3 vright_edge_cart[];
 flat in vec4 vphi_plane_le[];
 flat in vec4 vphi_plane_re[];
 
 flat out vec3 dx;
 flat out vec3 left_edge;
 flat out vec3 right_edge;
+flat out vec3 left_edge_cart;
+flat out vec3 right_edge_cart;
 flat out mat4 inverse_proj; // always cartesian
 flat out mat4 inverse_mvm; // always cartesian
 flat out mat4 inverse_pmvm; // always cartesian
@@ -62,16 +66,19 @@ void main() {
 
     vec4 center = gl_in[0].gl_Position;  // always cartesian
 
-    vec3 width = vright_edge[0] - vleft_edge[0];
+    vec3 width = vright_edge_cart[0] - vleft_edge_cart[0];
 
     vec4 newPos;
     vec3 current_pos;
 
     for (int i = 0; i < 14; i++) {
         // walks through each vertex of the triangle strip, emit it. need to
         // emit gl_Position in cartesian, but pass native coords out in v_model
-        current_pos = vec3(vleft_edge[0] + width * arrangement[aindex[i]]);
-        newPos = vec4(transform_vec3(current_pos), 1.0); // cartesian
+
+        // hm, this seems wrong. maybe should use the cartesian bounding box
+        // nodes for building the triangle strip primitive?
+        current_pos = vec3(vleft_edge_cart[0] + width * arrangement[aindex[i]]);
+        newPos = vec4(current_pos, 1.0); // cartesian
         gl_Position = projection * modelview * newPos;  // cartesian
         left_edge = vleft_edge[0];
         right_edge = vright_edge[0];
@@ -80,9 +87,13 @@ void main() {
         inverse_mvm = vinverse_mvm[0];
         phi_plane_le = vphi_plane_le[0];
         phi_plane_re = vphi_plane_re[0];
+        left_edge_cart = vleft_edge_cart[0];
+        right_edge_cart = vright_edge_cart[0];
         dx = vdx[0];
         v_model = vec4(current_pos, 1.0);
         texture_offset = ivec3(0);
         EmitVertex();
     }
+
+    // why no endprimitive?
 }

yt_idv/shaders/grid_position.vert.glsl

@@ -9,6 +9,10 @@ in vec3 in_right_edge;
 in vec4 phi_plane_le;  // first 3 elements are the normal vector, final is constant
 in vec4 phi_plane_re;
 
+// pre-computed cartesian le, re
+in vec3 le_cart;
+in vec3 re_cart;
+
 flat out vec4 vv_model;
 flat out mat4 vinverse_proj;  // always cartesian
 flat out mat4 vinverse_mvm;  // always cartesian
@@ -19,41 +23,51 @@ flat out vec3 vright_edge;
 
 flat out vec4 vphi_plane_le;
 flat out vec4 vphi_plane_re;
+flat out vec3 vleft_edge_cart;
+flat out vec3 vright_edge_cart;
 
-vec3 transform_vec3(vec3 v) {
-    if (is_spherical) {
-        // in yt, phi is the polar angle from (0, 2pi), theta is the azimuthal
-        // angle (0, pi). the id_ values below are uniforms that depend on the
-        // yt dataset coordinate ordering
-        return vec3(
-            v[id_r] * sin(v[id_theta]) * cos(v[id_phi]),
-            v[id_r] * sin(v[id_theta]) * sin(v[id_phi]),
-            v[id_r] * cos(v[id_theta])
-        );
-    } else {
-        return v;
-    }
-}
-
-vec4 transform_vec4(vec4 v) {
-    vec3 v3 = transform_vec3(vec3(v));
-    return vec4(v3[0], v3[1], v3[2], v[3]);
-}
+//vec3 transform_vec3(vec3 v) {
+//    if (is_spherical) {
+//        // in yt, phi is the polar angle from (0, 2pi), theta is the azimuthal
+//        // angle (0, pi). the id_ values below are uniforms that depend on the
+//        // yt dataset coordinate ordering
+//        return vec3(
+//            v[id_r] * sin(v[id_theta]) * cos(v[id_phi]),
+//            v[id_r] * sin(v[id_theta]) * sin(v[id_phi]),
+//            v[id_r] * cos(v[id_theta])
+//        );
+//    } else {
+//        return v;
+//    }
+//}
+//
+//vec4 transform_vec4(vec4 v) {
+////    vec3 v3 = transform_vec3(vec3(v));
+//    return vec4(v3[0], v3[1], v3[2], v[3]);
+//}
 
 void main()
 {
     // camera uniforms:
     // projection, modelview
     vv_model = model_vertex;
     vinverse_proj = inverse(projection);
+
     // inverse model-view-matrix
     vinverse_mvm = inverse(modelview);
     vinverse_pmvm = inverse(projection * modelview);
-    gl_Position = projection * modelview * transform_vec4(model_vertex);
+    gl_Position = projection * modelview * vv_model; //transform_vec4(model_vertex);
+
+    // spherical
     vdx = vec3(in_dx);
     vleft_edge = vec3(in_left_edge);
     vright_edge = vec3(in_right_edge);
 
+    // pre-computed phi-plane info
     vphi_plane_le = vec4(phi_plane_le);
     vphi_plane_re = vec4(phi_plane_re);
+
+    // cartesian bounding boxes
+    vleft_edge_cart = vec3(le_cart);
+    vright_edge_cart = vec3(re_cart);
 }