Add functions to apply an area opening and to compute granulometry.
authorGreg Burri <greg.burri@gmail.com>
Fri, 18 Dec 2015 19:27:24 +0000 (20:27 +0100)
committerGreg Burri <greg.burri@gmail.com>
Fri, 18 Dec 2015 19:27:24 +0000 (20:27 +0100)
Parasitemia/Parasitemia/Classifier.fs
Parasitemia/Parasitemia/Ellipse.fs
Parasitemia/Parasitemia/Granulometry.fs [new file with mode: 0644]
Parasitemia/Parasitemia/ImgTools.fs
Parasitemia/Parasitemia/KMedians.fs
Parasitemia/Parasitemia/MainAnalysis.fs
Parasitemia/Parasitemia/MatchingEllipses.fs
Parasitemia/Parasitemia/Parasitemia.fsproj
Parasitemia/Parasitemia/Program.fs
Parasitemia/Parasitemia/Types.fs

index 8467d24..3ed42f0 100644 (file)
@@ -40,13 +40,13 @@ 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 a = int (e.A + 0.5)
             cx - a, cy - a, cx + a, cy + a
 
-        // 1) Remove ellipses touching the edges.
         let w = float fg.Width
         let h = float fg.Height
         let w = float fg.Width
         let h = float fg.Height
-        let ellipsesInside = ellipses |> List.map (fun e -> EllipseFlaggedKd (e, Removed = e.isOutside w h))
 
 
-        // 2) Associate touching ellipses with each ellipses.
-        let tree = KdTree.Tree.BuildTree ellipsesInside
+        let ellipses = ellipses |> List.map EllipseFlaggedKd
+
+        // 1) Associate touching ellipses with each ellipses.
+        let tree = KdTree.Tree.BuildTree ellipses
         let neighbors (e: Ellipse) : (EllipseFlaggedKd * PointD * PointD) list =
             tree.Search (searchRegion e)
                 // We only keep the ellipses touching 'e'.
         let neighbors (e: Ellipse) : (EllipseFlaggedKd * PointD * PointD) list =
             tree.Search (searchRegion e)
                 // We only keep the ellipses touching 'e'.
@@ -57,10 +57,10 @@ let findCells (ellipses: Ellipse list) (parasites: ParasitesMarker.Result) (fg:
                     | _ ->
                         None )
 
                     | _ ->
                         None )
 
-        let ellipsesWithNeigbors = ellipsesInside |> List.choose (fun e -> if e.Removed then None else Some (e, neighbors e))
+        let ellipsesWithNeigbors = ellipses |> List.choose (fun e -> if e.Removed then None else Some (e, neighbors e))
 
 
-        // 3) Remove ellipses with a lower percentage of foreground.
-        for e, neighbors in ellipsesWithNeigbors do
+        // 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
             let minX, minY, maxX, maxY = ellipseWindow e
 
             let mutable totalElement = 0
             let minX, minY, maxX, maxY = ellipseWindow e
 
             let mutable totalElement = 0
