From 5ac2dedf8ead5275ac216e0b41829ab39c843800 Mon Sep 17 00:00:00 2001 From: Greg Burri Date: Sun, 20 Dec 2015 14:06:49 +0100 Subject: [PATCH] * Remove the 'DoG' filter. * Fix the are opening. --- Parasitemia/Parasitemia/Config.fs | 4 +- Parasitemia/Parasitemia/EEOver.fs | 66 +++--- Parasitemia/Parasitemia/Granulometry.fs | 8 +- Parasitemia/Parasitemia/ImgTools.fs | 228 +++++++++++--------- Parasitemia/Parasitemia/MainAnalysis.fs | 50 ++--- Parasitemia/Parasitemia/MatchingEllipses.fs | 2 +- Parasitemia/Parasitemia/Program.fs | 8 +- Parasitemia/Parasitemia/Types.fs | 3 + 8 files changed, 194 insertions(+), 175 deletions(-) diff --git a/Parasitemia/Parasitemia/Config.fs b/Parasitemia/Parasitemia/Config.fs index 464fd7e..48a6334 100644 --- a/Parasitemia/Parasitemia/Config.fs +++ b/Parasitemia/Parasitemia/Config.fs @@ -12,9 +12,7 @@ type Config = { minRBCSize: float maxRBCSize: float - doGSigma1: float - doGSigma2: float - doGLowFreqPercentageReduction: float + preFilterSigma: float // Ellipse. factorNbPick: float diff --git a/Parasitemia/Parasitemia/EEOver.fs b/Parasitemia/Parasitemia/EEOver.fs index 52b5963..1a6bb5f 100644 --- a/Parasitemia/Parasitemia/EEOver.fs +++ b/Parasitemia/Parasitemia/EEOver.fs @@ -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 5 let r = Array2D.zeroCreate 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 diff --git a/Parasitemia/Parasitemia/Granulometry.fs b/Parasitemia/Parasitemia/Granulometry.fs index a9a102a..db299dc 100644 --- a/Parasitemia/Parasitemia/Granulometry.fs +++ b/Parasitemia/Parasitemia/Granulometry.fs @@ -10,7 +10,7 @@ open Utils // 'range': a minimum and maximum radius. // 'scale': <= 1.0, to speed up the process. -let findRadius (img: Image) (range: int * int) (scale: float) : int = +let findRadius (img: Image) (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) (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 diff --git a/Parasitemia/Parasitemia/ImgTools.fs b/Parasitemia/Parasitemia/ImgTools.fs index e9cbbe9..24c6998 100644 --- a/Parasitemia/Parasitemia/ImgTools.fs +++ b/Parasitemia/Parasitemia/ImgTools.fs @@ -26,7 +26,13 @@ let gaussianFilter (img : Image<'TColor, 'TDepth>) (standardDeviation : float) : img.SmoothGaussian(size, size, standardDeviation, standardDeviation) -let findMaxima (img: Image) : IEnumerable> = +type Points = HashSet + +let drawPoints (img: Image) (points: Points) (intensity: byte) = + for p in points do + img.Data.[p.Y, p.X, 0] <- intensity + +let findMaxima (img: Image) : IEnumerable = use suppress = new Image(img.Size) let w = img.Width let h = img.Height @@ -90,159 +96,189 @@ let findMaxima (img: Image) : IEnumerable> = if maxima.Count > 0 then result.AddRange(maxima) - result.Select(fun l -> HashSet(l)) + result.Select(fun l -> Points(l)) type PriorityQueue () = - let q = List>() // 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([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 - mutable intensity: byte option - mutable state: MaximaState } - +[] +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) (area: int) = let w = img.Width let h = img.Height + let imgData = img.Data + + let areas = List(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) + 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() - 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() + 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) (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) (startPoints: List) : Li let w = img.Width let h = img.Height - let pointChecked = HashSet() + let pointChecked = Points() let pointToCheck = List(startPoints); let data = img.Data diff --git a/Parasitemia/Parasitemia/MainAnalysis.fs b/Parasitemia/Parasitemia/MainAnalysis.fs index 155bfcc..39383fc 100644 --- a/Parasitemia/Parasitemia/MainAnalysis.fs +++ b/Parasitemia/Parasitemia/MainAnalysis.fs @@ -20,46 +20,25 @@ let doAnalysis (img: Image) (name: string) (config: Config) : Cell li use green = scaledImg.Item(1) use red = scaledImg.Item(2) - let greenFloat = green.Convert() - 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() // 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() - use yEdges = filteredGreen.Convolution(sobelKernel.Transpose()).Convert() + use xEdges = filteredGreenFloat.Convolution(sobelKernel).Convert() + use yEdges = filteredGreenFloat.Convolution(sobelKernel.Transpose()).Convert() let xEdgesData = xEdges.Data let yEdgesData = yEdges.Data @@ -91,9 +70,9 @@ let doAnalysis (img: Image) (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) (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 diff --git a/Parasitemia/Parasitemia/MatchingEllipses.fs b/Parasitemia/Parasitemia/MatchingEllipses.fs index 92a767a..99e3cb5 100644 --- a/Parasitemia/Parasitemia/MatchingEllipses.fs +++ b/Parasitemia/Parasitemia/MatchingEllipses.fs @@ -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 diff --git a/Parasitemia/Parasitemia/Program.fs b/Parasitemia/Parasitemia/Program.fs index fdbc908..86de435 100644 --- a/Parasitemia/Parasitemia/Program.fs +++ b/Parasitemia/Parasitemia/Program.fs @@ -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 diff --git a/Parasitemia/Parasitemia/Types.fs b/Parasitemia/Parasitemia/Types.fs index a551142..634a7b7 100644 --- a/Parasitemia/Parasitemia/Types.fs +++ b/Parasitemia/Parasitemia/Types.fs @@ -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 -- 2.43.0