@@ 49,8 49,7 @@
(square (Z vec3))))))
(defun normalize (vec3)
- (let ((m (magnitude vec3)))
- (scale-vec3 vec3 (/ 1 m))))
+ (scale-vec3 vec3 (/ 1 (magnitude vec3))))
(defun dot-product (v1 v2)
(reduce #'+ (vec3* v1 v2)))
@@ 63,9 62,7 @@
(defun blend (c1 c2)
(vec3* c1 c2))
-(defstruct (canvas
- (:constructor
- make-canvas (&key width height)))
+(defstruct (canvas (:constructor make-canvas (&key width height)))
width
height
(body (make-array (list width height) :initial-element (vector 0 0 0))))
@@ 137,9 134,12 @@
(aref b k j))))))))
(defun determinant (m)
- ;; only operates on 2x2 matrices
- (- (* (aref m 0 0) (aref m 1 1))
- (* (aref m 0 1) (aref m 1 0))))
+ (if (equalp (array-dimensions m) '(2 2))
+ (- (* (aref m 0 0) (aref m 1 1))
+ (* (aref m 0 1) (aref m 1 0)))
+ (let ((result 0))
+ (dotimes (i (array-dimension m 0) result)
+ (incf result (* (cofactor m 0 i) (aref m 0 i)))))))
(defun submatrix (m row column)
(let ((result (make-array (mapcar #'1- (array-dimensions m)))))
@@ 152,8 152,28 @@
(if (> i row) (setf ii (1- i)))
(if (> j column) (setf jj (1- j)))
(setf (aref result ii jj)
- (aref m i j)))))
- result))
+ (aref m i j))))
+ finally (return result))))
+
+(defun minor (m row column)
+ (determinant (submatrix m row column)))
+
+(defun cofactor (m row column)
+ (if (oddp (+ row column))
+ (- (minor m row column))
+ (minor m row column)))
+
+(defun invertible? (m)
+ (not (eq (determinant m) 0)))
+
+(defun inverse (m)
+ (when (invertible? m)
+ (let ((d (determinant m))
+ (result (make-array (array-dimensions m))))
+ (dotimes (row (array-dimension m 0) result)
+ (dotimes (column (array-dimension m 1))
+ (setf (aref result column row)
+ (/ (cofactor m row column) d)))))))
;;; since we're never taking the identity but for a 4x4 matrix...
(defvar identity-matrix #2A((1 0 0 0)
@@ 161,97 181,12 @@
(0 0 1 0)
(0 0 0 1)))
-(defun transpose-matrix (a)
- (let ((result (make-array (array-dimensions a))))
- (dotimes (i (array-dimension a 0) result)
- (dotimes (j (array-dimension a 1))
+(defun transpose-matrix (m)
+ (let ((result (make-array (array-dimensions m))))
+ (dotimes (i (array-dimension m 0) result)
+ (dotimes (j (array-dimension m 1))
(setf (aref result i j)
- (aref a j i))))))
-
-(deftest multiply-4x4-matrices
- (let ((a #2A((1 2 3 4)
- (5 6 7 8)
- (9 8 7 6)
- (5 4 3 2)))
- (b #2A((-2 1 2 3)
- ( 3 2 1 -1)
- ( 4 3 6 5)
- ( 1 2 7 8))))
- (check (equalp (matrix*matrix a b)
- #2A((20 22 50 48)
- (44 54 114 108)
- (40 58 110 102)
- (16 26 46 42))))))
-
-(deftest multiply-matrix-with-vector
- (let ((a #2A((1 2 3 4)
- (2 4 4 2)
- (8 6 4 1)
- (0 0 0 1)))
- (b (vector 1 2 3 1)))
- (check (equalp (matrix*vector a b) (vector 18 24 33 1)))))
-
-(deftest identity-matrix-returns-original-matrix
- (check (equalp
- (matrix*matrix identity-matrix #2A((1 2 3 4)
- (2 4 4 2)
- (8 6 4 1)
- (0 0 0 1)))
- #2A((1 2 3 4)
- (2 4 4 2)
- (8 6 4 1)
- (0 0 0 1)))))
-
-(deftest identity-matrix-times-tuple
- (check (equalp
- (matrix*vector identity-matrix (vector 1 2 3 4))
- (vector 1 2 3 4))))
-
-(deftest transpose-a-matrix
- (let ((a #2A((0 9 3 0)
- (9 8 0 8)
- (1 8 5 3)
- (0 0 5 6))))
- (check (equalp (transpose-matrix a) #2A((0 9 1 0)
- (9 8 8 0)
- (3 0 5 5)
- (0 8 3 6))))))
-
-(deftest transpose-identity-matrix
- (check (equalp (transpose-matrix identity-matrix)
- identity-matrix)))
-
-(deftest verify-2x2-determinant
- (let ((a #2A(( 1 5)
- (-3 2))))
- (check (equalp (determinant a) 17))))
-
-(deftest submatrix-of-3x3-is-2x2
- (let ((m #2A(( 1 5 0)
- (-3 2 7)
- ( 0 6 -3))))
- (check (equalp (submatrix m 0 2) #2A((-3 2)
- ( 0 6))))))
-
-(deftest submatrix-of-4x4-is-3x3
- (let ((m #2A((-6 1 1 6)
- (-8 5 8 6)
- (-1 0 8 2)
- (-7 1 -1 1))))
- (check (equalp (submatrix m 2 1) #2A((-6 1 6)
- (-8 8 6)
- (-7 -1 1))))))
-
-(deftest matrix-basics
- (multiply-4x4-matrices)
- (multiply-matrix-with-vector)
- (identity-matrix-returns-original-matrix)
- (identity-matrix-times-tuple)
- (transpose-a-matrix)
- (transpose-identity-matrix)
- (verify-2x2-determinant)
- (submatrix-of-3x3-is-2x2)
- (submatrix-of-4x4-is-3x3))
+ (aref m j i))))))
;;; tests:
(deftest adding-vec3
@@ 377,7 312,7 @@
(let ((canvas (make-canvas :width 10 :height 20))
(color (vector 1 0 0)))
(write-pixel canvas 2 3 color)
- (check (equalp (aref (canvas-body canvas) 3 2) (vector 1 0 0)))))
+ (check (equalp (aref (canvas-body canvas) 2 3) (vector 1 0 0)))))
(deftest check-canvas-pixel-string
(let ((c (make-canvas :width 1 :height 2)))
@@ 410,9 345,172 @@ 0 0 0 0 0 0 0 0 0 0
(output-canvas-to-ppm)
(ppm-includes-trailing-newline))
+(deftest multiply-4x4-matrices
+ (let ((a #2A((1 2 3 4)
+ (5 6 7 8)
+ (9 8 7 6)
+ (5 4 3 2)))
+ (b #2A((-2 1 2 3)
+ ( 3 2 1 -1)
+ ( 4 3 6 5)
+ ( 1 2 7 8))))
+ (check (equalp (matrix*matrix a b)
+ #2A((20 22 50 48)
+ (44 54 114 108)
+ (40 58 110 102)
+ (16 26 46 42))))))
+
+(deftest multiply-matrix-with-vector
+ (let ((a #2A((1 2 3 4)
+ (2 4 4 2)
+ (8 6 4 1)
+ (0 0 0 1)))
+ (b (vector 1 2 3 1)))
+ (check (equalp (matrix*vector a b) (vector 18 24 33 1)))))
+
+(deftest identity-matrix-returns-original-matrix
+ (check (equalp
+ (matrix*matrix identity-matrix #2A((1 2 3 4)
+ (2 4 4 2)
+ (8 6 4 1)
+ (0 0 0 1)))
+ #2A((1 2 3 4)
+ (2 4 4 2)
+ (8 6 4 1)
+ (0 0 0 1)))))
+
+(deftest identity-matrix-times-tuple
+ (check (equalp
+ (matrix*vector identity-matrix (vector 1 2 3 4))
+ (vector 1 2 3 4))))
+
+(deftest transpose-a-matrix
+ (let ((a #2A((0 9 3 0)
+ (9 8 0 8)
+ (1 8 5 3)
+ (0 0 5 6))))
+ (check (equalp (transpose-matrix a) #2A((0 9 1 0)
+ (9 8 8 0)
+ (3 0 5 5)
+ (0 8 3 6))))))
+
+(deftest transpose-identity-matrix
+ (check (equalp (transpose-matrix identity-matrix)
+ identity-matrix)))
+
+(deftest verify-2x2-determinant
+ (let ((a #2A(( 1 5)
+ (-3 2))))
+ (check (equalp (determinant a) 17))))
+
+(deftest submatrix-of-3x3-is-2x2
+ (let ((m #2A(( 1 5 0)
+ (-3 2 7)
+ ( 0 6 -3))))
+ (check (equalp (submatrix m 0 2) #2A((-3 2)
+ ( 0 6))))))
+
+(deftest submatrix-of-4x4-is-3x3
+ (let ((m #2A((-6 1 1 6)
+ (-8 5 8 6)
+ (-1 0 8 2)
+ (-7 1 -1 1))))
+ (check (equalp (submatrix m 2 1) #2A((-6 1 6)
+ (-8 8 6)
+ (-7 -1 1))))))
+
+(deftest 3x3-matrix-minor
+ (let* ((a #2A((3 5 0)
+ (2 -1 7)
+ (6 -1 5)))
+ (b (submatrix a 1 0)))
+ (check (equalp (determinant b) 25)
+ (equalp (minor a 1 0) 25))))
+
+(deftest 3x3-matrix-cofactor
+ (let ((a #2A((3 5 0)
+ (2 -1 -7)
+ (6 -1 5))))
+ (check (equalp (minor a 0 0) -12)
+ (equalp (cofactor a 0 0) -12)
+ (equalp (minor a 1 0) 25)
+ (equalp (cofactor a 1 0) -25))))
+
+(deftest 3x3-matrix-determinant
+ (let ((a #2A(( 1 2 6)
+ (-5 8 -4)
+ ( 2 6 4))))
+ (check (equalp (cofactor a 0 0) 56)
+ (equalp (cofactor a 0 1) 12)
+ (equalp (cofactor a 0 2) -46)
+ (equalp (determinant a) -196))))
+
+(deftest 4x4-matrix-determinant
+ (let ((a #2A((-2 -8 3 5)
+ (-3 1 7 3)
+ ( 1 2 -9 6)
+ (-6 7 7 -9))))
+ (check (equalp (cofactor a 0 0) 690)
+ (equalp (cofactor a 0 1) 447)
+ (equalp (cofactor a 0 2) 210)
+ (equalp (cofactor a 0 3) 51)
+ (equalp (determinant a) -4071))))
+
+(deftest invertibility-predicate
+ (let ((a #2A((6 4 4 4)
+ (5 5 7 6)
+ (4 -9 3 -7)
+ (9 1 7 -6))))
+ (check (equalp (determinant a) -2120)
+ (equalp (invertible? a) t))))
+
+(deftest invertibility-negative-case
+ (let ((a #2A((-4 2 -2 -3)
+ ( 9 6 2 6)
+ ( 0 -5 1 -5)
+ ( 0 0 0 0))))
+ (check (equalp (determinant a) 0)
+ (equalp (invertible? a) nil))))
+
+(deftest calculate-matrix-inverse
+ (let* ((a #2A((-5 2 6 -8)
+ ( 1 -5 1 8)
+ ( 7 7 -6 -7)
+ ( 1 -3 7 4)))
+ (b (inverse a)))
+ (check (equalp (determinant a) 532)
+ (equalp (cofactor a 2 3) -160)
+ (equalp (aref b 3 2) -160/532)
+ (equalp (cofactor a 3 2) 105)
+ (equalp (aref b 2 3) 105/532)
+
+ (equalp b #2A(( 29/133 60/133 32/133 -6/133)
+ (-215/266 -775/532 -59/133 277/532)
+ ( -3/38 -17/76 -1/19 15/76)
+ (-139/266 -433/532 -40/133 163/532))))))
+
+(deftest matrix-basics
+ (multiply-4x4-matrices)
+ (multiply-matrix-with-vector)
+ (identity-matrix-returns-original-matrix)
+ (identity-matrix-times-tuple)
+ (transpose-a-matrix)
+ (transpose-identity-matrix)
+ (verify-2x2-determinant)
+ (submatrix-of-3x3-is-2x2)
+ (submatrix-of-4x4-is-3x3)
+ (3x3-matrix-minor)
+ (3x3-matrix-cofactor)
+ (3x3-matrix-determinant)
+ (4x4-matrix-determinant)
+ (invertibility-predicate)
+ (invertibility-negative-case)
+ (calculate-matrix-inverse))
+
(deftest suite
(vector-basics)
- (canvas-and-visuals))
+ (canvas-and-visuals)
+ (matrix-basics))
;;; sample
(defstruct projectile position velocity)