Change the parasite detection method.
[master-thesis.git] / Parasitemia / Parasitemia / ImgTools.fs
index 24c6998..2f021c8 100644 (file)
@@ -32,13 +32,18 @@ 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)
+
+type ExtremumType =
+    | Maxima = 1
+    | Minima = 2
+
+let findExtremum (img: Image<Gray, byte>) (extremumType: ExtremumType) : IEnumerable<Points> =
     let w = img.Width
     let h = img.Height
+    let se = [| -1, 0; 0, -1; 1, 0; 0, 1 |]
 
     let imgData = img.Data
-    let suppressData = suppress.Data
+    let suppress: bool[,] = Array2D.zeroCreate h w
 
     let result = List<List<Point>>()
 
@@ -51,9 +56,9 @@ let findMaxima (img: Image<Gray, byte>) : IEnumerable<Points> =
 
         while betterLevelToCheck.Count > 0 do
             let p = betterLevelToCheck.Pop()
-            if suppressData.[p.Y, p.X, 0] = 0uy
+            if not suppress.[p.Y, p.X]
             then
-                suppressData.[p.Y, p.X, 0] <- 1uy
+                suppress.[p.Y, p.X] <- true
                 sameLevelToCheck.Push(p)
                 let current = List<Point>()
 
@@ -63,27 +68,24 @@ let findMaxima (img: Image<Gray, byte>) : IEnumerable<Points> =
                     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
+                    for i, j in se do
+                        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 = not suppress.[ni, nj]
+
+                            if level = currentLevel && notSuppressed
                             then
-                                let ni = i + p'.Y
-                                let nj = j + p'.X
-                                if ni >= 0 && ni < h && nj >= 0 && nj < w
+                                suppress.[ni, nj] <- true
+                                sameLevelToCheck.Push(Point(nj, ni))
+                            elif if extremumType = ExtremumType.Maxima then level > currentLevel else level < currentLevel
+                            then
+                                betterExists <- true
+                                if notSuppressed
                                 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))
+                                    betterLevelToCheck.Push(Point(nj, ni))
 
                 if not betterExists
                 then
@@ -94,18 +96,26 @@ let findMaxima (img: Image<Gray, byte>) : IEnumerable<Points> =
         for j in 0 .. w - 1 do
             let maxima = flood (Point(j, i))
             if maxima.Count > 0
-            then result.AddRange(maxima)
+            then
+                result.AddRange(maxima)
 
     result.Select(fun l -> Points(l))
 
 
+let findMaxima (img: Image<Gray, byte>) : IEnumerable<Points> =
+    findExtremum img ExtremumType.Maxima
+
+let findMinima (img: Image<Gray, byte>) : IEnumerable<Points> =
+    findExtremum img ExtremumType.Minima
+
+
 type PriorityQueue () =
     let size = 256
-    let q: Points[] = Array.init size (fun i -> Points()) // TODO: Check performance with an HasSet
+    let q: Points[] = Array.init size (fun i -> Points())
     let mutable highest = -1 // Value of the first elements of 'q'.
     let mutable lowest = size
 
-    member this.Next () : byte * Point =
+    member this.NextMax () : byte * Point =
         if this.IsEmpty
         then
             invalidOp "Queue is empty"
@@ -127,9 +137,34 @@ type PriorityQueue () =
 
             value, next
 
+    member this.NextMin () : byte * Point =
+        if this.IsEmpty
+        then
+            invalidOp "Queue is empty"
+        else
+            let l = q.[lowest + 1]
+            let next = l.First()
+            l.Remove(next) |> ignore
+            let value = byte (lowest + 1)
+
+            if l.Count = 0
+            then
+                lowest <- lowest + 1
+                while lowest < highest && q.[lowest + 1].Count = 0 do
+                    lowest <- lowest + 1
+                if highest = lowest
+                then
+                    highest <- -1
+                    lowest <- size
+
+            value, next
+
     member this.Max =
         highest |> byte
 
+    member this.Min =
+        lowest + 1 |> byte
+
     member this.Add (value: byte) (p: Point) =
         let vi = int value
 
@@ -173,23 +208,28 @@ type PriorityQueue () =
         lowest <- size
 
 
-type AreaState =
+type private AreaState =
     | Removed = 1
     | Unprocessed = 2
     | Validated = 3
 
+type private AreaOperation =
+    | Opening = 1
+    | Closing = 2
+
 [<AllowNullLiteral>]
-type Area (elements: Points) =
+type private 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 private areaOperation (img: Image<Gray, byte>) (area: int) (op: AreaOperation) =
     let w = img.Width
     let h = img.Height
     let imgData = img.Data
+    let se = [| -1, 0; 0, -1; 1, 0; 0, 1 |]
 
