Look, Ma:
No Matrices

AGACSE 2024

Steven De Keninck
Computer Vision Group
University of Amsterdam

Since the SIGGRAPH'19 course :

• Marc Ten Bosch : nD rigid body physics
• David Ruhe & co : GCAN, Clifford Group Equivariant Networks, ...
• Pim De Haan & co : Geometric Algebra Transformers
• Jeremy Ong : Klein Library
• Where's the simple 3D rendering stuff?

Look, Ma: No Matrices

https://github.com/enkimute/LookMaNoMatrices

A glTF compliant forward renderer in PGA \( \mathbb R_{3,0,1} \)!

Why No Matrices?

We love matrices.

We also love GA.

Speed vs Elegance?

PGA recap.

PGA recap.

Classic Matrix/Vector Approach

\(\stackrel {\text{transformations}} {\begin{pmatrix} & & \\ & & \\ & & \end{pmatrix}} \qquad \stackrel {\text{elements}} {\begin{pmatrix} \\ \\ \\ \end{pmatrix}} \)

PGA recap.

Geometric Algebra Approach

Model only the symmetry group.

Instead of matrices use multivectors (vectors, bivectors, ...)

PGA recap.

PGA models nD Geometry by modeling reflections in planes and their compositions

Orthogonal reflections, called blades, double as the elements of geometry.

PGA recap.

PGA models nD Geometry by modeling reflections in planes and their compositions

\[ \mathbb R_{3,0,1} \]

\({\mathbf e_1}^2 = 1\), \({\mathbf e_2}^2 = 1\), \({\mathbf e_3}^2 = 1\), \({\mathbf e_0}^2 = 0\)

plane(-reflections) as vectors. \( ax + by + cz + d = 0 \Leftrightarrow a\mathbf e_1 + b\mathbf e_2 + c\mathbf e_3 + d\mathbf e_0 \)

line(-reflections) as bivectors.

point(-reflections) as trivectors. \( [x,y,z]^\top \Leftrightarrow x\mathbf e_1^* + y\mathbf e_2^* + z\mathbf e_3^* + \mathbf e_0^* \)

PGA recap.

PGA models nD Geometry by modeling reflections in planes and their compositions

PGA recap.

PGA models nD Geometry by modeling reflections in planes and their compositions

PGA recap.

PGA models nD Geometry by modeling reflections in planes and their compositions

PGA recap.

PGA models nD Geometry by modeling reflections in planes and their compositions

PGA recap.

PGA models nD Geometry by modeling reflections in planes and their compositions

PGA recap.

PGA models nD Geometry by modeling reflections in planes and their compositions

Universal coordinate, dimension and basis agnostic expressions that 'just work' for any point, line, plane, ...

FPGA - Fast PGA?

The composition motors - aka dual quaternions.

function gp_mm (a,b,res) {
    const a0=a[0],a1=a[1],a2=a[2],a3=a[3],a4=a[4],a5=a[5],a6=a[6],a7=a[7],
            b0=b[0],b1=b[1],b2=b[2],b3=b[3],b4=b[4],b5=b[5],b6=b[6],b7=b[7];
    res[0] = a0*b0-a1*b1-a2*b2-a3*b3;
    res[1] = a0*b1+a1*b0+a3*b2-a2*b3;
    res[2] = a0*b2+a1*b3+a2*b0-a3*b1;
    res[3] = a0*b3+a2*b1+a3*b0-a1*b2;
    res[4] = a0*b4+a3*b5+a4*b0+a6*b2-a1*b7-a2*b6-a5*b3-a7*b1;
    res[5] = a0*b5+a1*b6+a4*b3+a5*b0-a2*b7-a3*b4-a6*b1-a7*b2;
    res[6] = a0*b6+a2*b4+a5*b1+a6*b0-a1*b5-a3*b7-a4*b2-a7*b3;
    res[7] = a0*b7+a1*b4+a2*b5+a3*b6+a4*b1+a5*b2+a6*b3+a7*b0;
    return res;
}

