clean up geometry

This commit is contained in:
water111
2026-08-08 00:07:40 -07:00
parent 0667197a05
commit 6a9f5e540f
5 changed files with 699 additions and 864 deletions
+2 -2
View File
@@ -3315,11 +3315,11 @@ and max-speed. Leave the velocity unchanged otherwise."
(define-extern circle-test "Print a diagnostic example of xz circle intersection." (function none))
(define-extern vector-circle-tangent-new "Write the two external tangent contact points for a pair of xz circles." (function vector vector vector vector none))
(define-extern vector-circle-tangent "Write the two xz tangent points from an external point to a circle." (function vector vector vector vector none))
(define-extern find-knot-span "Find the knot span containing a clamped parameter value." (function int int float (inline-array vector) int))
(define-extern find-knot-span "Find the knot span containing a clamped parameter value." (function int int float (pointer float) int))
(define-extern calculate-basis-functions-vector! "Compute the four nonzero cubic B-spline basis weights for a knot span." (function vector int float (pointer float) vector))
(define-extern curve-closest-point "Refine a normalized cubic B-spline parameter toward the point nearest a target." (function curve vector float float int float float))
(define-extern vector-plane-distance "Return the signed distance from a point to a plane and write the plane normal." (function vector plane vector float))
(define-extern vector-plane-distance "Return the signed distance from a point to a plane and write the plane normal." (function vector vector vector float))
(define-extern curve-get-pos! "Evaluate a curve at normalized input and write its position." (function vector float curve vector))
(define-extern curve-evaluate! "Evaluate a clamped cubic nonuniform B-spline at normalized input." (function vector float (inline-array vector) int (pointer float) int vector))
+21 -3
View File
@@ -1421,10 +1421,21 @@
"args": ["dst", "from", "to", "max-angle"]
},
"matrix-from-two-vectors-max-angle-partial!": {
"args": ["dst", "from", "to", "max-angle", "fraction"]
"args": ["dst", "from", "to", "max-angle", "fraction"],
"vars": {
"s4-1":"axis",
"f28-0":"cos-angle",
"f30-0":"max-cos",
"f0-2":"blended-cos"
}
},
"matrix-from-two-vectors-partial-linear!": {
"args": ["dst", "from", "to", "fraction"]
"args": ["dst", "from", "to", "fraction"],
"vars": {
"gp-1":"axis",
"f0-1":"cos-angle",
"f0-4":"blended-cos"
}
},
"matrix-remove-z-rot": {
"args": ["rotation", "reference"]
@@ -1442,7 +1453,14 @@
"args": ["dst", "from", "to", "t"]
},
"vector-vector-deg-slerp!": {
"args": ["dst", "from", "to", "t", "up"]
"args": ["dst", "from", "to", "t", "up"],
"vars": {
"s0-0":"from-unit",
"s1-0":"to-unit",
"s0-1":"from-rot",
"a2-5":"to-rot",
"a1-6":"blended-rot"
}
},
"normal-of-plane": {
"args": ["dst", "point-a", "point-b", "point-c"]
@@ -0,0 +1,514 @@
;; PS2-specific assembly functions from geometry.gc
(defun circle-circle-xz-intersect ((circle-a sphere) (circle-b sphere) (intersection-a vector) (intersection-b vector))
"Intersect two circles in the xz plane.
Returns -1 for coincident centers, 0 for no real intersection, 1 for a tangent,
or 2 for two intersections, writing the available points to intersection-a and intersection-b.
Each sphere argument supplies x, z and w only; y is ignored. The outputs receive x and z, and
their y and w are left as they were. On a tangent both outputs get the same point, and on 0 or
-1 neither is written at all, so the return value must be checked before the points are used."
;; Solve the radical-axis equation along the larger center delta so the division uses the
;; better-conditioned coordinate.
;;
;; Subtracting the two circle equations cancels the quadratic terms and leaves the radical axis,
;; a straight line: dx*x + dz*z = K/2, where dx and dz are the center deltas and K is
;; ar^2 + bx^2 + bz^2 - br^2 - ax^2 - az^2, built by the opening accumulator chain. Solving that
;; for the coordinate with the larger delta and substituting into circle A gives an ordinary
;; quadratic in the other coordinate. For the branch that divides by dx, the coefficients are
;; A = 1 + (dz/dx)^2
;; B = 2*(ax*dz/dx - (K/2)*dz/dx^2 - az)
;; C = (K/2)^2/dx^2 - K*ax/dx + ax^2 + az^2 - ar^2
;; and the sign of B^2 - 4AC selects between two roots, one root and none. The other branch is
;; the same expression with x and z exchanged. The two roots come out as the pair
;; (sqrt(D) - B)/2A and -(sqrt(D) + B)/2A, and the remaining coordinate is recovered from the
;; radical axis rather than from a second square root.
;;
;; Coincidence is decided on the raw bit patterns of the two deltas, so it is exact but not
;; symmetric: a delta of -0.0 has a nonzero bit pattern and takes the ordinary path, where the
;; division by it produces an infinity.
;;
;; work0 through work12 are f0 through f12. They alias the named registers above rather than
;; being separate storage, so a-x and work3 are one register, half and delta-x-squared and
;; work7 are another, and equation-right and dominance and work6 are a third. Each named value
;; is dead by the time the workN spelling appears for it.
(rlet ((circle-a-reg :reg a0)
(circle-b-reg :reg a1)
(intersection-a-reg :reg a2)
(intersection-b-reg :reg a3)
(result :reg v0)
(bits :reg v1)
(a-x :reg f3)
(a-z :reg f4)
(a-radius :reg f5)
(b-x :reg f1)
(b-z :reg f0)
(b-radius :reg f2)
(equation-right :reg f6)
(half :reg f7)
(delta-x-squared :reg f7)
(delta-z-squared :reg f8)
(dominance :reg f6)
(work0 :reg f0)
(work1 :reg f1)
(work2 :reg f2)
(work3 :reg f3)
(work4 :reg f4)
(work5 :reg f5)
(work6 :reg f6)
(work7 :reg f7)
(work8 :reg f8)
(work9 :reg f9)
(work10 :reg f10)
(work11 :reg f11)
(work12 :reg f12))
(l.s a-x circle-a-reg)
(l.s a-z circle-a-reg 8)
(l.s a-radius circle-a-reg 12)
(l.s b-x circle-b-reg)
(l.s b-z circle-b-reg 8)
(l.s b-radius circle-b-reg 12)
(m! bits (the-as uint 0.5))
(mula.s a-radius a-radius)
(madda.s b-x b-x)
(madda.s b-z b-z)
(msuba.s b-radius b-radius)
(msuba.s a-x a-x)
(msub.s equation-right a-z a-z)
(m half bits)
(sub.s work1 b-x a-x)
(sub.s work2 b-z a-z)
(mul.s work0 equation-right half)
(m circle-a-reg work1)
(abs.s work6 work1)
(m bits work2)
(abs.s work8 work2)
(nop!)
(mul.s delta-x-squared work1 work1)
(nop!)
(sub.s dominance work6 work8)
(b.nz circle-a-reg circle-centers-differ :delay (mul.s delta-z-squared work2 work2))
(b.z bits circle-centers-coincident :delay (nop!))
(label circle-centers-differ)
(m bits dominance)
(nop!)
(b.lt bits 0 solve-with-z :delay (nop!))
;; Substitute z = (equation-right / 2 - delta-z*x) / delta-x.
(m! bits (the-as uint 1.0))
(m work7 bits)
(add.s work6 work7 work7)
(div.s work1 work7 work1)
(mul.s work12 work6 work0)
(mul.s work10 work0 work2)
(mul.s work11 a-x work2)
(mul.s work9 work0 work0)
(mul.s work12 work12 a-x)
(mula.s a-x a-x)
(madda.s a-z a-z)
(msub.s work3 a-radius a-radius)
(mul.s work5 work12 work1)
(mul.s work12 work1 work1)
(mul.s work11 work11 work1)
(mul.s work10 work10 work12)
(mul.s work8 work8 work12)
(mul.s work9 work9 work12)
(sub.s work5 work3 work5)
(sub.s work10 work11 work10)
(add.s work3 work8 work7)
(sub.s work4 work10 a-z)
(add.s work5 work9 work5)
(mul.s work4 work4 work6)
(mul.s work5 work3 work5)
(add.s work6 work6 work6)
(mula.s work4 work4)
(msub.s work5 work6 work5)
(m bits work5)
(sqrt.s work5 work5)
(b.lt bits 0 circles-disjoint :delay (nop!))
(b.z bits circle-x-tangent :delay (nop!))
(m! bits (the-as uint 0.5))
(m work7 bits)
(m! work6 0.0)
(div.s work7 work7 work3)
(add.s work3 work5 work4)
(sub.s work4 work5 work4)
(sub.s work5 work6 work3)
(mul.s work3 work4 work7)
(mul.s work4 work5 work7)
(mul.s work5 work3 work2)
(mul.s work2 work4 work2)
(sub.s work5 work0 work5)
(sub.s work0 work0 work2)
(mul.s work2 work1 work5)
(mul.s work0 work1 work0)
(s.s work2 intersection-a-reg)
(s.s work3 intersection-a-reg 8)
(s.s work0 intersection-b-reg)
(b two-intersections :delay (s.s work4 intersection-b-reg 8))
(label circle-x-tangent)
(m! bits (the-as uint -2.0))
(m work5 bits)
(mul.s work3 work5 work3)
(div.s work3 work4 work3)
(mul.s work2 work3 work2)
(sub.s work0 work0 work2)
(mul.s work0 work0 work1)
(s.s work0 intersection-a-reg)
(s.s work3 intersection-a-reg 8)
(s.s work0 intersection-b-reg)
(b tangent-intersection :delay (s.s work3 intersection-b-reg 8))
(label solve-with-z)
;; The same quadratic with x and z exchanged.
(m! bits (the-as uint 1.0))
(m work8 bits)
(add.s work6 work8 work8)
(div.s work2 work8 work2)
(mul.s work12 work6 work0)
(mul.s work10 work0 work1)
(mul.s work11 a-z work1)
(mul.s work9 work0 work0)
(mul.s work12 work12 a-z)
(mula.s a-z a-z)
(madda.s a-x a-x)
(msub.s work4 a-radius a-radius)
(mul.s work5 work12 work2)
(mul.s work12 work2 work2)
(mul.s work11 work11 work2)
(mul.s work10 work10 work12)
(mul.s work7 delta-x-squared work12)
(mul.s work9 work9 work12)
(sub.s work5 work4 work5)
(sub.s work10 work11 work10)
(add.s work4 work7 work8)
(sub.s work3 work10 a-x)
(add.s work5 work9 work5)
(mul.s work3 work3 work6)
(mul.s work5 work4 work5)
(add.s work6 work6 work6)
(mula.s work3 work3)
(msub.s work5 work6 work5)
(m bits work5)
(sqrt.s work5 work5)
(b.lt bits 0 circles-disjoint :delay (nop!))
(b.z bits circle-z-tangent :delay (nop!))
(m! bits (the-as uint 0.5))
(m work7 bits)
(m! work6 0.0)
(div.s work4 work7 work4)
(add.s work7 work5 work3)
(sub.s work3 work5 work3)
(sub.s work5 work6 work7)
(mul.s work3 work3 work4)
(mul.s work4 work5 work4)
(mul.s work5 work3 work1)
(mul.s work1 work4 work1)
(sub.s work5 work0 work5)
(sub.s work0 work0 work1)
(mul.s work1 work2 work5)
(mul.s work0 work2 work0)
(s.s work3 intersection-a-reg)
(s.s work1 intersection-a-reg 8)
(s.s work4 intersection-b-reg)
(b two-intersections :delay (s.s work0 intersection-b-reg 8))
(label circle-z-tangent)
(m! bits (the-as uint -2.0))
(m work5 bits)
(mul.s work4 work5 work4)
(div.s work3 work3 work4)
(mul.s work1 work3 work1)
(sub.s work0 work0 work1)
(mul.s work0 work0 work2)
(s.s work3 intersection-a-reg)
(s.s work0 intersection-a-reg 8)
(s.s work3 intersection-b-reg)
(b tangent-intersection :delay (s.s work0 intersection-b-reg 8))
(label circles-disjoint)
(b circle-intersection-done :delay (m! result 0))
(label tangent-intersection)
(b circle-intersection-done :delay (m! result 1))
(label two-intersections)
(b circle-intersection-done :delay (m! result 2))
(label circle-centers-coincident)
(b circle-intersection-done :delay (m! result -1))
(label circle-intersection-done)
(j ra :delay (m sp sp))
(nop!)
(nop!)))
(defun calculate-basis-functions-vector! ((dst vector) (span int) (value float) (knots (pointer float)))
"Compute the four nonzero cubic B-spline basis weights for value in span, writing them to dst's
x, y, z and w in that order. They apply to control points span - 3 through span and sum to one
when knots[span] <= value < knots[span + 1]. knots[span - 2] through knots[span + 3] are all
read, so span must be at least 3 and at most knot-count - 5. Any two of those six knots being
equal makes a denominator zero and divides by it."
;; This is the cubic Cox-de Boor recurrence unrolled in the same scalar instruction order.
;; left-j is value - knot[span+1-j], and right-j is knot[span+j] - value.
(rlet ((dst-reg :reg a0)
(span-reg :reg a1)
(value-bits :reg a2)
(knots-reg :reg a3)
(knot-offset)
(span-knot)
(one)
(u)
(left-1)
(left-2)
(left-3)
(right-1)
(right-2)
(right-3)
(denominator)
(degree-1-scale)
(degree-1-temp)
(degree-2-scale-0)
(degree-2-temp-0)
(degree-2-scale-1)
(degree-2-temp-1)
(degree-3-scale-0)
(degree-3-temp-0)
(degree-3-scale-1)
(degree-3-temp-1)
(degree-3-scale-2)
(degree-3-temp-2)
(saved)
(basis-0)
(basis-1)
(basis-2)
(basis-3))
(.sll knot-offset span-reg 2)
(m! one 1.0)
(m u value-bits)
(.addu span-knot knots-reg knot-offset)
;; Degree zero begins with N[0] = 1. The two overwritten zero writes are retained
;; before the reciprocal multiply can use it.
(m basis-0 one)
(m! left-1 0.0)
(m! left-1 0.0)
;; Degree one.
(l.s left-1 span-knot)
(l.s right-1 span-knot 4)
(sub.s left-1 u left-1)
(sub.s right-1 right-1 u)
(add.s denominator right-1 left-1)
(div.s degree-1-scale one denominator)
;; Load degree-two and degree-three distances early to match the FPU schedule.
(l.s left-2 span-knot -4)
(l.s right-2 span-knot 8)
(sub.s left-2 u left-2)
(add.s denominator right-1 left-2)
(l.s left-3 span-knot -8)
(l.s right-3 span-knot 12)
(div.s degree-2-scale-0 one denominator)
;; Finish degree one.
(mul.s degree-1-temp basis-0 degree-1-scale)
(mul.s basis-0 right-1 degree-1-temp)
(mul.s basis-1 left-1 degree-1-temp)
;; Degree two.
(sub.s right-2 right-2 u)
(mul.s degree-2-temp-0 basis-0 degree-2-scale-0)
(add.s denominator right-2 left-1)
(div.s degree-2-scale-1 one denominator)
(mul.s basis-0 right-1 degree-2-temp-0)
(mul.s saved left-2 degree-2-temp-0)
(mul.s degree-2-temp-1 basis-1 degree-2-scale-1)
(mul.s basis-1 right-2 degree-2-temp-1)
(add.s basis-1 basis-1 saved)
(mul.s basis-2 left-1 degree-2-temp-1)
;; Degree three.
(sub.s left-3 u left-3)
(sub.s right-3 right-3 u)
(add.s denominator right-1 left-3)
(div.s degree-3-scale-0 one denominator)
(mul.s degree-3-temp-0 basis-0 degree-3-scale-0)
(mul.s basis-0 right-1 degree-3-temp-0)
(mul.s saved left-3 degree-3-temp-0)
(add.s denominator right-2 left-2)
(div.s degree-3-scale-1 one denominator)
(mul.s degree-3-temp-1 basis-1 degree-3-scale-1)
(mul.s basis-1 right-2 degree-3-temp-1)
(add.s basis-1 basis-1 saved)
(mul.s saved left-2 degree-3-temp-1)
(add.s denominator right-3 left-1)
(div.s degree-3-scale-2 one denominator)
(mul.s degree-3-temp-2 basis-2 degree-3-scale-2)
(mul.s basis-2 right-3 degree-3-temp-2)
(add.s basis-2 basis-2 saved)
(mul.s basis-3 left-1 degree-3-temp-2)
(s.s basis-0 dst-reg)
(s.s basis-1 dst-reg 4)
(s.s basis-2 dst-reg 8)
(s.s basis-3 dst-reg 12)
(m v0 dst-reg)
(j ra :delay (m sp sp))
(nop!)
(nop!)
(nop!)))
(defun curve-evaluate! ((dst vector)
(input float)
(control-points (inline-array vector))
(control-point-count int)
(knots (pointer float))
(knot-count int))
"Evaluate a clamped cubic nonuniform B-spline at normalized input.
Requires at least four controls, knot-count = control-point-count + 4, nondecreasing knots,
and a zero first knot for normalized input. The selected span supplies four controls and their
Cox-de Boor weights; input is parameter progress, not guaranteed distance progress.
input is clamped to the knot domain, so values outside zero to one give the endpoints rather
than extrapolating. dst.w is 1.0, since only xyz accumulate the weighted controls, and dst may
not overlap the control points. control-point-count is accepted but never used."
;; The span lookup first tries the integer part
;; of u as a likely uniform-knot span, then falls back to binary search.
(rlet ((dst-reg :reg a0)
(input-bits :reg a1)
(controls-reg :reg a2)
(control-count-reg :reg a3)
(knots-reg :reg t0)
(knot-count-reg :reg t1)
(basis)
(span :reg s3)
(saved-controls :reg s5)
(saved-dst :reg gp)
(first-knot)
(last-knot)
(scaled-input)
(u)
(u-bits)
(last-index)
(last-offset)
(last-span)
(end-knot-offset)
(end-knot-address)
(end-knot)
(rounded-u-float)
(rounded-u)
(candidate-span)
(candidate-offset)
(candidate-address)
(candidate-low)
(candidate-high)
(low-span)
(high-span)
(mid-span)
(mid-offset)
(mid-address)
(mid-low)
(mid-high)
(basis-fn)
(first-control)
(control-offset)
(control-address)
(result :reg v0)
(curve-point-bits :reg v1)
(weights :class vf)
(control-0 :class vf)
(control-1 :class vf)
(control-2 :class vf)
(control-3 :class vf)
(curve-point :class vf)
(vf0 :class vf)
(acc :class vf))
(init-vf0-vector)
(m saved-dst dst-reg)
(m saved-controls controls-reg)
(m! basis (new 'static 'vector))
;; The curve representation expects first-knot = 0, making this a normalized-domain map.
(l.s first-knot knots-reg)
(set! last-index (- knot-count-reg 1))
(.dsll last-offset last-index 2)
(.daddu end-knot-address knots-reg last-offset)
(l.s last-knot end-knot-address)
(m scaled-input input-bits)
(mul.s scaled-input scaled-input last-knot)
(m u-bits scaled-input)
(m scaled-input u-bits)
(min.s scaled-input scaled-input last-knot)
(max.s u scaled-input first-knot)
(m u-bits u)
(set! last-span (- knot-count-reg 5))
(set! low-span 3)
;; The last endpoint belongs to the final span even though its upper bound is closed.
(set! end-knot-offset (* (+ last-span 1) 4))
(.daddu end-knot-address knots-reg end-knot-offset)
(l.s end-knot end-knot-address)
(c.eq.s u end-knot)
(b.fpt curve-basis-ready :delay (m span last-span))
;; For [0 0 0 0 1 2 3 ...], u in [k, k+1) belongs to span int(u) + 3.
;; Validate that candidate so fractional or nonuniform knots fall back to binary search.
(cvt.w.s rounded-u-float u)
(m rounded-u rounded-u-float)
(set! candidate-span (+ rounded-u 3))
(.dsll candidate-offset candidate-span 2)
(.daddu candidate-address knots-reg candidate-offset)
(l.s candidate-low candidate-address)
(l.s candidate-high candidate-address 4)
(c.le.s candidate-low u)
(b.fpf curve-span-binary-search :delay (nop!))
(c.lt.s u candidate-high)
(b.fpf curve-span-binary-search :delay (m span candidate-span))
(b curve-basis-ready :delay (nop!))
(label curve-span-binary-search)
(m low-span low-span)
(set! high-span (+ last-span 1))
(label curve-span-search-loop)
(set! mid-span (sar (+ low-span high-span) 1))
(.dsll mid-offset mid-span 2)
(.daddu mid-address knots-reg mid-offset)
(l.s mid-low mid-address)
(c.lt.s u mid-low)
(b.fpf curve-span-check-upper :delay (nop!))
(b curve-span-search-loop :delay (m high-span mid-span))
(label curve-span-check-upper)
(l.s mid-high mid-address 4)
(c.lt.s u mid-high)
(b.fpt curve-span-found :delay (nop!))
(b curve-span-search-loop :delay (m low-span mid-span))
(label curve-span-found)
(m span mid-span)
(b.eq span rounded-u curve-basis-ready :delay (nop!))
(nop!)
(nop!)
(label curve-basis-ready)
(m! basis-fn calculate-basis-functions-vector!)
(m a0 basis)
(m a1 span)
(m a2 u-bits)
(m a3 knots-reg)
(jalr ra basis-fn :delay (.sll v0 ra 0))
;; Four basis weights select controls span-3 through span. The no-ops retain the
;; VU load and multiply-accumulate order.
(.addiu first-control span -3)
(.lvf weights basis)
(.sll control-offset first-control 4)
(.add.x.vf.w curve-point vf0 vf0)
(.daddu control-address saved-controls control-offset)
(nop!)
(nop!)
(.lvf control-0 control-address)
(nop!)
(.lvf control-1 control-address 16)
(nop!)
(.lvf control-2 control-address 32)
(nop!)
(.lvf control-3 control-address 48)
(.mul.x.vf acc control-0 weights)
(nop!)
(.add.mul.y.vf.xyz acc control-1 weights acc)
(nop!)
(.add.mul.z.vf.xyz acc control-2 weights acc)
(nop!)
(.add.mul.w.vf.xyz curve-point control-3 weights acc)
(nop!)
(nop!)
(nop!)
(nop!)
(nop!)
(nop!)
(nop!)
(nop!)
(.svf saved-dst curve-point)
(m curve-point-bits curve-point)
(m result saved-dst)
(j ra :delay (m sp sp))
(nop!)
(nop!)))
File diff suppressed because it is too large Load Diff
+2 -2
View File
@@ -19,9 +19,9 @@
;; normally.
;; There is a bug in some of the cosine functions that can be fixed by toggling this flag.
;; The shared scalar sine/cosine routine constructs its x^2 cosine coefficient with only
;; The shared scalar sine/cosine function constructs its x^2 cosine coefficient with only
;; `lui #xbeff`, producing -0.498046875 instead of the intended #xbefffd62 value used by
;; cos-rad. The standalone cosine routines are unaffected.
;; cos-rad. The standalone cosine functions are unaffected.
(defglobalconstant FIX_COSINE_BUG #f)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;