faster bounding rectangle for imposter spheres

CHANGELOG.md

@@ -5,7 +5,9 @@ Note that since we don't clearly distinguish between a public and private interf
 
 
 ## [Unreleased]
+
 - Add some elements support for `guessElementSymbolString` function
+- Faster bounding rectangle calculation for imposter spheres
 
 ## [v3.38.3] - 2023-07-29
 

src/mol-gl/shader/spheres.vert.ts

@@ -17,6 +17,7 @@ precision highp int;
 
 uniform mat4 uModelView;
 uniform mat4 uInvProjection;
+uniform float uIsOrtho;
 
 uniform vec2 uTexDim;
 uniform sampler2D tPositionGroup;
@@ -30,45 +31,27 @@ varying vec3 vPointViewPosition;
 
 #include matrix_scale
 
-const mat4 D = mat4(
-    1.0, 0.0, 0.0, 0.0,
-    0.0, 1.0, 0.0, 0.0,
-    0.0, 0.0, 1.0, 0.0,
-    0.0, 0.0, 0.0, -1.0
-);
-
 /**
- * Compute point size and center using the technique described in:
- * "GPU-Based Ray-Casting of Quadratic Surfaces" http://dl.acm.org/citation.cfm?id=2386396
- * by Christian Sigg, Tim Weyrich, Mario Botsch, Markus Gross.
+ * Bounding rectangle of a clipped, perspective-projected 3D Sphere.
+ * Michael Mara, Morgan McGuire. 2013
+ *
+ * Specialization by Arseny Kapoulkine, MIT License Copyright (c) 2018
+ * https://github.com/zeux/niagara
  */
-void quadraticProjection(const in float radius, const in vec3 position, const in vec2 mapping){
-    vec2 xbc, ybc;
-
-    mat4 T = mat4(
-        radius, 0.0, 0.0, 0.0,
-        0.0, radius, 0.0, 0.0,
-        0.0, 0.0, radius, 0.0,
-        position.x, position.y, position.z, 1.0
-    );
-
-    mat4 R = transpose4(uProjection * uModelView * aTransform * T);
-    float A = dot(R[3], D * R[3]);
-    float B = -2.0 * dot(R[0], D * R[3]);
-    float C = dot(R[0], D * R[0]);
-    xbc[0] = (-B - sqrt(B * B - 4.0 * A * C)) / (2.0 * A);
-    xbc[1] = (-B + sqrt(B * B - 4.0 * A * C)) / (2.0 * A);
-    float sx = abs(xbc[0] - xbc[1]) * 0.5;
-
-    A = dot(R[3], D * R[3]);
-    B = -2.0 * dot(R[1], D * R[3]);
-    C = dot(R[1], D * R[1]);
-    ybc[0] = (-B - sqrt(B * B - 4.0 * A * C)) / (2.0 * A);
-    ybc[1] = (-B + sqrt(B * B - 4.0 * A * C)) / (2.0 * A);
-    float sy = abs(ybc[0] - ybc[1]) * 0.5;
-
-    gl_Position.xy = vec2(0.5 * (xbc.x + xbc.y), 0.5 * (ybc.x + ybc.y));
-    gl_Position.xy -= mapping * vec2(sx, sy);
+void sphereProjection(const in vec3 p, const in float r, const in vec2 mapping) {
+    vec3 pr = p * r;
+    float pzr2 = p.z * p.z - r * r;
+
+    float vx = sqrt(p.x * p.x + pzr2);
+    float minx = ((vx * p.x - pr.z) / (vx * p.z + pr.x)) * uProjection[0][0];
+    float maxx = ((vx * p.x + pr.z) / (vx * p.z - pr.x)) * uProjection[0][0];
+
+    float vy = sqrt(p.y * p.y + pzr2);
+    float miny = ((vy * p.y - pr.z) / (vy * p.z + pr.y)) * uProjection[1][1];
+    float maxy = ((vy * p.y + pr.z) / (vy * p.z - pr.y)) * uProjection[1][1];
+
+    gl_Position.xy = vec2(maxx + minx, maxy + miny) * -0.5;
+    gl_Position.xy -= mapping * vec2(maxx - minx, maxy - miny) * 0.5;
     gl_Position.xy *= gl_Position.w;
 }
 
@@ -106,8 +89,14 @@ void main(void){
         mvCorner.xy += mapping * vRadius;
         gl_Position = uProjection * mvCorner;
     #else
-        gl_Position = uProjection * vec4(mvPosition.xyz, 1.0);
-        quadraticProjection(vRadius, position, mapping);
+        if (uIsOrtho == 1.0) {
+            vec4 mvCorner = vec4(mvPosition.xyz, 1.0);
+            mvCorner.xy += mapping * vRadius;
+            gl_Position = uProjection * mvCorner;
+        } else {
+            gl_Position = uProjection * vec4(mvPosition.xyz, 1.0);
+            sphereProjection(mvPosition.xyz, vRadius, mapping);
+        }
     #endif
 
     vec4 vPoint4 = uInvProjection * gl_Position;