kaze/orx
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

实现等高线的绘制

Browse Source
This commit is contained in:
tabidachinokaze
2025-11-25 19:43:51 +08:00
parent de15029b2b
commit 816e954ed8
31 changed files with 3553 additions and 257 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4
android/build.gradle.kts
View File

@@ -57,6 +57,10 @@ dependencies {
implementation(libs.androidx.lifecycle.runtime.ktx)
implementation(project(":icegps-common"))
implementation(project(":icegps-shared"))
implementation(project(":orx-marching-squares"))
implementation(project(":orx-palette")) {
exclude(group = "org.openrndr", module = "openrndr-draw")
}
testImplementation(libs.junit)
androidTestImplementation(libs.ext.junit)

401
android/src/main/java/com/icegps/orx/ContoursManager.kt Normal file
View File

@@ -0,0 +1,401 @@
package com.icegps.orx
import ColorBrewer2Type
import android.content.Context
import android.util.Log
import colorBrewer2Palettes
import com.icegps.math.geometry.Vector3D
import com.icegps.orx.ktx.area
import com.icegps.orx.ktx.toColorInt
import com.icegps.orx.ktx.toMapboxPoint
import com.icegps.orx.ktx.toast
import com.icegps.shared.ktx.TAG
import com.mapbox.geojson.Feature
import com.mapbox.geojson.FeatureCollection
import com.mapbox.geojson.LineString
import com.mapbox.geojson.Polygon
import com.mapbox.maps.MapView
import com.mapbox.maps.Style
import com.mapbox.maps.extension.style.expressions.generated.Expression
import com.mapbox.maps.extension.style.layers.addLayer
import com.mapbox.maps.extension.style.layers.generated.fillLayer
import com.mapbox.maps.extension.style.layers.generated.lineLayer
import com.mapbox.maps.extension.style.layers.properties.generated.LineCap
import com.mapbox.maps.extension.style.layers.properties.generated.LineJoin
import com.mapbox.maps.extension.style.sources.addSource
import com.mapbox.maps.extension.style.sources.generated.geoJsonSource
import kotlinx.coroutines.CoroutineScope
import kotlinx.coroutines.Dispatchers
import kotlinx.coroutines.async
import kotlinx.coroutines.awaitAll
import kotlinx.coroutines.launch
import kotlinx.coroutines.withContext
import com.icegps.orx.triangulation.DelaunayTriangulation3D
import com.icegps.orx.triangulation.Triangle3D
import org.openrndr.math.Vector2
import org.openrndr.math.Vector3
import org.openrndr.shape.Rectangle
import org.openrndr.shape.ShapeContour
import kotlin.math.max
class ContoursManager(
private val context: Context,
private val mapView: MapView,
private val scope: CoroutineScope
) {
private val sourceId: String = "contours-source-id-10"
private val layerId: String = "contours-layer-id-10"
private val fillSourceId: String = "contours-fill-source-id-10"
private val fillLayerId: String = "contours-fill-layer-id-10"
private val gridSourceId: String = "grid-polygon-source-id"
private val gridLayerId: String = "grid-polygon-layer-id"
private var contourSize: Int = 6
private var heightRange: ClosedFloatingPointRange<Double> = 0.0..100.0
private var cellSize: Double? = 10.0
private val simplePalette = SimplePalette(
range = 0.0..100.0
)
private var colors = colorBrewer2Palettes(
numberOfColors = contourSize,
paletteType = ColorBrewer2Type.Any
).first().colors.reversed()
private var points: List<Vector3D> = emptyList()
private val polylineManager = PolylineManager(mapView)
fun updateContourSize(contourSize: Int) {
this.contourSize = contourSize
colors = colorBrewer2Palettes(
numberOfColors = contourSize,
paletteType = ColorBrewer2Type.Any
).first().colors.reversed()
}
fun updateCellSize(value: Double) {
cellSize = value
}
fun updatePoints(
points: List<Vector3D>,
) {
this.points = points
}
fun updateHeightRange(
heightRange: ClosedFloatingPointRange<Double>? = null
) {
if (heightRange == null) {
if (points.isEmpty()) {
return
}
val height = points.map { it.z }
val range = height.min()..height.max()
this.heightRange = range
simplePalette.setRange(range)
} else {
this.heightRange = heightRange
simplePalette.setRange(heightRange)
}
}
private var isGridVisible: Boolean = true
private var gridModel: GridModel? = null
fun setGridVisible(visible: Boolean) {
if (visible != isGridVisible) {
isGridVisible = visible
if (visible) {
if (gridModel != null) mapView.displayGridModel(
grid = gridModel!!,
sourceId = gridSourceId,
layerId = gridLayerId,
palette = simplePalette::palette
)
} else {
mapView.mapboxMap.getStyle { style ->
try {
style.removeStyleLayer(gridLayerId)
} catch (_: Exception) {
}
if (style.styleSourceExists(gridSourceId)) {
style.removeStyleSource(gridSourceId)
}
}
}
}
}
private var triangles: List<Triangle3D> = listOf()
private var isTriangleVisible: Boolean = true
fun setTriangleVisible(visible: Boolean) {
if (visible != isTriangleVisible) {
isTriangleVisible = visible
if (visible) {
polylineManager.update(
triangles.map {
listOf(it.x1, it.x2, it.x3)
.map { Vector3D(it.x, it.y, it.z) }
}
)
} else {
polylineManager.clearContours()
}
}
}
fun refresh() {
val points = points
if (points.size <= 3) {
context.toast("points size ${points.size}")
return
}
scope.launch {
mapView.mapboxMap.getStyle { style ->
val step = heightRange.endInclusive / contourSize
val zip = (0..contourSize).map { index ->
heightRange.start + index * step
}.zipWithNext { a, b -> a..b }
val points = points.map { Vector3(it.x, it.y, it.z) }
val area = points.area
val triangulation = DelaunayTriangulation3D(points)
val triangles: MutableList<Triangle3D> = triangulation.triangles()
val cellSize: Double = if (cellSize == null || cellSize!! < 0.1) {
(max(triangulation.points.area.width, triangulation.points.area.height) / 50)
} else {
cellSize!!
}
scope.launch {
val gridModel = triangulationToGrid(
delaunator = triangulation,
cellSize = cellSize,
)
this@ContoursManager.gridModel = gridModel
if (isGridVisible) mapView.displayGridModel(
grid = gridModel,
sourceId = gridSourceId,
layerId = gridLayerId,
palette = simplePalette::palette
)
}
scope.launch(Dispatchers.Default) {
val lineFeatures = mutableListOf<List<Feature>>()
val features = zip.mapIndexed { index, range ->
async {
val contours = findContours(
triangles = triangles,
range = range,
area = area,
cellSize = cellSize
)
val color = colors[index].toColorInt()
lineFeatures.add(contoursToLineFeatures(contours, color).flatten())
contoursToPolygonFeatures(contours, color)
}
}.awaitAll()
withContext(Dispatchers.Main) {
if (false) setupLineLayer(
style = style,
sourceId = sourceId,
layerId = layerId,
features = lineFeatures.flatten()
)
setupFillLayer(
style = style,
sourceId = fillSourceId,
layerId = fillLayerId,
features = features.filterNotNull(),
)
Log.d(TAG, "refresh: 刷新完成")
}
}
}
}
}
fun findContours(
triangles: MutableList<Triangle3D>,
range: ClosedFloatingPointRange<Double>,
area: Rectangle,
cellSize: Double
): List<ShapeContour> {
return org.openrndr.extra.marchingsquares.findContours(
f = { v ->
val triangle = triangles.firstOrNull { triangle ->
isPointInTriangle3D(v, listOf(triangle.x1, triangle.x2, triangle.x3))
}
(triangle?.let { triangle ->
val interpolate = interpolateHeight(
point = v,
triangle = listOf(
triangle.x1,
triangle.x2,
triangle.x3,
)
)
if (interpolate.z in range) -1.0
else 1.0
} ?: 1.0).also {
Log.d(TAG, "findContours: ${v} -> ${it}")
}
},
area = area,
cellSize = cellSize,
)
}
private fun setupLineLayer(
style: Style,
sourceId: String,
layerId: String,
features: List<Feature>
) {
style.removeStyleLayer(layerId)
style.removeStyleSource(sourceId)
val source = geoJsonSource(sourceId) {
featureCollection(FeatureCollection.fromFeatures(features))
}
style.addSource(source)
val layer = lineLayer(layerId, sourceId) {
lineColor(Expression.Companion.toColor(Expression.Companion.get("color"))) // 从属性获取颜色
lineWidth(1.0)
lineCap(LineCap.Companion.ROUND)
lineJoin(LineJoin.Companion.ROUND)
lineOpacity(0.8)
}
style.addLayer(layer)
}
private fun setupFillLayer(
style: Style,
sourceId: String,
layerId: String,
features: List<Feature>
) {
style.removeStyleLayer(layerId)
style.removeStyleSource(sourceId)
val source = geoJsonSource(sourceId) {
featureCollection(FeatureCollection.fromFeatures(features))
}
style.addSource(source)
val layer = fillLayer(layerId, sourceId) {
fillColor(Expression.Companion.toColor(Expression.Companion.get("color"))) // 从属性获取颜色
fillOpacity(0.5)
fillAntialias(true)
}
style.addLayer(layer)
}
fun contoursToLineFeatures(contours: List<ShapeContour>, color: Int): List<List<Feature>> {
return contours.drop(1).map { contour ->
contour.segments.map { segment ->
LineString.fromLngLats(listOf(segment.start.toMapboxPoint(), segment.end.toMapboxPoint()))
}.map { lineString ->
Feature.fromGeometry(lineString).apply {
// 将颜色Int转换为十六进制字符串
addStringProperty("color", color.toHexColorString())
}
}
}
}
fun contoursToPolygonFeatures(contours: List<ShapeContour>, color: Int): Feature? {
val lists = contours.drop(0).filter { it.segments.isNotEmpty() }.map { contour ->
val start = contour.segments[0].start
listOf(start) + contour.segments.map { it.end }
}.map { points -> points.map { it.toMapboxPoint() } }
if (lists.isEmpty()) {
Log.w(TAG, "contoursToPolygonFeatures: 没有有效的轮廓数据")
return null
}
val polygon = Polygon.fromLngLats(lists)
return Feature.fromGeometry(polygon).apply {
// 将颜色Int转换为十六进制字符串
addStringProperty("color", color.toHexColorString())
}
}
fun Int.toHexColorString(): String {
return String.format("#%06X", 0xFFFFFF and this)
}
fun clearContours() {
mapView.mapboxMap.getStyle { style ->
try {
style.removeStyleLayer(layerId)
} catch (_: Exception) {
}
try {
style.removeStyleSource(sourceId)
} catch (_: Exception) {
}
}
}
}
fun isPointInTriangle3D(point: Vector2, triangle: List<Vector3>): Boolean {
require(triangle.size == 3) { "三角形必须有3个顶点" }
val (v1, v2, v3) = triangle
// 计算重心坐标
val denominator = (v2.y - v3.y) * (v1.x - v3.x) + (v3.x - v2.x) * (v1.y - v3.y)
if (denominator == 0.0) return false // 退化三角形
val alpha = ((v2.y - v3.y) * (point.x - v3.x) + (v3.x - v2.x) * (point.y - v3.y)) / denominator
val beta = ((v3.y - v1.y) * (point.x - v3.x) + (v1.x - v3.x) * (point.y - v3.y)) / denominator
val gamma = 1.0 - alpha - beta
// 点在三角形内当且仅当所有重心坐标都在[0,1]范围内
return alpha >= 0 && beta >= 0 && gamma >= 0 &&
alpha <= 1 && beta <= 1 && gamma <= 1
}
/**
* 使用重心坐标计算点在三角形上的高度
*
* @param point 二维点 (x, y)
* @param triangle 三角形的三个顶点
* @return 三维点 (x, y, z)
*/
fun interpolateHeight(point: Vector2, triangle: List<Vector3>): Vector3 {
/**
* 计算点在三角形中的重心坐标
*/
fun calculateBarycentricCoordinates(
point: Vector2,
v1: Vector3,
v2: Vector3,
v3: Vector3
): Triple<Double, Double, Double> {
val denom = (v2.y - v3.y) * (v1.x - v3.x) + (v3.x - v2.x) * (v1.y - v3.y)
val alpha = ((v2.y - v3.y) * (point.x - v3.x) + (v3.x - v2.x) * (point.y - v3.y)) / denom
val beta = ((v3.y - v1.y) * (point.x - v3.x) + (v1.x - v3.x) * (point.y - v3.y)) / denom
val gamma = 1.0 - alpha - beta
return Triple(alpha, beta, gamma)
}
require(triangle.size == 3) { "三角形必须有3个顶点" }
val (v1, v2, v3) = triangle
// 计算重心坐标
val (alpha, beta, gamma) = calculateBarycentricCoordinates(point, v1, v2, v3)
// 使用重心坐标插值z值
val z = alpha * v1.z + beta * v2.z + gamma * v3.z
return Vector3(point.x, point.y, z)
}

53
android/src/main/java/com/icegps/orx/CoordinateGenerator.kt Normal file
View File

@@ -0,0 +1,53 @@
package com.icegps.orx
import com.icegps.math.geometry.Angle
import com.icegps.math.geometry.Vector3D
import com.icegps.math.geometry.degrees
import kotlin.math.cos
import kotlin.math.sin
import kotlin.random.Random
/**
* @author tabidachinokaze
* @date 2025/11/25
*/
fun coordinateGenerate(): List<Vector3D> {
val minX = -20.0
val maxX = 20.0
val minY = -20.0
val maxY = 20.0
val minZ = -20.0
val maxZ = 20.0
val x: () -> Double = { Random.nextDouble(minX, maxX) }
val y: () -> Double = { Random.nextDouble(minY, maxY) }
val z: () -> Double = { Random.nextDouble(minZ, maxZ) }
val dPoints = (0..60).map {
Vector3D(x(), y(), z())
}
return dPoints
}
fun coordinateGenerate1(): List<List<Vector3D>> {
/**
* 绕 Z 轴旋转指定角度（弧度）
*/
fun Vector3D.rotateAroundZ(angle: Angle): Vector3D {
val cosAngle = cos(angle.radians)
val sinAngle = sin(angle.radians)
return Vector3D(
x = x * cosAngle - y * sinAngle,
y = x * sinAngle + y * cosAngle,
z = z
)
}
val center = Vector3D()
val direction = Vector3D(0.0, 1.0, -1.0)
return (0..360).step(10).map {
val nowDirection = direction.rotateAroundZ(it.degrees)
listOf(2, 6, 10).map {
center + nowDirection * it
}
}
}

83
android/src/main/java/com/icegps/orx/GridDisplay.kt Normal file
View File

@@ -0,0 +1,83 @@
package com.icegps.orx
import com.icegps.common.helper.GeoHelper
import com.icegps.math.geometry.Vector2D
import com.mapbox.geojson.Feature
import com.mapbox.geojson.FeatureCollection
import com.mapbox.geojson.Point
import com.mapbox.geojson.Polygon
import com.mapbox.maps.MapView
import com.mapbox.maps.extension.style.expressions.generated.Expression
import com.mapbox.maps.extension.style.layers.addLayer
import com.mapbox.maps.extension.style.layers.generated.FillLayer
import com.mapbox.maps.extension.style.sources.addSource
import com.mapbox.maps.extension.style.sources.generated.geoJsonSource
/**
* @author tabidachinokaze
* @date 2025/11/25
*/
fun MapView.displayGridModel(
grid: GridModel,
sourceId: String,
layerId: String,
palette: (Double?) -> String,
) {
val geoHelper = GeoHelper.getSharedInstance()
mapboxMap.getStyle { style ->
val polygonFeatures = mutableListOf<Feature>()
val minX = grid.minX
val maxY = grid.maxY
val cellSize = grid.cellSize
for (r in 0 until grid.rows) {
for (c in 0 until grid.cols) {
val idx = r * grid.cols + c
val v = grid.cells[idx] ?: continue
val x0 = minX + c * cellSize
val y0 = maxY - r * cellSize
val x1 = x0 + cellSize
val y1 = y0 - cellSize
val ring = listOf(
Vector2D(x0, y0),
Vector2D(x1, y0),
Vector2D(x1, y1),
Vector2D(x0, y1),
Vector2D(x0, y0),
).map {
geoHelper.enuToWGS84Object(GeoHelper.ENU(it.x, it.y))
}.map {
Point.fromLngLat(it.lon, it.lat)
}
val poly = Polygon.fromLngLats(listOf(ring))
val polyFeature = Feature.fromGeometry(poly)
polyFeature.addStringProperty("color", palette(v))
polyFeature.addNumberProperty("value", v ?: -9999.0)
polygonFeatures.add(polyFeature)
}
}
try {
style.removeStyleLayer(layerId)
} catch (_: Exception) {
}
if (style.styleSourceExists(sourceId)) {
style.removeStyleSource(sourceId)
}
val polygonSource = geoJsonSource(sourceId) {
featureCollection(FeatureCollection.fromFeatures(polygonFeatures))
}
style.addSource(polygonSource)
val fillLayer = FillLayer(layerId, sourceId).apply {
fillColor(Expression.toColor(Expression.get("color")))
fillOpacity(0.5)
}
style.addLayer(fillLayer)
}
}

132
android/src/main/java/com/icegps/orx/GridModel.kt Normal file
View File

@@ -0,0 +1,132 @@
package com.icegps.orx
import com.icegps.math.geometry.Vector2D
import com.icegps.orx.triangulation.DelaunayTriangulation3D
import org.openrndr.math.Vector3
import kotlin.math.absoluteValue
import kotlin.math.ceil
/**
* @author tabidachinokaze
* @date 2025/11/25
*/
data class GridModel(
val minX: Double,
val maxX: Double,
val minY: Double,
val maxY: Double,
val rows: Int,
val cols: Int,
val cellSize: Double,
val cells: Array<Double?>
)
fun triangulationToGrid(
delaunator: DelaunayTriangulation3D,
cellSize: Double = 50.0,
maxSidePixels: Int = 5000
): GridModel {
fun pointInTriangle(pt: Vector2D, a: Vector3, b: Vector3, c: Vector3): Boolean {
val v0x = c.x - a.x
val v0y = c.y - a.y
val v1x = b.x - a.x
val v1y = b.y - a.y
val v2x = pt.x - a.x
val v2y = pt.y - a.y
val dot00 = v0x * v0x + v0y * v0y
val dot01 = v0x * v1x + v0y * v1y
val dot02 = v0x * v2x + v0y * v2y
val dot11 = v1x * v1x + v1y * v1y
val dot12 = v1x * v2x + v1y * v2y
val denom = dot00 * dot11 - dot01 * dot01
if (denom == 0.0) return false
val invDenom = 1.0 / denom
val u = (dot11 * dot02 - dot01 * dot12) * invDenom
val v = (dot00 * dot12 - dot01 * dot02) * invDenom
return u >= 0 && v >= 0 && u + v <= 1
}
fun barycentricInterpolateLegacy(pt: Vector2D, a: Vector3, b: Vector3, c: Vector3, values: DoubleArray): Double {
val area = { p1: Vector2D, p2: Vector3, p3: Vector3 ->
((p2.x - p1.x) * (p3.y - p1.y) - (p3.x - p1.x) * (p2.y - p1.y)).absoluteValue / 2.0
}
val area2 = { p1: Vector3, p2: Vector3, p3: Vector3 ->
((p2.x - p1.x) * (p3.y - p1.y) - (p3.x - p1.x) * (p2.y - p1.y)).absoluteValue / 2.0
}
val areaTotal = area2(a, b, c)
if (areaTotal == 0.0) return values[0]
val wA = area(pt, b, c) / areaTotal
val wB = area(pt, c, a) / areaTotal
val wC = area(pt, a, b) / areaTotal
return values[0] * wA + values[1] * wB + values[2] * wC
}
val pts = delaunator.points
require(pts.isNotEmpty()) { "points empty" }
val x = pts.map { it.x }
val y = pts.map { it.y }
val minX = x.min()
val maxX = x.max()
val minY = y.min()
val maxY = y.max()
val width = maxX - minX
val height = maxY - minY
var cols = ceil(width / cellSize).toInt()
var rows = ceil(height / cellSize).toInt()
// 防止过大
if (cols > maxSidePixels) cols = maxSidePixels
if (rows > maxSidePixels) rows = maxSidePixels
val cells = Array<Double?>(rows * cols) { null }
val triangles = delaunator.triangles()
for (ti in 0 until triangles.size) {
val (a, b, c) = triangles[ti]
val tminX = minOf(a.x, b.x, c.x)
val tmaxX = maxOf(a.x, b.x, c.x)
val tminY = minOf(a.y, b.y, c.y)
val tmaxY = maxOf(a.y, b.y, c.y)
val colMin = ((tminX - minX) / cellSize).toInt().coerceIn(0, cols - 1)
val colMax = ((tmaxX - minX) / cellSize).toInt().coerceIn(0, cols - 1)
val rowMin = ((maxY - tmaxY) / cellSize).toInt().coerceIn(0, rows - 1)
val rowMax = ((maxY - tminY) / cellSize).toInt().coerceIn(0, rows - 1)
val triVertexVals = doubleArrayOf(a.z, b.z, c.z)
for (r in rowMin..rowMax) {
for (cIdx in colMin..colMax) {
val centerX = minX + (cIdx + 0.5) * cellSize
val centerY = maxY - (r + 0.5) * cellSize
val pt = Vector2D(centerX, centerY)
if (pointInTriangle(pt, a, b, c)) {
val idx = r * cols + cIdx
val valInterp = barycentricInterpolateLegacy(pt, a, b, c, triVertexVals)
cells[idx] = valInterp
}
}
}
}
val grid = GridModel(
minX = minX,
minY = minY,
maxX = maxX,
maxY = maxY,
rows = rows,
cols = cols,
cellSize = cellSize,
cells = cells
)
return grid
}

318
android/src/main/java/com/icegps/orx/MainActivity.kt
View File

@@ -1,56 +1,31 @@
package com.icegps.orx
import android.graphics.Color
import android.os.Bundle
import android.util.Log
import androidx.activity.enableEdgeToEdge
import androidx.appcompat.app.AppCompatActivity
import androidx.core.view.ViewCompat
import androidx.core.view.WindowInsetsCompat
import androidx.lifecycle.ViewModel
import androidx.lifecycle.viewModelScope
import androidx.lifecycle.ViewModelProvider
import androidx.lifecycle.lifecycleScope
import com.google.android.material.slider.RangeSlider
import com.google.android.material.slider.Slider
import com.icegps.common.helper.GeoHelper
import com.icegps.math.geometry.Angle
import com.icegps.math.geometry.Line3D
import com.icegps.math.geometry.Vector3D
import com.icegps.math.geometry.degrees
import com.icegps.orx.databinding.ActivityMainBinding
import com.icegps.shared.SharedHttpClient
import com.icegps.shared.SharedJson
import com.icegps.shared.api.OpenElevation
import com.icegps.shared.api.OpenElevationApi
import com.icegps.shared.ktx.TAG
import com.icegps.shared.model.GeoPoint
import com.mapbox.geojson.Feature
import com.mapbox.geojson.FeatureCollection
import com.mapbox.geojson.LineString
import com.mapbox.geojson.Point
import com.mapbox.geojson.Polygon
import com.mapbox.maps.CameraOptions
import com.mapbox.maps.MapView
import com.mapbox.maps.Style
import com.mapbox.maps.extension.style.layers.addLayer
import com.mapbox.maps.extension.style.layers.generated.fillLayer
import com.mapbox.maps.extension.style.layers.generated.lineLayer
import com.mapbox.maps.extension.style.layers.properties.generated.LineCap
import com.mapbox.maps.extension.style.layers.properties.generated.LineJoin
import com.mapbox.maps.extension.style.sources.addSource
import com.mapbox.maps.extension.style.sources.generated.geoJsonSource
import kotlinx.coroutines.flow.MutableStateFlow
import kotlinx.coroutines.flow.catch
import com.mapbox.maps.plugin.gestures.addOnMapClickListener
import kotlinx.coroutines.flow.launchIn
import kotlinx.coroutines.flow.map
import kotlinx.coroutines.flow.onEach
import kotlinx.coroutines.flow.update
import org.openrndr.extra.triangulation.DelaunayTriangulation
import org.openrndr.math.Vector2
import org.openrndr.math.YPolarity
import kotlin.math.cos
import kotlin.math.sin
class MainActivity : AppCompatActivity() {
private lateinit var binding: ActivityMainBinding
private lateinit var mapView: MapView
private val viewModel: MainViewModel by lazy {
ViewModelProvider(this)[MainViewModel::class.java]
}
private lateinit var contoursManager: ContoursManager
init {
initGeoHelper()
@@ -78,15 +53,96 @@ class MainActivity : AppCompatActivity() {
)
val points = coordinateGenerate1()
val polygonTest = PolygonTest(mapView)
polygonTest.clear()
val innerPoints = points.map { it[0] }
val outerPoints = points.map { it[1] }
polygonTest.update(
if (false) polygonTest.update(
outer = outerPoints,
inner = innerPoints,
other = points.map { it[2] }
)
// divider
contoursManager = ContoursManager(
context = this,
mapView = mapView,
scope = lifecycleScope
)
val points2 = points.flatten()
contoursManager.updateContourSize(6)
contoursManager.updatePoints(points2)
val height = points2.map { it.z }
val min = height.min()
val max = height.max()
contoursManager.updateHeightRange((min / 2)..max)
binding.heightRange.values = listOf(min.toFloat() / 2, max.toFloat())
binding.heightRange.valueFrom = min.toFloat()
binding.heightRange.valueTo = max.toFloat()
contoursManager.refresh()
binding.sliderTargetHeight.addOnSliderTouchListener(
object : Slider.OnSliderTouchListener {
override fun onStartTrackingTouch(p0: Slider) {
}
override fun onStopTrackingTouch(p0: Slider) {
val present = p0.value / p0.valueTo
// val targetHeight = ((valueRange.endInclusive - valueRange.start) * present) + valueRange.start
// val contours = findContours(triangles, targetHeight)
// contoursTest.clearContours()
// if (false) contoursTest.updateContours(contours)
}
}
)
binding.heightRange.addOnSliderTouchListener(
object : RangeSlider.OnSliderTouchListener {
override fun onStartTrackingTouch(slider: RangeSlider) {
}
override fun onStopTrackingTouch(slider: RangeSlider) {
contoursManager.updateHeightRange((slider.values.min().toDouble() - 1.0)..(slider.values.max().toDouble() + 1.0))
contoursManager.refresh()
}
}
)
binding.switchGrid.setOnCheckedChangeListener { _, isChecked ->
contoursManager.setGridVisible(isChecked)
}
binding.switchTriangle.setOnCheckedChangeListener { _, isChecked ->
contoursManager.setTriangleVisible(isChecked)
}
binding.update.setOnClickListener {
contoursManager.refresh()
}
binding.cellSize.addOnSliderTouchListener(
object : Slider.OnSliderTouchListener {
override fun onStartTrackingTouch(slider: Slider) {
}
override fun onStopTrackingTouch(slider: Slider) {
contoursManager.updateCellSize(slider.value.toDouble())
}
}
)
mapView.mapboxMap.addOnMapClickListener {
viewModel.addPoint(it)
true
}
binding.clearPoints.setOnClickListener {
viewModel.clearPoints()
}
initData()
}
private fun initData() {
viewModel.points.onEach {
contoursManager.updatePoints(it)
contoursManager.updateHeightRange()
}.launchIn(lifecycleScope)
}
}
@@ -100,195 +156,3 @@ fun initGeoHelper(base: GeoPoint = home) {
hgt = base.altitude
)
}
fun fromPoints(
points: List<Vector3D>,
closed: Boolean,
polarity: YPolarity = YPolarity.CW_NEGATIVE_Y
) = if (points.isEmpty()) {
emptyList()
} else {
if (!closed) {
(0 until points.size - 1).map {
Line3D(
points[it],
points[it + 1]
)
}
} else {
val d = (points.last() - points.first()).length
val usePoints = if (d > 1E-6) points else points.dropLast(1)
(usePoints.indices).map {
Line3D(
usePoints[it],
usePoints[(it + 1) % usePoints.size]
)
}
}
}
fun coordinateGenerate1(): List<List<Vector3D>> {
/**
* 绕 Z 轴旋转指定角度（弧度）
*/
fun Vector3D.rotateAroundZ(angle: Angle): Vector3D {
val cosAngle = cos(angle.radians)
val sinAngle = sin(angle.radians)
return Vector3D(
x = x * cosAngle - y * sinAngle,
y = x * sinAngle + y * cosAngle,
z = z
)
}
val center = Vector3D()
val direction = Vector3D(0.0, 1.0, -1.0)
return (0..360).step(10).map {
val nowDirection = direction.rotateAroundZ(it.degrees)
listOf(2, 6, 10).map {
center + nowDirection * it
}
}
}
class PolygonTest(
private val mapView: MapView
) {
private val geoHelper = GeoHelper.getSharedInstance()
private val contourSourceId = "contour-source-id-0"
private val contourLayerId = "contour-layer-id-0"
private val fillSourceId = "fill-source-id-0"
private val fillLayerId = "fill-layer-id-0"
fun Vector3D.toMapboxPoint(): Point {
return geoHelper.enuToWGS84Object(GeoHelper.ENU(x, y, z)).run {
Point.fromLngLat(lon, lat, hgt)
}
}
fun update(
outer: List<Vector3D>,
inner: List<Vector3D>,
other: List<Vector3D>
) {
val lineFeatures = mutableListOf<Feature>()
val fillFeatures = mutableListOf<Feature>()
val outerPoints = outer.map { it.toMapboxPoint() }
val innerPoints = inner.map { it.toMapboxPoint() }
val otherPoints = other.map { it.toMapboxPoint() }
val outerLine = LineString.fromLngLats(outerPoints)
Feature.fromGeometry(outerLine).also {
lineFeatures.add(it)
}
val innerLine = LineString.fromLngLats(innerPoints)
Feature.fromGeometry(innerLine).also {
lineFeatures.add(it)
}
//val polygon = Polygon.fromOuterInner(outerLine, innerLine)
val polygon = Polygon.fromLngLats(listOf(outerPoints, otherPoints, innerPoints))
mapView.mapboxMap.getStyle { style ->
setupLineLayer(
style = style,
sourceId = contourSourceId,
layerId = contourLayerId,
features = lineFeatures
)
setupFillLayer(
style = style,
sourceId = fillSourceId,
layerId = fillLayerId,
features = listOf(Feature.fromGeometry(polygon))
)
}
}
private fun setupLineLayer(
style: Style,
sourceId: String,
layerId: String,
features: List<Feature>
) {
style.removeStyleLayer(layerId)
style.removeStyleSource(sourceId)
val source = geoJsonSource(sourceId) {
featureCollection(FeatureCollection.fromFeatures(features))
}
style.addSource(source)
val layer = lineLayer(layerId, sourceId) {
lineColor(Color.RED)
lineWidth(2.0)
lineCap(LineCap.ROUND)
lineJoin(LineJoin.ROUND)
lineOpacity(0.8)
}
style.addLayer(layer)
}
private fun setupFillLayer(
style: Style,
sourceId: String,
layerId: String,
features: List<Feature>
) {
style.removeStyleLayer(layerId)
style.removeStyleSource(sourceId)
val source = geoJsonSource(sourceId) {
featureCollection(FeatureCollection.fromFeatures(features))
}
style.addSource(source)
val layer = fillLayer(fillLayerId, fillSourceId) {
fillColor(Color.YELLOW)
fillOpacity(0.3)
fillAntialias(true)
}
style.addLayer(layer)
}
fun clear() {
}
}
class MainViewModel : ViewModel() {
private val geoHelper = GeoHelper.getSharedInstance()
private val openElevation: OpenElevationApi = OpenElevation(SharedHttpClient(SharedJson()))
private val _points = MutableStateFlow<List<Point>>(emptyList())
init {
_points.map {
openElevation.lookup(it.map { GeoPoint(it.longitude(), it.latitude(), it.altitude()) })
}.catch {
Log.e(TAG, "高程请求失败", it)
}.map {
it.map {
val enu =
geoHelper.wgs84ToENU(lon = it.longitude, lat = it.latitude, hgt = it.altitude)
Vector2(enu.x, enu.y)
}
}.onEach {
val triangulation = DelaunayTriangulation(it)
triangulation.triangles().map {
it.contour
}
}.launchIn(viewModelScope)
}
fun addPoint(point: Point) {
_points.update {
it.toMutableList().apply {
add(point)
}
}
}
}

59
android/src/main/java/com/icegps/orx/MainViewModel.kt Normal file
View File

@@ -0,0 +1,59 @@
package com.icegps.orx
import android.app.Application
import android.util.Log
import androidx.lifecycle.AndroidViewModel
import androidx.lifecycle.viewModelScope
import com.icegps.common.helper.GeoHelper
import com.icegps.math.geometry.Vector3D
import com.icegps.orx.ktx.toast
import com.icegps.shared.SharedHttpClient
import com.icegps.shared.SharedJson
import com.icegps.shared.api.OpenElevation
import com.icegps.shared.api.OpenElevationApi
import com.icegps.shared.ktx.TAG
import com.icegps.shared.model.GeoPoint
import com.mapbox.geojson.Point
import kotlinx.coroutines.flow.MutableStateFlow
import kotlinx.coroutines.flow.SharingStarted
import kotlinx.coroutines.flow.catch
import kotlinx.coroutines.flow.debounce
import kotlinx.coroutines.flow.filter
import kotlinx.coroutines.flow.map
import kotlinx.coroutines.flow.stateIn
import kotlinx.coroutines.flow.update
class MainViewModel(private val context: Application) : AndroidViewModel(context) {
private val geoHelper = GeoHelper.Companion.getSharedInstance()
private val openElevation: OpenElevationApi = OpenElevation(SharedHttpClient(SharedJson()))
private val _points = MutableStateFlow<List<Point>>(emptyList())
val points = _points.filter { it.size > 3 }.debounce(1000).map {
openElevation.lookup(it.map { GeoPoint(it.longitude(), it.latitude(), it.altitude()) })
}.catch {
Log.e(TAG, "高程请求失败", it)
context.toast("高程请求失败")
}.map {
it.map {
val enu = geoHelper.wgs84ToENU(lon = it.longitude, lat = it.latitude, hgt = it.altitude)
Vector3D(enu.x, enu.y, enu.z)
}
}.stateIn(
scope = viewModelScope,
started = SharingStarted.Companion.Eagerly,
initialValue = emptyList()
)
fun addPoint(point: Point) {
context.toast("${point.longitude()}, ${point.latitude()}")
_points.update {
it.toMutableList().apply {
add(point)
}
}
}
fun clearPoints() {
_points.value = emptyList()
}
}

123
android/src/main/java/com/icegps/orx/PolygonTest.kt Normal file
View File

@@ -0,0 +1,123 @@
package com.icegps.orx
import android.graphics.Color
import com.icegps.common.helper.GeoHelper
import com.icegps.math.geometry.Vector3D
import com.icegps.orx.ktx.toMapboxPoint
import com.mapbox.geojson.Feature
import com.mapbox.geojson.FeatureCollection
import com.mapbox.geojson.LineString
import com.mapbox.geojson.Polygon
import com.mapbox.maps.MapView
import com.mapbox.maps.Style
import com.mapbox.maps.extension.style.layers.addLayer
import com.mapbox.maps.extension.style.layers.generated.fillLayer
import com.mapbox.maps.extension.style.layers.generated.lineLayer
import com.mapbox.maps.extension.style.layers.properties.generated.LineCap
import com.mapbox.maps.extension.style.layers.properties.generated.LineJoin
import com.mapbox.maps.extension.style.sources.addSource
import com.mapbox.maps.extension.style.sources.generated.geoJsonSource
class PolygonTest(
private val mapView: MapView
) {
private val geoHelper = GeoHelper.Companion.getSharedInstance()
private val contourSourceId = "contour-source-id-0"
private val contourLayerId = "contour-layer-id-0"
private val fillSourceId = "fill-source-id-0"
private val fillLayerId = "fill-layer-id-0"
fun update(
outer: List<Vector3D>,
inner: List<Vector3D>,
other: List<Vector3D>
) {
val lineFeatures = mutableListOf<Feature>()
val fillFeatures = mutableListOf<Feature>()
val outerPoints = outer.map { it.toMapboxPoint() }
val innerPoints = inner.map { it.toMapboxPoint() }
val otherPoints = other.map { it.toMapboxPoint() }
val outerLine = LineString.fromLngLats(outerPoints)
Feature.fromGeometry(outerLine).also {
lineFeatures.add(it)
}
val innerLine = LineString.fromLngLats(innerPoints)
Feature.fromGeometry(innerLine).also {
lineFeatures.add(it)
}
Feature.fromGeometry(LineString.fromLngLats(otherPoints)).also {
lineFeatures.add(it)
}
//val polygon = Polygon.fromOuterInner(outerLine, innerLine)
val polygon = Polygon.fromLngLats(listOf(outerPoints, otherPoints, innerPoints))
mapView.mapboxMap.getStyle { style ->
if (false) setupLineLayer(
style = style,
sourceId = contourSourceId,
layerId = contourLayerId,
features = lineFeatures
)
setupFillLayer(
style = style,
sourceId = fillSourceId,
layerId = fillLayerId,
features = listOf(Feature.fromGeometry(polygon))
)
}
}
private fun setupLineLayer(
style: Style,
sourceId: String,
layerId: String,
features: List<Feature>
) {
style.removeStyleLayer(layerId)
style.removeStyleSource(sourceId)
val source = geoJsonSource(sourceId) {
featureCollection(FeatureCollection.fromFeatures(features))
}
style.addSource(source)
val layer = lineLayer(layerId, sourceId) {
lineColor(Color.RED)
lineWidth(2.0)
lineCap(LineCap.Companion.ROUND)
lineJoin(LineJoin.Companion.ROUND)
lineOpacity(0.8)
}
style.addLayer(layer)
}
private fun setupFillLayer(
style: Style,
sourceId: String,
layerId: String,
features: List<Feature>
) {
style.removeStyleLayer(layerId)
style.removeStyleSource(sourceId)
val source = geoJsonSource(sourceId) {
featureCollection(FeatureCollection.fromFeatures(features))
}
style.addSource(source)
val layer = fillLayer(fillLayerId, fillSourceId) {
fillColor(Color.YELLOW)
fillOpacity(0.3)
fillAntialias(true)
}
style.addLayer(layer)
}
fun clear() {
}
}

123
android/src/main/java/com/icegps/orx/PolylineManager.kt Normal file
View File

@@ -0,0 +1,123 @@
package com.icegps.orx
import android.graphics.Color
import com.icegps.math.geometry.Line3D
import com.icegps.math.geometry.Vector3D
import com.icegps.orx.ktx.toMapboxPoint
import com.mapbox.geojson.Feature
import com.mapbox.geojson.FeatureCollection
import com.mapbox.geojson.LineString
import com.mapbox.maps.MapView
import com.mapbox.maps.Style
import com.mapbox.maps.extension.style.layers.addLayer
import com.mapbox.maps.extension.style.layers.generated.lineLayer
import com.mapbox.maps.extension.style.layers.properties.generated.LineCap
import com.mapbox.maps.extension.style.layers.properties.generated.LineJoin
import com.mapbox.maps.extension.style.sources.addSource
import com.mapbox.maps.extension.style.sources.generated.geoJsonSource
import org.openrndr.math.YPolarity
class PolylineManager(
private val mapView: MapView
) {
private val sourceId: String = "polyline-source-id-0"
private val layerId: String = "polyline-layer-id-0"
fun update(
points: List<List<Vector3D>>
) {
val lineStrings: List<List<Feature>> = points.map {
val lines = fromPoints(it, true)
lines.map {
LineString.fromLngLats(listOf(it.a.toMapboxPoint(), it.b.toMapboxPoint()))
}
}.map {
it.map { Feature.fromGeometry(it) }
}
mapView.mapboxMap.getStyle { style ->
setupLineLayer(
style = style,
sourceId = sourceId,
layerId = layerId,
features = lineStrings.flatten()
)
}
}
fun updateFeatures(
features: List<Feature>
) {
mapView.mapboxMap.getStyle { style ->
setupLineLayer(
style = style,
sourceId = sourceId,
layerId = layerId,
features = features
)
}
}
private fun setupLineLayer(
style: Style,
sourceId: String,
layerId: String,
features: List<Feature>
) {
style.removeStyleLayer(layerId)
style.removeStyleSource(sourceId)
val source = geoJsonSource(sourceId) {
featureCollection(FeatureCollection.fromFeatures(features))
}
style.addSource(source)
val layer = lineLayer(layerId, sourceId) {
lineColor(Color.RED)
lineWidth(2.0)
lineCap(LineCap.Companion.ROUND)
lineJoin(LineJoin.Companion.ROUND)
lineOpacity(0.8)
}
style.addLayer(layer)
}
fun clearContours() {
mapView.mapboxMap.getStyle { style ->
try {
style.removeStyleLayer(layerId)
} catch (_: Exception) {
}
try {
style.removeStyleSource(sourceId)
} catch (_: Exception) {
}
}
}
}
fun fromPoints(
points: List<Vector3D>,
closed: Boolean,
polarity: YPolarity = YPolarity.CW_NEGATIVE_Y
) = if (points.isEmpty()) {
emptyList()
} else {
if (!closed) {
(0 until points.size - 1).map {
Line3D(
points[it],
points[it + 1]
)
}
} else {
val d = (points.last() - points.first()).length
val usePoints = if (d > 1E-6) points else points.dropLast(1)
(usePoints.indices).map {
Line3D(
usePoints[it],
usePoints[(it + 1) % usePoints.size]
)
}
}
}

123
android/src/main/java/com/icegps/orx/SimplePalette.kt Normal file
View File

@@ -0,0 +1,123 @@
package com.icegps.orx
import android.util.Log
/**
* @author tabidachinokaze
* @date 2025/11/25
*/
class SimplePalette(
private var range: ClosedFloatingPointRange<Double>
) {
fun setRange(range: ClosedFloatingPointRange<Double>) {
this.range = range
}
private val colors: Map<Int, String>
init {
colors = generateTerrainColorMap()
}
fun palette(value: Double?): String {
if (value == null) return "#00000000"
val minH = range.start
val maxH = range.endInclusive
val normalized = ((value - minH) / (maxH - minH)).coerceIn(0.0, 1.0)
return colors[(normalized * 255).toInt()] ?: "#00000000"
}
fun palette1(value: Double?): String {
return if (value == null) "#00000000" else {
// 假设您已经知道高度范围，或者动态计算
val minH = range.start
val maxH = range.endInclusive
val normalized = ((value - minH) / (maxH - minH)).coerceIn(0.0, 1.0)
val alpha = (normalized * 255).toInt()
String.format("#%02X%02X%02X", alpha, 0, 0)
}.also {
Log.d("simplePalette", "$value -> $it")
}
}
fun generateGrayscaleColorMap2(): MutableMap<Int, String> {
val colorMap = mutableMapOf<Int, String>()
// 定义关键灰度点
val black = Color(0, 0, 0) // 低地势 - 黑色
val darkGray = Color(64, 64, 64) // 过渡
val midGray = Color(128, 128, 128) // 中间
val lightGray = Color(192, 192, 192) // 过渡
val white = Color(255, 255, 255) // 高地势 - 白色
for (i in 0..255) {
val position = i / 255.0
val color = when {
position < 0.25 -> interpolateColor(black, darkGray, position / 0.25)
position < 0.5 -> interpolateColor(darkGray, midGray, (position - 0.25) / 0.25)
position < 0.75 -> interpolateColor(midGray, lightGray, (position - 0.5) / 0.25)
else -> interpolateColor(lightGray, white, (position - 0.75) / 0.25)
}
colorMap[i] = color.toHex()
}
return colorMap
}
fun generateGrayscaleColorMap(): MutableMap<Int, String> {
val colorMap = mutableMapOf<Int, String>()
for (i in 0..255) {
// 从黑色到白色的线性渐变
val grayValue = i
val color = Color(grayValue, grayValue, grayValue)
colorMap[i] = color.toHex()
}
return colorMap
}
fun generateTerrainColorMap(): MutableMap<Int, String> {
val colorMap = mutableMapOf<Int, String>()
// 定义关键颜色点
val blue = Color(0, 0, 255) // 低地势 - 蓝色
val cyan = Color(0, 255, 255) // 中间过渡
val green = Color(0, 255, 0) // 中间过渡
val yellow = Color(255, 255, 0) // 中间过渡
val red = Color(255, 0, 0) // 高地势 - 红色
for (i in 0..255) {
val position = i / 255.0
val color = when {
position < 0.25 -> interpolateColor(blue, cyan, position / 0.25)
position < 0.5 -> interpolateColor(cyan, green, (position - 0.25) / 0.25)
position < 0.75 -> interpolateColor(green, yellow, (position - 0.5) / 0.25)
else -> interpolateColor(yellow, red, (position - 0.75) / 0.25)
}
colorMap[i] = color.toHex()
}
return colorMap
}
fun interpolateColor(start: Color, end: Color, fraction: Double): Color {
val r = (start.red + (end.red - start.red) * fraction).toInt()
val g = (start.green + (end.green - start.green) * fraction).toInt()
val b = (start.blue + (end.blue - start.blue) * fraction).toInt()
return Color(r, g, b)
}
// Color类简化实现
class Color(val red: Int, val green: Int, val blue: Int) {
fun toArgb(): Int {
return (0xFF shl 24) or (red shl 16) or (green shl 8) or blue
}
fun toHex(): String {
return String.format("#%06X", toArgb() and 0xFFFFFF)
}
}
}

23
android/src/main/java/com/icegps/orx/ktx/ColorRGBa.kt Normal file
View File

@@ -0,0 +1,23 @@
package com.icegps.orx.ktx
import org.openrndr.color.ColorRGBa
/**
* @author tabidachinokaze
* @date 2025/11/25
*/
fun ColorRGBa.toColorInt(): Int {
val clampedR = r.coerceIn(0.0, 1.0)
val clampedG = g.coerceIn(0.0, 1.0)
val clampedB = b.coerceIn(0.0, 1.0)
val clampedAlpha = alpha.coerceIn(0.0, 1.0)
return ((clampedAlpha * 255).toInt() shl 24) or
((clampedR * 255).toInt() shl 16) or
((clampedG * 255).toInt() shl 8) or
((clampedB * 255).toInt())
}
fun ColorRGBa.toColorHex(): String {
return String.format("#%06X", 0xFFFFFF and toColorInt())
}

12
android/src/main/java/com/icegps/orx/ktx/Context.kt Normal file
View File

@@ -0,0 +1,12 @@
package com.icegps.orx.ktx
import android.content.Context
import android.widget.Toast
/**
* @author tabidachinokaze
* @date 2025/11/25
*/
fun Context.toast(text: String, duration: Int = Toast.LENGTH_SHORT) {
Toast.makeText(this, text, duration).show()
}

22
android/src/main/java/com/icegps/orx/ktx/Vector2.kt Normal file
View File

@@ -0,0 +1,22 @@
package com.icegps.orx.ktx
import com.icegps.common.helper.GeoHelper
import com.mapbox.geojson.Point
import org.openrndr.math.Vector2
fun Vector2.niceStr(): String {
return "[$x, $y, 0.0]".format(this)
}
fun List<Vector2>.niceStr(): String {
return joinToString(", ", "[", "]") {
it.niceStr()
}
}
fun Vector2.toMapboxPoint(): Point {
val geoHelper = GeoHelper.getSharedInstance()
return geoHelper.enuToWGS84Object(GeoHelper.ENU(x, y)).run {
Point.fromLngLat(lon, lat, hgt)
}
}

22
android/src/main/java/com/icegps/orx/ktx/Vector3.kt Normal file
View File

@@ -0,0 +1,22 @@
package com.icegps.orx.ktx
import org.openrndr.math.Vector3
fun Vector3.niceStr(): String {
return "[$x, $y, $z]".format(this)
}
fun List<Vector3>.niceStr(): String {
return joinToString(", ", "[", "]") {
it.niceStr()
}
}
val List<Vector3>.area: org.openrndr.shape.Rectangle
get() {
val minX = minOf { it.x }
val maxX = maxOf { it.x }
val minY = minOf { it.y }
val maxY = maxOf { it.y }
return org.openrndr.shape.Rectangle(x = minX, y = minY, width = maxX - minX, height = maxY - minY)
}

32
android/src/main/java/com/icegps/orx/ktx/Vector3D.kt Normal file
View File

@@ -0,0 +1,32 @@
package com.icegps.orx.ktx
import com.icegps.common.helper.GeoHelper
import com.icegps.math.geometry.Rectangle
import com.icegps.math.geometry.Vector3D
import com.mapbox.geojson.Point
fun Vector3D.niceStr(): String {
return "[$x, $y, $z]".format(this)
}
fun List<Vector3D>.niceStr(): String {
return joinToString(", ", "[", "]") {
it.niceStr()
}
}
fun Vector3D.toMapboxPoint(): Point {
val geoHelper = GeoHelper.getSharedInstance()
return geoHelper.enuToWGS84Object(GeoHelper.ENU(x, y, z)).run {
Point.fromLngLat(lon, lat, hgt)
}
}
val List<Vector3D>.area: Rectangle
get() {
val minX = minOf { it.x }
val maxX = maxOf { it.x }
val minY = minOf { it.y }
val maxY = maxOf { it.y }
return Rectangle(x = minX, y = minY, width = maxX - minX, height = maxY - minY)
}

91
android/src/main/java/com/icegps/orx/triangulation/DelaunayTriangulation3D.kt Normal file
View File

@@ -0,0 +1,91 @@
package com.icegps.orx.triangulation
import org.openrndr.extra.triangulation.Delaunay
import org.openrndr.math.Vector3
import org.openrndr.shape.path3D
/**
* Kotlin/OPENRNDR idiomatic interface to `Delaunay`
*/
class DelaunayTriangulation3D(val points: List<Vector3>) {
val delaunay: Delaunay = Delaunay.Companion.from(points.map { it.xy })
fun neighbors(pointIndex: Int): Sequence<Int> {
return delaunay.neighbors(pointIndex)
}
fun neighborPoints(pointIndex: Int): List<Vector3> {
return neighbors(pointIndex).map { points[it] }.toList()
}
fun triangleIndices(): List<IntArray> {
val list = mutableListOf<IntArray>()
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 }): MutableList<Triangle3D> {
val list = mutableListOf<Triangle3D>()
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]
// originally they are defined *counterclockwise*
if (filterPredicate(t2, t1, t0)) {
val p1 = points[t0]
val p2 = points[t1]
val p3 = points[t2]
list.add(Triangle3D(p1, p2, p3))
}
}
return list
}
// Inner edges of the delaunay triangulation (without hull)
fun halfedges() = path3D {
for (i in delaunay.halfedges.indices) {
val j = delaunay.halfedges[i]
if (j < i) continue
val ti = delaunay.triangles[i]
val tj = delaunay.triangles[j]
moveTo(points[ti])
lineTo(points[tj])
}
}
fun hull() = path3D {
for (h in delaunay.hull) {
moveOrLineTo(points[h])
}
close()
}
fun nearest(query: Vector3): Int = delaunay.find(query.x, query.y)
fun nearestPoint(query: Vector3): Vector3 = points[nearest(query)]
}
/**
* Computes the Delaunay triangulation for the list of 2D points.
*
* The Delaunay triangulation is a triangulation of a set of points such that
* no point is inside the circumcircle of any triangle. It maximizes the minimum
* angle of all the angles in the triangles, avoiding skinny triangles.
*
* @return A DelaunayTriangulation object representing the triangulation of the given points.
*/
fun List<Vector3>.delaunayTriangulation(): DelaunayTriangulation3D {
return DelaunayTriangulation3D(this)
}

