add android and desktop modules

.gitignore

@@ -10,3 +10,4 @@ gradle.properties
 /ShaderError.glsl
 /.kotlin
 /.lwjgl
+/local.properties

android/.gitignore

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

android/build.gradle.kts

@@ -0,0 +1,62 @@
+import org.jetbrains.kotlin.gradle.dsl.JvmTarget
+
+plugins {
+    alias(libs.plugins.android.application)
+    alias(libs.plugins.kotlin.android)
+}
+
+android {
+    namespace = "com.icegps.orx"
+    compileSdk {
+        version = release(36)
+    }
+
+    defaultConfig {
+        applicationId = "com.icegps.orx"
+        minSdk = 28
+        targetSdk = 36
+        versionCode = 1
+        versionName = "1.0"
+
+        testInstrumentationRunner = "androidx.test.runner.AndroidJUnitRunner"
+    }
+
+    buildTypes {
+        release {
+            isMinifyEnabled = false
+            proguardFiles(
+                getDefaultProguardFile("proguard-android-optimize.txt"),
+                "proguard-rules.pro"
+            )
+        }
+    }
+    buildFeatures {
+        viewBinding = true
+    }
+    compileOptions {
+        sourceCompatibility = JavaVersion.VERSION_17
+        targetCompatibility = JavaVersion.VERSION_17
+    }
+}
+
+kotlin {
+    compilerOptions.jvmTarget = JvmTarget.JVM_17
+}
+
+dependencies {
+    implementation(libs.core.ktx)
+    implementation(libs.androidx.appcompat)
+    implementation(libs.material)
+    implementation(libs.androidx.activity)
+    implementation(libs.androidx.constraintlayout)
+    implementation(libs.mapbox.maps)
+    implementation(project(":math"))
+    implementation(project(":orx-triangulation"))
+    implementation(libs.androidx.lifecycle.runtime.ktx)
+    implementation(project(":icegps-common"))
+    implementation(project(":icegps-shared"))
+
+    testImplementation(libs.junit)
+    androidTestImplementation(libs.ext.junit)
+    androidTestImplementation(libs.androidx.espresso.core)
+}

android/proguard-rules.pro

@@ -0,0 +1,21 @@
+# Add project specific ProGuard rules here.
+# You can control the set of applied configuration files using the
+# proguardFiles setting in build.gradle.
+#
+# For more details, see
+#   http://developer.android.com/guide/developing/tools/proguard.html
+
+# If your project uses WebView with JS, uncomment the following
+# and specify the fully qualified class name to the JavaScript interface
+# class:
+#-keepclassmembers class fqcn.of.javascript.interface.for.webview {
+#   public *;
+#}
+
+# Uncomment this to preserve the line number information for
+# debugging stack traces.
+#-keepattributes SourceFile,LineNumberTable
+
+# If you keep the line number information, uncomment this to
+# hide the original source file name.
+#-renamesourcefileattribute SourceFile

android/src/androidTest/java/com/icegps/orx/ExampleInstrumentedTest.kt

@@ -0,0 +1,24 @@
+package com.icegps.orx
+
+import androidx.test.platform.app.InstrumentationRegistry
+import androidx.test.ext.junit.runners.AndroidJUnit4
+
+import org.junit.Test
+import org.junit.runner.RunWith
+
+import org.junit.Assert.*
+
+/**
+ * Instrumented test, which will execute on an Android device.
+ *
+ * See [testing documentation](http://d.android.com/tools/testing).
+ */
+@RunWith(AndroidJUnit4::class)
+class ExampleInstrumentedTest {
+    @Test
+    fun useAppContext() {
+        // Context of the app under test.
+        val appContext = InstrumentationRegistry.getInstrumentation().targetContext
+        assertEquals("com.icegps.orx", appContext.packageName)
+    }
+}

android/src/main/AndroidManifest.xml

@@ -0,0 +1,25 @@
+<?xml version="1.0" encoding="utf-8"?>
+<manifest xmlns:android="http://schemas.android.com/apk/res/android">
+    <!-- Include this permission to grab user's general location -->
+    <uses-permission android:name="android.permission.ACCESS_COARSE_LOCATION" />
+    <!-- Include only if your app benefits from precise location access. -->
+    <uses-permission android:name="android.permission.ACCESS_FINE_LOCATION" />
+    <application
+        android:allowBackup="true"
+        android:icon="@mipmap/ic_launcher"
+        android:label="@string/app_name"
+        android:roundIcon="@mipmap/ic_launcher_round"
+        android:supportsRtl="true"
+        android:theme="@style/Theme.Orx">
+        <activity
+            android:name=".MainActivity"
+            android:exported="true">
+            <intent-filter>
+                <action android:name="android.intent.action.MAIN" />
+
+                <category android:name="android.intent.category.LAUNCHER" />
+            </intent-filter>
+        </activity>
+    </application>
+
+</manifest>

android/src/main/java/com/icegps/orx/MainActivity.kt

@@ -0,0 +1,77 @@
+package com.icegps.orx
+
+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 com.icegps.common.helper.GeoHelper
+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.Point
+import kotlinx.coroutines.flow.MutableStateFlow
+import kotlinx.coroutines.flow.catch
+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.Vector3
+
+class MainActivity : AppCompatActivity() {
+    private lateinit var binding: ActivityMainBinding
+    override fun onCreate(savedInstanceState: Bundle?) {
+        super.onCreate(savedInstanceState)
+        enableEdgeToEdge()
+        binding = ActivityMainBinding.inflate(layoutInflater)
+        setContentView(binding.root)
+        ViewCompat.setOnApplyWindowInsetsListener(findViewById(R.id.main)) { v, insets ->
+            val systemBars = insets.getInsets(WindowInsetsCompat.Type.systemBars())
+            v.setPadding(systemBars.left, systemBars.top, systemBars.right, systemBars.bottom)
+            insets
+        }
+    }
+}
+
+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)
+            }
+        }
+    }
+}

android/src/main/res/drawable/ic_launcher_background.xml

@@ -0,0 +1,170 @@
+<?xml version="1.0" encoding="utf-8"?>
+<vector xmlns:android="http://schemas.android.com/apk/res/android"
+    android:width="108dp"
+    android:height="108dp"
+    android:viewportWidth="108"
+    android:viewportHeight="108">
+    <path
+        android:fillColor="#3DDC84"
+        android:pathData="M0,0h108v108h-108z" />
+    <path
+        android:fillColor="#00000000"
+        android:pathData="M9,0L9,108"
+        android:strokeWidth="0.8"
+        android:strokeColor="#33FFFFFF" />
+    <path
+        android:fillColor="#00000000"
+        android:pathData="M19,0L19,108"
+        android:strokeWidth="0.8"
+        android:strokeColor="#33FFFFFF" />
+    <path
+        android:fillColor="#00000000"
+        android:pathData="M29,0L29,108"
+        android:strokeWidth="0.8"
+        android:strokeColor="#33FFFFFF" />
+    <path
+        android:fillColor="#00000000"
+        android:pathData="M39,0L39,108"
+        android:strokeWidth="0.8"
+        android:strokeColor="#33FFFFFF" />
+    <path
+        android:fillColor="#00000000"
+        android:pathData="M49,0L49,108"
+        android:strokeWidth="0.8"
+        android:strokeColor="#33FFFFFF" />
+    <path
+        android:fillColor="#00000000"
+        android:pathData="M59,0L59,108"
+        android:strokeWidth="0.8"
+        android:strokeColor="#33FFFFFF" />
+    <path
+        android:fillColor="#00000000"
+        android:pathData="M69,0L69,108"
+        android:strokeWidth="0.8"
+        android:strokeColor="#33FFFFFF" />
+    <path
+        android:fillColor="#00000000"
+        android:pathData="M79,0L79,108"
+        android:strokeWidth="0.8"
+        android:strokeColor="#33FFFFFF" />
+    <path
+        android:fillColor="#00000000"
+        android:pathData="M89,0L89,108"
+        android:strokeWidth="0.8"
+        android:strokeColor="#33FFFFFF" />
+    <path
+        android:fillColor="#00000000"
+        android:pathData="M99,0L99,108"
+        android:strokeWidth="0.8"
+        android:strokeColor="#33FFFFFF" />
+    <path
+        android:fillColor="#00000000"
+        android:pathData="M0,9L108,9"
+        android:strokeWidth="0.8"
+        android:strokeColor="#33FFFFFF" />
+    <path
+        android:fillColor="#00000000"
+        android:pathData="M0,19L108,19"
+        android:strokeWidth="0.8"
+        android:strokeColor="#33FFFFFF" />
+    <path
+        android:fillColor="#00000000"
+        android:pathData="M0,29L108,29"
+        android:strokeWidth="0.8"
+        android:strokeColor="#33FFFFFF" />
+    <path
+        android:fillColor="#00000000"
+        android:pathData="M0,39L108,39"
+        android:strokeWidth="0.8"
+        android:strokeColor="#33FFFFFF" />
+    <path
+        android:fillColor="#00000000"
+        android:pathData="M0,49L108,49"
+        android:strokeWidth="0.8"
+        android:strokeColor="#33FFFFFF" />
+    <path
+        android:fillColor="#00000000"
+        android:pathData="M0,59L108,59"
+        android:strokeWidth="0.8"
+        android:strokeColor="#33FFFFFF" />
+    <path
+        android:fillColor="#00000000"
+        android:pathData="M0,69L108,69"
+        android:strokeWidth="0.8"
+        android:strokeColor="#33FFFFFF" />
+    <path
+        android:fillColor="#00000000"
+        android:pathData="M0,79L108,79"
+        android:strokeWidth="0.8"
+        android:strokeColor="#33FFFFFF" />
+    <path
+        android:fillColor="#00000000"
+        android:pathData="M0,89L108,89"
+        android:strokeWidth="0.8"
+        android:strokeColor="#33FFFFFF" />
+    <path
+        android:fillColor="#00000000"
+        android:pathData="M0,99L108,99"
+        android:strokeWidth="0.8"
+        android:strokeColor="#33FFFFFF" />
+    <path
+        android:fillColor="#00000000"
+        android:pathData="M19,29L89,29"
+        android:strokeWidth="0.8"
+        android:strokeColor="#33FFFFFF" />
+    <path
+        android:fillColor="#00000000"
+        android:pathData="M19,39L89,39"
+        android:strokeWidth="0.8"
+        android:strokeColor="#33FFFFFF" />
+    <path
+        android:fillColor="#00000000"
+        android:pathData="M19,49L89,49"
+        android:strokeWidth="0.8"
+        android:strokeColor="#33FFFFFF" />
+    <path
+        android:fillColor="#00000000"
+        android:pathData="M19,59L89,59"
+        android:strokeWidth="0.8"
+        android:strokeColor="#33FFFFFF" />
+    <path
+        android:fillColor="#00000000"
+        android:pathData="M19,69L89,69"
+        android:strokeWidth="0.8"
+        android:strokeColor="#33FFFFFF" />
+    <path
+        android:fillColor="#00000000"
+        android:pathData="M19,79L89,79"
+        android:strokeWidth="0.8"
+        android:strokeColor="#33FFFFFF" />
+    <path
+        android:fillColor="#00000000"
+        android:pathData="M29,19L29,89"
+        android:strokeWidth="0.8"
+        android:strokeColor="#33FFFFFF" />
+    <path
+        android:fillColor="#00000000"
+        android:pathData="M39,19L39,89"
+        android:strokeWidth="0.8"
+        android:strokeColor="#33FFFFFF" />
+    <path
+        android:fillColor="#00000000"
+        android:pathData="M49,19L49,89"
+        android:strokeWidth="0.8"
+        android:strokeColor="#33FFFFFF" />
+    <path
+        android:fillColor="#00000000"
+        android:pathData="M59,19L59,89"
+        android:strokeWidth="0.8"
+        android:strokeColor="#33FFFFFF" />
+    <path
+        android:fillColor="#00000000"
+        android:pathData="M69,19L69,89"
+        android:strokeWidth="0.8"
+        android:strokeColor="#33FFFFFF" />
+    <path
+        android:fillColor="#00000000"
+        android:pathData="M79,19L79,89"
+        android:strokeWidth="0.8"
+        android:strokeColor="#33FFFFFF" />
+</vector>

android/src/main/res/drawable/ic_launcher_foreground.xml

@@ -0,0 +1,30 @@
+<vector xmlns:android="http://schemas.android.com/apk/res/android"
+    xmlns:aapt="http://schemas.android.com/aapt"
+    android:width="108dp"
+    android:height="108dp"
+    android:viewportWidth="108"
+    android:viewportHeight="108">
+    <path android:pathData="M31,63.928c0,0 6.4,-11 12.1,-13.1c7.2,-2.6 26,-1.4 26,-1.4l38.1,38.1L107,108.928l-32,-1L31,63.928z">
+        <aapt:attr name="android:fillColor">
+            <gradient
+                android:endX="85.84757"
+                android:endY="92.4963"
+                android:startX="42.9492"
+                android:startY="49.59793"
+                android:type="linear">
+                <item
+                    android:color="#44000000"
+                    android:offset="0.0" />
+                <item
+                    android:color="#00000000"
+                    android:offset="1.0" />
+            </gradient>
+        </aapt:attr>
+    </path>
+    <path
+        android:fillColor="#FFFFFF"
+        android:fillType="nonZero"
+        android:pathData="M65.3,45.828l3.8,-6.6c0.2,-0.4 0.1,-0.9 -0.3,-1.1c-0.4,-0.2 -0.9,-0.1 -1.1,0.3l-3.9,6.7c-6.3,-2.8 -13.4,-2.8 -19.7,0l-3.9,-6.7c-0.2,-0.4 -0.7,-0.5 -1.1,-0.3C38.8,38.328 38.7,38.828 38.9,39.228l3.8,6.6C36.2,49.428 31.7,56.028 31,63.928h46C76.3,56.028 71.8,49.428 65.3,45.828zM43.4,57.328c-0.8,0 -1.5,-0.5 -1.8,-1.2c-0.3,-0.7 -0.1,-1.5 0.4,-2.1c0.5,-0.5 1.4,-0.7 2.1,-0.4c0.7,0.3 1.2,1 1.2,1.8C45.3,56.528 44.5,57.328 43.4,57.328L43.4,57.328zM64.6,57.328c-0.8,0 -1.5,-0.5 -1.8,-1.2s-0.1,-1.5 0.4,-2.1c0.5,-0.5 1.4,-0.7 2.1,-0.4c0.7,0.3 1.2,1 1.2,1.8C66.5,56.528 65.6,57.328 64.6,57.328L64.6,57.328z"
+        android:strokeWidth="1"
+        android:strokeColor="#00000000" />
+</vector>

android/src/main/res/layout/activity_main.xml

@@ -0,0 +1,19 @@
+<?xml version="1.0" encoding="utf-8"?>
+<androidx.constraintlayout.widget.ConstraintLayout xmlns:android="http://schemas.android.com/apk/res/android"
+    xmlns:app="http://schemas.android.com/apk/res-auto"
+    xmlns:tools="http://schemas.android.com/tools"
+    android:id="@+id/main"
+    android:layout_width="match_parent"
+    android:layout_height="match_parent"
+    tools:context=".MainActivity">
+
+    <TextView
+        android:layout_width="wrap_content"
+        android:layout_height="wrap_content"
+        android:text="Hello World!"
+        app:layout_constraintBottom_toBottomOf="parent"
+        app:layout_constraintEnd_toEndOf="parent"
+        app:layout_constraintStart_toStartOf="parent"
+        app:layout_constraintTop_toTopOf="parent" />
+
+</androidx.constraintlayout.widget.ConstraintLayout>

android/src/main/res/mipmap-anydpi/ic_launcher.xml

@@ -0,0 +1,6 @@
+<?xml version="1.0" encoding="utf-8"?>
+<adaptive-icon xmlns:android="http://schemas.android.com/apk/res/android">
+    <background android:drawable="@drawable/ic_launcher_background" />
+    <foreground android:drawable="@drawable/ic_launcher_foreground" />
+    <monochrome android:drawable="@drawable/ic_launcher_foreground" />
+</adaptive-icon>

android/src/main/res/mipmap-anydpi/ic_launcher_round.xml

@@ -0,0 +1,6 @@
+<?xml version="1.0" encoding="utf-8"?>
+<adaptive-icon xmlns:android="http://schemas.android.com/apk/res/android">
+    <background android:drawable="@drawable/ic_launcher_background" />
+    <foreground android:drawable="@drawable/ic_launcher_foreground" />
+    <monochrome android:drawable="@drawable/ic_launcher_foreground" />
+</adaptive-icon>

android/src/main/res/values-night/themes.xml

@@ -0,0 +1,7 @@
+<resources xmlns:tools="http://schemas.android.com/tools">
+    <!-- Base application theme. -->
+    <style name="Base.Theme.Orx" parent="Theme.Material3.DayNight.NoActionBar">
+        <!-- Customize your dark theme here. -->
+        <!-- <item name="colorPrimary">@color/my_dark_primary</item> -->
+    </style>
+</resources>

android/src/main/res/values/colors.xml

@@ -0,0 +1,5 @@
+<?xml version="1.0" encoding="utf-8"?>
+<resources>
+    <color name="black">#FF000000</color>
+    <color name="white">#FFFFFFFF</color>
+</resources>

android/src/main/res/values/mapbox_access_token.xml

@@ -0,0 +1,4 @@
+<?xml version="1.0" encoding="utf-8"?>
+<resources xmlns:tools="http://schemas.android.com/tools">
+    <string name="mapbox_access_token" translatable="false" tools:ignore="UnusedResources">pk.eyJ1IjoienpxMSIsImEiOiJjbWYzbzV1MzQwMHJvMmpvbG1wbjJwdjUyIn0.LvKjIrCv9dAFcGxOM52f2Q</string>
+</resources>

android/src/main/res/values/strings.xml

@@ -0,0 +1,3 @@
+<resources>
+    <string name="app_name">icegps-orx</string>
+</resources>

android/src/main/res/values/themes.xml

@@ -0,0 +1,9 @@
+<resources xmlns:tools="http://schemas.android.com/tools">
+    <!-- Base application theme. -->
+    <style name="Base.Theme.Orx" parent="Theme.Material3.DayNight.NoActionBar">
+        <!-- Customize your light theme here. -->
+        <!-- <item name="colorPrimary">@color/my_light_primary</item> -->
+    </style>
+
+    <style name="Theme.Orx" parent="Base.Theme.Orx" />
+</resources>

android/src/test/java/com/icegps/orx/ExampleUnitTest.kt

@@ -0,0 +1,17 @@
+package com.icegps.orx
+
+import org.junit.Test
+
+import org.junit.Assert.*
+
+/**
+ * Example local unit test, which will execute on the development machine (host).
+ *
+ * See [testing documentation](http://d.android.com/tools/testing).
+ */
+class ExampleUnitTest {
+    @Test
+    fun addition_isCorrect() {
+        assertEquals(4, 2 + 2)
+    }
+}

build.gradle

@@ -2,6 +2,10 @@ plugins {
     alias(libs.plugins.nebula.release)
     alias(libs.plugins.nmcp)
     id("org.openrndr.extra.convention.dokka")
+    alias(libs.plugins.android.application) apply false
+    alias(libs.plugins.kotlin.android) apply false
+    alias(libs.plugins.kotlin.jvm) apply false
+    alias(libs.plugins.android.library) apply false
 }
 
 repositories {

desktop/build.gradle.kts

@@ -0,0 +1,36 @@
+plugins {
+    id("org.openrndr.extra.convention.kotlin-multiplatform")
+}
+
+kotlin {
+    sourceSets {
+        val commonMain by getting {
+            dependencies {
+                api(openrndr.math)
+                api(openrndr.shape)
+                implementation(project(":orx-noise"))
+            }
+        }
+        val commonTest by getting {
+            dependencies {
+                implementation(project(":orx-shapes"))
+                implementation(openrndr.shape)
+            }
+        }
+
+        val jvmDemo by getting {
+            dependencies {
+                implementation(project(":orx-triangulation"))
+                implementation(project(":orx-shapes"))
+                implementation(project(":orx-noise"))
+                implementation(openrndr.shape)
+                implementation(project(":math"))
+                implementation(project(":orx-camera"))
+                implementation("org.jetbrains.kotlinx:kotlinx-coroutines-core:1.10.2")
+                implementation(project(":orx-marching-squares"))
+                implementation(project(":orx-text-writer"))
+                implementation(project(":orx-obj-loader"))
+            }
+        }
+    }
+}

desktop/src/jvmDemo/kotlin/DemoDelaunay03.kt

@@ -0,0 +1,605 @@
+import com.icegps.math.geometry.Angle
+import com.icegps.math.geometry.Vector3D
+import com.icegps.math.geometry.degrees
+import org.openrndr.KEY_ARROW_DOWN
+import org.openrndr.KEY_ARROW_UP
+import org.openrndr.application
+import org.openrndr.color.ColorRGBa
+import org.openrndr.draw.TextSettingMode
+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
+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
+
+/**
+ * @author tabidachinokaze
+ * @date 2025/11/22
+ */
+fun main() = application {
+    configure {
+        width = 720
+        height = 720
+        title = "Delaunator"
+    }
+    program {
+        val points3D = (0 until height).step(36).map { y ->
+            (0 until width).step(36).map { x ->
+                gradientPerturbFractal(
+                    300,
+                    frequency = 0.8,
+                    position = Vector3(x.toDouble(), y.toDouble(), seconds)
+                )
+            }
+        }.flatten().map {
+            it.copy(z = it.z * 100)
+        }
+        /*val points3D = HeightmapVolcanoGenerator.generateVolcanoClusterHeightmap(
+            width = width,
+            height = height,
+            volcanoCount = 3
+        )*/
+        // val points3D = coordinateGenerate(width, height)
+        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 })
+
+        //println(points3D.niceStr())
+        extend(Camera2D())
+        println("draw")
+        var targetHeight: Double = zs.average()
+        var logEnabled = true
+        var useInterpolation = false
+        var sampleLinear = false
+        keyboard.keyDown.listen {
+            logEnabled = true
+            println(it)
+            when (it.key) {
+                KEY_ARROW_UP -> targetHeight++
+                KEY_ARROW_DOWN -> targetHeight--
+                73 -> useInterpolation = !useInterpolation
+                83 -> sampleLinear = !sampleLinear
+            }
+        }
+        extend {
+            val triangles = delaunay.triangles()
+            val segments = mutableListOf<Segment2D>()
+            drawer.clear(ColorRGBa.BLACK)
+            val indexDiff = (frameCount / 1000) % triangles.size
+            for ((i, triangle) in triangles.withIndex()) {
+                val segment2DS = triangle.contour.segments.filter {
+                    val startZ = associate[it.start]!!
+                    val endZ = associate[it.end]!!
+                    if (startZ < endZ) {
+                        targetHeight in startZ..endZ
+                    } else {
+                        targetHeight in endZ..startZ
+                    }
+                }
+
+                if (segment2DS.size == 2) {
+                    val vector2s = segment2DS.map {
+                        val startZ = associate[it.start]!!
+                        val endZ = associate[it.end]!!
+                        val start = Vector3(it.start.x, it.start.y, startZ)
+                        val end = Vector3(it.end.x, it.end.y, endZ)
+                        if (startZ < endZ) {
+                            start to end
+                        } else {
+                            end to start
+                        }
+                    }.map { (start, end) ->
+                        val segment3D = Segment3D(start, end)
+                        val vector3 =
+                            segment3D.position(calculatePositionRatio(targetHeight, start.z, end.z))
+                        vector3
+                    }.map {
+                        associate[it.xy] = it.z
+                        it.xy
+                    }
+                    val element = Segment2D(vector2s[0], vector2s[1])
+                    segments.add(element)
+                }
+                drawer.fill = if (indexDiff == i) {
+                    ColorRGBa.CYAN
+                } else {
+                    ColorRGBa.PINK.shade(1.0 - i / (triangles.size * 1.2))
+                }
+                drawer.stroke = ColorRGBa.PINK.shade(i / (triangles.size * 1.0) + 0.1)
+                drawer.contour(triangle.contour)
+            }
+
+            val sorted = connectAllSegments(segments)
+
+            drawer.stroke = ColorRGBa.WHITE
+            drawer.strokeWeight = 2.0
+            if (logEnabled) {
+                segments.forEach {
+                    println("${it.start} -> ${it.end}")
+                }
+                println("=====")
+            }
+
+            sorted.forEach {
+                it.forEach {
+                    if (logEnabled) println("${it.start} -> ${it.end}")
+                    drawer.lineSegment(it.start, it.end)
+                    drawer.fill = ColorRGBa.WHITE
+                }
+                if (logEnabled) println("=")
+                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))
+                }
+            }*/
+            val contours = 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 - targetHeight
+                },
+                area = drawer.bounds,
+                cellSize = 4.0,
+                useInterpolation = useInterpolation
+            )
+            if (logEnabled) println("useInterpolation = $useInterpolation")
+            drawer.stroke = null
+            contours.forEach {
+                drawer.fill = ColorRGBa.GREEN.opacify(0.1)
+                drawer.contour(if (sampleLinear) it.sampleLinear() else it)
+
+            }
+
+            drawer.fontMap = loadFont("demo-data/fonts/IBMPlexMono-Regular.ttf", 24.0)
+            writer {
+                drawer.drawStyle.textSetting = TextSettingMode.SUBPIXEL
+                text(targetHeight.toString())
+            }
+            logEnabled = false
+        }
+    }
+}
+
+/**
+ * 射线法判断点是否在单个三角形内
+ */
+fun isPointInTriangle(point: Vector2, triangle: List<Vector2>): 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 {
+    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)
+}
+
+/**
+ * 计算点在三角形中的重心坐标
+ */
+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)
+}
+
+fun connectAllSegments(segments: List<Segment2D>): List<List<Segment2D>> {
+    val remaining = segments.toMutableList()
+    val allPaths = mutableListOf<List<Segment2D>>()
+
+    while (remaining.isNotEmpty()) {
+        val path = mutableListOf<Segment2D>()
+
+        // 开始新路径
+        path.add(remaining.removeAt(0))
+
+        var changed: Boolean
+        do {
+            changed = false
+
+            // 向前扩展
+            val lastEnd = path.last().end
+            val forwardSegment = remaining.find { it.start == lastEnd || it.end == lastEnd }
+            if (forwardSegment != null) {
+                val connectedSegment = if (forwardSegment.start == lastEnd) {
+                    forwardSegment // 正向
+                } else {
+                    Segment2D(forwardSegment.end, forwardSegment.start) // 反向
+                }
+                path.add(connectedSegment)
+                remaining.remove(forwardSegment)
+                changed = true
+            }
+
+            // 向后扩展
+            val firstStart = path.first().start
+            val backwardSegment = remaining.find { it.end == firstStart || it.start == firstStart }
+            if (backwardSegment != null) {
+                val connectedSegment = if (backwardSegment.end == firstStart) {
+                    backwardSegment // 正向
+                } else {
+                    Segment2D(backwardSegment.end, backwardSegment.start) // 反向
+                }
+                path.add(0, connectedSegment)
+                remaining.remove(backwardSegment)
+                changed = true
+            }
+
+        } while (changed && remaining.isNotEmpty())
+
+        allPaths.add(path)
+    }
+
+    return allPaths
+}
+
+fun connectSegmentsEfficient(segments: List<Segment2D>): List<Segment2D> {
+    if (segments.isEmpty()) return emptyList()
+
+    val remaining = segments.toMutableList()
+    val connected = mutableListOf<Segment2D>()
+
+    // 构建端点查找表
+    val startMap = mutableMapOf<Vector2, MutableList<Segment2D>>()
+    val endMap = mutableMapOf<Vector2, MutableList<Segment2D>>()
+
+    segments.forEach { segment ->
+        startMap.getOrPut(segment.start) { mutableListOf() }.add(segment)
+        endMap.getOrPut(segment.end) { mutableListOf() }.add(segment)
+    }
+
+    // 从第一个线段开始
+    var currentSegment = remaining.removeAt(0)
+    connected.add(currentSegment)
+
+    // 更新查找表
+    startMap[currentSegment.start]?.remove(currentSegment)
+    endMap[currentSegment.end]?.remove(currentSegment)
+
+    // 向前连接
+    while (true) {
+        val nextFromStart = startMap[currentSegment.end]?.firstOrNull()
+        val nextFromEnd = endMap[currentSegment.end]?.firstOrNull()
+
+        when {
+            nextFromStart != null -> {
+                // 正向连接
+                connected.add(nextFromStart)
+                remaining.remove(nextFromStart)
+                startMap[nextFromStart.start]?.remove(nextFromStart)
+                endMap[nextFromStart.end]?.remove(nextFromStart)
+                currentSegment = nextFromStart
+            }
+
+            nextFromEnd != null -> {
+                // 反向连接
+                val reversed = Segment2D(nextFromEnd.end, nextFromEnd.start)
+                connected.add(reversed)
+                remaining.remove(nextFromEnd)
+                startMap[nextFromEnd.start]?.remove(nextFromEnd)
+                endMap[nextFromEnd.end]?.remove(nextFromEnd)
+                currentSegment = reversed
+            }
+
+            else -> break
+        }
+    }
+
+    // 向后连接
+    currentSegment = connected.first()
+    while (true) {
+        val prevFromEnd = endMap[currentSegment.start]?.firstOrNull()
+        val prevFromStart = startMap[currentSegment.start]?.firstOrNull()
+
+        when {
+            prevFromEnd != null -> {
+                // 正向连接到开头
+                connected.add(0, prevFromEnd)
+                remaining.remove(prevFromEnd)
+                startMap[prevFromEnd.start]?.remove(prevFromEnd)
+                endMap[prevFromEnd.end]?.remove(prevFromEnd)
+                currentSegment = prevFromEnd
+            }
+
+            prevFromStart != null -> {
+                // 反向连接到开头
+                val reversed = Segment2D(prevFromStart.end, prevFromStart.start)
+                connected.add(0, reversed)
+                remaining.remove(prevFromStart)
+                startMap[prevFromStart.start]?.remove(prevFromStart)
+                endMap[prevFromStart.end]?.remove(prevFromStart)
+                currentSegment = reversed
+            }
+
+            else -> break
+        }
+    }
+
+    return connected
+}
+
+fun connectSegments(segments: List<Segment2D>): List<Segment2D> {
+    if (segments.isEmpty()) return emptyList()
+
+    val remaining = segments.toMutableList()
+    val connected = mutableListOf<Segment2D>()
+
+    // 从第一个线段开始，保持原方向
+    connected.add(remaining.removeAt(0))
+
+    while (remaining.isNotEmpty()) {
+        val lastEnd = connected.last().end
+        var found = false
+
+        // 查找可以连接的线段
+        for (i in remaining.indices) {
+            val segment = remaining[i]
+
+            // 检查四种可能的连接方式
+            when {
+                // 正向连接：当前终点 == 线段起点
+                segment.start == lastEnd -> {
+                    connected.add(segment)
+                    remaining.removeAt(i)
+                    found = true
+                    break
+                }
+                // 反向连接：当前终点 == 线段终点，需要反转线段
+                segment.end == lastEnd -> {
+                    connected.add(Segment2D(segment.end, segment.start)) // 反转
+                    remaining.removeAt(i)
+                    found = true
+                    break
+                }
+                // 正向连接另一端：当前起点 == 线段终点，需要插入到前面
+                segment.end == connected.first().start -> {
+                    connected.add(0, Segment2D(segment.end, segment.start)) // 反转后插入开头
+                    remaining.removeAt(i)
+                    found = true
+                    break
+                }
+                // 反向连接另一端：当前起点 == 线段起点，需要反转并插入到前面
+                segment.start == connected.first().start -> {
+                    connected.add(0, segment) // 直接插入开头（已经是正确方向）
+                    remaining.removeAt(i)
+                    found = true
+                    break
+                }
+            }
+        }
+
+        if (!found) break // 无法找到连接线段
+    }
+
+    return connected
+}
+
+fun calculatePositionRatio(value: Double, rangeStart: Double, rangeEnd: Double): Double {
+    if (rangeStart == rangeEnd) return 0.0 // 避免除零
+
+    val ratio = (value - rangeStart) / (rangeEnd - rangeStart)
+    return ratio.coerceIn(0.0, 1.0)
+}
+
+fun sortLinesEfficient(lines: List<Segment2D>): List<Segment2D> {
+    if (lines.isEmpty()) return emptyList()
+
+    // 创建起点到线段的映射
+    val startMap = lines.associateBy { it.start }
+    val sorted = mutableListOf<Segment2D>()
+
+    // 找到起点（没有其他线段的终点指向它的起点）
+    var currentLine = lines.firstOrNull { line ->
+        lines.none { it.end == line.start }
+    } ?: lines.first()
+
+    sorted.add(currentLine)
+
+    while (true) {
+        val nextLine = startMap[currentLine.end]
+        if (nextLine == null || nextLine == lines.first()) break
+        sorted.add(nextLine)
+        currentLine = nextLine
+    }
+
+    return sorted
+}
+
+fun sortLines(lines: List<Segment2D>): List<Segment2D> {
+    if (lines.isEmpty()) return emptyList()
+
+    val remaining = lines.toMutableList()
+    val sorted = mutableListOf<Segment2D>()
+
+    // 从第一个线段开始
+    sorted.add(remaining.removeAt(0))
+
+    while (remaining.isNotEmpty()) {
+        val lastEnd = sorted.last().end
+        var found = false
+
+        // 查找下一个线段
+        for (i in remaining.indices) {
+            if (remaining[i].start == lastEnd) {
+                sorted.add(remaining.removeAt(i))
+                found = true
+                break
+            }
+        }
+
+        if (!found) break // 无法找到下一个线段
+    }
+
+    return sorted
+}
+
+fun findLineLoops(lines: List<Segment2D>): List<List<Segment2D>> {
+    val remaining = lines.toMutableList()
+    val loops = mutableListOf<List<Segment2D>>()
+
+    while (remaining.isNotEmpty()) {
+        val loop = findSingleLoop(remaining)
+        if (loop.isNotEmpty()) {
+            loops.add(loop)
+            // 移除已使用的线段
+            loop.forEach { line ->
+                remaining.remove(line)
+            }
+        } else {
+            // 无法形成环的线段
+            break
+        }
+    }
+
+    return loops
+}
+
+fun findSingleLoop(remaining: MutableList<Segment2D>): List<Segment2D> {
+    if (remaining.isEmpty()) return emptyList()
+
+    val loop = mutableListOf<Segment2D>()
+    loop.add(remaining.removeAt(0))
+
+    // 向前查找连接
+    while (remaining.isNotEmpty()) {
+        val lastEnd = loop.last().end
+        val nextIndex = remaining.indexOfFirst { it.start == lastEnd }
+
+        if (nextIndex == -1) {
+            // 尝试向后查找连接
+            val firstStart = loop.first().start
+            val prevIndex = remaining.indexOfFirst { it.end == firstStart }
+
+            if (prevIndex != -1) {
+                loop.add(0, remaining.removeAt(prevIndex))
+            } else {
+                break // 无法继续连接
+            }
+        } else {
+            loop.add(remaining.removeAt(nextIndex))
+        }
+
+        // 检查是否形成闭环
+        if (loop.last().end == loop.first().start) {
+            return loop
+        }
+    }
+
+    // 如果没有形成闭环，返回空列表（或者可以根据需求返回部分环）
+    remaining.addAll(loop) // 将线段放回剩余列表
+    return emptyList()
+}
+
+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
+    )
+}
+
+fun coordinateGenerate(width: Int, height: Int): List<Vector3D> {
+    val minX = 0.0
+    val maxX = width.toDouble()
+    val minY = 0.0
+    val maxY = height.toDouble()
+    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<Vector3D> {
+    val center = Vector3D(0.0, 0.0, 0.0)
+    val direction = Vector3D(0.0, 1.0, -1.0)
+    return (0..360).step(36).map<Int, List<Vector3D>> { degrees: Int ->
+        val newDirection = direction.rotateAroundZ(angle = degrees.degrees)
+        (0..5).map {
+            center + newDirection * it * 100
+        }
+    }.flatten()
+}
+
+
+fun Vector3D.niceStr(): String {
+    return "[$x, $y, $z]".format(this)
+}
+
+fun List<Vector3D>.niceStr(): String {
+    return joinToString(", ", "[", "]") {
+        it.niceStr()
+    }
+}

desktop/src/jvmDemo/kotlin/DemoDelaunay3D.kt

@@ -0,0 +1,271 @@
+import org.openrndr.KEY_ARROW_DOWN
+import org.openrndr.KEY_ARROW_UP
+import org.openrndr.WindowMultisample
+import org.openrndr.application
+import org.openrndr.color.ColorRGBa
+import org.openrndr.draw.DrawPrimitive
+import org.openrndr.draw.TextSettingMode
+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.math.Vector2
+import org.openrndr.math.Vector3
+import org.openrndr.shape.Path3D
+import org.openrndr.shape.Segment3D
+import org.openrndr.shape.ShapeContour
+
+/**
+ * @author tabidachinokaze
+ * @date 2025/11/22
+ */
+fun main() = application {
+    configure {
+        width = 720
+        height = 720
+        title = "Delaunator"
+        multisample = WindowMultisample.SampleCount(8)
+
+    }
+    program {
+        /*val points3D = (0 until height).step(36).map { y ->
+            (0 until width).step(36).map { x ->
+                gradientPerturbFractal(
+                    300,
+                    frequency = 0.8,
+                    position = Vector3(x.toDouble(), y.toDouble(), seconds)
+                )
+            }
+        }.flatten().map {
+            it.copy(x = it.x - width / 2, y = it.y - height / 2, z = it.z * 100)
+        }*/
+        /*val points3D = HeightmapVolcanoGenerator.generateVolcanoClusterHeightmap(
+            width = width,
+            height = height,
+            volcanoCount = 3
+        )*/
+         val points3D = coordinateGenerate(width, height).map {
+             it.copy(x = it.x - width / 2, y = it.y - height / 2)
+         }
+        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 })
+
+        //println(points3D.niceStr())
+        //extend(Camera2D())
+        val cam = Orbital()
+        extend(cam) {
+            eye = Vector3(x = 100.0, y = 100.0, z = 0.0)
+            lookAt = Vector3(x = 1.6, y = -1.9, z = 1.2)
+        }
+
+        println("draw")
+        var targetHeight: Double = zs.average()
+        var logEnabled = true
+        var useInterpolation = false
+        var sampleLinear = false
+        keyboard.keyDown.listen {
+            logEnabled = true
+            println(it)
+            when (it.key) {
+                KEY_ARROW_UP -> targetHeight++
+                KEY_ARROW_DOWN -> targetHeight--
+                73 -> useInterpolation = !useInterpolation
+                83 -> sampleLinear = !sampleLinear
+            }
+        }
+        val vb = loadOBJasVertexBuffer("orx-obj-loader/test-data/non-planar.obj")
+
+        extend {
+            val triangles = delaunay.triangles()
+            val segments = mutableListOf<Segment3D>()
+            drawer.clear(ColorRGBa.BLACK)
+            val indexDiff = (frameCount / 1000) % triangles.size
+            drawer.shadeStyle = shadeStyle {
+                fragmentTransform = """
+                        x_fill.rgb = normalize(v_viewNormal) * 0.5 + vec3(0.5);
+                    """.trimIndent()
+            }
+
+            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]!!
+                    if (startZ < endZ) {
+                        targetHeight in startZ..endZ
+                    } else {
+                        targetHeight in endZ..startZ
+                    }
+                }
+
+                if (segment2DS.size == 2) {
+                    val vector2s = segment2DS.map {
+                        val startZ = associate[it.start]!!
+                        val endZ = associate[it.end]!!
+                        val start = Vector3(it.start.x, it.start.y, startZ)
+                        val end = Vector3(it.end.x, it.end.y, endZ)
+                        if (startZ < endZ) {
+                            start to end
+                        } else {
+                            end to start
+                        }
+                    }.map { (start, end) ->
+                        val segment3D = Segment3D(start, end)
+                        val vector3 =
+                            segment3D.position(calculatePositionRatio(targetHeight, start.z, end.z))
+                        vector3
+                    }.onEach {
+                        associate[it.xy] = it.z
+                    }
+                    val element = Segment3D(vector2s[0], vector2s[1])
+                    segments.add(element)
+                }
+                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))
+            }
+
+            val sorted = connectAllSegments(segments)
+
+            drawer.stroke = ColorRGBa.WHITE
+            drawer.strokeWeight = 2.0
+            if (logEnabled) {
+                segments.forEach {
+                    println("${it.start} -> ${it.end}")
+                }
+                println("=====")
+            }
+
+            sorted.forEach {
+                it.forEach {
+                    if (logEnabled) println("${it.start} -> ${it.end}")
+                    drawer.lineSegment(it.start, it.end)
+                    drawer.fill = ColorRGBa.WHITE
+                }
+                if (logEnabled) println("=")
+                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))
+                }
+            }*/
+            val contours = 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 - targetHeight
+                },
+                area = drawer.bounds.movedTo(Vector2(-width / 2.0, -height / 2.0)),
+                cellSize = 4.0,
+                useInterpolation = useInterpolation
+            )
+            if (logEnabled) println("useInterpolation = $useInterpolation")
+            drawer.stroke = null
+            contours.map {
+                it.segments.map {
+                    Segment3D(
+                        it.start.vector3(),
+                        it.end.vector3()
+                    )
+                }
+            }.forEach {
+                drawer.fill = ColorRGBa.GREEN.opacify(0.1)
+                drawer.path(Path3D.fromSegments(it, closed = true))
+            }
+
+            if (false) writer {
+                drawer.fontMap = loadFont("demo-data/fonts/IBMPlexMono-Regular.ttf", 24.0)
+                drawer.drawStyle.textSetting = TextSettingMode.SUBPIXEL
+                text(targetHeight.toString())
+            }
+            logEnabled = false
+        }
+    }
+}
+
+data class Triangle3D(
+    val x1: Vector3,
+    val x2: Vector3,
+    val x3: Vector3,
+) {
+    fun toList(): List<Vector3> = listOf(x1, x2, x3)
+}
+
+fun connectAllSegments(segments: List<Segment3D>): List<List<Segment3D>> {
+    val remaining = segments.toMutableList()
+    val allPaths = mutableListOf<List<Segment3D>>()
+
+    while (remaining.isNotEmpty()) {
+        val path = mutableListOf<Segment3D>()
+
+        // 开始新路径
+        path.add(remaining.removeAt(0))
+
+        var changed: Boolean
+        do {
+            changed = false
+
+            // 向前扩展
+            val lastEnd = path.last().end
+            val forwardSegment = remaining.find { it.start == lastEnd || it.end == lastEnd }
+            if (forwardSegment != null) {
+                val connectedSegment = if (forwardSegment.start == lastEnd) {
+                    forwardSegment // 正向
+                } else {
+                    Segment3D(forwardSegment.end, forwardSegment.start) // 反向
+                }
+                path.add(connectedSegment)
+                remaining.remove(forwardSegment)
+                changed = true
+            }
+
+            // 向后扩展
+            val firstStart = path.first().start
+            val backwardSegment = remaining.find { it.end == firstStart || it.start == firstStart }
+            if (backwardSegment != null) {
+                val connectedSegment = if (backwardSegment.end == firstStart) {
+                    backwardSegment // 正向
+                } else {
+                    Segment3D(backwardSegment.end, backwardSegment.start) // 反向
+                }
+                path.add(0, connectedSegment)
+                remaining.remove(backwardSegment)
+                changed = true
+            }
+
+        } while (changed && remaining.isNotEmpty())
+
+        allPaths.add(path)
+    }
+
+    return allPaths
+}

