Божин обнови решението на 18.11.2018 23:41 (преди 9 месеца)
+package main
+
+import (
+ "math"
+
+ "github.com/fmi/go-homework/geom"
+)
+
+func main() {}
+
+func equal(a, b float64) bool {
+ return math.Abs(a-b) < 1e-6
+}
+
+func Det(m [][]float64) float64 {
+ diag1 := m[0][0]*m[1][1]*m[2][2] +
+ m[0][1]*m[1][2]*m[2][0] +
+ m[0][2]*m[1][0]*m[2][1]
+
+ diag2 := m[2][0]*m[1][1]*m[0][2] +
+ m[2][1]*m[1][2]*m[0][0] +
+ m[2][2]*m[1][0]*m[0][1]
+
+ return diag1 - diag2
+}
+
+func Len(v geom.Vector) float64 {
+ return math.Sqrt(v.X*v.X + v.Y*v.Y + v.Z*v.Z)
+}
+
+func Add(v, w geom.Vector) geom.Vector {
+ return geom.Vector{v.X + w.X, v.Y + w.Y, v.Z + w.Z}
+}
+
+func Sub(v, w geom.Vector) geom.Vector {
+ return geom.Vector{v.X - w.X, v.Y - w.Y, v.Z - w.Z}
+}
+
+func Mul(v geom.Vector, s float64) geom.Vector {
+ return geom.Vector{s * v.X, s * v.Y, s * v.Z}
+}
+
+func Dot(v, w geom.Vector) float64 {
+ return v.X*w.X + v.Y*w.Y + v.Z*w.Z
+}
+
+func Cross(v, w geom.Vector) geom.Vector {
+ return geom.Vector{
+ v.Y*w.Z - v.Z*w.Y,
+ v.Z*w.X - v.X*w.Z,
+ v.X*w.Y - v.Y*w.X,
+ }
+}
+
+func sameSign(a, b float64) bool {
+ return a*b >= 0
+}
+
+func RayContainsPoint(ray geom.Ray, point geom.Vector) bool {
+ // The point must lie on the ray.
+
+ return sameSign(ray.Direction.X, point.X-ray.Origin.X) &&
+ sameSign(ray.Direction.Y, point.Y-ray.Origin.Y) &&
+ sameSign(ray.Direction.Z, point.Z-ray.Origin.Z)
+}
+
+type Plane struct {
+ triangle Triangle
+ A, B, C, D float64
+ offset geom.Vector
+}
+
+func NewPlane(t Triangle) Plane {
+ D := Det([][]float64{
+ {t.A.X, t.A.Y, t.A.Z},
+ {t.B.X, t.B.Y, t.B.Z},
+ {t.C.X, t.C.Y, t.C.Z},
+ })
+
+ var offset geom.Vector
+ if D == 0 {
+ offset = geom.NewVector(100, 100, 100)
+ t = NewTriangle(
+ Add(t.A, offset),
+ Add(t.B, offset),
+ Add(t.C, offset),
+ )
+ D = Det([][]float64{
+ {t.A.X, t.A.Y, t.A.Z},
+ {t.B.X, t.B.Y, t.B.Z},
+ {t.C.X, t.C.Y, t.C.Z},
+ })
+ } else {
+ offset = geom.NewVector(0, 0, 0)
+ }
+
+ d := float64(1)
+
+ a := -d / D * Det([][]float64{
+ {1, t.A.Y, t.A.Z},
+ {1, t.B.Y, t.B.Z},
+ {1, t.C.Y, t.C.Z},
+ })
+
+ b := -d / D * Det([][]float64{
+ {t.A.X, 1, t.A.Z},
+ {t.B.X, 1, t.B.Z},
+ {t.C.X, 1, t.C.Z},
+ })
+
+ c := -d / D * Det([][]float64{
+ {t.A.X, t.A.Y, 1},
+ {t.B.X, t.B.Y, 1},
+ {t.C.X, t.C.Y, 1},
+ })
+
+ return Plane{
+ triangle: t,
+ A: a,
+ B: b,
+ C: c,
+ D: d,
+ offset: offset,
+ }
+}
+
+func (p Plane) Normal() geom.Vector {
+ return Cross(
+ Sub(p.triangle.C, p.triangle.A),
+ Sub(p.triangle.B, p.triangle.A),
+ )
+}
+
+func (p Plane) Contains(point geom.Vector) bool {
+ point = Add(point, p.offset)
+ return equal(
+ 0,
+ p.A*point.X+
+ p.B*point.Y+
+ p.C*point.Z+
+ p.D,
+ )
+}
+
+// Intersections returns the number of points of intersection.
+// If the intersection is exactly one, it is also returned.
+func (p Plane) Intersections(ray geom.Ray) (int, geom.Vector) {
+ ray = geom.NewRay(
+ Add(ray.Origin, p.offset),
+ ray.Direction,
+ )
+
+ normal := p.Normal()
+
+ diff := Sub(ray.Origin, p.triangle.A)
+ prod1 := Dot(diff, normal)
+ prod2 := Dot(ray.Direction, normal)
+
+ if equal(prod2, 0) {
+ if p.Contains(Sub(ray.Origin, p.offset)) {
+ return 2, geom.Vector{} // The ray lies on the plane.
+ }
+
+ return 0, geom.Vector{} // The ray is parallel to the plane.
+ }
+
+ intersection := Sub(ray.Origin, Mul(ray.Direction, prod1/prod2))
+ return 1, Sub(intersection, p.offset)
+}
+
+type Triangle struct {
+ A, B, C geom.Vector
+}
+
+func NewTriangle(a, b, c geom.Vector) Triangle {
+ return Triangle{a, b, c}
+}
+
+func (t Triangle) Area() float64 {
+ return Len(
+ Cross(
+ Sub(t.B, t.A),
+ Sub(t.C, t.A),
+ ),
+ ) / 2
+}
+
+func (t Triangle) Plane() Plane {
+ return NewPlane(t)
+}
+
+func (t Triangle) Contains(point geom.Vector) bool {
+ if !t.Plane().Contains(point) {
+ return false
+ }
+
+ ta := NewTriangle(point, t.B, t.C)
+ tb := NewTriangle(point, t.A, t.C)
+ tc := NewTriangle(point, t.A, t.B)
+ return equal(t.Area(), ta.Area()+tb.Area()+tc.Area())
+}
+
+func (t Triangle) Intersect(ray geom.Ray) bool {
+ count, intersection := t.Plane().Intersections(ray)
+
+ if count == 0 {
+ return false
+ }
+
+ if count == 1 {
+ return RayContainsPoint(ray, intersection) && t.Contains(intersection)
+ }
+
+ // Don't hate the player.
+ for quot := float64(0); !equal(1, quot); quot += 1e-4 {
+ if t.Contains(Add(
+ ray.Origin,
+ Mul(ray.Direction, quot),
+ )) {
+ return true
+ }
+ }
+
+ return false
+}
+
+type Quad struct {
+ a, b, c, d geom.Vector
+}
+
+func NewQuad(a, b, c, d geom.Vector) Quad {
+ return Quad{a, b, c, d}
+}
+
+func (q Quad) Intersect(ray geom.Ray) bool {
+ return NewTriangle(q.a, q.b, q.c).Intersect(ray) ||
+ NewTriangle(q.a, q.b, q.d).Intersect(ray) ||
+ NewTriangle(q.a, q.c, q.d).Intersect(ray) ||
+ NewTriangle(q.b, q.c, q.d).Intersect(ray)
+}
+
+type Sphere struct {
+ Origin geom.Vector
+ R float64
+}
+
+func NewSphere(origin geom.Vector, r float64) Sphere {
+ return Sphere{origin, r}
+}
+
+func (s Sphere) Contains(point geom.Vector) bool {
+ return Len(Sub(point, s.Origin)) <= s.R
+}
+
+// Closest returns the point from `ray`'s line, closest to `s`'s origin.
+func (s Sphere) Closest(ray geom.Ray) geom.Vector {
+ rayToSphere := Sub(s.Origin, ray.Origin)
+ quot := Dot(rayToSphere, ray.Direction) / Dot(ray.Direction, ray.Direction)
+ return Add(ray.Origin, Mul(ray.Direction, quot))
+}
+
+func (s Sphere) Intersect(ray geom.Ray) bool {
+ h := s.Closest(ray)
+
+ if !s.Contains(h) {
+ return false
+ }
+
+ return RayContainsPoint(ray, h) ||
+ s.Contains(ray.Origin) ||
+ s.Contains(Add(ray.Origin, ray.Direction))
+}
Трябва да знаеш, че лъчите са безкрайно дълги. Твоята имплементация тук всъщност проверява за пресичане между сфера и права.
Всъщност сега виждам че 178-я ред е излишен, но ми се струва, че и с него и без него имплементацията е коректна