Change the parasite detection method.
authorGreg Burri <greg.burri@gmail.com>
Mon, 21 Dec 2015 22:34:17 +0000 (23:34 +0100)
committerGreg Burri <greg.burri@gmail.com>
Mon, 21 Dec 2015 22:34:17 +0000 (23:34 +0100)
Parasitemia/Parasitemia/Classifier.fs
Parasitemia/Parasitemia/Config.fs
Parasitemia/Parasitemia/Ellipse.fs
Parasitemia/Parasitemia/Heap.fs [new file with mode: 0644]
Parasitemia/Parasitemia/ImgTools.fs
Parasitemia/Parasitemia/MainAnalysis.fs
Parasitemia/Parasitemia/MatchingEllipses.fs
Parasitemia/Parasitemia/Parasitemia.fsproj
Parasitemia/Parasitemia/ParasitesMarker.fs
Parasitemia/Parasitemia/ParasitesMarker2.fs [new file with mode: 0644]
Parasitemia/Parasitemia/Program.fs

index 3ed42f0..4c6cc56 100644 (file)
@@ -21,7 +21,7 @@ type private EllipseFlaggedKd (e: Ellipse) =
         member this.Y = this.Cy
 
 
-let findCells (ellipses: Ellipse list) (parasites: ParasitesMarker.Result) (fg: Image<Gray, byte>) (config: Config.Config) : Cell list =
+let findCells (ellipses: Ellipse list) (parasites: ParasitesMarker2.Result) (img: Image<Gray, byte>) (config: Config.Config) : Cell list =
     if ellipses.IsEmpty
     then
         []
