* Remove the 'DoG' filter.
authorGreg Burri <greg.burri@gmail.com>
Sun, 20 Dec 2015 13:06:49 +0000 (14:06 +0100)
committerGreg Burri <greg.burri@gmail.com>
Sun, 20 Dec 2015 13:06:49 +0000 (14:06 +0100)
* Fix the are opening.

Parasitemia/Parasitemia/Config.fs
Parasitemia/Parasitemia/EEOver.fs
Parasitemia/Parasitemia/Granulometry.fs
Parasitemia/Parasitemia/ImgTools.fs
Parasitemia/Parasitemia/MainAnalysis.fs
Parasitemia/Parasitemia/MatchingEllipses.fs
Parasitemia/Parasitemia/Program.fs
Parasitemia/Parasitemia/Types.fs

index 464fd7e..48a6334 100644 (file)
@@ -12,9 +12,7 @@ type Config = {
     minRBCSize: float
     maxRBCSize: float
 
-    doGSigma1: float
-    doGSigma2: float
-    doGLowFreqPercentageReduction: float
+    preFilterSigma: float
 
     // Ellipse.
     factorNbPick: float
index 52b5963..1a6bb5f 100644 (file)
@@ -7,7 +7,7 @@ let private EPS = 1.0e-5
 let inline private ellipse2tr (x: float) (y: float) (aa: float) (bb: float) (cc: float) (dd: float) (ee: float) (ff: float) : float =
     aa * x * x + bb * x * y + cc * y * y + dd * x + ee * y + ff
 
-let private nointpts (a1: float) (b1: float) (a2: float) (b2: float) (h1: float) (k1: float) (h2: float) (k2: float) (phi_1: float) (phi_2: float) (h2_tr: float) (k2_tr: float) (aa: float) (bb: float) (cc: float) (dd: float) (ee: float) (ff: float) = 
+let private nointpts (a1: float) (b1: float) (a2: float) (b2: float) (h1: float) (k1: float) (h2: float) (k2: float) (phi_1: float) (phi_2: float) (h2_tr: float) (k2_tr: float) (aa: float) (bb: float) (cc: float) (dd: float) (ee: float) (ff: float) =
     let a1b1 = a1 * b1
     let a2b2 = a2 * b2
     let area_1 = Math.PI * a1b1
@@ -34,13 +34,13 @@ let private nointpts (a1: float) (b1: float) (a2: float) (b2: float) (h1: float)
 type private PointType = TANGENT_POINT | INTERSECTION_POINT
 
 let private istanpt (x: float) (y: float) (a1: float) (b1: float) (aa: float) (bb: float) (cc: float) (dd: float) (ee: float) (ff: float) : PointType =
-    let x = 
+    let x =
         if abs x > a1
         then
             if x < 0.0 then -a1 else a1
         else x
 
-    let theta = 
+    let theta =
         if y < 0.0
         then 2.0 * Math.PI - acos (x / a1)
         else acos (x / a1)
@@ -105,7 +105,7 @@ let private twointpts (x: float[]) (y: float[]) (a1: float) (b1: float) (phi_1:
     if theta1 > theta2
     then
         theta1 <- theta1 - 2.0 * Math.PI
-    
+
     let trsign = if (theta2 - theta1) > Math.PI then 1.0 else -1.0
 
     let mutable area1 = 0.5 * (a1 * b1 * (theta2 - theta1) + trsign * abs (x.[0] * y.[1] - x.[1] * y.[0]))
@@ -182,7 +182,7 @@ let private twointpts (x: float[]) (y: float[]) (a1: float) (b1: float) (phi_1:
     area1 + area2
 
 
-let private threeintpts (xint: float[]) (yint: float[]) (a1: float) (b1: float) (phi_1: float) (a2: float) (b2: float) (h2_tr: float) (k2_tr: float) (phi_2: float) (aa: float) (bb: float) (cc: float) (dd: float) (ee: float) (ff: float) : float = 
+let private threeintpts (xint: float[]) (yint: float[]) (a1: float) (b1: float) (phi_1: float) (a2: float) (b2: float) (h2_tr: float) (k2_tr: float) (phi_2: float) (aa: float) (bb: float) (cc: float) (dd: float) (ee: float) (ff: float) : float =
     let mutable tanpts = 0
     let mutable tanindex = 0
     for i in 0..2 do
@@ -205,7 +205,7 @@ let private threeintpts (xint: float[]) (yint: float[]) (a1: float) (b1: float)
         | 1 ->
             xint.[1] <- xint.[2]
             yint.[1] <- yint.[2]
-        | _ -> 
+        | _ ->
             ()
         twointpts xint yint a1 b1 phi_1 a2 b2 h2_tr k2_tr phi_2 aa bb cc dd ee ff
 
@@ -244,8 +244,8 @@ let private fourintpts (xint: float[]) (yint: float[]) (a1: float) (b1: float) (
             else
                 theta.[k+1] <- theta.[k]
                 xint.[k+1] <- xint.[k]
-                yint.[k+1] <- yint.[k]            
-                k <- k - 1                
+                yint.[k+1] <- yint.[k]
+                k <- k - 1
                 k2 <- k + 1
 
         theta.[k2] <- tmp0
@@ -258,7 +258,7 @@ let private fourintpts (xint: float[]) (yint: float[]) (a1: float) (b1: float) (
     for k in 0..3 do
         printf "k=%d:  Theta = %f, xint=%f, yint=%f\n" k theta.[k] xint.[k] yint.[k]
 #endif
-                
+
     let area1 = 0.5 * abs ((xint.[2] - xint.[0]) * (yint.[3] - yint.[1]) - (xint.[3] - xint.[1]) * (yint.[2] - yint.[0]))
 
     let cosphi = cos (phi_1 - phi_2)
@@ -275,7 +275,7 @@ let private fourintpts (xint: float[]) (yint: float[]) (a1: float) (b1: float) (
         if abs xint_tr.[i] > a2
         then
             xint_tr.[i] <- if xint_tr.[i] < 0.0 then -a2 else a2
-        
+
         theta_tr.[i] <- if yint_tr.[i] < 0.0 then 2.0 * Math.PI - acos (xint_tr.[i] / a2) else acos (xint_tr.[i] / a2)
 
     let xmid = a1 * cos ((theta.[0] + theta.[1]) / 2.0)
@@ -298,9 +298,9 @@ let private fourintpts (xint: float[]) (yint: float[]) (a1: float) (b1: float) (
         area3 <- 0.5 * (a1b1 * (theta.[0] - (theta.[3] - 2.0 * Math.PI)) - abs (xint.[3] * yint.[0] - xint.[0] * yint.[3]))
         area4 <- 0.5 * (a2b2 * (theta_tr.[1] - theta_tr.[0]) - abs (xint_tr.[0] * yint_tr.[1] - xint_tr.[1] * yint_tr.[0]))
         area5 <- 0.5 * (a2b2 * (theta_tr.[3] - theta_tr.[2]) - abs (xint_tr.[2] * yint_tr.[3] - xint_tr.[3] * yint_tr.[2]))
-        
+
     if area5 < 0.0
-    then 
+    then
 #if DEBUG_LOG
         printf "\n\t\t-------------> area5 is negativ (%f). Add: pi*A2*B2=%f <------------\n" area5 area_2
 #endif
@@ -321,7 +321,7 @@ let private fourintpts (xint: float[]) (yint: float[]) (a1: float) (b1: float) (
         area3 <- area3 + area_1
 
     if area2 < 0.0
-    then    
+    then
 #if DEBUG_LOG
         printf "\n\t\t-------------> area2 is negativ (%f). Add: pi*A2*B2=%f <------------\n" area2 area_1
 #endif
@@ -334,7 +334,7 @@ let private fourintpts (xint: float[]) (yint: float[]) (a1: float) (b1: float) (
     area1 + area2 + area3 + area4 + area5
 
 
-let private quadroots (p: float[]) (r: float[,]) = 
+let private quadroots (p: float[]) (r: float[,]) =
     let mutable b = -p.[1] / (2.0 * p.[0])
     let c = p.[2] / p.[0]
     let mutable d = b * b - c
@@ -342,7 +342,7 @@ let private quadroots (p: float[]) (r: float[,]) =
     if d >= 0.0
     then
         if b > 0.0
-        then 
+        then
             b <- sqrt d + b
             r.[1, 2] <- b
         else
@@ -370,7 +370,7 @@ let private cubicroots (p: float[]) (r: float[,]) =
     let mutable c = t * t * t
     let mutable d = b * b - c
 
-    if d >= 0.0 
+    if d >= 0.0
     then
         d <- ((sqrt d) + (abs b)) ** (1.0 / 3.0)
         if d <> 0.0
@@ -427,7 +427,7 @@ let private cubicroots (p: float[]) (r: float[,]) =
         r.[1, 1] <- t
         for k in 1..3 do
             r.[2, k] <- 0.0
-        
+
 
 let private biquadroots (p: float[]) (r: float[,]) =
     if p.[0] <> 1.0
@@ -445,7 +445,7 @@ let private biquadroots (p: float[]) (r: float[,]) =
     a <- a - d
 
     let quadExecuted = ref false
-    let quad () = 
+    let quad () =
         if not !quadExecuted
         then
             p.[2] <- !c / !b
@@ -482,7 +482,7 @@ let private biquadroots (p: float[]) (r: float[,]) =
                     p.[1] <- -(sqrt d)
                 b := 0.5 * (a + !b / p.[1])
                 quad ()
-                k <- 4            
+                k <- 4
             k <- k + 1
 
     if not !quadExecuted && p.[2] < 0.0
@@ -517,17 +517,17 @@ let EEOverlapArea (e1: Types.Ellipse) (e2: Types.Ellipse) : (float * float[] * f
     let h1, k1, a1, b1, phi_1 = e1.Cx, e1.Cy, e1.A, e1.B, e1.Alpha
     let h2, k2, a2, b2, phi_2 = e2.Cx, e2.Cy, e2.A, e2.B, e2.Alpha
 
-    if a1 <= EPS || b1 <= EPS || a2 <= EPS || b2 <= EPS 
+    if a1 <= EPS || b1 <= EPS || a2 <= EPS || b2 <= EPS
     then
         None
     else
-        let phi_1 = phi_1 % Math.PI
-        let phi_2 = phi_2 % Math.PI
+        let phi_1 = phi_1 % Math.PI //(if phi_1 > Math.PI / 2.0 then phi_1 - Math.PI else phi_1) % Math.PI
+        let phi_2 = phi_2 % Math.PI //(if phi_2 > Math.PI / 2.0 then phi_2 - Math.PI else phi_2) % Math.PI
         let h2_tr, k2_tr, phi_2r =
             let cosphi = cos phi_1
             let sinphi = sin phi_1
             (h2 - h1) * cosphi + (k2 - k1) * sinphi, (h1 - h2) * sinphi + (k2 - k1) * cosphi, (phi_2 - phi_1) % (2.0 * Math.PI)
-        
+
 #if DEBUG_LOG
         printf "H2_TR=%f, K2_TR=%f, PHI_2R=%f\n" h2_tr k2_tr phi_2r
 #endif
@@ -562,7 +562,7 @@ let EEOverlapArea (e1: Types.Ellipse) (e2: Types.Ellipse) : (float * float[] * f
 #if DEBUG_LOG
         for i in 0..4 do
             printf "cy[%d]=%f\n" i cy.[i]
-#endif 
+#endif
 
         let py = Array.zeroCreate<float> 5
         let r = Array2D.zeroCreate<float> 3 5
@@ -587,7 +587,7 @@ let EEOverlapArea (e1: Types.Ellipse) (e2: Types.Ellipse) : (float * float[] * f
                 py.[0] <- 1.0
                 cubicroots py r
                 3
-        
+
             elif abs cy.[2] > EPS
             then
                 for i in 0..1 do
@@ -652,7 +652,7 @@ let EEOverlapArea (e1: Types.Ellipse) (e2: Types.Ellipse) : (float * float[] * f
 #if DEBUG_LOG
                 printf "\tx1=%f, y1=%f, A=%f. B=%f ---> ellipse2tr(x1)= %f\n" x1 ychk.[i] a1 b1 (ellipse2tr x1 ychk.[i] aa bb cc dd ee ff)
                 printf "\tx2=%f, y1=%f, A=%f. B=%f ---> ellipse2tr(x2) %f\n" x2 ychk.[i] a1 b1 (ellipse2tr x2 ychk.[i] aa bb cc dd ee ff)
-#endif         
+#endif
 
                 if abs (ellipse2tr x1 ychk.[i] aa bb cc dd ee ff) < EPS
                 then
@@ -661,7 +661,7 @@ let EEOverlapArea (e1: Types.Ellipse) (e2: Types.Ellipse) : (float * float[] * f
                     printf "first if x1. acc nintps=%d\n" nintpts
 #endif
                     if nintpts > 4
-                    then 
+                    then
                         returnValue <- -1.0
                     else
                         xint.[nintpts-1] <- x1
@@ -677,7 +677,7 @@ let EEOverlapArea (e1: Types.Ellipse) (e2: Types.Ellipse) : (float * float[] * f
                     printf "first if x2. nintps=%d, Dx=%f (eps2=%f) \n" nintpts (abs (x2 - x1)) EPS
 #endif
                     if nintpts > 4
-                    then 
+                    then
                         returnValue <- -1.0
                     else
                         xint.[nintpts-1] <- x2
@@ -692,8 +692,8 @@ let EEOverlapArea (e1: Types.Ellipse) (e2: Types.Ellipse) : (float * float[] * f
                 printf "i=%d, multiple roots: %f  <--------> %f. continue\n" i ychk.[i] ychk.[i-1]
 #endif
             i <- i + 1
-                    
-        
+
+
         if returnValue = -1.0
         then
             None
@@ -702,13 +702,13 @@ let EEOverlapArea (e1: Types.Ellipse) (e2: Types.Ellipse) : (float * float[] * f
                 match nintpts with
                 | 0 | 1 -> nointpts a1 b1 a2 b2 h1 k1 h2 k2 phi_1 phi_2 h2_tr k2_tr aa bb cc dd ee ff
                 | 2 -> match istanpt xint.[0] yint.[0] a1 b1 aa bb cc dd ee ff with
-                       | TANGENT_POINT -> 
+                       | TANGENT_POINT ->
 #if DEBUG_LOG
                             printf "one point is tangent\n"
 #endif
                             nointpts a1 b1 a2 b2 h1 k1 h2 k2 phi_1 phi_2 h2_tr k2_tr aa bb cc dd ee ff
 
-                       | INTERSECTION_POINT -> 
+                       | INTERSECTION_POINT ->
 #if DEBUG_LOG
                             printf "check twointpts\n"
 #endif
@@ -718,7 +718,7 @@ let EEOverlapArea (e1: Types.Ellipse) (e2: Types.Ellipse) : (float * float[] * f
                 | _ -> -1.0
             if nintpts = 0
             then Some (area, [||], [||])
-            else 
+            else
                 let xTransform = Array.zeroCreate nintpts
                 let yTransform = Array.zeroCreate nintpts
                 for i in 0 .. (nintpts - 1) do
index a9a102a..db299dc 100644 (file)
@@ -10,7 +10,7 @@ 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 =
+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
@@ -24,15 +24,15 @@ let findRadius (img: Image<Gray, byte>)  (range: int * int) (scale: float) : int
         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
+        let n = closed.GetSum().Intensity
 
         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)
+    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
+    float (max + r1') / scale |> roundInt
 
 
index e9cbbe9..24c6998 100644 (file)
@@ -26,7 +26,13 @@ let gaussianFilter (img : Image<'TColor, 'TDepth>) (standardDeviation : float) :
     img.SmoothGaussian(size, size, standardDeviation, standardDeviation)
 
 
-let findMaxima (img: Image<Gray, byte>) : IEnumerable<HashSet<Point>> =
+type Points = HashSet<Point>
+
+let drawPoints (img: Image<Gray, byte>) (points: Points) (intensity: byte) =
+    for p in points do
+        img.Data.[p.Y, p.X, 0] <- intensity
+
+let findMaxima (img: Image<Gray, byte>) : IEnumerable<Points> =
     use suppress = new Image<Gray, byte>(img.Size)
     let w = img.Width
     let h = img.Height
@@ -90,159 +96,189 @@ let findMaxima (img: Image<Gray, byte>) : IEnumerable<HashSet<Point>> =
             if maxima.Count > 0
             then result.AddRange(maxima)
 
-    result.Select(fun l -> HashSet<Point>(l))
+    result.Select(fun l -> Points(l))
 
 
 type PriorityQueue () =
-    let q = List<HashSet<Point>>() // TODO: Check performance with an HasSet
+    let size = 256
+    let q: Points[] = Array.init size (fun i -> Points()) // TODO: Check performance with an HasSet
     let mutable highest = -1 // Value of the first elements of 'q'.
+    let mutable lowest = size
 
     member this.Next () : byte * Point =
         if this.IsEmpty
         then
             invalidOp "Queue is empty"
         else
-            let l = q.[0]
+            let l = q.[highest]
             let next = l.First()
             l.Remove(next) |> ignore
             let value = byte highest
+
             if l.Count = 0
             then
-                q.RemoveAt(0)
                 highest <- highest - 1
-                while q.Count > 0 && q.[0] = null do
-                    q.RemoveAt(0)
+                while highest > lowest && q.[highest].Count = 0 do
                     highest <- highest - 1
+                if highest = lowest
+                then
+                    highest <- -1
+                    lowest <- size
+
             value, next
 
     member this.Max =
         highest |> byte
 
-    (*member this.UnionWith (other: PriorityQueue) =
-        while not other.IsEmpty do
-            let p, v = other.Next
-            this.Add p v*)
-
     member this.Add (value: byte) (p: Point) =
         let vi = int value
 
-        if this.IsEmpty
-        then
-            highest <- int value
-            q.Insert(0, null)
-        elif vi > highest
+        if vi > highest
         then
-            for i in highest .. vi - 1  do
-                q.Insert(0, null)
             highest <- vi
-        elif highest - vi >= q.Count
+        if vi <= lowest
         then
-            for i in 0 .. highest - vi - q.Count do
-                q.Add(null)
+            lowest <- vi - 1
 
-        let pos = highest - vi
-        if q.[pos] = null
+        q.[vi].Add(p) |> ignore
+
+    member this.Remove (value: byte) (p: Point) =
+        let vi = int value
+        if q.[vi].Remove(p) && q.[vi].Count = 0
         then
-            q.[pos] <- HashSet<Point>([p])
-        else
-            q.[pos].Add(p) |> ignore
+            if vi = highest
+            then
+                highest <- highest - 1
+                while highest > lowest && q.[highest].Count = 0 do
+                    highest <- highest - 1
+            elif vi - 1 = lowest
+            then
+                lowest <- lowest + 1
+                while lowest < highest && q.[lowest + 1].Count = 0 do
+                    lowest <- lowest + 1
 
+            if highest = lowest // The queue is now empty.
+            then
+                highest <- -1
+                lowest <- size
 
     member this.IsEmpty =
-        q.Count = 0
+        highest = -1
 
     member this.Clear () =
-        while highest >= 0 do
-            q.[highest] <- null
+        while highest > lowest  do
+            q.[highest].Clear()
             highest <- highest - 1
+        highest <- -1
+        lowest <- size
 
 
+type AreaState =
+    | Removed = 1
+    | Unprocessed = 2
+    | Validated = 3
 
-type MaximaState =  Uncertain | Validated | TooBig
-type Maxima = {
-    elements : HashSet<Point>
-    mutable intensity: byte option
-    mutable state: MaximaState }
-
+[<AllowNullLiteral>]
+type Area (elements: Points) =
+    member this.Elements = elements
+    member val Intensity = None with get, set
+    member val State = AreaState.Unprocessed with get, set
 
 let areaOpen (img: Image<Gray, byte>) (area: int) =
     let w = img.Width
     let h = img.Height
+    let imgData = img.Data
+
+    let areas = List<Area>(findMaxima img |> Seq.map Area)
+
+    let pixels: Area[,] = Array2D.create h w null
+    for m in areas do
+        for e in m.Elements do
+            pixels.[e.Y, e.X] <- m
 
-    let maxima = findMaxima img |> Seq.map (fun m -> { elements = m; intensity = None; state = Uncertain }) |> List.ofSeq
-    let toValidated = Stack<Maxima>(maxima)
+    let queue = PriorityQueue()
 
-    while toValidated.Count > 0 do
-        let m = toValidated.Pop()
-        if m.elements.Count <= area
+    let addEdgeToQueue (elements: Points) =
+        for p in elements do
+            for i in -1 .. 1 do
+                for j in -1 .. 1 do
+                    if i <> 0 || j <> 0
+                    then
+                        let ni = i + p.Y
+                        let nj = j + p.X
+                        let p' = Point(nj, ni)
+                        if ni >= 0 && ni < h && nj >= 0 && nj < w && not (elements.Contains(p'))
+                        then
+                            queue.Add (imgData.[ni, nj, 0]) p'
+
+    // Reverse order is quicker.
+    for i in areas.Count - 1 .. -1 .. 0 do
+        let m = areas.[i]
+        if m.Elements.Count <= area && m.State <> AreaState.Removed
         then
-            let queue =
-                let q = PriorityQueue()
-                let firstElements = HashSet<Point>()
-                for p in m.elements do
-                    for i in -1 .. 1 do
-                        for j in -1 .. 1 do
-                            if i <> 0 || j <> 0
-                            then
-                                let ni = i + p.Y
-                                let nj = j + p.X
-                                let p' = Point(nj, ni)
-                                if ni >= 0 && ni < h && nj >= 0 && nj < w && not (m.elements.Contains(p')) && not (firstElements.Contains(p'))
-                                then
-                                    firstElements.Add(p') |> ignore
-                                    q.Add (img.Data.[ni, nj, 0]) p'
-                q
+            queue.Clear()
+            addEdgeToQueue m.Elements
 
             let mutable intensity = queue.Max
-            let nextElements = HashSet<Point>()
+            let nextElements = Points()
 
             let mutable stop = false
             while not stop do
                 let intensity', p = queue.Next ()
+                let mutable merged = false
 
                 if intensity' = intensity // The intensity doesn't change.
                 then
-                    if m.elements.Count + nextElements.Count + 1 > area
+                    if m.Elements.Count + nextElements.Count + 1 > area
                     then
-                        m.state <- Validated
-                        m.intensity <- Some intensity
+                        m.State <- AreaState.Validated
+                        m.Intensity <- Some intensity
                         stop <- true
                     else
                         nextElements.Add(p) |> ignore
+
                 elif intensity' < intensity
                 then
-                    m.elements.UnionWith(nextElements)
-                    if m.elements.Count = area
+                    m.Elements.UnionWith(nextElements)
+                    for e in nextElements do
+                        pixels.[e.Y, e.X] <- m
+
+                    if m.Elements.Count = area
                     then
-                        m.state <- Validated
-                        m.intensity <- Some (intensity')
+                        m.State <- AreaState.Validated
+                        m.Intensity <- Some (intensity')
                         stop <- true
                     else
                         intensity <- intensity'
                         nextElements.Clear()
                         nextElements.Add(p) |> ignore
+
                 else // i' > i
-                    seq {
-                        for m' in maxima do
-                            if m' <> m && m'.elements.Contains(p) then
-                                if m'.elements.Count + m.elements.Count <= area
-                                then
-                                    m'.state <- Uncertain
-                                    m'.elements.UnionWith(m.elements)
-                                    if not <| toValidated.Contains m' // FIXME: Maybe use state instead of scanning the whole list.
-                                    then
-                                        toValidated.Push(m')
-                                    stop <- true
-                                yield false
-                    } |> Seq.forall id |> ignore
+                    let m' = pixels.[p.Y, p.X]
+                    if m' <> null
+                    then
+                        if m'.Elements.Count + m.Elements.Count <= area
+                        then
+                            m'.State <- AreaState.Removed
+                            for e in m'.Elements do
+                                pixels.[e.Y, e.X] <- m
+                                queue.Remove imgData.[e.Y, e.X, 0] e
+                            addEdgeToQueue m'.Elements
+                            m.Elements.UnionWith(m'.Elements)
+                            let intensityMax = queue.Max
+                            if intensityMax <> intensity
+                            then
+                                intensity <- intensityMax
+                                nextElements.Clear()
+                            merged <- true
 
-                    if not stop
+                    if not merged
                     then
-                        m.state <- Validated
-                        m.intensity <- Some (intensity)
+                        m.State <- AreaState.Validated
+                        m.Intensity <- Some (intensity)
                         stop <- true
 
-                if not stop
+                if not stop && not merged
                 then
                     for i in -1 .. 1 do
                         for j in -1 .. 1 do
@@ -253,29 +289,29 @@ let areaOpen (img: Image<Gray, byte>) (area: int) =
                                 let p' = Point(nj, ni)
                                 if ni < 0 || ni >= h || nj < 0 || nj >= w
                                 then
-                                    m.state <- Validated
-                                    m.intensity <- Some (intensity)
+                                    m.State <- AreaState.Validated
+                                    m.Intensity <- Some (intensity)
                                     stop <- true
-                                elif not (m.elements.Contains(p')) && not (nextElements.Contains(p'))
+                                elif not (m.Elements.Contains(p')) && not (nextElements.Contains(p'))
                                 then
-                                    queue.Add (img.Data.[ni, nj, 0]) p'
+                                    queue.Add (imgData.[ni, nj, 0]) p'
 
                 if queue.IsEmpty
                 then
-                    if m.elements.Count + nextElements.Count <= area
+                    if m.Elements.Count + nextElements.Count <= area
                     then
-                        m.state <- Validated
-                        m.intensity <- Some intensity'
-                        m.elements.UnionWith(nextElements)
+                        m.State <- AreaState.Validated
+                        m.Intensity <- Some intensity'
+                        m.Elements.UnionWith(nextElements)
                     stop <- true
 
-    for m in maxima do
-        if m.state = Validated
+    for m in areas do
+        if m.State = AreaState.Validated
         then
-            match m.intensity with
+            match m.Intensity with
             | Some i ->
-                for p in m.elements do
-                    img.Data.[p.Y, p.X, 0] <- i
+                for p in m.Elements do
+                    imgData.[p.Y, p.X, 0] <- i
             | _ -> ()
     ()
 
@@ -386,7 +422,7 @@ let connectedComponents (img: Image<Gray, byte>) (startPoints: List<Point>) : Li
     let w = img.Width
     let h = img.Height
 
-    let pointChecked = HashSet<Point>()
+    let pointChecked = Points()
     let pointToCheck = List<Point>(startPoints);
 
     let data = img.Data
index 155bfcc..39383fc 100644 (file)
@@ -20,46 +20,25 @@ let doAnalysis (img: Image<Bgr, byte>) (name: string) (config: Config) : Cell li
     use green = scaledImg.Item(1)
     use red = scaledImg.Item(2)
 
-
     let greenFloat = green.Convert<Gray, float32>()
 
-    let green = gaussianFilter green 1.5
-
-    // let RBCSize = Granulometry.findRadius green (10, 100) 0.5
+    let filteredGreen = gaussianFilter green config.preFilterSigma
+    logTime "areaOpen" (fun () -> ImgTools.areaOpen filteredGreen 2000)
 
-    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)
+    let RBCSize = Granulometry.findRadius filteredGreen (10, 100) 0.5 |> float
+    let RBCSizeDelta = RBCSize / 3.0
 
-        saveImg green (buildFileName " - green.png")
+    let config = { config with minRBCSize = RBCSize - RBCSizeDelta ; maxRBCSize = RBCSize + RBCSizeDelta }
 
-        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)
+    let filteredGreenFloat = filteredGreen.Convert<Gray, float32>() // Is it neccessary?
 
     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 = filteredGreen.Convolution(sobelKernel).Convert<Gray, float>()
-    use yEdges = filteredGreen.Convolution(sobelKernel.Transpose()).Convert<Gray, float>()
+    use xEdges = filteredGreenFloat.Convolution(sobelKernel).Convert<Gray, float>()
+    use yEdges = filteredGreenFloat.Convolution(sobelKernel.Transpose()).Convert<Gray, float>()
 
     let xEdgesData = xEdges.Data
     let yEdgesData = yEdges.Data
@@ -91,9 +70,9 @@ 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 filteredGreen 1.0)
+    let kmediansResults = logTime "Finding foreground (k-medians)" (fun () -> KMedians.kmedians filteredGreenFloat 1.0)
 
-    let parasites = ParasitesMarker.find greenFloat filteredGreen kmediansResults config
+    let parasites = ParasitesMarker.find greenFloat filteredGreenFloat kmediansResults config
 
     let allEllipses, ellipses = logTime "Finding ellipses" (fun () ->
         let matchingEllipses = Ellipse.find edges xEdges yEdges config
@@ -133,10 +112,15 @@ let doAnalysis (img: Image<Bgr, byte>) (name: string) (config: Config) : Cell li
         drawCells imgCells' true cells
         saveImg imgCells' (buildFileName " - cells - full.png")
 
-        saveImg (normalizeAndConvert filteredGreen) (buildFileName " - filtered.png")
+        let filteredGreenMaxima = gaussianFilter green config.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 green (buildFileName " - green.png")
         saveImg red (buildFileName " - red.png")
     | _ -> ()
 
-    cells*)
+    cells
index 92a767a..99e3cb5 100644 (file)
@@ -13,7 +13,7 @@ open Utils
 let matchingScoreThreshold1 = 0.6
 
 // All ellipsee with a score below this is removed.
-let matchingScoreThreshold2 = 1. / 50.
+let matchingScoreThreshold2 = 1. / 100.
 
 type private EllipseScoreFlaggedKd (matchingScore: float, e: Ellipse) =
     let mutable matchingScore = matchingScore
index fdbc908..86de435 100644 (file)
@@ -67,11 +67,9 @@ let main args =
             minRBCSize = 20.
             maxRBCSize = 42.
 
-            doGSigma1 = 1.5
-            doGSigma2 = 30.
-            doGLowFreqPercentageReduction = 0.75
+            preFilterSigma = 1.5
 
-            factorNbPick = 1.0
+            factorNbPick = 1.0 // 1.0
             factorWindowSize = 2.0
 
             darkStainLevel = 0.4 // Lower -> more sensitive.
@@ -81,7 +79,7 @@ let main args =
             stainSpreadRequired = 3.0
 
             infectionSigma = 2.2
-            infectionLevel = 0.85
+            infectionLevel = 0.87
             infectionPixelsRequired = 1
 
             percentageOfFgValidCell = 0.4
index a551142..634a7b7 100644 (file)
@@ -37,6 +37,9 @@ type Ellipse (cx: float, cy: float, a: float, b: float, alpha: float) =
     member this.Perimeter =
         Math.PI * (3.0 * (this.A + this.B) - sqrt ((3.0 * this.A + this.B) * (this.A + 3.0 * this.B)))
 
+    override this.ToString () =
+        sprintf "(cx: %A, cy: %A, a: %A, b: %A, alpha: %A)" this.Cx this.Cy this.A this.B this.Alpha
+
 
 type CellClass = HealthyRBC | InfectedRBC | Peculiar