@@ -68,7 +68,7 @@ let findCells (ellipses: Ellipse list) (parasites: ParasitesMarker.Result) (fg:
             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
                     let yf, xf = float y, float x
             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
                     let yf, xf = float y, float x
-                    if e.Contains xf yf && neighbors |> List.forall (fun (otherE, _, _) -> not <| otherE.Contains xf yf)
+                    if e.Contains xf yf && neighbors |> List.forall (fun (otherE, _, _) -> otherE.Removed || not <| otherE.Contains xf yf)
                         then
                             totalElement <- totalElement + 1
                             if fg.Data.[y, x, 0] > 0uy
                         then
                             totalElement <- totalElement + 1
                             if fg.Data.[y, x, 0] > 0uy
@@ -79,6 +79,10 @@ let findCells (ellipses: Ellipse list) (parasites: ParasitesMarker.Result) (fg:
             then
                 e.Removed <- true
 
             then
                 e.Removed <- true
 
+        // 3) Remove ellipses touching the edges.
+        for e in ellipses do
+            if e.isOutside w h then e.Removed <- true
+
         // 4) Remove ellipses with little area.
         for e, neighbors in ellipsesWithNeigbors do
             if not e.Removed
         // 4) Remove ellipses with little area.
         for e, neighbors in ellipsesWithNeigbors do
             if not e.Removed
index 128c823..139b7b5 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 =
 // 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 = 4
-    let accuracy_extremum_search_2 = 3
+    let accuracy_extremum_search_1 = 7 // 3
+    let accuracy_extremum_search_2 = 7 // 4
 
     // p3 as the referencial.
     let p1x = p1x - p3x
 
     // p3 as the referencial.
     let p1x = p1x - p3x
@@ -178,18 +178,31 @@ let private areVectorsValid (p1x: float) (p1y: float) (p2x: float) (p2y: float)
 
     let b1 = -m1 * p1x + p1y
     let b2 = -m2 * p2x + p2y
 
     let b1 = -m1 * p1x + p1y
     let b2 = -m2 * p2x + p2y
-    let px = -((b1 - b2)/(m1 - m2))
-    let py = -((m2 * b1 - m1 * b2)/(m1 - m2))
+    let px = -((b1 - b2) / (m1 - m2))
+    let py = -((m2 * b1 - m1 * b2) / (m1 - m2))
 
     let rot1 = vectorRotation p1x p1y v1x v1y px py
     let rot2 = vectorRotation p2x p2y v2x v2y px py
 
 
     let rot1 = vectorRotation p1x p1y v1x v1y px py
     let rot2 = vectorRotation p2x p2y v2x v2y px py
 
-    if rot1 = rot2 || rot1 * atan2 (p1y - py) (p1x - px) + rot2 * atan2 (p2y - py) (p2x - px) <= 0.0
+    if rot1 = rot2
     then
         None
     else
     then
         None
     else
+        let alpha1 = atan2 (p1y - py) (p1x - px)
+        let alpha2 = atan2 (p2y - py) (p2x - px)
+
+        let alpha1' = if alpha1 < 0.0 then 2.0 * Math.PI + alpha1 else alpha1
+        let alpha2' = if alpha2 < 0.0 then 2.0 * Math.PI + alpha2 else alpha2
+
+        let diff = rot1 * alpha1' + rot2 * alpha2'
+
+        if diff > Math.PI || (diff < 0.0 && diff > -Math.PI)
+        then
+            None
+        else
         Some (m1, m2)
 
         Some (m1, m2)
 
+
 let find (edges: Matrix<byte>)
          (xDir: Image<Gray, float>)
          (yDir: Image<Gray, float>)
 let find (edges: Matrix<byte>)
          (xDir: Image<Gray, float>)
          (yDir: Image<Gray, float>)
@@ -262,6 +275,8 @@ let find (edges: Matrix<byte>)
                                 match ellipse p1xf p1yf m1 p2xf p2yf m2 p3xf p3yf with
                                 | Some e when e.Cx > 0.0 && e.Cx < (float w) - 1.0 && e.Cy > 0.0 && e.Cy < (float h) - 1.0 &&
                                               e.A >= r1 - radiusTolerance && e.A <= r2 + radiusTolerance && e.B >= r1 - radiusTolerance && e.B <= r2 + radiusTolerance ->
                                 match ellipse p1xf p1yf m1 p2xf p2yf m2 p3xf p3yf with
                                 | Some e when e.Cx > 0.0 && e.Cx < (float w) - 1.0 && e.Cy > 0.0 && e.Cy < (float h) - 1.0 &&
                                               e.A >= r1 - radiusTolerance && e.A <= r2 + radiusTolerance && e.B >= r1 - radiusTolerance && e.B <= r2 + radiusTolerance ->
+
+                                     let prout = areVectorsValid p1xf p1yf p2xf p2yf -xDirData.[p1y, p1x, 0] -yDirData.[p1y, p1x, 0] -xDirData.[p2y, p2x, 0] -yDirData.[p2y, p2x, 0]
                                      ellipses.Add e
                                 | _ -> ()
                             | _ -> ()
                                      ellipses.Add e
                                 | _ -> ()
                             | _ -> ()
diff --git a/Parasitemia/Parasitemia/Granulometry.fs b/Parasitemia/Parasitemia/Granulometry.fs
new file mode 100644 (file)
index 0000000..a9a102a
--- /dev/null
@@ -0,0 +1,38 @@
+module Granulometry
+
+open System
+open System.Drawing
+
+open Emgu.CV
+open Emgu.CV.Structure
+
+open Utils
+
+// 'range': a minimum and maximum radius.
+// 'scale': <= 1.0, to speed up the process.
+let findRadius (img: Image<Gray, byte>)  (range: int * int) (scale: float) : int =
+    use scaledImg = if scale = 1.0 then img else img.Resize(scale, CvEnum.Inter.Area)
+
+    let r1, r2 = range
+    let r1', r2' = roundInt (float r1 * scale), roundInt (float r2 * scale)
+
+    let patternSpectrum = Array.zeroCreate (r2' - r1')
+    let intensityImg = scaledImg.GetSum().Intensity
+
+    let mutable previous_n = Double.NaN
+    for r in r1' .. r2' do
+        let se = CvInvoke.GetStructuringElement(CvEnum.ElementShape.Ellipse, Size(2 * r, 2 * r), Point(-1, -1))
+        use closed = scaledImg.MorphologyEx(CvEnum.MorphOp.Close, se, Point(-1, -1), 1, CvEnum.BorderType.Replicate, MCvScalar(0.0))
+
+        let n = 1.0 - closed.GetSum().Intensity / intensityImg
+
+        if not (Double.IsNaN previous_n)
+        then
+            patternSpectrum.[r - r1' - 1] <- abs (n - previous_n)
+        previous_n <- n
+
+    let max = patternSpectrum |> Array.indexed |> Array.fold (fun (iMax, sMax) (i, s) -> if s > sMax then (i, s) else (iMax, sMax)) (0, Double.MinValue)
+
+    0
+
+
index 7253390..e9cbbe9 100644 (file)
@@ -3,11 +3,13 @@
 open System
 open System.Drawing
 open System.Collections.Generic
 open System
 open System.Drawing
 open System.Collections.Generic
+open System.Linq
 
 open Emgu.CV
 open Emgu.CV.Structure
 
 open Utils
 
 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> =
 
 // Normalize image values between 0uy and 255uy.
 let normalizeAndConvert (img: Image<Gray, float32>) : Image<Gray, byte> =
@@ -18,10 +20,266 @@ 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>()
 
     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)
 
 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)
 
+
+let findMaxima (img: Image<Gray, byte>) : IEnumerable<HashSet<Point>> =
+    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 -> HashSet<Point>(l))
+
+
+type PriorityQueue () =
+    let q = List<HashSet<Point>>() // TODO: Check performance with an HasSet
+    let mutable highest = -1 // Value of the first elements of 'q'.
+
+    member this.Next () : byte * Point =
+        if this.IsEmpty
+        then
+            invalidOp "Queue is empty"
+        else
+            let l = q.[0]
+            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)
+                    highest <- highest - 1
+            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
+        then
+            for i in highest .. vi - 1  do
+                q.Insert(0, null)
+            highest <- vi
+        elif highest - vi >= q.Count
+        then
+            for i in 0 .. highest - vi - q.Count do
+                q.Add(null)
+
+        let pos = highest - vi
+        if q.[pos] = null
+        then
+            q.[pos] <- HashSet<Point>([p])
+        else
+            q.[pos].Add(p) |> ignore
+
+
+    member this.IsEmpty =
+        q.Count = 0
+
+    member this.Clear () =
+        while highest >= 0 do
+            q.[highest] <- null
+            highest <- highest - 1
+
+
+
+type MaximaState =  Uncertain | Validated | TooBig
+type Maxima = {
+    elements : HashSet<Point>
+    mutable intensity: byte option
+    mutable state: MaximaState }
+
+
+let areaOpen (img: Image<Gray, byte>) (area: int) =
+    let w = img.Width
+    let h = img.Height
+
+    let maxima = findMaxima img |> Seq.map (fun m -> { elements = m; intensity = None; state = Uncertain }) |> List.ofSeq
+    let toValidated = Stack<Maxima>(maxima)
+
+    while toValidated.Count > 0 do
+        let m = toValidated.Pop()
+        if m.elements.Count <= area
+        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
+
+            let mutable intensity = queue.Max
+            let nextElements = HashSet<Point>()
+
+            let mutable stop = false
+            while not stop do
+                let intensity', p = queue.Next ()
+
+                if intensity' = intensity // The intensity doesn't change.
+                then
+                    if m.elements.Count + nextElements.Count + 1 > area
+                    then
+                        m.state <- 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
+                    then
+                        m.state <- 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
+
+                    if not stop
+                    then
+                        m.state <- Validated
+                        m.intensity <- Some (intensity)
+                        stop <- true
+
+                if not stop
+                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 <- Validated
+                                    m.intensity <- Some (intensity)
+                                    stop <- true
+                                elif not (m.elements.Contains(p')) && not (nextElements.Contains(p'))
+                                then
+                                    queue.Add (img.Data.[ni, nj, 0]) p'
+
+                if queue.IsEmpty
+                then
+                    if m.elements.Count + nextElements.Count <= area
+                    then
+                        m.state <- Validated
+                        m.intensity <- Some intensity'
+                        m.elements.UnionWith(nextElements)
+                    stop <- true
+
+    for m in maxima do
+        if m.state = Validated
+        then
+            match m.intensity with
+            | Some i ->
+                for p in m.elements do
+                    img.Data.[p.Y, p.X, 0] <- i
+            | _ -> ()
+    ()
+
+
 // Zhang and Suen algorithm.
 // Modify 'mat' in place.
 let thin (mat: Matrix<byte>) =
 // Zhang and Suen algorithm.
 // Modify 'mat' in place.
 let thin (mat: Matrix<byte>) =
