[orx-triangulation] Small fixes

This commit is contained in:
Edwin Jakobs
2022-12-19 23:49:29 +01:00
parent 608ef6e33a
commit 163218378d
4 changed files with 137 additions and 122 deletions

View File

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

View File

@@ -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,11 +77,11 @@ class Delaunay(val points: DoubleArray) {
init()
}
fun neighbors(i: Int) = sequence {
val e0 = inedges.getOrNull(i) ?: return@sequence
fun neighbors(i:Int) = sequence<Int> {
val e0 = inedges[i]
if (e0 != -1) {
var e = e0
var p0 : Int
var p0 = -1
loop@do {
p0 = triangles[e]
@@ -105,28 +109,26 @@ 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 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
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)
}
fun init() {
if (hull.size > 2 && collinear()) {
println("warning: triangulation is collinear")
val r = 1E-8
for (i in points.indices step 2) {
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]
@@ -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)
}

View File

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

View File

@@ -154,7 +154,11 @@ class Voronoi(val delaunay: Delaunay, val bounds: Rectangle) {
}
private fun cell(i: Int): MutableList<Double>? {
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()