Spherical volume rendering

yt_idv/scene_data/block_collection.py

@@ -41,6 +41,39 @@ def add_data(self, field, no_ghost=False):
         # Every time we change our data source, we wipe all existing ones.
         # We now set up our vertices into our current data source.
         vert, dx, le, re = [], [], [], []
+
+        # spherical-only: SHOULD ONLY DO THIS IF DATASET IS SPHERICAL
+        # normal vectors and plane constants for min/max phi face
+        phi_plane_le = []
+        phi_plane_re = []
+
+        axis_id = self.data_source.ds.coordinates.axis_id
+
+        def calculate_phi_normal_plane(le_or_re):
+            # calculates the parameters for a plane parallel to the z-axis and
+            # passing through the x-y values defined by the polar angle phi
+            # the result plane defines the boundary on +/-phi faces.
+            #
+            # eq of plane:
+            # X dot N = d where N is a normal vector, d is a constant,
+            #             X is a position vector
+            # this function returns N and d.
+            phi = le_or_re[axis_id["phi"]]
+            theta = le_or_re[axis_id["theta"]]
+            r = le_or_re[axis_id["r"]]
+
+            # prob can/should use a yt func here to ensure consistency....
+            z = r * np.cos(theta)
+            r_xy = r * np.sin(theta)
+            x = r_xy * np.cos(phi)
+            y = r_xy * np.sin(phi)
+            pt = np.array([x, y, z])
+
+            z_hat = np.array([0, 0, 1])
+            normal_vec = np.cross(pt, z_hat)
+            d = np.dot(normal_vec, pt)
+            return normal_vec, d
+
         self.min_val = +np.inf
         self.max_val = -np.inf
         if self.scale:
@@ -65,6 +98,12 @@ def add_data(self, field, no_ghost=False):
             le.append(block.LeftEdge.tolist())
             re.append(block.RightEdge.tolist())
 
+            normal, const = calculate_phi_normal_plane(le[-1])
+            phi_plane_le.append(np.array([*normal, const]))
+
+            normal, const = calculate_phi_normal_plane(re[-1])
+            phi_plane_re.append(np.array([*normal, const]))
+
         for (g, node, (sl, _dims, _gi)) in self.data_source.tiles.slice_traverse():
             block = node.data
             self.blocks_by_grid[g.id - g._id_offset].append((id(block), i))
@@ -95,6 +134,15 @@ def add_data(self, field, no_ghost=False):
             VertexAttribute(name="in_right_edge", data=re)
         )
 
+        phi_plane_re = np.array(phi_plane_re, dtype="f4")
+        phi_plane_le = np.array(phi_plane_le, dtype="f4")
+        self.vertex_array.attributes.append(
+            VertexAttribute(name="phi_plane_le", data=phi_plane_le)
+        )
+        self.vertex_array.attributes.append(
+            VertexAttribute(name="phi_plane_re", data=phi_plane_re)
+        )
+
         # Now we set up our textures
         self._load_textures()
 

yt_idv/shaders/grid_expand.geom.glsl

@@ -6,6 +6,8 @@ layout ( triangle_strip, max_vertices = 14 ) out;
 flat in vec3 vdx[];
 flat in vec3 vleft_edge[];
 flat in vec3 vright_edge[];
+flat in vec4 vphi_plane_le[];
+flat in vec4 vphi_plane_re[];
 
 flat out vec3 dx;
 flat out vec3 left_edge;
@@ -22,6 +24,9 @@ flat in vec4 vv_model[];
 
 flat out ivec3 texture_offset;
 
+flat out vec4 phi_plane_le;
+flat out vec4 phi_plane_re;
+
 // https://stackoverflow.com/questions/28375338/cube-using-single-gl-triangle-strip
 // suggests that the triangle strip we want for the cube is
 
