* Remove the 'DoG' filter.
[master-thesis.git] / Parasitemia / Parasitemia / ImgTools.fs
index 7253390..24c6998 100644 (file)
@@ -3,11 +3,13 @@
 open System
 open System.Drawing
 open System.Collections.Generic
+open System.Linq
 
 open Emgu.CV
 open Emgu.CV.Structure
 
 open Utils
+open Heap
 
 // Normalize image values between 0uy and 255uy.
 let normalizeAndConvert (img: Image<Gray, float32>) : Image<Gray, byte> =
@@ -18,10 +20,302 @@ let normalizeAndConvert (img: Image<Gray, float32>) : Image<Gray, byte> =
     img.MinMax(min, max, minLocation, maxLocation)
     ((img - (!min).[0]) / ((!max).[0] - (!min).[0]) * 255.0).Convert<Gray, byte>()
 
+
 let gaussianFilter (img : Image<'TColor, 'TDepth>) (standardDeviation : float) : Image<'TColor, 'TDepth> =
     let size = 2 * int (ceil (4.0 * standardDeviation)) + 1
     img.SmoothGaussian(size, size, standardDeviation, standardDeviation)
 
+
+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
+
+    let imgData = img.Data
+    let suppressData = suppress.Data
+
+    let result = List<List<Point>>()
+
+    let flood (start: Point) : List<List<Point>> =
+        let sameLevelToCheck = Stack<Point>()
+        let betterLevelToCheck = Stack<Point>()
+        betterLevelToCheck.Push(start)
+
+        let result' = List<List<Point>>()
+
+        while betterLevelToCheck.Count > 0 do
+            let p = betterLevelToCheck.Pop()
+            if suppressData.[p.Y, p.X, 0] = 0uy
+            then
+                suppressData.[p.Y, p.X, 0] <- 1uy
+                sameLevelToCheck.Push(p)
+                let current = List<Point>()
+
+                let mutable betterExists = false
+
+                while sameLevelToCheck.Count > 0 do
+                    let p' = sameLevelToCheck.Pop()
+                    let currentLevel = imgData.[p'.Y, p'.X, 0]
+                    current.Add(p') |> ignore
+                    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
+                                if ni >= 0 && ni < h && nj >= 0 && nj < w
+                                then
+                                    let level = imgData.[ni, nj, 0]
+                                    let notSuppressed = suppressData.[ni, nj, 0] = 0uy
+
+                                    if level = currentLevel && notSuppressed
+                                    then
+                                        suppressData.[ni, nj, 0] <- 1uy
+                                        sameLevelToCheck.Push(Point(nj, ni))
+                                    elif level > currentLevel
+                                    then
+                                        betterExists <- true
+                                        if notSuppressed
+                                        then
+                                            betterLevelToCheck.Push(Point(nj, ni))
+
+                if not betterExists
+                then
+                    result'.Add(current)
+        result'
+
+    for i in 0 .. h - 1 do
+        for j in 0 .. w - 1 do
+            let maxima = flood (Point(j, i))
+            if maxima.Count > 0
+            then result.AddRange(maxima)
+
+    result.Select(fun l -> Points(l))
+
+
+type PriorityQueue () =
+    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.[highest]
+            let next = l.First()
+            l.Remove(next) |> ignore
+            let value = byte highest
+
+            if l.Count = 0
+            then
+                highest <- highest - 1
+                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.Add (value: byte) (p: Point) =
+        let vi = int value
+
+        if vi > highest
+        then
+            highest <- vi
+        if vi <= lowest
+        then
+            lowest <- vi - 1
+
+        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
+            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 =
+        highest = -1
+
+    member this.Clear () =
+        while highest > lowest  do
+            q.[highest].Clear()
+            highest <- highest - 1
+        highest <- -1
+        lowest <- size
+
+
+type AreaState =
+    | Removed = 1
+    | Unprocessed = 2
+    | Validated = 3
+
+[<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 queue = PriorityQueue()
+
+    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
+            queue.Clear()
+            addEdgeToQueue m.Elements
+
+            let mutable intensity = queue.Max
+            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
+                    then
+                        m.State <- AreaState.Validated
+                        m.Intensity <- Some intensity
+                        stop <- true
+                    else
+                        nextElements.Add(p) |> ignore
+
+                elif intensity' < intensity
+                then
+                    m.Elements.UnionWith(nextElements)
+                    for e in nextElements do
+                        pixels.[e.Y, e.X] <- m
+
+                    if m.Elements.Count = area
+                    then
+                        m.State <- AreaState.Validated
+                        m.Intensity <- Some (intensity')
+                        stop <- true
+                    else
+                        intensity <- intensity'
+                        nextElements.Clear()
+                        nextElements.Add(p) |> ignore
+
+                else // i' > i
+                    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 merged
+                    then
+                        m.State <- AreaState.Validated
+                        m.Intensity <- Some (intensity)
+                        stop <- true
+
+                if not stop && not merged
+                then
+                    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
+                                then
+                                    m.State <- AreaState.Validated
+                                    m.Intensity <- Some (intensity)
+                                    stop <- true
+                                elif not (m.Elements.Contains(p')) && not (nextElements.Contains(p'))
+                                then
+                                    queue.Add (imgData.[ni, nj, 0]) p'
+
+                if queue.IsEmpty
+                then
+                    if m.Elements.Count + nextElements.Count <= area
+                    then
+                        m.State <- AreaState.Validated
+                        m.Intensity <- Some intensity'
+                        m.Elements.UnionWith(nextElements)
+                    stop <- true
+
+    for m in areas do
+        if m.State = AreaState.Validated
+        then
+            match m.Intensity with
+            | Some i ->
+                for p in m.Elements do
+                    imgData.[p.Y, p.X, 0] <- i
+            | _ -> ()
+    ()
+
+
 // Zhang and Suen algorithm.
 // Modify 'mat' in place.
 let thin (mat: Matrix<byte>) =
@@ -73,6 +367,7 @@ let thin (mat: Matrix<byte>) =
         data2 <- tmp
 
 
+// FIXME: replace by a queue or stack.
 let pop (l: List<'a>) : 'a =
     let n = l.[l.Count - 1]
     l.RemoveAt(l.Count - 1)
@@ -127,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