Cleaning, micro-optimizations.
authorGreg Burri <greg.burri@gmail.com>
Mon, 4 Jan 2016 12:09:54 +0000 (13:09 +0100)
committerGreg Burri <greg.burri@gmail.com>
Mon, 4 Jan 2016 12:09:54 +0000 (13:09 +0100)
Parasitemia/Parasitemia/Classifier.fs
Parasitemia/Parasitemia/EEOver.fs
Parasitemia/Parasitemia/Heap.fs
Parasitemia/Parasitemia/ImgTools.fs
Parasitemia/Parasitemia/MainAnalysis.fs
Parasitemia/Parasitemia/MatchingEllipses.fs
Parasitemia/Parasitemia/Program.fs

index f17c0b9..0d21375 100644 (file)
@@ -93,7 +93,12 @@ let findCells (ellipses: Ellipse list) (parasites: ParasitesMarker.Result) (img:
         // We reverse the list to get the lower score ellipses first.
         let ellipsesWithNeigbors = ellipses |> List.map (fun e -> e, neighbors e) |> List.rev
 
-        // 2) Remove ellipses with a high standard deviation (high contrast).
+
+        // 2) Remove ellipses touching the edges.
+        for e in ellipses do
+            if e.isOutside w_f h_f then e.Removed <- true
+
+        // 3) Remove ellipses with a high standard deviation (high contrast).
 
         // CvInvoke.Normalize(img, img, 0.0, 255.0, CvEnum.NormType.MinMax) // Not needed.
 
@@ -103,21 +108,20 @@ let findCells (ellipses: Ellipse list) (parasites: ParasitesMarker.Result) (img:
                     yield float img.Data.[y, x, 0] })
 
         for e in ellipses do
-            let minX, minY, maxX, maxY = ellipseWindow e
-
-            let stdDeviation = MathNet.Numerics.Statistics.Statistics.StandardDeviation (seq {
-                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
-                        if e.Contains (float x) (float y)
-                        then
-                            yield float img.Data.[y, x, 0] })
+            if not e.Removed
+            then
+                let shrinkedE = e.Scale 0.9
+                let minX, minY, maxX, maxY = ellipseWindow shrinkedE
 
-            if stdDeviation > globalStdDeviation * config.Parameters.standardDeviationMaxRatio then
-                e.Removed <- true
+                let stdDeviation = MathNet.Numerics.Statistics.Statistics.StandardDeviation (seq {
+                    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
+                            if shrinkedE.Contains (float x) (float y)
+                            then
+                                yield float img.Data.[y, x, 0] })
 
-        // 3) Remove ellipses touching the edges.
-        for e in ellipses do
-            if e.isOutside w_f h_f then e.Removed <- true
+                if stdDeviation > globalStdDeviation * config.Parameters.standardDeviationMaxRatio then
+                    e.Removed <- true
 
         // 4) Remove ellipses with little area.
         let minArea = config.RBCMinArea
index ef21051..528632a 100644 (file)
@@ -437,33 +437,33 @@ let private biquadroots (p: float[]) (r: float[,]) =
             p.[k] <- p.[k] / p.[0]
         p.[0] <- 1.0
     let e = 0.25 * p.[1]
-    let b = ref (2.0 * e)
-    let c = ref (!b ** 2.0)
-    let mutable d = 0.75 * !c
-    b := p.[3] + !b *(!c - p.[2])
+    let mutable b = 2.0 * e
+    let mutable c = b ** 2.0
+    let mutable d = 0.75 * c
+    b <- p.[3] + b *(c - p.[2])
     let mutable a = p.[2] - d
-    c := p.[4] + e * (e * a - p.[3])
+    c <- p.[4] + e * (e * a - p.[3])
     a <- a - d
 
-    let quadExecuted = ref false
+    let mutable quadExecuted = false
     let quad () =
-        if not !quadExecuted
+        if not quadExecuted
         then
-            p.[2] <- !c / !b
+            p.[2] <- c / b
             quadroots p r
             for k in 1..2 do
                 for j in 1..2 do
                     r.[j, k+2] <- r.[j, k]
             p.[1] <- -p.[1]
