make points and vectors differentiable

add rotational transformations (x, y, and z)
1 files changed, 244 insertions(+), 88 deletions(-)

M tracer.lisp
M tracer.lisp +244 -88
@@ 14,53 14,60 @@ 
   (< (abs (- f1 f2)) single-float-epsilon))
 
 ;;; code:
-(defmacro defun-vec3-binop (binop)
-  `(defun ,(intern (concatenate 'string "VEC3" (string binop))) (v1 v2)
+(defmacro defun-vector-binop (binop)
+  `(defun ,(intern (concatenate 'string "VECTOR" (string binop))) (v1 v2)
      (vector (,binop (aref v1 0) (aref v2 0))
              (,binop (aref v1 1) (aref v2 1))
-             (,binop (aref v1 2) (aref v2 2)))))
+             (,binop (aref v1 2) (aref v2 2))
+             (,binop (aref v1 3) (aref v2 3)))))
 
-(defun-vec3-binop +)
-(defun-vec3-binop -)
-(defun-vec3-binop *)
-(defun-vec3-binop /)
+(defun-vector-binop +)
+(defun-vector-binop -)
+(defun-vector-binop *)
+(defun-vector-binop /)
 
-(defmacro vec3-getters (aliases)
+(defmacro vector-getters (aliases)
   `(progn
      ,@(loop for name in aliases
           for index = 0 then (1+ index)
           collect
-            `(defmacro ,name (vec3)
-               `(aref ,vec3 ,,index)))))
+            `(defmacro ,name (some-vector)
+               `(aref ,some-vector ,,index)))))
+
+(defmacro point (x y z)
+  `(vector ,x ,y ,z 1))
 
-(vec3-getters (RED GREEN BLUE))
-(vec3-getters (X Y Z))
+(defmacro vec3 (x y z)
+  `(vector ,x ,y ,z 0))
 
-(defun negate-vec3 (v)
-  (VEC3- (vector 0 0 0) v))
+(vector-getters (RED GREEN BLUE))
+(vector-getters (X Y Z W))
 
-(defun scale-vec3 (v scalar)
-  (VEC3* v (vector scalar scalar scalar)))
+(defun negate-vector (v)
+  (VECTOR- (vec3 0 0 0) v))
 
-(defun magnitude (vec3)
+(defun scale-vector (v scalar)
+  (VECTOR* v (vec3 scalar scalar scalar)))
+
+(defun magnitude (some-vector)
   (labels ((square (n) (expt n 2)))
-    (sqrt (+ (square (X vec3))
-             (square (Y vec3))
-             (square (Z vec3))))))
+    (sqrt (+ (square (X some-vector))
+             (square (Y some-vector))
+             (square (Z some-vector))))))
 
-(defun normalize (vec3)
-  (scale-vec3 vec3 (/ 1 (magnitude vec3))))
+(defun normalize (some-vector)
+  (scale-vector some-vector (/ 1 (magnitude some-vector))))
 
 (defun dot-product (v1 v2)
-  (reduce #'+ (vec3* v1 v2)))
+  (reduce #'+ (VECTOR* v1 v2)))
 
 (defun cross-product (v1 v2)
-  (vector (- (* (Y v1) (Z v2)) (* (Z v1) (Y v2)))
-          (- (* (Z v1) (X v2)) (* (X v1) (Z v2)))
-          (- (* (X v1) (Y v2)) (* (Y v1) (X v2)))))
+  (vec3 (- (* (Y v1) (Z v2)) (* (Z v1) (Y v2)))
+        (- (* (Z v1) (X v2)) (* (X v1) (Z v2)))
+        (- (* (X v1) (Y v2)) (* (Y v1) (X v2)))))
 
 (defun blend (c1 c2)
-  (vec3* c1 c2))
+  (VECTOR* c1 c2))
 
 (defstruct (canvas (:constructor make-canvas (&key width height)))
   width

          
@@ 78,8 85,9 @@ 
         (lines (list))
         (offset 0)
         (previous 0 next)
-        (next (position #\Space text) (when (< (1+ previous) len)
-                                        (position #\Space text :start (1+ previous)))))
+        (next (position #\Space text)
+              (when (< (1+ previous) len)
+                (position #\Space text :start (1+ previous)))))
        ((null next) (progn
                       (push (subseq text offset (1- len)) lines)
                       (nreverse lines)))

          
@@ 189,89 197,89 @@ 
               (aref m j i))))))
 
 ;;; tests:
-(deftest adding-vec3
-  (let ((a (vector 3 -2 5))
-        (b (vector -2 3 1)))
-    (check (equalp (VEC3+ a b)
-                   (vector 1 1 6)))))
+(deftest adding-vectors
+  (let ((a (vec3 3 -2 5))
+        (b (vec3 -2 3 1)))
+    (check (equalp (VECTOR+ a b)
+                   (vec3 1 1 6)))))
 
-(deftest subtracting-vec3
-  (let ((a (vector 3 2 1))
-        (b (vector 5 6 7)))
-    (check (equalp (VEC3- a b)
-                   (vector -2 -4 -6)))))
+(deftest subtracting-vectors
+  (let ((a (vec3 3 2 1))
+        (b (vec3 5 6 7)))
+    (check (equalp (VECTOR- a b)
+                   (vec3 -2 -4 -6)))))
 
 (deftest point-minus-a-vector
-  (let ((p (vector 3 2 1))
-        (v (vector 5 6 7)))
-    (check (equalp (VEC3- p v)
-                   (vector -2 -4 -6)))))
+  (let ((p (point 3 2 1))
+        (v (vec3 5 6 7)))
+    (check (equalp (VECTOR- p v)
+                   (point -2 -4 -6)))))
 
 (deftest subtracting-zero-vector
-  (let ((zero (vector 0 0 0))
-        (v (vector 1 -2 3)))
-    (check (equalp (VEC3- zero v)
-                   (vector -1 2 -3)))))
+  (let ((zero (vec3 0 0 0))
+        (v (vec3 1 -2 3)))
+    (check (equalp (VECTOR- zero v)
+                   (vec3 -1 2 -3)))))
 
 (deftest negating-a-vector
-  (let ((v (vector 1 -2 3)))
-    (check (equalp (negate-vec3 v)
-                   (vector -1 2 -3)))))
+  (let ((v (vec3 1 -2 3)))
+    (check (equalp (negate-vector v)
+                   (vec3 -1 2 -3)))))
 
 (deftest multiply-by-scalar
   (let ((t1 (vector 1 -2 3 -4))
         (scalar 3.5))
-    (check (equalp (scale-vec3 t1 scalar)
-                   (vector 3.5 -7 10.5)))))
+    (check (equalp (scale-vector t1 scalar)
+                   (vec3 3.5 -7 10.5)))))
 
 (deftest multiply-by-scalar-fraction
-  (let ((v (vector 1 -2 3))
+  (let ((v (vec3 1 -2 3))
         (scalar (/ 1 2)))
-    (check (equalp (scale-vec3 v scalar)
-                   (vector 0.5 -1 1.5)))))
+    (check (equalp (scale-vector v scalar)
+                   (vec3 0.5 -1 1.5)))))
 
 (deftest dividing-scalar
-  (let ((v (vector 1 -2 3))
+  (let ((v (vec3 1 -2 3))
         (scalar 2))
-    (check (equalp (scale-vec3 v (/ 1 scalar))
-                   (vector 0.5 -1 1.5)))))
+    (check (equalp (scale-vector v (/ 1 scalar))
+                   (vec3 0.5 -1 1.5)))))
 
 (deftest magnitude-of-vector
-  (let ((v1 (vector 0 1 0))
-        (v2 (vector 0 0 1))
-        (v3 (vector 1 2 3))
-        (v4 (vector -1 -2 -3)))
+  (let ((v1 (vec3 0 1 0))
+        (v2 (vec3 0 0 1))
+        (v3 (vec3 1 2 3))
+        (v4 (vec3 -1 -2 -3)))
     (check (equalp (magnitude v1) 1)
            (equalp (magnitude v2) 1)
            (equalp (magnitude v3) (sqrt 14))
            (equalp (magnitude v4) (sqrt 14)))))
 
 (deftest normalize-vector
-  (let ((v1 (vector 4 0 0))
-        (v2 (vector 1 2 3)))
-    (check (equalp (normalize v1) (vector 1 0 0))
-           (equalp (normalize v2) (vector (/ 1 (sqrt 14))
+  (let ((v1 (vec3 4 0 0))
+        (v2 (vec3 1 2 3)))
+    (check (equalp (normalize v1) (vec3 1 0 0))
+           (equalp (normalize v2) (vec3 (/ 1 (sqrt 14))
                                           (/ 2 (sqrt 14))
                                           (/ 3 (sqrt 14)))))))
 
 (deftest normalize-magnitude-of-vector
-  (let ((v1 (vector 1 2 3)))
+  (let ((v1 (vec3 1 2 3)))
     (check (float= (magnitude (normalize v1)) 1))))
 
 (deftest test-dot-product
-  (let ((v1 (vector 1 2 3))
-        (v2 (vector 2 3 4)))
+  (let ((v1 (vec3 1 2 3))
+        (v2 (vec3 2 3 4)))
     (check (equalp (dot-product v1 v2) 20))))
 
 (deftest cross-product-results
-  (let ((v1 (vector 1 2 3))
-        (v2 (vector 2 3 4)))
-    (check (equalp (cross-product v1 v2) (vector -1 2 -1))
-           (equalp (cross-product v2 v1) (vector 1 -2 1)))))
+  (let ((v1 (vec3 1 2 3))
+        (v2 (vec3 2 3 4)))
+    (check (equalp (cross-product v1 v2) (vec3 -1 2 -1))
+           (equalp (cross-product v2 v1) (vec3 1 -2 1)))))
 
 (deftest vector-basics
-  (adding-vec3)
-  (subtracting-vec3)
+  (adding-vectors)
+  (subtracting-vectors)
   (point-minus-a-vector)
   (subtracting-zero-vector)
   (negating-a-vector)

          
@@ 291,8 299,8 @@ 
            (equalp (BLUE c) 1.7))))
 
 (deftest blending-colors
-  (let* ((c1 (vector 1 0.2 0.4))
-         (c2 (vector 0.9 1 0.1))
+  (let* ((c1 (vec3 1 0.2 0.4))
+         (c2 (vec3 0.9 1 0.1))
          (blended-color (blend c1 c2)))
     (check (float= (RED blended-color) 0.9)
            (float= (GREEN blended-color) 0.2)

          
@@ 529,35 537,183 @@ 0 0 0 0 0 0 0 0 0 0
   (inverse-is-reversible)
   (inverse-undoes-multiplication))
 
+(defun translation (x y z)
+  (make-array '(4 4) :initial-contents `((1 0 0 ,x)
+                                         (0 1 0 ,y)
+                                         (0 0 1 ,z)
+                                         (0 0 0 1))))
+
+(defun scaling (x y z)
+  (make-array '(4 4) :initial-contents `((,x  0  0 0)
+                                         ( 0 ,y  0 0)
+                                         ( 0  0 ,z 0)
+                                         ( 0  0  0  1))))
+
+(defun rotation-x (r)
+  (make-array '(4 4) :initial-contents `((1     0         0      0)
+                                         (0 ,(cos r) ,(-(sin r)) 0)
+                                         (0 ,(sin r)   ,(cos r)  0)
+                                         (0     0         0      1))))
+
+(defun rotation-y (r)
+  (make-array '(4 4) :initial-contents `((,(cos r)    0 ,(sin r) 0)
+                                         (    0       1     0    0)
+                                         (,(-(sin r)) 0 ,(cos r) 0)
+                                         (    0       0     0    1))))
+
+(defun rotation-z (r)
+  (make-array '(4 4) :initial-contents `((,(cos r) ,(- (sin r)) 0 0)
+                                         (,(sin r) ,(cos r)     0 0)
+                                         (    0        0        1 0)
+                                         (    0        0        0 1))))
+
+(deftest multiplying-translation-matrix
+  (let ((transform (translation 5 -3 2))
+        (p (point -3 4 5)))
+    (check (equalp (matrix*vector transform p)
+                   (point 2 1 7)))))
+
+(deftest multiply-inverse-of-translation
+  (let* ((transform (translation 5 -3 2))
+         (inv (inverse transform))
+         (p (point -3 4 5)))
+    (check (equalp (matrix*vector inv p)
+                   (point -8 7 3)))))
+
+(deftest translation-of-a-vector
+  (let ((transform (translation 5 -3 2))
+        (v (vec3 -3 4 5)))
+    (check (equalp (matrix*vector transform v) v))))
+
+(deftest scaling-matrix-to-point
+  (let ((transform (scaling 2 3 4))
+        (p (point -4 6 8)))
+    (check (equalp (matrix*vector transform p)
+                   (point -8 18 32)))))
+
+(deftest scaling-matrix-to-vector
+  (let ((transform (scaling 2 3 4))
+        (v (vec3 -4 6 8)))
+    (check (equalp (matrix*vector transform v)
+                   (vec3 -8 18 32)))))
+
+(deftest multiplying-inverse-of-scaling-matrix
+  (let* ((transform (scaling 2 3 4))
+         (inv (inverse transform))
+         (v (vec3 -4 6 8)))
+    (check (equalp (matrix*vector inv v)
+                   (vec3 -2 2 2)))))
+
+(deftest scaling-negatively-is-reflection
+  (let ((transform (scaling -1 1 1))
+        (p (point 2 3 4)))
+    (check (equalp (matrix*vector transform p)
+                   (point -2 3 4)))))
+
+(deftest rotation-around-x-axis
+  (let* ((p (point 0 1 0))
+         (half-quarter (rotation-x (/ pi 4)))
+         (full-quarter (rotation-x (/ pi 2)))
+         (hq-values (loop for i across (matrix*vector half-quarter p)
+                       collect i))
+         (hq-actuals (loop for i across (point 0 (/ (sqrt 2) 2) (/ (sqrt 2) 2))
+                        collect i))
+         (fq-values (loop for i across (matrix*vector full-quarter p)
+                       collect i))
+         (fq-actuals (loop for i across (point 0 0 1)
+                        collect i)))
+    (check (not (member nil (mapcar #'float= hq-values hq-actuals)))
+           (not (member nil (mapcar #'float= fq-values fq-actuals))))))
+
+(deftest inverse-x-rotation-rotates-oppositely
+  (let* ((p (point 0 1 0))
+         (half-quarter (rotation-x (/ pi 4)))
+         (inv (inverse half-quarter))
+         (inverse-values (loop for i across (matrix*vector inv p)
+                            collect i))
+         (inverse-actuals (loop for i across (point 0 (/ (sqrt 2) 2) (-(/ (sqrt 2) 2)))
+                             collect i)))
+    (check
+      ;; is there a better way to do float= over two arrays?
+      (not (member nil (mapcar #'float= inverse-values inverse-actuals))))))
+
+(deftest rotation-around-y-axis
+  (let* ((p (point 0 0 1))
+         (half-quarter (rotation-y (/ pi 4)))
+         (full-quarter (rotation-y (/ pi 2)))
+         (hq-values (loop for i across (matrix*vector half-quarter p)
+                       collect i))
+         (hq-actuals (loop for i across (point (/ (sqrt 2) 2)
+                                               0
+                                               (/ (sqrt 2) 2))
+                        collect i))
+         (fq-values (loop for i across (matrix*vector full-quarter p)
+                       collect i))
+         (fq-actuals (loop for i across (point 1 0 0)
+                        collect i)))
+    (check (not (member nil (mapcar #'float= hq-values hq-actuals)))
+           (not (member nil (mapcar #'float= fq-values fq-actuals))))))
+
+(deftest rotation-around-z-axis
+  (let* ((p (point 0 1 0))
+         (half-quarter (rotation-z (/ pi 4)))
+         (full-quarter (rotation-z (/ pi 2)))
+         (hq-values (loop for i across
+                         (matrix*vector half-quarter p)
+                       collect i))
+         (hq-actuals (loop for i across
+                          (point (- (/ (sqrt 2) 2)) (/ (sqrt 2) 2) 0)
+                        collect i))
+         (fq-values (loop for i across (matrix*vector full-quarter p)
+                       collect i))
+         (fq-actuals (loop for i across (point -1 0 0)
+                        collect i)))
+    (check (not (member nil (mapcar #'float= hq-values hq-actuals)))
+           (not (member nil (mapcar #'float= fq-values fq-actuals))))))
+
+(deftest matrix-transformations
+  (multiplying-translation-matrix)
+  (multiply-inverse-of-translation)
+  (translation-of-a-vector)
+  (scaling-matrix-to-point)
+  (scaling-matrix-to-vector)
+  (multiplying-inverse-of-scaling-matrix)
+  (scaling-negatively-is-reflection)
+  (rotation-around-x-axis)
+  (inverse-x-rotation-rotates-oppositely)
+  (rotation-around-y-axis)
+  (rotation-around-z-axis))
+
 (deftest suite
   (vector-basics)
   (canvas-and-visuals)
-  (matrix-basics))
+  (matrix-basics)
+  (matrix-transformations))
 
 ;;; sample
 (defstruct projectile position velocity)
 (defstruct environment gravity wind)
 
 (defun tick (env proj)
-  (let ((p (VEC3+ (projectile-position proj) (projectile-velocity proj)))
-        (v (VEC3+ (projectile-velocity proj) (VEC3+ (environment-gravity env)
+  (let ((p (VECTOR+ (projectile-position proj) (projectile-velocity proj)))
+        (v (VECTOR+ (projectile-velocity proj) (VECTOR+ (environment-gravity env)
                                                     (environment-wind env)))))
     (make-projectile :position p :velocity v)))
 
 (defun projectile-tracking ()
-  (do* ((env (make-environment :gravity (vector 0 -0.1 0)
-                               :wind (vector -0.01 0 0)))
-        (p (make-projectile :position (vector 0 1 0)
-                            :velocity (normalize (vector 1 1 0)))
+  (do* ((env (make-environment :gravity (vec3 0 -0.1 0)
+                               :wind (vec3 -0.01 0 0)))
+        (p (make-projectile :position (point 0 1 0)
+                            :velocity (normalize (vec3 1 1 0)))
            (tick env p)))
        ((if (<= (Y (projectile-position p)) 0) (return p)))
     (format t "~a~%" (projectile-position p))))
 
 (defun projectile-visualization ()
-  (do* ((env (make-environment :gravity (vector 0 -0.093 0)
-                               :wind (vector -0.01 0 0)))
-        (p (make-projectile :position (vector 0 1 0)
-                            :velocity (scale-vec3 (normalize (vector 1 1.5 0)) 10))
+  (do* ((env (make-environment :gravity (vec3 0 -0.093 0)
+                               :wind (vec3 -0.01 0 0)))
+        (p (make-projectile :position (point 0 1 0)
+                            :velocity (scale-vector (normalize (vec3 1 1.5 0)) 10))
            (tick env p))
         (c (make-canvas :width 900 :height 550)))
        ((<= (Y (projectile-position p)) 0)