@@ -73,6 +331,7 @@ let thin (mat: Matrix<byte>) =
         data2 <- tmp
 
 
         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)
 let pop (l: List<'a>) : 'a =
     let n = l.[l.Count - 1]
     l.RemoveAt(l.Count - 1)
index 1822c75..9b5d50f 100644 (file)
@@ -10,7 +10,7 @@ type Result = {
     fg: Image<Gray, byte>
     median_bg: float
     median_fg: float
     fg: Image<Gray, byte>
     median_bg: float
     median_fg: float
-    d_fg: Image<Gray, float32> } // Distances to median_fg.
+    d_fg: Image<Gray, float32> } // Euclidean distances of the foreground to median_fg.
 
 let kmedians (img: Image<Gray, float32>) (fgFactor: float) : Result =
     let nbIteration = 3
 
 let kmedians (img: Image<Gray, float32>) (fgFactor: float) : Result =
     let nbIteration = 3
@@ -32,17 +32,19 @@ let kmedians (img: Image<Gray, float32>) (fgFactor: float) : Result =
     for i in 1 .. nbIteration do
         CvInvoke.Pow(img - median_bg, 2.0, d_bg)
         CvInvoke.Pow(img - median_fg, 2.0, d_fg)
     for i in 1 .. nbIteration do
         CvInvoke.Pow(img - median_bg, 2.0, d_bg)
         CvInvoke.Pow(img - median_fg, 2.0, d_fg)
