Finding ellipses and parasites.
[master-thesis.git] / Parasitemia / Parasitemia / KdTree.fs
diff --git a/Parasitemia/Parasitemia/KdTree.fs b/Parasitemia/Parasitemia/KdTree.fs
new file mode 100644 (file)
index 0000000..dd0247a
--- /dev/null
@@ -0,0 +1,154 @@
+module KdTree
+
+open System
+
+type I2DCoords =
+    abstract X : float
+    abstract Y : float
+
+// Compare 'e1' and 'e2' by X.
+let cmpX (e1: I2DCoords) (e2: I2DCoords) : int =
+    match e1.X.CompareTo(e2.X) with
+    | 0 -> match e1.Y.CompareTo(e2.Y) with 
+           | 0 -> e1.GetHashCode().CompareTo(e2.GetHashCode())
+           | v -> v
+    | v -> v
+    
+// Compare 'e1' and 'e2' by Y.
+let cmpY (e1: I2DCoords) (e2: I2DCoords) : int =
+    match e1.Y.CompareTo(e2.Y) with
+    | 0 -> match e1.X.CompareTo(e2.X) with 
+           | 0 -> e1.GetHashCode().CompareTo(e2.GetHashCode())
+           | v -> v
+    | v -> v
+
+type Region = { minX: float; maxX: float; minY: float; maxY: float } with
+    member this.Contains px py : bool =
+        px >= this.minX && px <= this.maxX && 
+        py >= this.minY && py <= this.maxY
+
+    member this.IsSub otherRegion : bool =
+        this.minX >= otherRegion.minX && this.maxX <= otherRegion.maxX && 
+        this.minY >= otherRegion.minY && this.maxY <= otherRegion.maxY
+
+    member this.Intersects otherRegion : bool = 
+        this.minX < otherRegion.maxX && this.maxX >= otherRegion.minX &&
+        this.minY < otherRegion.maxY && this.maxY >= otherRegion.minY
+
+type Tree<'a when 'a :> I2DCoords> = 
+    | Node of float * Tree<'a> * Tree<'a>
+    | Leaf of 'a
+
+    static member buildTree (l: 'a list) : Tree<'a> =
+        let xSorted = List.toArray l
+        let ySorted = List.toArray l
+        Array.sortInPlaceWith cmpX xSorted
+        Array.sortInPlaceWith cmpY ySorted
+
+        let rec buildTreeFromSortedArray (pXSorted: 'a[]) (pYSorted: 'a[]) (depth: int) : Tree<'a> =
+            if pXSorted.Length = 1
+            then 
+                Leaf pXSorted.[0]
+            else                
+                if depth % 2 = 1 // 'depth' is odd -> vertical splitting else horizontal splitting.
+                then                
+                    let leftX, rightX = Array.splitAt ((pXSorted.Length + 1) / 2) pXSorted
+                    let splitElement = Array.last leftX
+                    let leftY, rightY = Array.partition (fun (e: 'a) -> cmpX e splitElement <= 0) pYSorted // FIXME: Maybe this operation can be optimized.
+                    Node (splitElement.X, buildTreeFromSortedArray leftX leftY (depth + 1), buildTreeFromSortedArray rightX rightY (depth + 1))
+                else
+                    let downY, upY = Array.splitAt ((pYSorted.Length + 1) / 2) pYSorted
+                    let splitElement = Array.last downY
+                    let downX, upX = Array.partition (fun (e: 'a) -> cmpY e splitElement <= 0) pXSorted // FIXME: Maybe this operation can be optimized.
+                    Node (splitElement.Y, buildTreeFromSortedArray downX downY (depth + 1), buildTreeFromSortedArray upX upY (depth + 1))
+
+        buildTreeFromSortedArray xSorted ySorted 1
+
+    static member search (tree: Tree<'a>) (searchRegion: Region) : 'a list =
+        let rec valuesFrom (tree: Tree<'a>) : 'a list =
+            match tree with
+            | Leaf v -> [v]
+            | Node (_, part1, part2) -> (valuesFrom part1) @ (valuesFrom part2)
+
+        let rec searchWithRegion (tree: Tree<'a>) (currentRegion: Region) (depth: int) : 'a list =
+            match tree with
+            | Leaf v -> if searchRegion.Contains v.X v.Y then [v] else []
+            | Node (splitValue, part1, part2) ->
+                let valuesInRegion (region: Region) (treeRegion: Tree<'a>) =  
+                    if region.IsSub searchRegion
+                    then
+                        valuesFrom treeRegion
+                    elif region.Intersects searchRegion
+                    then
+                        searchWithRegion treeRegion region (depth + 1)
+                    else 
+                        []
+
+                if depth % 2 = 1 // Vertical splitting.
+                then                
+                    let leftRegion = { currentRegion with maxX = splitValue }  
+                    let rightRegion = { currentRegion with minX = splitValue }
+                    (valuesInRegion leftRegion part1) @ (valuesInRegion rightRegion part2)
+                else // Horizontal splitting.
+                    let downRegion = { currentRegion with maxY = splitValue }  
+                    let upRegion = { currentRegion with minY = splitValue }
+                    (valuesInRegion downRegion part1) @ (valuesInRegion upRegion part2)
+                
+        searchWithRegion tree { minX = Double.MinValue; maxX = Double.MaxValue; minY = Double.MinValue; maxY = Double.MaxValue } 1
+            
+        
+///// Tests. TODO: to put in a unit test.
+
+type Point (x: float, y: float) =
+    interface I2DCoords with
+        member this.X = x
+        member this.Y = y
+
+    override this.ToString () =
+        sprintf "(%.1f, %.1f)" x y
+
+// TODO: test with identical X or Y coords
+let test () = 
+    let pts = [
+        Point(1.0, 1.0)
+        Point(2.0, 2.0)
+        Point(1.5, 3.6)
+        Point(3.0, 3.2)
+        Point(4.0, 4.0)
+        Point(3.5, 1.5)
+        Point(2.5, 0.5) ]
+    
+    let tree = Tree.buildTree pts
+    Utils.dprintfn "Tree: %A" tree
+
+    let s1 = Tree.search tree { minX = 0.0; maxX = 5.0; minY = 0.0; maxY = 5.0 } // All points.
+    Utils.dprintfn "s1: %A" s1
+
+    let s2 = Tree.search tree { minX = 2.8; maxX = 4.5; minY = 3.0; maxY = 4.5 }
+    Utils.dprintfn "s2: %A" s2
+    
+    let s3 = Tree.search tree { minX = 2.0; maxX = 2.0; minY = 2.0; maxY = 2.0 }
+    Utils.dprintfn "s3: %A" s3
+
+let test2 () = 
+    let pts = [
+        Point(1.0, 1.0)
+        Point(1.0, 2.0)
+        Point(1.0, 3.0) ]
+        
+    let tree = Tree.buildTree pts
+    Utils.dprintfn "Tree: %A" tree
+        
+    let s1 = Tree.search tree { minX = 1.0; maxX = 1.0; minY = 1.0; maxY = 1.0 }
+    Utils.dprintfn "s1: %A" s1
+
+    let s2 = Tree.search tree { minX = 1.0; maxX = 1.0; minY = 2.0; maxY = 2.0 }
+    Utils.dprintfn "s2: %A" s2
+    
+    // This case result is wrong: FIXME
+    let s3 = Tree.search tree { minX = 1.0; maxX = 1.0; minY = 3.0; maxY = 3.0 }
+    Utils.dprintfn "s3: %A" s3
+
+    let s4 = Tree.search tree { minX = 0.0; maxX = 2.0; minY = 0.0; maxY = 4.0 }
+    Utils.dprintfn "s4: %A" s4
+