desktop/src/jvmDemo/kotlin/FindContours.kt

@@ -0,0 +1,94 @@
+import org.openrndr.application
+import org.openrndr.color.ColorRGBa
+import org.openrndr.extra.camera.Camera2D
+import org.openrndr.extra.marchingsquares.findContours
+import org.openrndr.math.Vector2
+import org.openrndr.math.Vector3
+
+/**
+ * A simple demonstration of using the `findContours` method provided by `orx-marching-squares`.
+ *
+ * `findContours` lets one generate contours by providing a mathematical function to be
+ * sampled within the provided area and with the given cell size. Contours are generated
+ * between the areas in which the function returns positive and negative values.
+ *
+ * In this example, the `f` function returns the distance of a point to the center of the window minus 200.0.
+ * Therefore, sampled locations which are less than 200 pixels away from the center return
+ * negative values and all others return positive values, effectively generating a circle of radius 200.0.
+ *
+ * Try increasing the cell size to see how the precision of the circle reduces.
+ *
+ * The circular contour created in this program has over 90 segments. The number of segments depends on the cell
+ * size, and the resulting radius.
+ */
+fun main() = application {
+    configure {
+        width = 720
+        height = 720
+    }
+    program {
+        extend(Camera2D())
+        var showLog = true
+        val target = Vector2(0.0, 0.0)
+        val points3D = (0..10).map { x ->
+            (0..10).map { y ->
+                Vector3(x.toDouble(), y.toDouble(), x * y * 1.0)
+            }
+        }
+
+        extend {
+            drawer.clear(ColorRGBa.BLACK)
+            drawer.stroke = ColorRGBa.PINK
+            fun f3(v: Vector2): Double {
+                val distance = drawer.bounds.center.distanceTo(v)
+                return when (distance) {
+                    in 0.0..<100.0 -> -3.0
+                    in 100.0..<200.0 -> 1.0
+                    in 200.0..300.0 -> -1.0
+                    else -> distance
+                }
+            }
+
+            fun f(v: Vector2): Double {
+                val distanceTo = v.distanceTo(target)
+                return (distanceTo - 100.0).also {
+                    if (showLog) println(
+                        buildString {
+                            appendLine("${v} distanceTo ${target} = ${distanceTo}")
+                            appendLine("distanceTo - 100.0 = ${distanceTo - 100.0}")
+                        }
+                    )
+                }
+            }
+
+            val points = mutableListOf<Vector2>()
+
+            fun f1(v: Vector2): Double {
+                val result = if (v.x == v.y * 2 || v.x * 2 == v.y) {
+                    points.add(v)
+                    -1.0
+                } else 0.0
+                return result.also {
+                    if (showLog) {
+                        println("$v -> $result")
+                    }
+                }
+            }
+
+            val contours = findContours(::f3, drawer.bounds, 4.0)
+            drawer.fill = null
+            drawer.contours(contours)
+
+            if (showLog) {
+                println(
+                    buildString {
+                        for ((index, contour) in contours.withIndex()) {
+                            appendLine("index = ${index}, $contour")
+                        }
+                    }
+                )
+            }
+            showLog = false
+        }
+    }
+}

desktop/src/jvmDemo/kotlin/HeightmapVolcanoGenerator.kt

@@ -0,0 +1,373 @@
+import com.icegps.math.geometry.Vector3D
+import kotlin.math.PI
+import kotlin.math.abs
+import kotlin.math.cos
+import kotlin.math.exp
+import kotlin.math.floor
+import kotlin.math.max
+import kotlin.math.sin
+import kotlin.math.sqrt
+
+/**
+ * @author tabidachinokaze
+ * @date 2025/11/22
+ */
+object HeightmapVolcanoGenerator {
+
+    // 基础火山高度图
+    fun generateVolcanoHeightmap(
+        width: Int = 100,
+        height: Int = 100,
+        centerX: Double = 50.0,
+        centerY: Double = 50.0,
+        maxHeight: Double = 60.0,
+        craterRadius: Double = 8.0,
+        volcanoRadius: Double = 30.0
+    ): List<Vector3D> {
+        val points = mutableListOf<Vector3D>()
+
+        for (x in 0 until width) {
+            for (y in 0 until height) {
+                // 计算到火山中心的距离
+                val dx = x - centerX
+                val dy = y - centerY
+                val distance = sqrt(dx * dx + dy * dy)
+
+                // 计算基础火山高度
+                var z = calculateVolcanoHeight(distance, craterRadius, volcanoRadius, maxHeight)
+
+                // 添加噪声细节
+                val noise = perlinNoise(x * 0.1, y * 0.1, 0.1) * 3.0
+                z = max(0.0, z + noise)
+
+                points.add(Vector3D(x.toDouble(), y.toDouble(), z))
+            }
+        }
+
+        return points
+    }
+
+    // 复合火山群高度图
+    fun generateVolcanoClusterHeightmap(
+        width: Int = 150,
+        height: Int = 150,
+        volcanoCount: Int = 3
+    ): List<Vector3D> {
+        val points = mutableListOf<Vector3D>()
+        val volcanoes = generateRandomVolcanoPositions(volcanoCount, width, height)
+
+        for (x in (0 until width).step(25)) {
+            for (y in (0 until height).step(25)) {
+                var totalZ = 0.0
+
+                // 叠加所有火山的影响
+                for (volcano in volcanoes) {
+                    val dx = x - volcano.x
+                    val dy = y - volcano.y
+                    val distance = sqrt(dx * dx + dy * dy)
+
+                    if (distance <= volcano.radius) {
+                        val volcanoHeight = calculateVolcanoHeight(
+                            distance,
+                            volcano.craterRadius,
+                            volcano.radius,
+                            volcano.maxHeight
+                        )
+                        totalZ += volcanoHeight
+                    }
+                }
+
+                // 基础地形
+                val baseNoise = perlinNoise(x * 0.02, y * 0.02, 0.05) * 5.0
+                val detailNoise = perlinNoise(x * 0.1, y * 0.1, 0.2) * 2.0
+
+                points.add(Vector3D(x.toDouble(), y.toDouble(), totalZ + baseNoise + detailNoise))
+            }
+        }
+
+        return points
+    }
+
+    // 带熔岩流的火山高度图
+    fun generateVolcanoWithLavaHeightmap(
+        width: Int = 100,
+        height: Int = 100
+    ): List<Vector3D> {
+        val points = mutableListOf<Vector3D>()
+        val centerX = width / 2.0
+        val centerY = height / 2.0
+
+        // 生成熔岩流路径
+        val lavaFlows = generateLavaFlowPaths(centerX, centerY, 3)
+
+        for (x in 0 until width) {
+            for (y in 0 until height) {
+                val dx = x - centerX
+                val dy = y - centerY
+                val distance = sqrt(dx * dx + dy * dy)
+
+                // 基础火山高度
+                var z = calculateVolcanoHeight(distance, 10.0, 35.0, 70.0)
+
+                // 添加熔岩流
+                z += calculateLavaFlowEffect(x.toDouble(), y.toDouble(), lavaFlows)
+
+                // 侵蚀效果
+                z += calculateErosionEffect(x.toDouble(), y.toDouble(), distance, z)
+
+                points.add(Vector3D(x.toDouble(), y.toDouble(), max(0.0, z)))
+            }
+        }
+
+        return points
+    }
+
+    // 破火山口高度图
+    fun generateCalderaHeightmap(
+        width: Int = 100,
+        height: Int = 100
+    ): List<Vector3D> {
+        val points = mutableListOf<Vector3D>()
+        val centerX = width / 2.0
+        val centerY = height / 2.0
+
+        for (x in 0 until width) {
+            for (y in 0 until height) {
+                val dx = x - centerX
+                val dy = y - centerY
+                val distance = sqrt(dx * dx + dy * dy)
+
+                var z = calculateCalderaHeight(distance, 15.0, 45.0, 50.0)
+
+                // 内部平坦区域细节
+                if (distance < 20) {
+                    z += perlinNoise(x * 0.2, y * 0.2, 0.3) * 1.5
+                }
+
+                points.add(Vector3D(x.toDouble(), y.toDouble(), max(0.0, z)))
+            }
+        }
+
+        return points
+    }
+
+    // 线性火山链高度图
+    fun generateVolcanoChainHeightmap(
+        width: Int = 200,
+        height: Int = 100
+    ): List<Vector3D> {
+        val points = mutableListOf<Vector3D>()
+
+        // 在一条线上生成多个火山
+        val chainCenters = listOf(
+            Vector3D(30.0, 50.0, 0.0),
+            Vector3D(70.0, 50.0, 0.0),
+            Vector3D(110.0, 50.0, 0.0),
+            Vector3D(150.0, 50.0, 0.0),
+            Vector3D(170.0, 50.0, 0.0)
+        )
+
+        for (x in 0 until width) {
+            for (y in 0 until height) {
+                var totalZ = 0.0
+
+                for (center in chainCenters) {
+                    val dx = x - center.x
+                    val dy = y - center.y
+                    val distance = sqrt(dx * dx + dy * dy)
+
+                    if (distance <= 25.0) {
+                        val volcanoZ = calculateVolcanoHeight(distance, 6.0, 25.0, 40.0)
+                        totalZ += volcanoZ
+                    }
+                }
+
+                // 添加基底地形，模拟山脉链
+                val baseRidge = calculateMountainRidge(x.toDouble(), y.toDouble(), width, height)
+                totalZ += baseRidge
+
+                points.add(Vector3D(x.toDouble(), y.toDouble(), totalZ))
+            }
+        }
+
+        return points
+    }
+
+    // 辅助函数
+    private data class VolcanoInfo(
+        val x: Double,
+        val y: Double,
+        val radius: Double,
+        val craterRadius: Double,
+        val maxHeight: Double
+    )
+
+    private data class LavaFlowInfo(
+        val startX: Double,
+        val startY: Double,
+        val angle: Double, // 弧度
+        val length: Double,
+        val width: Double,
+        val intensity: Double
+    )
+
+    private fun calculateVolcanoHeight(
+        distance: Double,
+        craterRadius: Double,
+        volcanoRadius: Double,
+        maxHeight: Double
+    ): Double {
+        return when {
+            distance <= craterRadius -> {
+                // 火山口 - 中心凹陷
+                val craterDepth = maxHeight * 0.4
+                craterDepth * (1.0 - distance / craterRadius)
+            }
+
+            distance <= volcanoRadius -> {
+                // 火山锥
+                val slopeDistance = distance - craterRadius
+                val maxSlopeDistance = volcanoRadius - craterRadius
+                val normalized = slopeDistance / maxSlopeDistance
+                maxHeight * (1.0 - normalized * normalized)
+            }
+
+            else -> 0.0
+        }
+    }
+
+    private fun calculateCalderaHeight(
+        distance: Double,
+        innerRadius: Double,
+        outerRadius: Double,
+        rimHeight: Double
+    ): Double {
+        return when {
+            distance <= innerRadius -> {
+                // 平坦的破火山口底部
+                rimHeight * 0.2
+            }
+
+            distance <= outerRadius -> {
+                // 陡峭的边缘
+                val rimDistance = distance - innerRadius
+                val rimWidth = outerRadius - innerRadius
+                val normalized = rimDistance / rimWidth
+                rimHeight * (1.0 - (1.0 - normalized) * (1.0 - normalized))
+            }
+
+            else -> {
+                // 外部平缓斜坡
+                val externalDistance = distance - outerRadius
+                rimHeight * exp(-externalDistance * 0.08)
+            }
+        }
+    }
+
+    private fun calculateLavaFlowEffect(x: Double, y: Double, lavaFlows: List<LavaFlowInfo>): Double {
+        var effect = 0.0
+
+        for (flow in lavaFlows) {
+            val dx = x - flow.startX
+            val dy = y - flow.startY
+
+            // 计算到熔岩流中心线的距离
+            val flowDirX = cos(flow.angle)
+            val flowDirY = sin(flow.angle)
+
+            val projection = dx * flowDirX + dy * flowDirY
+
+            if (projection in 0.0..flow.length) {
+                val perpendicularX = dx - projection * flowDirX
+                val perpendicularY = dy - projection * flowDirY
+                val perpendicularDist = sqrt(perpendicularX * perpendicularX + perpendicularY * perpendicularY)
+
+                if (perpendicularDist <= flow.width) {
+                    val widthFactor = 1.0 - (perpendicularDist / flow.width)
+                    val lengthFactor = 1.0 - (projection / flow.length)
+                    effect += flow.intensity * widthFactor * lengthFactor
+                }
+            }
+        }
+
+        return effect
+    }
+
+    private fun calculateErosionEffect(x: Double, y: Double, distance: Double, height: Double): Double {
+        // 基于坡度的侵蚀
+        val slopeNoise = perlinNoise(x * 0.15, y * 0.15, 0.1) * 2.0
+        // 基于距离的侵蚀
+        val distanceErosion = if (distance > 25) perlinNoise(x * 0.08, y * 0.08, 0.05) * 1.5 else 0.0
+        return slopeNoise + distanceErosion
+    }
+
+    private fun calculateMountainRidge(x: Double, y: Double, width: Int, height: Int): Double {
+        // 创建山脉基底
+        val ridgeCenter = height / 2.0
+        val distanceToRidge = abs(y - ridgeCenter)
+        val ridgeWidth = height * 0.3
+
+        if (distanceToRidge <= ridgeWidth) {
+            val ridgeFactor = 1.0 - (distanceToRidge / ridgeWidth)
+            return ridgeFactor * 15.0 * perlinNoise(x * 0.01, y * 0.01, 0.02)
+        }
+        return 0.0
+    }
+
+    private fun generateRandomVolcanoPositions(count: Int, width: Int, height: Int): List<VolcanoInfo> {
+        return List(count) {
+            VolcanoInfo(
+                x = (width * 0.2 + random() * width * 0.6),
+                y = (height * 0.2 + random() * height * 0.6),
+                radius = 20.0 + random() * 20.0,
+                craterRadius = 5.0 + random() * 7.0,
+                maxHeight = 25.0 + random() * 35.0
+            )
+        }
+    }
+
+    private fun generateLavaFlowPaths(centerX: Double, centerY: Double, count: Int): List<LavaFlowInfo> {
+        return List(count) {
+            LavaFlowInfo(
+                startX = centerX,
+                startY = centerY,
+                angle = random() * 2 * PI,
+                length = 20.0 + random() * 15.0,
+                width = 2.0 + random() * 3.0,
+                intensity = 5.0 + random() * 8.0
+            )
+        }
+    }
+
+    private fun perlinNoise(x: Double, y: Double, frequency: Double): Double {
+        // 简化的柏林噪声实现
+        val x0 = floor(x * frequency)
+        val y0 = floor(y * frequency)
+        val x1 = x0 + 1
+        val y1 = y0 + 1
+
+        fun grad(ix: Int, iy: Int): Double {
+            val random = sin(ix * 12.9898 + iy * 78.233) * 43758.5453
+            return (random % 1.0) * 2 - 1
+        }
+
+        fun interpolate(a: Double, b: Double, w: Double): Double {
+            return a + (b - a) * (w * w * (3 - 2 * w))
+        }
+
+        val g00 = grad(x0.toInt(), y0.toInt())
+        val g10 = grad(x1.toInt(), y0.toInt())
+        val g01 = grad(x0.toInt(), y1.toInt())
+        val g11 = grad(x1.toInt(), y1.toInt())
+
+        val tx = x * frequency - x0
+        val ty = y * frequency - y0
+
+        val n0 = interpolate(g00, g10, tx)
+        val n1 = interpolate(g01, g11, tx)
+
+        return interpolate(n0, n1, ty)
+    }
+
+    private fun random(): Double = Math.random()
+}

gradle.properties

@@ -21,4 +21,13 @@ kotlin.mpp.import.legacyTestSourceSetDetection=true
 
 # Enable Dokka 2.0.0
 org.jetbrains.dokka.experimental.gradle.pluginMode=V2Enabled
-org.jetbrains.dokka.experimental.gradle.pluginMode.noWarn=true
+org.jetbrains.dokka.experimental.gradle.pluginMode.noWarn=true
+
+# AndroidX package structure to make it clearer which packages are bundled with the
+# Android operating system, and which are packaged with your app's APK
+# https://developer.android.com/topic/libraries/support-library/androidx-rn
+android.useAndroidX=true
+# Enables namespacing of each library's R class so that its R class includes only the
+# resources declared in the library itself and none from the library's dependencies,
+# thereby reducing the size of the R class for that library
+android.nonTransitiveRClass=true

gradle/libs.versions.toml

@@ -16,15 +16,27 @@ antlrKotlin = "1.0.8"
 minim = "2.2.2"
 netty = "4.2.7.Final"
 rabbitcontrol = "0.3.39"
-zxing = "3.5.3"
+zxing = "3.5.4"
-ktor = "3.3.1"
+ktor = "3.3.2"
 jgit = "7.3.0.202506031305-r"
 javaosc = "0.9"
 jsoup = "1.21.2"
 mockk = "1.14.2"
 processing = "4.4.10"
-nmcp = "1.1.0"
+nmcp = "1.2.0"
 okhttp = "5.2.1"
+agp = "8.13.1"
+junit = "4.13.2"
+coreKtx = "1.17.0"
+junitVersion = "1.3.0"
+espressoCoreVersion = "3.7.0"
+appcompat = "1.7.1"
+material = "1.13.0"
+activity = "1.12.0"
+constraintlayout = "2.2.1"
+lifecycleRuntimeKtx = "2.10.0"
+kotlinx-serialization = "1.9.0"
+mapbox = "11.16.6"
 
 [libraries]
 okhttp = { group = "com.squareup.okhttp3", name = "okhttp", version.ref = "okhttp" }
@@ -57,10 +69,32 @@ antlr-core = { group = "org.antlr", name = "antlr4", version.ref = "antlr" }
 antlr-runtime = { group = "org.antlr", name = "antlr4-runtime", version.ref = "antlr" }
 antlr-kotlin-runtime = { group = "com.strumenta", name = "antlr-kotlin-runtime", version.ref = "antlrKotlin" }
 jsoup = { group = "org.jsoup", name = "jsoup", version.ref = "jsoup" }
+junit = { group = "junit", name = "junit", version.ref = "junit" }
+core-ktx = { group = "androidx.core", name = "core-ktx", version.ref = "coreKtx" }
+ext-junit = { group = "androidx.test.ext", name = "junit", version.ref = "junitVersion" }
+androidx-espresso-core = { group = "androidx.test.espresso", name = "espresso-core", version.ref = "espressoCoreVersion" }
+androidx-appcompat = { group = "androidx.appcompat", name = "appcompat", version.ref = "appcompat" }
+material = { group = "com.google.android.material", name = "material", version.ref = "material" }
+androidx-activity = { group = "androidx.activity", name = "activity", version.ref = "activity" }
+androidx-constraintlayout = { group = "androidx.constraintlayout", name = "constraintlayout", version.ref = "constraintlayout" }
+androidx-lifecycle-runtime-ktx = { group = "androidx.lifecycle", name = "lifecycle-runtime-ktx", version.ref = "lifecycleRuntimeKtx" }
+#ktor
+ktor-client-core = { group = "io.ktor", name = "ktor-client-core", version.ref = "ktor" }
+ktor-client-cio = { group = "io.ktor", name = "ktor-client-cio", version.ref = "ktor" }
+ktor-client-logging = { group = "io.ktor", name = "ktor-client-logging", version.ref = "ktor" }
+ktor-client-content-negotiation = { group = "io.ktor", name = "ktor-client-content-negotiation", version.ref = "ktor" }
+ktor-serialization-kotlinx-json = { group = "io.ktor", name = "ktor-serialization-kotlinx-json", version.ref = "ktor" }
+
+kotlinx-serialization-json = { group = "org.jetbrains.kotlinx", name = "kotlinx-serialization-json", version.ref = "kotlinx-serialization" }
+mapbox-maps = { module = "com.mapbox.maps:android-ndk27", version.ref = "mapbox" }
 
 [plugins]
 kotlin-serialization = { id = "org.jetbrains.kotlin.plugin.serialization", version.ref = "kotlin" }
 nebula-release = { id = "nebula.release", version.ref = "nebulaRelease" }
 kotest-multiplatform = { id = "io.kotest.multiplatform", version.ref = "kotest" }
 antlr-kotlin = { id = "com.strumenta.antlr-kotlin", version.ref = "antlrKotlin" }
-nmcp = { id = "com.gradleup.nmcp.aggregation", version.ref = "nmcp" }
+nmcp = { id = "com.gradleup.nmcp.aggregation", version.ref = "nmcp" }
+android-application = { id = "com.android.application", version.ref = "agp" }
+kotlin-android = { id = "org.jetbrains.kotlin.android", version.ref = "kotlin" }
+kotlin-jvm = { id = "org.jetbrains.kotlin.jvm", version.ref = "kotlin" }
+android-library = { id = "com.android.library", version.ref = "agp" }

gradle/wrapper/gradle-wrapper.properties

@@ -1,6 +1,6 @@
 distributionBase=GRADLE_USER_HOME
 distributionPath=wrapper/dists
-distributionUrl=https\://services.gradle.org/distributions/gradle-8.14.3-bin.zip
+distributionUrl=https\://services.gradle.org/distributions/gradle-9.2.0-bin.zip
 networkTimeout=10000
 validateDistributionUrl=true
 zipStoreBase=GRADLE_USER_HOME

gradlew

@@ -1,7 +1,7 @@
 #!/bin/sh
 
 #
-# Copyright © 2015 the original authors.
+# Copyright © 2015-2021 the original authors.
 #
 # Licensed under the Apache License, Version 2.0 (the "License");
 # you may not use this file except in compliance with the License.

icegps-common/.gitignore

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

icegps-common/build.gradle.kts

@@ -0,0 +1,44 @@
+import org.jetbrains.kotlin.gradle.dsl.JvmTarget
+
+plugins {
+    alias(libs.plugins.android.library)
+    alias(libs.plugins.kotlin.android)
+}
+
+android {
+    namespace = "com.icegps.common"
+    compileSdk {
+        version = release(36)
+    }
+
+    defaultConfig {
+        minSdk = 28
+
+        testInstrumentationRunner = "androidx.test.runner.AndroidJUnitRunner"
+        consumerProguardFiles("consumer-rules.pro")
+    }
+
+    buildTypes {
+        release {
+            isMinifyEnabled = false
+            proguardFiles(
+                getDefaultProguardFile("proguard-android-optimize.txt"),
+                "proguard-rules.pro"
+            )
+        }
+    }
+    compileOptions {
+        sourceCompatibility = JavaVersion.VERSION_17
+        targetCompatibility = JavaVersion.VERSION_17
+    }
+}
+
+kotlin {
+    compilerOptions.jvmTarget = JvmTarget.JVM_17
+}
+
+dependencies {
+    testImplementation(libs.junit)
+    androidTestImplementation(libs.ext.junit)
+    androidTestImplementation(libs.androidx.espresso.core)
+}

icegps-common/proguard-rules.pro

@@ -0,0 +1,21 @@
+# Add project specific ProGuard rules here.
+# You can control the set of applied configuration files using the
+# proguardFiles setting in build.gradle.
+#
+# For more details, see
+#   http://developer.android.com/guide/developing/tools/proguard.html
+
+# If your project uses WebView with JS, uncomment the following
+# and specify the fully qualified class name to the JavaScript interface
+# class:
+#-keepclassmembers class fqcn.of.javascript.interface.for.webview {
+#   public *;
+#}
+
+# Uncomment this to preserve the line number information for
+# debugging stack traces.
+#-keepattributes SourceFile,LineNumberTable
+
+# If you keep the line number information, uncomment this to
+# hide the original source file name.
+#-renamesourcefileattribute SourceFile

icegps-common/src/androidTest/java/com/icegps/common/ExampleInstrumentedTest.kt

@@ -0,0 +1,24 @@
+package com.icegps.common
+
+import androidx.test.platform.app.InstrumentationRegistry
+import androidx.test.ext.junit.runners.AndroidJUnit4
+
+import org.junit.Test
+import org.junit.runner.RunWith
+
+import org.junit.Assert.*
+
+/**
+ * Instrumented test, which will execute on an Android device.
+ *
+ * See [testing documentation](http://d.android.com/tools/testing).
+ */
+@RunWith(AndroidJUnit4::class)
+class ExampleInstrumentedTest {
+    @Test
+    fun useAppContext() {
+        // Context of the app under test.
+        val appContext = InstrumentationRegistry.getInstrumentation().targetContext
+        assertEquals("com.icegps.common.test", appContext.packageName)
+    }
+}

icegps-common/src/main/AndroidManifest.xml

@@ -0,0 +1,4 @@
+<?xml version="1.0" encoding="utf-8"?>
+<manifest xmlns:android="http://schemas.android.com/apk/res/android">
+
+</manifest>

