diff --git a/orx-triangulation/build.gradle.kts b/orx-triangulation/build.gradle.kts index 64844493..02702077 100644 --- a/orx-triangulation/build.gradle.kts +++ b/orx-triangulation/build.gradle.kts @@ -33,6 +33,7 @@ kotlin { dependencies { api(libs.openrndr.math) api(libs.openrndr.shape) + implementation(project(":orx-noise")) } } diff --git a/orx-triangulation/src/commonMain/kotlin/Delaunay.kt b/orx-triangulation/src/commonMain/kotlin/Delaunay.kt index 10fec43c..726d2ac6 100644 --- a/orx-triangulation/src/commonMain/kotlin/Delaunay.kt +++ b/orx-triangulation/src/commonMain/kotlin/Delaunay.kt @@ -2,10 +2,6 @@ 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 @@ -77,13 +73,13 @@ class Delaunay(val points: DoubleArray) { init() } - fun neighbors(i:Int) = sequence { + fun neighbors(i: Int) = sequence { val e0 = inedges.getOrNull(i) ?: return@sequence if (e0 != -1) { var e = e0 - var p0 = -1 + var p0 : Int - loop@do { + loop@ do { p0 = triangles[e] yield(p0) e = if (e % 3 == 2) e - 2 else e + 1 @@ -109,29 +105,31 @@ class Delaunay(val points: DoubleArray) { } fun collinear(): Boolean { - for (i in 0 until triangles.size step 3) { + for (i in triangles.indices 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 0 until points.size step 2) { - val p = jitter(points[i], points[i+1], r) + for (i in points.indices 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) diff --git a/orx-triangulation/src/commonMain/kotlin/DelaunayTriangulation.kt b/orx-triangulation/src/commonMain/kotlin/DelaunayTriangulation.kt index 73a0b314..25c5e23b 100644 --- a/orx-triangulation/src/commonMain/kotlin/DelaunayTriangulation.kt +++ b/orx-triangulation/src/commonMain/kotlin/DelaunayTriangulation.kt @@ -14,28 +14,43 @@ class DelaunayTriangulation(val points: List) { fun voronoiDiagram(bounds: Rectangle) = VoronoiDiagram(this, bounds) - fun neighbors(pointIndex: Int) : Sequence { + fun neighbors(pointIndex: Int): Sequence { return delaunay.neighbors(pointIndex) } - fun neighborPoints(pointIndex: Int) : List { + fun neighborPoints(pointIndex: Int): List { return neighbors(pointIndex).map { points[it] }.toList() } - fun triangles(): List { + fun triangleIndices(): List { + val list = mutableListOf() + for (i in delaunay.triangles.indices step 3) { + list.add( + intArrayOf( + delaunay.triangles[i], + delaunay.triangles[i + 1], + delaunay.triangles[i + 2] + ) + ) + } + return list + } + + fun triangles(filterPredicate: (Int, Int, Int) -> Boolean = { _, _, _ -> true }): List { val list = mutableListOf() - for (i in delaunay.triangles.indices step 3 ) { + for (i in delaunay.triangles.indices step 3) { val t0 = delaunay.triangles[i] val t1 = delaunay.triangles[i + 1] val t2 = delaunay.triangles[i + 2] - val p1 = points[t0] - val p2 = points[t1] - val p3 = points[t2] - // originally they are defined *counterclockwise* - list.add(Triangle(p3, p2, p1)) + if (filterPredicate(t2, t1, t0)) { + val p1 = points[t0] + val p2 = points[t1] + val p3 = points[t2] + list.add(Triangle(p3, p2, p1)) + } } return list } @@ -61,11 +76,11 @@ class DelaunayTriangulation(val points: List) { close() } - fun nearest(query: Vector2) : Int = delaunay.find(query.x, query.y) + fun nearest(query: Vector2): Int = delaunay.find(query.x, query.y) - fun nearestPoint(query: Vector2) : Vector2 = points[nearest(query)] + fun nearestPoint(query: Vector2): Vector2 = points[nearest(query)] } -fun List.delaunayTriangulation() : DelaunayTriangulation { +fun List.delaunayTriangulation(): DelaunayTriangulation { return DelaunayTriangulation(this) } \ No newline at end of file diff --git a/orx-triangulation/src/commonMain/kotlin/DoubleDouble.kt b/orx-triangulation/src/commonMain/kotlin/DoubleDouble.kt index d8650e56..58defed3 100644 --- a/orx-triangulation/src/commonMain/kotlin/DoubleDouble.kt +++ b/orx-triangulation/src/commonMain/kotlin/DoubleDouble.kt @@ -17,8 +17,8 @@ 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; + val x = a - b + val y = (a - x) - b return doubleArrayOf(y, x) } @@ -60,11 +60,11 @@ internal fun reduceSignificand( 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; +private const val f = 134217729 // 2**27 + 1; /** @@ -90,9 +90,9 @@ private const val f = 134217729; // 2**27 + 1; * @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; + val c = f * a + val a_h = c - (c - a) + val a_l = a - a_h return doubleArrayOf(a_h, a_l) } @@ -104,9 +104,9 @@ private fun split(a: Double): DoubleArray { * @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); + val x = a - b + val bvirt = a - x + val y = (a - (x + bvirt)) + (bvirt - b) return doubleArrayOf(y, x) } @@ -130,15 +130,15 @@ internal 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); @@ -149,14 +149,14 @@ internal fun twoProduct(a: Double, b: Double): DoubleArray { } internal fun twoSquare(a: Double): DoubleArray { - val x = a * a; + 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) } @@ -174,8 +174,8 @@ 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; + val x = a + b + val bv = x - a return doubleArrayOf((a - (x - bv)) + (b - bv), x) } @@ -194,21 +194,21 @@ internal fun twoSum(a: Double, b: Double): DoubleArray { * @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]; + 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) } @@ -229,19 +229,19 @@ internal 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) } @@ -261,21 +261,21 @@ internal fun ddMultDd(x: DoubleArray, y: DoubleArray): DoubleArray { * @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]; + 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) } @@ -297,24 +297,24 @@ internal fun ddAddDd(x: DoubleArray, y: DoubleArray): DoubleArray { * @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]; + 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/Predicates.kt b/orx-triangulation/src/commonMain/kotlin/Predicates.kt index 265c9f83..7b5e129d 100644 --- a/orx-triangulation/src/commonMain/kotlin/Predicates.kt +++ b/orx-triangulation/src/commonMain/kotlin/Predicates.kt @@ -1,6 +1,6 @@ package org.openrndr.extra.triangulation -internal fun orient2d(bx: Double, by: Double, ax: Double, ay: Double, cx: Double, cy: Double): Double { +fun orient2d(bx: Double, by: Double, ax: Double, ay: Double, cx: Double, cy: Double): Double { // (ax,ay) (bx,by) are swapped such that the sign of the determinant is flipped. which is what Delaunator.kt expects. /* @@ -8,12 +8,12 @@ internal fun orient2d(bx: Double, by: Double, ax: Double, ay: Double, cx: Double | c d | | bx - cx by - cy | */ - val a = ax - cx - val b = ay - cy - val c = bx - cx - val d = by - cy + val a = twoDiff(ax, cx) + val b = twoDiff(ay, cy) + val c = twoDiff(bx, cx) + val d = twoDiff(by, cy) - val determinant = ddDiffDd(twoProduct(a, d), twoProduct(b, c)) + val determinant = ddDiffDd(ddMultDd(a, d), ddMultDd(b, c)) return determinant[1] } \ No newline at end of file diff --git a/orx-triangulation/src/commonMain/kotlin/SmoothScatter.kt b/orx-triangulation/src/commonMain/kotlin/SmoothScatter.kt new file mode 100644 index 00000000..2c28613f --- /dev/null +++ b/orx-triangulation/src/commonMain/kotlin/SmoothScatter.kt @@ -0,0 +1,120 @@ +package org.openrndr.extra.triangulation + +import org.openrndr.extra.noise.scatter +import org.openrndr.math.Vector2 +import org.openrndr.shape.ShapeProvider +import org.openrndr.shape.bounds +import kotlin.random.Random + +fun ShapeProvider.smoothScatterSeq( + placementRadius: Double, + distanceToEdge: Double = placementRadius * 2.0, + smoothing: Double = 0.5, + random: Random = Random.Default +) = sequence { + val boundaryPointSets = this@smoothScatterSeq.shape.contours.map { + it.equidistantPositions((it.length / placementRadius).toInt()) + } + + val boundaryPoints = boundaryPointSets.flatten() + val interiorPoints = this@smoothScatterSeq.shape.scatter( + placementRadius = placementRadius, distanceToEdge = distanceToEdge, random = random + ) + + val bounds = interiorPoints.bounds.offsetEdges(100.0) + var relaxedPoints = interiorPoints + + while (true) { + val dt = (relaxedPoints + boundaryPoints) + val v = dt.voronoiDiagram(bounds) + + relaxedPoints = relaxedPoints.mapIndexed { index, it -> + val c = v.cellCentroid(index) + if (c.x == c.x && c.y == c.y) { + it * smoothing + c * (1.0 - smoothing) + } else { + it + } + } + yield(relaxedPoints) + } +} + +fun ShapeProvider.smoothScatterWeightedSeq( + placementRadius: Double, + distanceToEdge: Double = placementRadius * 2.0, + smoothing: Double = 0.5, + random: Random = Random.Default +) = sequence { + val boundaryPointSets = this@smoothScatterWeightedSeq.shape.contours.map { + it.equidistantPositions((it.length / placementRadius).toInt()) + } + + val boundaryPoints = boundaryPointSets.flatten() + val interiorPoints = this@smoothScatterWeightedSeq.shape.scatter( + placementRadius = placementRadius, distanceToEdge = distanceToEdge, random = random + ) + + val bounds = interiorPoints.bounds.offsetEdges(100.0) + var relaxedPoints = interiorPoints + + fun isBoundaryPoint(i: Int) = i >= interiorPoints.size + + val targetAreas = interiorPoints.map { if (random.nextDouble() < 0.1) 450.0 else null } + + while (true) { + val dt = (relaxedPoints + boundaryPoints) + val v = dt.voronoiDiagram(bounds) + + + relaxedPoints = relaxedPoints.mapIndexed { index, it -> + val c = v.cellCentroid(index) + if (c.x == c.x && c.y == c.y) { + it * smoothing + c * (1.0 - smoothing) + } else { + it + } + } + val resolvedPoints = relaxedPoints.map { it }.toMutableList() + + + for (i in interiorPoints.indices) { + + if (targetAreas[i] != null) { + + val targetArea = targetAreas[i]!! + val cellArea = v.cellArea(i) + val cellCentroid = v.cellCentroid(i) + val areaDiff = targetArea - cellArea + + val ns = v.neighbors(i).filter { !isBoundaryPoint(it) }.toList() + + var force: Vector2 + val scale = 1.0 / ns.size + for (n in ns) { + force = v.cellCentroid(n) - cellCentroid + resolvedPoints[n] += force.normalized * (areaDiff * 0.01) * scale + } + } + relaxedPoints = resolvedPoints + } + yield(relaxedPoints) + } +} + + +fun ShapeProvider.smoothScatter( + placementRadius: Double, + distanceToEdge: Double = placementRadius * 2.0, + iterations: Int = 10, + smoothing: Double = 0.5, + random: Random = Random.Default +): List { + + val seq = smoothScatterSeq(placementRadius, distanceToEdge, smoothing, random).iterator() + + for (i in 0 until iterations - 1) { + seq.next() + } + return seq.next() +} \ 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 b665cca8..7f3a873d 100644 --- a/orx-triangulation/src/commonMain/kotlin/Voronoi.kt +++ b/orx-triangulation/src/commonMain/kotlin/Voronoi.kt @@ -154,11 +154,7 @@ 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 @@ -195,8 +191,8 @@ class Voronoi(val delaunay: Delaunay, val bounds: Rectangle) { if (cj != null) { val li = ci.size val lj = cj.size - loop@ for (ai in 0 until ci.size step 2) { - for (aj in 0 until cj.size step 2) { + loop@ for (ai in ci.indices step 2) { + for (aj in cj.indices step 2) { if (ci[ai] == cj[aj] && ci[ai + 1] == cj[aj + 1] && ci[(ai + 2) % li] == cj[(aj + lj - 2) % lj] @@ -256,10 +252,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 @@ -298,7 +294,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? = null + var e0: Int? var e1: Int? = 0 for (j in 0 until n step 2) { @@ -360,7 +356,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, @@ -398,19 +394,22 @@ 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 } } @@ -458,33 +457,40 @@ 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 diff --git a/orx-triangulation/src/commonMain/kotlin/VoronoiDiagram.kt b/orx-triangulation/src/commonMain/kotlin/VoronoiDiagram.kt index bafb164d..a00b14fe 100644 --- a/orx-triangulation/src/commonMain/kotlin/VoronoiDiagram.kt +++ b/orx-triangulation/src/commonMain/kotlin/VoronoiDiagram.kt @@ -20,6 +20,37 @@ class VoronoiDiagram(val delaunayTriangulation: DelaunayTriangulation, val bound } } + fun cellArea(i: Int, contour: ShapeContour = cellPolygon(i)): Double { + val segments = contour.segments + var sum = 0.0 + for (j in segments.indices) { + val v0 = segments[j].start + val v1 = segments[(j + 1).mod(segments.size)].start + sum += v0.x * v1.y - v1.x * v0.y + } + return sum / 2.0 + } + + fun cellCentroid(i: Int, contour: ShapeContour = cellPolygon(i)): Vector2 { + val segments = cellPolygon(i).segments + var cx = 0.0 + var cy = 0.0 + for (j in segments.indices) { + val v0 = segments[j].start + val v1 = segments[(j + 1).mod(segments.size)].start + cx += (v0.x + v1.x) * (v0.x * v1.y - v1.x * v0.y) + cy += (v0.y + v1.y) * (v0.x * v1.y - v1.x * v0.y) + } + val a = cellArea(i, contour) * 6.0 + cx /= a + cy /= a + return Vector2(cx, cy) + } + + fun cellCentroids() = (delaunayTriangulation.points.indices).map { + cellCentroid(it) + } + fun cellPolygon(i: Int): ShapeContour { val points = voronoi.clip(i)