48 muls and 40 adds

FPGA - Fast PGA?

The composition motors - aka dual quaternions.

motor gp_mm( motor a, motor b ) {
  return motor(
    a[0].x*b[0].x   - dot(a[0].yzw, b[0].yzw), 
    a[0].x*b[0].yzw + b[0].x*a[0].yzw + cross(b[0].yzw, a[0].yzw),
    a[0].x*b[1].xyz + b[0].x*a[1].xyz + cross(b[0].yzw, a[1].xyz) 
    + cross(b[1].xyz, a[0].yzw) - b[1].w*a[0].yzw - a[1].w*b[0].yzw, 
    a[0].x*b[1].w + b[0].x*a[1].w 
    + dot(a[0].yzw, b[1].xyz) + dot(a[1].xyz, b[0].yzw));
}

48 muls and 40 adds

FPGA - Fast PGA?

The composition motors - aka dual quaternions.

FPGA - Fast PGA?

// 21 mul, 18 add
point sw_mp( motor a, point b ) {
    direction t = cross(b, a[0].yzw)  - a[1].xyz;
    return  (a[0].x * t + cross(t, a[0].yzw) - a[0].yzw * a[1].w) * 2. + b;  
}

FPGA - Fast PGA?

// 21 mul, 18 add
point sw_mp( motor a, point b ) {
    direction t = cross(b, a[0].yzw)  - a[1].xyz;
    return  (a[0].x * t + cross(t, a[0].yzw) - a[0].yzw * a[1].w) * 2. + b;  
}

// 18 mul, 12 add
direction sw_md( motor a, direction b ) {
    direction t = cross(b, a[0].yzw);
    return  (a[0].x * t + cross(t, a[0].yzw)) * 2. + b;                      
}

recall, transform points with a 4x4 matrix is 16/12, and with 3x4 it is 12/9.

FPGA - Fast PGA?

// 6 mul, 4 add
direction sw_mx( motor a ) {
    return direction(
        0.5 - a[0].w*a[0].w - a[0].z*a[0].z,                                 
        a[0].z*a[0].y - a[0].x*a[0].w, 
        a[0].w*a[0].y + a[0].x*a[0].z
    );
}

That's cheap.

FPGA - Fast PGA?

Many more useful and fun bits for various PGA functions in GLSL are in the writeup and code.

• normalization
• exponential/logarithm
• inverses
• factorizations
• conversion to/from matrices

Matrix-free Forward Rendering

• PGA toolbox on CPU
  
  • Convert matrices → motors
  • Handle scaling
• PGA toolbox on GPU

  #define motor mat2x4   // [[s,e23,e31,e12],[e01,e02,e03,e0123]] 
  #define line  mat2x3   // [[e23,e31,e12],[e01,e02,e03]]
  #define point vec3     // [e032,e013,e021 ] 1 e123
  #define direction vec3 // [e032,e013,e021 ] 0 e123

• Model and View motors
• Projection function

Tangent Space Normalmapping

image from wikimedia.

• Classic setup : use normal and tangent to construct orthonormal TBN matrix in FS.
• Storage cost : normal + tangent + handedness
• Use PGA versors instead?

Tangent Space Normalmapping

The naive approach of simply converting from matrices to motors.

normal_out      = normal_matrix * normal_in;                                 
tangent_out.xyz = model_matrix * tangent_in.xyz;
tangent_out.w   = tangent_in.w;

normal_out      = sw_md( toWorld, normal_in );                               
tangent_out.xyz = sw_md( toWorld, tangent_in );
tangent_out.w   = tangent_in.w;

Tangent Space Normalmapping

The naive approach of simply converting from matrices to motors.

normal_out      = normal_matrix * normal_in;                                 
tangent_out.xyz = normal_matrix * tangent_in.xyz;
tangent_out.w   = tangent_in.w;