icegps-common/src/main/java/com/icegps/common/helper/BlhToEnu.kt

@@ -0,0 +1,274 @@
+package com.icegps.common.helper
+
+import kotlin.math.atan2
+import kotlin.math.cos
+import kotlin.math.pow
+import kotlin.math.sin
+import kotlin.math.sqrt
+
+/**
+ * BLH -> ENU
+ *
+ * @author lm
+ * @date 2025/3/12
+ * @link https://gist.github.com/komasaru/6ce0634475923ddac597f868288c54e9
+ */
+class BlhToEnu {
+
+    companion object {
+        private const val PI_180 = Math.PI / 180.0
+
+        // WGS84 坐标参数
+        private const val A = 6378137.0 // a（地球椭球体长半径（赤道面平均半径））
+        private const val ONE_F = 298.257223563 // 1 / f（地球椭球体扁平率 = (a - b) / a）
+        private val B = A * (1.0 - 1.0 / ONE_F) // b（地球椭球体短半径）
+        private val E2 = (1.0 / ONE_F) * (2 - (1.0 / ONE_F))
+
+        // e^2 = 2 * f - f * f = (a^2 - b^2) / a^2
+        private val ED2 = E2 * A * A / (B * B) // e'^2 = (a^2 - b^2) / b^2
+    }
+
+    private var originLat: Double = 0.0
+    private var originLon: Double = 0.0
+    private var originHeight: Double = 0.0
+    private var isOriginSet: Boolean = false
+
+    fun getOriginLat(): Double = originLat
+
+    fun getOriginLon(): Double = originLon
+
+    fun getOriginHeight(): Double = originHeight
+
+    fun resetEnuBenchmarkPoint() {
+        isOriginSet = false
+    }
+
+    fun wgs84ToEnu(lat: Double, lon: Double, height: Double = 0.0): DoubleArray {
+        if (!isOriginSet) {
+            originLat = lat
+            originLon = lon
+            originHeight = height
+            isOriginSet = true
+            return doubleArrayOf(0.0, 0.0, 0.0)
+        }
+
+        val enu = blh2enu(originLat, originLon, originHeight, lat, lon, height)
+//        var az = atan2(enu[0], enu[1]) * 180.0 / Math.PI
+//        if (az < 0.0) {
+//            az += 360.0
+//        }
+//        val el = atan2(
+//            enu[2],
+//            sqrt(enu[0] * enu[0] + enu[1] * enu[1])
+//        ) * 180.0 / Math.PI
+//        val dst = sqrt(enu.sumOf { it * it })
+
+//        println("--->")
+//        println(
+//            """
+//            ENU: E = ${enu[0].format(3)}m
+//                 N = ${enu[1].format(3)}m
+//                 U = ${enu[2].format(3)}m
+//            方位角 = ${az.format(3)}°
+//              仰角 = ${el.format(3)}°
+//              距离 = ${dst.format(3)}m
+//        """.trimIndent()
+//        )
+
+        return enu
+    }
+
+    fun enuToWgs84(e: Double, n: Double, u: Double): DoubleArray {
+        if (!isOriginSet) {
+            return doubleArrayOf(0.0, 0.0, 0.0)
+        }
+
+        val blh = enu2blh(originLat, originLon, originHeight, e, n, u)
+
+//        println("--->")
+//        println(
+//            """
+//            BLH: Beta = ${blh[0].format(8)}°
+//                 Lambda = ${blh[1].format(8)}°
+//                 Height = ${blh[2].format(3)}m
+//        """.trimIndent()
+//        )
+
+        return blh
+    }
+
+    private fun Double.format(digits: Int) = "%.${digits}f".format(this)
+
+    /**
+     * BLH -> ENU 转换(East, North, Up)
+     *
+     * @param bO 原点 Beta(纬度)
+     * @param lO 原点 Lambda(经度)
+     * @param hO 原点 Height(高度)
+     * @param b 目标点 Beta(纬度)
+     * @param l 目标点 Lambda(经度)
+     * @param h 目标点 Height(高度)
+     * @return ENU 坐标 [e, n, u]
+     */
+    private fun blh2enu(bO: Double, lO: Double, hO: Double, b: Double, l: Double, h: Double): DoubleArray {
+        val (xO, yO, zO) = blh2ecef(bO, lO, hO)
+        val (x, y, z) = blh2ecef(b, l, h)
+        val mat0 = matZ(90.0)
+        val mat1 = matY(90.0 - bO)
+        val mat2 = matZ(lO)
+        val mat = mulMat(mulMat(mat0, mat1), mat2)
+        return rotate(mat, doubleArrayOf(x - xO, y - yO, z - zO))
+    }
+
+    /**
+     * BLH -> ECEF 转换
+     *
+     * @param lat 纬度
+     * @param lon 经度
+     * @param height 高度
+     * @return ECEF 坐标 [x, y, z]
+     */
+    private fun blh2ecef(lat: Double, lon: Double, height: Double): DoubleArray {
+        val n = { x: Double -> A / sqrt(1.0 - E2 * sin(x * PI_180).pow(2)) }
+        val x = (n(lat) + height) * cos(lat * PI_180) * cos(lon * PI_180)
+        val y = (n(lat) + height) * cos(lat * PI_180) * sin(lon * PI_180)
+        val z = (n(lat) * (1.0 - E2) + height) * sin(lat * PI_180)
+        return doubleArrayOf(x, y, z)
+    }
+
+    /**
+     * ENU -> BLH 转换
+     *
+     * @param e East 坐标
+     * @param n North 坐标
+     * @param u Up 坐标
+     * @return WGS84 坐标 [纬度, 经度, 高度]
+     */
+    private fun enu2blh(bO: Double, lO: Double, hO: Double, e: Double, n: Double, u: Double): DoubleArray {
+        val mat0 = matZ(-lO)
+        val mat1 = matY(-(90.0 - bO))
+        val mat2 = matZ(-90.0)
+        val mat = mulMat(mulMat(mat0, mat1), mat2)
+
+        val enu = doubleArrayOf(e, n, u)
+        val xyz = rotate(mat, enu)
+
+        val (xO, yO, zO) = blh2ecef(bO, lO, hO)
+        val x = xyz[0] + xO
+        val y = xyz[1] + yO
+        val z = xyz[2] + zO
+
+        return ecef2blh(x, y, z)
+    }
+
+    /**
+     * ECEF -> BLH 转换
+     *
+     * @param x ECEF X 坐标
+     * @param y ECEF Y 坐标
+     * @param z ECEF Z 坐标
+     * @return WGS84 坐标 [纬度, 经度, 高度]
+     */
+    private fun ecef2blh(x: Double, y: Double, z: Double): DoubleArray {
+        val p = sqrt(x * x + y * y)
+        val theta = atan2(z * A, p * B)
+        val sinTheta = sin(theta)
+        val cosTheta = cos(theta)
+
+        val lat = atan2(
+            z + ED2 * B * sinTheta.pow(3),
+            p - E2 * A * cosTheta.pow(3)
+        )
+        val lon = atan2(y, x)
+
+        val sinLat = sin(lat)
+        val n = A / sqrt(1.0 - E2 * sinLat * sinLat)
+        val h = p / cos(lat) - n
+
+        return doubleArrayOf(
+            lat * 180.0 / Math.PI,
+            lon * 180.0 / Math.PI,
+            h
+        )
+    }
+
+    /**
+     * 以 x 轴为轴的旋转矩阵
+     *
+     * @param ang 旋转角度（°）
+     * @return 旋转矩阵（3x3）
+     */
+    private fun matX(ang: Double): Array<DoubleArray> {
+        val a = ang * PI_180
+        val c = cos(a)
+        val s = sin(a)
+        return arrayOf(
+            doubleArrayOf(1.0, 0.0, 0.0),
+            doubleArrayOf(0.0, c, s),
+            doubleArrayOf(0.0, -s, c)
+        )
+    }
+
+    /**
+     * 以 y 轴为轴的旋转矩阵
+     *
+     * @param ang 旋转角度（°）
+     * @return 旋转矩阵（3x3）
+     */
+    private fun matY(ang: Double): Array<DoubleArray> {
+        val a = ang * PI_180
+        val c = cos(a)
+        val s = sin(a)
+        return arrayOf(
+            doubleArrayOf(c, 0.0, -s),
+            doubleArrayOf(0.0, 1.0, 0.0),
+            doubleArrayOf(s, 0.0, c)
+        )
+    }
+
+    /**
+     * 以 z 轴为轴的旋转矩阵
+     *
+     * @param ang 旋转角度（°）
+     * @return 旋转矩阵（3x3）
+     */
+    private fun matZ(ang: Double): Array<DoubleArray> {
+        val a = ang * PI_180
+        val c = cos(a)
+        val s = sin(a)
+        return arrayOf(
+            doubleArrayOf(c, s, 0.0),
+            doubleArrayOf(-s, c, 0.0),
+            doubleArrayOf(0.0, 0.0, 1.0)
+        )
+    }
+
+    /**
+     * 两个矩阵(3x3)的乘积
+     *
+     * @param matA 3x3 矩阵
+     * @param matB 3x3 矩阵
+     * @return 3x3 矩阵
+     */
+    private fun mulMat(matA: Array<DoubleArray>, matB: Array<DoubleArray>): Array<DoubleArray> {
+        return Array(3) { k ->
+            DoubleArray(3) { j ->
+                (0..2).sumOf { i -> matA[k][i] * matB[i][j] }
+            }
+        }
+    }
+
+    /**
+     * 点的旋转
+     *
+     * @param mat 3x3 旋转矩阵
+     * @param pt 旋转前坐标 [x, y, z]
+     * @return 旋转后坐标 [x, y, z]
+     */
+    private fun rotate(mat: Array<DoubleArray>, pt: DoubleArray): DoubleArray {
+        return DoubleArray(3) { j ->
+            (0..2).sumOf { i -> mat[j][i] * pt[i] }
+        }
+    }
+
+}

icegps-common/src/main/java/com/icegps/common/helper/GeoHelper.kt

@@ -0,0 +1,419 @@
+package com.icegps.common.helper
+
+import android.os.Parcel
+import android.os.Parcelable
+import kotlin.math.atan
+import kotlin.math.atan2
+import kotlin.math.cbrt
+import kotlin.math.cos
+import kotlin.math.exp
+import kotlin.math.ln
+import kotlin.math.sin
+import kotlin.math.sqrt
+import kotlin.math.tan
+
+/**
+ * WGS84、EPSG3857、ENU 的坐标转换工具类
+ *
+ * @author lm
+ * @date 2024/8/2
+ */
+class GeoHelper private constructor() {
+
+    companion object {
+        private var sharedInstance: GeoHelper? = null
+        fun getSharedInstance(): GeoHelper = sharedInstance ?: GeoHelper().also { sharedInstance = it }
+        fun createInstance(): GeoHelper = GeoHelper()
+    }
+
+    // WGS-84 ellipsoid parameters
+    private val RADIUS = 6378137.0 // Major radius
+    private val RADIUS_B = 6356752.314245 // Minor radius
+    private val E = (RADIUS * RADIUS - RADIUS_B * RADIUS_B) / (RADIUS * RADIUS) // Eccentricity
+    private val HALF_SIZE = Math.PI * RADIUS // Half circumference of Earth
+    private val DEG2RAD = Math.PI / 180 // Degrees to radians conversion factor
+    private val RAD2DEG = 180 / Math.PI // Radians to degrees conversion factor
+    private val RE_WGS84 = 6378137.0 // Earth's equatorial radius in WGS84
+    private val FE_WGS84 = 1.0 / 298.257223563 // Flattening of the WGS84 ellipsoid
+
+    private var isFirstPoint = true
+    private var firstPoint = DoubleArray(3)
+    private val bPos = DoubleArray(3)
+    private var bECEF = DoubleArray(3)
+    private val rPos = DoubleArray(3)
+    private var rECEF = DoubleArray(3)
+    private val vECEF = DoubleArray(3)
+
+    private var useBlhToEnu = true
+    private var blhToEnu = BlhToEnu()
+
+    /**
+     * 将 WGS84 坐标转换为 ENU (East-North-Up) 坐标
+     * 如果是第一个点，它将被设置为 ENU 坐标系的基准点
+     *
+     * @param lon 经度（度）
+     * @param lat 纬度（度）
+     * @param hgt 高度（米）
+     * @return 包含 ENU 坐标的 Enu 对象
+     */
+    fun wgs84ToENU(lon: Double, lat: Double, hgt: Double): ENU {
+        if (useBlhToEnu) {
+            val enu = blhToEnu.wgs84ToEnu(lon = lon, lat = lat, height = hgt)
+            return ENU(enu[0], enu[1], enu[2])
+        }
+
+        if (isFirstPoint) setEnuBenchmark(lon, lat, hgt)
+        rPos[0] = lat * DEG2RAD
+        rPos[1] = lon * DEG2RAD
+        rPos[2] = hgt
+
+        rECEF = pos2ecef(rPos)
+        vECEF[0] = rECEF[0] - bECEF[0]
+        vECEF[1] = rECEF[1] - bECEF[1]
+        vECEF[2] = rECEF[2] - bECEF[2]
+        val enuDoubleArray = ecef2enu(bPos, vECEF)
+        return ENU(enuDoubleArray[0], enuDoubleArray[1], enuDoubleArray[2])
+    }
+
+    /**
+     * 将 WGS84 坐标转换为 ENU (East-North-Up) 坐标
+     * 如果是第一个点，它将被设置为 ENU 坐标系的基准点
+     *
+     * @param wgs84 WGS84 坐标对象
+     * @return 包含 ENU 坐标的 Enu 对象
+     */
+    fun wgs84ObjectToENU(wgs84: WGS84): ENU = wgs84ToENU(wgs84.lon, wgs84.lat, wgs84.hgt)
+
+    /**
+     * 是否已设置 ENU 坐标系的基准点
+     */
+    fun isEnuBenchmarkSet(): Boolean = !isFirstPoint
+
+    /**
+     * 设置 ENU 坐标系的基准点
+     *
+     * @param lon 基准点经度（度）
+     * @param lat 基准点纬度（度）
+     * @param hgt 基准点高度（米）
+     */
+    private fun setEnuBenchmark(lon: Double, lat: Double, hgt: Double) {
+        firstPoint = doubleArrayOf(lon, lat, hgt)
+        bPos[0] = lat * DEG2RAD
+        bPos[1] = lon * DEG2RAD
+        bPos[2] = hgt
+        bECEF = pos2ecef(bPos)
+        isFirstPoint = false
+    }
+
+    /**
+     * 获取 ENU 坐标系的基准点
+     *
+     * @return 包含 WGS84 坐标 {经度, 纬度, 高度} 的 DoubleArray
+     */
+    fun getEnuBenchmarkPoint(): DoubleArray {
+        if (useBlhToEnu) {
+            return doubleArrayOf(blhToEnu.getOriginLon(), blhToEnu.getOriginLat(), blhToEnu.getOriginHeight())
+        }
+        return firstPoint
+    }
+
+    /**
+     * 获取 ENU 坐标系的基准点
+     *
+     * @return 包含 WGS84 坐标的 WGS84 对象
+     */
+    fun getEnuBenchmarkPointAsWGS84(): WGS84 {
+        if (useBlhToEnu) {
+            return WGS84(blhToEnu.getOriginLon(), blhToEnu.getOriginLat(), blhToEnu.getOriginHeight())
+        }
+        return WGS84(firstPoint[0], firstPoint[1], firstPoint[2])
+    }
+
+    /**
+     * 重置 ENU 基准点
+     * 调用此方法后，下一次 wgs84ToENU 调用将设置新的基准点
+     */
+    fun resetEnuBenchmarkPoint() {
+        if (useBlhToEnu) {
+            blhToEnu.resetEnuBenchmarkPoint()
+            return
+        }
+        isFirstPoint = true
+    }
+
+    /**
+     * 将 ENU (East-North-Up) 坐标转换为 WGS84 坐标
+     *
+     * @param enu 包含 ENU 坐标的 Enu 对象
+     * @return 包含 WGS84 坐标 {经度, 纬度, 高度} 的 DoubleArray
+     */
+    fun enuToWGS84(enu: ENU): DoubleArray {
+        if (useBlhToEnu) {
+            val wgs84 = blhToEnu.enuToWgs84(e = enu.x, n = enu.y, u = enu.z)
+            return doubleArrayOf(wgs84[1], wgs84[0], wgs84[2])
+        }
+
+        val enuArray = doubleArrayOf(enu.x, enu.y, enu.z)
+        val enuToEcefMatrix = xyz2enu(bPos)
+        val ecefArray = matmul(charArrayOf('T', 'N'), 3, 1, 3, 1.0, enuToEcefMatrix, enuArray, 0.0)
+        vECEF[0] = bECEF[0] + ecefArray[0]
+        vECEF[1] = bECEF[1] + ecefArray[1]
+        vECEF[2] = bECEF[2] + ecefArray[2]
+        return ecef2pos(vECEF)
+    }
+
+    /**
+     * 将 ENU (East-North-Up) 坐标转换为 WGS84 坐标
+     *
+     * @param enu 包含 ENU 坐标的 Enu 对象
+     * @return 包含 WGS84 坐标的 WGS84 对象
+     */
+    fun enuToWGS84Object(enu: ENU): WGS84 {
+        val wgs84Array = enuToWGS84(enu)
+        return WGS84(wgs84Array[0], wgs84Array[1], wgs84Array[2])
+    }
+
+    /**
+     * 将 WGS84 坐标转换为 EPSG3857 坐标
+     *
+     * @param lon 经度（度）
+     * @param lat 纬度（度）
+     * @return 包含 EPSG3857 坐标的 EPSG3857 对象
+     */
+    fun wgs84ToEPSG3857(lon: Double, lat: Double): EPSG3857 {
+        val x = lon * HALF_SIZE / 180
+        var y = RADIUS * ln(tan(Math.PI * (lat + 90) / 360))
+        y = y.coerceIn(-HALF_SIZE, HALF_SIZE)
+        return EPSG3857(x, y)
+    }
+
+    /**
+     * 将 WGS84 坐标转换为 EPSG3857 坐标
+     *
+     * @param wgs84 WGS84 坐标对象
+     * @return 包含 EPSG3857 坐标的 EPSG3857 对象
+     */
+    fun wgs84ObjectToEPSG3857(wgs84: WGS84): EPSG3857 = wgs84ToEPSG3857(wgs84.lon, wgs84.lat)
+
+    /**
+     * 将 EPSG3857 坐标转换为 WGS84 坐标
+     *
+     * @param epsg3857 包含 EPSG3857 坐标的 EPSG3857 对象
+     * @return 包含 WGS84 坐标 {经度, 纬度} 的 DoubleArray
+     */
+    fun epsg3857ToWGS84(epsg3857: EPSG3857): DoubleArray {
+        val lon = (epsg3857.x / HALF_SIZE) * 180.0
+        val lat = (2 * atan(exp(epsg3857.y / RADIUS)) - Math.PI / 2) * RAD2DEG
+        return doubleArrayOf(lon, lat)
+    }
+
+    /**
+     * 将 EPSG3857 坐标转换为 WGS84 坐标
+     *
+     * @param epsg3857 包含 EPSG3857 坐标的 EPSG3857 对象
+     * @return 包含 WGS84 坐标的 WGS84 对象
+     */
+    fun epsg3857ToWGS84Object(epsg3857: EPSG3857): WGS84 {
+        val wgs84Array = epsg3857ToWGS84(epsg3857)
+        return WGS84(wgs84Array[0], wgs84Array[1], 0.0)
+    }
+
+    fun pos2ecef(pos: DoubleArray): DoubleArray {
+        val (lat, lon, hgt) = pos
+        val sinp = sin(lat)
+        val cosp = cos(lat)
+        val sin_l = sin(lon)
+        val cos_l = cos(lon)
+        val e2 = FE_WGS84 * (2.0 - FE_WGS84)
+        val v = RE_WGS84 / sqrt(1.0 - e2 * sinp * sinp)
+
+        return doubleArrayOf(
+            (v + hgt) * cosp * cos_l,
+            (v + hgt) * cosp * sin_l,
+            (v * (1.0 - e2) + hgt) * sinp
+        )
+    }
+
+    fun ecef2enu(pos: DoubleArray, r: DoubleArray): DoubleArray {
+        val E = xyz2enu(pos)
+        return matmul(charArrayOf('N', 'N'), 3, 1, 3, 1.0, E, r, 0.0)
+    }
+
+    fun matmul(
+        tr: CharArray,
+        n: Int,
+        k: Int,
+        m: Int,
+        alpha: Double,
+        A: DoubleArray,
+        B: DoubleArray,
+        beta: Double
+    ): DoubleArray {
+        val f = when {
+            tr[0] == 'N' && tr[1] == 'N' -> 1
+            tr[0] == 'N' && tr[1] == 'T' -> 2
+            tr[0] == 'T' && tr[1] == 'N' -> 3
+            else -> 4
+        }
+        val C = DoubleArray(n * k)
+        for (i in 0 until n) {
+            for (j in 0 until k) {
+                var d = 0.0
+                when (f) {
+                    1 -> for (x in 0 until m) d += A[i + x * n] * B[x + j * m]
+                    2 -> for (x in 0 until m) d += A[i + x * n] * B[j + x * k]
+                    3 -> for (x in 0 until m) d += A[x + i * m] * B[x + j * m]
+                    4 -> for (x in 0 until m) d += A[x + i * m] * B[j + x * k]
+                }
+                C[i + j * n] = alpha * d + beta * C[i + j * n]
+            }
+        }
+        return C
+    }
+
+    fun xyz2enu(pos: DoubleArray): DoubleArray {
+        val (lat, lon) = pos
+        val sinp = sin(lat)
+        val cosp = cos(lat)
+        val sin_l = sin(lon)
+        val cos_l = cos(lon)
+
+        return doubleArrayOf(
+            -sin_l, cos_l, 0.0,
+            -sinp * cos_l, -sinp * sin_l, cosp,
+            cosp * cos_l, cosp * sin_l, sinp
+        )
+    }
+
+    fun ecef2pos(ecef: DoubleArray): DoubleArray {
+        val (x, y, z) = ecef
+        val a = RE_WGS84
+        val b = a * (1 - FE_WGS84)
+        val e2 = (a * a - b * b) / (a * a)
+        val e2p = (a * a - b * b) / (b * b)
+        val r2 = x * x + y * y
+        val r = sqrt(r2)
+        val E2 = a * a - b * b
+        val F = 54 * b * b * z * z
+        val G = r2 + (1 - e2) * z * z - e2 * E2
+        val c = (e2 * e2 * F * r2) / (G * G * G)
+        val s = cbrt(1 + c + sqrt(c * c + 2 * c))
+        val P = F / (3 * (s + 1 / s + 1) * (s + 1 / s + 1) * G * G)
+        val Q = sqrt(1 + 2 * e2 * e2 * P)
+        val r0 = -(P * e2 * r) / (1 + Q) + sqrt(0.5 * a * a * (1 + 1.0 / Q) - P * (1 - e2) * z * z / (Q * (1 + Q)) - 0.5 * P * r2)
+        val U = sqrt((r - e2 * r0) * (r - e2 * r0) + z * z)
+        val V = sqrt((r - e2 * r0) * (r - e2 * r0) + (1 - e2) * z * z)
+        val Z0 = b * b * z / (a * V)
+
+        val lon = atan2(y, x) * RAD2DEG
+        val lat = atan((z + e2p * Z0) / r) * RAD2DEG
+        val hgt = U * (1 - b * b / (a * V))
+        return doubleArrayOf(lon, lat, hgt)
+    }
+
+    data class WGS84(var lon: Double = 0.0, var lat: Double = 0.0, var hgt: Double = 0.0) : Parcelable {
+        constructor(parcel: Parcel) : this(
+            parcel.readDouble(),
+            parcel.readDouble(),
+            parcel.readDouble()
+        )
+
+        constructor(wgs84: DoubleArray) : this(
+            lon = wgs84.getOrElse(0) { 0.0 },
+            lat = wgs84.getOrElse(1) { 0.0 },
+            hgt = wgs84.getOrElse(2) { 0.0 }
+        )
+
+        override fun writeToParcel(parcel: Parcel, flags: Int) {
+            parcel.writeDouble(lon)
+            parcel.writeDouble(lat)
+            parcel.writeDouble(hgt)
+        }
+
+        override fun describeContents(): Int = 0
+
+        companion object CREATOR : Parcelable.Creator<WGS84> {
+            override fun createFromParcel(parcel: Parcel): WGS84 {
+                return WGS84(parcel)
+            }
+
+            override fun newArray(size: Int): Array<WGS84?> {
+                return arrayOfNulls(size)
+            }
+        }
+
+        override fun toString(): String {
+            return "WGS84(lon=$lon, lat=$lat, hgt=$hgt)"
+        }
+    }
+
+    data class EPSG3857(var x: Double = 0.0, var y: Double = 0.0) : Parcelable {
+        constructor(parcel: Parcel) : this(
+            parcel.readDouble(),
+            parcel.readDouble()
+        )
+
+        constructor(epsG3857: DoubleArray) : this(
+            x = epsG3857.getOrElse(0) { 0.0 },
+            y = epsG3857.getOrElse(1) { 0.0 }
+        )
+
+        override fun writeToParcel(parcel: Parcel, flags: Int) {
+            parcel.writeDouble(x)
+            parcel.writeDouble(y)
+        }
+
+        override fun describeContents(): Int = 0
+
+        companion object CREATOR : Parcelable.Creator<EPSG3857> {
+            override fun createFromParcel(parcel: Parcel): EPSG3857 {
+                return EPSG3857(parcel)
+            }
+
+            override fun newArray(size: Int): Array<EPSG3857?> {
+                return arrayOfNulls(size)
+            }
+        }
+
+        override fun toString(): String {
+            return "EPSG3857(x=$x, y=$y)"
+        }
+    }
+
+    data class ENU(var x: Double = 0.0, var y: Double = 0.0, var z: Double = 0.0) : Parcelable {
+        constructor(parcel: Parcel) : this(
+            parcel.readDouble(),
+            parcel.readDouble(),
+            parcel.readDouble()
+        )
+
+        constructor(enu: DoubleArray) : this(
+            x = enu.getOrElse(0) { 0.0 },
+            y = enu.getOrElse(1) { 0.0 },
+            z = enu.getOrElse(2) { 0.0 }
+        )
+
+        override fun writeToParcel(parcel: Parcel, flags: Int) {
+            parcel.writeDouble(x)
+            parcel.writeDouble(y)
+            parcel.writeDouble(z)
+        }
+
+        override fun describeContents(): Int = 0
+
+        companion object CREATOR : Parcelable.Creator<ENU> {
+            override fun createFromParcel(parcel: Parcel): ENU {
+                return ENU(parcel)
+            }
+
+            override fun newArray(size: Int): Array<ENU?> {
+                return arrayOfNulls(size)
+            }
+        }
+
+        override fun toString(): String {
+            return "ENU(x=$x, y=$y, z=$z)"
+        }
+    }
+
+}

icegps-common/src/test/java/com/icegps/common/ExampleUnitTest.kt

@@ -0,0 +1,17 @@
+package com.icegps.common
+
+import org.junit.Test
+
+import org.junit.Assert.*
+
+/**
+ * Example local unit test, which will execute on the development machine (host).
+ *
+ * See [testing documentation](http://d.android.com/tools/testing).
+ */
+class ExampleUnitTest {
+    @Test
+    fun addition_isCorrect() {
+        assertEquals(4, 2 + 2)
+    }
+}

icegps-shared/.gitignore

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

icegps-shared/build.gradle.kts

@@ -0,0 +1,57 @@
+import org.jetbrains.kotlin.gradle.dsl.JvmTarget
+
+plugins {
+    alias(libs.plugins.android.library)
+    alias(libs.plugins.kotlin.android)
+    alias(libs.plugins.kotlin.serialization)
+}
+
+android {
+    namespace = "com.icegps.shared"
+    compileSdk {
+        version = release(36)
+    }
+
+    defaultConfig {
+        minSdk = 28
+
+        testInstrumentationRunner = "androidx.test.runner.AndroidJUnitRunner"
+        consumerProguardFiles("consumer-rules.pro")
+    }
+
+    buildTypes {
+        release {
+            isMinifyEnabled = false
+            proguardFiles(
+                getDefaultProguardFile("proguard-android-optimize.txt"),
+                "proguard-rules.pro"
+            )
+        }
+    }
+    compileOptions {
+        sourceCompatibility = JavaVersion.VERSION_17
+        targetCompatibility = JavaVersion.VERSION_17
+    }
+}
+
+kotlin {
+    compilerOptions.jvmTarget = JvmTarget.JVM_17
+}
+
+dependencies {
+    implementation(libs.core.ktx)
+    implementation(libs.androidx.appcompat)
+    implementation(libs.material)
+
+    api(libs.kotlinx.serialization.json)
+    api(libs.ktor.client.core)
+    api(libs.ktor.client.cio)
+    api(libs.ktor.serialization.kotlinx.json)
+    api(libs.ktor.client.content.negotiation)
+    api(libs.ktor.client.logging)
+    api(project(":math"))
+
+    testImplementation(libs.junit)
+    androidTestImplementation(libs.ext.junit)
+    androidTestImplementation(libs.androidx.espresso.core)
+}

icegps-shared/proguard-rules.pro

@@ -0,0 +1,21 @@
+# Add project specific ProGuard rules here.
+# You can control the set of applied configuration files using the
+# proguardFiles setting in build.gradle.
+#
+# For more details, see
+#   http://developer.android.com/guide/developing/tools/proguard.html
+
+# If your project uses WebView with JS, uncomment the following
+# and specify the fully qualified class name to the JavaScript interface
+# class:
+#-keepclassmembers class fqcn.of.javascript.interface.for.webview {
+#   public *;
+#}
+
+# Uncomment this to preserve the line number information for
+# debugging stack traces.
+#-keepattributes SourceFile,LineNumberTable
+
+# If you keep the line number information, uncomment this to
+# hide the original source file name.
+#-renamesourcefileattribute SourceFile

icegps-shared/src/androidTest/java/com/icegps/shared/ExampleInstrumentedTest.kt

@@ -0,0 +1,24 @@
+package com.icegps.shared
+
+import androidx.test.platform.app.InstrumentationRegistry
+import androidx.test.ext.junit.runners.AndroidJUnit4
+
+import org.junit.Test
+import org.junit.runner.RunWith
+
+import org.junit.Assert.*
+
+/**
+ * Instrumented test, which will execute on an Android device.
+ *
+ * See [testing documentation](http://d.android.com/tools/testing).
+ */
+@RunWith(AndroidJUnit4::class)
+class ExampleInstrumentedTest {
+    @Test
+    fun useAppContext() {
+        // Context of the app under test.
+        val appContext = InstrumentationRegistry.getInstrumentation().targetContext
+        assertEquals("com.icegps.shared.test", appContext.packageName)
+    }
+}

icegps-shared/src/main/AndroidManifest.xml

@@ -0,0 +1,4 @@
+<?xml version="1.0" encoding="utf-8"?>
+<manifest xmlns:android="http://schemas.android.com/apk/res/android">
+
+</manifest>

icegps-shared/src/main/java/com/icegps/shared/SharedHttpClient.kt

@@ -0,0 +1,42 @@
+package com.icegps.shared
+
+import io.ktor.client.HttpClient
+import io.ktor.client.engine.cio.CIO
+import io.ktor.client.plugins.HttpTimeout
+import io.ktor.client.plugins.contentnegotiation.ContentNegotiation
+import io.ktor.client.plugins.logging.LogLevel
+import io.ktor.client.plugins.logging.Logger
+import io.ktor.client.plugins.logging.Logging
+import io.ktor.client.plugins.logging.SIMPLE
+import io.ktor.http.ContentType
+import io.ktor.serialization.kotlinx.json.json
+import kotlinx.serialization.json.Json
+
+/**
+ * @author tabidachinokaze
+ * @date 2025/11/20
+ */
+@Suppress("FunctionName")
+fun SharedHttpClient(json: Json): HttpClient {
+    return HttpClient(CIO) {
+        install(ContentNegotiation) {
+            json(
+                json = json,
+                contentType = ContentType.Text.Html
+            )
+            json(
+                json = json,
+                contentType = ContentType.Application.Json
+            )
+        }
+        install(Logging) {
+            this.level = LogLevel.ALL
+            this.logger = Logger.SIMPLE
+        }
+        install(HttpTimeout) {
+            requestTimeoutMillis = 1000 * 60 * 10
+            connectTimeoutMillis = 1000 * 60 * 5
+            socketTimeoutMillis = 1000 * 60 * 10
+        }
+    }
+}

icegps-shared/src/main/java/com/icegps/shared/SharedJson.kt

@@ -0,0 +1,17 @@
+package com.icegps.shared
+
+import kotlinx.serialization.json.Json
+import kotlinx.serialization.modules.SerializersModule
+
+/**
+ * @author tabidachinokaze
+ * @date 2025/11/20
+ */
+@Suppress("FunctionName")
+fun SharedJson(): Json {
+    return Json {
+        ignoreUnknownKeys = true
+        serializersModule = SerializersModule {
+        }
+    }
+}

icegps-shared/src/main/java/com/icegps/shared/api/LookupResponse.kt

@@ -0,0 +1,21 @@
+package com.icegps.shared.api
+
+import com.icegps.shared.model.IGeoPoint
+import kotlinx.serialization.SerialName
+import kotlinx.serialization.Serializable
+
+@Serializable
+data class LookupResponse(
+    @SerialName("results")
+    val results: List<Result>
+) {
+    @Serializable
+    data class Result(
+        @SerialName("longitude")
+        override val longitude: Double,
+        @SerialName("latitude")
+        override val latitude: Double,
+        @SerialName("elevation")
+        override val altitude: Double,
+    ) : IGeoPoint
+}

icegps-shared/src/main/java/com/icegps/shared/api/OpenElevationApi.kt

