[orx-shapes] Add 3d Hobby curves

This commit is contained in:
Edwin Jakobs
2024-04-03 08:15:14 +02:00
parent c9f9f67dec
commit 70ffc7dc14
3 changed files with 249 additions and 70 deletions

View File

@@ -21,6 +21,12 @@ kotlin {
}
val commonTest by getting {
dependencies {
implementation(project(":orx-noise"))
}
}
val jvmTest by getting {
dependencies {
implementation(libs.kotest.assertions)

View File

@@ -1,15 +1,10 @@
package org.openrndr.extra.shapes.hobbycurve
// Code adapted from http://weitz.de/hobby/
import org.openrndr.math.Vector2
import org.openrndr.shape.Segment2D
import org.openrndr.shape.Shape
import org.openrndr.shape.ShapeContour
import kotlin.math.atan2
import kotlin.math.cos
import kotlin.math.sin
import kotlin.math.sqrt
import org.openrndr.math.*
import org.openrndr.math.transforms.buildTransform
import org.openrndr.shape.*
import kotlin.math.*
fun ShapeContour.hobbyCurve(curl: Double = 0.0): ShapeContour {
val vertices = if (closed)
@@ -19,12 +14,29 @@ fun ShapeContour.hobbyCurve(curl: Double = 0.0): ShapeContour {
return hobbyCurve(vertices, closed, curl)
}
fun Shape.hobbyCurve(curl: Double = 0.0) : Shape {
fun Shape.hobbyCurve(curl: Double = 0.0): Shape {
return Shape(contours.map {
it.hobbyCurve(curl)
})
}
private fun Vector2.atan22(other: Vector2): Double {
val u = this.normalized
val v = other.normalized
val x = u.cross(v)
val y = u.dot(v)
return atan2(x, y)
}
private fun Vector3.atan22(other: Vector3): Double {
val u = this.normalized
val v = other.normalized
val x = u.cross(v).length
val y = u.dot(v)
return atan2(x, y)
}
/**
* Uses Hobby's algorithm to construct a [ShapeContour] through a given list of points.
* @param points The list of points through which the curve should go.
@@ -38,109 +50,222 @@ fun hobbyCurve(points: List<Vector2>, closed: Boolean = false, curl: Double = 0.
val m = points.size
val n = if (closed) m else m - 1
val diffs = Array(n) { points[(it+1) % m] - points[it] }
val distances = Array(n) { diffs[it].length }
val chords = Array(n) { points[(it + 1) % m] - points[it] }
val distances = Array(n) { chords[it].length }
val gamma = arrayOfNulls<Double>(m)
for (i in (if (closed) 0 else 1) until n){
val k = (i + m - 1) % m
val n1 = diffs[k].normalized
val s = n1.y
val c = n1.x
val v = rotate(diffs[i], -s, c)
gamma[i] = atan2(v.y, v.x)
require(distances.all { it > 0.0 })
val gamma = DoubleArray(m)
for (i in (if (closed) 0 else 1) until n) {
gamma[i] = chords[(i - 1).mod(m)].atan22(chords[(i).mod(m)])
}
if (!closed) gamma[n] = 0.0
val a = arrayOfNulls<Double>(m)
val b = arrayOfNulls<Double>(m)
val c = arrayOfNulls<Double>(m)
val d = arrayOfNulls<Double>(m)
val a = DoubleArray(n) { 0.0 }
val b = DoubleArray(n) { 0.0 }
val c = DoubleArray(n) { 0.0 }
val d = DoubleArray(n) { 0.0 }
for (i in (if (closed) 0 else 1) until n){
val j = (i + 1) % m
val k = (i + m - 1) % m
for (i in (if (closed) 0 else 1) until n) {
val j = (i + 1).mod(m)
val k = (i - 1).mod(m)
a[i] = 1 / distances[k]
b[i] = (2 * distances[k] + 2 * distances[i]) / (distances[k] * distances[i])
c[i] = 1 / distances[i]
d[i] = -(2 * gamma[i]!! * distances[i] + gamma[j]!! * distances[k]) / (distances[k] * distances[i])
d[i] = -(2 * gamma[i] * distances[i] + gamma[j] * distances[k]) / (distances[k] * distances[i])
}
lateinit var alpha: Array<Double>
lateinit var beta: Array<Double?>
val alpha: DoubleArray
val beta: DoubleArray
if (!closed) {
a[0] = 0.0
b[0] = 2 + curl
c[0] = 2 * curl + 1
d[0] = -c[0]!! * gamma[1]!!
d[0] = -c[0] * gamma[1]
a[n] = 2 * curl + 1
b[n] = 2 + curl
c[n] = 0.0
d[n] = 0.0
alpha = thomas(a.requireNoNulls(), b.requireNoNulls(), c.requireNoNulls(), d.requireNoNulls())
beta = arrayOfNulls(n)
for (i in 0 until n-1){
beta[i] = -gamma[i+1]!! - alpha[i+1]
alpha = thomas(a, b, c, d)
beta = DoubleArray(n) { 0.0 }
for (i in 0 until n - 1) {
beta[i] = -gamma[i + 1] - alpha[i + 1]
}
beta[n-1] = -alpha[n]
beta[n - 1] = -alpha[n]
} else {
val s = a[0]!!
val s = a[0]
a[0] = 0.0
val t = c[n-1]!!
c[n-1] = 0.0
alpha = sherman(a.requireNoNulls(), b.requireNoNulls(), c.requireNoNulls(), d.requireNoNulls(), s, t)
beta = arrayOfNulls(n)
for (i in 0 until n){
val j = (i+1) % n
beta[i] = -gamma[j]!! - alpha[j]
val t = c[n - 1]
c[n - 1] = 0.0
alpha = sherman(a, b, c, d, s, t)
beta = DoubleArray(n) { 0.0 }
for (i in 0 until n) {
val j = (i + 1) % n
beta[i] = -gamma[j] - alpha[j]
}
}
val c1s = mutableListOf<Vector2>()
val c2s = mutableListOf<Vector2>()
for (i in 0 until n){
val v1 = rotateAngle(diffs[i], alpha[i]).normalized
val v2 = rotateAngle(diffs[i], -beta[i]!!).normalized
c1s.add(points[i % m] + v1 * rho(alpha[i], beta[i]!!) * distances[i] / 3.0)
c2s.add(points[(i+1) % m] - v2 * rho(beta[i]!!, alpha[i]) * distances[i] / 3.0)
for (i in 0 until n) {
val v1 = rotateAngle(chords[i], alpha[i]).normalized
val v2 = rotateAngle(chords[i], -beta[i]).normalized
c1s.add(points[i % m] + v1 * rho(alpha[i], beta[i]) * distances[i] / 3.0)
c2s.add(points[(i + 1) % m] - v2 * rho(beta[i], alpha[i]) * distances[i] / 3.0)
}
return ShapeContour(List(n) { Segment2D(points[it], c1s[it], c2s[it], points[(it+1)%m]) }, closed=closed)
return ShapeContour(List(n) { Segment2D(points[it], c1s[it], c2s[it], points[(it + 1) % m]) }, closed = closed)
}
private fun thomas(a: Array<Double>, b: Array<Double>, c: Array<Double>, d: Array<Double>): Array<Double> {
/**
* Uses Hobby's algorithm to construct a [Path3D] through a given list of points.
* @param points The list of points through which the curve should go.
* @param closed Whether to construct a closed or open curve.
* @param curl The 'curl' at the endpoints of the curve; this is only applicable when [closed] is false. Best results for values in [-1, 1], where a higher value makes segments closer to circular arcs.
* @param tensions A function that returns the in and out tensions given a chord index.
* @return A [Path3D] through [points].
*/
fun hobbyCurve(
points: List<Vector3>,
closed: Boolean = false,
curl: Double = 0.0,
tensions: (chordIndex: Int) -> Pair<Double, Double> = { _ -> Pair(1.0, 1.0) }
): Path3D {
if (points.size <= 1) return Path3D.EMPTY
val m = points.size
val n = if (closed) m else m - 1
val chords = Array(n) { points[(it + 1) % m] - points[it] }
val distances = Array(n) { chords[it].length }
val normals = Array(n) { Vector3.ZERO }
require(distances.all { it > 0.0 })
val gamma = DoubleArray(m)
for (i in (if (closed) 0 else 1) until n) {
val zc = chords[i].normalized
val zp = chords[(i - 1).mod(m)].normalized
val normal = zc.cross(-zp)
gamma[i] = zp.atan22(zc) * sign(normal.z)
normals[i] = normal * sign(normal.z)
}
if (!closed) {
gamma[n] = 0.0
}
val a = DoubleArray(m) { 0.0 }
val b = DoubleArray(m) { 0.0 }
val c = DoubleArray(m) { 0.0 }
val d = DoubleArray(m) { 0.0 }
for (i in (if (closed) 0 else 1) until n) {
val j = (i + 1).mod(m)
val k = (i - 1).mod(m)
a[i] = 1 / distances[k]
b[i] = (2 * distances[k] + 2 * distances[i]) / (distances[k] * distances[i])
c[i] = 1 / distances[i]
d[i] = -(2 * gamma[i] * distances[i] + gamma[j] * distances[k]) / (distances[k] * distances[i])
}
val alpha: DoubleArray
val beta: DoubleArray
if (!closed) {
a[0] = 0.0
b[0] = 2 + curl
c[0] = 2 * curl + 1
d[0] = -c[0] * gamma[1]
a[n] = 2 * curl + 1
b[n] = 2 + curl
c[n] = 0.0
d[n] = 0.0
alpha = thomas(a, b, c, d)
beta = DoubleArray(n) { 0.0 }
for (i in 0 until n - 1) {
beta[i] = -gamma[i + 1] - alpha[i + 1]
}
beta[n - 1] = -alpha[n]
} else {
val s = a[0]
a[0] = 0.0
val t = c[n - 1]
c[n - 1] = 0.0
alpha = sherman(a, b, c, d, s, t)
beta = DoubleArray(n) { 0.0 }
for (i in 0 until n) {
val j = (i + 1) % n
beta[i] = -gamma[j] - alpha[j]
}
}
val c1s = mutableListOf<Vector3>()
val c2s = mutableListOf<Vector3>()
for (i in 0 until n) {
val r1 = buildTransform { rotate(normals[i], alpha[i].asDegrees) }
val r2 = buildTransform { rotate(normals[(i + 1).mod(normals.size)], -beta[i].asDegrees) }
val v1 = (r1 * chords[i].xyz0).xyz.normalized
val v2 = (r2 * chords[i].xyz0).xyz.normalized
val t = tensions(i)
c1s.add(points[i % m] + v1 * rho(alpha[i], beta[i]) * distances[i] * t.first / 3.0)
c2s.add(points[(i + 1) % m] - v2 * rho(beta[i], alpha[i]) * distances[i] * t.second / 3.0)
}
return Path3D(List(n) {
Segment3D(points[it], c1s[it], c2s[it], points[(it + 1) % m])
}, closed = closed)
}
/** The Thomas algorithm: solve a system of linear equations encoded in a tridiagonal matrix.
https://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm
*/
private fun thomas(a: DoubleArray, b: DoubleArray, c: DoubleArray, d: DoubleArray): DoubleArray {
val n = a.size
val cc = arrayOfNulls<Double>(n)
val dd = arrayOfNulls<Double>(n)
val cc = DoubleArray(n) { 0.0 }
val dd = DoubleArray(n) { 0.0 }
cc[0] = c[0] / b[0]
dd[0] = d[0] / b[0]
for (i in 1 until n){
val den = b[i] - cc[i-1]!! * a[i]
for (i in 1 until n) {
val den = b[i] - cc[i - 1] * a[i]
cc[i] = c[i] / den
dd[i] = (d[i] - dd[i-1]!!*a[i]) / den
dd[i] = (d[i] - dd[i - 1] * a[i]) / den
}
val x = arrayOfNulls<Double>(n)
x[n-1] = dd[n-1]
for (i in n-2 downTo 0){
x[i] = dd[i]!! - cc[i]!! * x[i+1]!!
val x = DoubleArray(n) { 0.0 }
x[n - 1] = dd[n - 1]
for (i in n - 2 downTo 0) {
x[i] = dd[i] - cc[i] * x[i + 1]
}
return x.requireNoNulls()
return x
}
private fun sherman(a: Array<Double>, b: Array<Double>, c: Array<Double>, d: Array<Double>, s: Double, t: Double): Array<Double> {
private fun sherman(
a: DoubleArray,
b: DoubleArray,
c: DoubleArray,
d: DoubleArray,
s: Double,
t: Double
): DoubleArray {
val n = a.size
val u = Array(n) { if (it == 0 || it == n-1) 1.0 else 0.0 }
val v = Array(n) { when (it){ 0 -> t; n-1 -> s; else -> 0.0 } }
val u = DoubleArray(n) { if (it == 0 || it == n - 1) 1.0 else 0.0 }
val v = DoubleArray(n) {
when (it) {
0 -> t; n - 1 -> s; else -> 0.0
}
}
b[0] -= t
b[n-1] -= s
b[n - 1] -= s
val Td = thomas(a, b, c, d)
val Tu = thomas(a, b, c, u)
val factor = (t * Td[0] + s*Td[n-1]) / (1 + t * Tu[0] + s*Tu[n-1])
return Array(n) {
val factor = (t * Td[0] + s * Td[n - 1]) / (1 + t * Tu[0] + s * Tu[n - 1])
return DoubleArray(n) {
Td[it] - factor * Tu[it]
}
}
@@ -151,9 +276,9 @@ private fun rho(a: Double, b: Double): Double {
val ca = cos(a)
val cb = cos(b)
val s5 = sqrt(5.0)
val num = 4 + sqrt(8.0) * (sa - sb/16) * (sb - sa/16) * (ca - cb)
val num = 4 + sqrt(8.0) * (sa - sb / 16) * (sb - sa / 16) * (ca - cb)
val den = 2 + (s5 - 1) * ca + (3 - s5) * cb
return num/den
return num / den
}
private fun rotate(v: Vector2, s: Double, c: Double) = Vector2(v.x * c - v.y * s, v.x * s + v.y * c)

View File

@@ -0,0 +1,48 @@
package hobbycurve
import org.openrndr.WindowMultisample
import org.openrndr.application
import org.openrndr.color.ColorRGBa
import org.openrndr.draw.isolated
import org.openrndr.extra.noise.scatter
import org.openrndr.extra.noise.uniform
import org.openrndr.extra.shapes.hobbycurve.hobbyCurve
import org.openrndr.math.Vector3
import kotlin.math.cos
import kotlin.random.Random
fun main() = application {
configure {
width = 720
height = 720
multisample = WindowMultisample.SampleCount(4)
}
program {
val pts = drawer.bounds.scatter(30.0, distanceToEdge = 200.0, random = Random(3000))
extend {
drawer.stroke = ColorRGBa.PINK
drawer.strokeWeight = 4.0
drawer.fill = null
val r = Random(3000)
val hobby3D = hobbyCurve(
pts.map { it.xy0 + Vector3(0.0, 0.0, Double.uniform(-360.0, 360.0, r)) },
true,
tensions = { chordIndex: Int ->
Pair(
cos(seconds + chordIndex * 0.1) * 0.5 + 0.5,
cos(seconds + (1.0 + chordIndex) * 0.1) * 0.5 + 0.5
)
})
drawer.isolated {
drawer.ortho(0.0, width.toDouble(), height.toDouble(), 0.0, -4000.0, 4000.0)
drawer.translate(width / 2.0, height / 2.0)
drawer.rotate(Vector3.UNIT_Y, seconds * 16.0)
drawer.translate(-width / 2.0, -height / 2.0)
drawer.path(hobby3D)
}
}
}
}