diff options
Diffstat (limited to 'src/Objects.scala')
-rw-r--r-- | src/Objects.scala | 359 |
1 files changed, 356 insertions, 3 deletions
diff --git a/src/Objects.scala b/src/Objects.scala index b4e6435..a74ccc2 100644 --- a/src/Objects.scala +++ b/src/Objects.scala @@ -31,10 +31,26 @@ object Objects { 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) { + 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) @@ -117,10 +133,344 @@ object Objects { case (x: ArrayIndexOutOfBoundsException) => { /* well, we're not properly rendering a region so uh just ignore the failure i guess lol */ } } } - def normal: Segment = { + def normal(t: Double): Ray = { val normalMag = Math.sqrt(x * x + y * y) val finalNormMult = 1.5 / normalMag - Ray(-y * finalNormMult, x * finalNormMult, (at(0) + at(1)) / 2).toSegment + 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 { @@ -147,5 +497,8 @@ object Objects { def mag: Double = { Math.sqrt(x * x + y * y) } + def at(t: Double): Point = { + initial + Point(x * t, y * t) + } } } |