24
android/src/main/java/com/icegps/orx/triangulation/Triangle3D.kt Normal file
View File

@@ -0,0 +1,24 @@
package com.icegps.orx.triangulation
import org.openrndr.math.Vector3
import org.openrndr.shape.BezierSegment
import org.openrndr.shape.Path
import org.openrndr.shape.Path3D
/**
* @author tabidachinokaze
* @date 2025/11/24
*/
data class Triangle3D(
val x1: Vector3,
val x2: Vector3,
val x3: Vector3,
) : Path<Vector3> {
val path = Path3D.fromPoints(points = listOf(x1, x2, x3), closed = true)
override fun sub(t0: Double, t1: Double): Path<Vector3> = path.sub(t0, t1)
override val closed: Boolean get() = path.closed
override val empty: Boolean get() = path.empty
override val infinity: Vector3 get() = path.infinity
override val segments: List<BezierSegment<Vector3>> get() = path.segments
}

105
android/src/main/res/layout/activity_main.xml
View File

@@ -4,11 +4,110 @@
android:id="@+id/main"
android:layout_width="match_parent"
android:layout_height="match_parent"
android:orientation="vertical"
android:orientation="horizontal"
tools:context=".MainActivity">
<LinearLayout
android:layout_width="0dp"
android:layout_height="match_parent"
android:layout_weight="1"
android:orientation="vertical"
android:paddingTop="32dp">
<com.google.android.material.slider.Slider
android:id="@+id/slider_target_height"
android:layout_width="match_parent"
android:layout_height="wrap_content"
android:value="0"
android:valueFrom="0"
android:valueTo="100" />
<LinearLayout
android:layout_width="match_parent"
android:layout_height="wrap_content"
android:gravity="center_vertical"
android:orientation="horizontal">
<TextView
android:layout_width="wrap_content"
android:layout_height="wrap_content"
android:text="栅格大小：" />
<com.google.android.material.slider.Slider
android:id="@+id/cell_size"
android:layout_width="0dp"
android:layout_height="wrap_content"
android:layout_weight="1"
android:value="1"
android:valueFrom="1"
android:valueTo="100" />
</LinearLayout>
<LinearLayout
android:layout_width="match_parent"
android:layout_height="wrap_content"
android:gravity="center_vertical"
android:orientation="horizontal">
<TextView
android:layout_width="wrap_content"
android:layout_height="wrap_content"
android:text="高度范围：" />
<com.google.android.material.slider.RangeSlider
android:id="@+id/height_range"
android:layout_width="0dp"
android:layout_height="wrap_content"
android:layout_weight="1"
android:valueFrom="0"
android:valueTo="100" />
</LinearLayout>
<Switch
android:id="@+id/switch_grid"
android:layout_width="wrap_content"
android:layout_height="wrap_content"
android:switchPadding="16dp"
android:text="栅格网" />
<Switch
android:id="@+id/switch_triangle"
android:layout_width="wrap_content"
android:layout_height="wrap_content"
android:switchPadding="16dp"
android:text="三角网" />
<LinearLayout
android:layout_width="match_parent"
android:layout_height="wrap_content">
<TextView
android:layout_width="wrap_content"
android:layout_height="wrap_content"
android:text="坐标数量：" />
<TextView
android:id="@+id/point_count"
android:layout_width="wrap_content"
android:layout_height="wrap_content" />
</LinearLayout>
<Button
android:id="@+id/update"
android:layout_width="wrap_content"
android:layout_height="wrap_content"
android:text="刷新界面" />
<Button
android:id="@+id/clear_points"
android:layout_width="wrap_content"
android:layout_height="wrap_content"
android:text="清除所有点" />
</LinearLayout>
<com.mapbox.maps.MapView
android:id="@+id/map_view"
android:layout_width="match_parent"
android:layout_height="match_parent" />
android:layout_width="0dp"
android:layout_height="match_parent"
android:layout_weight="3" />
</LinearLayout>

