Browse code

Add transformations

Robert Cranston authored on 15/07/2024 00:23:13
Showing 3 changed files

1 1
new file mode 100644
2 2
Binary files /dev/null and b/doc/03_transformations.tga differ
3 3
new file mode 100755
... ...
@@ -0,0 +1,153 @@
1
+#!/usr/bin/env python3
2
+
3
+from sympy import *
4
+
5
+
6
+def mat(*args):
7
+    return Matrix(list(map(list, args))).T
8
+
9
+def vec(*args):
10
+    return Matrix(args)
11
+
12
+def normalize(v):
13
+    return v / sqrt(v.dot(v))
14
+
15
+def elementwise(a, b, f):
16
+    return Matrix([
17
+        [
18
+            f(a[r, c], b[r, c])
19
+            for c in range(a.cols)
20
+        ]
21
+        for r in range(a.rows)
22
+    ])
23
+
24
+def mul(a, b):
25
+    return elementwise(a, b, lambda a, b: a * b)
26
+
27
+def div(a, b):
28
+    return elementwise(a, b, lambda a, b: a / b)
29
+
30
+
31
+def window(viewport, depth_range):
32
+    size = vec(*viewport[2:4], depth_range[1]-depth_range[0])
33
+    offs = vec(*viewport[0:2],                depth_range[0])
34
+    s = 0.5 * size
35
+    t = 0.5 * size + offs
36
+    return Matrix([
37
+        [s[0], 0,    0,    t[0]],
38
+        [0,    s[1], 0,    t[1]],
39
+        [0,    0,    s[2], t[2]],
40
+        [0,    0,    0,    1],
41
+    ])
42
+
43
+def ortho(left, right, bottom, top, near, far):
44
+    size = vec(right-left, top-bottom, far-near)
45
+    offs = vec(     -left,    -bottom,    +near)
46
+    s = div(vec(2, 2, -2), size)
47
+    t = mul(s, offs) - vec(1, 1, 1)
48
+    return Matrix([
49
+        [s[0], 0,    0,    t[0]],
50
+        [0,    s[1], 0,    t[1]],
51
+        [0,    0,    s[2], t[2]],
52
+        [0,    0,    0,    1],
53
+    ])
54
+
55
+def frustum(left, right, bottom, top, near, far):
56
+    z = vec(0, 0, near+far, -1)
57
+    w = vec(0, 0, near*far,  0)
58
+    return ortho(left, right, bottom, top, near, far) * Matrix([
59
+        [near, 0,    z[0], z[0]],
60
+        [0,    near, z[1], w[1]],
61
+        [0,    0,    z[2], w[2]],
62
+        [0,    0,    z[3], w[3]],
63
+    ])
64
+
65
+def perspective(fovy, aspect, near, far):
66
+    y = near * tan(0.5 * fovy)
67
+    x = y * aspect
68
+    return frustum(-x, +x, -y, +y, near, far)
69
+
70
+def lookat(position, target, up):
71
+    z = normalize(position - target)
72
+    x = normalize(up.cross(z))
73
+    y = z.cross(x)
74
+    R_inv = mat(x, y, z).T
75
+    t_inv = -position
76
+    return mat(R_inv.row(0), R_inv.row(0), R_inv.row(0), R_inv * t_inv)
77
+
78
+
79
+# https://registry.khronos.org/OpenGL-Refpages/gl4/html/glViewport.xhtml
80
+# https://registry.khronos.org/OpenGL-Refpages/gl4/html/glDepthRange.xhtml
81
+# https://www.khronos.org/opengl/wiki/Vertex_Post-Processing#Viewport_transform
82
+# https://www.songho.ca/opengl/gl_viewport.html
83
+viewport    = symbols('x, y, width, height')
84
+depth_range = symbols('near, far')
85
+pprint(simplify(window(viewport, depth_range)))
86
+# ⎡0.5⋅width      0               0             0.5⋅width + x   ⎤
87
+# ⎢                                                             ⎥
88
+# ⎢    0      0.5⋅height          0             0.5⋅height + y  ⎥
89
+# ⎢                                                             ⎥
90
+# ⎢    0          0       0.5⋅far - 0.5⋅near  0.5⋅far + 0.5⋅near⎥
91
+# ⎢                                                             ⎥
92
+# ⎣    0          0               0                   1         ⎦
93
+
94
+# https://registry.khronos.org/OpenGL-Refpages/gl2.1/xhtml/glOrtho.xml
95
+# https://www.songho.ca/opengl/gl_projectionmatrix.html#ortho
96
+left, right, bottom, top, near, far = symbols('left, right, bottom, top, near, far')
97
+pprint(simplify(ortho(left, right, bottom, top, near, far)))
98
+# ⎡    -2                                  left + right⎤
99
+# ⎢────────────       0            0       ────────────⎥
100
+# ⎢left - right                            left - right⎥
101
+# ⎢                                                    ⎥
102
+# ⎢                  -2                    bottom + top⎥
103
+# ⎢     0        ────────────      0       ────────────⎥
104
+# ⎢              bottom - top              bottom - top⎥
105
+# ⎢                                                    ⎥
106
+# ⎢                               -2       -far - near ⎥
107
+# ⎢     0             0        ──────────  ─────────── ⎥
108
+# ⎢                            far - near   far - near ⎥
109
+# ⎢                                                    ⎥
110
+# ⎣     0             0            0            1      ⎦
111
+
112
+# https://registry.khronos.org/OpenGL-Refpages/gl2.1/xhtml/glFrustum.xml
113
+# https://www.songho.ca/opengl/gl_projectionmatrix.html#perspective
114
+left, right, bottom, top, near, far = symbols('left, right, bottom, top, near, far')
115
+pprint(simplify(frustum(left, right, bottom, top, near, far)))
116
+# ⎡  -2⋅near                   -left - right              ⎤
117
+# ⎢────────────       0        ─────────────       0      ⎥
118
+# ⎢left - right                 left - right              ⎥
119
+# ⎢                                                       ⎥
120
+# ⎢                -2⋅near     -bottom - top              ⎥
121
+# ⎢     0        ────────────  ─────────────       0      ⎥
122
+# ⎢              bottom - top   bottom - top              ⎥
123
+# ⎢                                                       ⎥
124
+# ⎢                             -far - near   -2⋅far⋅near ⎥
125
+# ⎢     0             0         ───────────   ────────────⎥
126
+# ⎢                              far - near    far - near ⎥
127
+# ⎢                                                       ⎥
128
+# ⎣     0             0             -1             0      ⎦
129
+
130
+# https://registry.khronos.org/OpenGL-Refpages/gl2.1/xhtml/gluPerspective.xml
131
+# https://www.songho.ca/opengl/gl_projectionmatrix.html#fov
132
+fovy, aspect, near, far = symbols('fovy, aspect, near, far')
133
+pprint(simplify(perspective(fovy, aspect, near, far)))
134
+# ⎡         1                                                    ⎤
135
+# ⎢────────────────────        0             0            0      ⎥
136
+# ⎢aspect⋅tan(0.5⋅fovy)                                          ⎥
137
+# ⎢                                                              ⎥
138
+# ⎢                            1                                 ⎥
139
+# ⎢         0            ─────────────       0            0      ⎥
140
+# ⎢                      tan(0.5⋅fovy)                           ⎥
141
+# ⎢                                                              ⎥
142
+# ⎢                                     -far - near  -2⋅far⋅near ⎥
143
+# ⎢         0                  0        ───────────  ────────────⎥
144
+# ⎢                                      far - near   far - near ⎥
145
+# ⎢                                                              ⎥
146
+# ⎣         0                  0            -1            0      ⎦
147
+
148
+# https://registry.khronos.org/OpenGL-Refpages/gl2.1/xhtml/gluLookAt.xml
149
+# https://www.songho.ca/opengl/gl_camera.html#lookat
150
+position = Matrix(MatrixSymbol('position', 3, 1))
151
+target   = Matrix(MatrixSymbol('target',   3, 1))
152
+up       = Matrix(MatrixSymbol('up',       3, 1))
153
+# pprint(simplify(lookat(position, target, up)))
... ...
@@ -85,17 +85,24 @@ struct Scene
85 85
 //// Camera
