1+ #version 330
2+ #pragma vscode_glsllint_stage : vert
3+
4+ /* *
5+ Compute covariance relatec computation & projection etc
6+ This emits the 2 2D basis vector for geometry shader's quad construction
7+ The RGBA color of the GS is also passed as is to the geometry shader
8+ */
9+
10+ // uniform int H;
11+ // uniform int W;
12+ // uniform float n;
13+ // uniform float f;
14+ // uniform mat3x3 K;
15+
16+ uniform mat4x4 P;
17+ uniform mat4x4 VM;
18+ uniform vec2 focal;
19+ uniform vec2 principal;
20+ uniform vec2 basisViewport;
21+
22+ // uniform float discardAlpha = 0.1;
23+ uniform float discardAlpha = 0.0001 ;
24+ // uniform float discardAlpha = 0.15;
25+ uniform float maxScreenSpaceSplatSize = 2048.0 ;
26+ uniform float sqrt8 = sqrt (8 );
27+ // uniform float sqrt8 = 3;
28+
29+ layout (location = 0 ) in vec3 aPos; // xyz
30+ layout (location = 1 ) in vec3 aCov0_3; // cov6
31+ layout (location = 2 ) in vec3 aCov3_6; // cov6
32+ layout (location = 3 ) in vec4 aColor; // rgba
33+
34+ out vec2 basisVector0;
35+ out vec2 basisVector1;
36+ out vec4 gColor; // pass through
37+ out float gDepth; // pass through
38+
39+ void main() {
40+ if (aColor.a < discardAlpha) {
41+ gColor.a = 0.0 ; // will not emit things in geometry shader
42+ return ;
43+ }
44+
45+ // Compute the view and clip space coordinates of the center of the ellipse
46+ vec4 viewCenter = VM * vec4 (aPos, 1.0 );
47+ vec4 clipCenter = P * vec4 (aPos, 1.0 );
48+ clipCenter = clipCenter / clipCenter.w; // perspective division
49+
50+ // Construct the 3D covariance matrix
51+ mat3 Vrk = mat3 (
52+ aCov0_3[0 ], aCov0_3[1 ], aCov0_3[2 ],
53+ aCov0_3[1 ], aCov3_6[0 ], aCov3_6[1 ],
54+ aCov0_3[2 ], aCov3_6[1 ], aCov3_6[2 ]);
55+
56+ // Construct the Jacobian of the affine approximation of the projection matrix. It will be used to transform the
57+ float width = 1 / basisViewport[0 ];
58+ float height = 1 / basisViewport[1 ];
59+ float fx = focal[0 ];
60+ float fy = focal[1 ];
61+ float cx = principal[0 ];
62+ float cy = principal[1 ];
63+ float x = viewCenter.x;
64+ float y = viewCenter.y;
65+ float z = viewCenter.z;
66+
67+ float tan_fovx = 0.5 * width / fx;
68+ float tan_fovy = 0.5 * height / fy;
69+ float lim_x_pos = (width - cx) / fx + 0.3 * tan_fovx;
70+ float lim_x_neg = cx / fx + 0.3 * tan_fovx;
71+ float lim_y_pos = (height - cy) / fy + 0.3 * tan_fovy;
72+ float lim_y_neg = cy / fy + 0.3 * tan_fovy;
73+
74+ float rz = 1.0 / z;
75+ float rz2 = rz * rz;
76+ float tx = z * min (lim_x_pos, max (- lim_x_neg, x * rz));
77+ float ty = z * min (lim_y_pos, max (- lim_y_neg, y * rz));
78+
79+ mat3 J = mat3 (
80+ fx * rz, 0 ., - fx * tx * rz2,
81+ 0 ., fy * rz, - fy * ty * rz2,
82+ 0 ., 0 ., 0 .);
83+
84+ // Concatenate the projection approximation with the model-view transformation
85+ mat3 W = transpose (mat3 (VM));
86+ mat3 T = W * J;
87+
88+ // Transform the 3D covariance matrix (Vrk) to compute the 2D covariance matrix
89+ mat3 cov2Dm = transpose (T) * Vrk * T;
90+ cov2Dm[0 ][0 ] += 0.3 ;
91+ cov2Dm[1 ][1 ] += 0.3 ;
92+
93+ // We are interested in the upper-left 2x2 portion of the projected 3D covariance matrix because
94+ // we only care about the X and Y values. We want the X-diagonal, cov2Dm[0][0],
95+ // the Y-diagonal, cov2Dm[1][1], and the correlation between the two cov2Dm[0][1]. We don't
96+ // need cov2Dm[1][0] because it is a symetric matrix.
97+ vec3 cov2Dv = vec3 (cov2Dm[0 ][0 ], cov2Dm[0 ][1 ], cov2Dm[1 ][1 ]);
98+
99+ // We now need to solve for the eigen-values and eigen vectors of the 2D covariance matrix
100+ // so that we can determine the 2D basis for the splat. This is done using the method described
101+ // here: https://people.math.harvard.edu/~knill/teaching/math21b2004/exhibits/2dmatrices/index.html
102+ // After calculating the eigen-values and eigen-vectors, we calculate the basis for rendering the splat
103+ // by normalizing the eigen-vectors and then multiplying them by (sqrt(8) * eigen-value), which is
104+ // equal to scaling them by sqrt(8) standard deviations.
105+ //
106+ // This is a different approach than in the original work at INRIA. In that work they compute the
107+ // max extents of the projected splat in screen space to form a screen-space aligned bounding rectangle
108+ // which forms the geometry that is actually rasterized. The dimensions of that bounding box are 3.0
109+ // times the maximum eigen-value, or 3 standard deviations. They then use the inverse 2D covariance
110+ // matrix (called 'conic') in the CUDA rendering thread to determine fragment opacity by calculating the
111+ // full gaussian: exp(-0.5 * (X - mean) * conic * (X - mean)) * splat opacity
112+ float a = cov2Dv.x;
113+ float d = cov2Dv.z;
114+ float b = cov2Dv.y;
115+ float D = a * d - b * b;
116+
117+ if (D <= 0.0 || cov2Dv.x <= 0.0 || cov2Dv.z <= 0.0 ) {
118+ // Illegal cov matrix, this point should be pruned with zero gradients
119+ gColor.a = 0.0 ; // will not emit things
120+ return ;
121+ }
122+
123+ float trace = a + d;
124+ float traceOver2 = 0.5 * trace;
125+ float term2 = sqrt (max (0 .1f, traceOver2 * traceOver2 - D));
126+ float eigenValue0 = traceOver2 + term2;
127+ float eigenValue1 = traceOver2 - term2;
128+
129+ // if (eigenValue1 < 0) {
130+ // gColor.a = 0.0; // will not emit things
131+ // return;
132+ // }
133+ if (eigenValue0 <= 0.01 || eigenValue1 <= 0.01 || eigenValue0 < eigenValue1 || (eigenValue0 / eigenValue1) > 10000.0 ) {
134+ gColor.a = 0.0 ; // will not emit things
135+ return ;
136+ }
137+
138+ vec2 eigenVector0 = normalize (vec2 (b, eigenValue0 - a));
139+ // since the eigen vectors are orthogonal, we derive the second one from the first
140+ vec2 eigenVector1 = vec2 (eigenVector0.y, - eigenVector0.x);
141+
142+ // We use sqrt(8) standard deviations instead of 3 to eliminate more of the splat with a very low opacity.
143+ basisVector0 = eigenVector0 * min (sqrt8 * sqrt (eigenValue0), maxScreenSpaceSplatSize);
144+ basisVector1 = eigenVector1 * min (sqrt8 * sqrt (eigenValue1), maxScreenSpaceSplatSize);
145+
146+ gl_Position = clipCenter; // doing a perspective projection to ndc space
147+ gColor = aColor; // passing through
148+ gDepth = clipCenter.z; // passing through
149+ }
0 commit comments