-        fg <- (d_fg * fgFactor).Cmp(d_bg, CvEnum.CmpType.LessThan)
+        CvInvoke.Compare(d_fg, d_bg, fg, CvEnum.CmpType.LessThan)
 
 
-        median_fg <- MathNet.Numerics.Statistics.Statistics.Median(seq {
-            for i in 0 .. h - 1 do
-                for j in 0 .. w - 1 do
-                    if fg.Data.[i, j, 0] > 0uy then yield img.Data.[i, j, 0] |> float })
+        let bg_values = List<float>()
+        let fg_values = List<float>()
 
 
-        median_bg <- MathNet.Numerics.Statistics.Statistics.Median(seq {
-            for i in 0 .. h - 1 do
-                for j in 0 .. w - 1 do
-                    if fg.Data.[i, j, 0] = 0uy then yield img.Data.[i, j, 0] |> float })
+        for i in 0 .. h - 1 do
+            for j in 0 .. w - 1 do
+                if fg.Data.[i, j, 0] > 0uy
+                then fg_values.Add(float img.Data.[i, j, 0])
+                else bg_values.Add(float img.Data.[i, j, 0])
+
+        median_bg <- MathNet.Numerics.Statistics.Statistics.Median(bg_values)
+        median_fg <- MathNet.Numerics.Statistics.Statistics.Median(fg_values)
 
     CvInvoke.Sqrt(d_fg, d_fg)
 
 
     CvInvoke.Sqrt(d_fg, d_fg)
 
