Projecting 3D gaussian's transform matrix to a 2D viewspace transform matrix

Question

Essentially I'm wanting to project:
A 3D transformation matrix of a Gaussian (or in this case for simplicity a unit sphere) that includes scale, rotation and translation in world space TO
A 2D transformation matrix for a unit ellipse in screen space

I'm attempting to render Gaussians efficiently within the limitations of a game's API. It's a Unity game, and all I really have access to is meshes and shaders.

How do I go about doing this?

I'm also wondering if this might be a question for mathematics.

---

Currently after caressing ChatGPT for awhile I managed to get it to spit out some snippets I managed to piece together, it almost works. It works perfectly with orthographic projection or perspective at distance, but perspective has major error up close. Here's my code:

float2x2 Project3DGaussianTo2D(
    float3x3 gaussianLocalTransform,
    float3x3 cameraRotationWorldToCamera,
    float3 cameraRelativeDirection
) {
    // Step 1: 3D Covariance Matrix
    float3x3 G = gaussianLocalTransform;
    float3x3 C = mul(G, transpose(G));

    // Step 2: Transform covariance to camera/view space
    float3x3 R = cameraRotationWorldToCamera;
    float3x3 C_view = mul(R, mul(C, transpose(R)));

    // Step 3: Projection Jacobian J
    float3 center = mul(R, cameraRelativeDirection);
    float x = center.x;
    float y = center.y;
    float z = center.z;
    float invZ = 1.0 / z;
    float invZ2 = invZ * invZ;

    // Corrected Jacobian using perspective projection ∂proj/∂pos
    float2x3 J;
    J[0] = float3(invZ, 0, -x * invZ2);
    J[1] = float3(0, invZ, -y * invZ2);

    // Step 4: Project 3D covariance to 2D
    float scaleZ = 1.0 / (z * z);
    float2x2 C_2D = mul(J, mul(C_view * scaleZ, transpose(J)));

    // Step 5: Eigen-decomposition to get ellipse axes
    float a = C_2D[0][0], b = C_2D[0][1];
    float c = C_2D[1][0], d = C_2D[1][1];
    float T = a + d;
    float D = a * d - b * c;
    float disc = sqrt(max(T*T - 4*D, 0));
    float lambda1 = 0.5 * (T + disc);
    float lambda2 = 0.5 * (T - disc);

    // Handle eigenvectors
    float2 majorAxis = normalize(float2(b, lambda1 - a));
    if (dot(majorAxis, majorAxis) < 1e-5)
        majorAxis = normalize(float2(lambda1 - d, c));
    float2 minorAxis = float2(majorAxis.y, -majorAxis.x);

    float2 scale = float2(sqrt(max(lambda1, 0)), sqrt(max(lambda2, 0)));
    float2x2 basis = float2x2(majorAxis, minorAxis);

    return mul(basis, float2x2(scale.x, 0, 0, scale.y));
}

I don't really know how it all works if I'm honest. I just need this magic sauce to work.

And rendered in unity using a plane mesh with UV matching plane, scaled down to zero in Blender, and using shader:

v2f vert(appdata v) {
    v2f o;
    float2 uv = float2(saturate(v.uv.x), v.uv.y);
    uv = uv - 0.5;
    o.uv = uv*2;
    
    float3x3 sphereMatRot = ...;
    float3x3 sphereMatScale = ...;
    
    float3 cameraTangentX = mul(UNITY_MATRIX_V[0].xyz, transpose((float3x3)unity_WorldToObject));
    float3 cameraTangentY = mul(UNITY_MATRIX_V[1].xyz, transpose((float3x3)unity_WorldToObject));
    float3 cameraDir = normalize(mul(unity_WorldToObject, float4(_WorldSpaceCameraPos.xyz, 1)) - v.vertex);
    float3x3 cameraRot = ExtractRotation(UNITY_MATRIX_V);
    
    // ===== PROJECT SPHERE TO ECLIPSE
    float2x2 ellipseTrans = Project3DGaussianTo2D(
        mul(sphereMatRot, sphereMatScale),
        cameraRot,
        cameraDir
    );
    uv = mul(ellipseTrans, uv.xy);
    
    float3 offset = uv.x * cameraTangentX + uv.y * cameraTangentY;
    o.pos = UnityObjectToClipPos(v.vertex + float4(offset, 0));
    return o;
}
float4 frag(v2f i) : SV_Target {
    float circle = saturate(1-length(i.uv));
    if (circle < 0.02) discard;
    return float4(0,0,1,1);
}
```