1
desktop/build.gradle.kts
View File

@@ -30,6 +30,7 @@ kotlin {
implementation(project(":orx-marching-squares"))
implementation(project(":orx-text-writer"))
implementation(project(":orx-obj-loader"))
implementation(project(":orx-palette"))
}
}
}

35
desktop/src/jvmDemo/kotlin/DemoColorBrewer2.kt Normal file
View File

@@ -0,0 +1,35 @@
import org.openrndr.application
import org.openrndr.shape.Rectangle
/**
* Demonstrates how to use a ColorBrewer2 palette.
* Finds the first available palette with 5 colors,
* then draws concentric circles filled with those colors.
*/
fun main() = application {
configure {
width = 720
height = 720
}
program {
val palette = colorBrewer2Palettes(
numberOfColors = 6,
paletteType = ColorBrewer2Type.Any
).first().colors
val cellSize = 50.0
extend {
palette.forEachIndexed { i, color ->
drawer.fill = color
drawer.rectangle(
Rectangle(
x = 0.0,
y = cellSize * i,
width = cellSize,
height = cellSize
)
)
// drawer.circle(drawer.bounds.center, 300.0 - i * 40.0)
}
}
}
}

59
desktop/src/jvmDemo/kotlin/DemoDelaunay03.kt
View File