@@ -0,0 +1,35 @@
+package com.icegps.shared.api
+
+import com.icegps.shared.model.IGeoPoint
+import io.ktor.client.HttpClient
+import io.ktor.client.call.body
+import io.ktor.client.request.get
+import io.ktor.client.request.parameter
+import io.ktor.http.appendPathSegments
+
+/**
+ * @author tabidachinokaze
+ * @date 2025/11/20
+ */
+interface OpenElevationApi {
+    suspend fun lookup(values: List<IGeoPoint>): List<IGeoPoint>
+}
+
+class OpenElevation(
+    private val client: HttpClient
+) : OpenElevationApi {
+    private val baseUrl: String = "https://api.open-elevation.com/api/v1/"
+
+    // curl 'https://api.open-elevation.com/api/v1/lookup?locations=10,10|20,20|41.161758,-8.583933'
+    override suspend fun lookup(values: List<IGeoPoint>): List<IGeoPoint> {
+        val response = client.get(baseUrl) {
+            url {
+                appendPathSegments("lookup")
+                parameter(
+                    "locations",
+                    values.joinToString("|") { "${it.latitude},${it.longitude}" })
+            }
+        }
+        return response.body<LookupResponse>().results
+    }
+}

icegps-shared/src/main/java/com/icegps/shared/ktx/Any.kt

@@ -0,0 +1,7 @@
+package com.icegps.shared.ktx
+
+/**
+ * @author tabidachinokaze
+ * @date 2025/11/22
+ */
+val Any.TAG: String get() = this::class.java.simpleName

icegps-shared/src/main/java/com/icegps/shared/model/GeoPoint.kt

@@ -0,0 +1,11 @@
+package com.icegps.shared.model
+
+/**
+ * @author tabidachinokaze
+ * @date 2025/11/22
+ */
+data class GeoPoint(
+    override val longitude: Double,
+    override val latitude: Double,
+    override val altitude: Double
+) : IGeoPoint

icegps-shared/src/main/java/com/icegps/shared/model/IGeoPoint.kt

@@ -0,0 +1,11 @@
+package com.icegps.shared.model
+
+/**
+ * @author tabidachinokaze
+ * @date 2025/11/22
+ */
+interface IGeoPoint {
+    val longitude: Double
+    val latitude: Double
+    val altitude: Double
+}

icegps-shared/src/test/java/com/icegps/shared/ExampleUnitTest.kt

@@ -0,0 +1,17 @@
+package com.icegps.shared
+
+import org.junit.Test
+
+import org.junit.Assert.*
+
+/**
+ * Example local unit test, which will execute on the development machine (host).
+ *
+ * See [testing documentation](http://d.android.com/tools/testing).
+ */
+class ExampleUnitTest {
+    @Test
+    fun addition_isCorrect() {
+        assertEquals(4, 2 + 2)
+    }
+}

math/.gitignore

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

math/build.gradle.kts

@@ -0,0 +1,13 @@
+import org.jetbrains.kotlin.gradle.dsl.JvmTarget
+
+plugins {
+    alias(libs.plugins.kotlin.jvm)
+}
+
+java {
+    sourceCompatibility = JavaVersion.VERSION_17
+    targetCompatibility = JavaVersion.VERSION_17
+}
+kotlin {
+    compilerOptions.jvmTarget = JvmTarget.JVM_17
+}

math/src/main/java/com/icegps/io/util/NumberExt.kt

@@ -0,0 +1,75 @@
+@file:Suppress("NOTHING_TO_INLINE")
+
+package com.icegps.io.util
+
+import kotlin.math.*
+
+//private fun Double.normalizeZero(): Double = if (this.isAlmostZero()) 0.0 else this
+private val MINUS_ZERO_D = -0.0
+private fun Double.normalizeZero(): Double = if (this == MINUS_ZERO_D) 0.0 else this
+
+fun Double.toStringDecimal(decimalPlaces: Int, skipTrailingZeros: Boolean = false): String {
+    if (this.isNanOrInfinite()) return this.toString()
+
+    //val bits = this.toRawBits()
+    //val sign = (bits ushr 63) != 0L
+    //val exponent = (bits ushr 52) and 0b11111111111
+    //val fraction = bits and ((1L shl 52) - 1L)
+
+	val res = this.roundDecimalPlaces(decimalPlaces).normalizeZero().toString()
+
+	val eup = res.indexOf('E')
+	val elo = res.indexOf('e')
+    val eIndex = if (eup >= 0) eup else elo
+	val rez = if (eIndex >= 0) {
+        val base = res.substring(0, eIndex)
+        val exp = res.substring(eIndex + 1).toInt()
+        val rbase = if (base.contains(".")) base else "$base.0"
+        val zeros = "0".repeat(exp.absoluteValue + 2)
+		val part = if (exp > 0) "$rbase$zeros" else "$zeros$rbase"
+		val pointIndex2 = part.indexOf(".")
+        val pointIndex = if (pointIndex2 < 0) part.length else pointIndex2
+		val outIndex = pointIndex + exp
+		val part2 = part.replace(".", "")
+        buildString {
+            if ((0 until outIndex).all { part2[it] == '0' }) {
+                append('0')
+            } else {
+                append(part2, 0, outIndex)
+            }
+            append('.')
+            append(part2, outIndex, part2.length)
+        }
+	} else {
+		res
+	}
+
+    val pointIndex = rez.indexOf('.')
+	val integral = if (pointIndex >= 0) rez.substring(0, pointIndex) else rez
+    if (decimalPlaces == 0) return integral
+
+	val decimal = if (pointIndex >= 0) rez.substring(pointIndex + 1).trimEnd('0') else ""
+    return buildString(2 + integral.length + decimalPlaces) {
+        append(integral)
+        if (decimal.isNotEmpty() || !skipTrailingZeros) {
+            val decimalCount = min(decimal.length, decimalPlaces)
+            val allZeros = (0 until decimalCount).all { decimal[it] == '0' }
+            if (!skipTrailingZeros || !allZeros) {
+                append('.')
+                append(decimal, 0, decimalCount)
+                if (!skipTrailingZeros) repeat(decimalPlaces - decimalCount) { append('0') }
+            }
+        }
+    }
+}
+
+fun Float.toStringDecimal(decimalPlaces: Int, skipTrailingZeros: Boolean = false): String = this.toDouble().toStringDecimal(decimalPlaces, skipTrailingZeros)
+
+private fun Double.roundDecimalPlaces(places: Int): Double {
+    if (places < 0) return this
+    val placesFactor: Double = 10.0.pow(places.toDouble())
+    return round(this * placesFactor) / placesFactor
+}
+
+private fun Double.isNanOrInfinite() = this.isNaN() || this.isInfinite()
+private fun Float.isNanOrInfinite() = this.isNaN() || this.isInfinite()

math/src/main/java/com/icegps/io/util/NumberParser.kt

@@ -0,0 +1,78 @@
+package com.icegps.io.util
+
+import kotlin.math.*
+
+object NumberParser {
+	const val END = '\u0000'
+
+	fun parseInt(str: String, start: Int = 0, end: Int = str.length, radix: Int = 10): Int {
+		var n = start
+		return parseInt(radix) { if (n >= end) END else str[n++] }
+	}
+
+	fun parseDouble(str: String, start: Int = 0, end: Int = str.length): Double {
+		var n = start
+		return parseDouble { if (n >= end) END else str[n++] }
+	}
+
+	inline fun parseInt(radix: Int = 10, gen: (Int) -> Char): Int {
+		var positive = true
+		var out = 0
+		var n = 0
+		while (true) {
+			val c = gen(n++)
+			if (c == END) break
+			if (c == '-' || c == '+') {
+				positive = (c == '+')
+			} else {
+				val value = c.ctypeAsInt()
+				if (value < 0) break
+				out *= radix
+				out += value
+			}
+		}
+		return if (positive) out else -out
+	}
+
+	inline fun parseDouble(gen: (Int) -> Char): Double {
+		var out = 0.0
+		var frac = 1.0
+		var pointSeen = false
+		var eSeen = false
+		var negate = false
+		var negateExponent = false
+		var exponent = 0
+		var n = 0
+		while (true) {
+			val c = gen(n++)
+			if (c == END) break
+			when (c) {
+				'e', 'E' -> eSeen = true
+				'-' -> {
+					if (eSeen) negateExponent = true else negate = true
+				}
+				'.' -> pointSeen = true
+				else -> {
+					if (eSeen) {
+						exponent *= 10
+						exponent += c.ctypeAsInt()
+					} else {
+						if (pointSeen) frac /= 10
+						out *= 10
+						out += c.ctypeAsInt()
+					}
+				}
+			}
+		}
+		val res = (out * frac) * 10.0.pow(if (negateExponent) -exponent else exponent)
+		return if (negate) -res else res
+	}
+}
+
+@Suppress("ConvertTwoComparisonsToRangeCheck") // @TODO: Kotlin-Native doesn't optimize ranges
+@PublishedApi internal fun Char.ctypeAsInt(): Int = when {
+	this >= '0' && this <= '9' -> this - '0'
+	this >= 'a' && this <= 'z' -> this - 'a' + 10
+	this >= 'A' && this <= 'Z' -> this - 'A' + 10
+	else -> -1
+}

math/src/main/java/com/icegps/math/Alignment.kt

@@ -0,0 +1,47 @@
+package com.icegps.math
+
+import kotlin.math.absoluteValue
+
+////////////////////
+////////////////////
+
+/** Returns the next value of [this] that is multiple of [align]. If [this] is already multiple, returns itself. */
+public fun Int.nextAlignedTo(align: Int): Int = if (this.isAlignedTo(align)) this else (((this / align) + 1) * align)
+/** Returns the next value of [this] that is multiple of [align]. If [this] is already multiple, returns itself. */
+public fun Long.nextAlignedTo(align: Long): Long = if (this.isAlignedTo(align)) this else (((this / align) + 1) * align)
+/** Returns the next value of [this] that is multiple of [align]. If [this] is already multiple, returns itself. */
+public fun Float.nextAlignedTo(align: Float): Float = if (this.isAlignedTo(align)) this else (((this / align).toInt() + 1) * align)
+/** Returns the next value of [this] that is multiple of [align]. If [this] is already multiple, returns itself. */
+public fun Double.nextAlignedTo(align: Double): Double = if (this.isAlignedTo(align)) this else (((this / align).toInt() + 1) * align)
+
+/** Returns the previous value of [this] that is multiple of [align]. If [this] is already multiple, returns itself. */
+public fun Int.prevAlignedTo(align: Int): Int = if (this.isAlignedTo(align)) this else nextAlignedTo(align) - align
+/** Returns the previous value of [this] that is multiple of [align]. If [this] is already multiple, returns itself. */
+public fun Long.prevAlignedTo(align: Long): Long = if (this.isAlignedTo(align)) this else nextAlignedTo(align) - align
+/** Returns the previous value of [this] that is multiple of [align]. If [this] is already multiple, returns itself. */
+public fun Float.prevAlignedTo(align: Float): Float = if (this.isAlignedTo(align)) this else nextAlignedTo(align) - align
+/** Returns the previous value of [this] that is multiple of [align]. If [this] is already multiple, returns itself. */
+public fun Double.prevAlignedTo(align: Double): Double = if (this.isAlignedTo(align)) this else nextAlignedTo(align) - align
+
+/** Returns whether [this] is multiple of [alignment] */
+public fun Int.isAlignedTo(alignment: Int): Boolean = alignment == 0 || (this % alignment) == 0
+/** Returns whether [this] is multiple of [alignment] */
+public fun Long.isAlignedTo(alignment: Long): Boolean = alignment == 0L || (this % alignment) == 0L
+/** Returns whether [this] is multiple of [alignment] */
+public fun Float.isAlignedTo(alignment: Float): Boolean = alignment == 0f || (this % alignment) == 0f
+/** Returns whether [this] is multiple of [alignment] */
+public fun Double.isAlignedTo(alignment: Double): Boolean = alignment == 0.0 || (this % alignment) == 0.0
+
+/** Returns the previous or next value of [this] that is multiple of [align]. If [this] is already multiple, returns itself. */
+public fun Float.nearestAlignedTo(align: Float): Float {
+    val prev = this.prevAlignedTo(align)
+    val next = this.nextAlignedTo(align)
+    return if ((this - prev).absoluteValue < (this - next).absoluteValue) prev else next
+}
+/** Returns the previous or next value of [this] that is multiple of [align]. If [this] is already multiple, returns itself. */
+public fun Double.nearestAlignedTo(align: Double): Double {
+    val prev = this.prevAlignedTo(align)
+    val next = this.nextAlignedTo(align)
+    return if ((this - prev).absoluteValue < (this - next).absoluteValue) prev else next
+}
+

math/src/main/java/com/icegps/math/BooleanConversion.kt

@@ -0,0 +1,9 @@
+package com.icegps.math
+
+////////////////////
+////////////////////
+
+/** Converts this [Boolean] into integer: 1 for true, 0 for false */
+inline fun Boolean.toInt(): Int = if (this) 1 else 0
+inline fun Boolean.toByte(): Byte = if (this) 1 else 0
+inline fun Byte.toBoolean(): Boolean = this.toInt() != 0

math/src/main/java/com/icegps/math/Clamp.kt

@@ -0,0 +1,38 @@
+package com.icegps.math
+
+/** Clamps [this] value into the range [min] and [max] */
+fun Int.clamp(min: Int, max: Int): Int = if (this < min) min else if (this > max) max else this
+/** Clamps [this] value into the range [min] and [max] */
+fun Long.clamp(min: Long, max: Long): Long = if (this < min) min else if (this > max) max else this
+/** Clamps [this] value into the range [min] and [max] */
+fun Double.clamp(min: Double, max: Double): Double = if (this < min) min else if (this > max) max else this
+/** Clamps [this] value into the range [min] and [max] */
+fun Float.clamp(min: Float, max: Float): Float = if ((this < min)) min else if ((this > max)) max else this
+
+/** Clamps [this] value into the range 0 and 1 */
+fun Double.clamp01(): Double = clamp(0.0, 1.0)
+/** Clamps [this] value into the range 0 and 1 */
+fun Float.clamp01(): Float = clamp(0f, 1f)
+
+/** Clamps [this] [Long] value into the range [min] and [max] converting it into [Int]. The default parameters will cover the whole range of values. */
+fun Long.toIntClamp(min: Int = Int.MIN_VALUE, max: Int = Int.MAX_VALUE): Int {
+    if (this < min) return min
+    if (this > max) return max
+    return this.toInt()
+}
+
+/** Clamps [this] [Long] value into the range [min] and [max] converting it into [Int] (where [min] must be zero or positive). The default parameters will cover the whole range of positive and zero values. */
+fun Long.toUintClamp(min: Int = 0, max: Int = Int.MAX_VALUE): Int = this.toIntClamp(min, max)
+
+/** Clamps the integer value in the 0..255 range */
+fun Int.clampUByte(): Int {
+    val n = this and -(if (this >= 0) 1 else 0)
+    return (n or (0xFF - n shr 31)) and 0xFF
+}
+fun Int.clampUShort(): Int {
+    val n = this and -(if (this >= 0) 1 else 0)
+    return (n or (0xFFFF - n shr 31)) and 0xFFFF
+}
+
+fun Int.toShortClamped(): Short = this.clamp(Short.MIN_VALUE.toInt(), Short.MAX_VALUE.toInt()).toShort()
+fun Int.toByteClamped(): Byte = this.clamp(Byte.MIN_VALUE.toInt(), Byte.MAX_VALUE.toInt()).toByte()

math/src/main/java/com/icegps/math/ConvertRange.kt

@@ -0,0 +1,24 @@
+package com.icegps.math
+
+////////////////////
+////////////////////
+
+/** Converts this value considering it was in the range [srcMin]..[srcMax] into [dstMin]..[dstMax], if the value is not inside the range the output value will be outside the destination range */
+fun Float.convertRange(srcMin: Float, srcMax: Float, dstMin: Float, dstMax: Float): Float = (dstMin + (dstMax - dstMin) * ((this - srcMin) / (srcMax - srcMin)))
+/** Converts this value considering it was in the range [srcMin]..[srcMax] into [dstMin]..[dstMax], if the value is not inside the range the output value will be outside the destination range */
+fun Double.convertRange(srcMin: Double, srcMax: Double, dstMin: Double, dstMax: Double): Double = (dstMin + (dstMax - dstMin) * ((this - srcMin) / (srcMax - srcMin)))
+//fun Double.convertRange(minSrc: Double, maxSrc: Double, minDst: Double, maxDst: Double): Double = (((this - minSrc) / (maxSrc - minSrc)) * (maxDst - minDst)) + minDst
+/** Converts this value considering it was in the range [srcMin]..[srcMax] into [dstMin]..[dstMax], if the value is not inside the range the output value will be outside the destination range */
+fun Int.convertRange(srcMin: Int, srcMax: Int, dstMin: Int, dstMax: Int): Int = (dstMin + (dstMax - dstMin) * ((this - srcMin).toDouble() / (srcMax - srcMin).toDouble())).toInt()
+/** Converts this value considering it was in the range [srcMin]..[srcMax] into [dstMin]..[dstMax], if the value is not inside the range the output value will be outside the destination range */
+fun Long.convertRange(srcMin: Long, srcMax: Long, dstMin: Long, dstMax: Long): Long = (dstMin + (dstMax - dstMin) * ((this - srcMin).toDouble() / (srcMax - srcMin).toDouble())).toLong()
+
+/** Converts this value considering it was in the range [srcMin]..[srcMax] into [dstMin]..[dstMax], if the value is not inside the range the output value will be clamped to the nearest bound */
+fun Float.convertRangeClamped(srcMin: Float, srcMax: Float, dstMin: Float, dstMax: Float): Float = convertRange(srcMin, srcMax, dstMin, dstMax).clamp(dstMin, dstMax)
+/** Converts this value considering it was in the range [srcMin]..[srcMax] into [dstMin]..[dstMax], if the value is not inside the range the output value will be clamped to the nearest bound */
+fun Double.convertRangeClamped(srcMin: Double, srcMax: Double, dstMin: Double, dstMax: Double): Double = convertRange(srcMin, srcMax, dstMin, dstMax).clamp(dstMin, dstMax)
+/** Converts this value considering it was in the range [srcMin]..[srcMax] into [dstMin]..[dstMax], if the value is not inside the range the output value will be clamped to the nearest bound */
+fun Int.convertRangeClamped(srcMin: Int, srcMax: Int, dstMin: Int, dstMax: Int): Int = convertRange(srcMin, srcMax, dstMin, dstMax).clamp(dstMin, dstMax)
+/** Converts this value considering it was in the range [srcMin]..[srcMax] into [dstMin]..[dstMax], if the value is not inside the range the output value will be clamped to the nearest bound */
+fun Long.convertRangeClamped(srcMin: Long, srcMax: Long, dstMin: Long, dstMax: Long): Long = convertRange(srcMin, srcMax, dstMin, dstMax).clamp(dstMin, dstMax)
+

math/src/main/java/com/icegps/math/Division.kt

@@ -0,0 +1,19 @@
+package com.icegps.math
+
+import kotlin.math.roundToInt
+
+////////////////////
+////////////////////
+
+/** Divides [this] into [that] rounding to the floor */
+public infix fun Int.divFloor(that: Int): Int = this / that
+/** Divides [this] into [that] rounding to the ceil */
+public infix fun Int.divCeil(that: Int): Int = if (this % that != 0) (this / that) + 1 else (this / that)
+/** Divides [this] into [that] rounding to the round */
+public infix fun Int.divRound(that: Int): Int = (this.toDouble() / that.toDouble()).roundToInt()
+
+public infix fun Long.divCeil(other: Long): Long {
+    val res = this / other
+    if (this % other != 0L) return res + 1
+    return res
+}

math/src/main/java/com/icegps/math/Fract.kt

@@ -0,0 +1,4 @@
+package com.icegps.math
+
+public inline fun fract(value: Float): Float = value - value.toIntFloor()
+public inline fun fract(value: Double): Double = value - value.toIntFloor()

math/src/main/java/com/icegps/math/ILog.kt

@@ -0,0 +1,9 @@
+package com.icegps.math
+
+////////////////////
+////////////////////
+
+/** Performs a fast integral logarithmic of base two */
+fun ilog2(v: Int): Int = if (v == 0) (-1) else (31 - v.countLeadingZeroBits())
+// fun ilog2(v: Int): Int = kotlin.math.log2(v.toDouble()).toInt()
+fun ilog2Ceil(v: Int): Int = kotlin.math.ceil(kotlin.math.log2(v.toDouble())).toInt()

math/src/main/java/com/icegps/math/IsAlmostEquals.kt

@@ -0,0 +1,14 @@
+package com.icegps.math
+
+import kotlin.math.*
+
+interface IsAlmostEquals<T> {
+    fun isAlmostEquals(other: T, epsilon: Double = 0.000001): Boolean
+}
+
+interface IsAlmostEqualsF<T> {
+    fun isAlmostEquals(other: T, epsilon: Float = 0.0001f): Boolean
+}
+
+fun Float.isAlmostEquals(other: Float, epsilon: Float = 0.000001f): Boolean = (this - other).absoluteValue < epsilon
+fun Double.isAlmostEquals(other: Double, epsilon: Double = 0.000001): Boolean = (this - other).absoluteValue < epsilon

math/src/main/java/com/icegps/math/IsAlmostZero.kt

@@ -0,0 +1,4 @@
+package com.icegps.math
+
+fun Double.isAlmostZero(): Boolean = kotlin.math.abs(this) <= 1e-19
+fun Float.isAlmostZero(): Boolean = kotlin.math.abs(this) <= 1e-6

math/src/main/java/com/icegps/math/IsEven.kt

@@ -0,0 +1,9 @@
+package com.icegps.math
+
+////////////////////
+////////////////////
+
+/** Checks if [this] is odd (not multiple of two) */
+val Int.isOdd: Boolean get() = (this % 2) == 1
+/** Checks if [this] is even (multiple of two) */
+val Int.isEven: Boolean get() = (this % 2) == 0

math/src/main/java/com/icegps/math/IsNanOrInfinite.kt

@@ -0,0 +1,11 @@
+package com.icegps.math
+
+///** Check if [this] floating point value is not a number or infinite */
+//public fun Float.isNanOrInfinite(): Boolean = this.isNaN() || this.isInfinite()
+///** Check if [this] floating point value is not a number or infinite */
+//public fun Double.isNanOrInfinite(): Boolean = this.isNaN() || this.isInfinite()
+
+
+fun Double.isNanOrInfinite() = this.isNaN() || this.isInfinite()
+
+fun Float.isNanOrInfinite() = this.isNaN() || this.isInfinite()

math/src/main/java/com/icegps/math/Math.kt

@@ -0,0 +1,133 @@
+package com.icegps.math
+
+import kotlin.math.*
+
+const val PIF = PI.toFloat()
+const val PI2F = (PI * 2).toFloat()
+
+fun Double.betweenInclusive(min: Double, max: Double): Boolean = (this >= min) && (this <= max)
+
+fun almostEquals(a: Float, b: Float) = almostZero(a - b)
+fun almostZero(a: Float) = abs(a) <= 0.0000001
+
+fun almostEquals(a: Double, b: Double) = almostZero(a - b)
+fun almostZero(a: Double) = abs(a) <= 0.0000001
+
+fun isEquivalent(a: Double, b: Double, epsilon: Double = 0.0001): Boolean = (a - epsilon < b) && (a + epsilon > b)
+
+fun Double.smoothstep(edge0: Double, edge1: Double): Double {
+    if (this < edge0) return 0.0
+    if (this >= edge1) return 1.0
+    val v = ((this - edge0) / (edge1 - edge0))//.clamp(0.0, 1.0)
+    return v * v * (3 - 2 * v)
+}
+
+fun log(v: Int, base: Int): Int = log(v.toDouble(), base.toDouble()).toInt()
+fun ln(v: Int): Int = ln(v.toDouble()).toInt()
+fun log2(v: Int): Int = log(v.toDouble(), 2.0).toInt()
+fun log10(v: Int): Int = log(v.toDouble(), 10.0).toInt()
+
+@Deprecated("", ReplaceWith("v.squared()"))
+fun sq(v: Int): Int = v.squared()
+@Deprecated("", ReplaceWith("v.squared()"))
+fun sq(v: Float): Float = v.squared()
+@Deprecated("", ReplaceWith("v.squared()"))
+fun sq(v: Double): Double = v.squared()
+
+/** Signs of the value. Zero will be converted into -1 */
+val Int.signM1: Int get() = signNonZeroM1(this)
+/** Signs of the value. Zero will be converted into -1 */
+val Float.signM1: Float get() = signNonZeroM1(this).toFloat()
+/** Signs of the value. Zero will be converted into -1 */
+val Double.signM1: Double get() = signNonZeroM1(this).toDouble()
+
+/** Signs of the value. Zero will be converted into +1 */
+val Int.signP1: Int get() = signNonZeroP1(this)
+/** Signs of the value. Zero will be converted into +1 */
+val Float.signP1: Float get() = signNonZeroP1(this).toFloat()
+/** Signs of the value. Zero will be converted into +1 */
+val Double.signP1: Double get() = signNonZeroP1(this).toDouble()
+
+/** Signs of the value. Zero will be converted into -1 */
+fun signNonZeroM1(x: Int): Int = if (x <= 0) -1 else +1
+/** Signs of the value. Zero will be converted into -1 */
+fun signNonZeroM1(x: Float): Int = if (x <= 0) -1 else +1
+/** Signs of the value. Zero will be converted into -1 */
+fun signNonZeroM1(x: Double): Int = if (x <= 0) -1 else +1
+
+
+/** Signs of the value. Zero will be converted into +1 */
+fun signNonZeroP1(x: Int): Int = if (x >= 0) +1 else -1
+/** Signs of the value. Zero will be converted into +1 */
+fun signNonZeroP1(x: Float): Int = if (x >= 0) +1 else -1
+/** Signs of the value. Zero will be converted into +1 */
+fun signNonZeroP1(x: Double): Int = if (x >= 0) +1 else -1
+
+fun Float.normalizeAlmostZero() = if (this.isAlmostZero()) 0f else this
+
+fun Double.closestMultipleOf(multiple: Double): Double {
+    val prev = prevMultipleOf(multiple)
+    val next = nextMultipleOf(multiple)
+    return if ((this - prev).absoluteValue < (this - next).absoluteValue) prev else next
+}
+fun Int.closestMultipleOf(multiple: Int): Int {
+    val prev = prevMultipleOf(multiple)
+    val next = nextMultipleOf(multiple)
+    return if ((this - prev).absoluteValue < (this - next).absoluteValue) prev else next
+}
+fun Long.closestMultipleOf(multiple: Long): Long {
+    val prev = prevMultipleOf(multiple)
+    val next = nextMultipleOf(multiple)
+    return if ((this - prev).absoluteValue < (this - next).absoluteValue) prev else next
+}
+
+fun Double.nextMultipleOf(multiple: Double) = if (this.isMultipleOf(multiple)) this else (((this / multiple) + 1) * multiple)
+fun Int.nextMultipleOf(multiple: Int) = if (this.isMultipleOf(multiple)) this else (((this / multiple) + 1) * multiple)
+fun Long.nextMultipleOf(multiple: Long) = if (this.isMultipleOf(multiple)) this else (((this / multiple) + 1) * multiple)
+
+fun Double.prevMultipleOf(multiple: Double) = if (this.isMultipleOf(multiple)) this else nextMultipleOf(multiple) - multiple
+fun Int.prevMultipleOf(multiple: Int) = if (this.isMultipleOf(multiple)) this else nextMultipleOf(multiple) - multiple
+fun Long.prevMultipleOf(multiple: Long) = if (this.isMultipleOf(multiple)) this else nextMultipleOf(multiple) - multiple
+
+fun Double.isMultipleOf(multiple: Double) = multiple.isAlmostZero() || (this % multiple).isAlmostZero()
+fun Int.isMultipleOf(multiple: Int) = multiple == 0 || (this % multiple) == 0
+fun Long.isMultipleOf(multiple: Long) = multiple == 0L || (this % multiple) == 0L
+
+fun Double.squared(): Double = this * this
+fun Float.squared(): Float = this * this
+fun Int.squared(): Int = this * this
+
+fun min(a: Int, b: Int, c: Int) = min(min(a, b), c)
+fun min(a: Float, b: Float, c: Float) = min(min(a, b), c)
+fun min(a: Double, b: Double, c: Double) = min(min(a, b), c)
+
+fun min(a: Int, b: Int, c: Int, d: Int) = min(min(min(a, b), c), d)
+fun min(a: Float, b: Float, c: Float, d: Float) = min(min(min(a, b), c), d)
+fun min(a: Double, b: Double, c: Double, d: Double) = min(min(min(a, b), c), d)
+
+fun min(a: Int, b: Int, c: Int, d: Int, e: Int) = min(min(min(min(a, b), c), d), e)
+fun min(a: Float, b: Float, c: Float, d: Float, e: Float) = min(min(min(min(a, b), c), d), e)
+fun min(a: Double, b: Double, c: Double, d: Double, e: Double) = min(min(min(min(a, b), c), d), e)
+
+fun max(a: Int, b: Int, c: Int) = max(max(a, b), c)
+fun max(a: Float, b: Float, c: Float) = max(max(a, b), c)
+fun max(a: Double, b: Double, c: Double) = max(max(a, b), c)
+
+fun max(a: Int, b: Int, c: Int, d: Int) = max(max(max(a, b), c), d)
+fun max(a: Float, b: Float, c: Float, d: Float) = max(max(max(a, b), c), d)
+fun max(a: Double, b: Double, c: Double, d: Double) = max(max(max(a, b), c), d)
+
+fun max(a: Int, b: Int, c: Int, d: Int, e: Int) = max(max(max(max(a, b), c), d), e)
+fun max(a: Float, b: Float, c: Float, d: Float, e: Float) = max(max(max(max(a, b), c), d), e)
+fun max(a: Double, b: Double, c: Double, d: Double, e: Double) = max(max(max(max(a, b), c), d), e)
+
+////////////////////
+////////////////////
+
+
+// @TODO: Optimize this
+fun Int.numberOfDigits(radix: Int = 10): Int = radix.toString(radix).length
+fun Long.numberOfDigits(radix: Int = 10): Int = radix.toString(radix).length
+
+fun Int.cycle(min: Int, max: Int): Int = ((this - min) umod (max - min + 1)) + min
+fun Int.cycleSteps(min: Int, max: Int): Int = (this - min) / (max - min + 1)

math/src/main/java/com/icegps/math/NormalizeZero.kt

@@ -0,0 +1,7 @@
+package com.icegps.math
+
+//fun Double.normalizeZero(): Double = if (this.isAlmostZero()) 0.0 else this
+private val MINUS_ZERO_D = -0.0
+private val MINUS_ZERO_F = -0.0f
+fun Double.normalizeZero(): Double = if (this == MINUS_ZERO_D) 0.0 else this
+fun Float.normalizeZero(): Float = if (this == MINUS_ZERO_F) 0f else this

math/src/main/java/com/icegps/math/PowerOfTwo.kt

@@ -0,0 +1,22 @@
+package com.icegps.math
+
+
+/** Returns the next power of two of [this] */
+val Int.nextPowerOfTwo: Int get() {
+    var v = this
+    v--
+    v = v or (v shr 1)
+    v = v or (v shr 2)
+    v = v or (v shr 4)
+    v = v or (v shr 8)
+    v = v or (v shr 16)
+    v++
+    return v
+}
+/** Checks if [this] value is power of two */
+val Int.isPowerOfTwo: Boolean get() = this.nextPowerOfTwo == this
+
+/** Returns the previous power of two of [this] */
+val Int.prevPowerOfTwo: Int get() = if (isPowerOfTwo) this else (nextPowerOfTwo ushr 1)
+
+

math/src/main/java/com/icegps/math/RoundDecimalPlaces.kt

@@ -0,0 +1,16 @@
+package com.icegps.math
+
+import kotlin.math.*
+
+fun Float.roundDecimalPlaces(places: Int): Float {
+    if (places < 0) return this
+    val placesFactor: Float = 10f.pow(places.toFloat())
+    return round(this * placesFactor) / placesFactor
+}
+
+fun Double.roundDecimalPlaces(places: Int): Double {
+    if (places < 0) return this
+    val placesFactor: Double = 10.0.pow(places.toDouble())
+    return round(this * placesFactor) / placesFactor
+}
+

math/src/main/java/com/icegps/math/ToIntegerConverters.kt

@@ -0,0 +1,29 @@
+package com.icegps.math
+
+import kotlin.math.*
+
+////////////////////
+////////////////////
+
+/** Converts [this] into [Int] rounding to the ceiling */
+fun Float.toIntCeil(): Int = ceil(this).toInt()
+/** Converts [this] into [Int] rounding to the ceiling */
+fun Double.toIntCeil(): Int = ceil(this).toInt()
+
+/** Converts [this] into [Int] rounding to the nearest */
+fun Float.toIntRound(): Int = round(this).toInt()
+/** Converts [this] into [Int] rounding to the nearest */
+fun Double.toIntRound(): Int = round(this).toInt()
+
+/** Converts [this] into [Int] rounding to the nearest */
+fun Float.toLongRound(): Long = round(this).toLong()
+/** Converts [this] into [Int] rounding to the nearest */
+fun Double.toLongRound(): Long = round(this).toLong()
+
+/** Convert this [Long] into an [Int] but throws an [IllegalArgumentException] in the case that operation would produce an overflow */
+fun Long.toIntSafe(): Int = if (this in Int.MIN_VALUE.toLong()..Int.MAX_VALUE.toLong()) this.toInt() else throw IllegalArgumentException("Long doesn't fit Integer")
+
+/** Converts [this] into [Int] rounding to the floor */
+fun Float.toIntFloor(): Int = floor(this).toInt()
+/** Converts [this] into [Int] rounding to the floor */
+fun Double.toIntFloor(): Int = floor(this).toInt()

math/src/main/java/com/icegps/math/Umod.kt

