a49530da16f0 — Leonard Ritter a month ago
* checked in autoinverse
3 files changed, 229 insertions(+), 47 deletions(-)

M testing/bidimat.sc => lib/tukan/autoinverse.sc
A => testing/bidimat.cpp
A => testing/test_autoinverse.sc
M testing/bidimat.sc => lib/tukan/autoinverse.sc +114 -47
@@ 1,77 1,144 @@ 
 
 
+# auto-inverting 4x4 matrix, which transparently maintains its inverse through
+    all constructive and compositional operations.
+
 using import struct
 using import glm
-using import glsl
 
-import ..tukan.use
-using import tukan.rotation
-
-struct mat4bidi plain
-    right : mat4
-    left : mat4
+struct aimat4 plain
+    M : mat4
+    IM : mat4
 
     @@ memo
     inline __* (cls T)
         static-if (cls == T)
             fn (self other)
                 this-type
-                    self.right * other.right
-                    other.left * self.left
+                    self.M * other.M
+                    other.IM * self.IM
+        elseif (T == mat4.ColumnType)
+            inline (self other)
+                self.M * other
 
-    inline... scale
-    case (s : f32,)
-        this-function (vec3 s)
-    case (s : vec3,)
-        is := (/ s)
+    """"builds a matrix equivalent to the operations offset * rotate * scale
+    inline... transform (offset : vec3, rotate : mat3, scale : vec3)
+        iscale := (/ scale); ioffset := -offset
+        let rotx roty rotz = (unpack rotate)
+        let irotx iroty irotz =
+            rotx * iscale.x; roty * iscale.y; rotz * iscale.z
         this-type
-            mat4  s.x 0 0 0
-                \ 0 s.y 0 0
-                \ 0 0 s.z 0
-                \ 0 0 0 1
-            mat4  is.x 0 0 0
-                \ 0 is.y 0 0
-                \ 0 0 is.z 0
-                \ 0 0 0 1
-
-    inline... offset
-    case (o : vec3,)
-        io := -o
+            mat4  (rotx * scale.x) 0
+                \ (roty * scale.y) 0
+                \ (rotz * scale.z) 0
+                \ offset 1
+            mat4  irotx.x iroty.x irotz.x 0
+                \ irotx.y iroty.y irotz.y 0
+                \ irotx.z iroty.z irotz.z 0
+                \ (ioffset * iscale) 1
+    case (offset : vec3)
+        ioffset := -offset
         this-type
             mat4  1 0 0 0
                 \ 0 1 0 0
                 \ 0 0 1 0
-                \ o.x o.y o.z 1
+                \ offset 1
             mat4  1 0 0 0
                 \ 0 1 0 0
                 \ 0 0 1 0
-                \ io.x io.y io.z 1
+                \ ioffset 1
+    case (rotate : mat3)
+        irotate := (transpose rotate)
+        this-type
+            mat4 rotate
+            mat4 irotate
+    case (offset : vec3, rotate : mat3)
+        ioffset := -offset
+        let rotx roty rotz = (unpack rotate)
+        this-type
+            mat4  rotx 0
+                \ roty 0
+                \ rotz 0
+                \ offset 1
+            mat4  rotx.x roty.x rotz.x 0
+                \ rotx.y roty.y rotz.y 0
+                \ rotx.z roty.z rotz.z 0
+                \ ioffset 1
+    case (offset : vec3, scale : vec3)
+        iscale := (/ scale); ioffset := -offset
+        this-type
+            mat4  scale.x 0 0 0
+                \ 0 scale.y 0 0
+                \ 0 0 scale.z 0
+                \ offset 1
+            mat4  iscale.x 0 0 0
+                \ 0 iscale.y 0 0
+                \ 0 0 iscale.z 0
+                \ (ioffset * iscale) 1
+    case (rotate : mat3, scale : vec3)
+        iscale := (/ scale)
+        let rotx roty rotz = (unpack rotate)
+        let irotx iroty irotz =
+            rotx * iscale.x; roty * iscale.y; rotz * iscale.z
+        this-type
+            mat4  (rotx * scale.x) 0
+                \ (roty * scale.y) 0
+                \ (rotz * scale.z) 0
+                \ 0 0 0 1
+            mat4  irotx.x iroty.x irotz.x 0
+                \ irotx.y iroty.y irotz.y 0
+                \ irotx.z iroty.z irotz.z 0
+                \ 0 0 0 1
 
-    inline... affine
-    case (m : mat3,)
-        im := (transpose m)
+    # infinite far plane perspective transform
+        for this to work, the depth buffer must be reversed;
+        glDepthFunc must be set to GL_GREATER, glClearDepth set to 0,
+        and glDepthRanged set to -1, 1;
+    inline... perspective (aspect : vec2, near : f32)
+        iaspect := (/ aspect)
         this-type
-            mat4 m
-            mat4 im
+            mat4  aspect.x 0 0 0
+                \ 0 aspect.y 0 0
+                \ 0 0 0 1
+                \ 0 0 near 0
+            mat4  iaspect.x 0 0 0
+                \ 0 iaspect.y 0 0
+                \ 0 0 0 (/ near)
+                \ 0 0 1 0
+
+    # linear orthographic with near = 1 far = 0
+        for this to work, the depth buffer must be reversed;
+        glDepthFunc must be set to GL_GREATER, glClearDepth set to 0,
+        and glDepthRanged set to -1, 1;
+    inline... orthographic (aspect : vec2, near : f32, far : f32)
+        iaspect := (/ aspect)
+        izd := near - far
+        zd := (/ izd)
+        this-type
+            mat4  aspect.x 0 0 0
+                \ 0 aspect.y 0 0
+                \ 0 0 zd 0
+                \ 0 0 (-far * zd) 1
+            mat4  iaspect.x 0 0 0
+                \ 0 iaspect.y 0 0
+                \ 0 0 izd 0
+                \ 0 0 far 1
 
     fn complement (self)
-        self.right * self.left
+        """"should always be approximate to identity matrix
+        self.M * self.IM
+
+    inline inverse (self)
+        this-type self.IM self.M
 
     fn __repr (self)