86 86
 struct Camera
87 87
 {
88
-    float near;
89
-    float far;
90
-    uvec2 size;
91
-    Ray ray(vec4 const & frag_coord) const
88
+    mat4 camera_inverse;
89
+    explicit Camera(
90
+        mat4 const & view_,
91
+        mat4 const & projection_,
92
+        mat4 const & window_
93
+    )
94
+    :
95
+        camera_inverse{inverse(window_ * projection_ * view_)}
96
+    {}
97
+    vec3 unproject(vec3 const & frag_coord) const
92 98
     {
93
-        auto const point = vec3(
94
-            (2.0F * vec2(frag_coord.xy) - vec2(size)) / float(size.y),
95
-            -1.0F
96
-        );
97
-        auto const origin    = near * point;
98
-        auto const stop      = far  * point;
99
+            auto const world_coord = camera_inverse * vec4(frag_coord, 1.0F);
100
+            return vec3(world_coord) / world_coord.w;
101
+    }
102
+    Ray ray(vec4 const & frag_coord, vec2 const & depth_range) const
103
+    {
104
+        auto const origin    = unproject(vec3(frag_coord.xy, depth_range[0]));
105
+        auto const stop      = unproject(vec3(frag_coord.xy, depth_range[1]));
99 106
         auto const range     = length(stop - origin);
100 107
         auto const direction = (stop - origin) / range;
101 108
         return Ray{origin, direction, range};
... ...
@@ -113,15 +120,16 @@ struct Framebuffer
113 120
         color((std::size_t)(size_.x * size_.y))
114 121
     {}
115 122
     template<typename Shader>
116
-    void render(Shader const & shader)
123
+    void render(uvec4 viewport, vec2 depth_range, Shader const & shader)
117 124
     {
118
-        auto const begin = uvec2(0);
119
-        auto const end   = size;
125
+        auto const begin = max(uvec2(0),    uvec2(viewport.xy));
126
+        auto const end   = min(uvec2(size), uvec2(viewport.xy + viewport.zw));
127
+        auto const z     = 0.5F * (depth_range[0] + depth_range[1]);
120 128
         for (auto y = begin.y; y < end.y; ++y)
121 129
         for (auto x = begin.x; x < end.x; ++x)
122 130
         {
123 131
             auto const index      = size.x * y + x;
124
-            auto const frag_coord = vec4(vec2(x, y) + 0.5F, 0.0F, 1.0F);
132
+            auto const frag_coord = vec4(vec2(x, y) + 0.5F, z, 1.0F);
125 133
             auto const frag_color = vec4(shader(frag_coord).bgra);
126 134
             color[index] = clamp(frag_color, 0.0F, 1.0F) * 255.0F + 0.5F;
127 135
         }
... ...
@@ -150,6 +158,58 @@ struct Framebuffer
150 158
     }
151 159
 };
152 160
 
161
+/// Transformations
162
+
163
+//// window
164
+mat4 window(uvec4 viewport, vec2 depth_range)
165
+{
166
+    auto const size = vec3(uvec2(viewport.zw), depth_range[1]-depth_range[0]);
167
+    auto const offs = vec3(uvec2(viewport.xy),                depth_range[0]);
168
+    auto const s = 0.5F * size;
169
+    auto const t = 0.5F * size + offs;
170
+    auto const S = matrixCompMult(mat3(1.0F), outerProduct(vec3(1.0F), s));
171
+    return mat4(mat4x3(S[0], S[1], S[2], t));
172
+}
173
+
174
+//// ortho
175
+mat4 ortho(float left, float right, float bottom, float top, float near, float far)
176
+{
177
+    auto const size = vec3(right-left, top-bottom, far-near);
178
+    auto const offs = vec3(     -left,    -bottom,    +near);
179
+    auto const s = vec3(2.0F, 2.0F, -2.0F) / size;
180
+    auto const t = s * offs - 1.0F;
181
+    auto const S = matrixCompMult(mat3(1.0F), outerProduct(vec3(1.0F), s));
182
+    return mat4(mat4x3(S[0], S[1], S[2], t));
183
+}
184
+
185
+//// frustum
186
+mat4 frustum(float left, float right, float bottom, float top, float near, float far)
187
+{
188
+    auto const z = vec4(0.0F, 0.0F, near+far, -1.0F);
189
+    auto const w = vec4(0.0F, 0.0F, near*far,  0.0F);
190
+    auto const S = mat4(near);
191
+    return ortho(left, right, bottom, top, near, far) * mat4(S[0], S[1], z, w);
192
+}
193
+
194
+//// perspective
195
+mat4 perspective(float fovy, float aspect, float near, float far)
196
+{
197
+    auto const y = near * tan(0.5F * radians(fovy));
198
+    auto const x = y * aspect;
199
+    return frustum(-x, +x, -y, +y, near, far);
200
+}
201
+
202
+//// lookat
203
+mat4 lookat(vec3 position, vec3 target, vec3 up)
204
+{
205
+    auto const z = normalize(position - target);
206
+    auto const x = normalize(cross(up, z));
207
+    auto const y = cross(z, x);
208
+    auto const R_inv = transpose(mat3(x, y, z));
209
+    auto const t_inv = -position;
210
+    return mat4(mat4x3(R_inv[0], R_inv[1], R_inv[2], R_inv * t_inv));
211
+}
212
+
153 213
 /// Main
154 214
 
155 215
 int main(int argc, char const * argv[])
... ...
@@ -163,12 +223,24 @@ int main(int argc, char const * argv[])
163 223
     auto const * path = argv[1]; // NOLINT
164 224
     //// Configure
165 225
     auto size        = uvec2(640, 480);
226
+    auto aspect      = (float)size.x / (float)size.y;
227
+    auto viewport    = uvec4(uvec2(0), size);
228
+    auto depth_range = vec2(0.0F, 1.0F);
166 229
     auto framebuffer = Framebuffer(size);
167
-    auto camera = Camera{
168
-        1.0F,   // near
169
-        100.0F, // far
170
-        size,   // size
171
-    };
230
+    auto camera = Camera(
231
+        lookat( // view
232
+            {0.0F, 0.0F,  0.0F}, // position
233
+            {0.0F, 0.0F, -1.0F}, // target
234
+            {0.0F, 1.0F,  0.0F}  // up
235
+        ),
236
+        perspective( // projection
237
+            90.0F,  // fovy
238
+            aspect, // aspect
239
+            1.0F,   // near
240
+            100.0F  // far
241
+        ),
242
+        window(viewport, depth_range)
243
+    );
172 244
     auto scene = Scene{
173 245
         { // shapes
174 246
             Sphere{{   6.0F,  -4.0F,  -24.0F}, 12.0F},
... ...
@@ -180,9 +252,9 @@ int main(int argc, char const * argv[])
180 252
         },
181 253
     };
182 254
     //// Render
183
-    framebuffer.render([&](vec4 const & frag_coord)
255
+    framebuffer.render(viewport, depth_range, [&](vec4 const & frag_coord)
184 256
     {
185
-        auto const ray   = camera.ray(frag_coord);
257
+        auto const ray   = camera.ray(frag_coord, depth_range);
186 258
         auto const trace = scene.trace(ray);
187 259
         auto const rgb   = vec3(1.0F / (1.0F + trace.distance));
188 260
         auto const alpha = (bool)trace;