@@ -29,10 +29,10 @@ let findCells (ellipses: Ellipse list) (parasites: ParasitesMarker.Result) (fg:
         let infection = parasites.infection.Copy() // To avoid to modify the parameter.
 
         // This is the minimum window size to check if other ellipses touch 'e'.
-        let searchRegion (e: Ellipse) = { KdTree.minX = e.Cx - (e.A + config.maxRBCSize) * config.scale
-                                          KdTree.maxX = e.Cx + (e.A + config.maxRBCSize) * config.scale
-                                          KdTree.minY = e.Cy - (e.A + config.maxRBCSize) * config.scale
-                                          KdTree.maxY = e.Cy + (e.A + config.maxRBCSize) * config.scale }
+        let searchRegion (e: Ellipse) = { KdTree.minX = e.Cx - (e.A + config.RBCMax) * config.Parameters.scale
+                                          KdTree.maxX = e.Cx + (e.A + config.RBCMax) * config.Parameters.scale
+                                          KdTree.minY = e.Cy - (e.A + config.RBCMax) * config.Parameters.scale
+                                          KdTree.maxY = e.Cy + (e.A + config.RBCMax) * config.Parameters.scale }
 
         // The minimum window to contain a given ellipse.
         let ellipseWindow (e: Ellipse) =
@@ -40,8 +40,10 @@ let findCells (ellipses: Ellipse list) (parasites: ParasitesMarker.Result) (fg:
             let a = int (e.A + 0.5)
             cx - a, cy - a, cx + a, cy + a
 
-        let w = float fg.Width
-        let h = float fg.Height
+        let w = img.Width
+        let w_f = float w
+        let h = img.Height
+        let h_f = float h
 
         let ellipses = ellipses |> List.map EllipseFlaggedKd
 
@@ -60,7 +62,7 @@ let findCells (ellipses: Ellipse list) (parasites: ParasitesMarker.Result) (fg:
         let ellipsesWithNeigbors = ellipses |> List.choose (fun e -> if e.Removed then None else Some (e, neighbors e))
 
         // 2) Remove ellipses with a lower percentage of foreground. (taken from the lower score to the highest).
-        for e, neighbors in List.rev ellipsesWithNeigbors do
+        (*for e, neighbors in List.rev ellipsesWithNeigbors do
             let minX, minY, maxX, maxY = ellipseWindow e
 
             let mutable totalElement = 0
@@ -75,13 +77,13 @@ let findCells (ellipses: Ellipse list) (parasites: ParasitesMarker.Result) (fg:
                             then
                                 fgElement <- fgElement + 1
 
-            if (float fgElement) / (float totalElement) < config.percentageOfFgValidCell
+            if (float fgElement) / (float totalElement) < config.Parameters.percentageOfFgValidCell
             then
-                e.Removed <- true
+                e.Removed <- true*)
 
         // 3) Remove ellipses touching the edges.
         for e in ellipses do
-            if e.isOutside w h then e.Removed <- true
+            if e.isOutside w_f h_f then e.Removed <- true
 
         // 4) Remove ellipses with little area.
         for e, neighbors in ellipsesWithNeigbors do
@@ -90,16 +92,15 @@ let findCells (ellipses: Ellipse list) (parasites: ParasitesMarker.Result) (fg:
                 let minX, minY, maxX, maxY = ellipseWindow e
 
                 let mutable area = 0
-                for y in (if minY < 0 then 0 else minY) .. (if maxY >= fg.Height then fg.Height - 1 else maxY) do
-                    for x in (if minX < 0 then 0 else minX) .. (if maxX >= fg.Width then fg.Width - 1 else maxX) do
+                for y in (if minY < 0 then 0 else minY) .. (if maxY >= h then h - 1 else maxY) do
+                    for x in (if minX < 0 then 0 else minX) .. (if maxX >= w then w - 1 else maxX) do
                         let yf, xf = float y, float x
-                        if fg.Data.[y, x, 0] > 0uy &&
-                           e.Contains xf yf &&
+                        if e.Contains xf yf &&
                            neighbors |> List.forall (fun (otherE, _, _) -> otherE.Removed || not <| otherE.Contains xf yf)
                             then
                                 area <- area + 1
 
-                if area < config.minimumCellArea
+                if area < int config.RBCMinArea
                 then
                     e.Removed <- true
 
@@ -142,8 +143,7 @@ let findCells (ellipses: Ellipse list) (parasites: ParasitesMarker.Result) (fg:
                 for y in minY .. maxY do
                     for x in minX .. maxX do
                         let p = PointD(float x, float y)
-                        if pixelOwnedByE p e (neighbors |> List.choose (fun (otherE, p1, p2) -> if otherE.Removed then None else Some (Utils.lineFromTwoPoints p1 p2))) &&
-                           fg.Data.[y, x, 0] > 0uy
+                        if pixelOwnedByE p e (neighbors |> List.choose (fun (otherE, p1, p2) -> if otherE.Removed then None else Some (Utils.lineFromTwoPoints p1 p2)))
                         then
                             elements.[y-minY, x-minX] <- 1uy
                             nbElement <- nbElement + 1
@@ -167,11 +167,11 @@ let findCells (ellipses: Ellipse list) (parasites: ParasitesMarker.Result) (fg:
                                 darkStainPixels <- darkStainPixels + 1
 
                 let cellClass =
-                    if float darkStainPixels > config.MaxDarkStainRatio * (float nbElement) (* ||
+                    if float darkStainPixels > config.Parameters.maxDarkStainRatio * (float nbElement) (* ||
                         sqrt (((float sumCoords_x) / (float nbElement) - e.Cx) ** 2.0 + ((float sumCoords_y) / (float nbElement) - e.Cy) ** 2.0) > e.A * config.maxOffcenter *)
                     then
                         Peculiar
-                    elif infectedPixels.Count > config.infectionPixelsRequired
+                    elif infectedPixels.Count > config.Parameters.infectionPixelsRequired
                     then
                         let infectionToRemove = ImgTools.connectedComponents parasites.stain infectedPixels
                         for p in infectionToRemove do
index 48a6334..55e142c 100644 (file)
@@ -1,16 +1,16 @@
 module Config
 
+open System
+
 type Debug =
     | DebugOff
     | DebugOn of string // Output directory.
 
-type Config = {
-    debug: Debug
-
+type Parameters = {
     scale: float
 
-    minRBCSize: float
-    maxRBCSize: float
+    minRbcRadius: float
+    maxRbcRadius: float
 
     preFilterSigma: float
 
@@ -21,18 +21,27 @@ type Config = {
     // Parasites detection.
     darkStainLevel: float
 
-    stainSigma: float
-    stainLevel: float
-    stainSpreadRequired: float
+    stainArea: float // Factor of a RBC area. 0.5 means the half of RBC area.
+    stainLevel: float // [0, 1]
 
-    infectionSigma: float
-    infectionLevel: float
+    infectionArea: float // Factor of a RBC area. 0.5 means the half of RBC area.
+    infectionLevel: float // [0, 1]
     infectionPixelsRequired: int
 
-    percentageOfFgValidCell: float
+    maxDarkStainRatio: float
+
+    minimumCellArea: float // Factor of RBC area.
+}
+
+type Config (param: Parameters) =
+    member this.Parameters = param
+    member val Debug = DebugOff with get, set
+    member val RBCSize = 30. with get, set
+
+    member this.RBCMin = this.RBCSize + param.minRbcRadius * this.RBCSize
+    member this.RBCMax = this.RBCSize + param.maxRbcRadius * this.RBCSize
 
-    MaxDarkStainRatio: float
+    member this.RBCMinArea = param.minimumCellArea * Math.PI * this.RBCSize ** 2.0
 
-    minimumCellArea: int
-    maxOffcenter: float
-}
\ No newline at end of file
+    member this.ParasiteArea = param.infectionArea * Math.PI * this.RBCSize ** 2.0
+    member this.StainArea = param.stainArea * Math.PI * this.RBCSize ** 2.0
index 139b7b5..33333bf 100644 (file)
@@ -46,8 +46,8 @@ let private goldenSectionSearch (f: float -> float) (nbIter: int) (xmin: float)
 // Ellipse.A is always equal or greater than Ellipse.B.
 // Ellipse.Alpha is between 0 and Pi.
 let ellipse (p1x: float) (p1y: float) (m1: float) (p2x: float) (p2y: float) (m2: float) (p3x: float) (p3y: float) : Types.Ellipse option =
-    let accuracy_extremum_search_1 = 7 // 3
-    let accuracy_extremum_search_2 = 7 // 4
+    let accuracy_extremum_search_1 = 8 // 3
+    let accuracy_extremum_search_2 = 8 // 4
 
     // p3 as the referencial.
     let p1x = p1x - p3x
@@ -208,9 +208,9 @@ let find (edges: Matrix<byte>)
          (yDir: Image<Gray, float>)
          (config: Config) : MatchingEllipses =
 
-    let r1, r2 = config.scale * config.minRBCSize, config.scale * config.maxRBCSize
-    let windowSize = roundInt (config.factorWindowSize * r2)
-    let factorNbPick = config.factorNbPick
+    let r1, r2 = config.Parameters.scale * config.RBCMin, config.Parameters.scale * config.RBCMax // FIXME: scale factor should be applied in Config!?
+    let windowSize = roundInt (config.Parameters.factorWindowSize * r2)
+    let factorNbPick = config.Parameters.factorNbPick
 
     let increment = windowSize / 4
 
diff --git a/Parasitemia/Parasitemia/Heap.fs b/Parasitemia/Parasitemia/Heap.fs
new file mode 100644 (file)
index 0000000..c2ab69f
--- /dev/null
@@ -0,0 +1,41 @@
+module Heap
+
+open System.Collections.Generic
+
+let inline private parent (i: int) : int = (i - 1) / 2
+let inline private left (i: int) : int = 2 * (i + 1) - 1
+let inline private right (i: int) : int = 2 * (i + 1)
+
+type private Node<'k, 'v> = { key: 'k; value : 'v } with
+    override this.ToString () = sprintf "%A -> %A" this.key this.value
+
+type Heap<'k, 'v when 'k : comparison> () =
+    let a = List<Node<'k, 'v>>()
+
+    let rec heapUp (i: int) =
+        let l, r = left i, right i
+        let mutable max = if l < a.Count && a.[l].key > a.[i].key then l else i
+        if r < a.Count && a.[r].key > a.[max].key
+        then
+            max <- r
+        if max <> i
+        then
+            let tmp = a.[i]
+            a.[i] <- a.[max]
+            a.[max] <- tmp
+            heapUp max
+
+    member this.Next () : 'k * 'v =
+        let { key = key; value = value} = a.[0]
+        a.RemoveAt(0)
+        key, value
+
+    member this.Add (key: 'k) (value: 'v) =
+        a.Insert(0, { key = key; value = value})
+        heapUp 0
+
+    member this.IsEmpty = a.Count = 0
+    member this.Max : 'k = a.[0].key
+    member this.Clear () = a.Clear()
+
+
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>) =
index 39383fc..d959624 100644 (file)
@@ -14,31 +14,40 @@ open Types
 
 let doAnalysis (img: Image<Bgr, byte>) (name: string) (config: Config) : Cell list =
 
-    use scaledImg = if config.scale = 1.0 then img else img.Resize(config.scale, CvEnum.Inter.Area)
+    use scaledImg = if config.Parameters.scale = 1.0 then img else img.Resize(config.Parameters.scale, CvEnum.Inter.Area)
 
-    use blue  = scaledImg.Item(0)
     use green = scaledImg.Item(1)
-    use red = scaledImg.Item(2)
 
     let greenFloat = green.Convert<Gray, float32>()
 
-    let filteredGreen = gaussianFilter green config.preFilterSigma
-    logTime "areaOpen" (fun () -> ImgTools.areaOpen filteredGreen 2000)
+    let filteredGreen = gaussianFilter green config.Parameters.preFilterSigma
+
+    (*let maximaImg = filteredGreen.Copy()
+    let maxima = logTime "maxima" (fun () -> ImgTools.findMaxima maximaImg)
+    for m in maxima do
+        ImgTools.drawPoints maximaImg m 255uy
 
-    let RBCSize = Granulometry.findRadius filteredGreen (10, 100) 0.5 |> float
-    let RBCSizeDelta = RBCSize / 3.0
+    let greenOpen1 = filteredGreen.Copy()
+    logTime "areaOpen1" (fun () -> ImgTools.areaOpen greenOpen1 2000)*)
+
+    logTime "areaOpen" (fun () -> ImgTools.areaOpen filteredGreen 2000)
 
-    let config = { config with minRBCSize = RBCSize - RBCSizeDelta ; maxRBCSize = RBCSize + RBCSizeDelta }
+    config.RBCSize <- Granulometry.findRadius filteredGreen (10, 100) 0.5 |> float
 
     let filteredGreenFloat = filteredGreen.Convert<Gray, float32>() // Is it neccessary?
 
+    let kmediansResults = logTime "Finding foreground (k-medians)" (fun () -> KMedians.kmedians filteredGreenFloat 1.0)
+
+    let parasites, filteredGreenWhitoutParasites, filteredGreenWhitoutStain = ParasitesMarker2.find filteredGreen filteredGreenFloat kmediansResults config
+    let filteredGreenWhitoutParasitesFloat = filteredGreenWhitoutParasites.Convert<Gray, float32>()
+
     use sobelKernel =
         new ConvolutionKernelF(array2D [[ 1.0f; 0.0f; -1.0f ]
                                         [ 2.0f; 0.0f; -2.0f ]
                                         [ 1.0f; 0.0f; -1.0f ]], Point(0, 0))
 
-    use xEdges = filteredGreenFloat.Convolution(sobelKernel).Convert<Gray, float>()
-    use yEdges = filteredGreenFloat.Convolution(sobelKernel.Transpose()).Convert<Gray, float>()
+    use xEdges = filteredGreenWhitoutParasitesFloat.Convolution(sobelKernel).Convert<Gray, float>()
+    use yEdges = filteredGreenWhitoutParasitesFloat.Convolution(sobelKernel.Transpose()).Convert<Gray, float>()
 
     let xEdgesData = xEdges.Data
     let yEdgesData = yEdges.Data
@@ -70,18 +79,14 @@ let doAnalysis (img: Image<Bgr, byte>) (name: string) (config: Config) : Cell li
     logTime "Finding edges" (fun() -> thin edges)
     logTime "Removing small connected components from thinning" (fun () -> removeArea edges 12)
 
-    let kmediansResults = logTime "Finding foreground (k-medians)" (fun () -> KMedians.kmedians filteredGreenFloat 1.0)
-
-    let parasites = ParasitesMarker.find greenFloat filteredGreenFloat kmediansResults config
-
     let allEllipses, ellipses = logTime "Finding ellipses" (fun () ->
         let matchingEllipses = Ellipse.find edges xEdges yEdges config
         matchingEllipses.Ellipses, matchingEllipses.PrunedEllipses )
 
-    let cells = logTime "Classifier" (fun () -> Classifier.findCells ellipses parasites kmediansResults.fg config)
+    let cells = logTime "Classifier" (fun () -> Classifier.findCells ellipses parasites filteredGreenWhitoutParasites config)
 
     // Output pictures if debug flag is set.
-    match config.debug with
+    match config.Debug with
     | DebugOn output ->
         let dirPath = System.IO.Path.Combine(output, name)
         System.IO.Directory.CreateDirectory dirPath |> ignore
@@ -112,14 +117,23 @@ let doAnalysis (img: Image<Bgr, byte>) (name: string) (config: Config) : Cell li
         drawCells imgCells' true cells
         saveImg imgCells' (buildFileName " - cells - full.png")
 
-        let filteredGreenMaxima = gaussianFilter green config.preFilterSigma
+        let filteredGreenMaxima = gaussianFilter green config.Parameters.preFilterSigma
         for m in ImgTools.findMaxima filteredGreenMaxima do
             ImgTools.drawPoints filteredGreenMaxima m 255uy
         saveImg filteredGreenMaxima (buildFileName " - filtered - maxima.png")
 
         saveImg filteredGreen (buildFileName " - filtered.png")
-        saveImg blue (buildFileName " - blue.png")
+        saveImg filteredGreenWhitoutParasites (buildFileName " - filtered closed.png")
+
+        (*saveImg parasitesMarker (buildFileName " - parasites (area closing).png")
+        saveImg stainMarker (buildFileName " - stain (area closing).png")*)
+
         saveImg green (buildFileName " - green.png")
+
+        use blue  = scaledImg.Item(0)
+        saveImg blue (buildFileName " - blue.png")
+
+        use red = scaledImg.Item(2)
         saveImg red (buildFileName " - red.png")
     | _ -> ()
 
index 99e3cb5..2196122 100644 (file)
@@ -13,19 +13,17 @@ open Utils
 let matchingScoreThreshold1 = 0.6
 
 // All ellipsee with a score below this is removed.
-let matchingScoreThreshold2 = 1. / 100.
+let matchingScoreThreshold2 = 2.
 
 type private EllipseScoreFlaggedKd (matchingScore: float, e: Ellipse) =
     let mutable matchingScore = matchingScore
-    let perimeter = e.Perimeter
 
     member this.Ellipse = e
 
     member this.MatchingScore = matchingScore
 
-    // The score is proportional to the perimeter because large ellipse will receive more votes.
     member this.AddMatchingScore(score: float) =
-        matchingScore <- matchingScore + score / perimeter
+        matchingScore <- matchingScore + score
 
     member val Processed = false with get, set
     member val Removed = false with get, set
index 3e6609a..fe81290 100644 (file)
@@ -28,7 +28,7 @@
     <PlatformTarget>x64</PlatformTarget>
     <DocumentationFile>bin\Debug\Parasitemia.XML</DocumentationFile>
     <Prefer32Bit>false</Prefer32Bit>
-    <StartArguments>--folder "../../../Images/areaopen" --output "../../../Images/output" --debug</StartArguments>
+    <StartArguments>--folder "../../../Images/debug" --output "../../../Images/output" --debug</StartArguments>
   </PropertyGroup>
   <PropertyGroup Condition=" '$(Configuration)|$(Platform)' == 'Release|AnyCPU' ">
     <DebugType>pdbonly</DebugType>
@@ -72,6 +72,7 @@
     <Compile Include="Config.fs" />
     <Compile Include="KMedians.fs" />
     <Compile Include="EEOver.fs" />
+    <Compile Include="ParasitesMarker2.fs" />
     <Compile Include="ParasitesMarker.fs" />
     <Compile Include="KdTree.fs" />
     <Compile Include="MatchingEllipses.fs" />
index 9adf185..d1a51fa 100644 (file)
@@ -14,11 +14,11 @@ type Result = {
 // * 'Dark stain' corresponds to the colored pixel, it's independent of the size of the areas.
 // * 'Stain' corresponds to the stain around the parasites.
 // * 'Infection' corresponds to the parasite. It shouldn't contain thrombocytes.
-let find (green: Image<Gray, float32>) (filteredGreen: Image<Gray, float32>) (kmediansResult: KMedians.Result) (config: Config.Config) : Result =
+(*let find (green: Image<Gray, float32>) (filteredGreen: Image<Gray, float32>) (kmediansResult: KMedians.Result) (config: Config.Config) : Result =
 
     // We use the filtered image to find the dark stain.
     let { KMedians.fg = fg; KMedians.median_bg = median_bg; KMedians.median_fg = median_fg; KMedians.d_fg = d_fg } = kmediansResult
-    let darkStain = d_fg.Cmp(median_bg * config.darkStainLevel, CvEnum.CmpType.GreaterThan)
+    let darkStain = d_fg.Cmp(median_bg * config.Parameters.darkStainLevel, CvEnum.CmpType.GreaterThan)
     darkStain._And(filteredGreen.Cmp(median_fg, CvEnum.CmpType.LessThan))
     darkStain._And(fg)
 
@@ -34,8 +34,8 @@ let find (green: Image<Gray, float32>) (filteredGreen: Image<Gray, float32>) (km
         smears
 
     { darkStain = darkStain;
-      stain = findSmears config.stainSigma config.stainLevel
-      infection = findSmears config.infectionSigma config.infectionLevel }
+      stain = findSmears config.Parameters.stainSigma config.Parameters.stainLevel
+      infection = findSmears config.Parameters.infectionSigma config.Parameters.infectionLevel }*)
 
 
 
diff --git a/Parasitemia/Parasitemia/ParasitesMarker2.fs b/Parasitemia/Parasitemia/ParasitesMarker2.fs
new file mode 100644 (file)
index 0000000..03af2b6
--- /dev/null
@@ -0,0 +1,50 @@
+module ParasitesMarker2
+
+open System.Drawing
+
+open Emgu.CV
+open Emgu.CV.Structure
+
+type Result = {
+    darkStain: Image<Gray, byte>
+    stain: Image<Gray, byte>
+    infection: Image<Gray, byte> }
+
+// Create three binary markers :
+// * 'Dark stain' corresponds to the colored pixel, it's independent of the size of the areas.
+// * 'Stain' corresponds to the stain around the parasites.
+// * 'Infection' corresponds to the parasite. It shouldn't contain thrombocytes.
+let find (filteredGreen: Image<Gray, byte>) (filteredGreenFloat: Image<Gray, float32>) (kmediansResult: KMedians.Result) (config: Config.Config) : Result * Image<Gray, byte> * Image<Gray, byte> =
+
+    // We use the filtered image to find the dark stain.
+    let { KMedians.fg = fg; KMedians.median_bg = median_bg; KMedians.median_fg = median_fg; KMedians.d_fg = d_fg } = kmediansResult
+    let darkStain = d_fg.Cmp(median_bg * config.Parameters.darkStainLevel, CvEnum.CmpType.GreaterThan)
+    darkStain._And(filteredGreenFloat.Cmp(median_fg, CvEnum.CmpType.LessThan))
+
+    let marker (area: int) (threshold: float) : Image<Gray, byte> * Image<Gray, byte> =
+        let closed = filteredGreen.Copy()
+        ImgTools.areaClose closed area
+        let diff = closed - filteredGreen
+
+        let min = ref [| 0. |]
+        let minLocation = ref <| [| Point() |]
+        let max = ref [| 0. |]
+        let maxLocation = ref <| [| Point() |]
+        diff.MinMax(min, max, minLocation, maxLocation)
+
+        diff._ThresholdBinary((!max).[0] * threshold |> Gray, Gray(255.))
+        diff, closed
+
+    let parasitesMarker, filteredGreenWithoutParasites  = marker (int config.ParasiteArea) config.Parameters.infectionLevel
+    let stainMarker, filteredGreenWithoutStain = marker (int config.StainArea) config.Parameters.stainLevel
+
+    { darkStain = darkStain
+      stain = parasitesMarker
+      infection = stainMarker },
+    filteredGreenWithoutParasites,
+    filteredGreenWithoutStain
+
+
+
+
+
index 86de435..103596a 100644 (file)
@@ -57,41 +57,36 @@ let parseArgs (args: string[]) : Arguments =
 let main args =
     match parseArgs args with
     | mode, debug ->
-        let scale = 1.
-        let config = {
-            debug = if debug then DebugOn "." else DebugOff
+        let config =
+            Config(
+              { scale = 1.
 
-            scale = scale
+                minRbcRadius = -0.3
+                maxRbcRadius = 0.3
 
-            // RBC size range in px at scale = 1.0.
-            minRBCSize = 20.
-            maxRBCSize = 42.
+                preFilterSigma = 1.5
 
-            preFilterSigma = 1.5
+                factorNbPick = 1.0
+                factorWindowSize = 2.0
 
-            factorNbPick = 1.0 // 1.0
-            factorWindowSize = 2.0
+                darkStainLevel = 0.4 // Lower -> more sensitive.
 
-            darkStainLevel = 0.4 // Lower -> more sensitive.
+                stainArea = 0.02
+                stainLevel = 0.2
 
-            stainSigma = 10.
-            stainLevel = 0.9
-            stainSpreadRequired = 3.0
+                infectionArea = 0.06
+                infectionLevel = 0.2
+                infectionPixelsRequired = 1
 
-            infectionSigma = 2.2
-            infectionLevel = 0.87
-            infectionPixelsRequired = 1
+                maxDarkStainRatio = 0.1
 
-            percentageOfFgValidCell = 0.4
-
-            MaxDarkStainRatio = 0.1
-
-            minimumCellArea = 1000. * scale ** 2. |> int
-            maxOffcenter = 0.5 }
+                minimumCellArea = 0.3 })
 
         match mode with
         | CmdLine (input, output) ->
-            let config = if debug then { config with debug = DebugOn output } else config
+            if debug
+            then
+                config.Debug <- DebugOn output
 
             Directory.CreateDirectory output |> ignore
 
@@ -121,6 +116,10 @@ let main args =
             let app = new Application()
             let mainWindow = Views.MainWindow()
 
+            if debug
+            then
+                config.Debug <- DebugOn "."
+
             Utils.log <- (fun m -> log mainWindow m)
 
             //display mainWindow img