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) } } }