@@ -10,7 +10,6 @@ import org.openrndr.draw.loadFont
import org.openrndr.extra.camera.Camera2D
import org.openrndr.extra.marchingsquares.findContours
import org.openrndr.extra.noise.gradientPerturbFractal
import org.openrndr.extra.noise.simplex
import org.openrndr.extra.textwriter.writer
import org.openrndr.extra.triangulation.DelaunayTriangulation
import org.openrndr.math.Vector2
@@ -18,7 +17,6 @@ import org.openrndr.math.Vector3
import org.openrndr.shape.Segment2D
import org.openrndr.shape.Segment3D
import org.openrndr.shape.ShapeContour
import kotlin.math.absoluteValue
import kotlin.math.cos
import kotlin.math.sin
import kotlin.random.Random
@@ -30,7 +28,7 @@ import kotlin.random.Random
fun main() = application {
configure {
width = 720
height = 720
height = 480
title = "Delaunator"
}
program {
@@ -62,6 +60,10 @@ fun main() = application {
extend(Camera2D())
println("draw")
var targetHeight: Double = zs.average()
val step = zs.max() - zs.min() / 6
var heightList = (0..5).map { index ->
zs.min() + step * index
}
var logEnabled = true
var useInterpolation = false
var sampleLinear = false
@@ -171,12 +173,41 @@ fun main() = application {
cellSize = 4.0,
useInterpolation = useInterpolation
)
val associateWith: List<List<ShapeContour>> = heightList.map { height ->
findContours(
f = {
val triangle = triangles.firstOrNull { triangle ->
isPointInTriangle(it, listOf(triangle.x1, triangle.x2, triangle.x3))
}
triangle ?: return@findContours 0.0
val interpolate = interpolateHeight(
point = it,
triangle = listOf(
triangle.x1,
triangle.x2,
triangle.x3,
).map {
Vector3(it.x, it.y, associate[it]!!)
}
)
interpolate.z - height
},
area = drawer.bounds,
cellSize = 4.0,
useInterpolation = useInterpolation
)
}
if (logEnabled) println("useInterpolation = $useInterpolation")
drawer.stroke = null
contours.forEach {
if (true) contours.forEach {
drawer.fill = ColorRGBa.GREEN.opacify(0.1)
drawer.contour(if (sampleLinear) it.sampleLinear() else it)
}
if (false) associateWith.forEachIndexed { index, contours ->
contours.forEach {
drawer.fill = colorBrewer2[index].colors.first().opacify(0.1)
drawer.contour(it)
}
}
drawer.fontMap = loadFont("demo-data/fonts/IBMPlexMono-Regular.ttf", 24.0)
@@ -210,6 +241,24 @@ fun isPointInTriangle(point: Vector2, triangle: List<Vector2>): Boolean {
alpha <= 1 && beta <= 1 && gamma <= 1
}
fun isPointInTriangle3D(point: Vector2, triangle: List<Vector3>): Boolean {
require(triangle.size == 3) { "三角形必须有3个顶点" }
val (v1, v2, v3) = triangle
// 计算重心坐标
val denominator = (v2.y - v3.y) * (v1.x - v3.x) + (v3.x - v2.x) * (v1.y - v3.y)
if (denominator == 0.0) return false // 退化三角形
val alpha = ((v2.y - v3.y) * (point.x - v3.x) + (v3.x - v2.x) * (point.y - v3.y)) / denominator
val beta = ((v3.y - v1.y) * (point.x - v3.x) + (v1.x - v3.x) * (point.y - v3.y)) / denominator
val gamma = 1.0 - alpha - beta
// 点在三角形内当且仅当所有重心坐标都在[0,1]范围内
return alpha >= 0 && beta >= 0 && gamma >= 0 &&
alpha <= 1 && beta <= 1 && gamma <= 1
}
/**
* 使用重心坐标计算点在三角形上的高度
* @param point 二维点 (x, y)

40
desktop/src/jvmDemo/kotlin/DemoDelaunay3D.kt
View File

@@ -9,15 +9,14 @@ import org.openrndr.draw.loadFont
import org.openrndr.draw.shadeStyle
import org.openrndr.extra.camera.Orbital
import org.openrndr.extra.marchingsquares.findContours
import org.openrndr.extra.noise.gradientPerturbFractal
import org.openrndr.extra.objloader.loadOBJasVertexBuffer
import org.openrndr.extra.textwriter.writer
import org.openrndr.extra.triangulation.DelaunayTriangulation
import org.openrndr.extra.triangulation.DelaunayTriangulation3D
import org.openrndr.math.Vector2
import org.openrndr.math.Vector3
import org.openrndr.shape.Path3D
import org.openrndr.shape.Segment3D
import org.openrndr.shape.ShapeContour
/**
* @author tabidachinokaze
@@ -48,15 +47,16 @@ fun main() = application {
height = height,
volcanoCount = 3
)*/
val points3D = coordinateGenerate(width, height).map {
it.copy(x = it.x - width / 2, y = it.y - height / 2)
}
val points3D = coordinateGenerate(width, height).map {
it.copy(x = it.x - width / 2, y = it.y - height / 2, z = it.z * 10)
}
val zs = points3D.map { it.z }
println("zs = ${zs}")
val associate: MutableMap<Vector2, Double> = points3D.associate {
Vector2(it.x, it.y) to it.z
}.toMutableMap()
val delaunay = DelaunayTriangulation(associate.map { it.key })
val delaunay3D = DelaunayTriangulation3D(points3D.map { Vector3(it.x, it.y, it.z) })
//println(points3D.niceStr())
//extend(Camera2D())
@@ -84,7 +84,7 @@ fun main() = application {
val vb = loadOBJasVertexBuffer("orx-obj-loader/test-data/non-planar.obj")
extend {
val triangles = delaunay.triangles()
val triangles = delaunay3D.triangles()
val segments = mutableListOf<Segment3D>()
drawer.clear(ColorRGBa.BLACK)
val indexDiff = (frameCount / 1000) % triangles.size
@@ -95,10 +95,12 @@ fun main() = application {
}
drawer.vertexBuffer(vb, DrawPrimitive.TRIANGLES)
// 绘制等高线段区域
for ((i, triangle) in triangles.withIndex()) {
val segment2DS = triangle.contour.segments.filter {
val startZ = associate[it.start]!!
val endZ = associate[it.end]!!
val segment2DS = triangle.segments.filter {
val startZ = it.start.z
val endZ = it.end.z
if (startZ < endZ) {
targetHeight in startZ..endZ
} else {
@@ -108,8 +110,8 @@ fun main() = application {
if (segment2DS.size == 2) {
val vector2s = segment2DS.map {
val startZ = associate[it.start]!!
val endZ = associate[it.end]!!
val startZ = it.start.z
val endZ = it.end.z
val start = Vector3(it.start.x, it.start.y, startZ)
val end = Vector3(it.end.x, it.end.y, endZ)
if (startZ < endZ) {
@@ -130,14 +132,8 @@ fun main() = application {
}
drawer.strokeWeight = 20.0
drawer.stroke = ColorRGBa.PINK
val segment3DS = triangle.contour.segments.map {
val startZ = associate[it.start]!!
val endZ = associate[it.end]!!
Segment3D(it.start.vector3(z = startZ), it.end.vector3(z = endZ))
}
//drawer.contour(triangle.contour)
drawer.path(Path3D.fromSegments(segment3DS, closed = true))
drawer.path(triangle.path)
}
val sorted = connectAllSegments(segments)
@@ -161,6 +157,7 @@ fun main() = application {
drawer.fill = ColorRGBa.YELLOW
// if (false) drawer.contour(ShapeContour.fromSegments(it, closed = true))
}
// 结束绘制等高线
/*for (y in 0 until (area.height / cellSize).toInt()) {
for (x in 0 until (area.width / cellSize).toInt()) {
values[IntVector2(x, y)] = f(Vector2(x * cellSize + area.x, y * cellSize + area.y))
@@ -169,7 +166,7 @@ fun main() = application {
val contours = findContours(
f = {
val triangle = triangles.firstOrNull { triangle ->
isPointInTriangle(it, listOf(triangle.x1, triangle.x2, triangle.x3))
isPointInTriangle3D(it, listOf(triangle.x1, triangle.x2, triangle.x3))
}
triangle ?: return@findContours 0.0
val interpolate = interpolateHeight(
@@ -178,9 +175,7 @@ fun main() = application {
triangle.x1,
triangle.x2,
triangle.x3,
).map {
Vector3(it.x, it.y, associate[it]!!)
}
)
)
interpolate.z - targetHeight
},
@@ -191,6 +186,7 @@ fun main() = application {
if (logEnabled) println("useInterpolation = $useInterpolation")
drawer.stroke = null
contours.map {
if (false) drawer.contour(it)
it.segments.map {
Segment3D(
it.start.vector3(),

1
icegps-triangulation/.gitignore vendored Normal file
View File

@@ -0,0 +1 @@
/build

13
icegps-triangulation/build.gradle.kts Normal file
View File

@@ -0,0 +1,13 @@
plugins {
id("java-library")
alias(libs.plugins.kotlin.jvm)
}
java {
sourceCompatibility = JavaVersion.VERSION_11
targetCompatibility = JavaVersion.VERSION_11
}
kotlin {
compilerOptions {
jvmTarget = org.jetbrains.kotlin.gradle.dsl.JvmTarget.JVM_11
}
}

596
icegps-triangulation/src/main/java/com/icegps/triangulation/Delaunator.kt Normal file
View File

@@ -0,0 +1,596 @@
package org.openrndr.extra.triangulation
import kotlin.math.*
val EPSILON: Double = 2.0.pow(-52)
/**
* A Kotlin port of Mapbox's Delaunator incredibly fast JavaScript library for Delaunay triangulation of 2D points.
*
* @description Port of Mapbox's Delaunator (JavaScript) library - https://github.com/mapbox/delaunator
* @property coords flat positions' array - [x0, y0, x1, y1..]
*
* @since f0ed80d - commit
* @author Ricardo Matias
*/
@Suppress("unused")
class Delaunator(val coords: DoubleArray) {
val EDGE_STACK = IntArray(512)
private var count = coords.size shr 1
// arrays that will store the triangulation graph
val maxTriangles = (2 * count - 5).coerceAtLeast(0)
private val _triangles = IntArray(maxTriangles * 3)
private val _halfedges = IntArray(maxTriangles * 3)
lateinit var triangles: IntArray
lateinit var halfedges: IntArray
// temporary arrays for tracking the edges of the advancing convex hull
private var hashSize = ceil(sqrt(count * 1.0)).toInt()
private var hullPrev = IntArray(count) // edge to prev edge
private var hullNext = IntArray(count) // edge to next edge
private var hullTri = IntArray(count) // edge to adjacent triangle
private var hullHash = IntArray(hashSize) // angular edge hash
private var hullStart: Int = -1
// temporary arrays for sorting points
private var ids = IntArray(count)
private var dists = DoubleArray(count)
private var cx: Double = Double.NaN
private var cy: Double = Double.NaN
private var trianglesLen: Int = -1
lateinit var hull: IntArray
init {
update()
}
fun update() {
if (coords.size <= 2) {
halfedges = IntArray(0)
triangles = IntArray(0)
hull = IntArray(0)
return
}
// populate an array of point indices calculate input data bbox
var minX = Double.POSITIVE_INFINITY
var minY = Double.POSITIVE_INFINITY
var maxX = Double.NEGATIVE_INFINITY
var maxY = Double.NEGATIVE_INFINITY
// points -> points
// minX, minY, maxX, maxY
for (i in 0 until count) {
val x = coords[2 * i]
val y = coords[2 * i + 1]
if (x < minX) minX = x
if (y < minY) minY = y
if (x > maxX) maxX = x
if (y > maxY) maxY = y
ids[i] = i
}
val cx = (minX + maxX) / 2
val cy = (minY + maxY) / 2
var minDist = Double.POSITIVE_INFINITY
var i0: Int = -1
var i1: Int = -1
var i2: Int = -1
// pick a seed point close to the center
for (i in 0 until count) {
val d = dist(cx, cy, coords[2 * i], coords[2 * i + 1])
if (d < minDist) {
i0 = i
minDist = d
}
}
val i0x = coords[2 * i0]
val i0y = coords[2 * i0 + 1]
minDist = Double.POSITIVE_INFINITY
// Find the point closest to the seed
for(i in 0 until count) {
if (i == i0) continue
val d = dist(i0x, i0y, coords[2 * i], coords[2 * i + 1])
if (d < minDist && d > 0) {
i1 = i
minDist = d
}
}
var i1x = coords[2 * i1]
var i1y = coords[2 * i1 + 1]
var minRadius = Double.POSITIVE_INFINITY
// Find the third point which forms the smallest circumcircle with the first two
for (i in 0 until count) {
if(i == i0 || i == i1) continue
val r = circumradius(i0x, i0y, i1x, i1y, coords[2 * i], coords[2 * i + 1])
if(r < minRadius) {
i2 = i
minRadius = r
}
}
if (minRadius == Double.POSITIVE_INFINITY) {
// order collinear points by dx (or dy if all x are identical)
// and return the list as a hull
for (i in 0 until count) {
val a = (coords[2 * i] - coords[0])
val b = (coords[2 * i + 1] - coords[1])
dists[i] = if (a == 0.0) b else a
}
quicksort(ids, dists, 0, count - 1)
val nhull = IntArray(count)
var j = 0
var d0 = Double.NEGATIVE_INFINITY
for (i in 0 until count) {
val id = ids[i]
if (dists[id] > d0) {
nhull[j++] = id
d0 = dists[id]
}
}
hull = nhull.copyOf(j)
triangles = IntArray(0)
halfedges = IntArray(0)
return
}
var i2x = coords[2 * i2]
var i2y = coords[2 * i2 + 1]
// swap the order of the seed points for counter-clockwise orientation
if (orient2d(i0x, i0y, i1x, i1y, i2x, i2y) < 0.0) {
val i = i1
val x = i1x
val y = i1y
i1 = i2
i1x = i2x
i1y = i2y
i2 = i
i2x = x
i2y = y
}
val center = circumcenter(i0x, i0y, i1x, i1y, i2x, i2y)
this.cx = center[0]
this.cy = center[1]
for (i in 0 until count) {
dists[i] = dist(coords[2 * i], coords[2 * i + 1], center[0], center[1])
}
// sort the points by distance from the seed triangle circumcenter
quicksort(ids, dists, 0, count - 1)
// set up the seed triangle as the starting hull
hullStart = i0
var hullSize = 3
hullNext[i0] = i1
hullNext[i1] = i2
hullNext[i2] = i0
hullPrev[i2] = i1
hullPrev[i0] = i2
hullPrev[i1] = i0
hullTri[i0] = 0
hullTri[i1] = 1
hullTri[i2] = 2
hullHash.fill(-1)
hullHash[hashKey(i0x, i0y)] = i0
hullHash[hashKey(i1x, i1y)] = i1
hullHash[hashKey(i2x, i2y)] = i2
trianglesLen = 0
addTriangle(i0, i1, i2, -1, -1, -1)
var xp = 0.0
var yp = 0.0
for (k in ids.indices) {
val i = ids[k]
val x = coords[2 * i]
val y = coords[2 * i + 1]
// skip near-duplicate points
if (k > 0 && abs(x - xp) <= EPSILON && abs(y - yp) <= EPSILON) continue
xp = x
yp = y
// skip seed triangle points
if (i == i0 || i == i1 || i == i2) continue
// find a visible edge on the convex hull using edge hash
var start = 0
val key = hashKey(x, y)
for (j in 0 until hashSize) {
start = hullHash[(key + j) % hashSize]
if (start != -1 && start != hullNext[start]) break
}
start = hullPrev[start]
var e = start
var q = hullNext[e]
while (orient2d(x, y, coords[2 * e], coords[2 * e + 1], coords[2 * q], coords[2 * q + 1]) >= 0) {
e = q
if (e == start) {
e = -1
break
}
q = hullNext[e]
}
if (e == -1) continue // likely a near-duplicate point skip it
// add the first triangle from the point
var t = addTriangle(e, i, hullNext[e], -1, -1, hullTri[e])
// recursively flip triangles from the point until they satisfy the Delaunay condition
hullTri[i] = legalize(t + 2)
hullTri[e] = t // keep track of boundary triangles on the hull
hullSize++
// walk forward through the hull, adding more triangles and flipping recursively
var next = hullNext[e]
q = hullNext[next]
while (orient2d(x, y, coords[2 * next], coords[2 * next + 1], coords[2 * q], coords[2 * q + 1]) < 0) {
t = addTriangle(next, i, q, hullTri[i], -1, hullTri[next])
hullTri[i] = legalize(t + 2)
hullNext[next] = next // mark as removed
hullSize--
next = q
q = hullNext[next]
}
// walk backward from the other side, adding more triangles and flipping
if (e == start) {
q = hullPrev[e]
while (orient2d(x, y, coords[2 * q], coords[2 * q + 1], coords[2 * e], coords[2 * e + 1]) < 0) {
t = addTriangle(q, i, e, -1, hullTri[e], hullTri[q])
legalize(t + 2)
hullTri[q] = t
hullNext[e] = e // mark as removed
hullSize--
e = q
q = hullPrev[e]
}
}
// update the hull indices
hullStart = e
hullPrev[i] = e
hullNext[e] = i
hullPrev[next] = i
hullNext[i] = next
// save the two new edges in the hash table
hullHash[hashKey(x, y)] = i
hullHash[hashKey(coords[2 * e], coords[2 * e + 1])] = e
}
hull = IntArray(hullSize)
var e = hullStart
for (i in 0 until hullSize) {
hull[i] = e
e = hullNext[e]
}
// trim typed triangle mesh arrays
triangles = _triangles.copyOf(trianglesLen)
halfedges = _halfedges.copyOf(trianglesLen)
}
private fun legalize(a: Int): Int {
var i = 0
var na = a
var ar: Int
// recursion eliminated with a fixed-size stack
while (true) {
val b = _halfedges[na]
/* if the pair of triangles doesn't satisfy the Delaunay condition
* (p1 is inside the circumcircle of [p0, pl, pr]), flip them,
* then do the same check/flip recursively for the new pair of triangles
*
* pl pl
* /||\ / \
* al/ || \bl al/ \a
* / || \ / \
* / a||b \ flip /___ar___\
* p0\ || /p1 => p0\---bl---/p1
* \ || / \ /
* ar\ || /br b\ /br
* \||/ \ /
* pr pr
*/
val a0 = na - na % 3
ar = a0 + (na + 2) % 3
if (b == -1) { // convex hull edge
if (i == 0) break
na = EDGE_STACK[--i]
continue
}
val b0 = b - b % 3
val al = a0 + (na + 1) % 3
val bl = b0 + (b + 2) % 3
val p0 = _triangles[ar]
val pr = _triangles[na]
val pl = _triangles[al]
val p1 = _triangles[bl]
val illegal = inCircleRobust(
coords[2 * p0], coords[2 * p0 + 1],
coords[2 * pr], coords[2 * pr + 1],
coords[2 * pl], coords[2 * pl + 1],
coords[2 * p1], coords[2 * p1 + 1])
if (illegal) {
_triangles[na] = p1
_triangles[b] = p0
val hbl = _halfedges[bl]
// edge swapped on the other side of the hull (rare) fix the halfedge reference
if (hbl == -1) {
var e = hullStart
do {
if (hullTri[e] == bl) {
hullTri[e] = na
break
}
e = hullPrev[e]
} while (e != hullStart)
}
link(na, hbl)
link(b, _halfedges[ar])
link(ar, bl)
val br = b0 + (b + 1) % 3
// don't worry about hitting the cap: it can only happen on extremely degenerate input
if (i < EDGE_STACK.size) {
EDGE_STACK[i++] = br
}
} else {
if (i == 0) break
na = EDGE_STACK[--i]
}
}
return ar
}
private fun link(a:Int, b:Int) {
_halfedges[a] = b
if (b != -1) _halfedges[b] = a
}
// add a new triangle given vertex indices and adjacent half-edge ids
private fun addTriangle(i0: Int, i1: Int, i2: Int, a: Int, b: Int, c: Int): Int {
val t = trianglesLen
_triangles[t] = i0
_triangles[t + 1] = i1
_triangles[t + 2] = i2
link(t, a)
link(t + 1, b)
link(t + 2, c)
trianglesLen += 3
return t
}
private fun hashKey(x: Double, y: Double): Int {
return (floor(pseudoAngle(x - cx, y - cy) * hashSize) % hashSize).toInt()
}
}
fun circumradius(ax: Double, ay: Double,
bx: Double, by: Double,
cx: Double, cy: Double): Double {
val dx = bx - ax
val dy = by - ay
val ex = cx - ax
val ey = cy - ay
val bl = dx * dx + dy * dy
val cl = ex * ex + ey * ey
val d = 0.5 / (dx * ey - dy * ex)
val x = (ey * bl - dy * cl) * d
val y = (dx * cl - ex * bl) * d
return x * x + y * y
}
fun circumcenter(ax: Double, ay: Double,
bx: Double, by: Double,
cx: Double, cy: Double): DoubleArray {
val dx = bx - ax
val dy = by - ay
val ex = cx - ax
val ey = cy - ay
val bl = dx * dx + dy * dy
val cl = ex * ex + ey * ey
val d = 0.5 / (dx * ey - dy * ex)
val x = ax + (ey * bl - dy * cl) * d
val y = ay + (dx * cl - ex * bl) * d
return doubleArrayOf(x, y)
}
fun quicksort(ids: IntArray, dists: DoubleArray, left: Int, right: Int) {
if (right - left <= 20) {
for (i in (left + 1)..right) {
val temp = ids[i]
val tempDist = dists[temp]
var j = i - 1
while (j >= left && dists[ids[j]] > tempDist) ids[j + 1] = ids[j--]
ids[j + 1] = temp
}
} else {
val median = (left + right) shr 1
var i = left + 1
var j = right
swap(ids, median, i)
if (dists[ids[left]] > dists[ids[right]]) swap(ids, left, right)
if (dists[ids[i]] > dists[ids[right]]) swap(ids, i, right)
if (dists[ids[left]] > dists[ids[i]]) swap(ids, left, i)
val temp = ids[i]
val tempDist = dists[temp]
while (true) {
do i++ while (dists[ids[i]] < tempDist)
do j-- while (dists[ids[j]] > tempDist)
if (j < i) break
swap(ids, i, j)
}
ids[left + 1] = ids[j]
ids[j] = temp
if (right - i + 1 >= j - left) {
quicksort(ids, dists, i, right)
quicksort(ids, dists, left, j - 1)
} else {
quicksort(ids, dists, left, j - 1)
quicksort(ids, dists, i, right)
}
}
}
private fun swap(arr: IntArray, i: Int, j: Int) {
val tmp = arr[i]
arr[i] = arr[j]
arr[j] = tmp
}
// monotonically increases with real angle, but doesn't need expensive trigonometry
private fun pseudoAngle(dx: Double, dy: Double): Double {
val p = dx / (abs(dx) + abs(dy))
val a = if (dy > 0.0) 3.0 - p else 1.0 + p
return a / 4.0 // [0..1]
}
private fun inCircle(ax: Double, ay: Double,
bx: Double, by: Double,
cx: Double, cy: Double,
px: Double, py: Double): Boolean {
val dx = ax - px
val dy = ay - py
val ex = bx - px
val ey = by - py
val fx = cx - px
val fy = cy - py
val ap = dx * dx + dy * dy
val bp = ex * ex + ey * ey
val cp = fx * fx + fy * fy
return dx * (ey * cp - bp * fy) -
dy * (ex * cp - bp * fx) +
ap * (ex * fy - ey * fx) < 0
}
private fun inCircleRobust(
ax: Double, ay: Double,
bx: Double, by: Double,
cx: Double, cy: Double,
px: Double, py: Double
): Boolean {
val dx = twoDiff(ax, px)
val dy = twoDiff(ay, py)
val ex = twoDiff(bx, px)
val ey = twoDiff(by, py)
val fx = twoDiff(cx, px)
val fy = twoDiff(cy, py)
val ap = ddAddDd(ddMultDd(dx, dx), ddMultDd(dy, dy))
val bp = ddAddDd(ddMultDd(ex, ex), ddMultDd(ey, ey))
val cp = ddAddDd(ddMultDd(fx, fx), ddMultDd(fy, fy))
val dd = ddAddDd(
ddDiffDd(
ddMultDd(dx, ddDiffDd(ddMultDd(ey, cp), ddMultDd(bp, fy))),
ddMultDd(dy, ddDiffDd(ddMultDd(ex, cp), ddMultDd(bp, fx)))
),
ddMultDd(ap, ddDiffDd(ddMultDd(ex, fy), ddMultDd(ey, fx)))
)
return (dd[1]) <= 0
}
private fun dist(ax: Double, ay: Double, bx: Double, by: Double): Double {
//val dx = ax - bx
//val dy = ay - by
//return dx * dx + dy * dy
// double-double implementation but I think it is overkill.
val dx = twoDiff(ax, bx)
val dy = twoDiff(ay, by)
val dx2 = ddMultDd(dx, dx)
val dy2 = ddMultDd(dy, dy)
val d2 = ddAddDd(dx2, dy2)
return d2[0] + d2[1]
}

232
icegps-triangulation/src/main/java/com/icegps/triangulation/Delaunay.kt Normal file
View File

@@ -0,0 +1,232 @@
package org.openrndr.extra.triangulation
import org.openrndr.extra.triangulation.Delaunay.Companion.from
import org.openrndr.math.Vector2
import org.openrndr.shape.Rectangle
import kotlin.math.cos
import kotlin.math.pow
import kotlin.math.sin
/*
ISC License
Copyright 2021 Ricardo Matias.
Permission to use, copy, modify, and/or distribute this software for any purpose
with or without fee is hereby granted, provided that the above copyright notice
and this permission notice appear in all copies.
THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES WITH
REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND
FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY SPECIAL, DIRECT,
INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS
OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER
TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
THIS SOFTWARE.
*/
/**
* Use [from] static method to use the delaunay triangulation
*
* @description Port of d3-delaunay (JavaScript) library - https://github.com/d3/d3-delaunay
* @property points flat positions' array - [x0, y0, x1, y1..]
*
* @since 9258fa3 - commit
* @author Ricardo Matias
*/
@Suppress("unused")
class Delaunay(val points: DoubleArray) {
companion object {
/**
* Entry point for the delaunay triangulation
*
* @property points a list of 2D points
*/
fun from(points: List<Vector2>): Delaunay {
val n = points.size
val coords = DoubleArray(n * 2)
for (i in points.indices) {
val p = points[i]
coords[2 * i] = p.x
coords[2 * i + 1] = p.y
}
return Delaunay(coords)
}
}
private var delaunator: Delaunator = Delaunator(points)
val inedges = IntArray(points.size / 2)
private val hullIndex = IntArray(points.size / 2)
var halfedges: IntArray = delaunator.halfedges
var hull: IntArray = delaunator.hull
var triangles: IntArray = delaunator.triangles
init {
init()
}
fun update() {
delaunator.update()
init()
}
fun neighbors(i:Int) = sequence<Int> {
val e0 = inedges[i]
if (e0 != -1) {
var e = e0
var p0 = -1
loop@do {
p0 = triangles[e]
yield(p0)
e = if (e % 3 == 2) e - 2 else e + 1
if (e == -1) {
break@loop
}
if (triangles[e] != i) {
break@loop
//error("bad triangulation")
}
e = halfedges[e]
if (e == -1) {
val p = hull[(hullIndex[i] + 1) % hull.size]
if (p != p0) {
yield(p)
}
break@loop
}
} while (e != e0)
}
}
fun collinear(): Boolean {
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;
}
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 0 until points.size step 2) {
val p = jitter(points[i], points[i+1], r)
points[i] = p[0]
points[i+1] = p[1]
}
delaunator = Delaunator(points)
halfedges = delaunator.halfedges
hull = delaunator.hull
triangles = delaunator.triangles
}
inedges.fill(-1)
hullIndex.fill(-1)
// Compute an index from each point to an (arbitrary) incoming halfedge
// Used to give the first neighbor of each point for this reason,
// on the hull we give priority to exterior halfedges
for (e in halfedges.indices) {
val p = triangles[nextHalfedge(e)]
if (halfedges[e] == -1 || inedges[p] == -1) inedges[p] = e
}
for (i in hull.indices) {
hullIndex[hull[i]] = i
}
// degenerate case: 1 or 2 (distinct) points
if (hull.size in 1..2) {
triangles = IntArray(3) { -1 }
halfedges = IntArray(3) { -1 }
triangles[0] = hull[0]
inedges[hull[0]] = 1
if (hull.size == 2) {
inedges[hull[1]] = 0
triangles[1] = hull[1]
triangles[2] = hull[1]
}
}
}
fun find(x: Double, y: Double, i: Int = 0): Int {
var i1 = i
var c = step(i, x, y)
while (c >= 0 && c != i && c != i1) {
i1 = c
c = step(i1, x, y)
}
return c
}
fun nextHalfedge(e: Int) = if (e % 3 == 2) e - 2 else e + 1
fun prevHalfedge(e: Int) = if (e % 3 == 0) e + 2 else e - 1
fun step(i: Int, x: Double, y: Double): Int {
if (inedges[i] == -1 || points.isEmpty()) return (i + 1) % (points.size shr 1)
var c = i
var dc = (x - points[i * 2]).pow(2) + (y - points[i * 2 + 1]).pow(2)
val e0 = inedges[i]
var e = e0
do {
val t = triangles[e]
val dt = (x - points[t * 2]).pow(2) + (y - points[t * 2 + 1]).pow(2)
if (dt < dc) {
dc = dt
c = t
}
e = if (e % 3 == 2) e - 2 else e + 1
if (triangles[e] != i) {
//error("bad triangulation")
break
} // bad triangulation
e = halfedges[e]
if (e == -1) {
e = hull[(hullIndex[i] + 1) % hull.size]
if (e != t) {
if ((x - points[e * 2]).pow(2) + (y - points[e * 2 + 1]).pow(2) < dc) return e
}
break
}
} while (e != e0)
return c
}
/**
* Generates a Voronoi diagram based on the current Delaunay triangulation and the provided bounds.
*
* @param bounds A rectangle defining the boundaries within which the Voronoi diagram will be generated.
* @return A Voronoi instance representing the resulting Voronoi diagram.
*/
fun voronoi(bounds: Rectangle): Voronoi = Voronoi(this, bounds)
}

95
icegps-triangulation/src/main/java/com/icegps/triangulation/DelaunayTriangulation.kt Normal file
View File

@@ -0,0 +1,95 @@
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
/**
* Kotlin/OPENRNDR idiomatic interface to `Delaunay`
*/
class DelaunayTriangulation(val points: List<Vector2>) {
val delaunay: Delaunay = Delaunay.from(points)
fun voronoiDiagram(bounds: Rectangle) = VoronoiDiagram(this, bounds)
fun neighbors(pointIndex: Int): Sequence<Int> {
return delaunay.neighbors(pointIndex)
}
fun neighborPoints(pointIndex: Int): List<Vector2> {
return neighbors(pointIndex).map { points[it] }.toList()
}
fun triangleIndices(): List<IntArray> {
val list = mutableListOf<IntArray>()
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<Triangle> {
val list = mutableListOf<Triangle>()
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]
// originally they are defined *counterclockwise*
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
}
// Inner edges of the delaunay triangulation (without hull)
fun halfedges() = contours {
for (i in delaunay.halfedges.indices) {
val j = delaunay.halfedges[i]
if (j < i) continue
val ti = delaunay.triangles[i]
val tj = delaunay.triangles[j]
moveTo(points[ti])
lineTo(points[tj])
}
}
fun hull() = contour {
for (h in delaunay.hull) {
moveOrLineTo(points[h])
}
close()
}
fun nearest(query: Vector2): Int = delaunay.find(query.x, query.y)
fun nearestPoint(query: Vector2): Vector2 = points[nearest(query)]
}
/**
* Computes the Delaunay triangulation for the list of 2D points.
*
* The Delaunay triangulation is a triangulation of a set of points such that
* no point is inside the circumcircle of any triangle. It maximizes the minimum
* angle of all the angles in the triangles, avoiding skinny triangles.
*
* @return A DelaunayTriangulation object representing the triangulation of the given points.
*/
fun List<Vector2>.delaunayTriangulation(): DelaunayTriangulation {
return DelaunayTriangulation(this)
}

340
icegps-triangulation/src/main/java/com/icegps/triangulation/DoubleDouble.kt Normal file
View File

@@ -0,0 +1,340 @@
package org.openrndr.extra.triangulation
import kotlin.math.pow
// original code: https://github.com/FlorisSteenkamp/double-double/
/**
* Returns the difference and exact error of subtracting two floating point
* numbers.
* Uses an EFT (error-free transformation), i.e. `a-b === x+y` exactly.
* The returned result is a non-overlapping expansion (smallest value first!).
*
* * **precondition:** `abs(a) >= abs(b)` - A fast test that can be used is
* `(a > b) === (a > -b)`
*
* See https://people.eecs.berkeley.edu/~jrs/papers/robustr.pdf
*/
fun fastTwoDiff(a: Double, b: Double): DoubleArray {
val x = a - b;
val y = (a - x) - b;
return doubleArrayOf(y, x)
}
/**
* Returns the sum and exact error of adding two floating point numbers.
* Uses an EFT (error-free transformation), i.e. a+b === x+y exactly.
* The returned sum is a non-overlapping expansion (smallest value first!).
*
* Precondition: abs(a) >= abs(b) - A fast test that can be used is
* (a > b) === (a > -b)
*
* See https://people.eecs.berkeley.edu/~jrs/papers/robustr.pdf
*/
fun fastTwoSum(a: Double, b: Double): DoubleArray {
val x = a + b;
return doubleArrayOf(b - (x - a), x)
}
/**
* Truncates a floating point value's significand and returns the result.
* Similar to split, but with the ability to specify the number of bits to keep.
*
* **Theorem 17 (Veltkamp-Dekker)**: Let a be a p-bit floating-point number, where
* p >= 3. Choose a splitting point s such that p/2 <= s <= p-1. Then the
* following algorithm will produce a (p-s)-bit value a_hi and a
* nonoverlapping (s-1)-bit value a_lo such that abs(a_hi) >= abs(a_lo) and
* a = a_hi + a_lo.
*
* * see [Shewchuk](https://people.eecs.berkeley.edu/~jrs/papers/robustr.pdf)
*
* @param a a double
* @param bits the number of significand bits to leave intact
*/
fun reduceSignificand(
a: Double,
bits: Int
): Double {
val s = 53 - bits;
val f = 2.0.pow(s) + 1;
val c = f * a;
val r = c - (c - a);
return r;
}
/**
* === 2^Math.ceil(p/2) + 1 where p is the # of significand bits in a double === 53.
* @internal
*/
const val f = 134217729; // 2**27 + 1;
/**
* Returns the result of splitting a double into 2 26-bit doubles.
*
* Theorem 17 (Veltkamp-Dekker): Let a be a p-bit floating-point number, where
* p >= 3. Choose a splitting point s such that p/2 <= s <= p-1. Then the
* following algorithm will produce a (p-s)-bit value a_hi and a
* nonoverlapping (s-1)-bit value a_lo such that abs(a_hi) >= abs(a_lo) and
* a = a_hi + a_lo.
*
* see e.g. [Shewchuk](https://people.eecs.berkeley.edu/~jrs/papers/robustr.pdf)
* @param a A double floating point number
*/
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)
}
/**
* Returns the exact result of subtracting b from a.
*
* @param a minuend - a double-double precision floating point number
* @param b subtrahend - a double-double precision floating point number
*/
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)
}
/**
* Returns the exact result of multiplying two doubles.
*
* * the resulting array is the reverse of the standard twoSum in the literature.
*
* Theorem 18 (Shewchuk): Let a and b be p-bit floating-point numbers, where
* p >= 6. Then the following algorithm will produce a nonoverlapping expansion
* x + y such that ab = x + y, where x is an approximation to ab and y
* represents the roundoff error in the calculation of x. Furthermore, if
* round-to-even tiebreaking is used, x and y are non-adjacent.
*
* See https://people.eecs.berkeley.edu/~jrs/papers/robustr.pdf
* @param a A double
* @param b Another double
*/
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;
//const [bh, bl] = split(b);
val d = f * b;
val bh = d - (d - b);
val bl = b - bh;
val y = (al * bl) - ((x - (ah * bh)) - (al * bh) - (ah * bl));
//const err1 = x - (ah * bh);
//const err2 = err1 - (al * bh);
//const err3 = err2 - (ah * bl);
//const y = (al * bl) - err3;
return doubleArrayOf(y, x)
}
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 y = (al * al) - ((x - (ah * ah)) - 2 * (ah * al));
return doubleArrayOf(y, x)
}
/**
* Returns the exact result of adding two doubles.
*
* * the resulting array is the reverse of the standard twoSum in the literature.
*
* Theorem 7 (Knuth): Let a and b be p-bit floating-point numbers. Then the
* following algorithm will produce a nonoverlapping expansion x + y such that
* a + b = x + y, where x is an approximation to a + b and y is the roundoff
* error in the calculation of x.
*
* See https://people.eecs.berkeley.edu/~jrs/papers/robustr.pdf
*/
fun twoSum(a: Double, b: Double): DoubleArray {
val x = a + b;
val bv = x - a;
return doubleArrayOf((a - (x - bv)) + (b - bv), x)
}
/**
* Returns the result of subtracting the second given double-double-precision
* floating point number from the first.
*
* * relative error bound: 3u^2 + 13u^3, i.e. fl(a-b) = (a-b)(1+ϵ),
* where ϵ <= 3u^2 + 13u^3, u = 0.5 * Number.EPSILON
* * the error bound is not sharp - the worst case that could be found by the
* authors were 2.25u^2
*
* ALGORITHM 6 of https://hal.archives-ouvertes.fr/hal-01351529v3/document
* @param x a double-double precision floating point number
* @param y another double-double precision floating point number
*/
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);
//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;
//const [vl,vh] = fastTwoSum(sh,c)
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);
return doubleArrayOf(zl, zh)
}
/**
* Returns the product of two double-double-precision floating point numbers.
*
* * relative error bound: 7u^2, i.e. fl(a+b) = (a+b)(1+ϵ),
* where ϵ <= 7u^2, u = 0.5 * Number.EPSILON
* the error bound is not sharp - the worst case that could be found by the
* authors were 5u^2
*
* * ALGORITHM 10 of https://hal.archives-ouvertes.fr/hal-01351529v3/document
* @param x a double-double precision floating point number
* @param y another double-double precision floating point number
*/
fun ddMultDd(x: DoubleArray, y: DoubleArray): DoubleArray {
//const xl = x[0];
val xh = x[1];
//const yl = y[0];
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));
//return fastTwoSum(ch,cl1 + (xh*yl + xl*yh));
val b = cl1 + (xh * y[0] + x[0] * yh);
val xx = ch + b;
return doubleArrayOf(b - (xx - ch), xx)
}
/**
* Returns the result of adding two double-double-precision floating point
* numbers.
*
* * relative error bound: 3u^2 + 13u^3, i.e. fl(a+b) = (a+b)(1+ϵ),
* where ϵ <= 3u^2 + 13u^3, u = 0.5 * Number.EPSILON
* * the error bound is not sharp - the worst case that could be found by the
* authors were 2.25u^2
*
* ALGORITHM 6 of https://hal.archives-ouvertes.fr/hal-01351529v3/document
* @param x a double-double precision floating point number
* @param y another double-double precision floating point number
*/
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 [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 [vl,vh] = fastTwoSum(sh,c)
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);
return doubleArrayOf(zl, zh)
}
/**
* Returns the product of a double-double-precision floating point number and a
* double.
*
* * slower than ALGORITHM 8 (one call to fastTwoSum more) but about 2x more
* accurate
* * relative error bound: 1.5u^2 + 4u^3, i.e. fl(a+b) = (a+b)(1+ϵ),
* where ϵ <= 1.5u^2 + 4u^3, u = 0.5 * Number.EPSILON
* * the bound is very sharp
* * probably prefer `ddMultDouble2` due to extra speed
*
* * ALGORITHM 7 of https://hal.archives-ouvertes.fr/hal-01351529v3/document
* @param y a double
* @param x a double-double precision floating point number
*/
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 cl2 = xl * y;
//val [tl1,th] = fastTwoSum(ch,cl2);
val th = ch + cl2;
val tl1 = cl2 - (th - ch);
val tl2 = tl1 + cl1;
//val [zl,zh] = fastTwoSum(th,tl2);
val zh = th + tl2;
val zl = tl2 - (zh - th);
return doubleArrayOf(zl, zh);
}