-            p.[2] <- !b
+            p.[2] <- b
             quadroots p r
             for k in 1..4 do
                 r.[1,k] <- r.[1,k] - e
-            quadExecuted := true
+            quadExecuted <- true
 
     p.[1] <- 0.5 * a
-    p.[2] <- (p.[1] * p.[1] - !c) * 0.25
-    p.[3] <- !b * !b / -64.0
+    p.[2] <- (p.[1] * p.[1] - c) * 0.25
+    p.[3] <- b * b / -64.0
     if p.[3] < 0.0
     then
         cubicroots p r
@@ -473,43 +473,43 @@ let private biquadroots (p: float[]) (r: float[,]) =
             then
                 d <- r.[1, k] * 4.0
                 a <- a + d
-                if a >= 0.0 && !b >= 0.0
+                if a >= 0.0 && b >= 0.0
                 then
                     p.[1] <- sqrt d
-                elif a <= 0.0 && !b <= 0.0
+                elif a <= 0.0 && b <= 0.0
                 then
                     p.[1] <- sqrt d
                 else
                     p.[1] <- -(sqrt d)
-                b := 0.5 * (a + !b / p.[1])
+                b <- 0.5 * (a + b / p.[1])
                 quad ()
                 k <- 4
             k <- k + 1
 
-    if not !quadExecuted && p.[2] < 0.0
+    if not quadExecuted && p.[2] < 0.0
     then
-        b := sqrt !c
-        d <- !b + !b - a
+        b <- sqrt c
+        d <- b + b - a
         p.[1] <- 0.0
         if d > 0.0
         then
             p.[1] <- sqrt d
-    elif not !quadExecuted
+    elif not quadExecuted
     then
         if p.[1] > 0.0
         then
-            b := (sqrt p.[2]) * 2.0 + p.[1]
+            b <- (sqrt p.[2]) * 2.0 + p.[1]
         else
-            b := -(sqrt p.[2]) * 2.0 + p.[1]
+            b <- -(sqrt p.[2]) * 2.0 + p.[1]
 
-        if !b <> 0.0
+        if b <> 0.0
         then
             p.[1] <- 0.0
         else
             for k in 1..4 do
                 r.[1, k] <- -e
                 r.[2, k] <- 0.0
-            quadExecuted := true
+            quadExecuted <- true
 
     quad ()
 
