* Remove the 'DoG' filter.
[master-thesis.git] / Parasitemia / Parasitemia / ImgTools.fs
index e9cbbe9..24c6998 100644 (file)
@@ -26,7 +26,13 @@ let gaussianFilter (img : Image<'TColor, 'TDepth>) (standardDeviation : float) :
     img.SmoothGaussian(size, size, standardDeviation, standardDeviation)
 
 
-let findMaxima (img: Image<Gray, byte>) : IEnumerable<HashSet<Point>> =
+type Points = HashSet<Point>
+
+let drawPoints (img: Image<Gray, byte>) (points: Points) (intensity: byte) =
+    for p in points do
+        img.Data.[p.Y, p.X, 0] <- intensity
+
+let findMaxima (img: Image<Gray, byte>) : IEnumerable<Points> =
     use suppress = new Image<Gray, byte>(img.Size)
     let w = img.Width
     let h = img.Height
@@ -90,159 +96,189 @@ let findMaxima (img: Image<Gray, byte>) : IEnumerable<HashSet<Point>> =
             if maxima.Count > 0
             then result.AddRange(maxima)
 
-    result.Select(fun l -> HashSet<Point>(l))
+    result.Select(fun l -> Points(l))
 
 
 type PriorityQueue () =
-    let q = List<HashSet<Point>>() // TODO: Check performance with an HasSet
+    let size = 256
+    let q: Points[] = Array.init size (fun i -> Points()) // TODO: Check performance with an HasSet
     let mutable highest = -1 // Value of the first elements of 'q'.
+    let mutable lowest = size
 
     member this.Next () : byte * Point =
         if this.IsEmpty
         then
             invalidOp "Queue is empty"
         else
-            let l = q.[0]
+            let l = q.[highest]
             let next = l.First()
             l.Remove(next) |> ignore
             let value = byte highest
+
             if l.Count = 0
             then
-                q.RemoveAt(0)
                 highest <- highest - 1
-                while q.Count > 0 && q.[0] = null do
-                    q.RemoveAt(0)
+                while highest > lowest && q.[highest].Count = 0 do
                     highest <- highest - 1
+                if highest = lowest
+                then
+                    highest <- -1
+                    lowest <- size
+
             value, next
 
     member this.Max =
         highest |> byte
 
-    (*member this.UnionWith (other: PriorityQueue) =
-        while not other.IsEmpty do
-            let p, v = other.Next
-            this.Add p v*)
-
     member this.Add (value: byte) (p: Point) =
         let vi = int value
 
-        if this.IsEmpty
-        then
-            highest <- int value
-            q.Insert(0, null)
-        elif vi > highest
+        if vi > highest
         then
-            for i in highest .. vi - 1  do
-                q.Insert(0, null)
             highest <- vi
-        elif highest - vi >= q.Count
+        if vi <= lowest
         then
-            for i in 0 .. highest - vi - q.Count do
-                q.Add(null)
+            lowest <- vi - 1
 
-        let pos = highest - vi
-        if q.[pos] = null
+        q.[vi].Add(p) |> ignore
+
+    member this.Remove (value: byte) (p: Point) =
+        let vi = int value
+        if q.[vi].Remove(p) && q.[vi].Count = 0
         then
-            q.[pos] <- HashSet<Point>([p])
-        else
-            q.[pos].Add(p) |> ignore
+            if vi = highest
+            then
+                highest <- highest - 1
+                while highest > lowest && q.[highest].Count = 0 do
+                    highest <- highest - 1
+            elif vi - 1 = lowest
+            then
+                lowest <- lowest + 1
+                while lowest < highest && q.[lowest + 1].Count = 0 do
+                    lowest <- lowest + 1
 
+            if highest = lowest // The queue is now empty.
+            then
+                highest <- -1
+                lowest <- size
 
     member this.IsEmpty =
-        q.Count = 0
+        highest = -1
 
     member this.Clear () =