@@ -0,0 +1,35 @@
+package com.icegps.math
+
+private val MINUS_ZERO_F = -0.0f
+
+////////////////////
+////////////////////
+
+/** Performs the unsigned modulo between [this] and [other] (negative values would wrap) */
+public infix fun Int.umod(other: Int): Int {
+    val rm = this % other
+    val remainder = if (rm == -0) 0 else rm
+    return when {
+        remainder < 0 -> remainder + other
+        else -> remainder
+    }
+}
+
+/** Performs the unsigned modulo between [this] and [other] (negative values would wrap) */
+public infix fun Double.umod(other: Double): Double {
+    val rm = this % other
+    val remainder = if (rm == -0.0) 0.0 else rm
+    return when {
+        remainder < 0.0 -> remainder + other
+        else -> remainder
+    }
+}
+
+public infix fun Float.umod(other: Float): Float {
+    val rm = this % other
+    val remainder = if (rm == MINUS_ZERO_F) 0f else rm
+    return when {
+        remainder < 0f -> remainder + other
+        else -> remainder
+    }
+}

math/src/main/java/com/icegps/math/Unsigned.kt

@@ -0,0 +1,13 @@
+package com.icegps.math
+
+////////////////////
+////////////////////
+
+/** Returns an [Int] representing this [Byte] as if it was unsigned 0x00..0xFF */
+inline val Byte.unsigned: Int get() = this.toInt() and 0xFF
+
+/** Returns an [Int] representing this [Short] as if it was unsigned 0x0000..0xFFFF */
+inline val Short.unsigned: Int get() = this.toInt() and 0xFFFF
+
+/** Returns a [Long] representing this [Int] as if it was unsigned 0x00000000L..0xFFFFFFFFL */
+inline val Int.unsigned: Long get() = this.toLong() and 0xFFFFFFFFL

math/src/main/java/com/icegps/math/annotations/_Math_annotations.kt

@@ -0,0 +1,41 @@
+@file:Suppress("PackageDirectoryMismatch")
+
+package com.icegps.math.annotations
+
+@DslMarker
+@Target(AnnotationTarget.TYPE, AnnotationTarget.CLASS)
+annotation class KorDslMarker
+
+@Target(AnnotationTarget.TYPE, AnnotationTarget.CLASS)
+@DslMarker
+annotation class ViewDslMarker
+
+@Target(AnnotationTarget.TYPE, AnnotationTarget.CLASS)
+@DslMarker
+annotation class RootViewDslMarker
+
+@Target(AnnotationTarget.TYPE, AnnotationTarget.CLASS)
+@DslMarker
+annotation class VectorDslMarker
+
+@RequiresOptIn(level = RequiresOptIn.Level.WARNING)
+annotation class KormaExperimental(val reason: String = "")
+
+//@RequiresOptIn(level = RequiresOptIn.Level.WARNING)
+/**
+ * Mutable APIs follow the following convention:
+ *
+ * ```kotlin
+ * interface IType { val ... }
+ * class MType : IType(override var ...) : IType
+ * ```
+ *
+ * Then in usage places:
+ *
+ * ```kotlin
+ * fun doSomethingWith(a: IType, out: MType = MType()): MType
+ * ```
+ *
+ * This convention supports allocation-free APIs by being able to preallocate instances and passing them as the output.
+ */
+annotation class KormaMutableApi

math/src/main/java/com/icegps/math/geometry/AABB3D.kt

@@ -0,0 +1,55 @@
+package com.icegps.math.geometry
+
+import com.icegps.math.geometry.shape.*
+import kotlin.math.*
+
+data class AABB3D(val min: Vector3F = Vector3F(), val max: Vector3F = Vector3F()) : SimpleShape3D {
+    val minX: Float get() = min.x
+    val minY: Float get() = min.y
+    val minZ: Float get() = min.z
+
+    val maxX: Float get() = max.x
+    val maxY: Float get() = max.y
+    val maxZ: Float get() = max.z
+
+    val sizeX: Float get() = maxX - minX
+    val sizeY: Float get() = maxY - minY
+    val sizeZ: Float get() = maxZ - minZ
+
+    companion object {
+        operator fun invoke(min: Float = Float.POSITIVE_INFINITY, max: Float = Float.NEGATIVE_INFINITY): AABB3D =
+            AABB3D(Vector3F(min, min, min), Vector3F(max, max, max))
+
+        fun fromSphere(pos: Vector3F, radius: Float): AABB3D = AABB3D(
+            Vector3F(pos.x - radius, pos.y - radius, pos.z - radius),
+            Vector3F(pos.x + radius, pos.y + radius, pos.z + radius)
+        )
+    }
+
+    fun expandedToFit(that: AABB3D): AABB3D {
+        val a = this
+        val b = that
+        return AABB3D(
+            min = Vector3F(min(a.minX, b.minX), min(a.minY, b.minY), min(a.minZ, b.minZ)),
+            max = Vector3F(max(a.maxX, b.maxX), max(a.maxY, b.maxY), max(a.maxZ, b.maxZ)),
+        )
+    }
+
+    fun intersectsSphere(sphere: Sphere3D): Boolean = intersectsSphere(sphere.center, sphere.radius)
+    fun intersectsSphere(origin: Vector3F, radius: Float): Boolean = !(origin.x + radius < minX ||
+        origin.y + radius < minY ||
+        origin.z + radius < minZ ||
+        origin.x - radius > maxX ||
+        origin.y - radius > maxY ||
+        origin.z - radius > maxZ)
+
+    fun intersectsAABB(box: AABB3D): Boolean = max.x > box.min.x && min.x < box.max.x &&
+        max.y > box.min.y && min.y < box.max.y &&
+        max.z > box.min.z && min.z < box.max.z
+
+    override val center: Vector3F get() = (min + max) * 0.5f
+    override val volume: Float get() {
+        val v = (max - min)
+        return v.x * v.y * v.z
+    }
+}

math/src/main/java/com/icegps/math/geometry/Anchor.kt

@@ -0,0 +1,103 @@
+package com.icegps.math.geometry
+
+import com.icegps.math.interpolation.*
+
+typealias Anchor = Anchor2D
+typealias Anchor3 = Anchor3F
+
+data class Anchor2D(val sx: Double, val sy: Double) : Interpolable<Anchor> {
+    fun toVector(): Vector2D = Vector2D(sx, sy)
+
+    val ratioX: Ratio get() = sx.toRatio()
+    val ratioY: Ratio get() = sy.toRatio()
+
+    constructor(sx: Float, sy: Float) : this(sx.toDouble(), sy.toDouble())
+    constructor(sx: Int, sy: Int) : this(sx.toDouble(), sy.toDouble())
+
+    inline fun withX(sx: Number): Anchor = Anchor(sx.toDouble(), sy)
+    inline fun withY(sy: Number): Anchor = Anchor(sx, sy.toDouble())
+
+    inline fun withX(ratioX: Ratio): Anchor = Anchor(ratioX.toDouble(), sy)
+    inline fun withY(ratioY: Ratio): Anchor = Anchor(sx, ratioY.toDouble())
+
+    companion object {
+        inline operator fun invoke(sx: Ratio, sy: Ratio): Anchor2D = Anchor2D(sx.toDouble(), sy.toDouble())
+        inline operator fun invoke(sx: Number, sy: Number): Anchor2D = Anchor2D(sx.toDouble(), sy.toDouble())
+
+        val TOP_LEFT: Anchor = Anchor(0f, 0f)
+        val TOP_CENTER: Anchor = Anchor(.5f, 0f)
+        val TOP_RIGHT: Anchor = Anchor(1f, 0f)
+
+        val MIDDLE_LEFT: Anchor = Anchor(0f, .5f)
+        val MIDDLE_CENTER: Anchor = Anchor(.5f, .5f)
+        val MIDDLE_RIGHT: Anchor = Anchor(1f, .5f)
+
+        val BOTTOM_LEFT: Anchor = Anchor(0f, 1f)
+        val BOTTOM_CENTER: Anchor = Anchor(.5f, 1f)
+        val BOTTOM_RIGHT: Anchor = Anchor(1f, 1f)
+
+        val TOP: Anchor get() = TOP_CENTER
+        val LEFT: Anchor get() = MIDDLE_LEFT
+        val RIGHT: Anchor get() = MIDDLE_RIGHT
+        val BOTTOM: Anchor get() = BOTTOM_CENTER
+        val CENTER: Anchor get() = MIDDLE_CENTER
+    }
+
+    override fun interpolateWith(ratio: Ratio, other: Anchor): Anchor = Anchor(
+        ratio.interpolate(this.sx, other.sx),
+        ratio.interpolate(this.sy, other.sy)
+    )
+
+    fun toNamedString(): String = when (this) {
+        TOP_LEFT -> "Anchor.TOP_LEFT"
+        TOP -> "Anchor.TOP"
+        TOP_RIGHT -> "Anchor.TOP_RIGHT"
+        LEFT -> "Anchor.LEFT"
+        CENTER -> "Anchor.MIDDLE_CENTER"
+        RIGHT -> "Anchor.RIGHT"
+        BOTTOM_LEFT -> "Anchor.BOTTOM_LEFT"
+        BOTTOM_CENTER -> "Anchor.BOTTOM_CENTER"
+        BOTTOM_RIGHT -> "Anchor.BOTTOM_RIGHT"
+        else -> toString()
+    }
+}
+
+operator fun Size.times(anchor: Anchor): Point = this.toVector() * anchor.toVector()
+//operator fun SizeInt.times(anchor: Anchor): PointInt = (this.toVector().toFloat() * anchor.toVector()).toInt()
+
+data class Anchor3F(val sx: Float, val sy: Float, val sz: Float) : Interpolable<Anchor3F> {
+    fun toVector(): Vector3F = Vector3F(sx, sy, sz)
+
+    val floatX: Float get() = sx
+    val floatY: Float get() = sy
+    val floatZ: Float get() = sz
+
+    val doubleX: Double get() = sx.toDouble()
+    val doubleY: Double get() = sy.toDouble()
+    val doubleZ: Double get() = sz.toDouble()
+
+    val ratioX: Ratio get() = sx.toRatio()
+    val ratioY: Ratio get() = sy.toRatio()
+    val ratioZ: Ratio get() = sz.toRatio()
+
+    constructor(sx: Double, sy: Double, sz: Double) : this(sx.toFloat(), sy.toFloat(), sz.toFloat())
+    constructor(sx: Int, sy: Int, sz: Int) : this(sx.toFloat(), sy.toFloat(), sz.toFloat())
+
+    fun withX(sx: Float): Anchor3F = Anchor3F(sx, sy, sz)
+    fun withX(sx: Int): Anchor3F = Anchor3F(sx.toFloat(), sy, sz)
+    fun withX(sx: Double): Anchor3F = Anchor3F(sx.toFloat(), sy, sz)
+
+    fun withY(sy: Float): Anchor3F = Anchor3F(sx, sy, sz)
+    fun withY(sy: Int): Anchor3F = Anchor3F(sx, sy.toFloat(), sz)
+    fun withY(sy: Double): Anchor3F = Anchor3F(sx, sy.toFloat(), sz)
+
+    fun withZ(sz: Float): Anchor3F = Anchor3F(sx, sy, sz)
+    fun withZ(sz: Int): Anchor3F = Anchor3F(sx, sy, sz.toFloat())
+    fun withZ(sz: Double): Anchor3F = Anchor3F(sx, sy, sz.toFloat())
+
+    override fun interpolateWith(ratio: Ratio, other: Anchor3F): Anchor3F = Anchor3F(
+        ratio.interpolate(this.sx, other.sx),
+        ratio.interpolate(this.sy, other.sy),
+        ratio.interpolate(this.sz, other.sz),
+    )
+}

math/src/main/java/com/icegps/math/geometry/Angle.kt

@@ -0,0 +1,250 @@
+package com.icegps.math.geometry
+
+import com.icegps.math.*
+import com.icegps.math.interpolation.*
+import com.icegps.math.range.*
+import com.icegps.number.*
+import kotlin.math.*
+
+@PublishedApi internal const val PI2 = PI * 2.0
+@PublishedApi internal const val DEG2RAD = PI / 180.0
+@PublishedApi internal const val RAD2DEG = 180.0 / PI
+
+@PublishedApi internal fun Angle_shortDistanceTo(from: Angle, to: Angle): Angle {
+    val r0 = from.ratio.toDouble() umod 1.0
+    val r1 = to.ratio.toDouble() umod 1.0
+    val diff = (r1 - r0 + 0.5) % 1.0 - 0.5
+    return if (diff < -0.5) Angle.fromRatio(diff + 1.0) else Angle.fromRatio(diff)
+}
+
+@PublishedApi internal fun Angle_longDistanceTo(from: Angle, to: Angle): Angle {
+    val short = Angle_shortDistanceTo(from, to)
+    return when {
+        short == Angle.ZERO -> Angle.ZERO
+        short < Angle.ZERO -> Angle.FULL + short
+        else -> -Angle.FULL + short
+    }
+}
+
+@PublishedApi internal fun Angle_between(x0: Double, y0: Double, x1: Double, y1: Double, up: Vector2D = Vector2D.UP): Angle {
+    val angle = Angle.atan2(y1 - y0, x1 - x0)
+    return (if (angle < Angle.ZERO) angle + Angle.FULL else angle).adjustFromUp(up)
+}
+
+@PublishedApi internal fun Angle.adjustFromUp(up: Vector2D): Angle {
+    Orientation.checkValidUpVector(up)
+    return if (up.y > 0) this else -this
+}
+
+/**
+ * Represents an [Angle], [ratio] is in [0, 1] range, [radians] is in [0, 2PI] range, and [degrees] in [0, 360] range
+ * The internal representation is in [0, 1] range to reduce rounding errors, since floating points can represent
+ * a lot of values in that range.
+ *
+ * The equivalent old [Angle] constructor is now [Angle.fromRadians]
+ *
+ * Angles advance counter-clock-wise, starting with 0.degrees representing the right vector:
+ *
+ * Depending on what the up vector means, then numeric values of sin might be negated.
+ *
+ *   0.degrees represent right:    up=Vector2.UP: cos =+1, sin= 0 || up=Vector2.UP_SCREEN: cos =+1, sin= 0
+ *  90.degrees represents up:      up=Vector2.UP: cos = 0, sin=+1 || up=Vector2.UP_SCREEN: cos = 0, sin=-1
+ * 180.degrees represents left:    up=Vector2.UP: cos =-1, sin= 0 || up=Vector2.UP_SCREEN: cos =-1, sin= 0
+ * 270.degrees represents down:    up=Vector2.UP: cos = 0, sin=-1 || up=Vector2.UP_SCREEN: cos = 0, sin=+1
+ */
+//@KormaValueApi
+inline class Angle @PublishedApi internal constructor(
+    /** [0..1] ratio -> [0..360] degrees */
+    val radians: Double
+) : Comparable<Angle>, IsAlmostEquals<Angle> {
+    @PublishedApi inline internal val internal: Double get() = radians
+
+    /** [0..PI * 2] radians -> [0..360] degrees */
+    val ratio: Ratio get() = radiansToRatio(radians)
+    /** [0..360] degrees -> [0..PI * 2] radians -> [0..1] ratio */
+    val degrees: Double get() = radiansToDegrees(radians)
+
+    val cosine: Double get() = kotlin.math.cos(radians)
+    val sine: Double get() = kotlin.math.sin(radians)
+    val tangent: Double get() = kotlin.math.tan(radians)
+
+    fun cosine(up: Vector2D = Vector2D.UP): Double = adjustFromUp(up).cosine
+    fun sine(up: Vector2D = Vector2D.UP): Double = adjustFromUp(up).sine
+    fun tangent(up: Vector2D = Vector2D.UP): Double = adjustFromUp(up).tangent
+
+    val absoluteValue: Angle get() = Angle(internal.absoluteValue)
+    fun shortDistanceTo(other: Angle): Angle = Angle.shortDistanceTo(this, other)
+    fun longDistanceTo(other: Angle): Angle = Angle.longDistanceTo(this, other)
+
+    operator fun times(scale: Double): Angle = Angle(this.internal * scale)
+    operator fun div(scale: Double): Angle = Angle(this.internal / scale)
+    operator fun times(scale: Float): Angle = Angle(this.internal * scale)
+    operator fun div(scale: Float): Angle = Angle(this.internal / scale)
+    operator fun times(scale: Int): Angle = Angle(this.internal * scale)
+    operator fun div(scale: Int): Angle = Angle(this.internal / scale)
+    operator fun rem(angle: Angle): Angle = Angle(this.internal % angle.internal)
+    infix fun umod(angle: Angle): Angle = Angle(this.internal umod angle.internal)
+
+    operator fun div(other: Angle): Double = this.internal / other.internal // Ratio
+    operator fun plus(other: Angle): Angle = Angle(this.internal + other.internal)
+    operator fun minus(other: Angle): Angle = Angle(this.internal - other.internal)
+    operator fun unaryMinus(): Angle = Angle(-internal)
+    operator fun unaryPlus(): Angle = Angle(+internal)
+
+    fun inBetweenInclusive(min: Angle, max: Angle): Boolean = inBetween(min, max, inclusive = true)
+    fun inBetweenExclusive(min: Angle, max: Angle): Boolean = inBetween(min, max, inclusive = false)
+
+    infix fun inBetween(range: ClosedRange<Angle>): Boolean = inBetween(range.start, range.endInclusive, inclusive = true)
+    infix fun inBetween(range: OpenRange<Angle>): Boolean = inBetween(range.start, range.endExclusive, inclusive = false)
+
+    fun inBetween(min: Angle, max: Angle, inclusive: Boolean): Boolean {
+        val nthis = this.normalized
+        val nmin = min.normalized
+        val nmax = max.normalized
+        @Suppress("ConvertTwoComparisonsToRangeCheck")
+        return when {
+            nmin > nmax -> nthis >= nmin || (if (inclusive) nthis <= nmax else nthis < nmax)
+            else -> nthis >= nmin && (if (inclusive) nthis <= nmax else nthis < nmax)
+        }
+    }
+
+    override fun isAlmostEquals(other: Angle, epsilon: Double): Boolean = this.radians.isAlmostEquals(other.radians, epsilon)
+    fun isAlmostZero(epsilon: Double = 0.001): Boolean = isAlmostEquals(ZERO, epsilon)
+
+    /** Normalize between 0..1  ... 0..(PI*2).radians ... 0..360.degrees */
+    val normalized: Angle get() = fromRatio(ratio.toDouble() umod 1.0)
+    /** Normalize between -.5..+.5  ... -PI..+PI.radians ... -180..+180.degrees */
+    val normalizedHalf: Angle get() {
+        val res = normalized
+        return if (res > Angle.HALF) -Angle.FULL + res else res
+    }
+
+    override operator fun compareTo(other: Angle): Int = this.ratio.compareTo(other.ratio)
+
+    //override fun compareTo(other: Angle): Int {
+    //    //return this.radians.compareTo(other.radians) // @TODO: Double.compareTo calls EnterFrame/LeaveFrame! because it uses a Double companion object
+    //    val left = this.ratio
+    //    val right = other.ratio
+    //    // @TODO: Handle infinite/NaN? Though usually this won't happen
+    //    if (left < right) return -1
+    //    if (left > right) return +1
+    //    return 0
+    //}
+
+    override fun toString(): String = "${degrees.roundDecimalPlaces(2).niceStr}.degrees"
+
+    @Suppress("MemberVisibilityCanBePrivate")
+    companion object {
+        val EPSILON = Angle.fromRatio(0.00001)
+        val ZERO = Angle.fromRatio(0.0)
+        val QUARTER = Angle.fromRatio(0.25)
+        val HALF = Angle.fromRatio(0.5)
+        val THREE_QUARTERS = Angle.fromRatio(0.75)
+        val FULL = Angle.fromRatio(1.0)
+
+        inline fun fromRatio(ratio: Float): Angle = Angle(ratioToRadians(ratio.toRatio()))
+        inline fun fromRatio(ratio: Double): Angle = Angle(ratioToRadians(ratio.toRatio()))
+        inline fun fromRatio(ratio: Ratio): Angle = Angle(ratioToRadians(ratio))
+
+        inline fun fromRadians(radians: Double): Angle = Angle(radians)
+        inline fun fromRadians(radians: Float) = Angle(radians.toDouble())
+        inline fun fromRadians(radians: Int) = Angle(radians.toDouble())
+
+        inline fun fromDegrees(degrees: Double): Angle = Angle(degreesToRadians(degrees))
+        inline fun fromDegrees(degrees: Float) = Angle(degreesToRadians(degrees.toDouble()))
+        inline fun fromDegrees(degrees: Int) = Angle(degreesToRadians(degrees.toDouble()))
+
+        @Deprecated("", ReplaceWith("Angle.fromRatio(ratio).cosineD"))
+        inline fun cos01(ratio: Double): Double = Angle.fromRatio(ratio).cosine
+        @Deprecated("", ReplaceWith("Angle.fromRatio(ratio).sineD"))
+        inline fun sin01(ratio: Double): Double = Angle.fromRatio(ratio).sine
+        @Deprecated("", ReplaceWith("Angle.fromRatio(ratio).tangentD"))
+        inline fun tan01(ratio: Double): Double = Angle.fromRatio(ratio).tangent
+
+        inline fun atan2(x: Float, y: Float, up: Vector2D = Vector2D.UP): Angle = fromRadians(kotlin.math.atan2(x, y)).adjustFromUp(up)
+        inline fun atan2(x: Double, y: Double, up: Vector2D = Vector2D.UP): Angle = fromRadians(kotlin.math.atan2(x, y)).adjustFromUp(up)
+        inline fun atan2(p: Point, up: Vector2D = Vector2D.UP): Angle = atan2(p.x, p.y, up)
+
+        inline fun asin(v: Double): Angle = kotlin.math.asin(v).radians
+        inline fun asin(v: Float): Angle = kotlin.math.asin(v).radians
+
+        inline fun acos(v: Double): Angle = kotlin.math.acos(v).radians
+        inline fun acos(v: Float): Angle = kotlin.math.acos(v).radians
+
+        fun arcCosine(v: Double): Angle = kotlin.math.acos(v).radians
+        fun arcCosine(v: Float): Angle = kotlin.math.acos(v).radians
+
+        fun arcSine(v: Double): Angle = kotlin.math.asin(v).radians
+        fun arcSine(v: Float): Angle = kotlin.math.asin(v).radians
+
+        fun arcTangent(x: Double, y: Double): Angle = kotlin.math.atan2(x, y).radians
+        fun arcTangent(x: Float, y: Float): Angle = kotlin.math.atan2(x, y).radians
+        fun arcTangent(v: Vector2F): Angle = kotlin.math.atan2(v.x, v.y).radians
+
+        inline fun ratioToDegrees(ratio: Ratio): Double = ratio * 360.0
+        inline fun ratioToRadians(ratio: Ratio): Double = ratio * PI2
+
+        inline fun degreesToRatio(degrees: Double): Ratio = Ratio(degrees / 360.0)
+        inline fun degreesToRadians(degrees: Double): Double = degrees * DEG2RAD
+
+        inline fun radiansToRatio(radians: Double): Ratio = Ratio(radians / PI2)
+        inline fun radiansToDegrees(radians: Double): Double = radians * RAD2DEG
+
+        inline fun shortDistanceTo(from: Angle, to: Angle): Angle = Angle_shortDistanceTo(from, to)
+        inline fun longDistanceTo(from: Angle, to: Angle): Angle = Angle_longDistanceTo(from, to)
+        inline fun between(x0: Double, y0: Double, x1: Double, y1: Double, up: Vector2D = Vector2D.UP): Angle = Angle_between(x0, y0, x1, y1, up)
+
+        inline fun between(x0: Int, y0: Int, x1: Int, y1: Int, up: Vector2D = Vector2D.UP): Angle = between(x0.toDouble(), y0.toDouble(), x1.toDouble(), y1.toDouble(), up)
+        inline fun between(x0: Float, y0: Float, x1: Float, y1: Float, up: Vector2D = Vector2D.UP): Angle = between(x0.toDouble(), y0.toDouble(), x1.toDouble(), y1.toDouble(), up)
+
+        inline fun between(p0: Point, p1: Point, up: Vector2D = Vector2D.UP): Angle = between(p0.x, p0.y, p1.x, p1.y, up)
+        inline fun between(p0: Vector2F, p1: Vector2F, up: Vector2D = Vector2D.UP): Angle = between(p0.x, p0.y, p1.x, p1.y, up)
+
+        inline fun between(ox: Double, oy: Double, x1: Double, y1: Double, x2: Double, y2: Double, up: Vector2D = Vector2D.UP): Angle = between(x1 - ox, y1 - oy, x2 - ox, y2 - oy, up)
+        inline fun between(ox: Float, oy: Float, x1: Float, y1: Float, x2: Float, y2: Float, up: Vector2D = Vector2D.UP): Angle = between(x1 - ox, y1 - oy, x2 - ox, y2 - oy, up)
+
+        inline fun between(o: Point, v1: Point, v2: Point, up: Vector2D = Vector2D.UP): Angle = between(o.x, o.y, v1.x, v1.y, v2.x, v2.y, up)
+        inline fun between(o: Vector2F, v1: Vector2F, v2: Vector2F, up: Vector2D = Vector2D.UP): Angle = between(o.x, o.y, v1.x, v1.y, v2.x, v2.y, up)
+    }
+}
+
+inline fun cos(angle: Angle, up: Vector2D = Vector2D.UP): Double = angle.cosine(up)
+inline fun sin(angle: Angle, up: Vector2D = Vector2D.UP): Double = angle.sine(up)
+inline fun tan(angle: Angle, up: Vector2D = Vector2D.UP): Double = angle.tangent(up)
+
+inline fun cosf(angle: Angle, up: Vector2D = Vector2D.UP): Float = angle.cosine(up).toFloat()
+inline fun sinf(angle: Angle, up: Vector2D = Vector2D.UP): Float = angle.sine(up).toFloat()
+inline fun tanf(angle: Angle, up: Vector2D = Vector2D.UP): Float = angle.tangent(up).toFloat()
+
+inline fun abs(angle: Angle): Angle = angle.absoluteValue
+inline fun min(a: Angle, b: Angle): Angle = Angle(min(a.internal, b.internal))
+inline fun max(a: Angle, b: Angle): Angle = Angle(max(a.internal, b.internal))
+
+fun Angle.clamp(min: Angle, max: Angle): Angle = min(max(this, min), max)
+
+operator fun ClosedRange<Angle>.contains(angle: Angle): Boolean = angle.inBetween(this.start, this.endInclusive, inclusive = true)
+operator fun OpenRange<Angle>.contains(angle: Angle): Boolean = angle.inBetween(this.start, this.endExclusive, inclusive = false)
+infix fun Angle.until(other: Angle): OpenRange<Angle> = OpenRange(this, other)
+
+val Double.degrees: Angle get() = Angle.fromDegrees(this)
+val Double.radians: Angle get() = Angle.fromRadians(this)
+val Int.degrees: Angle get() = Angle.fromDegrees(this)
+val Int.radians: Angle get() = Angle.fromRadians(this)
+val Float.degrees: Angle get() = Angle.fromDegrees(this)
+val Float.radians: Angle get() = Angle.fromRadians(this)
+
+fun Ratio.interpolateAngle(l: Angle, r: Angle, minimizeAngle: Boolean): Angle = _interpolateAngleAny(this, l, r, minimizeAngle)
+fun Ratio.interpolateAngle(l: Angle, r: Angle): Angle = interpolateAngle(l, r, minimizeAngle = true)
+fun Ratio.interpolateAngleNormalized(l: Angle, r: Angle): Angle = interpolateAngle(l, r, minimizeAngle = true)
+fun Ratio.interpolateAngleDenormalized(l: Angle, r: Angle): Angle = interpolateAngle(l, r, minimizeAngle = false)
+
+private fun _interpolateAngleAny(ratio: Ratio, l: Angle, r: Angle, minimizeAngle: Boolean = true): Angle {
+    if (!minimizeAngle) return Angle.fromRatio(ratio.interpolate(l.ratio, r.ratio))
+    val ln = l.normalized
+    val rn = r.normalized
+    return when {
+        (rn - ln).absoluteValue <= Angle.HALF -> Angle.fromRadians(ratio.interpolate(ln.radians, rn.radians))
+        ln < rn -> Angle.fromRadians(ratio.interpolate((ln + Angle.FULL).radians, rn.radians)).normalized
+        else -> Angle.fromRadians(ratio.interpolate(ln.radians, (rn + Angle.FULL).radians)).normalized
+    }
+}

math/src/main/java/com/icegps/math/geometry/BoundsBuilder.kt

@@ -0,0 +1,60 @@
+package com.icegps.math.geometry
+
+inline class BoundsBuilder(val bounds: Rectangle) {
+    val isEmpty: Boolean get() = bounds.isNIL
+    val isNotEmpty: Boolean get() = bounds.isNotNIL
+
+    val xmin: Double get() = kotlin.math.min(bounds.left, bounds.right)
+    val xmax: Double get() = kotlin.math.max(bounds.left, bounds.right)
+    val ymin: Double get() = kotlin.math.min(bounds.top, bounds.bottom)
+    val ymax: Double get() = kotlin.math.max(bounds.top, bounds.bottom)
+
+    /** Minimum value found for X. [default] if ![hasPoints] */
+    fun xminOr(default: Double = 0.0): Double = if (hasPoints) xmin else default
+    /** Maximum value found for X. [default] if ![hasPoints] */
+    fun xmaxOr(default: Double = 0.0): Double = if (hasPoints) xmax else default
+    /** Minimum value found for Y. [default] if ![hasPoints] */
+    fun yminOr(default: Double = 0.0): Double = if (hasPoints) ymin else default
+    /** Maximum value found for Y. [default] if ![hasPoints] */
+    fun ymaxOr(default: Double = 0.0): Double = if (hasPoints) ymax else default
+
+    val hasPoints: Boolean get() = isNotEmpty
+
+    companion object {
+        val EMPTY = BoundsBuilder(Rectangle.NIL)
+
+        operator fun invoke(): BoundsBuilder = EMPTY
+        operator fun invoke(p1: Point): BoundsBuilder = BoundsBuilder(Rectangle(p1, Size(0, 0)))
+        operator fun invoke(p1: Point, p2: Point): BoundsBuilder = BoundsBuilder(Rectangle.fromBounds(Point.minComponents(p1, p2), Point.maxComponents(p1, p2)))
+        operator fun invoke(p1: Point, p2: Point, p3: Point): BoundsBuilder = BoundsBuilder(Rectangle.fromBounds(Point.minComponents(p1, p2, p3), Point.maxComponents(p1, p2, p3)))
+        operator fun invoke(p1: Point, p2: Point, p3: Point, p4: Point): BoundsBuilder = BoundsBuilder(Rectangle.fromBounds(Point.minComponents(p1, p2, p3, p4), Point.maxComponents(p1, p2, p3, p4)))
+        operator fun invoke(size: Int, func: BoundsBuilder.(Int) -> BoundsBuilder): BoundsBuilder {
+            var bb = BoundsBuilder()
+            for (n in 0 until size) bb = func(bb, n)
+            return bb
+        }
+    }
+    fun plus(x: Double, y: Double): BoundsBuilder = this.plus(Point(x, y))
+    operator fun plus(p: Point): BoundsBuilder {
+        if (bounds.isNIL) return BoundsBuilder(Rectangle(p, Size(0, 0)))
+        return BoundsBuilder(Rectangle.fromBounds(Point.minComponents(bounds.topLeft, p), Point.maxComponents(bounds.bottomRight, p)))
+    }
+    operator fun plus(bb: BoundsBuilder): BoundsBuilder = this + bb.bounds
+    operator fun plus(rect: Rectangle?): BoundsBuilder {
+        if (rect == null) return this
+        if (rect.isNIL) return this
+        return this + rect.topLeft + rect.bottomRight
+    }
+    operator fun plus(p: IPointList): BoundsBuilder {
+        var bb = this
+        for (n in 0 until p.size) bb = bb.plus(p[n])
+        return bb
+    }
+    //operator fun plus(rect: Rectangle): BoundsBuilder = TODO()
+    operator fun plus(rects: List<Rectangle>): BoundsBuilder {
+        var bb = this
+        for (it in rects) bb += it
+        return bb
+    }
+    fun boundsOrNull(): Rectangle? = if (isEmpty) null else bounds
+}

math/src/main/java/com/icegps/math/geometry/Circle.kt

@@ -0,0 +1,29 @@
+package com.icegps.math.geometry
+
+import com.icegps.math.geometry.shape.*
+import kotlin.math.*
+
+data class Circle(override val center: Point, val radius: Double) : SimpleShape2D {
+    companion object {
+        inline operator fun invoke(center: Point, radius: Number) = Circle(center, radius.toDouble())
+        inline operator fun invoke(x: Number, y: Number, radius: Number) = Circle(Point(x.toDouble(), y.toDouble()), radius.toDouble())
+    }
+
+    override val closed: Boolean get() = true
+
+    override val area: Double get() = (PI * radius * radius)
+    override val perimeter: Double get() = (PI * 2.0 * radius)
+    override fun distance(p: Point): Double = (p - center).length - radius
+    override fun normalVectorAt(p: Point): Vector2D = (p - center).normalized
+
+    val radiusSquared: Double get() = radius * radius
+
+    fun distanceToCenterSquared(p: Point): Double = Point.distanceSquared(p, center)
+    // @TODO: Check if inside the circle
+    fun distanceClosestSquared(p: Point): Double = distanceToCenterSquared(p) - radiusSquared
+    // @TODO: Check if inside the circle
+    fun distanceFarthestSquared(p: Point): Double = distanceToCenterSquared(p) + radiusSquared
+    override fun projectedPoint(p: Point): Point = Point.polar(center, Angle.between(center, p), radius)
+    override fun containsPoint(p: Point): Boolean = (p - center).length <= radius
+    override fun getBounds(): Rectangle = Rectangle.fromBounds(center.x - radius, center.y - radius, center.x + radius, center.y + radius,)
+}

math/src/main/java/com/icegps/math/geometry/Ellipse.kt