-        .. "<mat4bidi "
-            repr self.right
+        .. "<aimat4 "
+            repr self.M
             " "
-            repr self.left
+            repr self.IM
             ">"
 
-m :=
-    mat4bidi.offset (vec3 1 2 3)
 
-print
-    mat4bidi.affine (mat3 5)
-
-print
-    m.right * (vec4 0 0 0 1)
-
-
-;
  No newline at end of file
+do
+    let aimat4
+    locals;
  No newline at end of file

          
A => testing/bidimat.cpp +71 -0
@@ 0,0 1,71 @@ 
+// gcc bidimat.cpp -o bidimat -lstdc++ -lm && ./bidimat
+
+#include <stdio.h>
+#include <math.h>
+#include <glm/glm.hpp>
+#include <glm/gtc/quaternion.hpp>
+#include <glm/gtx/quaternion.hpp>
+
+using namespace glm;
+
+void print(const vec3 &v) {
+    printf("%f %f %f\n\n", v[0], v[1], v[2]);
+}
+
+void print(const vec4 &v) {
+    printf("%f %f %f %f\n\n", v[0], v[1], v[2], v[3]);
+}
+
+void print(const quat &q) {
+    printf("%f %f %f %f\n\n", q[0], q[1], q[2], q[3]);
+}
+
+void print(const mat3 &m) {
+    printf("%f %f %f\n%f %f %f\n%f %f %f\n\n",
+        m[0][0], m[0][1], m[0][2],
+        m[1][0], m[1][1], m[1][2],
+        m[2][0], m[2][1], m[2][2]);
+}
+
+void print(const mat4 &m) {
+    printf("%f %f %f %f\n%f %f %f %f\n%f %f %f %f\n%f %f %f %f\n\n",
+        m[0][0], m[0][1], m[0][2], m[0][3],
+        m[1][0], m[1][1], m[1][2], m[1][3],
+        m[2][0], m[2][1], m[2][2], m[2][3],
+        m[3][0], m[3][1], m[3][2], m[3][3]);
+}
+
+void print(const mat4x3 &m) {
+    printf("%f %f %f\n%f %f %f\n%f %f %f\n%f %f %f\n\n",
+        m[0][0], m[0][1], m[0][2],
+        m[1][0], m[1][1], m[1][2],
+        m[2][0], m[2][1], m[2][2],
+        m[3][0], m[3][1], m[3][2]);
+}
+
+void print(const mat3x4 &m) {
+    printf("%f %f %f %f\n%f %f %f %f\n%f %f %f %f\n\n",
+        m[0][0], m[0][1], m[0][2], m[0][3],
+        m[1][0], m[1][1], m[1][2], m[1][3],
+        m[2][0], m[2][1], m[2][2], m[2][3]);
+}
+
+int main(int argc, char **argv) {
+    mat4 t(
+        vec4(1,0,0,0),
+        vec4(0,1,0,0),
+        vec4(0,0,1,0),
+        vec4(1,2,3,1));
+    quat v = angleAxis(radians(37.0f), normalize(vec3(1,2,3)));
+    print(v);
+    mat4 r = toMat4(v);
+    print(r);
+    mat4 s(
+        vec4(2,0,0,0),
+        vec4(0,3,0,0),
+        vec4(0,0,4,0),
+        vec4(0,0,0,1));
+    print(t * r * s);
+    print(inverse(t * r * s));
+    return 0;
+}

          
A => testing/test_autoinverse.sc +44 -0
@@ 0,0 1,44 @@ 
+
+
+using import glm
+
+import ..lib.tukan.use
+
+using import tukan.autoinverse
+using import tukan.rotation
+
+t :=
+    aimat4.transform (vec3 1 2 3)
+v :=
+    versor (normalize (vec3 1 2 3)) (radians 37.0)
+r :=
+    aimat4.transform (rotation v)
+s :=
+    aimat4.transform (vec3 0) (vec3 2 3 4)
+A := t * r * s
+B :=
+    aimat4.transform
+        vec3 1 2 3
+        rotation v
+        vec3 2 3 4
+print "t * r * s =" A
+print "transform(t,r,s) =" B
+print (A * ('inverse B)) # identity
+
+M :=
+    aimat4.perspective (vec2 2 1.2) 0.1
+    #aimat4.orthographic (vec2 2 1.2) 50.0 150.0
+
+print M
+print
+    'complement M
+
+q := (vec4 1.5 2.0 100.0 1.0)
+q2 := M * q
+q3 := q2 / q2.w
+q4 := ('inverse M) * q3
+q5 := q4 / q4.w
+
+print q q3 q5
+print (dot (q - q5) (q - q5))
+;
  No newline at end of file