19
icegps-triangulation/src/main/java/com/icegps/triangulation/Predicates.kt Normal file
View File

@@ -0,0 +1,19 @@
package org.openrndr.extra.triangulation
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.
/*
| a b | = | ax - cx ay - cy |
| c d | | bx - cx 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(ddMultDd(a, d), ddMultDd(b, c))
return determinant[1]
}

622
icegps-triangulation/src/main/java/com/icegps/triangulation/Voronoi.kt Normal file
View File

@@ -0,0 +1,622 @@
package org.openrndr.extra.triangulation
import org.openrndr.math.Vector2
import org.openrndr.shape.Rectangle
import kotlin.math.abs
import kotlin.math.floor
import kotlin.math.sign
/*
ISC License
Copyright 2021 Ricardo Matias.
Permission to use, copy, modify, and/or distribute this software for any purpose
with or without fee is hereby granted, provided that the above copyright notice
and this permission notice appear in all copies.
THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES WITH
REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND
FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY SPECIAL, DIRECT,
INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS
OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER
TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
THIS SOFTWARE.
*/
/**
* This is a fast library for computing the Voronoi diagram of a set of two-dimensional points.
* The Voronoi diagram is constructed by connecting the circumcenters of adjacent triangles
* in the Delaunay triangulation.
*
* @description Port of d3-delaunay (JavaScript) library - https://github.com/d3/d3-delaunay
* @property points flat positions' array - [x0, y0, x1, y1..]
*
* @since 9258fa3 - commit
* @author Ricardo Matias
*/
class Voronoi(val delaunay: Delaunay, val bounds: Rectangle) {
private val _circumcenters = DoubleArray(delaunay.points.size * 2)
lateinit var circumcenters: DoubleArray
private set
val vectors = DoubleArray(delaunay.points.size * 2)
init {
init()
}
fun update() {
delaunay.update()
init()
}
fun init() {
val points = delaunay.points
if (delaunay.points.isEmpty()) {
return
}
val triangles = delaunay.triangles
val hull = delaunay.hull
if (points.size == 2) {
_circumcenters[0] = points[0]
_circumcenters[1] = points[1]
circumcenters = _circumcenters
return
}
circumcenters = _circumcenters.copyOf(delaunay.triangles.size / 3 * 2)
// Compute circumcenters
var i = 0
var j = 0
var x: Double
var y: Double
while (i < triangles.size) {
val t1 = triangles[i] * 2
val t2 = triangles[i + 1] * 2
val t3 = triangles[i + 2] * 2
val x1 = points[t1]
val y1 = points[t1 + 1]
val x2 = points[t2]
val y2 = points[t2 + 1]
val x3 = points[t3]
val y3 = points[t3 + 1]
val dx = x2 - x1
val dy = y2 - y1
val ex = x3 - x1
val ey = y3 - y1
val ab = (dx * ey - dy * ex) * 2
if (abs(ab) < 1e-9) {
var a = 1e9
val r = triangles[0] * 2
a *= sign((points[r] - x1) * ey - (points[r + 1] - y1) * ex)
x = (x1 + x3) / 2 - a * ey
y = (y1 + y3) / 2 + a * ex
} else {
val d = 1 / ab
val bl = dx * dx + dy * dy
val cl = ex * ex + ey * ey
x = x1 + (ey * bl - dy * cl) * d
y = y1 + (dx * cl - ex * bl) * d
}
circumcenters[j] = x
circumcenters[j + 1] = y
i += 3
j += 2
}
// Compute exterior cell rays.
var h = hull[hull.size - 1]
var p0: Int
var p1 = h * 4
var x0: Double
var x1 = points[2 * h]
var y0: Double
var y1 = points[2 * h + 1]
var y01: Double
var x10: Double
vectors.fill(0.0)
for (idx in hull.indices) {
h = hull[idx]
p0 = p1
x0 = x1
y0 = y1
p1 = h * 4
x1 = points[2 * h]
y1 = points[2 * h + 1]
y01 = y0 - y1
x10 = x1 - x0
vectors[p0 + 2] = y01
vectors[p1] = y01
vectors[p0 + 3] = x10
vectors[p1 + 1] = x10
}
}
private fun cell(i: Int): MutableList<Double>? {
val inedges = delaunay.inedges
val halfedges = delaunay.halfedges
val triangles = delaunay.triangles
val e0 = inedges[i]
if (e0 == -1) return null // coincident point
val points = mutableListOf<Double>()
var e = e0
do {
val t = floor(e / 3.0).toInt()
points.add(circumcenters[t * 2])
points.add(circumcenters[t * 2 + 1])
e = if (e % 3 == 2) e - 2 else e + 1 // next half edge
if (triangles[e] != i) break
e = halfedges[e]
} while (e != e0 && e != -1)
return points
}
fun neighbors(i: Int) = sequence {
val ci = clip(i)
if (ci != null) {
for (j in delaunay.neighbors(i)) {
val cj = clip(j)
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) {
if (ci[ai] == cj[aj]
&& ci[ai + 1] == cj[aj + 1]
&& ci[(ai + 2) % li] == cj[(aj + lj - 2) % lj]
&& ci[(ai + 3) % li] == cj[(aj + lj - 1) % lj]
) {
yield(j)
break@loop
}
}
}
}
}
}
}
internal fun clip(i: Int): List<Double>? {
// degenerate case (1 valid point: return the box)
if (i == 0 && delaunay.points.size == 2) {
return listOf(
bounds.xmax,
bounds.ymin,
bounds.xmax,
bounds.ymax,
bounds.xmin,
bounds.ymax,
bounds.xmin,
bounds.ymin
)
}
val points = cell(i) ?: return null
val clipVectors = vectors
val v = i * 4
val a = !clipVectors[v].isFalsy()
val b = !clipVectors[v + 1].isFalsy()
return if (a || b) {
this.clipInfinite(i, points, clipVectors[v], clipVectors[v + 1], clipVectors[v + 2], clipVectors[v + 3])
} else {
this.clipFinite(i, points)
}
}
private fun clipInfinite(
i: Int,
points: MutableList<Double>,
vx0: Double,
vy0: Double,
vxn: Double,
vyn: Double
): List<Double>? {
var P: MutableList<Double>? = points.mutableCopyOf()
P!!
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!!)
var n = 0
if (P != null) {
n = P!!.size
var c0 = -1
var c1 = edgeCode(P[n - 2], P[n - 1])
var j = 0
var n = P.size
while (j < n) {
c0 = c1
c1 = edgeCode(P[j], P[j + 1])
if (c0 != 0 && c1 != 0) {
j = edge(i, c0, c1, P, j)
n = P.size
}
j += 2
}
} else if (this.contains(i, (bounds.xmin + bounds.xmax) / 2.0, (bounds.ymin + bounds.ymax) / 2.0)) {
P = mutableListOf(
bounds.xmin,
bounds.ymin,
bounds.xmax,
bounds.ymin,
bounds.xmax,
bounds.ymax,
bounds.xmin,
bounds.ymax
)
}
return P
}
private fun clipFinite(i: Int, points: MutableList<Double>): MutableList<Double>? {
val n = points.size
val P = mutableListOf<Double>()
var x0: Double
var y0: Double
var x1 = points[n - 2]
var y1 = points[n - 1]
var c0: Int
var c1: Int = regionCode(x1, y1)
var e0: Int? = null
var e1: Int? = 0
for (j in 0 until n step 2) {
x0 = x1
y0 = y1
x1 = points[j]
y1 = points[j + 1]
c0 = c1
c1 = regionCode(x1, y1)
if (c0 == 0 && c1 == 0) {
e0 = e1
e1 = 0
P.add(x1)
P.add(y1)
} else {
var S: DoubleArray?
var sx0: Double
var sy0: Double
var sx1: Double
var sy1: Double
if (c0 == 0) {
S = clipSegment(x0, y0, x1, y1, c0, c1)
if (S == null) continue
sx0 = S[0]
sy0 = S[1]
sx1 = S[2]
sy1 = S[3]
} else {
S = clipSegment(x1, y1, x0, y0, c1, c0)
if (S == null) continue
sx1 = S[0]
sy1 = S[1]
sx0 = S[2]
sy0 = S[3]
e0 = e1
e1 = this.edgeCode(sx0, sy0)
if (e0 != 0 && e1 != 0) this.edge(i, e0!!, e1, P, P.size)
P.add(sx0)
P.add(sy0)
}
e0 = e1
e1 = this.edgeCode(sx1, sy1);
if (e0.isTruthy() && e1.isTruthy()) this.edge(i, e0!!, e1, P, P.size);
P.add(sx1)
P.add(sy1)
}
}
if (P.isNotEmpty()) {
e0 = e1
e1 = this.edgeCode(P[0], P[1])
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,
bounds.ymin,
bounds.xmax,
bounds.ymax,
bounds.xmin,
bounds.ymax,
bounds.xmin,
bounds.ymin
)
} else {
return null
}
return P
}
private fun clipSegment(x0: Double, y0: Double, x1: Double, y1: Double, c0: Int, c1: Int): DoubleArray? {
var nx0: Double = x0
var ny0: Double = y0
var nx1: Double = x1
var ny1: Double = y1
var nc0: Int = c0
var nc1: Int = c1
while (true) {
if (nc0 == 0 && nc1 == 0) return doubleArrayOf(nx0, ny0, nx1, ny1)
// SHAKY STUFF
if ((nc0 and nc1) != 0) return null
var x: Double
var y: Double
val c: Int = if (nc0 != 0) nc0 else nc1
when {
(c and 0b1000) != 0 -> {
x = nx0 + (nx1 - nx0) * (bounds.ymax - ny0) / (ny1 - ny0)
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;
}
}
if (nc0 != 0) {
nx0 = x
ny0 = y
nc0 = this.regionCode(nx0, ny0)
} else {
nx1 = x
ny1 = y
nc1 = this.regionCode(nx1, ny1)
}
}
}
private fun regionCode(x: Double, y: Double): Int {
val xcode = when {
x < bounds.xmin -> 0b0001
x > bounds.xmax -> 0b0010
else -> 0b0000
}
val ycode = when {
y < bounds.ymin -> 0b0100
y > bounds.ymax -> 0b1000
else -> 0b0000
}
return xcode or ycode
}
private fun contains(i: Int, x: Double, y: Double): Boolean {
if (x.isNaN() || y.isNaN()) return false
return this.delaunay.step(i, x, y) == i;
}
private fun edge(i: Int, e0: Int, e1: Int, p: MutableList<Double>, j: Int): Int {
var j = j
var e = e0
loop@ while (e != e1) {
var x: Double = Double.NaN
var y: Double = Double.NaN
when (e) {
0b0101 -> { // top-left
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
y = bounds.ymin
}
}
if (((j < p.size && (p[j] != x)) || ((j + 1) < p.size && p[j + 1] != y)) && contains(i, x, y)) {
require(!x.isNaN())
require(!y.isNaN())
p.add(j, y)
p.add(j, x)
j += 2
} else if (j >= p.size && contains(i, x, y)) {
require(!x.isNaN())
require(!y.isNaN())
p.add(x)
p.add(y)
j += 2
}
}
if (p.size > 4) {
var idx = 0
var n = p.size
while (idx < n) {
val j = (idx + 2) % p.size
val k = (idx + 4) % p.size
if ((p[idx] == p[j] && p[j] == p[k])
|| (p[idx + 1] == p[j + 1] && p[j + 1] == p[k + 1])
) {
// SHAKY
p.removeAt(j)
p.removeAt(j)
idx -= 2
n -= 2
}
idx += 2
}
}
return j
}
private fun project(x0: Double, y0: Double, vx: Double, vy: Double): Vector2? {
var t = Double.POSITIVE_INFINITY
var c: Double
var x = Double.NaN
var y = Double.NaN
// top
if (vy < 0) {
if (y0 <= bounds.ymin) return null
c = (bounds.ymin - y0) / vy
if (c < t) {
t = c
y = bounds.ymin
x = x0 + t * vx
}
} else if (vy > 0) { // bottom
if (y0 >= bounds.ymax) return null
c = (bounds.ymax - y0) / vy
if (c < t) {
t = c
y = bounds.ymax
x = x0 + t * vx
}
}
// right
if (vx > 0) {
if (x0 >= bounds.xmax) return null
c = (bounds.xmax - x0) / vx
if (c < t) {
t = c
x = bounds.xmax
y = y0 + t * vy
}
} else if (vx < 0) { // left
if (x0 <= bounds.xmin) return null
c = (bounds.xmin - x0) / vx
if (c < t) {
t = c
x = bounds.xmin
y = y0 + t * vy
}
}
if (x.isNaN() || y.isNaN()) return null
return Vector2(x, y)
}
private fun edgeCode(x: Double, y: Double): Int {
val xcode = when (x) {
bounds.xmin -> 0b0001
bounds.xmax -> 0b0010
else -> 0b0000
}
val ycode = when (y) {
bounds.ymin -> 0b0100
bounds.ymax -> 0b1000
else -> 0b0000
}
return xcode or ycode
}
}
private fun Int?.isTruthy(): Boolean {
return (this != null && this != 0)
}
private fun <T> List<T>.mutableCopyOf(): MutableList<T> {
val original = this
return mutableListOf<T>().apply { addAll(original) }
}
private val Rectangle.xmin: Double
get() = this.corner.x
private val Rectangle.xmax: Double
get() = this.corner.x + width
private val Rectangle.ymin: Double
get() = this.corner.y
private val Rectangle.ymax: Double
get() = this.corner.y + height
private fun Double?.isFalsy() = this == null || this == -0.0 || this == 0.0 || isNaN()

7
icegps-triangulation/src/main/java/com/icegps/triangulation/typealias.kt Normal file
View File

@@ -0,0 +1,7 @@
package com.icegps.triangulation
/**
* @author tabidachinokaze
* @date 2025/11/24
*/
typealias Vector2 = Vector