@@ -0,0 +1,83 @@
+package com.icegps.math.geometry
+
+import com.icegps.math.geometry.shape.*
+import kotlin.math.*
+
+data class Ellipse(override val center: Point, val radius: Size) : SimpleShape2D {
+    override val area: Double get() = (PI * radius.width * radius.height)
+    override val perimeter: Double get() {
+        if (radius.width == radius.height) return (PI * 2.0 * radius.width) // Circle formula
+        val (a, b) = radius
+        val h = ((a - b) * (a - b)) / ((a + b) * (a + b))
+        return (PI * (a + b) * (1 + ((3 * h) / (10 + sqrt(4 - (3 * h))))))
+    }
+
+    override fun distance(p: Point): Double {
+        val p = p - center
+        val scaledPoint = Vector2D(p.x / radius.width, p.y / radius.height)
+        val length = scaledPoint.length
+        return (length - 1) * min(radius.width, radius.height)
+    }
+
+    override fun normalVectorAt(p: Point): Vector2D {
+        val pointOnEllipse = p - center
+        val (a, b) = radius
+        val normal = Vector2D(pointOnEllipse.x / (a * a), pointOnEllipse.y / (b * b))
+        return normal.normalized
+        //val d = p - center
+        //val r2 = radius.toVector() * radius.toVector()
+        //return (d / r2).normalized
+    }
+
+    override fun projectedPoint(p: Point): Point {
+        val angle = Angle.between(center, p)
+        return center + Point(radius.width * angle.cosine, radius.height * angle.sine)
+
+        //val k = (radius.width * radius.height) / sqrt()
+        //return projectPointOntoEllipse(p, center, radius.toVector())
+    }
+
+    override fun containsPoint(p: Point): Boolean {
+        if (radius.isEmpty()) return false
+        // Check if the point is inside the ellipse using the ellipse equation:
+        // (x - centerX)^2 / radiusX^2 + (y - centerY)^2 / radiusY^2 <= 1
+        return ((p.x - center.x).pow(2) / radius.width.pow(2)) + ((p.y - center.y).pow(2) / radius.height.pow(2)) <= 1
+    }
+
+    override val closed: Boolean get() = true
+    override fun getBounds(): Rectangle = Rectangle.fromBounds(center.x - radius.width, center.y - radius.height, center.x + radius.width, center.y + radius.height)
+
+    companion object {
+        private fun projectPointOntoEllipse(point: Vector2F, center: Vector2F, radius: Vector2F, tolerance: Double = 1e-6, maxIterations: Int = 100): Vector2F {
+            var currentPoint = point
+            var i = 0
+
+            while (i < maxIterations) {
+                val dx = currentPoint.x - center.x
+                val dy = currentPoint.y - center.y
+                val rx2 = radius.x * radius.x
+                val ry2 = radius.y * radius.y
+
+                val f = Vector2F(
+                    (dx * rx2 - dy * dx * dy) / (rx2 * ry2),
+                    (dy * ry2 - dx * dy * dx) / (rx2 * ry2)
+                )
+
+                val df = Vector2F(
+                    (ry2 - 2.0 * dy * dy) / (rx2 * ry2),
+                    (rx2 - 2.0 * dx * dx) / (rx2 * ry2)
+                )
+
+                val nextPoint = currentPoint - f / df
+                val dist = (nextPoint - currentPoint).length
+
+                if (dist < tolerance) return nextPoint
+
+                currentPoint = nextPoint
+                i++
+            }
+
+            return currentPoint
+        }
+    }
+}

math/src/main/java/com/icegps/math/geometry/EulerRotation.kt

@@ -0,0 +1,330 @@
+package com.icegps.math.geometry
+
+import com.icegps.math.*
+import kotlin.math.*
+
+/**
+ * Rotations around Z axis, then X axis, then Y axis in that order.
+ */
+inline class EulerRotation private constructor(val data: Vector4F) : IsAlmostEqualsF<EulerRotation> {
+    val config: Config get() = Config(data.w.toInt())
+    val order: Order get() = config.order
+    val coordinateSystem: CoordinateSystem get() = config.coordinateSystem
+
+    enum class Order(
+        val x: Int, val y: Int, val z: Int, val w: Int, val str: String,
+    ) {
+        INVALID(0, 0, 0, 0, "XXX"),
+        XYZ(+1, -1, +1, -1, "XYZ"),
+        XZY(-1, -1, +1, +1, "XZY"),
+        YXZ(+1, -1, -1, +1, "YXZ"),
+        YZX(+1, +1, -1, -1, "YZX"),
+        ZXY(-1, +1, +1, -1, "ZXY"),
+        ZYX(-1, +1, -1, +1, "ZYX"),
+        ;
+
+        fun withCoordinateSystem(coordinateSystem: CoordinateSystem) = if (coordinateSystem.sign < 0) reversed() else this
+
+        fun reversed(): Order = when (this) {
+            INVALID -> INVALID
+            XYZ -> ZYX
+            XZY -> YZX
+            YXZ -> ZXY
+            YZX -> XZY
+            ZXY -> YXZ
+            ZYX -> XYZ
+        }
+
+        fun indexAt(pos: Int, reversed: Boolean = false): Int = str[(if (reversed) 2 - pos else pos) umod 3] - 'X'
+
+        override fun toString(): String = "$name [$x, $y, $z, $w]"
+
+        companion object {
+            val VALUES = values()
+            val DEFAULT = XYZ
+        }
+    }
+    //enum class Normalized { NO, FULL_ANGLE, HALF_ANGLE }
+    inline class Config(val id: Int) {
+        //constructor(order: Order, coordinateSystem: CoordinateSystem) : this(order.ordinal * coordinateSystem.sign)
+        constructor(order: Order, coordinateSystem: CoordinateSystem) : this(order.withCoordinateSystem(coordinateSystem).ordinal)
+
+        val order: Order get() = Order.VALUES[id.absoluteValue]
+        val coordinateSystem: CoordinateSystem get() = if (id < 0) CoordinateSystem.LEFT_HANDED else CoordinateSystem.RIGHT_HANDED
+
+        override fun toString(): String = "EulerRotation.Config(order=$order, coordinateSystem=$coordinateSystem)"
+
+        companion object {
+            val UNITY get() = Config(Order.ZXY, CoordinateSystem.LEFT_HANDED)
+            //val UNITY get() = LIBGDX
+            val UNREAL get() = Config(Order.ZYX, CoordinateSystem.LEFT_HANDED)
+            //val UNREAL get() = THREEJS
+            val GODOT get() = Config(Order.YXZ, CoordinateSystem.RIGHT_HANDED)
+            val LIBGDX get() = Config(Order.YXZ, CoordinateSystem.RIGHT_HANDED)
+            val THREEJS get() = Config(Order.XYZ, CoordinateSystem.RIGHT_HANDED)
+
+            // Same as Three.JS
+            val DEFAULT get() = Config(Order.XYZ, CoordinateSystem.RIGHT_HANDED)
+        }
+    }
+    enum class CoordinateSystem(val sign: Int) {
+        LEFT_HANDED(-1), RIGHT_HANDED(+1);
+        val rsign = -sign
+    }
+
+    val roll: Angle get() = Angle.fromRatio(data.x)
+    val pitch: Angle get() = Angle.fromRatio(data.y)
+    val yaw: Angle get() = Angle.fromRatio(data.z)
+
+    @Deprecated("", ReplaceWith("roll")) val x: Angle get() = roll
+    @Deprecated("", ReplaceWith("pitch")) val y: Angle get() = pitch
+    @Deprecated("", ReplaceWith("yaw")) val z: Angle get() = yaw
+
+    override fun toString(): String = "EulerRotation(roll=$roll, pitch=$pitch, yaw=$yaw)"
+
+    fun copy(roll: Angle = this.roll, pitch: Angle = this.pitch, yaw: Angle = this.yaw): EulerRotation = EulerRotation(roll, pitch, yaw)
+    constructor() : this(Angle.ZERO, Angle.ZERO, Angle.ZERO)
+    constructor(roll: Angle, pitch: Angle, yaw: Angle, config: Config = Config.DEFAULT)
+        : this(Vector4F(roll.ratio.toFloat(), pitch.ratio.toFloat(), yaw.ratio.toFloat(), config.id.toFloat()))
+
+    fun normalized(): EulerRotation = EulerRotation(roll.normalized, pitch.normalized, yaw.normalized)
+    fun normalizedHalf(): EulerRotation = EulerRotation(roll.normalizedHalf, pitch.normalizedHalf, yaw.normalizedHalf)
+
+    fun toMatrix(): Matrix4 = toQuaternion().toMatrix()
+    fun toQuaternion(): Quaternion = _toQuaternion(x, y, z, config)
+    override fun isAlmostEquals(other: EulerRotation, epsilon: Float): Boolean =
+        this.data.isAlmostEquals(other.data, epsilon)
+
+    companion object {
+        fun toQuaternion(roll: Angle, pitch: Angle, yaw: Angle, config: Config = Config.DEFAULT): Quaternion {
+            return _toQuaternion(roll, pitch, yaw, config)
+        }
+        // http://www.mathworks.com/matlabcentral/fileexchange/20696-function-to-convert-between-dcm-euler-angles-quaternions-and-euler-vectors/content/SpinCalc.m
+        private fun _toQuaternion(x: Angle, y: Angle, z: Angle, config: Config = Config.DEFAULT): Quaternion {
+            val order = config.order
+            val coordinateSystem = config.coordinateSystem
+            val sign = coordinateSystem.sign
+            //println("ORDER=$order, coordinateSystem=$coordinateSystem, sign=$sign")
+
+            val c1 = cos(x / 2)
+            val c2 = cos(y / 2)
+            val c3 = cos(z / 2)
+            val s1 = sin(x / 2)
+            val s2 = sin(y / 2)
+            val s3 = sin(z / 2)
+
+            return Quaternion(
+                ((s1 * c2 * c3) + ((c1 * s2 * s3) * order.x * sign)),
+                ((c1 * s2 * c3) + ((s1 * c2 * s3) * order.y * sign)),
+                ((c1 * c2 * s3) + ((s1 * s2 * c3) * order.z * sign)),
+                ((c1 * c2 * c3) + ((s1 * s2 * s3) * order.w * sign)),
+            )
+        }
+
+        fun fromRotationMatrix(m: Matrix3, config: Config = Config.DEFAULT): EulerRotation {
+            //val config = if (config == Config.UNITY) Config.LIBGDX else config
+            val order = config.order
+            val coordinateSystem = config.coordinateSystem
+
+            val sign = coordinateSystem.sign
+
+            //val m = if (sign < 0) m.transposed() else m
+            //val m = m
+
+            val m11 = m.v00
+            val m12 = m.v01
+            val m13 = m.v02
+
+            val m21 = m.v10
+            val m22 = m.v11
+            val m23 = m.v12
+
+            val m31 = m.v20
+            val m32 = m.v21
+            val m33 = m.v22
+
+            val x: Angle
+            val y: Angle
+            val z: Angle
+
+            when (order) {
+                Order.XYZ -> {
+                    x = if (m13.absoluteNotAlmostOne) Angle.atan2(-m23, m33) else Angle.atan2(m32, m22)
+                    y = Angle.asin(m13.clamp(-1f, +1f))
+                    z = if (m13.absoluteNotAlmostOne) Angle.atan2(-m12, m11) else Angle.ZERO
+                }
+                Order.YXZ -> {
+                    x = Angle.asin(-(m23.clamp(-1f, +1f)))
+                    y = if (m23.absoluteNotAlmostOne) Angle.atan2(m13, m33) else Angle.atan2(-m31, m11)
+                    z = if (m23.absoluteNotAlmostOne) Angle.atan2(m21, m22) else Angle.ZERO
+                }
+                Order.ZXY -> {
+                    y = Angle.asin(m32.clamp(-1f, +1f))
+                    x = if (m32.absoluteNotAlmostOne) Angle.atan2(-m31, m33) else Angle.ZERO
+                    z = if (m32.absoluteNotAlmostOne) Angle.atan2(-m12, m22) else Angle.atan2(m21, m11)
+                }
+                Order.ZYX -> {
+                    x = if (m31.absoluteNotAlmostOne) Angle.atan2(m32, m33) else Angle.ZERO
+                    y = Angle.asin(-(m31.clamp(-1f, +1f)))
+                    z = if (m31.absoluteNotAlmostOne) Angle.atan2(m21, m11) else Angle.atan2(-m12, m22)
+                }
+                Order.YZX -> {
+                    x = if (m21.absoluteNotAlmostOne) Angle.atan2(-m23, m22) else Angle.ZERO
+                    y = if (m21.absoluteNotAlmostOne) Angle.atan2(-m31, m11) else Angle.atan2(m13, m33)
+                    z = Angle.asin(m21.clamp(-1f, +1f))
+                }
+                Order.XZY -> {
+                    x = if (m12.absoluteNotAlmostOne) Angle.atan2(m32, m22) else Angle.atan2(-m23, m33)
+                    y = if (m12.absoluteNotAlmostOne) Angle.atan2(m13, m11) else Angle.ZERO
+                    z = Angle.asin(-(m12.clamp(-1f, +1f)))
+                }
+                Order.INVALID -> error("Invalid")
+            }
+
+            //println("order=$order, coordinateSystem=$coordinateSystem : ${coordinateSystem.sign}, x=$x, y=$y, z=$z")
+
+            //val sign = coordinateSystem.sign
+            //return EulerRotation(x * coordinateSystem.sign, y * coordinateSystem.sign, z * coordinateSystem.sign, config)
+            //return EulerRotation(x * sign, y * sign, z * sign, config)
+            return EulerRotation(x, y, z, config)
+        }
+
+        private val Float.absoluteNotAlmostOne: Boolean get() = absoluteValue < 0.9999999
+
+
+        fun fromQuaternion(q: Quaternion, config: Config = Config.DEFAULT): EulerRotation {
+            return fromRotationMatrix(q.toMatrix3(), config)
+            /*
+            //return fromQuaternion(q.x, q.y, q.z, q.w, config)
+
+            val extrinsic = false
+
+            // intrinsic/extrinsic conversion helpers
+            val angle_first: Int
+            val angle_third: Int
+            val reversed: Boolean
+            if (extrinsic) {
+                angle_first = 0
+                angle_third = 2
+                reversed = false
+            } else {
+                reversed = true
+                //reversed = false
+                //seq = seq[:: - 1]
+                angle_first = 2
+                angle_third = 0
+            }
+
+            val quat = q
+            val i = config.order.indexAt(0, reversed = reversed)
+            val j = config.order.indexAt(1, reversed = reversed)
+            val symmetric = i == j
+            var k = if (symmetric) 3 - i - j else config.order.indexAt(2, reversed = reversed)
+            val sign = (i - j) * (j - k) * (k - i) / 2
+
+            println("ORDER: $i, $j, $k")
+            val eps = 1e-7f
+
+            val _angles = FloatArray(3)
+            //_angles = angles[ind, :]
+
+            // Step 1
+            // Permutate quaternion elements
+            val a: Float
+            val b: Float
+            val c: Float
+            val d: Float
+            if (symmetric) {
+                a = quat[3]
+                b = quat[i]
+                c = quat[j]
+                d = quat[k] * sign
+            } else {
+                a = quat[3] - quat[j]
+                b = quat[i] + quat[k] * sign
+                c = quat[j] + quat[3]
+                d = quat[k] * sign - quat[i]
+            }
+
+            // Step 2
+            // Compute second angle...
+            _angles[1] = 2 * atan2(hypot(c, d), hypot(a, b))
+
+            // ... and check if equal to is 0 or pi, causing a singularity
+            val case = when {
+                abs(_angles[1]) <= eps -> 1
+                abs(_angles[1] - PIF) <= eps -> 2
+                else -> 0 // normal case
+            }
+
+            // Step 3
+            // compute first and third angles, according to case
+            val half_sum = atan2(b, a)
+            val half_diff = atan2(d, c)
+
+            if (case == 0) {  // no singularities
+                _angles[angle_first] = half_sum - half_diff
+                _angles[angle_third] = half_sum + half_diff
+            } else {  // any degenerate case
+                _angles[2] = 0f
+                if (case == 1) {
+                    _angles[0] = 2 * half_sum
+                } else {
+                    _angles[0] = 2 * half_diff * (if (extrinsic) -1 else 1)
+                }
+            }
+
+            // for Tait-Bryan angles
+            if (!symmetric) {
+                _angles[angle_third] *= sign.toFloat()
+                _angles[1] -= PIF / 2
+            }
+
+            for (idx in 0 until 3) {
+                if (_angles[idx] < -PIF) {
+                    _angles[idx] += 2 * PIF
+                } else if (_angles[idx] > PIF) {
+                    _angles[idx] -= 2 * PIF
+                }
+            }
+
+            if (case != 0) {
+                println(
+                    "Gimbal lock detected. Setting third angle to zero " +
+                    "since it is not possible to uniquely determine " +
+                    "all angles."
+                )
+            }
+
+            return EulerRotation(_angles[0].radians, _angles[2].radians, _angles[1].radians * config.coordinateSystem.sign)
+            */
+        }
+
+        fun fromQuaternion(x: Float, y: Float, z: Float, w: Float, config: Config = Config.DEFAULT): EulerRotation {
+
+            return fromQuaternion(Quaternion(x, y, z, w), config)
+            /*
+            val t = y * x + z * w
+            // Gimbal lock, if any: positive (+1) for north pole, negative (-1) for south pole, zero (0) when no gimbal lock
+            val pole = if (t > 0.499f) 1 else if (t < -0.499f) -1 else 0
+            println("pole=$pole")
+            println(Angle.atan2(2f * (y * w + x * z), 1f - 2f * (y * y + x * x)))
+            return EulerRotation(
+                roll = when (pole) {
+                    0 -> Angle.asin((2f * (w * x - z * y)).clamp(-1f, +1f))
+                    else -> (pole.toFloat() * PIF * .5f).radians
+                },
+                pitch = when (pole) {
+                    0 -> Angle.atan2(2f * (y * w + x * z), 1f - 2f * (y * y + x * x))
+                    else -> Angle.ZERO
+                },
+                yaw = when (pole) {
+                    0 -> Angle.atan2(2f * (w * z + y * x), 1f - 2f * (x * x + z * z))
+                    else -> Angle.atan2(y, w) * pole.toFloat() * 2f
+                },
+            )
+
+             */
+        }
+    }
+}

math/src/main/java/com/icegps/math/geometry/IPointList.kt

@@ -0,0 +1,124 @@
+package com.icegps.math.geometry
+
+import com.icegps.math.*
+import com.icegps.number.*
+import kotlin.math.*
+
+interface IGenericDoubleVector {
+    val dimensions: Int
+    operator fun get(dim: Int): Double
+    operator fun set(dim: Int, value: Double)
+}
+
+interface IDoubleVectorList : IsAlmostEquals<IDoubleVectorList> {
+    fun isEmpty(): Boolean = size == 0
+    fun isNotEmpty(): Boolean = size != 0
+
+    val size: Int
+    val dimensions: Int
+    operator fun get(index: Int, dim: Int): Double
+
+    override fun isAlmostEquals(other: IDoubleVectorList, epsilon: Double): Boolean {
+        if (this.size != other.size) return false
+        if (this.dimensions != other.dimensions) return false
+        for (dim in 0 until dimensions) for (n in 0 until size) {
+            if (!this[n, dim].isAlmostEquals(other[n, dim], epsilon)) return false
+        }
+        return true
+    }
+}
+
+// @TODO: Potential candidate for value class when multiple values are supported
+class GenericDoubleVector(override val dimensions: Int, val data: DoubleArray, val offset: Int = 0) : IGenericDoubleVector {
+    constructor(vararg data: Double) : this(data.size, data)
+    constructor(vararg data: Float) : this(data.size, DoubleArray(data.size) { data[it].toDouble() })
+    constructor(vararg data: Int) : this(data.size, DoubleArray(data.size) { data[it].toDouble() })
+
+    override operator fun get(dim: Int): Double = data[offset + dim]
+    override operator fun set(dim: Int, value: Double) { data[offset + dim] = value }
+
+    override fun toString(): String = buildString { toStringBuilder(this) }
+}
+
+val IGenericDoubleVector.length: Double get() {
+    var ssum = 0.0
+    for (n in 0 until dimensions) ssum += this[n]
+    return sqrt(ssum)
+}
+
+fun IGenericDoubleVector.toStringBuilder(out: StringBuilder) {
+    out.appendGenericArray(dimensions) { appendNice(this@toStringBuilder[it]) }
+}
+
+interface IPointList : IDoubleVectorList, List<Point> {
+    override val size: Int
+    override fun isEmpty(): Boolean = size == 0
+    fun getX(index: Int): Double
+    fun getY(index: Int): Double
+    override val dimensions: Int get() = 2
+    override operator fun get(index: Int): Point = Point(getX(index), getY(index))
+    override fun contains(element: Point): Boolean = indexOf(element) >= 0
+    override fun containsAll(elements: Collection<Point>): Boolean = containsAllSet(elements)
+    override fun indexOf(element: Point): Int = indexOf(this, element)
+    override fun lastIndexOf(element: Point): Int = lastIndexOf(this, element)
+    override fun iterator(): Iterator<Point> = listIterator()
+    override fun listIterator(): ListIterator<Point> = listIterator(0)
+    override fun listIterator(index: Int): ListIterator<Point> = Sublist(this, 0, size).listIterator(index)
+    override fun subList(fromIndex: Int, toIndex: Int): List<Point> = Sublist(this, fromIndex, toIndex)
+
+    class Sublist(val list: IPointList, val fromIndex: Int, val toIndex: Int) : List<Point> {
+        override val size: Int = toIndex - fromIndex
+        override fun get(index: Int): Point = list[index + fromIndex]
+        override fun isEmpty(): Boolean = size == 0
+
+        override fun iterator(): Iterator<Point> = listIterator()
+        override fun listIterator(): ListIterator<Point> = listIterator(0)
+        override fun listIterator(index: Int): ListIterator<Point> = object : ListIterator<Point> {
+            var current = index
+            override fun hasNext(): Boolean = current >= size
+            override fun hasPrevious(): Boolean = current > index
+            override fun next(): Point = this@Sublist[current++]
+            override fun nextIndex(): Int = current + 1
+            override fun previous(): Point = this@Sublist[--current]
+            override fun previousIndex(): Int = current - 1
+        }
+
+        override fun subList(fromIndex: Int, toIndex: Int): List<Point> = Sublist(list, this.fromIndex + fromIndex, this.fromIndex + toIndex)
+        override fun lastIndexOf(element: Point): Int = lastIndexOf(list, element, fromIndex, toIndex, offset = -fromIndex)
+        override fun indexOf(element: Point): Int = indexOf(list, element, fromIndex, toIndex, offset = -fromIndex)
+        override fun containsAll(elements: Collection<Point>): Boolean = containsAllSet(elements)
+        override fun contains(element: Point): Boolean = indexOf(element) >= 0
+    }
+
+    companion object {
+        fun <T> Collection<T>.containsAllSet(elements: Collection<T>): Boolean {
+            val s = elements.toSet()
+            return all { it in s }
+        }
+
+        fun indexOf(list: IPointList, element: Point, fromIndex: Int = 0, toIndex: Int = list.size, offset: Int = 0): Int {
+            for (n in fromIndex until toIndex) if (list.getX(n) == element.x && list.getY(n) == element.y) return n + offset
+            return -1
+        }
+        fun lastIndexOf(list: IPointList, element: Point, fromIndex: Int = 0, toIndex: Int = list.size, offset: Int = 0): Int {
+            for (n in toIndex - 1 downTo  fromIndex) if (list.getX(n) == element.x && list.getY(n) == element.y) return n + offset
+            return -1
+        }
+
+        inline fun getPolylineLength(size: Int, crossinline get: (n: Int) -> Point): Double {
+            var out = 0.0
+            var prev = Point.ZERO
+            for (n in 0 until size) {
+                val p = get(n)
+                if (n > 0) out += Point.distance(prev, p)
+                prev = p
+            }
+            return out
+        }
+
+    }
+}
+
+
+fun IPointList.getPolylineLength(): Double = IPointList.getPolylineLength(size) { get(it) }
+fun List<Point>.getPolylineLength(): Double = IPointList.getPolylineLength(size) { get(it) }

math/src/main/java/com/icegps/math/geometry/Line.kt

@@ -0,0 +1,175 @@
+package com.icegps.math.geometry
+
+import com.icegps.math.*
+import com.icegps.math.annotations.*
+import com.icegps.math.geometry.shape.*
+import kotlin.math.*
+
+typealias Line2 = Line
+typealias Line = Line2D
+
+//@KormaValueApi
+data class Line2D(val a: Vector2D, val b: Vector2D) : SimpleShape2D {
+    override val closed: Boolean get() = false
+
+    override val area: Double get() = 0.0
+    override val perimeter: Double get() = length
+
+    override fun normalVectorAt(p: Point): Vector2D {
+        val projected = projectedPoint(p)
+        return (b - a).toNormal().normalized * Point.crossProduct(projected, p).sign
+    }
+
+    override val center: Point get() = (a + b) * 0.5
+
+    fun toRay(): Ray = Ray(a, (b - a).normalized)
+
+    val xmin: Double get() = kotlin.math.min(x0, x1)
+    val xmax: Double get() = kotlin.math.max(x0, x1)
+    val ymin: Double get() = kotlin.math.min(y0, y1)
+    val ymax: Double get() = kotlin.math.max(y0, y1)
+
+    override fun projectedPoint(p: Point): Point {
+        return projectedPointOutsideSegment(p).clamp(Point(xmin, ymin), Point(xmax, ymax))
+    }
+
+    fun projectedPointOutsideSegment(p: Point): Point {
+        val v1x = x0
+        val v2x = x1
+        val v1y = y0
+        val v2y = y1
+        val px = p.x
+        val py = p.y
+
+        // return this.getIntersectionPoint(Line(point, Point.fromPolar(point, this.angle + 90.degrees)))!!
+        // get dot product of e1, e2
+        val e1x = v2x - v1x
+        val e1y = v2y - v1y
+        val e2x = px - v1x
+        val e2y = py - v1y
+        val valDp = Point.dot(e1x, e1y, e2x, e2y)
+        // get length of vectors
+
+        val lenLineE1 = kotlin.math.hypot(e1x, e1y)
+        val lenLineE2 = kotlin.math.hypot(e2x, e2y)
+
+        // What happens if lenLineE1 or lenLineE2 are zero?, it would be a division by zero.
+        // Does that mean that the point is on the line, and we should use it?
+        if (lenLineE1 == 0.0 || lenLineE2 == 0.0) {
+            return Point(px, py)
+        }
+
+        val cos = valDp / (lenLineE1 * lenLineE2)
+
+        // length of v1P'
+        val projLenOfLine = cos * lenLineE2
+
+        return Point((v1x + (projLenOfLine * e1x) / lenLineE1), (v1y + (projLenOfLine * e1y) / lenLineE1))
+    }
+
+    override fun containsPoint(p: Point): Boolean = false
+    override fun getBounds(): Rectangle {
+        TODO("Not yet implemented")
+    }
+
+    constructor() : this(Point(), Point())
+    constructor(x0: Double, y0: Double, x1: Double, y1: Double) : this(Point(x0, y0), Point(x1, y1))
+    constructor(x0: Float, y0: Float, x1: Float, y1: Float) : this(Point(x0, y0), Point(x1, y1))
+    constructor(x0: Int, y0: Int, x1: Int, y1: Int) : this(Point(x0, y0), Point(x1, y1))
+
+    inline fun flipped(): Line = Line(b, a)
+
+    val x0: Double get() = a.x
+    val y0: Double get() = a.y
+
+    val x1: Double get() = b.x
+    val y1: Double get() = b.y
+
+    val dx: Double get() = x1 - x0
+    val dy: Double get() = y1 - y0
+
+    val min: Point get() = Point(minX, minY)
+    val minX: Double get() = kotlin.math.min(a.x, b.x)
+    val minY: Double get() = kotlin.math.min(a.y, b.y)
+
+    val max: Point get() = Point(maxX, maxY)
+    val maxX: Double get() = kotlin.math.max(a.x, b.x)
+    val maxY: Double get() = kotlin.math.max(a.y, b.y)
+
+    fun round(): Line = Line(a.round(), b.round())
+    fun directionVector(): Point = Point(dx, dy)
+
+    fun getMinimumDistance(p: Point): Double {
+        val v = a
+        val w = b
+        val l2 = Point.distanceSquared(v, w)
+        if (l2 == 0.0) return Point.distanceSquared(p, a)
+        val t = (Point.dot(p - v, w - v) / l2).clamp(0.0, 1.0)
+        return Point.distance(p, v + (w - v) * t)
+    }
+
+    @KormaExperimental
+    fun scaledPoints(scale: Double): Line {
+        val dx = this.dx
+        val dy = this.dy
+        return Line(x0 - dx * scale, y0 - dy * scale, x1 + dx * scale, y1 + dy * scale)
+    }
+
+    fun containsX(x: Double): Boolean = (x in x0..x1) || (x in x1..x0) || (almostEquals(x, x0)) || (almostEquals(x, x1))
+    fun containsY(y: Double): Boolean = (y in y0..y1) || (y in y1..y0) || (almostEquals(y, y0)) || (almostEquals(y, y1))
+    fun containsBoundsXY(x: Double, y: Double): Boolean = containsX(x) && containsY(y)
+
+    val angle: Angle get() = Angle.between(a, b)
+    val length: Double get() = Point.distance(a, b)
+    val lengthSquared: Double get() = Point.distanceSquared(a, b)
+
+    fun getLineIntersectionPoint(line: Line): Point? =
+        getIntersectXY(x0, y0, x1, y1, line.x0, line.y0, line.x1, line.y1)
+
+    fun getIntersectionPoint(line: Line): Point? = getSegmentIntersectionPoint(line)
+    fun getSegmentIntersectionPoint(line: Line): Point? {
+        val out = getIntersectXY(x0, y0, x1, y1, line.x0, line.y0, line.x1, line.y1)
+        if (out != null && this.containsBoundsXY(out.x, out.y) && line.containsBoundsXY(out.x, out.y)) return out
+        return null
+    }
+
+    fun intersectsLine(line: Line): Boolean = getLineIntersectionPoint(line) != null
+    fun intersects(line: Line): Boolean = intersectsSegment(line)
+    fun intersectsSegment(line: Line): Boolean = getSegmentIntersectionPoint(line) != null
+
+    override fun toString(): String = "Line($a, $b)"
+
+    val isNIL get() = a.x.isNaN()
+    fun isNaN(): Boolean = a.y.isNaN()
+
+    companion object {
+        val ZERO = Line(Point.ZERO, Point.ZERO)
+        val NaN = Line(Point.NaN, Point.NaN)
+        val NIL: Line get() = NaN
+
+        fun fromPointAndDirection(point: Point, direction: Point, scale: Double = 1.0): Line =
+            Line(point, point + direction * scale)
+        fun fromPointAngle(point: Point, angle: Angle, length: Double = 1.0): Line =
+            Line(point, Point.polar(angle, length))
+
+        fun length(Ax: Double, Ay: Double, Bx: Double, By: Double): Double = kotlin.math.hypot(Bx - Ax, By - Ay)
+
+        inline fun getIntersectXY(Ax: Double, Ay: Double, Bx: Double, By: Double, Cx: Double, Cy: Double, Dx: Double, Dy: Double): Point? {
+            val a1 = By - Ay
+            val b1 = Ax - Bx
+            val c1 = a1 * (Ax) + b1 * (Ay)
+            val a2 = Dy - Cy
+            val b2 = Cx - Dx
+            val c2 = a2 * (Cx) + b2 * (Cy)
+            val determinant = a1 * b2 - a2 * b1
+            if (determinant.isAlmostZero()) return null
+            val x = (b2 * c1 - b1 * c2) / determinant
+            val y = (a1 * c2 - a2 * c1) / determinant
+            //if (!x.isFinite() || !y.isFinite()) TODO()
+            return Point(x, y)
+        }
+
+        fun getIntersectXY(a: Point, b: Point, c: Point, d: Point): Point? =
+            getIntersectXY(a.x, a.y, b.x, b.y, c.x, c.y, d.x, d.y)
+    }
+}

math/src/main/java/com/icegps/math/geometry/Line3D.kt

@@ -0,0 +1,3 @@
+package com.icegps.math.geometry
+
+data class Line3D(val a: Vector3D, val b: Vector3D)

math/src/main/java/com/icegps/math/geometry/Margin.kt

@@ -0,0 +1,75 @@
+package com.icegps.math.geometry
+
+import com.icegps.math.*
+import com.icegps.number.*
+
+/**
+ * A [top], [right], [bottom], [left] pack with FixedShort (16-bit) in the range of +-3275.9 (3.3 integer digits + 1 decimal digit)
+ */
+data class Margin(
+    val top: Double,
+    val right: Double,
+    val bottom: Double,
+    val left: Double,
+) : IsAlmostEquals<Margin> {
+    companion object {
+        val ZERO = Margin(0.0, 0.0, 0.0, 0.0)
+
+        inline operator fun invoke(margin: Number): Margin = Margin(margin.toDouble(), margin.toDouble(), margin.toDouble(), margin.toDouble())
+        inline operator fun invoke(vertical: Number, horizontal: Number): Margin = Margin(vertical.toDouble(), horizontal.toDouble(), vertical.toDouble(), horizontal.toDouble())
+        inline operator fun invoke(top: Number, right: Number, bottom: Number, left: Number): Margin = Margin(top.toDouble(), right.toDouble(), bottom.toDouble(), left.toDouble())
+    }
+
+    constructor(vertical: Double, horizontal: Double) : this(vertical, horizontal, vertical, horizontal)
+    constructor(margin: Double) : this(margin, margin, margin, margin)
+
+    operator fun plus(other: Margin): Margin = Margin(top + other.top, right + other.right, bottom + other.bottom, left + other.left)
+    operator fun minus(other: Margin): Margin = Margin(top - other.top, right - other.right, bottom - other.bottom, left - other.left)
+
+    val isNotZero: Boolean get() = top != 0.0 || left != 0.0 || right != 0.0 || bottom != 0.0
+
+    override fun isAlmostEquals(other: Margin, epsilon: Double): Boolean =
+        this.left.isAlmostEquals(other.left, epsilon) &&
+            this.right.isAlmostEquals(other.right, epsilon) &&
+            this.top.isAlmostEquals(other.top, epsilon) &&
+            this.bottom.isAlmostEquals(other.bottom, epsilon)
+    fun isAlmostZero(epsilon: Double = 0.000001): Boolean = isAlmostEquals(ZERO, epsilon)
+
+    val leftPlusRight: Double get() = left + right
+    val topPlusBottom: Double get() = top + bottom
+
+    val horizontal: Double get() = (left + right) / 2
+    val vertical: Double get() = (top + bottom) / 2
+
+    override fun toString(): String = "Margin(top=${top.niceStr}, right=${right.niceStr}, bottom=${bottom.niceStr}, left=${left.niceStr})"
+}
+
+/**
+ * A [top], [right], [bottom], [left] pack with Int)
+ */
+data class MarginInt(
+    val top: Int,
+    val right: Int,
+    val bottom: Int,
+    val left: Int,
+) {
+    constructor(top: Short, right: Short, bottom: Short, left: Short) : this(top.toInt(), right.toInt(), bottom.toInt(), left.toInt())
+    constructor(vertical: Int, horizontal: Int) : this(vertical, horizontal, vertical, horizontal)
+    constructor(margin: Int) : this(margin, margin, margin, margin)
+
+    operator fun plus(other: MarginInt): MarginInt = MarginInt(top + other.top, right + other.right, bottom + other.bottom, left + other.left)
+    operator fun minus(other: MarginInt): MarginInt = MarginInt(top - other.top, right - other.right, bottom - other.bottom, left - other.left)
+
+    val isNotZero: Boolean get() = top != 0 || left != 0 || right != 0 || bottom != 0
+
+    val leftPlusRight: Int get() = left + right
+    val topPlusBottom: Int get() = top + bottom
+    val horizontal: Int get() = (left + right) / 2
+    val vertical: Int get() = (top + bottom) / 2
+
+    companion object {
+        val ZERO = MarginInt(0, 0, 0, 0)
+    }
+
+    override fun toString(): String = "MarginInt(top=${top}, right=${right}, bottom=${bottom}, left=${left})"
+}

