From dec96d50e56e1932bbfa09e6bedf90d6b707ccbd Mon Sep 17 00:00:00 2001 From: Greg Burri Date: Fri, 11 Dec 2015 16:47:55 +0100 Subject: [PATCH] Finding ellipses and parasites. --- Parasitemia/Parasitemia/App.config | 2 +- Parasitemia/Parasitemia/Config.fs | 13 +- Parasitemia/Parasitemia/EEOver.fs | 145 ++++++++++++++++-- Parasitemia/Parasitemia/Ellipse.fs | 32 ++-- Parasitemia/Parasitemia/ImageAnalysis.fs | 31 ++-- Parasitemia/Parasitemia/ImgTools.fs | 22 ++- Parasitemia/Parasitemia/KMedians.fs | 73 +++++----- Parasitemia/Parasitemia/KdTree.fs | 154 ++++++++++++++++++++ Parasitemia/Parasitemia/MatchingEllipses.fs | 102 +++++++++---- Parasitemia/Parasitemia/Parasitemia.fsproj | 18 ++- Parasitemia/Parasitemia/ParasitesMarker.fs | 39 +++++ Parasitemia/Parasitemia/Program.fs | 35 ++++- Parasitemia/Parasitemia/Types.fs | 13 +- Parasitemia/Parasitemia/packages.config | 3 + 14 files changed, 554 insertions(+), 128 deletions(-) create mode 100644 Parasitemia/Parasitemia/KdTree.fs create mode 100644 Parasitemia/Parasitemia/ParasitesMarker.fs diff --git a/Parasitemia/Parasitemia/App.config b/Parasitemia/Parasitemia/App.config index 640d253..06fd2b0 100644 --- a/Parasitemia/Parasitemia/App.config +++ b/Parasitemia/Parasitemia/App.config @@ -1,7 +1,7 @@  - + diff --git a/Parasitemia/Parasitemia/Config.fs b/Parasitemia/Parasitemia/Config.fs index b22a4fd..75aa95f 100644 --- a/Parasitemia/Parasitemia/Config.fs +++ b/Parasitemia/Parasitemia/Config.fs @@ -1,9 +1,20 @@ module Config type Config = { + scale: float + doGSigma1: float doGSigma2: float doGLowFreqPercentageReduction: float - scale: float + // Parasites detection. + darkStainLevel: float + + stainSigma: float + stainLevel: float + stainSpreadRequired: float + + infectionSigma: float + infectionLevel: float + infectionPixelsRequired: int } \ No newline at end of file diff --git a/Parasitemia/Parasitemia/EEOver.fs b/Parasitemia/Parasitemia/EEOver.fs index e366b72..4891059 100644 --- a/Parasitemia/Parasitemia/EEOver.fs +++ b/Parasitemia/Parasitemia/EEOver.fs @@ -55,6 +55,16 @@ let private istanpt (x: float) (y: float) (a1: float) (b1: float) (aa: float) (b let test1 = ellipse2tr x1 y1 aa bb cc dd ee ff let test2 = ellipse2tr x2 y2 aa bb cc dd ee ff +#if DEBUG_LOG + printf "\t\t--- debug istanpt with (x,y)=(%f, %f), A1=%f, B1=%f\n" x y a1 b1 + printf "theta=%f\n" theta + printf "eps_Radian=%f\n" eps_radian + printf "(x1, y1)=(%f, %f)\n" x1 y1 + printf "(x2, y2)=(%f, %f)\n" x2 y2 + printf "test1=%f\n" test1 + printf "test2=%f\n" test2 +#endif + if test1 * test2 > 0.0 then TANGENT_POINT else INTERSECTION_POINT @@ -102,6 +112,9 @@ let private twointpts (x: float[]) (y: float[]) (a1: float) (b1: float) (phi_1: if area1 < 0.0 then +#if DEBUG_LOG + printf "TWO area1=%f\n" area1 +#endif area1 <- area1 + a1 * b1 let cosphi = cos (phi_1 - phi_2) @@ -161,6 +174,9 @@ let private twointpts (x: float[]) (y: float[]) (a1: float) (b1: float) (phi_1: let mutable area2 = 0.5 * (a2 * b2 * (theta2 - theta1) + trsign * abs (x1_tr * y2_tr - x2_tr * y1_tr)) if area2 < 0.0 then +#if DEBUG_LOG + printf "TWO area2=%f\n" area2 +#endif area2 <- area2 + a2 * b2 area1 + area2 @@ -172,9 +188,12 @@ let private threeintpts (xint: float[]) (yint: float[]) (a1: float) (b1: float) for i in 0..2 do if istanpt xint.[i] yint.[i] a1 b2 aa bb cc dd ee ff = TANGENT_POINT then - tanpts <- tanpts + 1; - tanindex <- i; - + tanpts <- tanpts + 1 + tanindex <- i +#if DEBUG_LOG + printf "tanindex=%d\n" tanindex +#endif + if tanpts <> 1 then -1.0 @@ -205,25 +224,40 @@ let private fourintpts (xint: float[]) (yint: float[]) (a1: float) (b1: float) ( xint.[i] <- if xint.[i] < 0.0 then -a1 else a1 theta.[i] <- if yint.[i] < 0.0 then 2.0 * Math.PI - acos (xint.[i] / a1) else acos (xint.[i] / a1) +#if DEBUG_LOG + for k in 0..3 do + printf "k=%d: Theta = %f, xint=%f, yint=%f\n" k theta.[k] xint.[k] yint.[k] +#endif + for j in 1..3 do let tmp0 = theta.[j] let tmp1 = xint.[j] let tmp2 = yint.[j] let mutable k = j - 1 + let mutable k2 = 0 while k >= 0 do if theta.[k] <= tmp0 then + k2 <- k + 1 k <- -1 else theta.[k+1] <- theta.[k] xint.[k+1] <- xint.[k] yint.[k+1] <- yint.[k] - k <- k - 1 + k <- k - 1 + k2 <- k + 1 + + theta.[k2] <- tmp0 + xint.[k2] <- tmp1 + yint.[k2] <- tmp2 - theta.[k+1] <- tmp0 - xint.[k+1] <- tmp1 - yint.[k+1] <- tmp2 + +#if DEBUG_LOG + printf "AFTER sorting\n" + 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])) @@ -267,20 +301,36 @@ let private fourintpts (xint: float[]) (yint: float[]) (a1: float) (b1: float) ( if area5 < 0.0 then +#if DEBUG_LOG + printf "\n\t\t-------------> area5 is negativ (%f). Add: pi*A2*B2=%f <------------\n" area5 area_2 +#endif area5 <- area5 + area_2 if area4 < 0.0 then +#if DEBUG_LOG + printf "\n\t\t-------------> area4 is negativ (%f). Add: pi*A2*B2=%f <------------\n" area4 area_2 +#endif area4 <- area4 + area_2 if area3 < 0.0 then +#if DEBUG_LOG + printf "\n\t\t-------------> area3 is negativ (%f). Add: pi*A2*B2=%f <------------\n" area3 area_1 +#endif 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 area2 <- area2 + area_1 +#if DEBUG_LOG + printf "\narea1=%f, area2=%f area3=%f, area4=%f, area5=%f\n\n" area1 area2 area3 area4 area5 +#endif + area1 + area2 + area3 + area4 + area5 @@ -463,11 +513,9 @@ let private biquadroots (p: float[]) (r: float[,]) = quad () -open Types - -let EEOverlapArea (e1: Ellipse) (e2: Ellipse) : float = - let { cx = h1; cy = k1; a = a1; b = b1; alpha = phi_1 } = e1 - let { cx = h2; cy = k2; a = a2; b = b2; alpha = phi_2 } = e2 +let EEOverlapArea (e1: Types.Ellipse) (e2: Types.Ellipse) : float = + 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 then @@ -480,6 +528,10 @@ let EEOverlapArea (e1: Ellipse) (e2: Ellipse) : float = 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 + let cosphi = cos phi_2r let cosphi2 = cosphi ** 2.0 let sinphi = sin phi_2r @@ -507,6 +559,11 @@ let EEOverlapArea (e1: Ellipse) (e2: Ellipse) : float = a1 * a1 * a1 * a1 * aa * aa + b1 * b1 * (a1 * a1 * (bb * bb - 2.0 * aa * cc) + b1 * b1 * cc * cc) |] +#if DEBUG_LOG + for i in 0..4 do + printf "cy[%d]=%f\n" i cy.[i] +#endif + let py = Array.zeroCreate 5 let r = Array2D.zeroCreate 3 5 @@ -516,6 +573,10 @@ let EEOverlapArea (e1: Ellipse) (e2: Ellipse) : float = for i in 0 .. 3 do py.[4-i] <- cy.[i] / cy.[4] py.[0] <- 1.0 +#if DEBUG_LOG + for i in 0..4 do + printf "py[%d]=%f\n" i py.[i] +#endif biquadroots py r 4 @@ -544,14 +605,27 @@ let EEOverlapArea (e1: Ellipse) (e2: Ellipse) : float = else 0 +#if DEBUG_LOG + 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 +#if DEBUG_LOG + printf "ROOT is Real, i=%d --> %f (B1=%f)\n" i r.[1, i] b1 +#endif |] Array.sortInPlace ychk +#if DEBUG_LOG + printf "nychk=%d\n" ychk.Length + for j in 0 .. ychk.Length - 1 do + printf "\t j=%d, ychk=%f\n" j ychk.[j] +#endif + let nychk = Array.length ychk let mutable nintpts = 0 @@ -562,30 +636,61 @@ let EEOverlapArea (e1: Ellipse) (e2: Ellipse) : float = let mutable i = 0 while returnValue = 0.0 && i < nychk do +#if DEBUG_LOG + printf "------------->i=%d (nychk=%d)\n" i nychk +#endif + if not (i < nychk - 1 && abs (ychk.[i] - ychk.[i+1]) < EPS / 2.0) then +#if DEBUG_LOG + printf "check intersecting points. nintps is %d" nintpts +#endif + let x1 = if abs ychk.[i] > b1 then 0.0 else a1 * sqrt (1.0 - (ychk.[i] * ychk.[i]) / (b1 * b1)) let x2 = -x1 +#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 + if abs (ellipse2tr x1 ychk.[i] aa bb cc dd ee ff) < EPS then nintpts <- nintpts + 1 +#if DEBUG_LOG + printf "first if x1. acc nintps=%d\n" nintpts +#endif if nintpts > 4 then returnValue <- -1.0 else xint.[nintpts-1] <- x1 yint.[nintpts-1] <- ychk.[i] +#if DEBUG_LOG + printf "nintpts=%d, xint=%f, x2=%f, i=%d, yint=%f\n" nintpts x1 x2 i ychk.[i] +#endif if returnValue <> -1.0 && abs (ellipse2tr x2 ychk.[i] aa bb cc dd ee ff) < EPS && abs (x2 - x1) > EPS then nintpts <- nintpts + 1 +#if DEBUG_LOG + printf "first if x2. nintps=%d, Dx=%f (eps2=%f) \n" nintpts (abs (x2 - x1)) EPS +#endif if nintpts > 4 then returnValue <- -1.0 else xint.[nintpts-1] <- x2 yint.[nintpts-1] <- ychk.[i] + +#if DEBUG_LOG + printf "nintpts=%d, x1=%f, xint=%f, i=%d, yint=%f\n" nintpts x1 x2 i ychk.[i] +#endif + +#if DEBUG_LOG + else + printf "i=%d, multiple roots: %f <--------> %f. continue\n" i ychk.[i] ychk.[i-1] +#endif i <- i + 1 @@ -596,9 +701,17 @@ let EEOverlapArea (e1: Ellipse) (e2: Ellipse) : float = 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 -> 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 -> twointpts xint yint a1 b1 phi_1 a2 b2 h2_tr k2_tr phi_2 aa bb cc dd ee ff + | 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 -> +#if DEBUG_LOG + printf "check twointpts\n" +#endif + twointpts xint yint a1 b1 phi_1 a2 b2 h2_tr k2_tr phi_2 aa bb cc dd ee ff | 3 -> threeintpts xint yint a1 b1 phi_1 a2 b2 h2_tr k2_tr phi_2 aa bb cc dd ee ff | 4 -> fourintpts xint yint a1 b1 phi_1 a2 b2 h2_tr k2_tr phi_2 aa bb cc dd ee ff | _ -> -1.0 - diff --git a/Parasitemia/Parasitemia/Ellipse.fs b/Parasitemia/Parasitemia/Ellipse.fs index 62273b5..5bfcc33 100644 --- a/Parasitemia/Parasitemia/Ellipse.fs +++ b/Parasitemia/Parasitemia/Ellipse.fs @@ -31,20 +31,20 @@ let private goldenSectionSearch (f: float -> float) (nbIter: int) (xmin: float) if fc < fd then - b <- d; - d <- c; - c <- b - gr * (b - a); + b <- d + d <- c + c <- b - gr * (b - a) else - a <- c; - c <- d; - d <- a + gr * (b - a); + a <- c + c <- d + d <- a + gr * (b - a) let x = (b + a) / 2.0 x, f x 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 = 4 + let accuracy_extremum_search_2 = 3 // p3 as the referencial. let p1x = p1x - p3x @@ -139,7 +139,7 @@ let ellipse (p1x: float) (p1y: float) (m1: float) (p2x: float) (p2y: float) (m2: let cx = cx + p3x let cy = cy + p3y - Some { cx = cx; cy = cy; a = r1e; b = r2e; alpha = alpha } + Some (Types.Ellipse(cx, cy, r1e, r2e, alpha)) else None @@ -195,13 +195,13 @@ let find (edges: Matrix) (windowSize: int) (factorNbPick: float) : Types.Ellipse list = - let increment = windowSize / 4; + let increment = windowSize / 4 let r1, r2 = radiusRange let radiusTolerance = (r2 - r1) * 0.2 - let minimumDistance = (r2 / 1.5) ** 2.0; - let squaredDistance x1 y1 x2 y2 = (x1 - x2) ** 2.0 + (y1 - y2) ** 2.0; + let minimumDistance = (r2 / 1.5) ** 2.0 + let squaredDistance x1 y1 x2 y2 = (x1 - x2) ** 2.0 + (y1 - y2) ** 2.0 let h = edges.Height let w = edges.Width @@ -214,9 +214,9 @@ let find (edges: Matrix) let xDirData = xDir.Data let yDirData = yDir.Data - let rng = Random() + let rng = Random(42) - let ellipses = MatchingEllipses () + let ellipses = MatchingEllipses(r1) for window_i in -windowSize + increment .. increment .. h - increment do for window_j in -windowSize + increment .. increment .. w - increment do @@ -257,8 +257,8 @@ let find (edges: Matrix) match areVectorsValid p1xf p1yf p2xf p2yf -xDirData.[p1y, p1x, 0] -yDirData.[p1y, p1x, 0] -xDirData.[p2y, p2x, 0] -yDirData.[p2y, p2x, 0] with | Some (m1, m2) -> 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 -> + | 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 -> ellipses.Add e | _ -> () | _ -> () diff --git a/Parasitemia/Parasitemia/ImageAnalysis.fs b/Parasitemia/Parasitemia/ImageAnalysis.fs index 1215965..0b96757 100644 --- a/Parasitemia/Parasitemia/ImageAnalysis.fs +++ b/Parasitemia/Parasitemia/ImageAnalysis.fs @@ -41,7 +41,7 @@ let doAnalysis (img: Image) (config: Config) : Result = //let test = greenMatrix.[10, 10] use filteredGreen = (gaussianFilter green config.doGSigma1) - config.doGLowFreqPercentageReduction * (gaussianFilter green config.doGSigma2) - + use sobelKernel = new ConvolutionKernelF(array2D [[ 1.0f; 0.0f; -1.0f ] [ 2.0f; 0.0f; -2.0f ] @@ -75,19 +75,32 @@ let doAnalysis (img: Image) (config: Config) : Result = use magnitudesByte = ((magnitudes / !max) * 255.0).Convert() // Otsu from OpenCV only support 'byte'. use edges = new Matrix(xEdges.Size) - let threshold = CvInvoke.Threshold(magnitudesByte, edges, 0.0, 1.0, CvEnum.ThresholdType.Otsu ||| CvEnum.ThresholdType.Binary) - thin edges - removeArea edges 12 + let threshold = CvInvoke.Threshold(magnitudesByte, edges, 0.0, 1.0, CvEnum.ThresholdType.Otsu ||| CvEnum.ThresholdType.Binary) + +// let filteredGreenMat = new Matrix(filteredGreen.Size) +// filteredGreen.CopyTo(filteredGreenMat) + let parasites = ParasitesMarker.find green filteredGreen config + + saveImg parasites.darkStain "parasites_dark_stain.png" + saveImg parasites.stain "parasites_stain.png" + saveImg parasites.infection "parasites_infection.png" + + logTime "Finding edges" (fun() -> + thin edges) + + logTime "Removing small connected components" (fun () -> + removeArea edges 12) saveMat (edges * 255.0) "edges.png" let radiusRange = config.scale * 20.0, config.scale * 40.0 let windowSize = roundInt (1.6 * (snd radiusRange)) - let factorNbPick = 1.0; - let ellipses = Ellipse.find edges xEdges yEdges radiusRange windowSize factorNbPick - - drawEllipse img (List.head ellipses) (Bgr(0.0, 255.0, 255.0)) - saveImg img "ellipses.png" + let factorNbPick = 1.5 + let ellipses = logTime "Finding ellipses" (fun () -> + Ellipse.find edges xEdges yEdges radiusRange windowSize factorNbPick) + + drawEllipses img ellipses (Bgr(0.0, 255.0, 255.0)) + //saveImg img "ellipses.png" { RBCPositions = []; infectedRBCPositions = []; img = img } diff --git a/Parasitemia/Parasitemia/ImgTools.fs b/Parasitemia/Parasitemia/ImgTools.fs index cc09405..5278313 100644 --- a/Parasitemia/Parasitemia/ImgTools.fs +++ b/Parasitemia/Parasitemia/ImgTools.fs @@ -9,6 +9,15 @@ open Emgu.CV.Structure open Utils +// Normalize image values between 0uy and 255uy. +let normalizeAndConvert (img: Image) : Image = + let min = ref [| 0.0 |] + let minLocation = ref <| [| Point() |] + let max = ref [| 0.0 |] + let maxLocation = ref <| [| Point() |] + img.MinMax(min, max, minLocation, maxLocation) + ((img - (!min).[0]) / ((!max).[0] - (!min).[0]) * 255.0).Convert() + 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) @@ -151,8 +160,8 @@ let drawLine (img: Image<'TColor, 'TDepth>) (color: 'TColor) (x0: float) (y0: fl img.Draw(LineSegment2D(Point(x0, y0), Point(x1, y1)), color, 1); let drawEllipse (img: Image<'TColor, 'TDepth>) (e: Types.Ellipse) (color: 'TColor) = - let cosAlpha = cos e.alpha - let sinAlpha = sin e.alpha + let cosAlpha = cos e.Alpha + let sinAlpha = sin e.Alpha let mutable x0 = 0.0 let mutable y0 = 0.0 @@ -164,8 +173,8 @@ let drawEllipse (img: Image<'TColor, 'TDepth>) (e: Types.Ellipse) (color: 'TColo for theta in 0.0 .. thetaIncrement .. 2.0 * Math.PI do let cosTheta = cos theta let sinTheta = sin theta - let x = e.cx + cosAlpha * e.a * cosTheta - sinAlpha * e.b * sinTheta - let y = e.cy + sinAlpha * e.a * cosTheta + cosAlpha * e.b * sinTheta + let x = e.Cx + cosAlpha * e.A * cosTheta - sinAlpha * e.B * sinTheta + let y = e.Cy + sinAlpha * e.A * cosTheta + cosAlpha * e.B * sinTheta if not first_iteration then @@ -174,4 +183,7 @@ let drawEllipse (img: Image<'TColor, 'TDepth>) (e: Types.Ellipse) (color: 'TColo first_iteration <- false x0 <- x - y0 <- y \ No newline at end of file + y0 <- y + +let drawEllipses (img: Image<'TColor, 'TDepth>) (ellipses: Types.Ellipse list) (color: 'TColor) = + List.iter (fun e -> drawEllipse img e color) ellipses \ No newline at end of file diff --git a/Parasitemia/Parasitemia/KMedians.fs b/Parasitemia/Parasitemia/KMedians.fs index 5ba15db..09d6c8f 100644 --- a/Parasitemia/Parasitemia/KMedians.fs +++ b/Parasitemia/Parasitemia/KMedians.fs @@ -1,45 +1,52 @@ module KMedians +open System.Collections.Generic +open System.Drawing + open Emgu.CV open Emgu.CV.Structure -open System.Drawing +type Result = { + fg: Image + median_bg: float + median_fg: float + d_fg: Image } // Distances to median_fg. -let kmedians (mat: Matrix) (fgFactor: float) : Matrix = +let kmedians (img: Image) (fgFactor: float) : Result = let nbIteration = 3 - - let min = ref 0.0 - let minLocation = ref <| Point() - let max = ref 0.0 - let maxLocation = ref <| Point() - mat.MinMax(min, max, minLocation, maxLocation) + let w = img.Width + let h = img.Height + + let min = ref [| 0.0 |] + let minLocation = ref <| [| Point() |] + let max = ref [| 0.0 |] + let maxLocation = ref <| [| Point() |] + img.MinMax(min, max, minLocation, maxLocation) - let mutable bgValue = !max - (!max - !min) / 4.0 - let mutable fgValue = !min + (!max - !min) / 4.0 - let mutable d_bg = new Matrix(mat.Size) - let mutable d_fg = new Matrix(mat.Size) - - for i in 1..3 do - CvInvoke.Pow(mat - bgValue, 2.0, d_bg) - CvInvoke.Pow(mat - fgValue, 2.0, d_fg) - let fg = (d_fg * fgFactor).Cmp(d_bg, CvEnum.CmpType.LessThan) - - printfn "test" - - - (*d_bg = (imgFloat - color_bg) .^ 2.0; - d_fg = (imgFloat - color_fg) .^ 2.0; - fg = d_fg * Background_weight < d_bg; - imgFilteredFg = imgFloat; - imgFilteredFg(~fg) = nan; - color_fg = median(reshape(imgFilteredFg, [], 1), 'omitnan'); - imgFilteredBg = imgFloat; - imgFilteredBg(fg) = nan; - color_bg = median(reshape(imgFilteredBg, [], 1), 'omitnan'); - *) - + let mutable median_bg = (!max).[0] - ((!max).[0] - (!min).[0]) / 4.0 + let mutable median_fg = (!min).[0] + ((!max).[0] - (!min).[0]) / 4.0 + let mutable d_bg = new Image(img.Size) + let mutable d_fg = new Image(img.Size) + let mutable fg = new Image(img.Size) - new Matrix(mat.Size) + 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) + + 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 }) + + 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 }) + + CvInvoke.Sqrt(d_fg, d_fg) + + { fg = fg; median_bg = median_bg; median_fg = median_fg; d_fg = d_fg } diff --git a/Parasitemia/Parasitemia/KdTree.fs b/Parasitemia/Parasitemia/KdTree.fs new file mode 100644 index 0000000..dd0247a --- /dev/null +++ b/Parasitemia/Parasitemia/KdTree.fs @@ -0,0 +1,154 @@ +module KdTree + +open System + +type I2DCoords = + abstract X : float + abstract Y : float + +// Compare 'e1' and 'e2' by X. +let cmpX (e1: I2DCoords) (e2: I2DCoords) : int = + match e1.X.CompareTo(e2.X) with + | 0 -> match e1.Y.CompareTo(e2.Y) with + | 0 -> e1.GetHashCode().CompareTo(e2.GetHashCode()) + | v -> v + | v -> v + +// Compare 'e1' and 'e2' by Y. +let cmpY (e1: I2DCoords) (e2: I2DCoords) : int = + match e1.Y.CompareTo(e2.Y) with + | 0 -> match e1.X.CompareTo(e2.X) with + | 0 -> e1.GetHashCode().CompareTo(e2.GetHashCode()) + | v -> v + | v -> v + +type Region = { minX: float; maxX: float; minY: float; maxY: float } with + member this.Contains px py : bool = + px >= this.minX && px <= this.maxX && + py >= this.minY && py <= this.maxY + + member this.IsSub otherRegion : bool = + this.minX >= otherRegion.minX && this.maxX <= otherRegion.maxX && + this.minY >= otherRegion.minY && this.maxY <= otherRegion.maxY + + member this.Intersects otherRegion : bool = + this.minX < otherRegion.maxX && this.maxX >= otherRegion.minX && + this.minY < otherRegion.maxY && this.maxY >= otherRegion.minY + +type Tree<'a when 'a :> I2DCoords> = + | Node of float * Tree<'a> * Tree<'a> + | Leaf of 'a + + static member buildTree (l: 'a list) : Tree<'a> = + let xSorted = List.toArray l + let ySorted = List.toArray l + Array.sortInPlaceWith cmpX xSorted + Array.sortInPlaceWith cmpY ySorted + + let rec buildTreeFromSortedArray (pXSorted: 'a[]) (pYSorted: 'a[]) (depth: int) : Tree<'a> = + if pXSorted.Length = 1 + then + Leaf pXSorted.[0] + else + if depth % 2 = 1 // 'depth' is odd -> vertical splitting else horizontal splitting. + then + let leftX, rightX = Array.splitAt ((pXSorted.Length + 1) / 2) pXSorted + let splitElement = Array.last leftX + let leftY, rightY = Array.partition (fun (e: 'a) -> cmpX e splitElement <= 0) pYSorted // FIXME: Maybe this operation can be optimized. + Node (splitElement.X, buildTreeFromSortedArray leftX leftY (depth + 1), buildTreeFromSortedArray rightX rightY (depth + 1)) + else + let downY, upY = Array.splitAt ((pYSorted.Length + 1) / 2) pYSorted + let splitElement = Array.last downY + let downX, upX = Array.partition (fun (e: 'a) -> cmpY e splitElement <= 0) pXSorted // FIXME: Maybe this operation can be optimized. + Node (splitElement.Y, buildTreeFromSortedArray downX downY (depth + 1), buildTreeFromSortedArray upX upY (depth + 1)) + + buildTreeFromSortedArray xSorted ySorted 1 + + static member search (tree: Tree<'a>) (searchRegion: Region) : 'a list = + let rec valuesFrom (tree: Tree<'a>) : 'a list = + match tree with + | Leaf v -> [v] + | Node (_, part1, part2) -> (valuesFrom part1) @ (valuesFrom part2) + + let rec searchWithRegion (tree: Tree<'a>) (currentRegion: Region) (depth: int) : 'a list = + match tree with + | Leaf v -> if searchRegion.Contains v.X v.Y then [v] else [] + | Node (splitValue, part1, part2) -> + let valuesInRegion (region: Region) (treeRegion: Tree<'a>) = + if region.IsSub searchRegion + then + valuesFrom treeRegion + elif region.Intersects searchRegion + then + searchWithRegion treeRegion region (depth + 1) + else + [] + + if depth % 2 = 1 // Vertical splitting. + then + let leftRegion = { currentRegion with maxX = splitValue } + let rightRegion = { currentRegion with minX = splitValue } + (valuesInRegion leftRegion part1) @ (valuesInRegion rightRegion part2) + else // Horizontal splitting. + let downRegion = { currentRegion with maxY = splitValue } + let upRegion = { currentRegion with minY = splitValue } + (valuesInRegion downRegion part1) @ (valuesInRegion upRegion part2) + + searchWithRegion tree { minX = Double.MinValue; maxX = Double.MaxValue; minY = Double.MinValue; maxY = Double.MaxValue } 1 + + +///// Tests. TODO: to put in a unit test. + +type Point (x: float, y: float) = + interface I2DCoords with + member this.X = x + member this.Y = y + + override this.ToString () = + sprintf "(%.1f, %.1f)" x y + +// TODO: test with identical X or Y coords +let test () = + let pts = [ + Point(1.0, 1.0) + Point(2.0, 2.0) + Point(1.5, 3.6) + Point(3.0, 3.2) + Point(4.0, 4.0) + Point(3.5, 1.5) + Point(2.5, 0.5) ] + + let tree = Tree.buildTree pts + Utils.dprintfn "Tree: %A" tree + + let s1 = Tree.search tree { minX = 0.0; maxX = 5.0; minY = 0.0; maxY = 5.0 } // All points. + Utils.dprintfn "s1: %A" s1 + + let s2 = Tree.search tree { minX = 2.8; maxX = 4.5; minY = 3.0; maxY = 4.5 } + Utils.dprintfn "s2: %A" s2 + + let s3 = Tree.search tree { minX = 2.0; maxX = 2.0; minY = 2.0; maxY = 2.0 } + Utils.dprintfn "s3: %A" s3 + +let test2 () = + let pts = [ + Point(1.0, 1.0) + Point(1.0, 2.0) + Point(1.0, 3.0) ] + + let tree = Tree.buildTree pts + Utils.dprintfn "Tree: %A" tree + + let s1 = Tree.search tree { minX = 1.0; maxX = 1.0; minY = 1.0; maxY = 1.0 } + Utils.dprintfn "s1: %A" s1 + + let s2 = Tree.search tree { minX = 1.0; maxX = 1.0; minY = 2.0; maxY = 2.0 } + Utils.dprintfn "s2: %A" s2 + + // This case result is wrong: FIXME + let s3 = Tree.search tree { minX = 1.0; maxX = 1.0; minY = 3.0; maxY = 3.0 } + Utils.dprintfn "s3: %A" s3 + + let s4 = Tree.search tree { minX = 0.0; maxX = 2.0; minY = 0.0; maxY = 4.0 } + Utils.dprintfn "s4: %A" s4 + diff --git a/Parasitemia/Parasitemia/MatchingEllipses.fs b/Parasitemia/Parasitemia/MatchingEllipses.fs index e4eae8d..afc4ab7 100644 --- a/Parasitemia/Parasitemia/MatchingEllipses.fs +++ b/Parasitemia/Parasitemia/MatchingEllipses.fs @@ -2,55 +2,95 @@ open System open System.Linq +open System.Collections open System.Collections.Generic open Types open Utils -type EllipseScore = { mutable matchingScore: float; e: Ellipse } +let matchingScoreThreshold1 = 0.6 +let matchingScoreThreshold2 = 1.0 -let matchingScoreThreshold1 = 0.5 -let matchingScoreThreshold2 = 2.0 +type EllipseScore (matchingScore: float, e: Ellipse) = + let mutable matchingScore = matchingScore -let ellipseArea e = e.a * e.b * Math.PI + member this.MatchingScore = matchingScore + member this.Ellipse = e + member val Processed = false with get, set + member val Removed = false with get, set -type MatchingEllipses () = + member this.AddMatchingScore(score: float) = + matchingScore <- matchingScore + score + + interface KdTree.I2DCoords with + member this.X = e.Cx + member this.Y = e.Cy + +type MatchingEllipses (radiusMin: float) = let ellipses = List() member this.Add (e: Ellipse) = - - // dprintfn "new ellipse: %A, nb: %A" e ellipses.Count - - let areaE = ellipseArea e - - let mutable matchingScoreE = 0.0 + ellipses.Add(EllipseScore(0.0, e)) - for other in ellipses do - let areaOther = ellipseArea other.e - let commonArea = EEOver.EEOverlapArea e other.e - let matchingScore = (2.0 * commonArea / (areaE + areaOther)) ** 2.0 - if matchingScore >= matchingScoreThreshold1 + member this.Ellipses : Ellipse list = + // 1) Create a kd-tree from the ellipses list. + let tree = KdTree.Tree.buildTree (List.ofSeq ellipses) + + // 2) Compute the matching score of each ellipses. + let windowSize = radiusMin + for e in ellipses do + e.Processed <- true + let areaE = e.Ellipse.Area + let window = { KdTree.minX = e.Ellipse.Cx - windowSize / 2.0 + KdTree.maxX = e.Ellipse.Cx + windowSize / 2.0 + KdTree.minY = e.Ellipse.Cy - windowSize / 2.0 + KdTree.maxY = e.Ellipse.Cy + windowSize / 2.0 } + for other in KdTree.Tree.search tree window do + if not other.Processed + then + let areaOther = other.Ellipse.Area + let commonArea = EEOver.EEOverlapArea e.Ellipse other.Ellipse + let matchingScore = 2.0 * commonArea / (areaE + areaOther) + if matchingScore >= matchingScoreThreshold1 + then + other.AddMatchingScore(matchingScore) + e.AddMatchingScore(matchingScore) + + // 3) Sort ellipses by their score. + ellipses.Sort(fun e1 e2 -> e2.MatchingScore.CompareTo(e1.MatchingScore)) + + // 4) Remove ellipses wich have a low score. + let i = ellipses.BinarySearch(EllipseScore(matchingScoreThreshold2, Ellipse(0.0, 0.0, 0.0, 0.0, 0.0)), + { new IComparer with + member this.Compare(e1, e2) = e2.MatchingScore.CompareTo(e1.MatchingScore) }) |> abs + let nbToRemove = ellipses.Count - i + if nbToRemove > 0 + then + for j in i .. nbToRemove - 1 do + ellipses.[j].Removed <- true + ellipses.RemoveRange(i, nbToRemove) + + // 5) Remove ellipses whose center is into an ellipse with a better score + for e in ellipses do + if not e.Removed then - other.matchingScore <- other.matchingScore + matchingScore - matchingScoreE <- matchingScoreE + matchingScore - printfn "Score" - - ellipses.Add({ matchingScore = matchingScoreE; e = e }) + let window = { KdTree.minX = e.Ellipse.Cx - e.Ellipse.A + KdTree.maxX = e.Ellipse.Cx + e.Ellipse.A + KdTree.minY = e.Ellipse.Cy - e.Ellipse.A + KdTree.maxY = e.Ellipse.Cy + e.Ellipse.A } + for other in KdTree.Tree.search tree window do + if not other.Removed && other.MatchingScore < e.MatchingScore + then + if e.Ellipse.Contains other.Ellipse.Cx other.Ellipse.Cy + then + other.Removed <- true + ellipses.RemoveAll(fun e -> e.Removed) |> ignore - member this.Ellipses : Ellipse list = dprintfn "Number of ellipse: %A" ellipses.Count - (*let sortedEllipses = - List.filter (fun e -> e.matchingScore >= matchingScoreThreshold2) (List.ofSeq ellipses) - |> List.sortWith (fun e1 e2 -> e2.matchingScore.CompareTo(e1.matchingScore))*) - - List.filter (fun e -> e.matchingScore >= matchingScoreThreshold2) (List.ofSeq ellipses) - |> List.sortWith (fun e1 e2 -> e2.matchingScore.CompareTo(e1.matchingScore)) - |> List.map (fun { e = e } -> e) + List.ofSeq ellipses |> List.map (fun e -> e.Ellipse) - // ellipses.Where(fun e -> e.matchingScore >= matchingScoreThreshold2) - // ([], fun acc { matchingScore = score; e = e } -> if score >= matchingScoreThreshold2 then e :: acc else acc) diff --git a/Parasitemia/Parasitemia/Parasitemia.fsproj b/Parasitemia/Parasitemia/Parasitemia.fsproj index d3a8fe4..388b059 100644 --- a/Parasitemia/Parasitemia/Parasitemia.fsproj +++ b/Parasitemia/Parasitemia/Parasitemia.fsproj @@ -9,7 +9,7 @@ WinExe Parasitemia Parasitemia - v4.6 + v4.6.1 true 4.4.0.0 Parasitemia @@ -67,9 +67,11 @@ + + @@ -82,7 +84,8 @@ ..\..\..\Emgu\emgucv-windows-universal 3.0.0.2157\bin\Emgu.Util.dll - + + ..\packages\FSharp.ViewModule.Core.0.9.9.1\lib\net45\FSharp.ViewModule.Core.Wpf.dll True @@ -99,16 +102,21 @@ ..\packages\log4net.1.2.10\lib\2.0\log4net.dll True - - + + ..\packages\MathNet.Numerics.3.9.0\lib\net40\MathNet.Numerics.dll True + + ..\packages\MathNet.Numerics.FSharp.3.9.0\lib\net40\MathNet.Numerics.FSharp.dll + True + + - + diff --git a/Parasitemia/Parasitemia/ParasitesMarker.fs b/Parasitemia/Parasitemia/ParasitesMarker.fs new file mode 100644 index 0000000..44da901 --- /dev/null +++ b/Parasitemia/Parasitemia/ParasitesMarker.fs @@ -0,0 +1,39 @@ +module ParasitesMarker + +open System.Drawing + +open Emgu.CV +open Emgu.CV.Structure + +type Result = { + darkStain: Image + stain: Image + infection: Image } + +let find (green: Image) (filteredGreen: Image) (config: Config.Config) : Result = + + // We use the filtered image to find the dark stain. + let { KMedians.fg = fg; KMedians.median_bg = median_bg; KMedians.median_fg = median_fg; KMedians.d_fg = d_fg } = KMedians.kmedians filteredGreen 1.0 + let darkStain = d_fg.Cmp(median_bg * config.darkStainLevel, CvEnum.CmpType.GreaterThan) + darkStain._And(filteredGreen.Cmp(median_fg, CvEnum.CmpType.LessThan)) + darkStain._And(fg) + + let fgFloat = (fg / 255.0).Convert() + let greenWithoutBg = green.Copy() + greenWithoutBg.SetValue(Gray(0.0), fg.Not()) + + let findSmears (sigma: float) (level: float) : Image = + let greenWithoutBgSmoothed = ImgTools.gaussianFilter greenWithoutBg sigma + let fgSmoothed = ImgTools.gaussianFilter fgFloat sigma + + let smears = (greenWithoutBg.Mul(fgSmoothed)).Cmp(greenWithoutBgSmoothed.Mul(level), CvEnum.CmpType.LessThan) + smears._And(fg) + smears + + { darkStain = darkStain; + stain = findSmears config.stainSigma config.stainLevel + infection = findSmears config.infectionSigma config.infectionLevel } + + + + diff --git a/Parasitemia/Parasitemia/Program.fs b/Parasitemia/Parasitemia/Program.fs index 56672bd..fd9a2f9 100644 --- a/Parasitemia/Parasitemia/Program.fs +++ b/Parasitemia/Parasitemia/Program.fs @@ -33,20 +33,39 @@ do Utils.log <- (fun m -> log mainWindow m) let config = { + scale = 1.0 + doGSigma1 = 1.5 doGSigma2 = 20.0 doGLowFreqPercentageReduction = 0.9 - scale = 1.0 - } + + darkStainLevel = 0.839 + + stainSigma = 10.0 + stainLevel = 0.9 + stainSpreadRequired = 3.0 + + infectionSigma = 2.0 + infectionLevel = 0.762 + infectionPixelsRequired = 1 } - //use img = new Image("../../../../src/Tests_hough/images/1508133543-7-4-120-0001.png") + ///// ELLIPSES ///// + //use img = new Image("../../../../src/Tests_hough/images/1508133543-7-4-120-0001.png") //use img = new Image("../../../../src/Tests_hough/images/rbc_single.png") - //use img = new Image("../../../../src/Tests_hough/images/rbc_single_oblong_2.png") - use img = new Image("../../../../src/Tests_hough/images/two_rbc_1.png") - - let result = ImageAnalysis.doAnalysis img config + //use img = new Image("../../../../src/Tests_hough/images/rbc_single_oblong_4.png") + //use img = new Image("../../../../src/Tests_hough/images/strange_rbc_1.png") + //use img = new Image("../../../../src/Tests_hough/images/rbc_single_blurred.png") + //use img = new Image("../../../../src/Tests_hough/images/lot.png") + + + ///// PARASITES ///// + use img = new Image("../../../../src/Parasites/images/1.png") + +// KdTree.test3 () +// Utils.dprintfn "area: %A" area + let result = ImageAnalysis.doAnalysis img config + display mainWindow result.img mainWindow.Root.Show() - app.Run() |> ignore diff --git a/Parasitemia/Parasitemia/Types.fs b/Parasitemia/Parasitemia/Types.fs index 308fc6b..589a735 100644 --- a/Parasitemia/Parasitemia/Types.fs +++ b/Parasitemia/Parasitemia/Types.fs @@ -2,7 +2,14 @@ open System -type Ellipse = { cx: float; cy: float; a: float; b: float; alpha: float } - +type Ellipse (cx: float, cy: float, a: float, b: float, alpha: float) = + member this.Cx = cx + member this.Cy = cy + member this.A = a + member this.B = b + member this.Alpha = alpha + member this.Area = a * b * Math.PI -//type PointImg = { x: int; y: int } \ No newline at end of file + // Does the ellipse contain the point (x, y)?. + member this.Contains x y = + ((x - cx) * cos alpha + (y - cy) * sin alpha) ** 2.0 / a ** 2.0 + ((x - cx) * sin alpha - (y - cy) * cos alpha) ** 2.0 / b ** 2.0 <= 1.0 \ No newline at end of file diff --git a/Parasitemia/Parasitemia/packages.config b/Parasitemia/Parasitemia/packages.config index 1e68f7b..20de8bc 100644 --- a/Parasitemia/Parasitemia/packages.config +++ b/Parasitemia/Parasitemia/packages.config @@ -2,7 +2,10 @@ + + + \ No newline at end of file -- 2.43.0