wip geometry

This commit is contained in:
water111
2026-08-08 01:44:11 -04:00
parent 6fd094fc17
commit 0667197a05
3 changed files with 78 additions and 122 deletions
+2 -2
View File
@@ -3224,8 +3224,8 @@ and max-speed. Leave the velocity unchanged otherwise."
(function matrix quaternion matrix))
(define-extern quaternion-left-mult-matrix! "Place quaternion coefficients into a matrix. Unused." (function matrix quaternion matrix))
(define-extern matrix-with-scale->quaternion "Remove independent row scale from src-mat and convert the remaining rotation to dst." (function quaternion matrix quaternion))
(define-extern quaternion-log! "Store the vector part of the logarithm of unit src in dst.xyz. Unused in Jak 1." (function quaternion quaternion quaternion))
(define-extern quaternion-exp! "Exponentiate a pure-vector quaternion src into unit dst. Unused in Jak 1." (function quaternion quaternion quaternion))
(define-extern quaternion-log! "Store the logarithm of unit src in dst.xyz. Unused in Jak 1." (function quaternion quaternion quaternion))
(define-extern quaternion-exp! "Exponentiate a quaternion src into unit dst. Unused in Jak 1." (function quaternion quaternion quaternion))
(define-extern quaternion-slerp! "Spherically interpolate from a to b, selecting the shortest quaternion hemisphere and using normalized lerp near zero angle." (function quaternion quaternion quaternion float quaternion))
(define-extern quaternion-pseudo-slerp! "Normalized linear interpolation with hemisphere correction. It is not constant-speed and degenerates for antipodal rotations; unused in Jak 1." (function quaternion quaternion quaternion float quaternion))
(define-extern quaternion-zxy! "Build dst from Z, then X, then Y rotations stored in angles." (function quaternion vector quaternion))
+17 -61
View File
@@ -5,72 +5,25 @@
(require "engine/math/quaternion.gc")
(require "engine/gfx/font.gc")
;; Geometry functions are common vector/plane utilities + the "curve" stuff
;; Conventions used throughout this file.
;;
;; Output first. Every function that produces a vector, matrix or quaternion takes it as the first
;; argument and returns it, so calls chain. Whether that output may share storage with an input
;; varies and each docstring says so; where it does not say so, assume it may not.
;;
;; The w lane is never incidental. init-vf0-vector loads vf0 with (0, 0, 0, 1), the value VU0's
;; hardwired vf0 holds, so (.mov.vf.w dst vf0) is how a function sets w to 1.0 and
;; (.add.x.vf.w dst vf0 vf0) is the same thing written as 0.0 + 1.0. An operation masked to xyz,
;; including every outer product, leaves the destination register's w alone, which is why that w
;; write usually appears well before the store rather than next to it.
;;
;; Units. Rotations are 65536 per revolution, so (degrees 1) is 182.044 and sin, cos, acos and
;; atan all take and return that scale. Distances are 4096 per meter, so (meters 1) is 4096.0. A
;; distance returned from this file is a true distance unless the name ends in -squared.
;;
;; A plane arrives in one of two forms and the argument type does not distinguish them. The plane
;; type is ax + by + cz = d with d in w, which is what plane-volume-intersect-dist and
;; point-in-vol? consume. vector-plane-distance instead reads its plane argument's xyz as a point
;; on the plane and takes the normal as a separate argument, ignoring w completely.
;;
;; Three VU0 idioms recur here.
;;
;; A cross product is a pair of instructions: (.outer.product.a.vf acc a b) followed by
;; (.outer.product.b.vf dst b a acc) leaves a x b in dst.xyz. Applying it twice computes a
;; projection, because n x (v x n) expands to v * (n . n) - n * (n . v), which is the component of
;; v lying in the plane through the origin with normal n as long as n is unit length. A non-unit
;; normal scales the result by n . n and the function has no way to notice. vector-flatten! and all
;; four vector-reflect* functions are built from that double cross, so every one of them requires a
;; unit normal.
;;
;; (.mov gpr vf) and (.mov vf gpr) move bits between the register files instead of converting
;; them. The move width comes from the destination's declared register class, not from the source:
;; a float destination moves 32 bits and so picks up lane x, while an int or uint destination moves
;; 64 bits and picks up lanes x and y together. closest-pt-in-triangle and point-in-triangle-cross
;; depend on the wide form. Each sums an edge test into lane y, moves the xy pair into a general
;; register, and reads the sign of the high word. Retyping such a local changes which lane the
;; sign test examines, so those declarations are load-bearing.
;;
;; Given that, testing a float's sign is an integer branch: shift the moved bits right by 63 to
;; get one bit per test, or or the carriers together and compare against zero to ask whether any
;; of them was negative. Testing a float for zero is the same trick and is exact, but it is not
;; symmetric about zero, since -0.0 has a nonzero bit pattern and reads as nonzero.
;; circle-circle-xz-intersect tests its two center deltas that way.
(defun vector-flatten! ((dst vector) (src vector) (plane-normal vector))
"Get the projection of src onto a plane with the given normal
The normal should have magnitude 1.0."
;; dst may be src or plane-normal, since both are in registers before the store. dst.w is 1.0.
;; A zero normal produces a zero result and a non-unit one scales the answer by its squared
;; length; neither is detected.
(rlet ((acc :class vf)
(vf0 :class vf)
(source :class vf) ;; src
(normal :class vf) ;; normal
(source :class vf)
(normal :class vf)
(in-plane :class vf))
(init-vf0-vector)
(.lvf source (&-> src quad))
(.lvf normal (&-> plane-normal quad))
(.mov.vf.w in-plane vf0)
;; has the right magnitude, but rotation is off by 90 degrees
;; use the triple-product to compute this without needing
;; a slow divide.
;; vec x normal gives us a vector with the right magnitude,
;; and in plane, but rotated 90 degrees.
(.outer.product.a.vf acc source normal)
(.outer.product.b.vf in-plane normal source acc)
;; rotate by 90 about normal of plane
;; rotate by 90 about normal of plane to get the final vector.
(.outer.product.a.vf acc normal in-plane)
(.outer.product.b.vf in-plane in-plane normal acc)
(.svf (&-> dst quad) in-plane)
@@ -93,13 +46,15 @@
(.lvf source (&-> src quad))
(.lvf normal (&-> plane-normal quad))
(.mov.vf.w in-plane vf0)
;; triple-product to project onto the plane
(.outer.product.a.vf acc source normal)
(.outer.product.b.vf in-plane normal source acc)
;; in-plane is the projection on the plane
(.outer.product.a.vf acc normal in-plane)
(.outer.product.b.vf in-plane in-plane normal acc)
(.add.vf.xyz acc in-plane in-plane) ;; double that part
(.sub.mul.w.vf.xyz in-plane source vf0 acc) ;; and subtract the original
;; add in-plane to itself to get 2 * T
(.add.vf.xyz acc in-plane in-plane)
;; subtract src.
(.sub.mul.w.vf.xyz in-plane source vf0 acc)
(.svf (&-> dst quad) in-plane)
dst))
@@ -117,9 +72,9 @@
(.lvf source (&-> src quad))
(.lvf normal (&-> plane-normal quad))
(.mov.vf.w in-plane vf0)
;; triple-product to project onto plane
(.outer.product.a.vf acc source normal)
(.outer.product.b.vf in-plane normal source acc)
;; part on the plane (requires normal to be unit)
(.outer.product.a.vf acc normal in-plane)
(.outer.product.b.vf in-plane in-plane normal acc)
(.add.vf.xyz in-plane in-plane normal) ;; add normal to that.
@@ -128,8 +83,7 @@
(defun vector-reflect-true-flat! ((dst vector) (src vector) (plane-normal vector))
"Project src onto the plane with unit normal plane-normal; equivalent to vector-flatten!."
;; Instruction for instruction the same function as vector-flatten!, kept as a separate symbol so
;; callers can name which of the reflect variants they mean. dst may be either input; dst.w is 1.0.
;; Instruction for instruction the same function as vector-flatten!
(rlet ((acc :class vf)
(vf0 :class vf)
(source :class vf)
@@ -149,7 +103,9 @@
(defun vector-reflect-flat-above! ((dst vector) (src vector) (plane-normal vector))
"Project src onto the plane with unit normal plane-normal, then lift the result off the plane by
0.32 times its own length, capped at four meters. dst is used as scratch, so it must not alias
plane-normal, though it may alias src. dst.w is 1.0."
plane-normal, though it may alias src."
;; this may just have a bug, normal-component below is always 0.
;; in any case, this is a bit of a softer "bounce" than reflect.
(rlet ((acc :class vf)
(vf0 :class vf)
(source :class vf)
+59 -59
View File
@@ -13,10 +13,7 @@
(defun set-eul! ((dst euler-angles) (angle0 float) (angle1 float) (angle2 float) (order int))
"Store three ordered Euler angles and their packed Shoemake order code."
(set! (-> dst data 0) angle0)
(set! (-> dst data 1) angle1)
(set! (-> dst data 2) angle2)
(set! (-> dst data 3) (the float order))
(set-vector! dst angle0 angle1 angle2 (the float order))
dst)
(defun eul->matrix ((dst-mat matrix) (src euler-angles))
@@ -81,67 +78,70 @@
(defun matrix->eul ((dst euler-angles) (src-mat matrix) (order int))
"Convert src-mat to the convention selected by packed Shoemake order. Repeated-axis and three-distinct-axis orders use separate formulas; near a coordinate singularity the underdetermined final angle is set to zero."
(let* ((parity (logand (sar order 2) 1))
(axis0 (-> EulSafe (logand (sar order 3) 3)))
0
0
0
(let* ((parity (logand (/ order 4) 1))
(axis0 (-> EulSafe (logand (/ order 8) 3)))
(axis1 (-> EulNext (+ axis0 parity)))
(axis2 (-> EulNext (+ (- 1 parity) axis0))))
(if (= (logand (sar order 1) 1) 1)
(let* ((f0-0 (-> (the-as (pointer float) (+ (+ (shl axis1 2) (shl axis0 4)) (the-as int src-mat)))))
(f0-2 (* f0-0 f0-0))
(f1-0 (-> (the-as (pointer float) (+ (+ (shl axis2 2) (shl axis0 4)) (the-as int src-mat)))))
(middle-sine-magnitude (sqrtf (+ f0-2 (* f1-0 f1-0)))))
(cond
((< 0.00001 middle-sine-magnitude)
(set! (-> dst data 0)
(atan (-> (the-as (pointer float) (+ (+ (shl axis1 2) (shl axis0 4)) (the-as int src-mat))))
(-> (the-as (pointer float) (+ (+ (shl axis2 2) (shl axis0 4)) (the-as int src-mat))))))
(set! (-> dst data 1)
(atan middle-sine-magnitude (-> (the-as (pointer float) (+ (+ (shl axis0 2) (shl axis0 4)) (the-as int src-mat))))))
(let ((f0-13 (atan (-> (the-as (pointer float) (+ (+ (shl axis0 2) (shl axis1 4)) (the-as int src-mat))))
(- (-> (the-as (pointer float) (+ (+ (shl axis0 2) (shl axis2 4)) (the-as int src-mat))))))))
(set! (-> dst data 2) f0-13)))
(else
(set! (-> dst data 0)
(atan (- (-> (the-as (pointer float) (+ (+ (shl axis2 2) (shl axis1 4)) (the-as int src-mat)))))
(-> (the-as (pointer float) (+ (+ (shl axis1 2) (shl axis1 4)) (the-as int src-mat))))))
(set! (-> dst data 1)
(atan middle-sine-magnitude (-> (the-as (pointer float) (+ (+ (shl axis0 2) (shl axis0 4)) (the-as int src-mat))))))
(let ((f0-20 0.0)) (set! (-> dst data 2) f0-20)))))
(let* ((f0-21 (-> (the-as (pointer float) (+ (+ (shl axis0 2) (shl axis0 4)) (the-as int src-mat)))))
(f0-23 (* f0-21 f0-21))
(f1-3 (-> (the-as (pointer float) (+ (+ (shl axis0 2) (shl axis1 4)) (the-as int src-mat)))))
(middle-cosine-magnitude (sqrtf (+ f0-23 (* f1-3 f1-3)))))
(cond
((< 0.00001 middle-cosine-magnitude)
(set! (-> dst data 0)
(atan (-> (the-as (pointer float) (+ (+ (shl axis1 2) (shl axis2 4)) (the-as int src-mat))))
(-> (the-as (pointer float) (+ (+ (shl axis2 2) (shl axis2 4)) (the-as int src-mat))))))
(set! (-> dst data 1)
(atan (- (-> (the-as (pointer float) (+ (+ (shl axis0 2) (shl axis2 4)) (the-as int src-mat)))))
middle-cosine-magnitude))
(let ((f0-34 (atan (-> (the-as (pointer float) (+ (+ (shl axis0 2) (shl axis1 4)) (the-as int src-mat))))
(-> (the-as (pointer float) (+ (+ (shl axis0 2) (shl axis0 4)) (the-as int src-mat)))))))
(set! (-> dst data 2) f0-34)))
(else
(set! (-> dst data 0)
(atan (- (-> (the-as (pointer float) (+ (+ (shl axis2 2) (shl axis1 4)) (the-as int src-mat)))))
(-> (the-as (pointer float) (+ (+ (shl axis1 2) (shl axis1 4)) (the-as int src-mat))))))
(set! (-> dst data 1)
(atan (- (-> (the-as (pointer float) (+ (+ (shl axis0 2) (shl axis2 4)) (the-as int src-mat)))))
middle-cosine-magnitude))
(let ((f0-42 0.0)) (set! (-> dst data 2) f0-42)))))))
(when (= (logand (sar order 2) 1) 1)
(set! (-> dst data 0) (- (-> dst data 0)))
(set! (-> dst data 1) (- (-> dst data 1)))
(let ((f0-48 (- (-> dst data 2)))) (set! (-> dst data 2) f0-48)))
(if (= (logand order 1) 1)
(let ((f0-49 (-> dst data 0))) (set! (-> dst data 0) (-> dst data 2)) (set! (-> dst data 2) f0-49)))
(set! (-> dst data 3) (the float order))
(cond
((= (logand (/ order 2) 1) 1)
(let ((middle-sine-magnitude (sqrtf (+ (square (-> (the-as (pointer float) (+ (+ (* axis1 4) (* axis0 16)) (the-as int src-mat)))))
(square (-> (the-as (pointer float) (+ (+ (* axis2 4) (* axis0 16)) (the-as int src-mat)))))))))
(cond
((< 0.00001 middle-sine-magnitude)
(set! (-> dst x)
(atan (-> (the-as (pointer float) (+ (+ (* axis1 4) (* axis0 16)) (the-as int src-mat))))
(-> (the-as (pointer float) (+ (+ (* axis2 4) (* axis0 16)) (the-as int src-mat))))))
(set! (-> dst y)
(atan middle-sine-magnitude (-> (the-as (pointer float) (+ (+ (* axis0 4) (* axis0 16)) (the-as int src-mat))))))
(set! (-> dst z)
(atan (-> (the-as (pointer float) (+ (+ (* axis0 4) (* axis1 16)) (the-as int src-mat))))
(- (-> (the-as (pointer float) (+ (+ (* axis0 4) (* axis2 16)) (the-as int src-mat))))))))
(else
(set! (-> dst x)
(atan (- (-> (the-as (pointer float) (+ (+ (* axis2 4) (* axis1 16)) (the-as int src-mat)))))
(-> (the-as (pointer float) (+ (+ (* axis1 4) (* axis1 16)) (the-as int src-mat))))))
(set! (-> dst y)
(atan middle-sine-magnitude (-> (the-as (pointer float) (+ (+ (* axis0 4) (* axis0 16)) (the-as int src-mat))))))
(set! (-> dst z) 0.0)))))
(else
(let ((middle-cosine-magnitude (sqrtf (+ (square (-> (the-as (pointer float) (+ (+ (* axis0 4) (* axis0 16)) (the-as int src-mat)))))
(square (-> (the-as (pointer float) (+ (+ (* axis0 4) (* axis1 16)) (the-as int src-mat)))))))))
(cond
((< 0.00001 middle-cosine-magnitude)
(set! (-> dst x)
(atan (-> (the-as (pointer float) (+ (+ (* axis1 4) (* axis2 16)) (the-as int src-mat))))
(-> (the-as (pointer float) (+ (+ (* axis2 4) (* axis2 16)) (the-as int src-mat))))))
(set! (-> dst y)
(atan (- (-> (the-as (pointer float) (+ (+ (* axis0 4) (* axis2 16)) (the-as int src-mat))))) middle-cosine-magnitude))
(set! (-> dst z)
(atan (-> (the-as (pointer float) (+ (+ (* axis0 4) (* axis1 16)) (the-as int src-mat))))
(-> (the-as (pointer float) (+ (+ (* axis0 4) (* axis0 16)) (the-as int src-mat)))))))
(else
(set! (-> dst x)
(atan (- (-> (the-as (pointer float) (+ (+ (* axis2 4) (* axis1 16)) (the-as int src-mat)))))
(-> (the-as (pointer float) (+ (+ (* axis1 4) (* axis1 16)) (the-as int src-mat))))))
(set! (-> dst y)
(atan (- (-> (the-as (pointer float) (+ (+ (* axis0 4) (* axis2 16)) (the-as int src-mat))))) middle-cosine-magnitude))
(set! (-> dst z) 0.0)))))))
(when (= (logand (/ order 4) 1) 1)
(set! (-> dst x) (- (-> dst x)))
(set! (-> dst y) (- (-> dst y)))
(set! (-> dst z) (- (-> dst z))))
(when (= (logand order 1) 1)
(let ((f0-49 (-> dst x)))
(set! (-> dst x) (-> dst z))
(set! (-> dst z) f0-49)))
(set! (-> dst w) (the float order))
dst)
(defun eul->quat ((dst quaternion) (src euler-angles))
"Convert src from its Euler convention to a quaternion."
(let ((rotation-mat (new 'stack-no-clear 'matrix))) (eul->matrix rotation-mat src) (matrix->quaternion dst rotation-mat))
(let ((rotation-mat (new 'stack-no-clear 'matrix)))
(eul->matrix rotation-mat src)
(matrix->quaternion dst rotation-mat))
dst)
(defun quat->eul ((dst euler-angles) (src quaternion) (order int))