implementation of geometry projections of vectors to a line/plane

kmath-geometry/src/commonMain/kotlin/space/kscience/kmath/geometry/Euclidean2DSpace.kt

@@ -18,8 +18,8 @@ public interface Vector2D : Point<Double>, Vector{
     override val size: Int get() = 2
 
     override operator fun get(index: Int): Double = when (index) {
-        1 -> x
+        0 -> x
-        2 -> y
+        1 -> y
         else -> error("Accessing outside of point bounds")
     }
 
@@ -27,7 +27,7 @@ public interface Vector2D : Point<Double>, Vector{
 }
 
 public val Vector2D.r: Double
-    get() = Euclidean2DSpace { sqrt(norm()) }
+    get() = Euclidean2DSpace { norm() }
 
 @Suppress("FunctionName")
 public fun Vector2D(x: Double, y: Double): Vector2D = Vector2DImpl(x, y)

kmath-geometry/src/commonMain/kotlin/space/kscience/kmath/geometry/Euclidean3DSpace.kt

@@ -19,9 +19,9 @@ public interface Vector3D : Point<Double>, Vector {
     override val size: Int get() = 3
 
     override operator fun get(index: Int): Double = when (index) {
-        1 -> x
+        0 -> x
-        2 -> y
+        1 -> y
-        3 -> z
+        2 -> z
         else -> error("Accessing outside of point bounds")
     }
 
@@ -31,7 +31,7 @@ public interface Vector3D : Point<Double>, Vector {
 @Suppress("FunctionName")
 public fun Vector3D(x: Double, y: Double, z: Double): Vector3D = Vector3DImpl(x, y, z)
 
-public val Vector3D.r: Double get() = Euclidean3DSpace { sqrt(norm()) }
+public val Vector3D.r: Double get() = Euclidean3DSpace { norm() }
 
 private data class Vector3DImpl(
     override val x: Double,

kmath-geometry/src/commonMain/kotlin/space/kscience/kmath/geometry/Projections.kt

@@ -0,0 +1,20 @@
+package space.kscience.kmath.geometry
+
+/**
+ * Project vector to a line
+ * @param vector to project
+ * @param line line to which vector should be projected
+ */
+public fun <V : Vector> GeometrySpace<V>.projectToLine(vector: V, line: Line<V>): V = with(line) {
+    base + (direction dot (vector - base)) / (direction dot direction) * direction
+}
+
+/**
+ * Project vector to a hyper-plane, which is defined by a normal and base
+ * In 2d case it is projection to a line, in 3d case it is projection to a plane
+ * @param vector to project
+ * @param normal normal (perpendicular) vector to a hyper-plane to which vector should be projected
+ * @param base point belonging to a hyper-plane to which vector should be projected
+ */
+public fun <V : Vector> GeometrySpace<V>.projectAlong(vector: V, normal: V, base: V): V =
+    vector + normal * ((base - vector) dot normal) / (normal dot normal)

kmath-geometry/src/commonTest/kotlin/space/kscience/kmath/geometry/Euclidean2DSpaceTest.kt

@@ -0,0 +1,62 @@
+package space.kscience.kmath.geometry
+
+import kotlin.math.sqrt
+import kotlin.test.Test
+import kotlin.test.assertEquals
+
+internal class Euclidean2DSpaceTest {
+    @Test
+    fun getZero() {
+        assertVectorEquals(Vector2D(0.0, 0.0), Euclidean2DSpace.zero)
+    }
+
+    @Test
+    fun norm() {
+        with(Euclidean2DSpace) {
+            assertEquals(0.0, zero.norm())
+            assertEquals(1.0, Vector2D(1.0, 0.0).norm())
+            assertEquals(sqrt(2.0), Vector2D(1.0, 1.0).norm())
+            assertEquals(sqrt(5.002001), Vector2D(-2.0, 1.001).norm())
+        }
+    }
+
+    @Test
+    fun dotProduct() {
+        with(Euclidean2DSpace) {
+            assertEquals(0.0, zero dot zero)
+            assertEquals(0.0, zero dot Vector2D(1.0, 0.0))
+            assertEquals(0.0, Vector2D(-2.0, 0.001) dot zero)
+            assertEquals(0.0, Vector2D(1.0, 0.0) dot Vector2D(0.0, 1.0))
+
+            assertEquals(1.0, Vector2D(1.0, 0.0) dot Vector2D(1.0, 0.0))
+            assertEquals(-2.0, Vector2D(0.0, 1.0) dot Vector2D(1.0, -2.0))
+            assertEquals(2.0, Vector2D(1.0, 1.0) dot Vector2D(1.0, 1.0))
+            assertEquals(4.001001, Vector2D(-2.0, 1.001) dot Vector2D(-2.0, 0.001))
+
+            assertEquals(-4.998, Vector2D(1.0, 2.0) dot Vector2D(-5.0, 0.001))
+        }
+    }
+
+    @Test
+    fun add() {
+        with(Euclidean2DSpace) {
+            assertVectorEquals(
+                Vector2D(-2.0, 0.001),
+                Vector2D(-2.0, 0.001) + zero
+            )
+            assertVectorEquals(
+                Vector2D(-3.0, 3.001),
+                Vector2D(2.0, 3.0) + Vector2D(-5.0, 0.001)
+            )
+        }
+    }
+
+    @Test
+    fun multiply() {
+        with(Euclidean2DSpace) {
+            assertVectorEquals(Vector2D(-4.0, 0.0), Vector2D(-2.0, 0.0) * 2)
+            assertVectorEquals(Vector2D(4.0, 0.0), Vector2D(-2.0, 0.0) * -2)
+            assertVectorEquals(Vector2D(300.0, 0.0003), Vector2D(100.0, 0.0001) * 3)
+        }
+    }
+}

kmath-geometry/src/commonTest/kotlin/space/kscience/kmath/geometry/Euclidean3DSpaceTest.kt

@@ -0,0 +1,74 @@
+package space.kscience.kmath.geometry
+
+import kotlin.test.Test
+import kotlin.test.assertEquals
+
+internal class Euclidean3DSpaceTest {
+    @Test
+    fun getZero() {
+        assertVectorEquals(Vector3D(0.0, 0.0, 0.0), Euclidean3DSpace.zero)
+    }
+
+    @Test
+    fun distance() {
+        with(Euclidean3DSpace) {
+            assertEquals(0.0, zero.distanceTo(zero))
+            assertEquals(1.0, zero.distanceTo(Vector3D(1.0, 0.0, 0.0)))
+            assertEquals(kotlin.math.sqrt(5.000001), Vector3D(1.0, -2.0, 0.001).distanceTo(zero))
+            assertEquals(0.0, Vector3D(1.0, -2.0, 0.001).distanceTo(Vector3D(1.0, -2.0, 0.001)))
+            assertEquals(0.0, Vector3D(1.0, 0.0, 0.0).distanceTo(Vector3D(1.0, 0.0, 0.0)))
+            assertEquals(kotlin.math.sqrt(2.0), Vector3D(1.0, 0.0, 0.0).distanceTo(Vector3D(1.0, 1.0, 1.0)))
+            assertEquals(3.1622778182822584, Vector3D(0.0, 1.0, 0.0).distanceTo(Vector3D(1.0, -2.0, 0.001)))
+            assertEquals(0.0, Vector3D(1.0, -2.0, 0.001).distanceTo(Vector3D(1.0, -2.0, 0.001)))
+            assertEquals(9.695050335093676, Vector3D(1.0, 2.0, 3.0).distanceTo(Vector3D(7.0, -5.0, 0.001)))
+        }
+    }
+
+    @Test
+    fun norm() {
+        with(Euclidean3DSpace) {
+            assertEquals(0.0, zero.norm())
+            assertEquals(1.0, Vector3D(1.0, 0.0, 0.0).norm())
+            assertEquals(kotlin.math.sqrt(3.0), Vector3D(1.0, 1.0, 1.0).norm())
+            assertEquals(kotlin.math.sqrt(5.000001), Vector3D(1.0, -2.0, 0.001).norm())
+        }
+    }
+
+    @Test
+    fun dotProduct() {
+        with(Euclidean3DSpace) {
+            assertEquals(0.0, zero dot zero)
+            assertEquals(0.0, zero dot Vector3D(1.0, 0.0, 0.0))
+            assertEquals(0.0, Vector3D(1.0, -2.0, 0.001) dot zero)
+
+            assertEquals(1.0, Vector3D(1.0, 0.0, 0.0) dot Vector3D(1.0, 0.0, 0.0))
+            assertEquals(1.0, Vector3D(1.0, 0.0, 0.0) dot Vector3D(1.0, 1.0, 1.0))
+            assertEquals(-2.0, Vector3D(0.0, 1.0, 0.0) dot Vector3D(1.0, -2.0, 0.001))
+            assertEquals(3.0, Vector3D(1.0, 1.0, 1.0) dot Vector3D(1.0, 1.0, 1.0))
+            assertEquals(5.000001, Vector3D(1.0, -2.0, 0.001) dot Vector3D(1.0, -2.0, 0.001))
+
+            assertEquals(-2.997, Vector3D(1.0, 2.0, 3.0) dot Vector3D(7.0, -5.0, 0.001))
+        }
+    }
+
+    @Test
+    fun add() {
+        with(Euclidean3DSpace) {
+            assertVectorEquals(
+                Vector3D(1.0, -2.0, 0.001),
+                Vector3D(1.0, -2.0, 0.001) + zero
+            )
+            assertVectorEquals(
+                Vector3D(8.0, -3.0, 3.001),
+                Vector3D(1.0, 2.0, 3.0) + Vector3D(7.0, -5.0, 0.001)
+            )
+        }
+    }
+
+    @Test
+    fun multiply() {
+        with(Euclidean3DSpace) {
+            assertVectorEquals(Vector3D(2.0, -4.0, 0.0), Vector3D(1.0, -2.0, 0.0) * 2)
+        }
+    }
+}

kmath-geometry/src/commonTest/kotlin/space/kscience/kmath/geometry/ProjectionAlongTest.kt

@@ -0,0 +1,60 @@
+package space.kscience.kmath.geometry
+
+import kotlin.test.Test
+import kotlin.test.assertTrue
+
+internal class ProjectionAlongTest {
+    @Test
+    fun projectionIntoYEqualsX() {
+        with(Euclidean2DSpace) {
+            val normal = Vector2D(-2.0, 2.0)
+            val base = Vector2D(2.3, 2.3)
+
+            assertVectorEquals(zero, projectAlong(zero, normal, base))
+
+            grid(-10.0..10.0, -10.0..10.0, 0.15).forEach { (x, y) ->
+                val d = (y - x) / 2.0
+                assertVectorEquals(Vector2D(x + d, y - d), projectAlong(Vector2D(x, y), normal, base))
+            }
+        }
+    }
+
+    @Test
+    fun projectionOntoLine() {
+        with(Euclidean2DSpace) {
+            val a = 5.0
+            val b = -3.0
+            val c = -15.0
+            val normal = Vector2D(-5.0, 3.0)
+            val base = Vector2D(3.0, 0.0)
+
+            grid(-10.0..10.0, -10.0..10.0, 0.15).forEach { (x, y) ->
+                val xProj = (b * (b * x - a * y) - a * c) / (a * a + b * b)
+                val yProj = (a * (-b * x + a * y) - b * c) / (a * a + b * b)
+                assertVectorEquals(Vector2D(xProj, yProj), projectAlong(Vector2D(x, y), normal, base))
+            }
+        }
+    }
+
+    @Test
+    fun projectOntoPlane() {
+        val normal = Vector3D(1.0, 3.5, 0.07)
+        val base = Vector3D(2.0, -0.0037, 11.1111)
+
+        with(Euclidean3DSpace) {
+            for (x in (-10.0..10.0).generateList(0.15)) {
+                for (y in (-10.0..10.0).generateList(0.15)) {
+                    for (z in (-10.0..10.0).generateList(0.15)) {
+                        val v = Vector3D(x, y, z)
+                        val result = projectAlong(v, normal, base)
+
+                        // assert that result is on plane
+                        assertTrue(isOrthogonal(result - base, normal))
+                        // assert that PV vector is collinear to normal vector
+                        assertTrue(isCollinear(v - result, normal))
+                    }
+                }
+            }
+        }
+    }
+}

kmath-geometry/src/commonTest/kotlin/space/kscience/kmath/geometry/ProjectionOntoLineTest.kt

@@ -0,0 +1,82 @@
+package space.kscience.kmath.geometry
+
+import kotlin.test.Test
+import kotlin.test.assertTrue
+
+internal class ProjectionOntoLineTest {
+    @Test
+    fun projectionIntoOx() {
+        with(Euclidean2DSpace) {
+            val ox = Line(zero, Vector2D(1.0, 0.0))
+
+            grid(-10.0..10.0, -10.0..10.0, 0.15).forEach { (x, y) ->
+                assertVectorEquals(Vector2D(x, 0.0), projectToLine(Vector2D(x, y), ox))
+            }
+        }
+    }
+
+    @Test
+    fun projectionIntoOy() {
+        with(Euclidean2DSpace) {
+            val line = Line(zero, Vector2D(0.0, 1.0))
+
+            grid(-10.0..10.0, -10.0..10.0, 0.15).forEach { (x, y) ->
+                assertVectorEquals(Vector2D(0.0, y), projectToLine(Vector2D(x, y), line))
+            }
+        }
+    }
+
+    @Test
+    fun projectionIntoYEqualsX() {
+        with(Euclidean2DSpace) {
+            val line = Line(zero, Vector2D(1.0, 1.0))
+
+            assertVectorEquals(zero, projectToLine(zero, line))
+
+            grid(-10.0..10.0, -10.0..10.0, 0.15).forEach { (x, y) ->
+                val d = (y - x) / 2.0
+                assertVectorEquals(Vector2D(x + d, y - d), projectToLine(Vector2D(x, y), line))
+            }
+        }
+    }
+
+    @Test
+    fun projectionOntoLine2d() {
+        with(Euclidean2DSpace) {
+            val a = 5.0
+            val b = -3.0
+            val c = -15.0
+            val line = Line(Vector2D(3.0, 0.0), Vector2D(3.0, 5.0))
+
+            grid(-10.0..10.0, -10.0..10.0, 0.15).forEach { (x, y) ->
+                val xProj = (b * (b * x - a * y) - a * c) / (a * a + b * b)
+                val yProj = (a * (-b * x + a * y) - b * c) / (a * a + b * b)
+                assertVectorEquals(Vector2D(xProj, yProj), projectToLine(Vector2D(x, y), line))
+            }
+        }
+    }
+
+    @Test
+    fun projectionOntoLine3d() {
+        val line = Line3D(
+            base = Vector3D(1.0, 3.5, 0.07),
+            direction = Vector3D(2.0, -0.0037, 11.1111)
+        )
+
+        with(Euclidean3DSpace) {
+            for (x in (-10.0..10.0).generateList(0.15)) {
+                for (y in (-10.0..10.0).generateList(0.15)) {
+                    for (z in (-10.0..10.0).generateList(0.15)) {
+                        val v = Vector3D(x, y, z)
+                        val result = projectToLine(v, line)
+
+                        // assert that result is on line
+                        assertTrue(isCollinear(result - line.base, line.direction))
+                        // assert that PV vector is orthogonal to direction vector
+                        assertTrue(isOrthogonal(v - result, line.direction))
+                    }
+                }
+            }
+        }
+    }
+}

kmath-geometry/src/commonTest/kotlin/space/kscience/kmath/geometry/Vector2DTest.kt

@@ -0,0 +1,36 @@
+package space.kscience.kmath.geometry
+
+import space.kscience.kmath.structures.asSequence
+import space.kscience.kmath.structures.toList
+import kotlin.test.assertEquals
+import kotlin.test.Test
+
+internal class Vector2DTest {
+    private val vector = Vector2D(1.0, -7.999)
+
+    @Test
+    fun getSize() {
+        assertEquals(2, vector.size)
+    }
+
+    @Test
+    fun get() {
+        assertEquals(1.0, vector[0])
+        assertEquals(-7.999, vector[1])
+    }
+
+    @Test
+    operator fun iterator() {
+        assertEquals(listOf(1.0, -7.999), vector.toList())
+    }
+
+    @Test
+    fun getX() {
+        assertEquals(1.0, vector.x)
+    }
+
+    @Test
+    fun getY() {
+        assertEquals(-7.999, vector.y)
+    }
+}

kmath-geometry/src/commonTest/kotlin/space/kscience/kmath/geometry/Vector3DTest.kt

@@ -0,0 +1,41 @@
+package space.kscience.kmath.geometry
+
+import space.kscience.kmath.structures.toList
+import kotlin.test.Test
+import kotlin.test.assertEquals
+
+internal class Vector3DTest {
+    private val vector = Vector3D(1.0, -7.999, 0.001)
+
+    @Test
+    fun getSize() {
+        assertEquals(3, vector.size)
+    }
+
+    @Test
+    fun get() {
+        assertEquals(1.0, vector[0])
+        assertEquals(-7.999, vector[1])
+        assertEquals(0.001, vector[2])
+    }
+
+    @Test
+    operator fun iterator() {
+        assertEquals(listOf(1.0, -7.999, 0.001), vector.toList())
+    }
+
+    @Test
+    fun getX() {
+        assertEquals(1.0, vector.x)
+    }
+
+    @Test
+    fun getY() {
+        assertEquals(-7.999, vector.y)
+    }
+
+    @Test
+    fun getZ() {
+        assertEquals(0.001, vector.z)
+    }
+}

kmath-geometry/src/commonTest/kotlin/space/kscience/kmath/geometry/testUtils.kt

@@ -0,0 +1,41 @@
+package space.kscience.kmath.geometry
+
+import kotlin.math.abs
+import kotlin.test.assertEquals
+
+fun ClosedRange<Double>.generateList(step: Double): List<Double> = generateSequence(start) { previous ->
+    if (previous == Double.POSITIVE_INFINITY) return@generateSequence null
+    val next = previous + step
+    if (next > endInclusive) null else next
+}.toList()
+
+fun grid(
+    xRange: ClosedRange<Double>,
+    yRange: ClosedRange<Double>,
+    step: Double
+): List<Pair<Double, Double>> {
+    val xs = xRange.generateList(step)
+    val ys = yRange.generateList(step)
+
+    return xs.flatMap { x -> ys.map { y -> x to y } }
+}
+
+fun assertVectorEquals(expected: Vector2D, actual: Vector2D, eps: Double = 1e-6) {
+    assertEquals(expected.x, actual.x, eps)
+    assertEquals(expected.y, actual.y, eps)
+}
+
+fun assertVectorEquals(expected: Vector3D, actual: Vector3D, eps: Double = 1e-6) {
+    assertEquals(expected.x, actual.x, eps)
+    assertEquals(expected.y, actual.y, eps)
+    assertEquals(expected.z, actual.z, eps)
+}
+
+fun <V : Vector> GeometrySpace<V>.isCollinear(a: V, b: V, eps: Double = 1e-6): Boolean {
+    val aDist = a.distanceTo(zero)
+    val bDist = b.distanceTo(zero)
+    return abs(aDist) < eps || abs(bDist) < eps || abs(abs((a dot b) / (aDist * bDist)) - 1) < eps
+}
+
+fun <V : Vector> GeometrySpace<V>.isOrthogonal(a: V, b: V, eps: Double = 1e-6): Boolean =
+    abs(a dot b) < eps