Траян обнови решението на 19.11.2018 10:24 (преди 9 месеца)
+package main
+
+import (
+ "github.com/fmi/go-homework/geom"
+ "math"
+)
+
+// type plane with some useful methods
+type plane struct {
+ A float64
+ B float64
+ C float64
+ D float64
+}
+
+func newPlane(triangle Triangle) plane {
+ vectAB := vectorByTwoPoints(triangle.A, triangle.B)
+ vectBC := vectorByTwoPoints(triangle.A, triangle.C)
+
+ normalizedVector := geom.NewVector(vectAB.Y*vectBC.Z-vectAB.Z*vectBC.Y,
+ vectAB.Z*vectBC.X-vectAB.X*vectBC.Z,
+ vectAB.X*vectBC.Y-vectAB.Y*vectBC.X)
+
+ dValue := -(normalizedVector.X*triangle.A.X +
+ normalizedVector.Y*triangle.A.Y +
+ normalizedVector.Z*triangle.A.Z)
+
+ plane := plane{
+ A: normalizedVector.X,
+ B: normalizedVector.Y,
+ C: normalizedVector.Z,
+ D: dValue,
+ }
+ return plane
+}
+
+func (plane *plane) rayParallelToPlane(ray geom.Ray) bool {
+ return (plane.A*ray.Direction.X +
+ plane.B*ray.Direction.Y +
+ plane.C*ray.Direction.Z) == 0.0
+}
+
+type Triangle struct {
+ A geom.Vector
+ B geom.Vector
+ C geom.Vector
+}
+
+func NewTriangle(a, b, c geom.Vector) Triangle {
+ return Triangle{
+ a,
+ b,
+ c,
+ }
+}
+
+func (triangle Triangle) Intersect(ray geom.Ray) bool {
+ trianglePlane := newPlane(triangle)
+ if trianglePlane.rayParallelToPlane(ray) {
+ return false
+ } else {
+ numerator := -(trianglePlane.A*ray.Origin.X +
+ trianglePlane.B*ray.Origin.Y +
+ trianglePlane.C*ray.Origin.Z) +
+ trianglePlane.D
+ denominator := trianglePlane.A*ray.Direction.X +
+ trianglePlane.B*ray.Direction.Y +
+ trianglePlane.C*ray.Direction.Z
+
+ t := numerator / denominator
+
+ if t < 0 {
+ return false
+ }
+ intersectionPoint := geom.NewVector(ray.Origin.X+t*ray.Direction.X,
+ ray.Origin.Y+t*ray.Direction.Y,
+ ray.Origin.Z+t*ray.Direction.Z)
+
+ triangleAPC := NewTriangle(triangle.A, intersectionPoint, triangle.C)
+ triangleABP := NewTriangle(triangle.A, triangle.B, intersectionPoint)
+ trianglePBC := NewTriangle(intersectionPoint, triangle.B, triangle.C)
+
+ baricentPBCtoABC := trianglePBC.area() / triangle.area()
+ baricentAPCtoABC := triangleAPC.area() / triangle.area()
+ baricentABPtoABC := triangleABP.area() / triangle.area()
+
+ if baricentPBCtoABC >= 0.0 && baricentPBCtoABC <= 1.0 {
+ if baricentAPCtoABC >= 0.0 && baricentAPCtoABC <= 1.0 {
+ if baricentABPtoABC >= 0.0 && baricentABPtoABC <= 1.0 {
+ return true
+ }
+ }
+ }
+
+ return false
+ }
+}
+
+func (triangle *Triangle) area() float64 {
+ vectAB := vectorByTwoPoints(triangle.A, triangle.B)
+ vectBC := vectorByTwoPoints(triangle.A, triangle.C)
+
+ normalizedVector := geom.NewVector(vectAB.Y*vectBC.Z-vectAB.Z*vectBC.Y,
+ vectAB.Z*vectBC.X-vectAB.X*vectBC.Z,
+ vectAB.X*vectBC.Y-vectAB.Y*vectBC.X)
+
+ return math.Sqrt(normalizedVector.X*normalizedVector.X+
+ normalizedVector.Y*normalizedVector.Y+
+ normalizedVector.Z*normalizedVector.Z) / 2
+}
+
+type Quadrangle struct {
+ triangle Triangle
+ D geom.Vector
+}
+
+func NewQuad(a, b, c, d geom.Vector) Quadrangle {
+ return Quadrangle{
+ NewTriangle(a, b, c),
+ d,
+ }
+}
+
+func (quadrangle Quadrangle) Intersect(ray geom.Ray) bool {
+ secondTriangle := NewTriangle(quadrangle.triangle.A, quadrangle.triangle.C, quadrangle.D)
+ return quadrangle.triangle.Intersect(ray) ||
+ secondTriangle.Intersect(ray)
+}
+
+// type quadratic equation with some useful methods
+type quadrEquation struct {
+ a float64
+ b float64
+ c float64
+}
+
+func newQuadrEqation(a, b, c float64) quadrEquation {
+ return quadrEquation{
+ a,
+ b,
+ c,
+ }
+}
+
+func (qdeq *quadrEquation) discriminant() float64 {
+ return qdeq.b*qdeq.b - 4*qdeq.a*qdeq.c
+}
+
+type Sphere struct {
+ center geom.Vector
+ radius float64
+}
+
+func NewSphere(origin geom.Vector, r float64) Sphere {
+ return Sphere{
+ origin,
+ r,
+ }
+}
+func (sphere Sphere) Intersect(ray geom.Ray) bool {
+ sphereQuadrEq := newQuadrEqation(
+ scalarMultiplication(ray.Direction, ray.Direction),
+ 2*scalarMultiplication(
+ ray.Direction,
+ subtractVectors(ray.Origin, sphere.center)),
+ scalarMultiplication(
+ subtractVectors(ray.Origin, sphere.center),
+ subtractVectors(ray.Origin, sphere.center))-
+ sphere.radius*sphere.radius)
+
+ det := sphereQuadrEq.discriminant()
+ if det < 0 {
+ return false
+ }
+ if det == 0 {
+ t := -sphereQuadrEq.b / 2
+ if t >= 0 {
+ return true
+ }
+ return false
+ }
+
+ t1 := (-sphereQuadrEq.b + math.Sqrt(det)) / 2
+ t2 := (-sphereQuadrEq.b - math.Sqrt(det)) / 2
+ if t1 >= 0 || t2 >= 0 {
+ return true
+ }
+
+ return false
+}
+
+// helper functions for vector operations
+func vectorByTwoPoints(point1 geom.Vector, point2 geom.Vector) geom.Vector {
+ return geom.NewVector(point2.X-point1.X, point2.Y-point1.Y, point2.Z-point1.Z)
+}
+
+func subtractVectors(v1 geom.Vector, v2 geom.Vector) geom.Vector {
+ return geom.NewVector(v1.X-v2.X, v1.Y-v2.Y, v1.Z-v2.Z)
+}
+
+func scalarMultiplication(v1 geom.Vector, v2 geom.Vector) float64 {
+ return v1.X*v2.X + v1.Y*v2.Y + v1.Z*v2.Z
+}