-        while highest >= 0 do
-            q.[highest] <- null
+        while highest > lowest  do
+            q.[highest].Clear()
             highest <- highest - 1
+        highest <- -1
+        lowest <- size
 
 
+type AreaState =
+    | Removed = 1
+    | Unprocessed = 2
+    | Validated = 3
 
-type MaximaState =  Uncertain | Validated | TooBig
-type Maxima = {
-    elements : HashSet<Point>
-    mutable intensity: byte option
-    mutable state: MaximaState }
-
+[<AllowNullLiteral>]
+type Area (elements: Points) =
+    member this.Elements = elements
+    member val Intensity = None with get, set
+    member val State = AreaState.Unprocessed with get, set
 
 let areaOpen (img: Image<Gray, byte>) (area: int) =
     let w = img.Width
     let h = img.Height
+    let imgData = img.Data
+
+    let areas = List<Area>(findMaxima img |> Seq.map Area)
+
+    let pixels: Area[,] = Array2D.create h w null
+    for m in areas do
+        for e in m.Elements do
+            pixels.[e.Y, e.X] <- m
 
-    let maxima = findMaxima img |> Seq.map (fun m -> { elements = m; intensity = None; state = Uncertain }) |> List.ofSeq
-    let toValidated = Stack<Maxima>(maxima)
+    let queue = PriorityQueue()
 
