Browse code

Add implementation

Robert Cranston authored on 21/12/2021 12:41:46
Showing 3 changed files

... ...
@@ -2,11 +2,106 @@
2 2
 
3 3
 A [GLSL][]/[OpenGL][] [\>=2.1][] [isometry][] matrix library.
4 4
 
5
+If a `mat4` matrix represents an [isometry][], which in this case means that it
6
+encodes only rotations and translations (no scaling or shearing, and a
7
+projective `w` entry equal to `1`), it is easy and performant to extract those
8
+rotations and translations, and to compute the [inverse][]. This library
9
+provides functions to do so.
10
+
5 11
 [`gliso`]: https://git.rcrnstn.net/rcrnstn/gliso
6 12
 [GLSL]: https://en.wikipedia.org/wiki/OpenGL_Shading_Language
7 13
 [OpenGL]: https://en.wikipedia.org/wiki/OpenGL
8 14
 [\>=2.1]: https://en.wikipedia.org/wiki/OpenGL#Version_history
9 15
 [isometry]: https://en.wikipedia.org/wiki/Isometry
16
+[inverse]: https://en.wikipedia.org/wiki/Invertible_matrix
17
+
18
+## Requirements
19
+
20
+Support for `#include` directives is required. This can be provided by e.g. the
21
+standardized [`ARB_shading_language_include`][] extension or some third party
22
+library.
23
+
24
+[`ARB_shading_language_include`]: https://www.khronos.org/registry/OpenGL/extensions/ARB/ARB_shading_language_include.txt
25
+
26
+## Usage
27
+
28
+Link shader programs with the provided `iso.glsl`.
29
+
30
+Include the provided `iso.h` from shaders. It declares the following functions:
31
+
32
+-   `mat{3,4} isorot{inv}{3,4}(mat4 iso)`
33
+
34
+    Extract (inverse) rotation.
35
+
36
+-   `vec{3,4} isotrans{inv}{3,4}(mat4 iso)`
37
+
38
+    Extract (inverse) translation.
39
+
40
+    The `z` components of the `4` versions are set to `0`.
41
+
42
+-   `mat4 isoinv(mat4 iso)`
43
+
44
+    Compute inverse.
45
+
46
+## Theory
47
+
48
+For simplicity, we use the same name for both the [projective][] and
49
+[Euclidean][] versions of vectors and transformations below.
50
+
51
+Assume $M$ is an isometry,
52
+
53
+$$
54
+M v
55
+=
56
+\begin{pmatrix}
57
+    r_{xx} & r_{yx} & r_{zx} & t_x \\
58
+    r_{xx} & r_{yx} & r_{zx} & t_y \\
59
+    r_{xx} & r_{yx} & r_{zx} & t_z \\
60
+         0 &      0 &      0 &   1 \\
61
+\end{pmatrix}
62
+\begin{pmatrix}
63
+    v_x \\
64
+    v_y \\
65
+    v_z \\
66
+      1 \\
67
+\end{pmatrix}
68
+=
69
+R v + t
70
+$$
71
+
72
+It is trivial to extract the rotation $R$ and translation $t$ directly from the
73
+components of $M$.
74
+
75
+Since the rotation part $R$ is [orthonormal][], its inverse is its own
76
+transpose. The inverse of the translation is its own negation.
77
+
78
+$$
79
+\begin{aligned}
80
+    R^{-1} &= R^T \\
81
+    t^{-1} &= -t \\
82
+\end{aligned}
83
+$$
84
+
85
+Assuming the $3 \mathsf{x} 3$ sub-matrix of $M^{-1}$ is $R^{-1}$, finding the
86
+full inverse is simply a matter of finding $M^{-1}$'s fourth column $m^{-1}$:
87
+
88
+$$
89
+\begin{aligned}
90
+    v
91
+    &= M^{-1} M v \\
92
+    &= R^{-1} (R v + t) + m^{-1} \\
93
+    &= v + R^{-1} t + m^{-1} \\
94
+    &\iff \\
95
+    m^{-1} &= -R^{-1} t = R^{-1} t^{-1} \\
96
+\end{aligned}
97
+$$
98
+
99
+Computing this is much faster than the built-in, general, [`inverse`][].
100
+
101
+[projective]: https://en.wikipedia.org/wiki/Homogeneous_coordinates
102
+[Euclidean]: https://en.wikipedia.org/wiki/Cartesian_coordinates
103
+[orthonormal]: https://en.wikipedia.org/wiki/Orthogonal_matrix
104
+[`inverse`]: https://www.khronos.org/registry/OpenGL-Refpages/gl4/html/inverse.xhtml
10 105
 
11 106
 ## Build system
12 107
 
13 108
new file mode 100644
... ...
@@ -0,0 +1,47 @@
1
+#version 120
2
+
3
+
4
+mat3 isorot3     (mat4 iso) { return mat3(iso); }
5
+vec3 isotrans3   (mat4 iso) { return iso[3].xyz; }
6
+
7
+mat3 isorotinv3  (mat4 iso) { return transpose(isorot3(iso)); }
8
+vec3 isotransinv3(mat4 iso) { return -isotrans3(iso); }
9
+
10
+mat4 isorot4     (mat4 iso) { return mat4(isorot3     (iso)); }
11
+mat4 isorotinv4  (mat4 iso) { return mat4(isorotinv3  (iso)); }
12
+vec4 isotrans4   (mat4 iso) { return vec4(isotrans3   (iso), 0); }
13
+vec4 isotransinv4(mat4 iso) { return vec4(isotransinv3(iso), 0); }
14
+
15
+mat4 isoinv(mat4 iso)
16
+{
17
+    // mat4 inv = isorotinv4(iso);
18
+    // inv[3] += inv * isotransinv4(iso);
19
+    mat3 rotinv = isorotinv3(iso);
20
+    mat4 inv = mat4(rotinv);
21
+    inv[3].xyz = rotinv * isotransinv3(iso);
22
+    return inv;
23
+}
24
+
25
+
26
+#ifdef GLSLRUN
27
+void main()
28
+{
29
+    float phi = 1;
30
+    mat4 iso = transpose(mat4(
31
+        +cos(phi), -sin(phi), 0, 1,
32
+        +sin(phi), +cos(phi), 0, 2,
33
+        0,         0,         1, 3,
34
+        0,         0,         0, 1
35
+    ));
36
+    vec4 vec = vec4(4, 5, 6, 1);
37
+
38
+    cout
39
+        << iso * vec << endl
40
+        << isorot4(iso) *      vec  + isotrans4(iso) << endl
41
+        << isorot3(iso) * vec3(vec) + isotrans3(iso) << endl;
42
+
43
+    cout
44
+        << inverse(iso) << endl
45
+        << isoinv (iso) << endl;
46
+}
47
+#endif // GLSLRUN
0 48
new file mode 100644
... ...
@@ -0,0 +1,13 @@
1
+mat3 isorot3(mat4 iso);
2
+mat4 isorot4(mat4 iso);
3
+
4
+mat3 isorotinv3(mat4 iso);
5
+mat4 isorotinv4(mat4 iso);
6
+
7
+vec3 isotrans3(mat4 iso);
8
+vec4 isotrans4(mat4 iso);
9
+
10
+vec3 isotransinv3(mat4 iso);
11
+vec4 isotransinv4(mat4 iso);
12
+
13
+mat4 isoinv(mat4 iso);