@@ -73,6 +78,8 @@ void main() {
         inverse_pmvm = vinverse_pmvm[0];
         inverse_proj = vinverse_proj[0];
         inverse_mvm = vinverse_mvm[0];
+        phi_plane_le = vphi_plane_le[0];
+        phi_plane_re = vphi_plane_re[0];
         dx = vdx[0];
         v_model = vec4(current_pos, 1.0);
         texture_offset = ivec3(0);

yt_idv/shaders/grid_position.vert.glsl

@@ -5,6 +5,10 @@ in vec4 model_vertex;
 in vec3 in_dx;
 in vec3 in_left_edge;
 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;
+
 flat out vec4 vv_model;
 flat out mat4 vinverse_proj;  // always cartesian
 flat out mat4 vinverse_mvm;  // always cartesian
@@ -13,6 +17,9 @@ flat out vec3 vdx;
 flat out vec3 vleft_edge;
 flat out vec3 vright_edge;
 
+flat out vec4 vphi_plane_le;
+flat out vec4 vphi_plane_re;
+
 vec3 transform_vec3(vec3 v) {
     if (is_spherical) {
         // in yt, phi is the polar angle from (0, 2pi), theta is the azimuthal
@@ -46,4 +53,7 @@ void main()
     vdx = vec3(in_dx);
     vleft_edge = vec3(in_left_edge);
     vright_edge = vec3(in_right_edge);
+
+    vphi_plane_le = vec4(phi_plane_le);
+    vphi_plane_re = vec4(phi_plane_re);
 }

yt_idv/shaders/ray_tracing.frag.glsl

@@ -8,6 +8,11 @@ flat in mat4 inverse_pmvm; // always cartesian
 flat in ivec3 texture_offset;
 out vec4 output_color;
 
+flat in vec4 phi_plane_le;
+flat in vec4 phi_plane_re;
+
+const float INFINITY = 1. / 0.;
+const float PI = 3.1415926535897932384626433832795;
 bool within_bb(vec3 pos)
 {
     bvec3 left =  greaterThanEqual(pos, left_edge);
@@ -45,6 +50,106 @@ vec4 cleanup_phase(in vec4 curr_color, in vec3 dir, in float t0, in float t1);
 //   void (vec3 tex_curr_pos, inout vec4 curr_color, float tdelta, float t,
 //         vec3 direction);
 
+float get_ray_plane_intersection(vec3 p_normal, float p_constant, vec3 ray_origin, vec3 ray_dir)
+{
+    float n_dot_u = dot(p_normal, ray_dir);
+    float n_dot_ro = dot(p_normal, ray_origin);
+    // add check for n_dot_u == 0 (ray is parallel to plane)
+    if (n_dot_u == float(0.0))
+    {
+        // the ray is parallel to the plane. there are either no intersections
+        // or infinite intersections. In the second case, we are guaranteed
+        // to intersect one of the other faces, so we can drop this plane.
+        return INFINITY;
+    }
+
+    return (p_constant - n_dot_ro) / n_dot_u;
+}
+
+vec2 get_ray_sphere_intersection(float r, vec3 ray_origin, vec3 ray_dir)
+{
+    float dir_dot_dir = dot(ray_dir, ray_dir);
+    float ro_dot_ro = dot(ray_origin, ray_origin);
+    float dir_dot_ro = dot(ray_dir, ray_origin);
+    float rsq = r * r; // could be calculated in vertex shader
+
+    float a_2 = 2.0 * dir_dot_dir;
+    float b = 2.0 * dir_dot_ro;
+    float c =  ro_dot_ro - rsq;
+    float determinate = b*b - 2.0 * a_2 * c;
+    float cutoff_val = 0.0;
+    if (determinate < cutoff_val)
+    {
+        return vec2(INFINITY, INFINITY);
+    }
+    else if (determinate == cutoff_val)
+    {
+        return vec2(-b / a_2, INFINITY);
+    }
+    else
+    {
+        return vec2((-b - sqrt(determinate))/ a_2, (-b + sqrt(determinate))/ a_2);
+    }
+
+}
+
+vec2 get_ray_cone_intersection(float theta, vec3 ray_origin, vec3 ray_dir)
+{
+    // note: cos(theta) and vhat could be calculated in vertex shader
+    float costheta;
+    vec3 vhat;
+    if (theta > PI/2.0)
+    {
+        // if theta is past PI/2, the cone will point in negative z and the
+        // half angle should be measured from the -z axis, not +z.
+        vhat = vec3(0.0, 0.0, -1.0);
+        costheta = cos(PI - theta);
+    }
+    else
+    {
+        vhat = vec3(0.0, 0.0, 1.0);
+        costheta = cos(theta);
+    }
+    float cos2t = pow(costheta, 2);
+    // note: theta = PI/2.0 is well defined. determinate = 0 in that case and
+    // the cone becomes a plane in x-y.
+
+    float dir_dot_vhat = dot(ray_dir, vhat);
+    float dir_dot_dir = dot(ray_dir, ray_dir);
+    float ro_dot_vhat = dot(ray_origin, vhat);
+    float ro_dot_dir = dot(ray_origin, ray_dir);
+    float ro_dot_ro = dot(ray_origin, ray_dir);
+
+    float a_2 = 2.0*(pow(dir_dot_vhat, 2) - dir_dot_dir * cos2t);
+    float b = 2.0 * (ro_dot_vhat * dir_dot_vhat - ro_dot_dir*cos2t);
+    float c = pow(ro_dot_vhat, 2) - ro_dot_ro*cos2t;
+    float determinate = b*b - 2.0 * a_2 * c;
+    if (determinate < 0.0)
+    {
+        return vec2(INFINITY, INFINITY);
+    }
+    else if (determinate == 0.0)
+    {
+        return vec2(-b / a_2, INFINITY);
+    }
+    else
+    {
+        vec2 t = vec2((-b - sqrt(determinate))/ a_2, (-b + sqrt(determinate))/ a_2);
+
+        // check that t values are not for the shadow cone
+        float cutoff_t = (-1.0 * ro_dot_vhat) / dir_dot_vhat;
+        if (t[0] < cutoff_t)
+        {
+            t[0] = INFINITY;
+        }
+        else if (t[1] < cutoff_t)
+        {
+            t[1] = INFINITY;
+        }
+        return t;
+    }
+}
+
 void spherical_coord_shortcircuit()
 {
     vec3 ray_position = v_model.xyz; // now spherical
@@ -53,6 +158,22 @@ void spherical_coord_shortcircuit()
     vec3 dir = -normalize(camera_pos.xyz - ray_position_xyz);
     vec4 curr_color = vec4(0.0);
 
+    // intersections
+    vec2 t_sphere_outer = get_ray_sphere_intersection(right_edge[id_r], ray_position_xyz, dir);
+    if (t_sphere_outer[0] == INFINITY && t_sphere_outer[1] == INFINITY)
+    {
+        // if there are no intersections with the outer sphere, then there
+        // will be no intersections with the remaining faces and we can stop
+        // looking.
+        discard;
+    }
+
+    vec2 t_sphere_inner = get_ray_sphere_intersection(left_edge[id_r], ray_position_xyz, dir);
+    float t_p_1 = get_ray_plane_intersection(vec3(phi_plane_le), phi_plane_le[3], ray_position_xyz, dir);
+    float t_p_2 = get_ray_plane_intersection(vec3(phi_plane_re), phi_plane_re[3], ray_position_xyz, dir);
+    vec2 t_cone_outer = get_ray_cone_intersection(right_edge[id_theta], ray_position_xyz, dir);
+    vec2 t_cone_inner= get_ray_cone_intersection(left_edge[id_theta], ray_position_xyz, dir);
+
     // do the texture evaluation in the native coordinate space
     vec3 range = (right_edge + dx/2.0) - (left_edge - dx/2.0);
     vec3 nzones = range / dx;