normal_out      = sw_md( toWorld, normal_in );                               
tangent_out.xyz = sw_md( toWorld, tangent_in );
tangent_out.w   = tangent_in.w;

18/12 \(\rightarrow\) 36/24

Please don't do this.

Tangent Rotors and Double Covers

• PGA \(\mathbb R_{3,0,1}\) encodes all k-reflections.
• Limit to origin : \( \mathbb R_{3,0,1} \rightarrow \mathbb R_3 \).
• Limit to 2-reflections : \( \mathbb R_{3} \rightarrow \mathbb R_3^+ \) + handedness.
• Use double cover to store handedness.

Tangent Rotors and Double Covers

Option 1 : Extract normal and tangent in vertexshader.

• From \(3 + 3 + 1\) floats to just \(4\).
• Compose tangent rotor with model rotor. (16/12).
• Extract tangent + normal in vertex shader. (9/8).
• Compare to transform normal and tangent. (18/12)
• No changes to the fragment shader.

Tangent Rotors and Double Covers

Option 1 : Extract normal and tangent in vertexshader.

// 9 muls, 8 adds
void extractNormalTangent(motor a,out direction normal,out direction tangent){
  float yw = a[0].y * a[0].w;
  float xz = a[0].x * a[0].z;
  float zz = a[0].z * a[0].z;
  normal = direction(yw-xz, a[0].z*a[0].w+a[0].y*a[0].x, 0.5-zz-a[0].y*a[0].y);
  tangent= direction(0.5-zz-a[0].w*a[0].w, a[0].z*a[0].y-a[0].x*a[0].w, yw+xz);
}

Tangent Rotors and Double Covers

Option 1 : Extract normal and tangent in vertexshader.

normal_out      = normal_matrix * normal_in;                                 
tangent_out.xyz = model_matrix * tangent_in.xyz;
tangent_out.w   = tangent_in.w;

motor tangentRotor = gp_rr( toWorld, motor(attrib_tangentRotor,vec4(0.)) );  
extractNormalTangent(tangentRotor, normal_out, tangent_out.xyz);
tangent_out.w = sign(1.0 / attrib_tangentRotor.x);

Tangent Rotors and Double Covers

Option 1 : Extract normal and tangent in vertexshader.

normal_out      = normal_matrix * normal_in;                                 
tangent_out.xyz = model_matrix * tangent_in.xyz;
tangent_out.w   = tangent_in.w;

motor tangentRotor = gp_rr( toWorld, motor(attrib_tangentRotor,vec4(0.)) );  
extractNormalTangent(tangentRotor, normal_out, tangent_out.xyz);
tangent_out.w = sign(1.0 / attrib_tangentRotor.x);

18/12 \(\rightarrow\) 25/18

Tangent Rotors and Double Covers

At this point in the project, discord user Criver pointed me towards a (video only) reference from the crytek team that me and the SIGGRAPH reviewers had missed.

"Spherical skinning with dual quaternions and QTangents" by Frey & Herzeg.

Luckily, their last slide contained some topics for future research ...

Tangent Rotors and Double Covers

Option 2 : Tangent Rotors in Fragment Shader.

• From \(7\) floats to \(4\) + one less varying.
• Compose tangent and model rotor. (16/12).
• Transform sampled normal with rotor, renorm. (22/15)
• Compare to transform in VS (18/12), re-orthonormalize + transform in FS. (33/18)

Tangent Rotors and Double Covers

Option 2 : Tangent Rotors in Fragment Shader.

• problem : incompatible tangent rotors.

Tangent Rotors and Double Covers

Option 2 : Tangent Rotors in Fragment Shader.

• problem : incompatible tangent rotors.
• solution : subdivide edges.

Summary of results

Storage and Compute requirements for the tangent frame.

*note : data just for tangent frame. for non-skinned vertex position transform, 4x3 matrices win. A hybrid solution that combines this with PGA tangent transform should be considered.

Thank you for your attention!

Join us on the https://bivector.net discord!