From 53440e757b4a4ab2a81b0f6a5dd1a2002c0133ba Mon Sep 17 00:00:00 2001 From: Greg Burri Date: Mon, 4 Jan 2016 13:09:54 +0100 Subject: [PATCH] Cleaning, micro-optimizations. --- Parasitemia/Parasitemia/Classifier.fs | 32 ++++++----- Parasitemia/Parasitemia/EEOver.fs | 62 ++++++++++----------- Parasitemia/Parasitemia/Heap.fs | 15 +++-- Parasitemia/Parasitemia/ImgTools.fs | 30 ++++------ Parasitemia/Parasitemia/MainAnalysis.fs | 2 +- Parasitemia/Parasitemia/MatchingEllipses.fs | 2 +- Parasitemia/Parasitemia/Program.fs | 4 +- 7 files changed, 74 insertions(+), 73 deletions(-) diff --git a/Parasitemia/Parasitemia/Classifier.fs b/Parasitemia/Parasitemia/Classifier.fs index f17c0b9..0d21375 100644 --- a/Parasitemia/Parasitemia/Classifier.fs +++ b/Parasitemia/Parasitemia/Classifier.fs @@ -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 diff --git a/Parasitemia/Parasitemia/EEOver.fs b/Parasitemia/Parasitemia/EEOver.fs index ef21051..528632a 100644 --- a/Parasitemia/Parasitemia/EEOver.fs +++ b/Parasitemia/Parasitemia/EEOver.fs @@ -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 diff --git a/Parasitemia/Parasitemia/Heap.fs b/Parasitemia/Parasitemia/Heap.fs index 82298e2..3d0879e 100644 --- a/Parasitemia/Parasitemia/Heap.fs +++ b/Parasitemia/Parasitemia/Heap.fs @@ -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 +[] +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 diff --git a/Parasitemia/Parasitemia/ImgTools.fs b/Parasitemia/Parasitemia/ImgTools.fs index cee21c7..b9a31f8 100644 --- a/Parasitemia/Parasitemia/ImgTools.fs +++ b/Parasitemia/Parasitemia/ImgTools.fs @@ -45,6 +45,7 @@ let suppressMConnections (img: Matrix) = then img.[i, j] <- 0uy + let findEdges (img: Image) : Matrix * Image * Image = let w = img.Width let h = img.Height @@ -76,10 +77,11 @@ let findEdges (img: Image) : Matrix * Image * 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() 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(xGradient.Size) @@ -106,28 +108,20 @@ let findEdges (img: Image) : Matrix * Image * 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 diff --git a/Parasitemia/Parasitemia/MainAnalysis.fs b/Parasitemia/Parasitemia/MainAnalysis.fs index 167e618..fdae997 100644 --- a/Parasitemia/Parasitemia/MainAnalysis.fs +++ b/Parasitemia/Parasitemia/MainAnalysis.fs @@ -22,7 +22,7 @@ let doAnalysis (img: Image) (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()) (10, 100) 0.5 |> float) + config.RBCRadius <- logTime "Granulometry" (fun() -> Granulometry.findRadius (filteredGreen.Convert()) (10, 100) 0.3 |> float) let secondAreaOpen = int <| config.RBCArea / 3. if secondAreaOpen > config.Parameters.initialAreaOpen diff --git a/Parasitemia/Parasitemia/MatchingEllipses.fs b/Parasitemia/Parasitemia/MatchingEllipses.fs index c39541e..2022c26 100644 --- a/Parasitemia/Parasitemia/MatchingEllipses.fs +++ b/Parasitemia/Parasitemia/MatchingEllipses.fs @@ -36,7 +36,7 @@ type MatchingEllipses (radiusMin: float) = let ellipses = List() // 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)) diff --git a/Parasitemia/Parasitemia/Program.fs b/Parasitemia/Parasitemia/Program.fs index 94d2999..6531da8 100644 --- a/Parasitemia/Parasitemia/Program.fs +++ b/Parasitemia/Parasitemia/Program.fs @@ -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 -- 2.43.0