-    while toValidated.Count > 0 do
-        let m = toValidated.Pop()
-        if m.elements.Count <= area
+    let addEdgeToQueue (elements: Points) =
+        for p in elements do
+            for i in -1 .. 1 do
+                for j in -1 .. 1 do
+                    if i <> 0 || j <> 0
+                    then
+                        let ni = i + p.Y
+                        let nj = j + p.X
+                        let p' = Point(nj, ni)
+                        if ni >= 0 && ni < h && nj >= 0 && nj < w && not (elements.Contains(p'))
+                        then
+                            queue.Add (imgData.[ni, nj, 0]) p'
+
+    // Reverse order is quicker.
+    for i in areas.Count - 1 .. -1 .. 0 do
+        let m = areas.[i]
+        if m.Elements.Count <= area && m.State <> AreaState.Removed
         then
-            let queue =
-                let q = PriorityQueue()
-                let firstElements = HashSet<Point>()
-                for p in m.elements do
-                    for i in -1 .. 1 do
-                        for j in -1 .. 1 do
-                            if i <> 0 || j <> 0
-                            then
-                                let ni = i + p.Y
-                                let nj = j + p.X
-                                let p' = Point(nj, ni)
-                                if ni >= 0 && ni < h && nj >= 0 && nj < w && not (m.elements.Contains(p')) && not (firstElements.Contains(p'))
-                                then
-                                    firstElements.Add(p') |> ignore
-                                    q.Add (img.Data.[ni, nj, 0]) p'
-                q
+            queue.Clear()
+            addEdgeToQueue m.Elements
 
             let mutable intensity = queue.Max
-            let nextElements = HashSet<Point>()
+            let nextElements = Points()
 
             let mutable stop = false
             while not stop do
                 let intensity', p = queue.Next ()
+                let mutable merged = false
 
                 if intensity' = intensity // The intensity doesn't change.
                 then
-                    if m.elements.Count + nextElements.Count + 1 > area
+                    if m.Elements.Count + nextElements.Count + 1 > area
                     then
-                        m.state <- Validated
-                        m.intensity <- Some intensity
+                        m.State <- AreaState.Validated
+                        m.Intensity <- Some intensity
                         stop <- true
                     else
                         nextElements.Add(p) |> ignore
+
                 elif intensity' < intensity
                 then
-                    m.elements.UnionWith(nextElements)
-                    if m.elements.Count = area
+                    m.Elements.UnionWith(nextElements)
+                    for e in nextElements do
+                        pixels.[e.Y, e.X] <- m
+
+                    if m.Elements.Count = area
                     then
-                        m.state <- Validated
-                        m.intensity <- Some (intensity')
+                        m.State <- AreaState.Validated
+                        m.Intensity <- Some (intensity')
                         stop <- true
                     else
                         intensity <- intensity'
                         nextElements.Clear()
                         nextElements.Add(p) |> ignore
+
                 else // i' > i
-                    seq {
-                        for m' in maxima do
-                            if m' <> m && m'.elements.Contains(p) then
-                                if m'.elements.Count + m.elements.Count <= area
-                                then
-                                    m'.state <- Uncertain
-                                    m'.elements.UnionWith(m.elements)
-                                    if not <| toValidated.Contains m' // FIXME: Maybe use state instead of scanning the whole list.
-                                    then
-                                        toValidated.Push(m')
-                                    stop <- true
-                                yield false
-                    } |> Seq.forall id |> ignore
+                    let m' = pixels.[p.Y, p.X]
+                    if m' <> null
+                    then
+                        if m'.Elements.Count + m.Elements.Count <= area
+                        then
+                            m'.State <- AreaState.Removed
+                            for e in m'.Elements do
+                                pixels.[e.Y, e.X] <- m
+                                queue.Remove imgData.[e.Y, e.X, 0] e
+                            addEdgeToQueue m'.Elements
+                            m.Elements.UnionWith(m'.Elements)
+                            let intensityMax = queue.Max
+                            if intensityMax <> intensity
+                            then
+                                intensity <- intensityMax
+                                nextElements.Clear()
+                            merged <- true
 
-                    if not stop
+                    if not merged
                     then
-                        m.state <- Validated
-                        m.intensity <- Some (intensity)
+                        m.State <- AreaState.Validated
+                        m.Intensity <- Some (intensity)
                         stop <- true
 
-                if not stop
+                if not stop && not merged
                 then
                     for i in -1 .. 1 do
                         for j in -1 .. 1 do
@@ -253,29 +289,29 @@ let areaOpen (img: Image<Gray, byte>) (area: int) =
                                 let p' = Point(nj, ni)
                                 if ni < 0 || ni >= h || nj < 0 || nj >= w
                                 then
-                                    m.state <- Validated
-                                    m.intensity <- Some (intensity)
+                                    m.State <- AreaState.Validated
+                                    m.Intensity <- Some (intensity)
                                     stop <- true
-                                elif not (m.elements.Contains(p')) && not (nextElements.Contains(p'))
+                                elif not (m.Elements.Contains(p')) && not (nextElements.Contains(p'))
                                 then
-                                    queue.Add (img.Data.[ni, nj, 0]) p'
+                                    queue.Add (imgData.[ni, nj, 0]) p'
 
                 if queue.IsEmpty
                 then
-                    if m.elements.Count + nextElements.Count <= area
+                    if m.Elements.Count + nextElements.Count <= area
                     then
-                        m.state <- Validated
-                        m.intensity <- Some intensity'
-                        m.elements.UnionWith(nextElements)
+                        m.State <- AreaState.Validated
+                        m.Intensity <- Some intensity'
+                        m.Elements.UnionWith(nextElements)
                     stop <- true
 
-    for m in maxima do
-        if m.state = Validated
+    for m in areas do
+        if m.State = AreaState.Validated
         then
-            match m.intensity with
+            match m.Intensity with
             | Some i ->
-                for p in m.elements do
-                    img.Data.[p.Y, p.X, 0] <- i
+                for p in m.Elements do
+                    imgData.[p.Y, p.X, 0] <- i
             | _ -> ()
     ()
 
@@ -386,7 +422,7 @@ let connectedComponents (img: Image<Gray, byte>) (startPoints: List<Point>) : Li
     let w = img.Width
     let h = img.Height
 
-    let pointChecked = HashSet<Point>()
+    let pointChecked = Points()
     let pointToCheck = List<Point>(startPoints);
 
     let data = img.Data