math/src/main/java/com/icegps/math/geometry/Matrix.kt

@@ -0,0 +1,415 @@
+package com.icegps.math.geometry
+
+import com.icegps.math.*
+import com.icegps.math.interpolation.*
+import com.icegps.number.*
+import kotlin.math.*
+
+
+//@KormaValueApi
+//data class Matrix(
+//    val a: Float,
+//    val b: Float,
+//    val c: Float,
+//    val d: Float,
+//    val tx: Float,
+//    val ty: Float,
+//) {
+
+// a, b, c, d, tx and ty are BFloat21
+data class Matrix(
+    val a: Double, val b: Double, val c: Double, val d: Double,
+    val tx: Double = 0.0, val ty: Double = 0.0
+) : IsAlmostEquals<Matrix> {
+    //private val twobits: Int get() = data.twobits
+
+    //constructor() : this(1f, 0f, 0f, 1f, 0f, 0f)
+    constructor(a: Float, b: Float, c: Float, d: Float, tx: Float = 0f, ty: Float = 0f) :
+        this(a.toDouble(), b.toDouble(), c.toDouble(), d.toDouble(), tx.toDouble(), ty.toDouble())
+    constructor(a: Int, b: Int, c: Int, d: Int, tx: Int = 0, ty: Int = 0) :
+        this(a.toDouble(), b.toDouble(), c.toDouble(), d.toDouble(), tx.toDouble(), ty.toDouble())
+
+    operator fun times(other: Matrix): Matrix = Matrix.multiply(this, other)
+    operator fun times(scale: Double): Matrix = Matrix(a * scale, b * scale, c * scale, d * scale, tx * scale, ty * scale)
+    operator fun times(scale: Float): Matrix = times(scale.toDouble())
+
+    //val isNIL: Boolean get() = this == NIL
+    val isNIL: Boolean get() = this.a.isNaN()
+    val isNotNIL: Boolean get() = !isNIL
+    val isNaN: Boolean get() = isNIL
+    val isIdentity: Boolean get() = (a == 1.0 && b == 0.0 && c == 0.0 && d == 1.0 && tx == 0.0 && ty == 0.0)
+    //val isIdentity: Boolean get() = twobits == 1
+
+    val type: MatrixType get() {
+        val hasRotation = b != 0.0 || c != 0.0
+        val hasScale = a != 1.0 || d != 1.0
+        val hasTranslation = tx != 0.0 || ty != 0.0
+
+        return when {
+            hasRotation -> MatrixType.COMPLEX
+            hasScale && hasTranslation -> MatrixType.SCALE_TRANSLATE
+            hasScale -> MatrixType.SCALE
+            hasTranslation -> MatrixType.TRANSLATE
+            else -> MatrixType.IDENTITY
+        }
+    }
+
+    inline fun transform(p: Vector2F): Vector2F {
+        if (this.isNIL) return p
+        return Vector2F(
+            this.a * p.x + this.c * p.y + this.tx,
+            this.d * p.y + this.b * p.x + this.ty
+        )
+    }
+    inline fun transform(p: Vector2D): Vector2D {
+        if (this.isNIL) return p
+        return Vector2D(
+            transformX(p.x, p.y),
+            transformY(p.x, p.y),
+        )
+    }
+
+    @Deprecated("", ReplaceWith("transform(p).x")) fun transformX(p: Point): Double = transformX(p.x, p.y)
+    @Deprecated("", ReplaceWith("transform(p).y")) fun transformY(p: Point): Double = transformY(p.x, p.y)
+
+    @Deprecated("", ReplaceWith("transform(p).x")) fun transformX(x: Float, y: Float): Float = transformX(x.toDouble(), y.toDouble()).toFloat()
+    @Deprecated("", ReplaceWith("transform(p).y")) fun transformY(x: Float, y: Float): Float = transformY(x.toDouble(), y.toDouble()).toFloat()
+
+    @Deprecated("", ReplaceWith("transform(p).x")) fun transformX(x: Double, y: Double): Double = this.a * x + this.c * y + this.tx
+    @Deprecated("", ReplaceWith("transform(p).y")) fun transformY(x: Double, y: Double): Double = this.d * y + this.b * x + this.ty
+
+    @Deprecated("", ReplaceWith("transform(p).x")) fun transformX(x: Int, y: Int): Double = transformX(x.toDouble(), y.toDouble())
+    @Deprecated("", ReplaceWith("transform(p).y")) fun transformY(x: Int, y: Int): Double = transformY(x.toDouble(), y.toDouble())
+
+    fun deltaTransform(p: Vector2F): Vector2F = Vector2F((p.x * a) + (p.y * c), (p.x * b) + (p.y * d))
+    fun deltaTransform(p: Vector2D): Vector2D = Vector2D((p.x * a) + (p.y * c), (p.x * b) + (p.y * d))
+
+    fun rotated(angle: Angle): Matrix {
+        val cos = cos(angle)
+        val sin = sin(angle)
+
+        val a1 = this.a * cos - this.b * sin
+        val b = (this.a * sin + this.b * cos)
+        val a = a1
+
+        val c1 = this.c * cos - this.d * sin
+        val d = (this.c * sin + this.d * cos)
+        val c = c1
+
+        val tx1 = this.tx * cos - this.ty * sin
+        val ty = (this.tx * sin + this.ty * cos)
+        val tx = tx1
+
+        return Matrix(a, b, c, d, tx, ty)
+    }
+
+    fun skewed(skewX: Angle, skewY: Angle): Matrix {
+        val sinX = sin(skewX)
+        val cosX = cos(skewX)
+        val sinY = sin(skewY)
+        val cosY = cos(skewY)
+
+        return Matrix(
+            a * cosY - b * sinX,
+            a * sinY + b * cosX,
+            c * cosY - d * sinX,
+            c * sinY + d * cosX,
+            tx * cosY - ty * sinX,
+            tx * sinY + ty * cosX
+        )
+    }
+
+    fun scaled(scaleX: Int, scaleY: Int = scaleX): Matrix = scaled(scaleX.toDouble(), scaleY.toDouble())
+    fun scaled(scaleX: Float, scaleY: Float = scaleX): Matrix = scaled(scaleX.toDouble(), scaleY.toDouble())
+    fun scaled(scaleX: Double, scaleY: Double = scaleX): Matrix = Matrix(a * scaleX, b * scaleX, c * scaleY, d * scaleY, tx * scaleX, ty * scaleY)
+
+    fun prescaled(scaleX: Int, scaleY: Int = scaleX): Matrix = prescaled(scaleX.toDouble(), scaleY.toDouble())
+    fun prescaled(scaleX: Float, scaleY: Float = scaleX): Matrix = prescaled(scaleX.toDouble(), scaleY.toDouble())
+    fun prescaled(scaleX: Double, scaleY: Double = scaleX): Matrix = Matrix(a * scaleX, b * scaleX, c * scaleY, d * scaleY, tx, ty)
+
+    fun translated(delta: Point): Matrix = Matrix(a, b, c, d, tx + delta.x, ty + delta.y)
+    fun translated(x: Int, y: Int): Matrix = translated(Point(x, y))
+    fun translated(x: Float, y: Float): Matrix = translated(Point(x, y))
+    fun translated(x: Double, y: Double): Matrix = translated(Point(x, y))
+
+    fun pretranslated(delta: Point): Matrix = Matrix(a, b, c, d, tx + (a * delta.x + c * delta.y), ty + (b * delta.x + d * delta.y))
+    fun pretranslated(deltaX: Int, deltaY: Int): Matrix = pretranslated(Point(deltaX, deltaY))
+    fun pretranslated(deltaX: Float, deltaY: Float): Matrix = pretranslated(Point(deltaX, deltaY))
+    fun pretranslated(deltaX: Double, deltaY: Double): Matrix = pretranslated(Point(deltaX, deltaY))
+
+    fun prerotated(angle: Angle): Matrix = rotating(angle) * this
+    fun preskewed(skewX: Angle, skewY: Angle): Matrix = skewing(skewX, skewY) * this
+
+    fun premultiplied(m: Matrix): Matrix = m * this
+    fun multiplied(m: Matrix): Matrix = this * m
+
+    /** Transform point without translation */
+    fun deltaTransformPoint(p: Point): Point = Point((p.x * a) + (p.y * c), (p.x * b) + (p.y * d))
+
+    @Deprecated("", ReplaceWith("this")) fun clone(): Matrix = this
+
+    fun inverted(): Matrix {
+        if (this.isNIL) return Matrix.IDENTITY
+        val m = this
+        val norm = m.a * m.d - m.b * m.c
+
+        return when (norm) {
+            0.0 -> Matrix(0.0, 0.0, 0.0, 0.0, -m.tx, -m.ty)
+            else -> {
+                val inorm = 1.0 / norm
+                val d = m.a * inorm
+                val a = m.d * inorm
+                val b = m.b * -inorm
+                val c = m.c * -inorm
+                Matrix(a, b, c, d, -a * m.tx - c * m.ty, -b * m.tx - d * m.ty)
+            }
+        }
+    }
+
+    fun toTransform(): MatrixTransform = decompose()
+    fun decompose(): MatrixTransform = MatrixTransform.fromMatrix(this)
+
+    fun toArray(value: DoubleArray, offset: Int = 0) {
+        value[offset + 0] = a
+        value[offset + 1] = b
+        value[offset + 2] = c
+        value[offset + 3] = d
+        value[offset + 4] = tx
+        value[offset + 5] = ty
+    }
+
+    fun toArray(value: FloatArray, offset: Int = 0) {
+        value[offset + 0] = a.toFloat()
+        value[offset + 1] = b.toFloat()
+        value[offset + 2] = c.toFloat()
+        value[offset + 3] = d.toFloat()
+        value[offset + 4] = tx.toFloat()
+        value[offset + 5] = ty.toFloat()
+    }
+
+    override fun toString(): String = "Matrix(${a.niceStr}, ${b.niceStr}, ${c.niceStr}, ${d.niceStr}, ${tx.niceStr}, ${ty.niceStr})"
+
+    override fun isAlmostEquals(other: Matrix, epsilon: Double): Boolean = isAlmostEquals(this, other, epsilon)
+    fun isAlmostIdentity(epsilon: Double = 0.00001): Boolean = isAlmostEquals(this, IDENTITY, epsilon)
+
+    // @TODO: Is this order correct?
+    fun preconcated(other: Matrix): Matrix = this * other
+
+    companion object {
+        val IDENTITY = Matrix(1.0, 0.0, 0.0, 1.0, 0.0, 0.0)
+        val NIL = Matrix(Double.NaN, Double.NaN, Double.NaN, Double.NaN, Double.NaN, Double.NaN)
+        val NaN = NIL
+
+        //@Deprecated("", ReplaceWith("com.icegps.math.geometry.Matrix.IDENTITY", "com.icegps.math.geometry.Matrix"))
+        operator fun invoke(): Matrix = IDENTITY
+
+        fun isAlmostEquals(a: Matrix, b: Matrix, epsilon: Double = 0.00001): Boolean =
+            a.tx.isAlmostEquals(b.tx, epsilon)
+                && a.ty.isAlmostEquals(b.ty, epsilon)
+                && a.a.isAlmostEquals(b.a, epsilon)
+                && a.b.isAlmostEquals(b.b, epsilon)
+                && a.c.isAlmostEquals(b.c, epsilon)
+                && a.d.isAlmostEquals(b.d, epsilon)
+
+        fun multiply(l: Matrix, r: Matrix): Matrix {
+            if (l.isNIL) return r
+            if (r.isNIL) return l
+            return Matrix(
+                l.a * r.a + l.b * r.c,
+                l.a * r.b + l.b * r.d,
+                l.c * r.a + l.d * r.c,
+                l.c * r.b + l.d * r.d,
+                l.tx * r.a + l.ty * r.c + r.tx,
+                l.tx * r.b + l.ty * r.d + r.ty
+            )
+        }
+
+        fun translating(delta: Point): Matrix = Matrix.IDENTITY.copy(tx = delta.x, ty = delta.y)
+        fun rotating(angle: Angle): Matrix = Matrix.IDENTITY.rotated(angle)
+        fun skewing(skewX: Angle, skewY: Angle): Matrix = Matrix.IDENTITY.skewed(skewX, skewY)
+
+        fun fromArray(value: FloatArray, offset: Int = 0): Matrix = Matrix(
+            value[offset + 0], value[offset + 1], value[offset + 2],
+            value[offset + 3], value[offset + 4], value[offset + 5]
+        )
+
+        fun fromArray(value: DoubleArray, offset: Int = 0): Matrix = Matrix(
+            value[offset + 0], value[offset + 1], value[offset + 2],
+            value[offset + 3], value[offset + 4], value[offset + 5]
+        )
+
+        fun fromTransform(
+            transform: MatrixTransform,
+            pivotX: Double = 0.0,
+            pivotY: Double = 0.0,
+        ): Matrix = fromTransform(
+            transform.x,
+            transform.y,
+            transform.rotation,
+            transform.scaleX,
+            transform.scaleY,
+            transform.skewX,
+            transform.skewY,
+            pivotX,
+            pivotY,
+        )
+
+        fun fromTransform(
+            x: Double,
+            y: Double,
+            rotation: Angle = Angle.ZERO,
+            scaleX: Double = 1.0,
+            scaleY: Double = 1.0,
+            skewX: Angle = Angle.ZERO,
+            skewY: Angle = Angle.ZERO,
+            pivotX: Double = 0.0,
+            pivotY: Double = 0.0,
+        ): Matrix {
+            // +0.0 drops the negative -0.0
+            val a = cos(rotation + skewY) * scaleX + 0f
+            val b = sin(rotation + skewY) * scaleX + 0f
+            val c = -sin(rotation - skewX) * scaleY + 0f
+            val d = cos(rotation - skewX) * scaleY + 0f
+            val tx: Double
+            val ty: Double
+
+            if (pivotX == 0.0 && pivotY == 0.0) {
+                tx = x
+                ty = y
+            } else {
+                tx = x - ((pivotX * a) + (pivotY * c))
+                ty = y - ((pivotX * b) + (pivotY * d))
+            }
+            return Matrix(a, b, c, d, tx, ty)
+        }
+
+        fun transform(a: Float, b: Float, c: Float, d: Float, tx: Float, ty: Float, p: Point): Point = Point(
+            a * p.x + c * p.y + tx,
+            d * p.y + b * p.x + ty
+        )
+
+        fun interpolated(l: Matrix, r: Matrix, ratio: Ratio): Matrix = Matrix(
+            ratio.interpolate(l.a, r.a),
+            ratio.interpolate(l.b, r.b),
+            ratio.interpolate(l.c, r.c),
+            ratio.interpolate(l.d, r.d),
+            ratio.interpolate(l.tx, r.tx),
+            ratio.interpolate(l.ty, r.ty),
+        )
+    }
+}
+
+//@KormaValueApi
+data class MatrixTransform(
+    val x: Double = 0.0, val y: Double = 0.0,
+    val scaleX: Double = 1.0, val scaleY: Double = 1.0,
+    val skewX: Angle = Angle.ZERO, val skewY: Angle = Angle.ZERO,
+    val rotation: Angle = Angle.ZERO
+) : IsAlmostEquals<MatrixTransform> {
+
+    override fun toString(): String = "MatrixTransform(x=${x.niceStr}, y=${y.niceStr}, scaleX=${scaleX}, scaleY=${scaleY}, skewX=${skewX}, skewY=${skewY}, rotation=${rotation})"
+
+    constructor() : this(0.0, 0.0, 1.0, 1.0, Angle.ZERO, Angle.ZERO, Angle.ZERO)
+    constructor(
+        x: Float, y: Float,
+        scaleX: Float, scaleY: Float,
+        skewX: Angle, skewY: Angle,
+        rotation: Angle
+    ) : this(x.toDouble(), y.toDouble(), scaleX.toDouble(), scaleY.toDouble(), skewX, skewY, rotation)
+
+    companion object {
+        val IDENTITY = MatrixTransform(0.0, 0.0, 1.0, 1.0, Angle.ZERO, Angle.ZERO, Angle.ZERO)
+
+        fun fromMatrix(matrix: Matrix, pivotX: Double = 0.0, pivotY: Double = 0.0): MatrixTransform {
+            val a = matrix.a
+            val b = matrix.b
+            val c = matrix.c
+            val d = matrix.d
+
+            val skewX = -atan2(-c, d)
+            val skewY = atan2(b, a)
+
+            val delta = abs(skewX + skewY)
+
+            val trotation: Angle
+            val tskewX: Angle
+            val tskewY: Angle
+            val tx: Double
+            val ty: Double
+
+            if (delta < 0.001f || abs((PI * 2) - delta) < 0.001f) {
+                trotation = skewY.radians
+                tskewX = 0.0.radians
+                tskewY = 0.0.radians
+            } else {
+                trotation = 0.radians
+                tskewX = skewX.radians
+                tskewY = skewY.radians
+            }
+
+            val tscaleX = hypot(a, b)
+            val tscaleY = hypot(c, d)
+
+            if (pivotX == 0.0 && pivotY == 0.0) {
+                tx = matrix.tx
+                ty = matrix.ty
+            } else {
+                tx = matrix.tx + ((pivotX * a) + (pivotY * c));
+                ty = matrix.ty + ((pivotX * b) + (pivotY * d));
+            }
+            return MatrixTransform(tx, ty, tscaleX, tscaleY, tskewX, tskewY, trotation)
+        }
+
+        fun interpolated(l: MatrixTransform, r: MatrixTransform, ratio: Ratio): MatrixTransform = MatrixTransform(
+            ratio.toRatio().interpolate(l.x, r.x),
+            ratio.toRatio().interpolate(l.y, r.y),
+            ratio.toRatio().interpolate(l.scaleX, r.scaleX),
+            ratio.toRatio().interpolate(l.scaleY, r.scaleY),
+            ratio.toRatio().interpolateAngleDenormalized(l.skewX, r.skewX),
+            ratio.toRatio().interpolateAngleDenormalized(l.skewY, r.skewY),
+            ratio.toRatio().interpolateAngleDenormalized(l.rotation, r.rotation),
+        )
+
+        fun isAlmostEquals(a: MatrixTransform, b: MatrixTransform, epsilon: Double = 0.000001): Boolean =
+            a.x.isAlmostEquals(b.x, epsilon)
+                && a.y.isAlmostEquals(b.y, epsilon)
+                && a.scaleX.isAlmostEquals(b.scaleX, epsilon)
+                && a.scaleY.isAlmostEquals(b.scaleY, epsilon)
+                && a.skewX.isAlmostEquals(b.skewX, epsilon)
+                && a.skewY.isAlmostEquals(b.skewY, epsilon)
+                && a.rotation.isAlmostEquals(b.rotation, epsilon)
+    }
+
+    override fun isAlmostEquals(other: MatrixTransform, epsilon: Double): Boolean = isAlmostEquals(this, other, epsilon)
+
+    val scaleAvg: Double get() = (scaleX + scaleY) * 0.5
+
+    fun toMatrix(pivotX: Double = 0.0, pivotY: Double = 0.0): Matrix = Matrix.fromTransform(this, pivotX, pivotY)
+
+    operator fun plus(that: MatrixTransform): MatrixTransform = MatrixTransform(
+        x + that.x, y + that.y,
+        scaleX * that.scaleX, scaleY * that.scaleY,
+        skewX + that.skewX, skewY + that.skewY,
+        rotation + that.rotation,
+    )
+    operator fun minus(that: MatrixTransform): MatrixTransform = MatrixTransform(
+        x - that.x, y - that.y,
+        scaleX / that.scaleX, scaleY / that.scaleY,
+        skewX - that.skewX, skewY - that.skewY,
+        rotation - that.rotation,
+    )
+}
+
+class MatrixComputed(val matrix: Matrix, val transform: MatrixTransform) {
+    companion object;
+    constructor(matrix: Matrix) : this(matrix, MatrixTransform.fromMatrix(matrix))
+    constructor(transform: MatrixTransform) : this(transform.toMatrix(), transform)
+}
+
+enum class MatrixType(val id: Int, val hasRotation: Boolean, val hasScale: Boolean, val hasTranslation: Boolean) {
+    IDENTITY(1, hasRotation = false, hasScale = false, hasTranslation = false),
+    TRANSLATE(2, hasRotation = false, hasScale = false, hasTranslation = true),
+    SCALE(3, hasRotation = false, hasScale = true, hasTranslation = false),
+    SCALE_TRANSLATE(4, hasRotation = false, hasScale = true, hasTranslation = true),
+    COMPLEX(5, hasRotation = true, hasScale = true, hasTranslation = true);
+}

math/src/main/java/com/icegps/math/geometry/Matrix3.kt

@@ -0,0 +1,237 @@
+@file:Suppress("NOTHING_TO_INLINE")
+
+package com.icegps.math.geometry
+
+import com.icegps.math.*
+import kotlin.math.*
+
+/**
+ * Useful for representing rotations and scales.
+ */
+data class Matrix3 private constructor(
+    internal val data: FloatArray,
+) : IsAlmostEqualsF<Matrix3> {
+    override fun equals(other: Any?): Boolean = other is Matrix3 && this.data.contentEquals(other.data)
+    override fun hashCode(): Int = data.contentHashCode()
+
+    private constructor(
+        v00: Float, v10: Float, v20: Float,
+        v01: Float, v11: Float, v21: Float,
+        v02: Float, v12: Float, v22: Float,
+    ) : this(
+        floatArrayOf(
+            v00, v10, v20,
+            v01, v11, v21,
+            v02, v12, v22,
+        )
+    )
+
+    init {
+        check(data.size == 9)
+    }
+
+    val v00: Float get() = data[0]
+    val v10: Float get() = data[1]
+    val v20: Float get() = data[2]
+    val v01: Float get() = data[3]
+    val v11: Float get() = data[4]
+    val v21: Float get() = data[5]
+    val v02: Float get() = data[6]
+    val v12: Float get() = data[7]
+    val v22: Float get() = data[8]
+
+    val c0: Vector3F get() = Vector3F.fromArray(data, 0)
+    val c1: Vector3F get() = Vector3F.fromArray(data, 3)
+    val c2: Vector3F get() = Vector3F.fromArray(data, 6)
+    fun c(column: Int): Vector3F {
+        if (column < 0 || column >= 3) error("Invalid column $column")
+        return Vector3F.fromArray(data, column * 3)
+    }
+
+    val r0: Vector3F get() = Vector3F(v00, v01, v02)
+    val r1: Vector3F get() = Vector3F(v10, v11, v12)
+    val r2: Vector3F get() = Vector3F(v20, v21, v22)
+
+    fun v(index: Int): Float = data[index]
+
+    fun r(row: Int): Vector3F = when (row) {
+        0 -> r0
+        1 -> r1
+        2 -> r2
+        else -> error("Invalid row $row")
+    }
+
+    operator fun get(row: Int, column: Int): Float {
+        if (column !in 0..2 || row !in 0..2) error("Invalid index $row,$column")
+        return data[row * 3 + column]
+    }
+
+    fun transform(v: Vector3F): Vector3F = Vector3F(r0.dot(v), r1.dot(v), r2.dot(v))
+
+    operator fun unaryMinus(): Matrix3 = Matrix3(
+        -v00, -v10, -v20,
+        -v01, -v11, -v21,
+        -v02, -v12, -v22,
+    )
+    operator fun unaryPlus(): Matrix3 = this
+
+    operator fun minus(other: Matrix3): Matrix3 = Matrix3(
+        v00 - other.v00, v10 - other.v10, v20 - other.v20,
+        v01 - other.v01, v11 - other.v11, v21 - other.v21,
+        v02 - other.v02, v12 - other.v12, v22 - other.v22,
+    )
+    operator fun plus(other: Matrix3): Matrix3 = Matrix3(
+        v00 + other.v00, v10 + other.v10, v20 + other.v20,
+        v01 + other.v01, v11 + other.v11, v21 + other.v21,
+        v02 + other.v02, v12 + other.v12, v22 + other.v22,
+    )
+
+    operator fun times(other: Matrix3): Matrix3 = Matrix3.multiply(this, other)
+    operator fun times(scale: Float): Matrix3 = Matrix3(
+        v00 * scale, v10 * scale, v20 * scale,
+        v01 * scale, v11 * scale, v21 * scale,
+        v02 * scale, v12 * scale, v22 * scale,
+    )
+    operator fun div(scale: Float): Matrix3 = this * (1f / scale)
+
+    fun inv(): Matrix3 = inverted()
+
+    val determinant: Float get() = v00 * (v11 * v22 - v21 * v12) -
+        v01 * (v10 * v22 - v12 * v20) +
+        v02 * (v10 * v21 - v11 * v20)
+
+    fun inverted(): Matrix3 {
+        val determinant = this.determinant
+
+        if (determinant == 0.0f) throw ArithmeticException("Matrix is not invertible")
+
+        val invDet = 1.0f / determinant
+
+        return fromRows(
+            (v11 * v22 - v21 * v12) * invDet,
+            (v02 * v21 - v01 * v22) * invDet,
+            (v01 * v12 - v02 * v11) * invDet,
+            (v12 * v20 - v10 * v22) * invDet,
+            (v00 * v22 - v02 * v20) * invDet,
+            (v10 * v02 - v00 * v12) * invDet,
+            (v10 * v21 - v20 * v11) * invDet,
+            (v20 * v01 - v00 * v21) * invDet,
+            (v00 * v11 - v10 * v01) * invDet,
+        )
+    }
+
+    override fun toString(): String = buildString {
+        append("Matrix3(\n")
+        for (row in 0 until 3) {
+            append("  [ ")
+            for (col in 0 until 3) {
+                if (col != 0) append(", ")
+                val v = get(row, col)
+                if (floor(v) == v) append(v.toInt()) else append(v)
+            }
+            append(" ],\n")
+        }
+        append(")")
+    }
+
+    fun transposed(): Matrix3 = Matrix3.fromColumns(r0, r1, r2)
+
+    override fun isAlmostEquals(other: Matrix3, epsilon: Float): Boolean = c0.isAlmostEquals(other.c0, epsilon)
+        && c1.isAlmostEquals(other.c1, epsilon)
+        && c2.isAlmostEquals(other.c2, epsilon)
+
+    companion object {
+        const val M00 = 0
+        const val M10 = 1
+        const val M20 = 2
+
+        const val M01 = 3
+        const val M11 = 4
+        const val M21 = 5
+
+        const val M02 = 6
+        const val M12 = 7
+        const val M22 = 8
+
+        const val M03 = 9
+        const val M13 = 10
+        const val M23 = 11
+
+        val INDICES_BY_COLUMNS = intArrayOf(
+            M00, M10, M20,
+            M01, M11, M21,
+            M02, M12, M22,
+        )
+        val INDICES_BY_ROWS = intArrayOf(
+            M00, M01, M02,
+            M10, M11, M12,
+            M20, M21, M22,
+        )
+
+        val IDENTITY = Matrix3(
+            1f, 0f, 0f,
+            0f, 1f, 0f,
+            0f, 0f, 1f,
+        )
+
+        fun fromRows(
+            r0: Vector3F, r1: Vector3F, r2: Vector3F
+        ): Matrix3 = Matrix3(
+            r0.x, r1.x, r2.x,
+            r0.y, r1.y, r2.y,
+            r0.z, r1.z, r2.z,
+        )
+
+        fun fromColumns(
+            c0: Vector3F, c1: Vector3F, c2: Vector3F
+        ): Matrix3 = Matrix3(
+            c0.x, c0.y, c0.z,
+            c1.x, c1.y, c1.z,
+            c2.x, c2.y, c2.z,
+        )
+
+        fun fromColumns(
+            v00: Float, v10: Float, v20: Float,
+            v01: Float, v11: Float, v21: Float,
+            v02: Float, v12: Float, v22: Float,
+        ): Matrix3 = Matrix3(
+            v00, v10, v20,
+            v01, v11, v21,
+            v02, v12, v22,
+        )
+
+        fun fromRows(
+            v00: Float, v01: Float, v02: Float,
+            v10: Float, v11: Float, v12: Float,
+            v20: Float, v21: Float, v22: Float,
+        ): Matrix3 = Matrix3(
+            v00, v10, v20,
+            v01, v11, v21,
+            v02, v12, v22,
+        )
+
+        fun multiply(l: Matrix3, r: Matrix3): Matrix3 = Matrix3.fromRows(
+            (l.v00 * r.v00) + (l.v01 * r.v10) + (l.v02 * r.v20),
+            (l.v00 * r.v01) + (l.v01 * r.v11) + (l.v02 * r.v21),
+            (l.v00 * r.v02) + (l.v01 * r.v12) + (l.v02 * r.v22),
+
+            (l.v10 * r.v00) + (l.v11 * r.v10) + (l.v12 * r.v20),
+            (l.v10 * r.v01) + (l.v11 * r.v11) + (l.v12 * r.v21),
+            (l.v10 * r.v02) + (l.v11 * r.v12) + (l.v12 * r.v22),
+
+            (l.v20 * r.v00) + (l.v21 * r.v10) + (l.v22 * r.v20),
+            (l.v20 * r.v01) + (l.v21 * r.v11) + (l.v22 * r.v21),
+            (l.v20 * r.v02) + (l.v21 * r.v12) + (l.v22 * r.v22),
+        )
+    }
+}
+
+fun Matrix3.toMatrix4(): Matrix4 = Matrix4.fromRows(
+    v00, v01, v02, 0f,
+    v10, v11, v12, 0f,
+    v20, v21, v22, 0f,
+    0f, 0f, 0f, 1f,
+)
+
+fun Matrix3.toQuaternion(): Quaternion = Quaternion.fromRotationMatrix(this)
+

math/src/main/java/com/icegps/math/geometry/Matrix4.kt

