From 7f26f5e5ded749100cff077f25058835f48614d6 Mon Sep 17 00:00:00 2001 From: Edwin Jakobs Date: Sat, 13 Jan 2024 23:11:17 +0100 Subject: [PATCH] [orx-shapes] Adopt code from openrndr-shape, openrndr-math --- .../src/jvmDemo/kotlin/DemoExtrude01.kt | 4 +- .../src/jvmDemo/kotlin/DemoExtrude02.kt | 4 +- .../src/commonMain/kotlin/offset/Offset.kt | 268 ++++++++++++++ .../src/commonMain/kotlin/simplify/Chaikin.kt | 93 +++++ .../kotlin/simplify/RamerDouglasPeucker.kt | 45 +++ .../commonMain/kotlin/splines/CatmullRom.kt | 333 ++++++++++++++++++ 6 files changed, 743 insertions(+), 4 deletions(-) create mode 100644 orx-shapes/src/commonMain/kotlin/offset/Offset.kt create mode 100644 orx-shapes/src/commonMain/kotlin/simplify/Chaikin.kt create mode 100644 orx-shapes/src/commonMain/kotlin/simplify/RamerDouglasPeucker.kt create mode 100644 orx-shapes/src/commonMain/kotlin/splines/CatmullRom.kt diff --git a/orx-mesh-generators/src/jvmDemo/kotlin/DemoExtrude01.kt b/orx-mesh-generators/src/jvmDemo/kotlin/DemoExtrude01.kt index c0ed6323..2eae562a 100644 --- a/orx-mesh-generators/src/jvmDemo/kotlin/DemoExtrude01.kt +++ b/orx-mesh-generators/src/jvmDemo/kotlin/DemoExtrude01.kt @@ -6,10 +6,10 @@ import org.openrndr.draw.shadeStyle import org.openrndr.extra.camera.Orbital import org.openrndr.extra.meshgenerators.buildTriangleMesh import org.openrndr.extra.meshgenerators.extrudeContourSteps +import org.openrndr.extra.shapes.splines.catmullRom +import org.openrndr.extra.shapes.splines.toPath3D import org.openrndr.math.Vector3 -import org.openrndr.math.catmullRom import org.openrndr.shape.Circle -import org.openrndr.shape.toPath3D fun main() { application { diff --git a/orx-mesh-generators/src/jvmDemo/kotlin/DemoExtrude02.kt b/orx-mesh-generators/src/jvmDemo/kotlin/DemoExtrude02.kt index dda3024d..64caa775 100644 --- a/orx-mesh-generators/src/jvmDemo/kotlin/DemoExtrude02.kt +++ b/orx-mesh-generators/src/jvmDemo/kotlin/DemoExtrude02.kt @@ -6,11 +6,11 @@ import org.openrndr.draw.shadeStyle import org.openrndr.extra.camera.Orbital import org.openrndr.extra.meshgenerators.buildTriangleMesh import org.openrndr.extra.meshgenerators.extrudeShapeSteps +import org.openrndr.extra.shapes.splines.catmullRom +import org.openrndr.extra.shapes.splines.toPath3D import org.openrndr.math.Vector3 -import org.openrndr.math.catmullRom import org.openrndr.shape.Circle import org.openrndr.shape.Shape -import org.openrndr.shape.toPath3D fun main() { application { diff --git a/orx-shapes/src/commonMain/kotlin/offset/Offset.kt b/orx-shapes/src/commonMain/kotlin/offset/Offset.kt new file mode 100644 index 00000000..fa443c79 --- /dev/null +++ b/orx-shapes/src/commonMain/kotlin/offset/Offset.kt @@ -0,0 +1,268 @@ +package offset + +import org.openrndr.math.Vector2 +import org.openrndr.math.YPolarity +import org.openrndr.math.times +import org.openrndr.shape.* +import kotlin.math.abs +import kotlin.math.sign +import kotlin.math.sqrt + + +private fun Segment.splitOnExtrema(): List { + var extrema = extrema().toMutableList() + + if (isStraight(0.05)) { + return listOf(this) + } + + if (simple && extrema.isEmpty()) { + return listOf(this) + } + + if (extrema.isEmpty()) { + return listOf(this) + } + if (extrema[0] <= 0.01) { + extrema[0] = 0.0 + } else { + extrema = (mutableListOf(0.0) + extrema).toMutableList() + } + + if (extrema.last() < 0.99) { + extrema = (extrema + listOf(1.0)).toMutableList() + } else if (extrema.last() >= 0.99) { + extrema[extrema.lastIndex] = 1.0 + } + + return extrema.zipWithNext().map { + sub(it.first, it.second) + } +} + +private fun Segment.splitToSimple(step: Double): List { + var t1 = 0.0 + var t2 = 0.0 + val result = mutableListOf() + while (t2 <= 1.0) { + t2 = t1 + step + while (t2 <= 1.0 + step) { + val segment = sub(t1, t2) + if (!segment.simple) { + t2 -= step + if (abs(t1 - t2) < step) { + return listOf(this) + } + val segment2 = sub(t1, t2) + result.add(segment2) + t1 = t2 + break + } + t2 += step + } + + } + if (t1 < 1.0) { + result.add(sub(t1, 1.0)) + } + if (result.isEmpty()) { + result.add(this) + } + return result +} + + +fun Segment.reduced(stepSize: Double = 0.01): List { + val pass1 = splitOnExtrema() + //return pass1 + return pass1.flatMap { it.splitToSimple(stepSize) } +} + +fun Segment.scale(scale: Double, polarity: YPolarity) = scale(polarity) { scale } + +fun Segment.scale(polarity: YPolarity, scale: (Double) -> Double): Segment { + if (control.size == 1) { + return cubic.scale(polarity, scale) + } + + val newStart = start + normal(0.0, polarity) * scale(0.0) + val newEnd = end + normal(1.0, polarity) * scale(1.0) + + val a = LineSegment(newStart, start) + val b = LineSegment(newEnd, end) + + val o = intersection(a, b, 1E7) + + if (o != Vector2.INFINITY) { + val newControls = control.mapIndexed { index, it -> + val d = it - o + val rc = scale((index + 1.0) / 3.0) + val s = normal(0.0, polarity).dot(d).sign + val nd = d.normalized * s + it + rc * nd + } + return copy(newStart, newControls.toTypedArray(), newEnd) + } else { + val newControls = control.mapIndexed { index, it -> + val rc = scale((index + 1.0) / 3.0) + it + rc * normal((index + 1.0), polarity) + } + return copy(newStart, newControls.toTypedArray(), newEnd) + } +} + +fun Segment.offset( + distance: Double, + stepSize: Double = 0.01, + yPolarity: YPolarity = YPolarity.CW_NEGATIVE_Y +): List { + return if (linear) { + val n = normal(0.0, yPolarity) + if (distance > 0.0) { + listOf(Segment(start + distance * n, end + distance * n)) + } else { + val d = direction() + val s = distance.coerceAtMost(length / 2.0) + val candidate = Segment( + start - s * d + distance * n, + end + s * d + distance * n + ) + if (candidate.length > 0.0) { + listOf(candidate) + } else { + emptyList() + } + } + } else { + reduced(stepSize).map { it.scale(distance, yPolarity) } + } +} + + + +/** + * Offsets a [ShapeContour]'s [Segment]s by given [distance]. + * + * [Segment]s are moved outwards if [distance] is > 0 or inwards if [distance] is < 0. + * + * @param joinType Specifies how to join together the moved [Segment]s. + */ +fun ShapeContour.offset(distance: Double, joinType: SegmentJoin = SegmentJoin.ROUND): ShapeContour { + val offsets = + segments.map { it.offset(distance, yPolarity = polarity) } + .filter { it.isNotEmpty() } + val tempContours = offsets.map { + ShapeContour.fromSegments(it, closed = false, distanceTolerance = 0.01) + } + val offsetContours = tempContours.map { it }.filter { it.length > 0.0 }.toMutableList() + + for (i in 0 until offsetContours.size) { + offsetContours[i] = offsetContours[i].removeLoops() + } + + for (i0 in 0 until if (this.closed) offsetContours.size else offsetContours.size - 1) { + val i1 = (i0 + 1) % (offsetContours.size) + val its = intersections(offsetContours[i0], offsetContours[i1]) + if (its.size == 1) { + offsetContours[i0] = offsetContours[i0].sub(0.0, its[0].a.contourT) + offsetContours[i1] = offsetContours[i1].sub(its[0].b.contourT, 1.0) + } + } + + if (offsets.isEmpty()) { + return ShapeContour(emptyList(), false) + } + + + val startPoint = if (closed) offsets.last().last().end else offsets.first().first().start + + val candidateContour = contour { + moveTo(startPoint) + for (offsetContour in offsetContours) { + val delta = (offsetContour.position(0.0) - cursor) + val joinDistance = delta.length + if (joinDistance > 10e-6) { + when (joinType) { + SegmentJoin.BEVEL -> lineTo(offsetContour.position(0.0)) + SegmentJoin.ROUND -> arcTo( + crx = joinDistance * 0.5 * sqrt(2.0), + cry = joinDistance * 0.5 * sqrt(2.0), + angle = 90.0, + largeArcFlag = false, + sweepFlag = true, + end = offsetContour.position(0.0) + ) + SegmentJoin.MITER -> { + val ls = lastSegment ?: offsetContours.last().segments.last() + val fs = offsetContour.segments.first() + val i = intersection( + ls.end, + ls.end + ls.direction(1.0), + fs.start, + fs.start - fs.direction(0.0), + eps = 10E8 + ) + if (i !== Vector2.INFINITY) { + lineTo(i) + lineTo(fs.start) + } else { + lineTo(fs.start) + } + } + } + } + for (offsetSegment in offsetContour.segments) { + segment(offsetSegment) + } + + } + if (this@offset.closed) { + close() + } + } + + val postProc = false + + var final = candidateContour.removeLoops() + + if (postProc && !final.empty) { + val head = Segment( + segments[0].start + segments[0].normal(0.0) + .perpendicular(polarity) * 1000.0, segments[0].start + ).offset(distance).firstOrNull()?.copy(end = final.segments[0].start)?.contour + + val tail = Segment( + segments.last().end, + segments.last().end - segments.last().normal(1.0) + .perpendicular(polarity) * 1000.0 + ).offset(distance).firstOrNull()?.copy(start = final.segments.last().end)?.contour + + if (head != null) { + val headInts = intersections(final, head) + if (headInts.size == 1) { + final = final.sub(headInts[0].a.contourT, 1.0) + } + if (headInts.size > 1) { + val sInts = headInts.sortedByDescending { it.a.contourT } + final = final.sub(sInts[0].a.contourT, 1.0) + } + } +// final = head + final +// + if (tail != null) { + val tailInts = intersections(final, tail) + if (tailInts.size == 1) { + final = final.sub(0.0, tailInts[0].a.contourT) + } + if (tailInts.size > 1) { + val sInts = tailInts.sortedBy { it.a.contourT } + final = final.sub(0.0, sInts[0].a.contourT) + } + } + +// final = final + tail + + } + + return final +} diff --git a/orx-shapes/src/commonMain/kotlin/simplify/Chaikin.kt b/orx-shapes/src/commonMain/kotlin/simplify/Chaikin.kt new file mode 100644 index 00000000..8bc10061 --- /dev/null +++ b/orx-shapes/src/commonMain/kotlin/simplify/Chaikin.kt @@ -0,0 +1,93 @@ +package org.openrndr.extra.shapes.simplify + +import org.openrndr.math.Vector2 + + +/** + * Chaikin's corner cutting algorithm generates an approximating curve from a [polyline] + * + * [Interactive Demo](https://observablehq.com/@infowantstobeseen/chaikins-curves) + * + * The code has been tweaked for performance + * instead of brevity or being idiomatic. + * + * @param polyline a list of vectors describing the polyline + * @param iterations the number of times to approximate + * @param closed when the polyline is supposed to be a closed shape + * @param bias a value above 0.0 and below 0.5 controlling + * where new vertices are located. Lower values produce vertices near + * existing vertices. Values near 0.5 biases new vertices towards + * the mid-point between existing vertices. + */ +tailrec fun chaikinSmooth( + polyline: List, + iterations: Int = 1, + closed: Boolean = false, + bias: Double = 0.25 +): List { + if (iterations <= 0 || polyline.size < 2) { + return polyline + } + + val biasInv = 1 - bias + val result = ArrayList(polyline.size * 2) + + if (closed) { + + val sz = polyline.size + for (i in 0 until sz) { + val p0 = polyline[i] // `if` is here faster than `%` + val p1 = polyline[if (i + 1 == sz) 0 else i + 1] + + val (p0x, p0y) = p0 + val (p1x, p1y) = p1 + + result.apply { + add( + Vector2( + biasInv * p0x + bias * p1x, + biasInv * p0y + bias * p1y + ) + ) + add( + Vector2( + bias * p0x + biasInv * p1x, + bias * p0y + biasInv * p1y + ) + ) + } + } + + } else { + + // make sure it starts at point 0 + result.add(polyline[0].copy()) + val sz = polyline.size - 1 + for (i in 0 until sz) { + val p0 = polyline[i] + val p1 = polyline[i + 1] + + val (p0x, p0y) = p0 + val (p1x, p1y) = p1 + + result.apply { + add( + Vector2( + biasInv * p0x + bias * p1x, + biasInv * p0y + bias * p1y + ) + ) + add( + Vector2( + bias * p0x + biasInv * p1x, + bias * p0y + biasInv * p1y + ) + ) + } + } + // make sure it ends at the last point + result.add(polyline[sz].copy()) + + } + return chaikinSmooth(result, iterations - 1, closed, bias) +} diff --git a/orx-shapes/src/commonMain/kotlin/simplify/RamerDouglasPeucker.kt b/orx-shapes/src/commonMain/kotlin/simplify/RamerDouglasPeucker.kt new file mode 100644 index 00000000..1ee80fcd --- /dev/null +++ b/orx-shapes/src/commonMain/kotlin/simplify/RamerDouglasPeucker.kt @@ -0,0 +1,45 @@ +package org.openrndr.extra.shapes.simplify + +import org.openrndr.shape.LineSegment +import org.openrndr.math.Vector2 + +/** + * The [Ramer–Douglas–Peucker algorithm](https://en.wikipedia.org/wiki/Ramer–Douglas–Peucker_algorithm), + * is an algorithm that decimates a curve composed of line segments to a similar curve with fewer points. + * + * When the epsilon is less than the distance between two points, [simplify] is applied recursively. + * + * @author Edwin Jakobs + * + * @param epsilon Maximum distance two points can have without simplifying them. + */ +fun simplify(points: List, epsilon: Double): List { + // Find the point with the maximum distance + + val startEndDistance = points.first().squaredDistanceTo(points.last()) + + val endIndex = if (startEndDistance < 1E-6) points.size - 2 else points.size - 1 + + var dMax = 0.0 + var index = 0 + val end = points.size + for (i in 1..(end - 2)) { + val ls = LineSegment(points[0], points[endIndex]).extend(1000000.0) + val d = ls.distance(points[i]) + if (d > dMax) { + index = i + dMax = d + } + + } + // If max distance is greater than epsilon, recursively org.openrndr.shape.simplify + return if (dMax > epsilon) { + // Recursive call + val recResults1 = simplify(points.subList(0, index + 1), epsilon) + val recResults2 = simplify(points.subList(index, end), epsilon) + // Build the result list + listOf(recResults1.subList(0, recResults1.lastIndex), recResults2).flatMap { it.toList() } + } else { + listOf(points[0], points[end - 1]) + } +} \ No newline at end of file diff --git a/orx-shapes/src/commonMain/kotlin/splines/CatmullRom.kt b/orx-shapes/src/commonMain/kotlin/splines/CatmullRom.kt new file mode 100644 index 00000000..59c7349e --- /dev/null +++ b/orx-shapes/src/commonMain/kotlin/splines/CatmullRom.kt @@ -0,0 +1,333 @@ +package org.openrndr.extra.shapes.splines + +import org.openrndr.math.Vector2 +import org.openrndr.math.Vector3 +import org.openrndr.math.mod +import org.openrndr.shape.* +import kotlin.jvm.JvmOverloads +import kotlin.math.abs +import kotlin.math.pow + +private const val almostZero = 0.00000001 +private const val almostOne = 0.99999999 + +/** + * Creates a 1D Catmull-Rom spline curve. + * + * @param p0 The first control point. + * @param p1 The starting anchor point. + * @param p2 The ending anchor point. + * @param p3 The second control point. + * @param alpha The *tension* of the curve. + * Use `0.0` for the uniform spline, `0.5` for the centripetal spline, `1.0` for the chordal spline. + */ +class CatmullRom1 @JvmOverloads constructor (val p0: Double, val p1: Double, val p2: Double, val p3: Double, val alpha: Double = 0.5) { + /** Value of t for p0. */ + val t0: Double = 0.0 + /** Value of t for p1. */ + val t1: Double = calculateT(t0, p0, p1) + /** Value of t for p2. */ + val t2: Double = calculateT(t1, p1, p2) + /** Value of t for p3. */ + val t3: Double = calculateT(t2, p2, p3) + + private fun f(x: Double): Double = if (abs(x) < almostZero) 1.0 else x + + /** + * @param rt segment parameter value in [0, 1] + * @return a position on the segment + */ + fun position(rt: Double): Double { + val t = (t2 - t1) * rt + t1 + + val a1 = p0 * ((t1 - t) / f(t1 - t0)) + p1 * ((t - t0) / f(t1 - t0)) + val a2 = p1 * ((t2 - t) / f(t2 - t1)) + p2 * ((t - t1) / f(t2 - t1)) + val a3 = p2 * ((t3 - t) / f(t3 - t2)) + p3 * ((t - t2) / f(t3 - t2)) + + val b1 = a1 * ((t2 - t) / f(t2 - t0)) + a2 * ((t - t0) / f(t2 - t0)) + val b2 = a2 * ((t3 - t) / f(t3 - t1)) + a3 * ((t - t1) / f(t3 - t1)) + + val c = b1 * ((t2 - t) / f(t2 - t1)) + b2 * ((t - t1) / f(t2 - t1)) + return c + } + + private fun calculateT(t: Double, p0: Double, p1: Double): Double { + val a = (p1 - p0).pow(2.0) + val b = a.pow(0.5) + val c = b.pow(alpha) + return c + t + } +} + +/** + * Calculates the 1D Catmull–Rom spline for a chain of points and returns the combined curve. + * + * For more details, see [CatmullRom1]. + * + * @param points The [List] of 1D points where [CatmullRom1] is applied in groups of 4. + * @param alpha The *tension* of the curve. + * Use `0.0` for the uniform spline, `0.5` for the centripetal spline, `1.0` for the chordal spline. + * @param loop Whether to connect the first and last point, such that it forms a closed shape. + */ +class CatmullRomChain1 @JvmOverloads constructor (points: List, alpha: Double = 0.5, val loop: Boolean = false) { + val segments = if (!loop) points.windowed(4, 1).map { + CatmullRom1(it[0], it[1], it[2], it[3], alpha) + } else { + val cleanPoints = if (loop && abs(points.first() - (points.last())) <= 1.0E-6) { + points.dropLast(1) + } else { + points + } + (cleanPoints + cleanPoints.take(3)).windowed(4, 1).map { + CatmullRom1(it[0], it[1], it[2], it[3], alpha) + } + } + + fun position(rt: Double): Double { + val st = if (loop) mod(rt, 1.0) else rt.coerceIn(0.0, 1.0) + val segmentIndex = (kotlin.math.min(almostOne, st) * segments.size).toInt() + val t = (kotlin.math.min(almostOne, st) * segments.size) - segmentIndex + return segments[segmentIndex].position(t) + } +} + +/** + * Creates a 2D Catmull-Rom spline curve. + * + * Can be represented as a segment drawn between [p1] and [p2], + * while [p0] and [p3] are used as control points. + * + * Under some circumstances alpha can have + * no perceptible effect, for example, + * when creating closed shapes with the vertices + * forming a regular 2D polygon. + * + * @param p0 The first control point. + * @param p1 The starting anchor point. + * @param p2 The ending anchor point. + * @param p3 The second control point. + * @param alpha The *tension* of the curve. + * Use `0.0` for the uniform spline, `0.5` for the centripetal spline, `1.0` for the chordal spline. + */ +class CatmullRom2 @JvmOverloads constructor (val p0: Vector2, val p1: Vector2, val p2: Vector2, val p3: Vector2, val alpha: Double = 0.5) { + /** Value of t for p0. */ + val t0: Double = 0.0 + /** Value of t for p1. */ + val t1: Double = calculateT(t0, p0, p1) + /** Value of t for p2. */ + val t2: Double = calculateT(t1, p1, p2) + /** Value of t for p3. */ + val t3: Double = calculateT(t2, p2, p3) + + fun position(rt: Double): Vector2 { + val t = t1 + rt * (t2 - t1) + val a1 = p0 * ((t1 - t) / (t1 - t0)) + p1 * ((t - t0) / (t1 - t0)) + val a2 = p1 * ((t2 - t) / (t2 - t1)) + p2 * ((t - t1) / (t2 - t1)) + val a3 = p2 * ((t3 - t) / (t3 - t2)) + p3 * ((t - t2) / (t3 - t2)) + + val b1 = a1 * ((t2 - t) / (t2 - t0)) + a2 * ((t - t0) / (t2 - t0)) + val b2 = a2 * ((t3 - t) / (t3 - t1)) + a3 * ((t - t1) / (t3 - t1)) + + val c = b1 * ((t2 - t) / (t2 - t1)) + b2 * ((t - t1) / (t2 - t1)) + return c + } + + private fun calculateT(t: Double, p0: Vector2, p1: Vector2): Double { + val a = (p1.x - p0.x).pow(2.0) + (p1.y - p0.y).pow(2.0) + val b = a.pow(0.5) + val c = b.pow(alpha) + return c + t + } +} + +/** + * Calculates the 2D Catmull–Rom spline for a chain of points and returns the combined curve. + * + * For more details, see [CatmullRom2]. + * + * @param points The [List] of 2D points where [CatmullRom2] is applied in groups of 4. + * @param alpha The *tension* of the curve. + * Use `0.0` for the uniform spline, `0.5` for the centripetal spline, `1.0` for the chordal spline. + * @param loop Whether to connect the first and last point, such that it forms a closed shape. + */ +class CatmullRomChain2 @JvmOverloads constructor (points: List, alpha: Double = 0.5, val loop: Boolean = false) { + val segments = if (!loop) { + val startPoints = points.take(2) + val endPoints = points.takeLast(2) + val mirrorStart = + startPoints.first() - (startPoints.last() - startPoints.first()).normalized + val mirrorEnd = endPoints.last() + (endPoints.last() - endPoints.first()).normalized + + (listOf(mirrorStart) + points + listOf(mirrorEnd)).windowed(4, 1).map { + CatmullRom2(it[0], it[1], it[2], it[3], alpha) + } + } else { + val cleanPoints = if (loop && points.first().distanceTo(points.last()) <= 1.0E-6) { + points.dropLast(1) + } else { + points + } + (cleanPoints + cleanPoints.take(3)).windowed(4, 1).map { + CatmullRom2(it[0], it[1], it[2], it[3], alpha) + } + } + + fun positions(steps: Int = segments.size * 4): List { + return (0..steps).map { + position(it.toDouble() / steps) + } + } + + fun position(rt: Double): Vector2 { + val st = if (loop) mod(rt, 1.0) else rt.coerceIn(0.0, 1.0) + val segmentIndex = (kotlin.math.min(almostOne, st) * segments.size).toInt() + val t = (kotlin.math.min(almostOne, st) * segments.size) - segmentIndex + return segments[segmentIndex].position(t) + } +} + +/** + * Creates a 3D Catmull-Rom spline curve. + * + * Can be represented as a segment drawn between [p1] and [p2], + * while [p0] and [p3] are used as control points. + * + * Under some circumstances alpha can have + * no perceptible effect, for example, + * when creating closed shapes with the vertices + * forming a regular 2D polygon (even on a 3D plane). + * + * @param p0 The first control point. + * @param p1 The starting anchor point. + * @param p2 The ending anchor point. + * @param p3 The second control point. + * @param alpha The *tension* of the curve. + * Use `0.0` for the uniform spline, `0.5` for the centripetal spline, `1.0` for the chordal spline. + */ +class CatmullRom3 @JvmOverloads constructor (val p0: Vector3, val p1: Vector3, val p2: Vector3, val p3: Vector3, val alpha: Double = 0.5) { + /** Value of t for p0. */ + val t0: Double = 0.0 + /** Value of t for p1. */ + val t1: Double = calculateT(t0, p0, p1) + /** Value of t for p2. */ + val t2: Double = calculateT(t1, p1, p2) + /** Value of t for p3. */ + val t3: Double = calculateT(t2, p2, p3) + + fun position(rt: Double): Vector3 { + val t = t1 + rt * (t2 - t1) + val a1 = p0 * ((t1 - t) / (t1 - t0)) + p1 * ((t - t0) / (t1 - t0)) + val a2 = p1 * ((t2 - t) / (t2 - t1)) + p2 * ((t - t1) / (t2 - t1)) + val a3 = p2 * ((t3 - t) / (t3 - t2)) + p3 * ((t - t2) / (t3 - t2)) + + val b1 = a1 * ((t2 - t) / (t2 - t0)) + a2 * ((t - t0) / (t2 - t0)) + val b2 = a2 * ((t3 - t) / (t3 - t1)) + a3 * ((t - t1) / (t3 - t1)) + + val c = b1 * ((t2 - t) / (t2 - t1)) + b2 * ((t - t1) / (t2 - t1)) + return c + } + + private fun calculateT(t: Double, p0: Vector3, p1: Vector3): Double { + val a = (p1.x - p0.x).pow(2.0) + (p1.y - p0.y).pow(2.0) + (p1.z - p0.z).pow(2.0) + val b = a.pow(0.5) + val c = b.pow(alpha) + return c + t + } +} + +/** + * Calculates the 3D Catmull–Rom spline for a chain of points and returns the combined curve. + * + * For more details, see [CatmullRom3]. + * + * @param points The [List] of 3D points where [CatmullRom3] is applied in groups of 4. + * @param alpha The *tension* of the curve. + * Use `0.0` for the uniform spline, `0.5` for the centripetal spline, `1.0` for the chordal spline. + * @param loop Whether to connect the first and last point, such that it forms a closed shape. + */ +class CatmullRomChain3 @JvmOverloads constructor (points: List, alpha: Double = 0.5, val loop: Boolean = false) { + val segments = if (!loop) { + val startPoints = points.take(2) + val endPoints = points.takeLast(2) + val mirrorStart = + startPoints.first() - (startPoints.last() - startPoints.first()).normalized + val mirrorEnd = endPoints.last() + (endPoints.last() - endPoints.first()).normalized + + (listOf(mirrorStart) + points + listOf(mirrorEnd)).windowed(4, 1).map { + CatmullRom3(it[0], it[1], it[2], it[3], alpha) + } + } else { + val cleanPoints = if (loop && points.first().distanceTo(points.last()) <= 1.0E-6) { + points.dropLast(1) + } else { + points + } + (cleanPoints + cleanPoints + cleanPoints.take(3)).windowed(4, 1).map { + CatmullRom3(it[0], it[1], it[2], it[3], alpha) + } + } + + fun positions(steps: Int = segments.size * 4): List { + return (0..steps).map { + position(it.toDouble() / steps) + } + } + + fun position(rt: Double): Vector3 { + val st = if (loop) mod(rt, 1.0) else rt.coerceIn(0.0, 1.0) + val segmentIndex = (kotlin.math.min(almostOne, st) * segments.size).toInt() + val t = (kotlin.math.min(almostOne, st) * segments.size) - segmentIndex + return segments[segmentIndex].position(t) + } +} + +fun List.catmullRom(alpha: Double = 0.5, closed: Boolean) = CatmullRomChain2(this, alpha, closed) + +fun List.catmullRom(alpha: Double = 0.5, closed: Boolean) = CatmullRomChain3(this, alpha, closed) + + +/** Converts spline to a [Segment]. */ +fun CatmullRom2.toSegment(): Segment { + val d1a2 = (p1 - p0).length.pow(2 * alpha) + val d2a2 = (p2 - p1).length.pow(2 * alpha) + val d3a2 = (p3 - p2).length.pow(2 * alpha) + val d1a = (p1 - p0).length.pow(alpha) + val d2a = (p2 - p1).length.pow(alpha) + val d3a = (p3 - p2).length.pow(alpha) + + val b0 = p1 + val b1 = (p2 * d1a2 - p0 * d2a2 + p1 * (2 * d1a2 + 3 * d1a * d2a + d2a2)) / (3 * d1a * (d1a + d2a)) + val b2 = (p1 * d3a2 - p3 * d2a2 + p2 * (2 * d3a2 + 3 * d3a * d2a + d2a2)) / (3 * d3a * (d3a + d2a)) + val b3 = p2 + + return Segment(b0, b1, b2, b3) +} + + + +/** + * Converts chain to a [ShapeContour]. + */ +@Suppress("unused") +fun CatmullRomChain2.toContour(): ShapeContour = + ShapeContour(segments.map { it.toSegment() }, this.loop) + + + +fun CatmullRom3.toSegment(): Segment3D { + val d1a2 = (p1 - p0).length.pow(2 * alpha) + val d2a2 = (p2 - p1).length.pow(2 * alpha) + val d3a2 = (p3 - p2).length.pow(2 * alpha) + val d1a = (p1 - p0).length.pow(alpha) + val d2a = (p2 - p1).length.pow(alpha) + val d3a = (p3 - p2).length.pow(alpha) + + val b0 = p1 + val b1 = (p2 * d1a2 - p0 * d2a2 + p1 * (2 * d1a2 + 3 * d1a * d2a + d2a2)) / (3 * d1a * (d1a + d2a)) + val b2 = (p1 * d3a2 - p3 * d2a2 + p2 * (2 * d3a2 + 3 * d3a * d2a + d2a2)) / (3 * d3a * (d3a + d2a)) + val b3 = p2 + + return Segment3D(b0, b1, b2, b3) +} + +fun CatmullRomChain3.toPath3D(): Path3D = Path3D(segments.map { it.toSegment() }, this.loop) \ No newline at end of file