-    let areas = List<Area>(findMaxima img |> Seq.map Area)
+    let areas = List<Area>((if op = AreaOperation.Opening then findMaxima img else findMinima img) |> Seq.map Area)
 
     let pixels: Area[,] = Array2D.create h w null
     for m in areas do
@@ -200,16 +240,13 @@ let areaOpen (img: Image<Gray, byte>) (area: int) =
 
     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'
+            for i, j in se do
+                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
@@ -219,12 +256,12 @@ let areaOpen (img: Image<Gray, byte>) (area: int) =
             queue.Clear()
             addEdgeToQueue m.Elements
 
-            let mutable intensity = queue.Max
+            let mutable intensity = if op = AreaOperation.Opening then queue.Max else queue.Min
             let nextElements = Points()
 
             let mutable stop = false
             while not stop do
-                let intensity', p = queue.Next ()
+                let intensity', p = if op = AreaOperation.Opening then queue.NextMax () else queue.NextMin ()
                 let mutable merged = false
 
                 if intensity' = intensity // The intensity doesn't change.
@@ -237,7 +274,7 @@ let areaOpen (img: Image<Gray, byte>) (area: int) =
                     else
                         nextElements.Add(p) |> ignore
 
-                elif intensity' < intensity
+                elif if op = AreaOperation.Opening then intensity' < intensity else intensity' > intensity
                 then
                     m.Elements.UnionWith(nextElements)
                     for e in nextElements do
@@ -253,7 +290,7 @@ let areaOpen (img: Image<Gray, byte>) (area: int) =
                         nextElements.Clear()
                         nextElements.Add(p) |> ignore
 
-                else // i' > i
+                else
                     let m' = pixels.[p.Y, p.X]
                     if m' <> null
                     then
@@ -265,7 +302,7 @@ let areaOpen (img: Image<Gray, byte>) (area: int) =
                                 queue.Remove imgData.[e.Y, e.X, 0] e
                             addEdgeToQueue m'.Elements
                             m.Elements.UnionWith(m'.Elements)
-                            let intensityMax = queue.Max
+                            let intensityMax = if op = AreaOperation.Opening then queue.Max else queue.Min
                             if intensityMax <> intensity
                             then
                                 intensity <- intensityMax
@@ -280,21 +317,18 @@ let areaOpen (img: Image<Gray, byte>) (area: int) =
 
                 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'
+                    for i, j in se do
+                        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
@@ -316,6 +350,67 @@ let areaOpen (img: Image<Gray, byte>) (area: int) =
     ()
 
 
+let areaOpen (img: Image<Gray, byte>) (area: int) =
+    areaOperation img area AreaOperation.Opening
+
+let areaClose (img: Image<Gray, byte>) (area: int) =
+    areaOperation img area AreaOperation.Closing
+
+let areaOpen2 (img: Image<Gray, byte>) (area: int) =
+    let w = img.Width
+    let h = img.Height
+    let imgData = img.Data
+    let se = [| -1, 0; 0, -1; 1, 0; 0, 1 |]
+
+    let histogram = Array.zeroCreate 256
+    for i in 0 .. h - 1 do
+        for j in 0 .. w - 1 do
+            let v = imgData.[i, j, 0] |> int
+            histogram.[v] <- histogram.[v] + 1
+
+    let flooded : bool[,] = Array2D.zeroCreate h w
+
+    let pointsChecked = HashSet<Point>()
+    let pointsToCheck = Stack<Point>()
+
+    for level in 255 .. -1 .. 0 do
+        let mutable n = histogram.[level]
+        if n > 0
+        then
+            for i in 0 .. h - 1 do
+                for j in 0 .. w - 1 do
+                    if not flooded.[i, j] && imgData.[i, j, 0] = byte level
+                    then
+                        let mutable maxNeighborValue = 0uy
+                        pointsChecked.Clear()
+                        pointsToCheck.Clear()
+                        pointsToCheck.Push(Point(j, i))
+
+                        while pointsToCheck.Count > 0 do
+                            let next = pointsToCheck.Pop()
+                            pointsChecked.Add(next) |> ignore
+                            flooded.[next.Y, next.X] <- true
+
+                            for nx, ny in se do
+                                let p = Point(next.X + nx, next.Y + ny)
+                                if p.X >= 0 && p.X < w && p.Y >= 0 && p.Y < h
+                                then
+                                    let v = imgData.[p.Y, p.X, 0]
+                                    if v = byte level
+                                    then
+                                        if not (pointsChecked.Contains(p))
+                                        then
+                                            pointsToCheck.Push(p)
+                                    elif v > maxNeighborValue
+                                    then
+                                        maxNeighborValue <- v
+
+                        if int maxNeighborValue < level && pointsChecked.Count <= area
+                        then
+                            for p in pointsChecked do
+                                imgData.[p.Y, p.X, 0] <- maxNeighborValue
+
+
 // Zhang and Suen algorithm.
 // Modify 'mat' in place.
 let thin (mat: Matrix<byte>) =