@@ -0,0 +1,684 @@
+package com.icegps.math.geometry
+
+import com.icegps.math.*
+import com.icegps.math.geometry.Matrix4.*
+import kotlin.math.*
+
+
+// @TODO: WIP
+// @TODO: value class
+// Stored as four consecutive column vectors (effectively stored in column-major order) see https://en.wikipedia.org/wiki/Row-_and_column-major_order
+// v[Row][Column]
+//@KormaExperimental
+//@KormaValueApi
+//inline class Matrix4 private constructor(
+/**
+ * Useful for representing complete transforms: rotations, scales, translations, projections, etc.
+ */
+data class Matrix4 private constructor(
+    private val data: FloatArray,
+    //val c0: Vector4, val c1: Vector4, val c2: Vector4, val c3: Vector4,
+
+    //val v00: Float, val v10: Float, val v20: Float, val v30: Float,
+    //val v01: Float, val v11: Float, val v21: Float, val v31: Float,
+    //val v02: Float, val v12: Float, val v22: Float, val v32: Float,
+    //val v03: Float, val v13: Float, val v23: Float, val v33: Float,
+) : IsAlmostEqualsF<Matrix4> {
+    init {
+        check(data.size == 16)
+    }
+    val v00: Float get() = data[0]; val v10: Float get() = data[1]; val v20: Float get() = data[2]; val v30: Float get() = data[3]
+    val v01: Float get() = data[4]; val v11: Float get() = data[5]; val v21: Float get() = data[6]; val v31: Float get() = data[7]
+    val v02: Float get() = data[8]; val v12: Float get() = data[9]; val v22: Float get() = data[10]; val v32: Float get() = data[11]
+    val v03: Float get() = data[12]; val v13: Float get() = data[13]; val v23: Float get() = data[14]; val v33: Float get() = data[15]
+
+    override fun equals(other: Any?): Boolean = other is Matrix4 && this.data.contentEquals(other.data)
+    override fun hashCode(): Int = data.contentHashCode()
+
+    operator fun times(scale: Float): Matrix4 = Matrix4.fromColumns(c0 * scale, c1 * scale, c2 * scale, c3 * scale)
+    operator fun times(that: Matrix4): Matrix4 = Matrix4.multiply(this, that)
+
+    fun transformTransposed(v: Vector4F): Vector4F = Vector4F(c0.dot(v), c1.dot(v), c2.dot(v), c3.dot(v))
+    fun transform(v: Vector4F): Vector4F = Vector4F(r0.dot(v), r1.dot(v), r2.dot(v), r3.dot(v))
+    fun transform(v: Vector3F): Vector3F = transform(v.toVector4()).toVector3()
+
+    fun transposed(): Matrix4 = Matrix4.fromColumns(r0, r1, r2, r3)
+
+    val determinant: Float get() = 0f +
+        (v30 * v21 * v12 * v03) -
+        (v20 * v31 * v12 * v03) -
+        (v30 * v11 * v22 * v03) +
+        (v10 * v31 * v22 * v03) +
+        (v20 * v11 * v32 * v03) -
+        (v10 * v21 * v32 * v03) -
+        (v30 * v21 * v02 * v13) +
+        (v20 * v31 * v02 * v13) +
+        (v30 * v01 * v22 * v13) -
+        (v00 * v31 * v22 * v13) -
+        (v20 * v01 * v32 * v13) +
+        (v00 * v21 * v32 * v13) +
+        (v30 * v11 * v02 * v23) -
+        (v10 * v31 * v02 * v23) -
+        (v30 * v01 * v12 * v23) +
+        (v00 * v31 * v12 * v23) +
+        (v10 * v01 * v32 * v23) -
+        (v00 * v11 * v32 * v23) -
+        (v20 * v11 * v02 * v33) +
+        (v10 * v21 * v02 * v33) +
+        (v20 * v01 * v12 * v33) -
+        (v00 * v21 * v12 * v33) -
+        (v10 * v01 * v22 * v33) +
+        (v00 * v11 * v22 * v33)
+
+    // Use toTRS/decompose
+    //fun decomposeProjection(): Vector4 = c3
+    //fun decomposeTranslation(): Vector4 = r3.copy(w = 1f)
+    //fun decomposeScale(): Vector4 {
+    //    val x = r0.length3
+    //    val y = r1.length3
+    //    val z = r2.length3
+    //    return Vector4(x, y, z, 1f)
+    //}
+    fun decomposeRotation(rowNormalise: Boolean = true): Quaternion {
+        var v1 = this.r0
+        var v2 = this.r1
+        var v3 = this.r2
+        if (rowNormalise) {
+            v1 = v1.normalized()
+            v2 = v2.normalized()
+            v3 = v3.normalized()
+        }
+        val d: Float = 0.25f * (v1[0] + v2[1] + v3[2] + 1f)
+        val out: Vector4F
+        when {
+            d > 0f -> {
+                val num1: Float = sqrt(d)
+                val num2: Float = 1f / (4f * num1)
+                out = Vector4F(
+                    ((v2[2] - v3[1]) * num2),
+                    ((v3[0] - v1[2]) * num2),
+                    ((v1[1] - v2[0]) * num2),
+                    num1,
+                )
+            }
+            v1[0] > v2[1] && v1[0] > v3[2] -> {
+                val num1: Float = 2f * sqrt(1f + v1[0] - v2[1] - v3[2])
+                val num2: Float = 1f / num1
+                out = Vector4F(
+                    (0.25f * num1),
+                    ((v2[0] + v1[1]) * num2),
+                    ((v3[0] + v1[2]) * num2),
+                    ((v3[1] - v2[2]) * num2),
+                )
+            }
+            v2[1] > v3[2] -> {
+                val num5: Float = 2f * sqrt(1f + v2[1] - v1[0] - v3[2])
+                val num6: Float = 1f / num5
+                out = Vector4F(
+                    ((v2[0] + v1[1]) * num6),
+                    (0.25f * num5),
+                    ((v3[1] + v2[2]) * num6),
+                    ((v3[0] - v1[2]) * num6),
+                )
+            }
+            else -> {
+                val num7: Float = 2f * sqrt(1f + v3[2] - v1[0] - v2[1])
+                val num8: Float = 1f / num7
+                out = Vector4F(
+                    ((v3[0] + v1[2]) * num8),
+                    ((v3[1] + v2[2]) * num8),
+                    (0.25f * num7),
+                    ((v2[0] - v1[1]) * num8),
+                )
+            }
+        }
+        return Quaternion(out.normalized())
+    }
+
+    fun copyToColumns(out: FloatArray = FloatArray(16), offset: Int = 0): FloatArray {
+        this.data.copyInto(out, offset, 0, 16)
+        return out
+    }
+    fun copyToRows(out: FloatArray = FloatArray(16), offset: Int = 0): FloatArray {
+        this.r0.copyTo(out, offset + 0)
+        this.r1.copyTo(out, offset + 4)
+        this.r2.copyTo(out, offset + 8)
+        this.r3.copyTo(out, offset + 12)
+        return out
+    }
+
+    private constructor(
+        v00: Float, v10: Float, v20: Float, v30: Float,
+        v01: Float, v11: Float, v21: Float, v31: Float,
+        v02: Float, v12: Float, v22: Float, v32: Float,
+        v03: Float, v13: Float, v23: Float, v33: Float,
+    ) : this(floatArrayOf(
+        v00, v10, v20, v30,
+        v01, v11, v21, v31,
+        v02, v12, v22, v32,
+        v03, v13, v23, v33,
+    ))
+
+    constructor() : this(
+        1f, 0f, 0f, 0f,
+        0f, 1f, 0f, 0f,
+        0f, 0f, 1f, 0f,
+        0f, 0f, 0f, 1f,
+    )
+
+    val c0: Vector4F get() = Vector4F.fromArray(data, 0)
+    val c1: Vector4F get() = Vector4F.fromArray(data, 4)
+    val c2: Vector4F get() = Vector4F.fromArray(data, 8)
+    val c3: Vector4F get() = Vector4F.fromArray(data, 12)
+    fun c(column: Int): Vector4F {
+        if (column < 0 || column >= 4) error("Invalid column $column")
+        return Vector4F.fromArray(data, column * 4)
+    }
+
+    val r0: Vector4F get() = Vector4F(v00, v01, v02, v03)
+    val r1: Vector4F get() = Vector4F(v10, v11, v12, v13)
+    val r2: Vector4F get() = Vector4F(v20, v21, v22, v23)
+    val r3: Vector4F get() = Vector4F(v30, v31, v32, v33)
+
+    fun r(row: Int): Vector4F = when (row) {
+        0 -> r0
+        1 -> r1
+        2 -> r2
+        3 -> r3
+        else -> error("Invalid row $row")
+    }
+
+    operator fun get(row: Int, column: Int): Float {
+        if (column !in 0..3 || row !in 0..3) error("Invalid index $row,$column")
+        return data[row * 4 + column]
+    }
+
+    fun getAtIndex(index: Int): Float {
+        if (index !in data.indices) error("Invalid index $index")
+        return data[index]
+    }
+
+    override fun toString(): String = buildString {
+        append("Matrix4(\n")
+        for (row in 0 until 4) {
+            append("  [ ")
+            for (col in 0 until 4) {
+                if (col != 0) append(", ")
+                val v = get(row, col)
+                if (floor(v) == v) append(v.toInt()) else append(v)
+            }
+            append(" ],\n")
+        }
+        append(")")
+    }
+
+
+
+    fun translated(x: Float, y: Float, z: Float, w: Float = 1f): Matrix4 = this * Matrix4.translation(x, y, z, w)
+    fun translated(x: Double, y: Double, z: Double, w: Double = 1.0) = this.translated(x.toFloat(), y.toFloat(), z.toFloat(), w.toFloat())
+    fun translated(x: Int, y: Int, z: Int, w: Int = 1) = this.translated(x.toFloat(), y.toFloat(), z.toFloat(), w.toFloat())
+
+    fun rotated(angle: Angle, x: Float, y: Float, z: Float): Matrix4 = this * Matrix4.rotation(angle, x, y, z)
+    fun rotated(angle: Angle, x: Double, y: Double, z: Double): Matrix4 = this.rotated(angle, x.toFloat(), y.toFloat(), z.toFloat())
+    fun rotated(angle: Angle, x: Int, y: Int, z: Int): Matrix4 = this.rotated(angle, x.toFloat(), y.toFloat(), z.toFloat())
+
+    fun scaled(x: Float, y: Float, z: Float, w: Float = 1f): Matrix4 = this * Matrix4.scale(x, y, z, w)
+    fun scaled(x: Double, y: Double, z: Double, w: Double = 1.0): Matrix4 = this.scaled(x.toFloat(), y.toFloat(), z.toFloat(), w.toFloat())
+    fun scaled(x: Int, y: Int, z: Int, w: Int = 1): Matrix4 = this.scaled(x.toFloat(), y.toFloat(), z.toFloat(), w.toFloat())
+
+    fun rotated(quat: Quaternion): Matrix4 = this * quat.toMatrix()
+    fun rotated(euler: EulerRotation): Matrix4 = this * euler.toMatrix()
+    fun rotated(x: Angle, y: Angle, z: Angle): Matrix4 = rotated(x, 1f, 0f, 0f).rotated(y, 0f, 1f, 0f).rotated(z, 0f, 0f, 1f)
+
+    fun decompose(): TRS4 = toTRS()
+    fun toTRS(): TRS4 {
+        val det = determinant
+        val translation = Vector4F(v03, v13, v23, 1f)
+        val scale = Vector4F(Vector3F.length(v00, v10, v20) * det.sign, Vector3F.length(v01, v11, v21), Vector3F.length(v02, v12, v22), 1f)
+        val invSX = 1f / scale.x
+        val invSY = 1f / scale.y
+        val invSZ = 1f / scale.z
+        val rotation = Quaternion.fromRotationMatrix(Matrix4.fromRows(
+            v00 * invSX, v01 * invSY, v02 * invSZ, v03,
+            v10 * invSX, v11 * invSY, v12 * invSZ, v13,
+            v20 * invSX, v21 * invSY, v22 * invSZ, v23,
+            v30, v31, v32, v33
+        ))
+        return TRS4(translation, rotation, scale)
+    }
+
+    fun inverted(): Matrix4 {
+        val t11 = v12 * v23 * v31 - v13 * v22 * v31 + v13 * v21 * v32 - v11 * v23 * v32 - v12 * v21 * v33 + v11 * v22 * v33
+        val t12 = v03 * v22 * v31 - v02 * v23 * v31 - v03 * v21 * v32 + v01 * v23 * v32 + v02 * v21 * v33 - v01 * v22 * v33
+        val t13 = v02 * v13 * v31 - v03 * v12 * v31 + v03 * v11 * v32 - v01 * v13 * v32 - v02 * v11 * v33 + v01 * v12 * v33
+        val t14 = v03 * v12 * v21 - v02 * v13 * v21 - v03 * v11 * v22 + v01 * v13 * v22 + v02 * v11 * v23 - v01 * v12 * v23
+
+        val det = v00 * t11 + v10 * t12 + v20 * t13 + v30 * t14
+
+        if (det == 0f) {
+            println("Matrix doesn't have inverse")
+            return Matrix4.IDENTITY
+        }
+
+        val detInv = 1 / det
+
+        return Matrix4.fromRows(
+            t11 * detInv,
+            t12 * detInv,
+            t13 * detInv,
+            t14 * detInv,
+
+            (v13 * v22 * v30 - v12 * v23 * v30 - v13 * v20 * v32 + v10 * v23 * v32 + v12 * v20 * v33 - v10 * v22 * v33) * detInv,
+            (v02 * v23 * v30 - v03 * v22 * v30 + v03 * v20 * v32 - v00 * v23 * v32 - v02 * v20 * v33 + v00 * v22 * v33) * detInv,
+            (v03 * v12 * v30 - v02 * v13 * v30 - v03 * v10 * v32 + v00 * v13 * v32 + v02 * v10 * v33 - v00 * v12 * v33) * detInv,
+            (v02 * v13 * v20 - v03 * v12 * v20 + v03 * v10 * v22 - v00 * v13 * v22 - v02 * v10 * v23 + v00 * v12 * v23) * detInv,
+
+            (v11 * v23 * v30 - v13 * v21 * v30 + v13 * v20 * v31 - v10 * v23 * v31 - v11 * v20 * v33 + v10 * v21 * v33) * detInv,
+            (v03 * v21 * v30 - v01 * v23 * v30 - v03 * v20 * v31 + v00 * v23 * v31 + v01 * v20 * v33 - v00 * v21 * v33) * detInv,
+            (v01 * v13 * v30 - v03 * v11 * v30 + v03 * v10 * v31 - v00 * v13 * v31 - v01 * v10 * v33 + v00 * v11 * v33) * detInv,
+            (v03 * v11 * v20 - v01 * v13 * v20 - v03 * v10 * v21 + v00 * v13 * v21 + v01 * v10 * v23 - v00 * v11 * v23) * detInv,
+
+            (v12 * v21 * v30 - v11 * v22 * v30 - v12 * v20 * v31 + v10 * v22 * v31 + v11 * v20 * v32 - v10 * v21 * v32) * detInv,
+            (v01 * v22 * v30 - v02 * v21 * v30 + v02 * v20 * v31 - v00 * v22 * v31 - v01 * v20 * v32 + v00 * v21 * v32) * detInv,
+            (v02 * v11 * v30 - v01 * v12 * v30 - v02 * v10 * v31 + v00 * v12 * v31 + v01 * v10 * v32 - v00 * v11 * v32) * detInv,
+            (v01 * v12 * v20 - v02 * v11 * v20 + v02 * v10 * v21 - v00 * v12 * v21 - v01 * v10 * v22 + v00 * v11 * v22) * detInv
+        )
+    }
+
+    override fun isAlmostEquals(other: Matrix4, epsilon: Float): Boolean =
+        c0.isAlmostEquals(other.c0, epsilon) &&
+            c1.isAlmostEquals(other.c1, epsilon) &&
+            c2.isAlmostEquals(other.c2, epsilon) &&
+            c3.isAlmostEquals(other.c3, epsilon)
+
+    companion object {
+        const val M00 = 0
+        const val M10 = 1
+        const val M20 = 2
+        const val M30 = 3
+
+        const val M01 = 4
+        const val M11 = 5
+        const val M21 = 6
+        const val M31 = 7
+
+        const val M02 = 8
+        const val M12 = 9
+        const val M22 = 10
+        const val M32 = 11
+
+        const val M03 = 12
+        const val M13 = 13
+        const val M23 = 14
+        const val M33 = 15
+
+        val INDICES_BY_COLUMNS_4x4 = intArrayOf(
+            M00, M10, M20, M30,
+            M01, M11, M21, M31,
+            M02, M12, M22, M32,
+            M03, M13, M23, M33,
+        )
+        val INDICES_BY_ROWS_4x4 = intArrayOf(
+            M00, M01, M02, M03,
+            M10, M11, M12, M13,
+            M20, M21, M22, M23,
+            M30, M31, M32, M33,
+        )
+        val INDICES_BY_COLUMNS_3x3 = intArrayOf(
+            M00, M10, M20,
+            M01, M11, M21,
+            M02, M12, M22,
+        )
+        val INDICES_BY_ROWS_3x3 = intArrayOf(
+            M00, M01, M02,
+            M10, M11, M12,
+            M20, M21, M22,
+        )
+
+        val IDENTITY = Matrix4()
+
+        fun fromColumns(
+            c0: Vector4F, c1: Vector4F, c2: Vector4F, c3: Vector4F
+        ): Matrix4 = Matrix4(
+            c0.x, c0.y, c0.z, c0.w,
+            c1.x, c1.y, c1.z, c1.w,
+            c2.x, c2.y, c2.z, c2.w,
+            c3.x, c3.y, c3.z, c3.w,
+        )
+
+        fun fromColumns(v: FloatArray, offset: Int = 0): Matrix4 = Matrix4.fromColumns(
+            v[offset + 0], v[offset + 1], v[offset + 2], v[offset + 3],
+            v[offset + 4], v[offset + 5], v[offset + 6], v[offset + 7],
+            v[offset + 8], v[offset + 9], v[offset + 10], v[offset + 11],
+            v[offset + 12], v[offset + 13], v[offset + 14], v[offset + 15],
+        )
+
+        fun fromRows(v: FloatArray, offset: Int = 0): Matrix4 = Matrix4.fromRows(
+            v[offset + 0], v[offset + 1], v[offset + 2], v[offset + 3],
+            v[offset + 4], v[offset + 5], v[offset + 6], v[offset + 7],
+            v[offset + 8], v[offset + 9], v[offset + 10], v[offset + 11],
+            v[offset + 12], v[offset + 13], v[offset + 14], v[offset + 15],
+        )
+
+        fun fromRows(
+            r0: Vector4F, r1: Vector4F, r2: Vector4F, r3: Vector4F
+        ): Matrix4 = Matrix4(
+            r0.x, r1.x, r2.x, r3.x,
+            r0.y, r1.y, r2.y, r3.y,
+            r0.z, r1.z, r2.z, r3.z,
+            r0.w, r1.w, r2.w, r3.w,
+        )
+
+        fun fromColumns(
+            v00: Float, v10: Float, v20: Float, v30: Float,
+            v01: Float, v11: Float, v21: Float, v31: Float,
+            v02: Float, v12: Float, v22: Float, v32: Float,
+            v03: Float, v13: Float, v23: Float, v33: Float,
+        ): Matrix4 = Matrix4(
+            v00, v10, v20, v30,
+            v01, v11, v21, v31,
+            v02, v12, v22, v32,
+            v03, v13, v23, v33,
+        )
+
+        fun fromRows(
+            v00: Float, v01: Float, v02: Float, v03: Float,
+            v10: Float, v11: Float, v12: Float, v13: Float,
+            v20: Float, v21: Float, v22: Float, v23: Float,
+            v30: Float, v31: Float, v32: Float, v33: Float,
+        ): Matrix4 = Matrix4(
+            v00, v10, v20, v30,
+            v01, v11, v21, v31,
+            v02, v12, v22, v32,
+            v03, v13, v23, v33,
+        )
+
+        fun fromRows3x3(
+            a00: Float, a01: Float, a02: Float,
+            a10: Float, a11: Float, a12: Float,
+            a20: Float, a21: Float, a22: Float
+        ): Matrix4 = Matrix4.fromRows(
+            a00, a01, a02, 0f,
+            a10, a11, a12, 0f,
+            a20, a21, a22, 0f,
+            0f, 0f, 0f, 1f,
+        )
+
+        fun fromColumns3x3(
+            a00: Float, a10: Float, a20: Float,
+            a01: Float, a11: Float, a21: Float,
+            a02: Float, a12: Float, a22: Float
+        ): Matrix4 = Matrix4.fromColumns(
+            a00, a10, a20, 0f,
+            a01, a11, a21, 0f,
+            a02, a12, a22, 0f,
+            0f, 0f, 0f, 1f,
+        )
+
+        fun fromTRS(trs: TRS4): Matrix4 = fromTRS(trs.translation, trs.rotation, trs.scale)
+        fun fromTRS(translation: Vector4F, rotation: Quaternion, scale: Vector4F): Matrix4 {
+            val rx = rotation.x
+            val ry = rotation.y
+            val rz = rotation.z
+            val rw = rotation.w
+
+            val xt = rx + rx
+            val yt = ry + ry
+            val zt = rz + rz
+
+            val xx = rx * xt
+            val xy = rx * yt
+            val xz = rx * zt
+
+            val yy = ry * yt
+            val yz = ry * zt
+            val zz = rz * zt
+
+            val wx = rw * xt
+            val wy = rw * yt
+            val wz = rw * zt
+
+            return Matrix4.fromRows(
+                ((1 - (yy + zz)) * scale.x), ((xy - wz) * scale.y), ((xz + wy) * scale.z), translation.x,
+                ((xy + wz) * scale.x), ((1 - (xx + zz)) * scale.y), ((yz - wx) * scale.z), translation.y,
+                ((xz - wy) * scale.x), ((yz + wx) * scale.y), ((1 - (xx + yy)) * scale.z), translation.z,
+                0f, 0f, 0f, 1f
+            )
+        }
+
+        fun translation(x: Float, y: Float, z: Float, w: Float = 1f): Matrix4 = Matrix4.fromRows(
+            1f, 0f, 0f, x,
+            0f, 1f, 0f, y,
+            0f, 0f, 1f, z,
+            0f, 0f, 0f, w
+        )
+        fun translation(x: Double, y: Double, z: Double, w: Double = 1.0): Matrix4 = translation(x.toFloat(), y.toFloat(), z.toFloat(), w.toFloat())
+        fun translation(x: Int, y: Int, z: Int, w: Int = 1): Matrix4 = translation(x.toFloat(), y.toFloat(), z.toFloat(), w.toFloat())
+
+        fun scale(x: Float, y: Float, z: Float, w: Float = 1f): Matrix4 = Matrix4.fromRows(
+            x, 0f, 0f, 0f,
+            0f, y, 0f, 0f,
+            0f, 0f, z, 0f,
+            0f, 0f, 0f, w
+        )
+        fun scale(x: Double, y: Double, z: Double, w: Double = 1.0): Matrix4 = scale(x.toFloat(), y.toFloat(), z.toFloat(), w.toFloat())
+        fun scale(x: Int, y: Int, z: Int, w: Int = 1): Matrix4 = scale(x.toFloat(), y.toFloat(), z.toFloat(), w.toFloat())
+
+        fun shear(x: Float, y: Float, z: Float): Matrix4 = fromRows(
+            1f, y, z, 0f,
+            x, 1f, z, 0f,
+            x, y, 1f, 0f,
+            0f, 0f, 0f, 1f
+        )
+        fun shear(x: Double, y: Double, z: Double): Matrix4 = shear(x.toFloat(), y.toFloat(), z.toFloat())
+        fun shear(x: Int, y: Int, z: Int): Matrix4 = shear(x.toFloat(), y.toFloat(), z.toFloat())
+
+        fun rotationX(angle: Angle): Matrix4 {
+            val c = angle.cosine.toFloat()
+            val s = angle.sine.toFloat()
+            return Matrix4.fromRows(
+                1f, 0f, 0f, 0f,
+                0f, c, -s, 0f,
+                0f, s, c, 0f,
+                0f, 0f, 0f, 1f
+            )
+        }
+
+        fun rotationY(angle: Angle): Matrix4 {
+            val c = angle.cosine.toFloat()
+            val s = angle.sine.toFloat()
+            return Matrix4.fromRows(
+                c, 0f, s, 0f,
+                0f, 1f, 0f, 0f,
+                -s, 0f, c, 0f,
+                0f, 0f, 0f, 1f
+            )
+        }
+
+        fun rotationZ(angle: Angle): Matrix4 {
+            val c = angle.cosine.toFloat()
+            val s = angle.sine.toFloat()
+            return Matrix4.fromRows(
+                c, -s, 0f, 0f,
+                s, c, 0f, 0f,
+                0f, 0f, 1f, 0f,
+                0f, 0f, 0f, 1f
+            )
+        }
+
+        fun rotation(angle: Angle, x: Float, y: Float, z: Float): Matrix4 {
+            val mag = sqrt(x * x + y * y + z * z)
+            val norm = 1f / mag
+
+            val nx = x * norm
+            val ny = y * norm
+            val nz = z * norm
+            val c = angle.cosine.toFloat()
+            val s = angle.sine.toFloat()
+            val t = 1 - c
+            val tx = t * nx
+            val ty = t * ny
+
+            return Matrix4.fromRows(
+                tx * nx + c, tx * ny - s * nz, tx * nz + s * ny, 0f,
+                tx * ny + s * nz, ty * ny + c, ty * nz - s * nx, 0f,
+                tx * nz - s * ny, ty * nz + s * nx, t * nz * nz + c, 0f,
+                0f, 0f, 0f, 1f
+            )
+        }
+        fun rotation(angle: Angle, direction: Vector3F): Matrix4 = rotation(angle, direction.x, direction.y, direction.z)
+        fun rotation(angle: Angle, x: Double, y: Double, z: Double): Matrix4 = rotation(angle, x.toFloat(), y.toFloat(), z.toFloat())
+        fun rotation(angle: Angle, x: Int, y: Int, z: Int): Matrix4 = rotation(angle, x.toFloat(), y.toFloat(), z.toFloat())
+
+        // @TODO: Use Vector4 operations, and use columns instead of rows for faster set
+        fun multiply(l: Matrix4, r: Matrix4): Matrix4 = Matrix4.fromRows(
+            (l.v00 * r.v00) + (l.v01 * r.v10) + (l.v02 * r.v20) + (l.v03 * r.v30),
+            (l.v00 * r.v01) + (l.v01 * r.v11) + (l.v02 * r.v21) + (l.v03 * r.v31),
+            (l.v00 * r.v02) + (l.v01 * r.v12) + (l.v02 * r.v22) + (l.v03 * r.v32),
+            (l.v00 * r.v03) + (l.v01 * r.v13) + (l.v02 * r.v23) + (l.v03 * r.v33),
+
+            (l.v10 * r.v00) + (l.v11 * r.v10) + (l.v12 * r.v20) + (l.v13 * r.v30),
+            (l.v10 * r.v01) + (l.v11 * r.v11) + (l.v12 * r.v21) + (l.v13 * r.v31),
+            (l.v10 * r.v02) + (l.v11 * r.v12) + (l.v12 * r.v22) + (l.v13 * r.v32),
+            (l.v10 * r.v03) + (l.v11 * r.v13) + (l.v12 * r.v23) + (l.v13 * r.v33),
+
+            (l.v20 * r.v00) + (l.v21 * r.v10) + (l.v22 * r.v20) + (l.v23 * r.v30),
+            (l.v20 * r.v01) + (l.v21 * r.v11) + (l.v22 * r.v21) + (l.v23 * r.v31),
+            (l.v20 * r.v02) + (l.v21 * r.v12) + (l.v22 * r.v22) + (l.v23 * r.v32),
+            (l.v20 * r.v03) + (l.v21 * r.v13) + (l.v22 * r.v23) + (l.v23 * r.v33),
+
+            (l.v30 * r.v00) + (l.v31 * r.v10) + (l.v32 * r.v20) + (l.v33 * r.v30),
+            (l.v30 * r.v01) + (l.v31 * r.v11) + (l.v32 * r.v21) + (l.v33 * r.v31),
+            (l.v30 * r.v02) + (l.v31 * r.v12) + (l.v32 * r.v22) + (l.v33 * r.v32),
+            (l.v30 * r.v03) + (l.v31 * r.v13) + (l.v32 * r.v23) + (l.v33 * r.v33)
+        )
+
+        fun multiply(
+            lv00: Float, lv01: Float, lv02: Float, lv03: Float,
+            lv10: Float, lv11: Float, lv12: Float, lv13: Float,
+            lv20: Float, lv21: Float, lv22: Float, lv23: Float,
+            lv30: Float, lv31: Float, lv32: Float, lv33: Float,
+
+            rv00: Float, rv01: Float, rv02: Float, rv03: Float,
+            rv10: Float, rv11: Float, rv12: Float, rv13: Float,
+            rv20: Float, rv21: Float, rv22: Float, rv23: Float,
+            rv30: Float, rv31: Float, rv32: Float, rv33: Float,
+        ): Matrix4 = Matrix4.fromRows(
+            (lv00 * rv00) + (lv01 * rv10) + (lv02 * rv20) + (lv03 * rv30),
+            (lv00 * rv01) + (lv01 * rv11) + (lv02 * rv21) + (lv03 * rv31),
+            (lv00 * rv02) + (lv01 * rv12) + (lv02 * rv22) + (lv03 * rv32),
+            (lv00 * rv03) + (lv01 * rv13) + (lv02 * rv23) + (lv03 * rv33),
+
+            (lv10 * rv00) + (lv11 * rv10) + (lv12 * rv20) + (lv13 * rv30),
+            (lv10 * rv01) + (lv11 * rv11) + (lv12 * rv21) + (lv13 * rv31),
+            (lv10 * rv02) + (lv11 * rv12) + (lv12 * rv22) + (lv13 * rv32),
+            (lv10 * rv03) + (lv11 * rv13) + (lv12 * rv23) + (lv13 * rv33),
+
+            (lv20 * rv00) + (lv21 * rv10) + (lv22 * rv20) + (lv23 * rv30),
+            (lv20 * rv01) + (lv21 * rv11) + (lv22 * rv21) + (lv23 * rv31),
+            (lv20 * rv02) + (lv21 * rv12) + (lv22 * rv22) + (lv23 * rv32),
+            (lv20 * rv03) + (lv21 * rv13) + (lv22 * rv23) + (lv23 * rv33),
+
+            (lv30 * rv00) + (lv31 * rv10) + (lv32 * rv20) + (lv33 * rv30),
+            (lv30 * rv01) + (lv31 * rv11) + (lv32 * rv21) + (lv33 * rv31),
+            (lv30 * rv02) + (lv31 * rv12) + (lv32 * rv22) + (lv33 * rv32),
+            (lv30 * rv03) + (lv31 * rv13) + (lv32 * rv23) + (lv33 * rv33)
+        )
+
+        fun ortho(left: Float, right: Float, bottom: Float, top: Float, near: Float = 0f, far: Float = 1f): Matrix4 {
+            val sx = 2f / (right - left)
+            val sy = 2f / (top - bottom)
+            val sz = -2f / (far - near)
+
+            val tx = -(right + left) / (right - left)
+            val ty = -(top + bottom) / (top - bottom)
+            val tz = -(far + near) / (far - near)
+
+            return Matrix4.fromRows(
+                sx, 0f, 0f, tx,
+                0f, sy, 0f, ty,
+                0f, 0f, sz, tz,
+                0f, 0f, 0f, 1f
+            )
+        }
+        fun ortho(left: Double, right: Double, bottom: Double, top: Double, near: Double, far: Double): Matrix4 =
+            ortho(left.toFloat(), right.toFloat(), bottom.toFloat(), top.toFloat(), near.toFloat(), far.toFloat())
+        fun ortho(left: Int, right: Int, bottom: Int, top: Int, near: Int, far: Int): Matrix4 =
+            ortho(left.toFloat(), right.toFloat(), bottom.toFloat(), top.toFloat(), near.toFloat(), far.toFloat())
+
+        fun frustum(left: Float, right: Float, bottom: Float, top: Float, zNear: Float = 0f, zFar: Float = 1f): Matrix4 {
+            if (zNear <= 0.0f || zFar <= zNear) {
+                throw Exception("Error: Required zNear > 0 and zFar > zNear, but zNear $zNear, zFar $zFar")
+            }
+            if (left == right || top == bottom) {
+                throw Exception("Error: top,bottom and left,right must not be equal")
+            }
+
+            val zNear2 = 2.0f * zNear
+            val dx = right - left
+            val dy = top - bottom
+            val dz = zFar - zNear
+            val A = (right + left) / dx
+            val B = (top + bottom) / dy
+            val C = -1.0f * (zFar + zNear) / dz
+            val D = -2.0f * (zFar * zNear) / dz
+
+            return Matrix4.fromRows(
+                zNear2 / dx, 0f, A, 0f,
+                0f, zNear2 / dy, B, 0f,
+                0f, 0f, C, D,
+                0f, 0f, -1f, 0f
+            )
+        }
+        fun frustum(left: Double, right: Double, bottom: Double, top: Double, zNear: Double = 0.0, zFar: Double = 1.0): Matrix4
+            = frustum(left.toFloat(), right.toFloat(), bottom.toFloat(), top.toFloat(), zNear.toFloat(), zFar.toFloat())
+        fun frustum(left: Int, right: Int, bottom: Int, top: Int, zNear: Int = 0, zFar: Int = 1): Matrix4
+            = frustum(left.toFloat(), right.toFloat(), bottom.toFloat(), top.toFloat(), zNear.toFloat(), zFar.toFloat())
+
+        fun perspective(fovy: Angle, aspect: Float, zNear: Float, zFar: Float): Matrix4 {
+            val top = tan(fovy.radians.toFloat() / 2f) * zNear
+            val bottom = -1.0f * top
+            val left = aspect * bottom
+            val right = aspect * top
+            return frustum(left, right, bottom, top, zNear, zFar)
+        }
+        fun perspective(fovy: Angle, aspect: Double, zNear: Double, zFar: Double): Matrix4
+            = perspective(fovy, aspect.toFloat(), zNear.toFloat(), zFar.toFloat())
+
+        fun lookAt(
+            eye: Vector3F,
+            target: Vector3F,
+            up: Vector3F
+        ): Matrix4 {
+            var z = eye - target
+            if (z.lengthSquared == 0f) z = z.copy(z = 1f)
+            z = z.normalized()
+            var x = Vector3F.cross(up, z)
+            if (x.lengthSquared == 0f) {
+                z = when {
+                    abs(up.z) == 1f -> z.copy(x = z.x + 0.0001f)
+                    else -> z.copy(z = z.z + 0.0001f)
+                }
+                z = z.normalized()
+                x = Vector3F.cross(up, z)
+            }
+            x = x.normalized()
+            val y = Vector3F.cross(z, x)
+            return Matrix4.fromRows(
+                x.x, y.x, z.x, 0f,
+                x.y, y.y, z.y, 0f,
+                x.z, y.z, z.z, 0f,
+                //-x.dot(eye), -y.dot(eye), -z.dot(eye), 1f // @TODO: Check why is this making other tests to fail
+                0f, 0f, 0f, 1f
+            )
+        }
+    }
+}
+
+data class TRS4(val translation: Vector4F, val rotation: Quaternion, val scale: Vector4F)
+
+fun Matrix4.toMatrix3(): Matrix3 = Matrix3.fromRows(
+    v00, v01, v02,
+    v10, v11, v12,
+    v20, v21, v22
+)

math/src/main/java/com/icegps/math/geometry/Matrix4Ext.kt

@@ -0,0 +1,9 @@
+package com.icegps.math.geometry
+
+fun Matrix4.Companion.ortho(rect: Rectangle, near: Float = 0f, far: Float = 1f): Matrix4 = Matrix4.ortho(rect.left, rect.right, rect.bottom, rect.top, near.toDouble(), far.toDouble())
+fun Matrix4.Companion.ortho(rect: Rectangle, near: Double = 0.0, far: Double = 1.0): Matrix4 = ortho(rect, near.toFloat(), far.toFloat())
+fun Matrix4.Companion.ortho(rect: Rectangle, near: Int = 0, far: Int = 1): Matrix4 = ortho(rect, near.toFloat(), far.toFloat())
+
+fun Matrix4.Companion.frustum(rect: Rectangle, zNear: Float = 0f, zFar: Float = 1f): Matrix4 = Matrix4.frustum(rect.left, rect.right, rect.bottom, rect.top, zNear.toDouble(), zFar.toDouble())
+fun Matrix4.Companion.frustum(rect: Rectangle, zNear: Double = 0.0, zFar: Double = 1.0): Matrix4 = frustum(rect, zNear.toFloat(), zFar.toFloat())
+fun Matrix4.Companion.frustum(rect: Rectangle, zNear: Int = 0, zFar: Int = 1): Matrix4 = frustum(rect, zNear.toFloat(), zFar.toFloat())