@@ -610,15 +610,16 @@ let EEOverlapArea (e1: Types.Ellipse) (e2: Types.Ellipse) : (float * float[] * f
         printf "nroots = %d\n" nroots
 #endif
 
-        let ychk = [|
-            for i in 1 .. nroots do
-                if abs r.[2, i] < EPS
-                then
-                    yield r.[1, i] * b1
+        let ychk = Array.init nroots (fun _ -> Double.MaxValue)
+        let mutable nychk = 0
+        for i in 1 .. nroots do
+            if abs r.[2, i] < EPS
+            then
+                ychk.[nychk] <- r.[1, i] * b1
+                nychk <- nychk + 1
 #if DEBUG_LOG
-                    printf "ROOT is Real,  i=%d --> %f (B1=%f)\n" i r.[1, i] b1
+                printf "ROOT is Real,  i=%d --> %f (B1=%f)\n" i r.[1, i] b1
 #endif
-        |]
         Array.sortInPlace ychk
 
 #if DEBUG_LOG
@@ -627,7 +628,6 @@ let EEOverlapArea (e1: Types.Ellipse) (e2: Types.Ellipse) : (float * float[] * f
             printf "\t j=%d, ychk=%f\n" j ychk.[j]
 #endif
 
-        let nychk = Array.length ychk
         let mutable nintpts = 0
 
         let xint = Array.zeroCreate 4
index 82298e2..3d0879e 100644 (file)
@@ -6,8 +6,11 @@ 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)
 
-// TODO: measure performance with a struct.
-type private Node<'k, 'v> = { mutable key: 'k; value : 'v } with
+[<Struct>]
+type private Node<'k, 'v> =
+    val mutable key: 'k
+    val value: 'v
+    new (k, v) = { key = k; value = v }
     override this.ToString () = sprintf "%A -> %A" this.key this.value
 
 type Heap<'k, 'v> (kComparer : IComparer<'k>) =
@@ -32,9 +35,9 @@ type Heap<'k, 'v> (kComparer : IComparer<'k>) =
             a.[max] <- tmp
             heapUp max
 
+    // Check the integrity of the heap, use for debug only.
     let rec checkIntegrity (i: int) : bool =
         let l, r = left i, right i
-
         let leftIntegrity =
             if l < a.Count
             then
@@ -62,11 +65,11 @@ type Heap<'k, 'v> (kComparer : IComparer<'k>) =
             (this :> IEnumerable<'k * 'v>).GetEnumerator() :> System.Collections.IEnumerator
 
     member this.Next () : 'k * 'v =
-        let { key = key; value = value } = a.[0]
+        let node = a.[0]
         a.[0] <- a.[a.Count - 1]
         a.RemoveAt(a.Count - 1)
         heapUp 0
-        key, value
+        node.key, node.value
 
     member this.RemoveNext () =
         a.[0] <- a.[a.Count - 1]
@@ -74,7 +77,7 @@ type Heap<'k, 'v> (kComparer : IComparer<'k>) =
         heapUp 0
 
     member this.Add (key: 'k) (value: 'v) =
-        a.Add({ key = key; value = value })
+        a.Add(Node(key, value))
 
         let mutable i = a.Count - 1
         while i > 0 && kComparer.Compare(a.[parent i].key, a.[i].key) < 0 do
index cee21c7..b9a31f8 100644 (file)
@@ -45,6 +45,7 @@ let suppressMConnections (img: Matrix<byte>) =
             then
                 img.[i, j] <- 0uy
 
+
 let findEdges (img: Image<Gray, float32>) : Matrix<byte> * Image<Gray, float> * Image<Gray, float> =
     let w = img.Width
     let h = img.Height
@@ -76,10 +77,11 @@ let findEdges (img: Image<Gray, float32>) : Matrix<byte> * Image<Gray, float> *
     CvInvoke.CartToPolar(xGradient, yGradient, magnitudes, angles) // Compute the magnitudes (without angles).
 
     let thresholdHigh, thresholdLow =
-        let sensibility = 0.1
+        let sensibilityHigh = 0.1
+        let sensibilityLow = 0.1
         use magnitudesByte = magnitudes.Convert<byte>()
         let threshold = CvInvoke.Threshold(magnitudesByte, magnitudesByte, 0.0, 1.0, CvEnum.ThresholdType.Otsu ||| CvEnum.ThresholdType.Binary)
-        threshold + (sensibility * threshold), threshold - (sensibility * threshold)
+        threshold + (sensibilityHigh * threshold), threshold - (sensibilityLow * threshold)
 
     // Non-maximum suppression.
     use nms = new Matrix<byte>(xGradient.Size)
@@ -106,28 +108,20 @@ let findEdges (img: Image<Gray, float32>) : Matrix<byte> * Image<Gray, float> *
 
                 let mNeigbors (sign: int) : float =
                     if angle < Math.PI / 4.
-                    then
-                        ratio1 * magnitudes.Data.[i, j + sign] + ratio2 * magnitudes.Data.[i + sign, j + sign]
+                    then ratio1 * magnitudes.Data.[i, j + sign] + ratio2 * magnitudes.Data.[i + sign, j + sign]
                     elif angle < Math.PI / 2.
-                    then
-                        ratio2 * magnitudes.Data.[i + sign, j + sign] + ratio1 * magnitudes.Data.[i + sign, j]
+                    then ratio2 * magnitudes.Data.[i + sign, j + sign] + ratio1 * magnitudes.Data.[i + sign, j]
                     elif angle < 3.0 * Math.PI / 4.
-                    then
-                        ratio1 * magnitudes.Data.[i + sign, j] + ratio2 * magnitudes.Data.[i + sign, j - sign]
+                    then ratio1 * magnitudes.Data.[i + sign, j] + ratio2 * magnitudes.Data.[i + sign, j - sign]
                     elif angle < Math.PI
-                    then
-                        ratio2 * magnitudes.Data.[i + sign, j - sign] + ratio1 * magnitudes.Data.[i, j - sign]
+                    then ratio2 * magnitudes.Data.[i + sign, j - sign] + ratio1 * magnitudes.Data.[i, j - sign]
                     elif angle < 5. * Math.PI / 4.
-                    then
-                        ratio1 * magnitudes.Data.[i, j - sign] + ratio2 * magnitudes.Data.[i - sign, j - sign]
+                    then ratio1 * magnitudes.Data.[i, j - sign] + ratio2 * magnitudes.Data.[i - sign, j - sign]
                     elif angle < 3. * Math.PI / 2.
-                    then
-                        ratio2 * magnitudes.Data.[i - sign, j - sign] + ratio1 * magnitudes.Data.[i - sign, j]
+                    then ratio2 * magnitudes.Data.[i - sign, j - sign] + ratio1 * magnitudes.Data.[i - sign, j]
                     elif angle < 7. * Math.PI / 4.
-                    then
-                        ratio1 * magnitudes.Data.[i - sign, j] + ratio2 * magnitudes.Data.[i - sign, j + sign]
-                    else
-                        ratio2 * magnitudes.Data.[i - sign, j + sign] + ratio1 * magnitudes.Data.[i, j + sign]
+                    then ratio1 * magnitudes.Data.[i - sign, j] + ratio2 * magnitudes.Data.[i - sign, j + sign]
+                    else ratio2 * magnitudes.Data.[i - sign, j + sign] + ratio1 * magnitudes.Data.[i, j + sign]
 
                 let m = magnitudes.Data.[i, j]
                 if m >= thresholdLow && m > mNeigbors 1 && m > mNeigbors -1
index 167e618..fdae997 100644 (file)
@@ -22,7 +22,7 @@ let doAnalysis (img: Image<Bgr, byte>) (name: string) (config: Config) : Cell li
 
     logTime "areaOpen 1" (fun () -> ImgTools.areaOpenF filteredGreen config.Parameters.initialAreaOpen)
 
-    config.RBCRadius <- logTime "Granulometry" (fun() -> Granulometry.findRadius (filteredGreen.Convert<Gray, byte>()) (10, 100) 0.5 |> float)
+    config.RBCRadius <- logTime "Granulometry" (fun() -> Granulometry.findRadius (filteredGreen.Convert<Gray, byte>()) (10, 100) 0.3 |> float)
 
     let secondAreaOpen = int <| config.RBCArea / 3.
     if secondAreaOpen > config.Parameters.initialAreaOpen
index c39541e..2022c26 100644 (file)
@@ -36,7 +36,7 @@ type MatchingEllipses (radiusMin: float) =
     let ellipses = List<EllipseScoreFlaggedKd>()
 
     // All ellipses with a score below this are removed.
-    let matchingScoreThreshold2 = 25. * radiusMin // 600.
+    let matchingScoreThreshold2 = 20. * radiusMin
 
     member this.Add (e: Ellipse) =
         ellipses.Add(EllipseScoreFlaggedKd(0.0, e))
index 94d2999..6531da8 100644 (file)
@@ -72,7 +72,7 @@ let main args =
                 factorWindowSize = 2.0
 
                 darkStainLevel = 0.25 // Lower -> more sensitive. 0.3
-                maxDarkStainRatio = 0.1
+                maxDarkStainRatio = 0.1 // 10 %
 
                 infectionArea = 0.012 // 1.2 %
                 infectionLevel = 1.12 // Lower -> more sensitive.
@@ -81,7 +81,7 @@ let main args =
                 stainLevel = 1.1 // Lower -> more sensitive.
                 maxStainRatio = 0.12 // 12 %
 
-                standardDeviationMaxRatio = 0.58 // 0.55
+                standardDeviationMaxRatio = 0.5 // 0.55
                 minimumCellArea = 0.5 })
 
         match mode with