path: root/src/Objects.scala

package net.iximeow.raytrace

import java.awt.image.BufferedImage

object Objects {
  /*
   * 3d nah
  case class Point(x: Double, y: Double, z: Double)
  case class Plane(pitch: Double, roll: Double, altitude: Double)
  case class BoundedPlane(pitch: Double, roll: Double, center: Point, h: Point, w: Point)
  */

  case class Point(x: Double, y: Double) {
    def +(other: Point): Point = Point(x + other.x, y + other.y)
    def -(other: Point): Point = Point(x - other.x, y - other.y)
    def /(scale: Double): Point = Point(x / scale, y / scale)
    def *(scale: Double): Point = Point(x * scale, y * scale)
    def magnitude: Double = distTo(Point.Zero)
    def distTo(other: Point) = Math.sqrt(sqDistTo(other))
    def sqDistTo(other: Point) = (other.x - x) * (other.x - x) + (other.y - y) * (other.y - y)
    def dot(other: Point): Double = x * other.x + y * other.y
  }
  object Point {
    val Zero = Point(0, 0)
  }
  case class Line(m: Double, b: Double)
  object Line {
    def fromPoints(p1: Point, p2: Point): Line = {
      val m = (p2.y - p1.y) / (p2.x - p1.x)
      val b = p1.y - m*p1.x
      Line(m, b)
    }
  }
  trait Surface {
    // returns a t where other first intersets with this surface
    def intersect(other: Segment): Double
    // intersect, but returns None if intersection (t) is outside [0,1]
    def intersectChecked(other: Segment): Option[Point]
    // surface is some parametric function Double => Point
    def at(t: Double): Point
    // aka at^-1, inverse of at(t). Option, for cases that p is not on this surface.
    def tFor(p: Point): Option[Double]
    // reflect a ray off this surface, assuming intersection at point
    //   returns a truncated form of source stopped at intersection
    //   and a continuation after interaction
    def scatter(source: Ray, intersection: Point): (Ray, Ray)
    def normal(t: Double): Ray
    def renderTo(buf: BufferedImage, scale: Double, xoff: Int, yoff: Int, color: Int): Unit
  }
  /*
   * Segments are defined for t in [0, 1]
   */
  case class Segment(x: Double, y: Double, initial: Point) extends Surface {
    def at(t: Double): Point = Point(x * t, y * t) + initial
    def apply = at _
    def length = Math.sqrt(x * x + y * y)
    def intersect(other: Segment): Double = {
      /*
       *  P1 = ai + t * a
       *  P2 = bi + u * b
       *  P1 == P2, so
       *  ai + t * a = bi + u * b
       *  ...
       *  ai.x + t * a.x = bi.x + u * b.x
       *  ai.y + t * a.y = bi.y + u * b.y
       *  t = (bi.x + u * b.x - ai.x) / a.x
       *  ai.y + a.y * (bi.x + u * b.x - ai.x) / a.x = bi.y + u * b.y
       *  ai.y + a.y * bi.x / a.x + u * b.x * a.y / a.x - ai.x * a.y / a.x = bi.y + u * b.y
       *  ai.y - bi.y + (a.y / a.x) (bi.x - ai.x) = u * b.y - u * b.x * a.y / a.x
       *  ai.y - bi.y + (a.y / a.x) (bi.x - ai.x) / (b.y - b.x * a.y / a.x) = u
       *
       *  (a.x * (ai.y - bi.y) + a.y * (bi.x - ai.x)) / (b.y - b.x * a.y)
       */
      val u = (
        x * (initial.y - other.initial.y) +
        y * (other.initial.x - initial.x)
      ) / (
        other.y * x - other.x * y
      )

      u
    }

    def intersectChecked(other: Segment): Option[Point] = {
      val u = intersect(other)
      if (u >= 0 && u <= 1) {
        //println("Intersection is at u=" + u)

        Some(other.at(u))
      } else {
        None
      }
    }
    def tFor(p: Point): Option[Double] = {
      (x, y) match {
        case (0, _) => Some((p.y - initial.y) / y)
        case (_, 0) => Some((p.x - initial.x) / x)
        case (_, _) => {
          val xT = (p.x - initial.x) / x
          val yT = (p.y - initial.y) / y
          if (Math.abs(xT - yT) < 0.000001) {
            Some(xT)
          } else {
            None
          }
        }
      }
    }
    def rotate(angle: Double, about: Point = Point(0, 0)): Segment = {
      val start = this.at(0)
      val end = this.at(1)
      val newStart = {
        val offset = start - about
        val m = offset.magnitude
        val newAngle = Math.atan2(offset.y, offset.x) + angle
        Point(Math.cos(newAngle) * m, Math.sin(newAngle) * m) + about
      }
      val newEnd = {
        val offset = end - about
        val m = offset.magnitude
        val newAngle = Math.atan2(offset.y, offset.x) + angle
        Point(Math.cos(newAngle) * m, Math.sin(newAngle) * m) + about
      }
      Segment.fromPoints(newStart, newEnd)
    }
    def renderTo(buf: BufferedImage, scale: Double = 1, xoff: Int = 0, yoff: Int = 0, color: Int = 0x808000): Unit = {
      try {
        for (i <- (0 to 100)) {
          val point = this.at(i / 100.0)
          buf.setRGB(Math.round((point.x * scale).toFloat) + xoff, Math.round((point.y * scale).toFloat) + yoff, color)
        }
      } catch {
        case (x: ArrayIndexOutOfBoundsException) => { /* well, we're not properly rendering a region so uh just ignore the failure i guess lol */ }
      }
    }
    def normal(t: Double): Ray = {
      val normalMag = Math.sqrt(x * x + y * y)
      val finalNormMult = 1.5 / normalMag
      Ray(-y * finalNormMult, x * finalNormMult, at(t))
    }
    def scatter(r: Ray, firstIntersection: Point): (Ray, Ray) = {
      def isBehind(start: Ray): Boolean = {
        val normal = Ray(this.normal(0.5).x, this.normal(0.5).y, Point(0, 0))
        val rebased = Ray(start.x, start.y, Point(0, 0))
        val cosAngle = normal.dot(rebased) / (normal.mag * rebased.mag)
        cosAngle > 0
      }
      if (isBehind(r)) { // as in, 'r is behind this'
        // stop.
        (r.endingAt(firstIntersection), Ray(0, 0, r.initial))
      } else {
        // reflect.
        val minAngle = {
          val fromStart = Raymath.angleBetween(
            r.initial,
            firstIntersection,
            this.at(0)
          )
          val fromEnd = Raymath.angleBetween(
            r.initial,
            firstIntersection,
            this.at(1)
          )

          if (Math.abs(fromStart) < Math.PI / 2) {
            fromStart
          } else {
            fromEnd
          }

          fromStart
        }

        val maxAngle = Math.PI - minAngle

        val baseAngle = Math.atan2(this.y, this.x)

        val reflectedAngle = baseAngle + minAngle

        if (minAngle < 0 || minAngle > Math.PI * 2) {
          (r.endingAt(firstIntersection), r.endingAt(firstIntersection)) //Ray(0, 0, firstIntersection._2))
        } else {
          val (x, y) = (
            Math.cos(reflectedAngle) * 3,
            Math.sin(reflectedAngle) * 3
          )

          // Sure hope this is right...
          (r.endingAt(firstIntersection), Ray(x, y, firstIntersection))
        }
      }
    }
  }
  case class ParabolicLens(center: Point, rotation: Double, radius: Double, rMinor: Double) extends Surface {
    // assume we can just use parabolic mirror equations here...
    // 4FD = R^2, F = focal length, D = depth, R = radius
    // so we know the intended radius and focal length, time to derive D..
    val depth = radius * radius / (4 * rMinor)
    val b = radius * radius / (4 * rMinor)
    val bi = 0 // -depth
    /*
     *  P2_x = u * cos(rotation) - (bi + u^2 * b) * sin(rotation)
     *  P2_y = u * sin(rotation) + (bi + u^2 * b) * cos(rotation)
     */
    def normal(t: Double): Ray =
      normalRaw(t - 0.5)

    def normalRaw(t: Double): Ray = {
      val cosRot = Math.cos(rotation)
      val sinRot = Math.sin(rotation)
      // p2_x' = -sin(rotation) - 2u*b*cos(rotation)
      // p2_y' = cos(rotation) - 2u*b*sin(rotation)
      val p2_xp = -sinRot - 2*(t) * b * cosRot
      val p2_yp = cosRot - 2*(t) * b * sinRot
      Ray(p2_xp, p2_yp, this.atRaw((t)))
    }
    def scatter(source: Ray, intersection: Point): (Ray, Ray) = {
      // next up...
      val rayEnd = intersection
      val t = tFor(intersection)
      //println(t)
      t.map(t => {
        val norm = normal(t)
        // so we rotate by the angle diff of source and normal? times refraction index?
        val angle = -(Math.PI - Raymath.angleBetween(
          source.initial,
          intersection,
          norm.at(1)
        )) * refractionIndex
        /*
         * rotation matrix:
         * cos(rot)  -sin(rot)
         * sin(rot)  +cos(rot)
         */
        val out = Ray(source.x * Math.cos(angle) - source.y * Math.sin(angle), source.x * Math.sin(angle) + source.y * Math.cos(angle), intersection)
        (source.endingAt(intersection), out)
      }).getOrElse((source, Ray(source.x, source.y, source.at(1))))
      //???
    }
    def tFor(p: Point): Option[Double] = {
      val cosRot = Math.cos(rotation)
      val sinRot = Math.sin(rotation)
      // P2_x u: = center.x + u * cos(rotation) - (bi + u^2 * b) * sin(rotation)
      //       0 = center.x - P2_x + u * cos(rotation) - (bi + u^2 * b) * sin(rotation)
      //       0 = center.x - P2_x - bi * sin(rotation) + u * cos(rotation) - u^2 * b * sin(rotation)
      // P2_y u: = center.y + u * sin(rotation) + (bi + u^2 * b) * cos(rotation)
      //       0 = center.y - P2_y + bi * cos(rotation) + u * sin(rotation) + u^2 * b * cos(rotation)
      val x_us = if (Math.abs(sinRot) > 0.0000001) {
        quadradicRoots(-b * sinRot, cosRot, center.x - p.x - bi * sinRot)
      } else {
        Seq(-(center.x - p.x) / cosRot / 2) // why does this appear to be off by a factor of 2?
      }
      val y_us = if (Math.abs(cosRot) > 0.0000001) {
        quadradicRoots(b * cosRot, sinRot, center.y - p.y + bi * cosRot)
      } else {
        Seq(-(center.y - p.y) / sinRot / 2)
      }
      //println(s"x_us: $x_us and y_us: $y_us")
      val matches = for {
        xu <- x_us
        yu <- y_us
      } yield {
        if (Math.abs(xu - yu) < 0.00001) {
          Some(xu)
        } else {
          None
        }
      }
      matches.flatten match {
        case Seq() => None
        case Seq(u) => Some(u + 0.5)
        case Seq(a, b) => {
          if (Math.abs(a - b) < 0.000001) {
            Some(a + 0.5)
          } else {
            throw new Exception("Too many t's for point " + p + " (" + a + ", " + b + ")")
          }
        }
        case x => throw new Exception("Too many t's for point " + p + " (" + x + ")")
      }
    }
    def renderTo(buf: BufferedImage, scale: Double = 1, xoff: Int = 0, yoff: Int = 0, color: Int = 0x205080): Unit = {
      var i = 0.0
      while (i <= 1.0) {
        val point = this.at(i)
        try {
          buf.setRGB(Math.round((point.x * scale).toFloat) + xoff, Math.round((point.y * scale).toFloat) + yoff, color)
        } catch {
          case e: Exception => { }
        }
        i = i + 0.01 // 100 points
      }
    }

    def at(t: Double): Point = atRaw(t - 0.5)
    def atRaw(t: Double): Point = {
      val cosRot = Math.cos(rotation)
      val sinRot = Math.sin(rotation)
      val t2 = t * t
      val bit2b = bi + t2 * b
      Point(t * cosRot * radius - bit2b * sinRot * radius, t * sinRot + bit2b * cosRot) + center
    }

    val refractionIndex = 1.52 // 1.52-1.75
    def intersect(other: Segment): Double = {
      /*
       *    WRONG:
       *  px = other.x * t + other.initial.x
       *  py = other.y * t + other.initial.y
       *  px = other.x 
       *  P1_x = t
       *  P1_y = ai + t * a
       *  P2_x = u * cos(rotation)
       *  p2_y = (bi + u^2 * b) * sin(rotation)
       *  ai + t * a = (bi + u^2 * b) * sin(rotation)
       *  0 = bi * sin(rotation) - ai - t * a + (t / cos(rotation))^2 * b * sin(rotation) // t == u by P1_x == P2_x
       *  0 = bi * sin(rotation) - ai - t * a + t ^ 2 * b * sin(rotation) / cos(rotation) ^2
       *
       *  t = a +- sqrt(a^2 - 4 (b * sin(rotation) / cos(rotation)^2  * (bi * sin(rotation) - ai))) / 2(b * sin(rotation) / cos(rotation) ^2)
       *
       *
       *    WRONG:
       *      definition of p2_{x,y} is wrong for rotation.
       *  px = other.x * t + other.initial.x
       *  py = other.y * t + other.initial.y
       *  P1_x = axi + t * ax
       *  P1_y = ayi + t * ay
       *  P2_x = u * cos(rotation)
       *  p2_y = (bi + u^2 * b) * sin(rotation)
       *
       *    P1_y = P2_y
       *  ayi + t * ay = (bi + u^2 * b) * sin(rotation)
       *
       *    P1_x = P2_x
       *  axi + t * ax = u * cos(rotation)
       *  (axi + t * ax) / cos(rotation) = u
       *
       *    sub u for t to have one variable to solve for
       *      reminder: t is parameter for `other` aka P1
       *  ayi + t * ay = (bi + ((axi + t * ax) / cos(rotation))^2 * b) * sin(rotation)
       *  ayi + t * ay = (bi + (axi + t * ax)^2 / cos(rotation)^2 * b) * sin(rotation)
       *  ayi / sin(rotation) + t * ay / sin(rotation) = bi + (axi + t * ax)^2 / cos(rotation)^2 * b
       *  ayi / (b * sin(rotation)) + t * a / (b * sin(rotation)) = bi / b + (axi + t * ax)^2 / cos(rotation)^2
       *  ayi / (b * sin(rotation)) - bi / b + t * ay / (b * sin(rotation)) = (axi + t * ax)^2 / cos(rotation)^2
       *  (ai - bi * sin(rotation)) / (b * sin(rotation)) + t * a / (b * sin(rotation)) = (axi + t * ax)^2 / cos(rotation)^2
       *  (cos(rotation)^2 * (ayi - bi * sin(rotation))) / (b * sin(rotation)) + t * ay * cos(rotation)^2 / (b * sin(rotation)) = (axi + t * ax)^2
       *  (cos(rotation)^2 * (ayi - bi * sin(rotation))) / (b * sin(rotation)) + t * ay * cos(rotation)^2 / (b * sin(rotation)) = axi^2 + 2 * ax * axi * t + (t * ax)^2
       *  (cos(rotation)^2 * (ayi - bi * sin(rotation))) / (b * sin(rotation)) - axi^2 - 2 * ax * axi * t + t * a * cos(rotation)^2 / (b * sin(rotation)) - t^2 * ax^2 = 0
       *  (cos(rotation)^2 * (ayi - bi * sin(rotation))) / (b * sin(rotation)) - axi^2 +
       *    - 2 * ax * axi * t + t * ay * cos(rotation)^2 / (b * sin(rotation)) +
       *    - t^2 * ax^2
       *  (cos(rotation)^2 * (ayi - bi * sin(rotation))) / (b * sin(rotation)) - axi^2 +
       *    t * (ay * cos(rotation)^2 - 2 * ax * axi * b * sin(rotation)) / (b * sin(rotation)) +
       *    - t^2 * ax^2
       *
       *  lol quadradic
       *
       * rotation matrix:
       * cos(rot)  -sin(rot)
       * sin(rot)  +cos(rot)
       *
       *    RIGHT..ish:
       *      doesn't account for possibility of non-centered parabola. FIX: subtract axi and ayi by center.x and center.y.
       *      this SHOULD be what the math shows anyway...
       *  px = other.x * t + other.initial.x
       *  py = other.y * t + other.initial.y
       *  P1_x = axi + t * ax
       *  P1_y = ayi + t * ay
       *  P2_x = u * cos(rotation) * r - (bi + u^2 * b) * sin(rotation) * r
       *  P2_y = u * sin(rotation) + (bi + u^2 * b) * cos(rotation)
       *
       *  so, intersect if P1_y(t1) == P2_y(u1) and P1_x(t1) == P2_x(u1)
       *    P1_x = P2_x:
       *  axi + t * ax = u * cos(rotation) * r - (bi + u^2 * b) * sin(rotation) * r
       *  t = (u * cos(rotation) * r - (bi + u^2 * b) * sin(rotation) * r - axi) / ax
       *    P1_y = P2_y:
       *  ayi + t * ay = u * sin(rotation) + (bi + u^2 * b) * cos(rotation)
       *    sub t to only solve for u
       *  ayi + (u * cos(rotation) * r - (bi + u^2 * b) * sin(rotation) * r - axi) / ax * ay = u * sin(rotation) + (bi + u^2 * b) * cos(rotation)
       *  ayi * ax / ay + u * cos(rotation) * r- (bi + u^2 * b) * sin(rotation) * r - axi = u * sin(rotation) * ax / ay + (bi + u^2 * b) * cos(rotation) * ax / ay
       *  ayi * ax / ay - axi + u * cos(rotation) * r - (bi + u^2 * b) * sin(rotation) * r = u * sin(rotation) * ax / ay + (bi + u^2 * b) * cos(rotation) * ax / ay
       *  ayi * ax / ay - axi + u * cos(rotation) * r - u * sin(rotation) * ax / ay - bi * sin(rotation) * r - u^2 * b * sin(rotation) * r = bi * cos(rotation) * ax / ay + u^2 * b * cos(rotation) * ax / ay
       *  ayi * ax / ay - axi - bi * cos(rotation) * ax / ay - bi * sin(rotation) * r + u * (cos(rotation) * r - sin(rotation) * ax / ay) - u^2 * b * sin(rotation) * r = u^2 * b * cos(rotation) * ax / ay
       *  (ayi - bi * cos(rotation)) * ax / ay - axi - bi * sin(rotation) * r  + u * (cos(rotation) * r - sin(rotation) * ax / ay) - u^2 * b * sin(rotation) * r - u^2 * b * cos(rotation) * ax / ay = 0
       *  (ayi - bi * cos(rotation)) * ax / ay - axi - bi * sin(rotation) * r + u * (cos(rotation) * r - sin(rotation) * ax / ay) - u^2 * b * (sin(rotation) * r + cos(rotation) * ax / ay) = 0
       *  (ayi - bi * cos(rotation)) * ax / ay - axi - bi * sin(rotation) * r +
       *    u * (cos(rotation) * r - sin(rotation) * ax / ay) +
       *    - u^2 * b * (sin(rotation) * r - cos(rotation) * ax / ay) = 0
       *
       *  but in cases like this where ax = 0...
       *  -axi - bi * sin(rotation) * r +
       *    u * cos(rotation) * r +
       *    - u^2 * b * sin(rotation) * r
       */
      val ax = other.x
      val axi = other.initial.x - center.x
      val ay = other.y
      val ayi = other.initial.y - center.y
      val cosRot = Math.cos(rotation)
      val cosRot2 = cosRot * cosRot
      val sinRot = Math.sin(rotation)

      val quad_a = b * (sinRot * radius + cosRot * ax / ay)
      val quad_b = cosRot * radius - sinRot * ax / ay
      val quad_c = (ayi - bi * cosRot) * ax / ay - axi - bi * sinRot * radius

      // potentially valid t for intersection
      val options: Seq[Double] = if (Math.abs(quad_a) > 0.000001) {
        quadradicRoots(quad_a, quad_b, quad_c)
      } else {
        // quad_a is basically 0, so...
        // quad_b * u + quad_c == 0
        Seq(-quad_c / quad_b)
      }

      /*
       * Returns all u such that P2_x(u) = p2_xi
       * p2_xi = u * cos(rotation) * r - (bi + u^2 * b) * sin(rotation) * r
       * p2_xi + bi * sin(rotation) * r - u * cos(rotation) * r + u*2 * b * sin(rotation) * r = 0
       */
      def `P2_x^-1`(p2_xi: Double): Seq[Double] = quadradicRoots(b * sinRot * radius, cosRot * radius, p2_xi + bi * sinRot * radius)
      def P1_x(t: Double): Double = other.at(t).x
      def `P1^-1`(pt: Point): Double = {
        val byX = (pt.x - other.initial.x) / other.x
        val byY = (pt.y - other.initial.y) / other.y
        if (Math.abs(other.x) <= 0.00001) {
          byY
        } else if (Math.abs(other.y) <= 0.00001) {
          byX
        } else {
          0
        }
      }

      /* returns Some(t) if there is an intersection there, None if there is not */
      // eg verifies that at other(t), u = P2_x^-1(other(t).x), that 
      val ts: Seq[Double] = options.filter(u => u >= -0.5 && u <= 0.5).flatMap(u => {
        val ourPoint = this.atRaw(u)
        val candidateT = `P1^-1`(ourPoint)
        //println("candidate t: " + candidateT)
        //println("ours: " + ourPoint + " their: " + other.at(candidateT))
        Some(candidateT).filter(t => ourPoint.distTo(other.at(t)) < 0.00001)
      })

      ts match {
        case Seq() => Double.NaN
        case Seq(t) => t
        case Seq(t1, t2) => Math.min(t1, t2)
      }
    }

    def quadradicRoots(a: Double, b: Double, c: Double): Seq[Double] = {
      val radical = b * b - 4 * a * c
      val bdiv2a = -b / (2 * a)
      if (radical < 0) {
        Nil
      } else if (radical == 0) {
        Seq(bdiv2a)
      } else {
        val sqrtRad = Math.sqrt(radical) / (2 * a)
        Seq(bdiv2a + sqrtRad, bdiv2a - sqrtRad)
      }
    }

    def intersectChecked(other: Segment): Option[Point] = {
      val u = intersect(other)
      //println("Intersection is at u=" + u)
      if (u >= 0 && u <= 1) {

        Some(other.at(u))
      } else {
        None
      }
    }
  }
  object Segment {
    def fromPoints(p1: Point, p2: Point): Segment = {
      val (x0, y0) = (p2.x - p1.x, p2.y - p1.y)
      Segment(x0, y0, p1)
    }
    def fromPointWithAngle(p: Point, angle: Double): Segment = {
      val (x0, y0) = (Math.cos(angle), Math.sin(angle))
      Segment(x0, y0, p)
    }
  }
  case class Ray(x: Double, y: Double, initial: Point) {
    def toSegment: Segment = Segment(x, y, initial)
    def endingAt(p: Point): Ray = {
      val intersectAt = this.toSegment.tFor(p)
      intersectAt.map(t =>
        Ray(x * t, y * t, initial)
      ).getOrElse(this)
    }
    def dot(other: Ray): Double = {
      x * other.x + y * other.y
    }
    def mag: Double = {
      Math.sqrt(x * x + y * y)
    }
    def at(t: Double): Point = {
      initial + Point(x * t, y * t)
    }
  }
}