diff --git a/orx-triangulation/src/commonMain/kotlin/Delaunator.kt b/orx-triangulation/src/commonMain/kotlin/Delaunator.kt index efa59df7..0958d1cc 100644 --- a/orx-triangulation/src/commonMain/kotlin/Delaunator.kt +++ b/orx-triangulation/src/commonMain/kotlin/Delaunator.kt @@ -2,7 +2,7 @@ package org.openrndr.extra.triangulation import kotlin.math.* -private val EPSILON: Double = 2.0.pow(-52) +val EPSILON: Double = 2.0.pow(-52) /** * A Kotlin port of Mapbox's Delaunator incredibly fast JavaScript library for Delaunay triangulation of 2D points. @@ -14,8 +14,8 @@ private val EPSILON: Double = 2.0.pow(-52) * @author Ricardo Matias */ @Suppress("unused") -internal class Delaunator(val coords: DoubleArray) { - private val EDGE_STACK = IntArray(512) +class Delaunator(val coords: DoubleArray) { + val EDGE_STACK = IntArray(512) private var count = coords.size shr 1 @@ -366,7 +366,7 @@ internal class Delaunator(val coords: DoubleArray) { val pl = _triangles[al] val p1 = _triangles[bl] - val illegal = inCircle( + val illegal = inCircleRobust( coords[2 * p0], coords[2 * p0 + 1], coords[2 * pr], coords[2 * pr + 1], coords[2 * pl], coords[2 * pl + 1], @@ -528,7 +528,7 @@ private fun pseudoAngle(dx: Double, dy: Double): Double { return a / 4.0 // [0..1] } -/* + private fun inCircle(ax: Double, ay: Double, bx: Double, by: Double, cx: Double, cy: Double, @@ -547,9 +547,9 @@ private fun inCircle(ax: Double, ay: Double, return dx * (ey * cp - bp * fy) - dy * (ex * cp - bp * fx) + ap * (ex * fy - ey * fx) < 0 -}*/ +} -private fun inCircle( +private fun inCircleRobust( ax: Double, ay: Double, bx: Double, by: Double, cx: Double, cy: Double, @@ -574,8 +574,7 @@ private fun inCircle( ), ddMultDd(ap, ddDiffDd(ddMultDd(ex, fy), ddMultDd(ey, fx))) ) - // add a small bias here, it seems to help - return (dd[1]) <= 1E-8 + return (dd[1]) <= 0 } diff --git a/orx-triangulation/src/commonMain/kotlin/Delaunay.kt b/orx-triangulation/src/commonMain/kotlin/Delaunay.kt index 726d2ac6..a312cfea 100644 --- a/orx-triangulation/src/commonMain/kotlin/Delaunay.kt +++ b/orx-triangulation/src/commonMain/kotlin/Delaunay.kt @@ -2,6 +2,10 @@ package org.openrndr.extra.triangulation import org.openrndr.math.Vector2 import org.openrndr.shape.Rectangle +import org.openrndr.shape.Triangle +import org.openrndr.shape.contour +import org.openrndr.shape.contours +import kotlin.js.JsName import kotlin.math.cos import kotlin.math.pow import kotlin.math.sin @@ -73,13 +77,13 @@ class Delaunay(val points: DoubleArray) { init() } - fun neighbors(i: Int) = sequence { - val e0 = inedges.getOrNull(i) ?: return@sequence + fun neighbors(i:Int) = sequence { + val e0 = inedges[i] if (e0 != -1) { var e = e0 - var p0 : Int + var p0 = -1 - loop@ do { + loop@do { p0 = triangles[e] yield(p0) e = if (e % 3 == 2) e - 2 else e + 1 @@ -105,31 +109,29 @@ class Delaunay(val points: DoubleArray) { } fun collinear(): Boolean { - for (i in triangles.indices step 3) { + for (i in 0 until triangles.size step 3) { val a = 2 * triangles[i] val b = 2 * triangles[i + 1] - val c = 2 * triangles[i + 2] + val c = 2 * triangles[i + 2] val coords = points val cross = (coords[c] - coords[a]) * (coords[b + 1] - coords[a + 1]) - -(coords[b] - coords[a]) * (coords[c + 1] - coords[a + 1]) - if (cross > 1e-10) return false + - (coords[b] - coords[a]) * (coords[c + 1] - coords[a + 1]) + if (cross > 1e-10) return false; } return true } - - private fun jitter(x: Double, y: Double, r: Double): DoubleArray { - return doubleArrayOf(x + sin(x + y) * r, y + cos(x - y) * r) + private fun jitter(x:Double, y:Double, r:Double): DoubleArray { + return doubleArrayOf(x + sin(x+y) * r, y + cos(x-y)*r) } - fun init() { if (hull.size > 2 && collinear()) { println("warning: triangulation is collinear") val r = 1E-8 - for (i in points.indices step 2) { - val p = jitter(points[i], points[i + 1], r) + for (i in 0 until points.size step 2) { + val p = jitter(points[i], points[i+1], r) points[i] = p[0] - points[i + 1] = p[1] + points[i+1] = p[1] } delaunator = Delaunator(points) @@ -139,6 +141,8 @@ class Delaunay(val points: DoubleArray) { } + + inedges.fill(-1) hullIndex.fill(-1) @@ -223,4 +227,3 @@ class Delaunay(val points: DoubleArray) { fun voronoi(bounds: Rectangle): Voronoi = Voronoi(this, bounds) } - diff --git a/orx-triangulation/src/commonMain/kotlin/DoubleDouble.kt b/orx-triangulation/src/commonMain/kotlin/DoubleDouble.kt index 58defed3..477adf04 100644 --- a/orx-triangulation/src/commonMain/kotlin/DoubleDouble.kt +++ b/orx-triangulation/src/commonMain/kotlin/DoubleDouble.kt @@ -16,9 +16,9 @@ import kotlin.math.pow * * See https://people.eecs.berkeley.edu/~jrs/papers/robustr.pdf */ -internal fun fastTwoDiff(a: Double, b: Double): DoubleArray { - val x = a - b - val y = (a - x) - b +fun fastTwoDiff(a: Double, b: Double): DoubleArray { + val x = a - b; + val y = (a - x) - b; return doubleArrayOf(y, x) } @@ -33,7 +33,7 @@ internal fun fastTwoDiff(a: Double, b: Double): DoubleArray { * * See https://people.eecs.berkeley.edu/~jrs/papers/robustr.pdf */ -internal fun fastTwoSum(a: Double, b: Double): DoubleArray { +fun fastTwoSum(a: Double, b: Double): DoubleArray { val x = a + b; return doubleArrayOf(b - (x - a), x) @@ -55,16 +55,16 @@ internal fun fastTwoSum(a: Double, b: Double): DoubleArray { * @param a a double * @param bits the number of significand bits to leave intact */ -internal fun reduceSignificand( +fun reduceSignificand( a: Double, bits: Int ): Double { - val s = 53 - bits - val f = 2.0.pow(s) + 1 + val s = 53 - bits; + val f = 2.0.pow(s) + 1; - val c = f * a - val r = c - (c - a) + val c = f * a; + val r = c - (c - a); return r; } @@ -74,7 +74,7 @@ internal fun reduceSignificand( * === 2^Math.ceil(p/2) + 1 where p is the # of significand bits in a double === 53. * @internal */ -private const val f = 134217729 // 2**27 + 1; +const val f = 134217729; // 2**27 + 1; /** @@ -89,10 +89,10 @@ private const val f = 134217729 // 2**27 + 1; * see e.g. [Shewchuk](https://people.eecs.berkeley.edu/~jrs/papers/robustr.pdf) * @param a A double floating point number */ -private fun split(a: Double): DoubleArray { - val c = f * a - val a_h = c - (c - a) - val a_l = a - a_h +fun split(a: Double): DoubleArray { + val c = f * a; + val a_h = c - (c - a); + val a_l = a - a_h; return doubleArrayOf(a_h, a_l) } @@ -103,10 +103,10 @@ private fun split(a: Double): DoubleArray { * @param a minuend - a double-double precision floating point number * @param b subtrahend - a double-double precision floating point number */ -internal fun twoDiff(a: Double, b: Double): DoubleArray { - val x = a - b - val bvirt = a - x - val y = (a - (x + bvirt)) + (bvirt - b) +fun twoDiff(a: Double, b: Double): DoubleArray { + val x = a - b; + val bvirt = a - x; + val y = (a - (x + bvirt)) + (bvirt - b); return doubleArrayOf(y, x) } @@ -126,19 +126,19 @@ internal fun twoDiff(a: Double, b: Double): DoubleArray { * @param a A double * @param b Another double */ -internal fun twoProduct(a: Double, b: Double): DoubleArray { +fun twoProduct(a: Double, b: Double): DoubleArray { val x = a * b; //const [ah, al] = split(a); - val c = f * a - val ah = c - (c - a) - val al = a - ah + val c = f * a; + val ah = c - (c - a); + val al = a - ah; //const [bh, bl] = split(b); - val d = f * b - val bh = d - (d - b) - val bl = b - bh + val d = f * b; + val bh = d - (d - b); + val bl = b - bh; - val y = (al * bl) - ((x - (ah * bh)) - (al * bh) - (ah * bl)) + val y = (al * bl) - ((x - (ah * bh)) - (al * bh) - (ah * bl)); //const err1 = x - (ah * bh); //const err2 = err1 - (al * bh); @@ -148,15 +148,15 @@ internal fun twoProduct(a: Double, b: Double): DoubleArray { return doubleArrayOf(y, x) } -internal fun twoSquare(a: Double): DoubleArray { - val x = a * a +fun twoSquare(a: Double): DoubleArray { + val x = a * a; //const [ah, al] = split(a); - val c = f * a - val ah = c - (c - a) - val al = a - ah + val c = f * a; + val ah = c - (c - a); + val al = a - ah; - val y = (al * al) - ((x - (ah * ah)) - 2 * (ah * al)) + val y = (al * al) - ((x - (ah * ah)) - 2 * (ah * al)); return doubleArrayOf(y, x) } @@ -173,9 +173,9 @@ internal fun twoSquare(a: Double): DoubleArray { * * See https://people.eecs.berkeley.edu/~jrs/papers/robustr.pdf */ -internal fun twoSum(a: Double, b: Double): DoubleArray { - val x = a + b - val bv = x - a +fun twoSum(a: Double, b: Double): DoubleArray { + val x = a + b; + val bv = x - a; return doubleArrayOf((a - (x - bv)) + (b - bv), x) } @@ -193,22 +193,28 @@ internal fun twoSum(a: Double, b: Double): DoubleArray { * @param x a double-double precision floating point number * @param y another double-double precision floating point number */ -internal fun ddDiffDd(x: DoubleArray, y: DoubleArray): DoubleArray { - val xl = x[0] - val xh = x[1] - val yl = y[0] - val yh = y[1] +fun ddDiffDd(x: DoubleArray, y: DoubleArray): DoubleArray { + val xl = x[0]; + val xh = x[1]; + val yl = y[0]; + val yh = y[1]; //const [sl,sh] = twoSum(xh,yh); - val sh = xh - yh; val _1 = sh - xh; val sl = (xh - (sh - _1)) + (-yh - _1) + val sh = xh - yh; + val _1 = sh - xh; + val sl = (xh - (sh - _1)) + (-yh - _1); //const [tl,th] = twoSum(xl,yl); - val th = xl - yl; val _2 = th - xl; val tl = (xl - (th - _2)) + (-yl - _2) - val c = sl + th + val th = xl - yl; + val _2 = th - xl; + val tl = (xl - (th - _2)) + (-yl - _2); + val c = sl + th; //const [vl,vh] = fastTwoSum(sh,c) - val vh = sh + c; val vl = c - (vh - sh) + val vh = sh + c; + val vl = c - (vh - sh); val w = tl + vl //const [zl,zh] = fastTwoSum(vh,w) - val zh = vh + w; val zl = w - (zh - vh) + val zh = vh + w; + val zl = w - (zh - vh); return doubleArrayOf(zl, zh) } @@ -225,23 +231,27 @@ internal fun ddDiffDd(x: DoubleArray, y: DoubleArray): DoubleArray { * @param x a double-double precision floating point number * @param y another double-double precision floating point number */ -internal fun ddMultDd(x: DoubleArray, y: DoubleArray): DoubleArray { +fun ddMultDd(x: DoubleArray, y: DoubleArray): DoubleArray { //const xl = x[0]; - val xh = x[1] + val xh = x[1]; //const yl = y[0]; - val yh = y[1] + val yh = y[1]; //const [cl1,ch] = twoProduct(xh,yh); - val ch = xh*yh - val c = f * xh; val ah = c - (c - xh); val al = xh - ah - val d = f * yh; val bh = d - (d - yh); val bl = yh - bh - val cl1 = (al*bl) - ((ch - (ah*bh)) - (al*bh) - (ah*bl)) + val ch = xh * yh; + val c = f * xh; + val ah = c - (c - xh); + val al = xh - ah; + val d = f * yh; + val bh = d - (d - yh); + val bl = yh - bh; + val cl1 = (al * bl) - ((ch - (ah * bh)) - (al * bh) - (ah * bl)); //return fastTwoSum(ch,cl1 + (xh*yl + xl*yh)); - val b = cl1 + (xh*y[0] + x[0]*yh) - val xx = ch + b + val b = cl1 + (xh * y[0] + x[0] * yh); + val xx = ch + b; return doubleArrayOf(b - (xx - ch), xx) } @@ -260,22 +270,28 @@ internal fun ddMultDd(x: DoubleArray, y: DoubleArray): DoubleArray { * @param x a double-double precision floating point number * @param y another double-double precision floating point number */ -internal fun ddAddDd(x: DoubleArray, y: DoubleArray): DoubleArray { - val xl = x[0] - val xh = x[1] - val yl = y[0] - val yh = y[1] +fun ddAddDd(x: DoubleArray, y: DoubleArray): DoubleArray { + val xl = x[0]; + val xh = x[1]; + val yl = y[0]; + val yh = y[1]; //const [sl,sh] = twoSum(xh,yh); - val sh = xh + yh; val _1 = sh - xh; val sl = (xh - (sh - _1)) + (yh - _1) + val sh = xh + yh; + val _1 = sh - xh; + val sl = (xh - (sh - _1)) + (yh - _1); //val [tl,th] = twoSum(xl,yl); - val th = xl + yl; val _2 = th - xl; val tl = (xl - (th - _2)) + (yl - _2) - val c = sl + th + val th = xl + yl; + val _2 = th - xl; + val tl = (xl - (th - _2)) + (yl - _2); + val c = sl + th; //val [vl,vh] = fastTwoSum(sh,c) - val vh = sh + c; val vl = c - (vh - sh) + val vh = sh + c; + val vl = c - (vh - sh); val w = tl + vl //val [zl,zh] = fastTwoSum(vh,w) - val zh = vh + w; val zl = w - (zh - vh) + val zh = vh + w; + val zl = w - (zh - vh); return doubleArrayOf(zl, zh) } @@ -296,25 +312,29 @@ internal fun ddAddDd(x: DoubleArray, y: DoubleArray): DoubleArray { * @param y a double * @param x a double-double precision floating point number */ -internal fun ddMultDouble1(y: Double, x: DoubleArray): DoubleArray { - val xl = x[0] - val xh = x[1] +fun ddMultDouble1(y: Double, x: DoubleArray): DoubleArray { + val xl = x[0]; + val xh = x[1]; //val [cl1,ch] = twoProduct(xh,y); - val ch = xh*y - val c = f * xh; val ah = c - (c - xh); val al = xh - ah - val d = f * y; val bh = d - (d - y); val bl = y - bh - val cl1 = (al*bl) - ((ch - (ah*bh)) - (al*bh) - (ah*bl)) + val ch = xh * y; + val c = f * xh; + val ah = c - (c - xh); + val al = xh - ah; + val d = f * y; + val bh = d - (d - y); + val bl = y - bh; + val cl1 = (al * bl) - ((ch - (ah * bh)) - (al * bh) - (ah * bl)); - val cl2 = xl*y + val cl2 = xl * y; //val [tl1,th] = fastTwoSum(ch,cl2); - val th = ch + cl2 - val tl1 = cl2 - (th - ch) + val th = ch + cl2; + val tl1 = cl2 - (th - ch); - val tl2 = tl1 + cl1 + val tl2 = tl1 + cl1; //val [zl,zh] = fastTwoSum(th,tl2); - val zh = th + tl2 - val zl = tl2 - (zh - th) + val zh = th + tl2; + val zl = tl2 - (zh - th); - return doubleArrayOf(zl,zh) + return doubleArrayOf(zl, zh); } \ No newline at end of file diff --git a/orx-triangulation/src/commonMain/kotlin/Voronoi.kt b/orx-triangulation/src/commonMain/kotlin/Voronoi.kt index 7f3a873d..3073c34d 100644 --- a/orx-triangulation/src/commonMain/kotlin/Voronoi.kt +++ b/orx-triangulation/src/commonMain/kotlin/Voronoi.kt @@ -154,7 +154,11 @@ class Voronoi(val delaunay: Delaunay, val bounds: Rectangle) { } + private fun cell(i: Int): MutableList? { + + + val inedges = delaunay.inedges val halfedges = delaunay.halfedges val triangles = delaunay.triangles @@ -191,8 +195,8 @@ class Voronoi(val delaunay: Delaunay, val bounds: Rectangle) { if (cj != null) { val li = ci.size val lj = cj.size - loop@ for (ai in ci.indices step 2) { - for (aj in cj.indices step 2) { + loop@ for (ai in 0 until ci.size step 2) { + for (aj in 0 until cj.size step 2) { if (ci[ai] == cj[aj] && ci[ai + 1] == cj[aj + 1] && ci[(ai + 2) % li] == cj[(aj + lj - 2) % lj] @@ -252,10 +256,10 @@ class Voronoi(val delaunay: Delaunay, val bounds: Rectangle) { project(P[0], P[1], vx0, vy0)?.let { p -> P!!.add(0, p[1]); P!!.add(0, p[0]) } project(P[P.size - 2], P[P.size - 1], vxn, vyn)?.let { p -> P!!.add(p[0]); P!!.add(p[1]) } - P = this.clipFinite(i, P) + P = this.clipFinite(i, P!!) var n = 0 if (P != null) { - n = P.size + n = P!!.size var c0 = -1 var c1 = edgeCode(P[n - 2], P[n - 1]) var j = 0 @@ -294,7 +298,7 @@ class Voronoi(val delaunay: Delaunay, val bounds: Rectangle) { var y1 = points[n - 1] var c0: Int var c1: Int = regionCode(x1, y1) - var e0: Int? + var e0: Int? = null var e1: Int? = 0 for (j in 0 until n step 2) { @@ -356,7 +360,7 @@ class Voronoi(val delaunay: Delaunay, val bounds: Rectangle) { e0 = e1 e1 = this.edgeCode(P[0], P[1]) - if (e0.isTruthy() && e1.isTruthy()) this.edge(i, e0!!, e1, P, P.size); + if (e0.isTruthy() && e1.isTruthy()) this.edge(i, e0!!, e1!!, P, P.size); } else if (this.contains(i, (bounds.xmin + bounds.xmax) / 2, (bounds.ymin + bounds.ymax) / 2)) { return mutableListOf( bounds.xmax, @@ -394,22 +398,19 @@ class Voronoi(val delaunay: Delaunay, val bounds: Rectangle) { when { (c and 0b1000) != 0 -> { x = nx0 + (nx1 - nx0) * (bounds.ymax - ny0) / (ny1 - ny0) - y = bounds.ymax + y = bounds.ymax; } - (c and 0b0100) != 0 -> { x = nx0 + (nx1 - nx0) * (bounds.ymin - ny0) / (ny1 - ny0) y = bounds.ymin } - (c and 0b0010) != 0 -> { y = ny0 + (ny1 - ny0) * (bounds.xmax - nx0) / (nx1 - nx0) x = bounds.xmax } - else -> { y = ny0 + (ny1 - ny0) * (bounds.xmin - nx0) / (nx1 - nx0) - x = bounds.xmin + x = bounds.xmin; } } @@ -457,40 +458,33 @@ class Voronoi(val delaunay: Delaunay, val bounds: Rectangle) { e = 0b0100 continue@loop } - 0b0100 -> { // top e = 0b0110 x = bounds.xmax y = bounds.ymin } - 0b0110 -> { // top-right e = 0b0010 continue@loop } - 0b0010 -> { // right e = 0b1010 x = bounds.xmax y = bounds.ymax } - 0b1010 -> { // bottom-right e = 0b1000 continue@loop } - 0b1000 -> { // bottom e = 0b1001 x = bounds.xmin y = bounds.ymax } - 0b1001 -> { // bottom-left e = 0b0001 continue@loop } - 0b0001 -> { // left e = 0b0101 x = bounds.xmin @@ -630,4 +624,3 @@ private val Rectangle.ymax: Double private fun Double?.isFalsy() = this == null || this == -0.0 || this == 0.0 || isNaN() -