index 3f35fb7..155bfcc 100644 (file)
@@ -14,12 +14,44 @@ open Types
 
 let doAnalysis (img: Image<Bgr, byte>) (name: string) (config: Config) : Cell list =
 
 
 let doAnalysis (img: Image<Bgr, byte>) (name: string) (config: Config) : Cell list =
 
-    let imgFloat = img.Convert<Bgr, float32>()
-    use scaledImg = if config.scale = 1.0 then imgFloat else imgFloat.Resize(config.scale, CvEnum.Inter.Area)
+    use scaledImg = if config.scale = 1.0 then img else img.Resize(config.scale, CvEnum.Inter.Area)
 
 
+    use blue  = scaledImg.Item(0)
     use green = scaledImg.Item(1)
     use green = scaledImg.Item(1)
+    use red = scaledImg.Item(2)
 
 
-    use filteredGreen = (gaussianFilter green config.doGSigma1) - config.doGLowFreqPercentageReduction * (gaussianFilter green config.doGSigma2)
+
+    let greenFloat = green.Convert<Gray, float32>()
+
+    let green = gaussianFilter green 1.5
+
+    // let RBCSize = Granulometry.findRadius green (10, 100) 0.5
+
+    match config.debug with
+    | DebugOn output ->
+        let dirPath = System.IO.Path.Combine(output, name)
+        System.IO.Directory.CreateDirectory dirPath |> ignore
+        let buildFileName postfix = System.IO.Path.Combine(dirPath, name + postfix)
+
+        saveImg green (buildFileName " - green.png")
+
+        let greenMaxima = green.Copy()
+        let maxima = ImgTools.findMaxima greenMaxima
+        for m in maxima do
+            for p in m do
+                greenMaxima.Data.[p.Y, p.X, 0] <- 255uy
+
+        saveImg greenMaxima (buildFileName " - maxima.png")
+
+        logTime "areaOpen" (fun () -> ImgTools.areaOpen green 800)
+        saveImg green (buildFileName " - green opened.png")
+
+    | _ -> ()
+
+    []
+    (*
+
+    use filteredGreen = (gaussianFilter greenFloat config.doGSigma1) - config.doGLowFreqPercentageReduction * (gaussianFilter greenFloat config.doGSigma2)
 
     use sobelKernel =
         new ConvolutionKernelF(array2D [[ 1.0f; 0.0f; -1.0f ]
 
     use sobelKernel =
         new ConvolutionKernelF(array2D [[ 1.0f; 0.0f; -1.0f ]
@@ -61,7 +93,7 @@ let doAnalysis (img: Image<Bgr, byte>) (name: string) (config: Config) : Cell li
 
     let kmediansResults = logTime "Finding foreground (k-medians)" (fun () -> KMedians.kmedians filteredGreen 1.0)
 
 
     let kmediansResults = logTime "Finding foreground (k-medians)" (fun () -> KMedians.kmedians filteredGreen 1.0)
 
-    let parasites = ParasitesMarker.find green filteredGreen kmediansResults config
+    let parasites = ParasitesMarker.find greenFloat filteredGreen kmediansResults config
 
     let allEllipses, ellipses = logTime "Finding ellipses" (fun () ->
         let matchingEllipses = Ellipse.find edges xEdges yEdges config
 
     let allEllipses, ellipses = logTime "Finding ellipses" (fun () ->
         let matchingEllipses = Ellipse.find edges xEdges yEdges config
@@ -72,7 +104,11 @@ let doAnalysis (img: Image<Bgr, byte>) (name: string) (config: Config) : Cell li
     // Output pictures if debug flag is set.
     match config.debug with
     | DebugOn output ->
     // Output pictures if debug flag is set.
     match config.debug with
     | DebugOn output ->
-        let buildFileName postfix = System.IO.Path.Combine(output, name + postfix)
+        let dirPath = System.IO.Path.Combine(output, name)
+        System.IO.Directory.CreateDirectory dirPath |> ignore
+
+        let buildFileName postfix = System.IO.Path.Combine(dirPath, name + postfix)
+
         saveMat (edges * 255.0) (buildFileName " - edges.png")
 
         saveImg parasites.darkStain (buildFileName " - parasites - dark stain.png")
         saveMat (edges * 255.0) (buildFileName " - edges.png")
 
         saveImg parasites.darkStain (buildFileName " - parasites - dark stain.png")
@@ -96,6 +132,11 @@ let doAnalysis (img: Image<Bgr, byte>) (name: string) (config: Config) : Cell li
         let imgCells' = img.Copy()
         drawCells imgCells' true cells
         saveImg imgCells' (buildFileName " - cells - full.png")
         let imgCells' = img.Copy()
         drawCells imgCells' true cells
         saveImg imgCells' (buildFileName " - cells - full.png")
+
+        saveImg (normalizeAndConvert filteredGreen) (buildFileName " - filtered.png")
+        saveImg blue (buildFileName " - blue.png")
+        saveImg green (buildFileName " - green.png")
+        saveImg red (buildFileName " - red.png")
     | _ -> ()
 
     | _ -> ()
 
-    cells
+    cells*)
index da74c65..92a767a 100644 (file)
@@ -9,17 +9,23 @@ open Types
 open Utils
 
 
 open Utils
 
 
+// Do not take in account matching score below this when two ellipses are matched.
 let matchingScoreThreshold1 = 0.6
 let matchingScoreThreshold1 = 0.6
-let matchingScoreThreshold2 = 1.
+
+// All ellipsee with a score below this is removed.
+let matchingScoreThreshold2 = 1. / 50.
 
 type private EllipseScoreFlaggedKd (matchingScore: float, e: Ellipse) =
     let mutable matchingScore = matchingScore
 
 type private EllipseScoreFlaggedKd (matchingScore: float, e: Ellipse) =
     let mutable matchingScore = matchingScore
+    let perimeter = e.Perimeter
 
     member this.Ellipse = e
 
     member this.MatchingScore = matchingScore
 
     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) =
     member this.AddMatchingScore(score: float) =
-        matchingScore <- matchingScore + score
+        matchingScore <- matchingScore + score / perimeter
 
     member val Processed = false with get, set
     member val Removed = false with get, set
 
     member val Processed = false with get, set
     member val Removed = false with get, set
@@ -38,6 +44,7 @@ type MatchingEllipses (radiusMin: float) =
     member this.Ellipses : Ellipse list =
         List.ofSeq ellipses |> List.map (fun e -> e.Ellipse)
 
     member this.Ellipses : Ellipse list =
         List.ofSeq ellipses |> List.map (fun e -> e.Ellipse)
 
+    // Process all ellipses and return ellipses ordered from the best score to the worst.
     member this.PrunedEllipses : Ellipse list =
         if ellipses.Count = 0
         then
     member this.PrunedEllipses : Ellipse list =
         if ellipses.Count = 0
         then
@@ -93,7 +100,7 @@ type MatchingEllipses (radiusMin: float) =
                     for other in tree.Search window do
                         if not other.Removed && other.MatchingScore < e.MatchingScore
                         then
                     for other in tree.Search window do
                         if not other.Removed && other.MatchingScore < e.MatchingScore
                         then
-                            if e.Ellipse.Contains other.Ellipse.Cx other.Ellipse.Cy
+                            if e.Ellipse.Scale(0.8).Contains other.Ellipse.Cx other.Ellipse.Cy
                             then
                                 other.Removed <- true
             ellipses.RemoveAll(fun e -> e.Removed) |> ignore
                             then
                                 other.Removed <- true
             ellipses.RemoveAll(fun e -> e.Removed) |> ignore
index b9316f1..3e6609a 100644 (file)
@@ -28,7 +28,7 @@
     <PlatformTarget>x64</PlatformTarget>
     <DocumentationFile>bin\Debug\Parasitemia.XML</DocumentationFile>
     <Prefer32Bit>false</Prefer32Bit>
     <PlatformTarget>x64</PlatformTarget>
     <DocumentationFile>bin\Debug\Parasitemia.XML</DocumentationFile>
     <Prefer32Bit>false</Prefer32Bit>
-    <StartArguments>--folder "../../../../Images/08.10.2015/debug" --output "../../../Images/output" --debug</StartArguments>
+    <StartArguments>--folder "../../../Images/areaopen" --output "../../../Images/output" --debug</StartArguments>
   </PropertyGroup>
   <PropertyGroup Condition=" '$(Configuration)|$(Platform)' == 'Release|AnyCPU' ">
     <DebugType>pdbonly</DebugType>
   </PropertyGroup>
   <PropertyGroup Condition=" '$(Configuration)|$(Platform)' == 'Release|AnyCPU' ">
     <DebugType>pdbonly</DebugType>
@@ -40,7 +40,7 @@
     <PlatformTarget>x64</PlatformTarget>
     <DocumentationFile>bin\Release\Parasitemia.XML</DocumentationFile>
     <Prefer32Bit>false</Prefer32Bit>
     <PlatformTarget>x64</PlatformTarget>
     <DocumentationFile>bin\Release\Parasitemia.XML</DocumentationFile>
     <Prefer32Bit>false</Prefer32Bit>
-    <StartArguments>--folder "../../../../Images/08.10.2015/release" --output "../../../Images/output" --debug</StartArguments>
+    <StartArguments>--folder "../../../Images/release" --output "../../../Images/output" --debug</StartArguments>
   </PropertyGroup>
   <PropertyGroup>
     <MinimumVisualStudioVersion Condition="'$(MinimumVisualStudioVersion)' == ''">11</MinimumVisualStudioVersion>
   </PropertyGroup>
   <PropertyGroup>
     <MinimumVisualStudioVersion Condition="'$(MinimumVisualStudioVersion)' == ''">11</MinimumVisualStudioVersion>
     <Compile Include="NumericUpDown.xaml.fs" />
     <Resource Include="MainWindow.xaml" />
     <Compile Include="MainWindow.xaml.fs" />
     <Compile Include="NumericUpDown.xaml.fs" />
     <Resource Include="MainWindow.xaml" />
     <Compile Include="MainWindow.xaml.fs" />
+    <Compile Include="Heap.fs" />
     <Compile Include="Types.fs" />
     <Compile Include="Utils.fs" />
     <Compile Include="ImgTools.fs" />
     <Compile Include="Types.fs" />
     <Compile Include="Utils.fs" />
     <Compile Include="ImgTools.fs" />
+    <Compile Include="Granulometry.fs" />
     <Compile Include="Config.fs" />
     <Compile Include="KMedians.fs" />
     <Compile Include="EEOver.fs" />
     <Compile Include="Config.fs" />
     <Compile Include="KMedians.fs" />
     <Compile Include="EEOver.fs" />
index d7348cf..fdbc908 100644 (file)
@@ -65,14 +65,14 @@ let main args =
 
             // RBC size range in px at scale = 1.0.
             minRBCSize = 20.
 
             // RBC size range in px at scale = 1.0.
             minRBCSize = 20.
-            maxRBCSize = 40.
+            maxRBCSize = 42.
 
             doGSigma1 = 1.5
 
             doGSigma1 = 1.5
-            doGSigma2 = 20.
+            doGSigma2 = 30.
             doGLowFreqPercentageReduction = 0.75
 
             doGLowFreqPercentageReduction = 0.75
 
-            factorNbPick = 2.0
-            factorWindowSize = 1.6
+            factorNbPick = 1.0
+            factorWindowSize = 2.0
 
             darkStainLevel = 0.4 // Lower -> more sensitive.
 
 
             darkStainLevel = 0.4 // Lower -> more sensitive.
 
@@ -88,18 +88,18 @@ let main args =
 
             MaxDarkStainRatio = 0.1
 
 
             MaxDarkStainRatio = 0.1
 
-            minimumCellArea = 1200. * scale ** 2. |> int
+            minimumCellArea = 1000. * scale ** 2. |> int
             maxOffcenter = 0.5 }
 
         match mode with
         | CmdLine (input, output) ->
             maxOffcenter = 0.5 }
 
         match mode with
         | CmdLine (input, output) ->
-            let config = { config with debug = DebugOn output }
+            let config = if debug then { config with debug = DebugOn output } else config
 
             Directory.CreateDirectory output |> ignore
 
             use logFile = new StreamWriter(new FileStream(Path.Combine(output, "log.txt"), FileMode.Append, FileAccess.Write))
             Utils.log <- (fun m -> logFile.WriteLine(m))
 
             Directory.CreateDirectory output |> ignore
 
             use logFile = new StreamWriter(new FileStream(Path.Combine(output, "log.txt"), FileMode.Append, FileAccess.Write))
             Utils.log <- (fun m -> logFile.WriteLine(m))
-            Utils.log (sprintf "=== New run : %A ===" DateTime.Now)
+            Utils.log (sprintf "=== New run : %A %A ===" DateTime.Now  (if debug then "[DEBUG]" else "[RELEASE]"))
 
             let files = match input with
                         | File file -> [ file ]
 
             let files = match input with
                         | File file -> [ file ]
@@ -111,7 +111,8 @@ let main args =
                 try
                     use img = new Image<Bgr, byte>(file)
                     Utils.log (sprintf "== File: %A" file)
                 try
                     use img = new Image<Bgr, byte>(file)
                     Utils.log (sprintf "== File: %A" file)
-                    let cells = Utils.logTime "Whole analyze" (fun () -> ImageAnalysis.doAnalysis img (FileInfo(file).Name) config)
+
+                    let cells = Utils.logTime "Whole analyze" (fun () -> ImageAnalysis.doAnalysis img (Path.GetFileNameWithoutExtension(FileInfo(file).Name)) config)
                     let total, infected = Utils.countCells cells
                     fprintf resultFile "File: %s %d %d %.2f\n" file total infected (100. * (float infected) / (float total))
                 with
                     let total, infected = Utils.countCells cells
                     fprintf resultFile "File: %s %d %d %.2f\n" file total infected (100. * (float infected) / (float total))
                 with
index d7bacbd..a551142 100644 (file)
@@ -30,6 +30,13 @@ type Ellipse (cx: float, cy: float, a: float, b: float, alpha: float) =
         this.CutAVericalLine 0.0 || this.CutAVericalLine width ||
         this.CutAnHorizontalLine 0.0 || this.CutAnHorizontalLine height
 
         this.CutAVericalLine 0.0 || this.CutAVericalLine width ||
         this.CutAnHorizontalLine 0.0 || this.CutAnHorizontalLine height
 
+    member this.Scale (factor: float) =
+        Ellipse(this.Cx, this.Cy, this.A * factor, this.B * factor, alpha)
+
+    // Approximation of Ramanujan.
+    member this.Perimeter =
+        Math.PI * (3.0 * (this.A + this.B) - sqrt ((3.0 * this.A + this.B) * (this.A + 3.0 * this.B)))
+
 
 type CellClass = HealthyRBC | InfectedRBC | Peculiar
 
 
 type CellClass = HealthyRBC | InfectedRBC | Peculiar