From e76da913cd58078ad2479357b2430ed62a6e0777 Mon Sep 17 00:00:00 2001 From: Greg Burri Date: Mon, 14 Dec 2015 02:29:31 +0100 Subject: [PATCH] The main process is now complete. --- Parasitemia/Parasitemia/Classifier.fs | 160 ++++++++++++++++++-- Parasitemia/Parasitemia/Config.fs | 18 ++- Parasitemia/Parasitemia/EEOver.fs | 8 +- Parasitemia/Parasitemia/Ellipse.fs | 15 +- Parasitemia/Parasitemia/ImageAnalysis.fs | 138 ----------------- Parasitemia/Parasitemia/ImgTools.fs | 130 +++++++++++----- Parasitemia/Parasitemia/KMedians.fs | 10 +- Parasitemia/Parasitemia/KdTree.fs | 24 +-- Parasitemia/Parasitemia/MainAnalysis.fs | 97 ++++++++++++ Parasitemia/Parasitemia/MatchingEllipses.fs | 132 ++++++++-------- Parasitemia/Parasitemia/Parasitemia.fsproj | 8 +- Parasitemia/Parasitemia/ParasitesMarker.fs | 14 +- Parasitemia/Parasitemia/Program.fs | 149 ++++++++++++------ Parasitemia/Parasitemia/Types.fs | 40 ++++- Parasitemia/Parasitemia/Utils.fs | 27 +++- 15 files changed, 629 insertions(+), 341 deletions(-) delete mode 100644 Parasitemia/Parasitemia/ImageAnalysis.fs create mode 100644 Parasitemia/Parasitemia/MainAnalysis.fs diff --git a/Parasitemia/Parasitemia/Classifier.fs b/Parasitemia/Parasitemia/Classifier.fs index 24ea4a7..37507c5 100644 --- a/Parasitemia/Parasitemia/Classifier.fs +++ b/Parasitemia/Parasitemia/Classifier.fs @@ -1,31 +1,163 @@ module Classifier open System +open System.Collections.Generic open System.Drawing open Emgu.CV open Emgu.CV.Structure +open Types +open Utils -type CellClass = HealthyRBC | InfectedRBC | Peculiar -type Cell = { - cellClass: CellClass - center: Point - elements: Matrix } +type private EllipseFlaggedKd (e: Ellipse) = + inherit Ellipse (e.Cx, e.Cy, e.A, e.B, e.Alpha) -type KdEllipse (e: Types.Ellipse) = - inherit Types.Ellipse (e.Cx, e.Cy, e.A, e.B, e.Alpha) - - interface KdTree.I2DCoords with + member val Removed = false with get, set + + interface KdTree.I2DCoords with member this.X = this.Cx member this.Y = this.Cy - - -let findCells (ellipses: Types.Ellipse list) (parasites: ParasitesMarker.Result) (fg: Image) : Cell list = + + +let findCells (ellipses: Ellipse list) (parasites: ParasitesMarker.Result) (fg: Image) (config: Config.Config) : Cell list = if ellipses.IsEmpty then [] else - let tree = KdTree.Tree.buildTree (ellipses |> List.map KdEllipse) - [] \ No newline at end of file + let infection = parasites.infection.Copy() // To avoid to modify the parameter. + + // This is the minimum window size to check if other ellipses touch 'e'. + let searchRegion (e: Ellipse) = { KdTree.minX = e.Cx - (e.A + config.maxRBCSize) * config.scale + KdTree.maxX = e.Cx + (e.A + config.maxRBCSize) * config.scale + KdTree.minY = e.Cy - (e.A + config.maxRBCSize) * config.scale + KdTree.maxY = e.Cy + (e.A + config.maxRBCSize) * config.scale } + + // The minimum window to contain a given ellipse. + let ellipseWindow (e: Ellipse) = + let cx, cy = roundInt e.Cx, roundInt e.Cy + let a = int (e.A + 0.5) + cx - a, cy - a, cx + a, cy + a + + // 1) Remove ellipses touching the edges. + let ellipsesInside = ellipses |> List.map (fun e -> + EllipseFlaggedKd (e, Removed = e.isOutside (float fg.Width) (float fg.Height))) + + // 2) Associate touching ellipses with each ellipses. + let tree = KdTree.Tree.BuildTree ellipsesInside + let neighbors (e: Ellipse) : (EllipseFlaggedKd * PointD * PointD) list = + tree.Search (searchRegion e) + // We only keep the ellipses touching 'e'. + |> List.choose (fun otherE -> + match EEOver.EEOverlapArea e otherE with + | Some (area, px, py) when area > 0.0 && px.Length >= 2 && py.Length >= 2 -> + Some (otherE, PointD(px.[0], py.[0]), PointD(px.[1], py.[1])) + | _ -> + None ) + + let ellipsesWithNeigbors = ellipsesInside |> List.choose (fun e -> if e.Removed then None else Some (e, neighbors e)) + + // 3) Remove ellipse with a lower percentage of foreground. + let fgData = fg.Data + for e, neighbors in ellipsesWithNeigbors do + let minX, minY, maxX, maxY = ellipseWindow e + + let mutable totalElement = 0 + let mutable fgElement = 0 + + for y in minY .. maxY do + for x in minX .. maxX do + let yf, xf = float y, float x + if e.Contains xf yf && neighbors |> List.forall (fun (otherE, _, _) -> not <| otherE.Contains xf yf) + then + totalElement <- totalElement + 1 + if fgData.[y, x, 0] > 0uy + then + fgElement <- fgElement + 1 + + if totalElement < config.minimumCellArea || (float fgElement) / (float totalElement) < config.percentageOfFgValidCell + then + e.Removed <- true + + // 3) Define pixels associated to each ellipse and create the cells. + + // Return 'true' if the point 'p' is owned by e. + // The lines represents all intersections with other ellipses. + let pixelOwnedByE (p: PointD) (e: Ellipse) (lines: Line list) = + e.Contains p.X p.Y && + seq { + let c = PointD(e.Cx, e.Cy) + for d1 in lines do + let d2 = Utils.lineFromTwoPoints c p + let p' = Utils.pointFromTwoLines d1 d2 + yield sign (c.X - p.X) <> sign (c.X - p'.X) || Utils.squaredDistanceTwoPoints c p' > Utils.squaredDistanceTwoPoints c p // 'false' -> the point is owned by another ellipse. + } |> Seq.forall id + + + ellipsesWithNeigbors + |> List.choose (fun (e, neighbors) -> + if e.Removed + then + None + else + let minX, minY, maxX, maxY = ellipseWindow e + + let infectedPixels = List() + let mutable stainPixels = 0 + let mutable darkStainPixels = 0 + let mutable nbElement = 0 + + let mutable stain_x = 0.0 + let mutable stain_x2 = 0.0 + let mutable stain_y = 0.0 + let mutable stain_y2 = 0.0 + + let mutable sumCoords_y = 0 + let mutable sumCoords_x = 0 + + let elements = new Matrix(maxY - minY + 1, maxX - minX + 1) + for y in minY .. maxY do + for x in minX .. maxX do + let p = PointD(float x, float y) + if pixelOwnedByE p e (neighbors |> List.choose (fun (otherE, p1, p2) -> if otherE.Removed then None else Some (Utils.lineFromTwoPoints p1 p2))) && + fg.Data.[y, x, 0] > 0uy + then + elements.[y-minY, x-minX] <- 1uy + nbElement <- nbElement + 1 + sumCoords_y <- sumCoords_y + y + sumCoords_x <- sumCoords_x + x + + if infection.Data.[y, x, 0] > 0uy + then + infectedPixels.Add(Point(x, y)) + + if parasites.stain.Data.[y, x, 0] > 0uy + then + stainPixels <- stainPixels + 1 + stain_x <- stain_x + (float x) + stain_x2 <- stain_x2 + (float x) ** 2.0 + stain_y <- stain_y + (float y) + stain_y2 <- stain_y2 + (float y) ** 2.0 + + if parasites.darkStain.Data.[y, x, 0] > 0uy + then + darkStainPixels <- darkStainPixels + 1 + + let cellClass = + if float darkStainPixels > config.MaxDarkStainRatio * (float nbElement) (* || + sqrt (((float sumCoords_x) / (float nbElement) - e.Cx) ** 2.0 + ((float sumCoords_y) / (float nbElement) - e.Cy) ** 2.0) > e.A * config.maxOffcenter *) + then + Peculiar + elif infectedPixels.Count > config.infectionPixelsRequired + then + let infectionToRemove = ImgTools.connectedComponents parasites.stain infectedPixels + for p in infectionToRemove do + infection.Data.[p.Y, p.X, 0] <- 0uy + InfectedRBC + else + HealthyRBC + + Some { cellClass = cellClass + center = Point(roundInt e.Cx, roundInt e.Cy) + elements = elements }) diff --git a/Parasitemia/Parasitemia/Config.fs b/Parasitemia/Parasitemia/Config.fs index 75aa95f..8a172b8 100644 --- a/Parasitemia/Parasitemia/Config.fs +++ b/Parasitemia/Parasitemia/Config.fs @@ -1,8 +1,17 @@ module Config +type Debug = + | DebugOff + | DebugOn of string // Output directory. + type Config = { + debug: Debug + scale: float + minRBCSize: float + maxRBCSize: float + doGSigma1: float doGSigma2: float doGLowFreqPercentageReduction: float @@ -13,8 +22,15 @@ type Config = { stainSigma: float stainLevel: float stainSpreadRequired: float - + infectionSigma: float infectionLevel: float infectionPixelsRequired: int + + percentageOfFgValidCell: float + + MaxDarkStainRatio: float + + minimumCellArea: int + maxOffcenter: float } \ No newline at end of file diff --git a/Parasitemia/Parasitemia/EEOver.fs b/Parasitemia/Parasitemia/EEOver.fs index 640d752..52b5963 100644 --- a/Parasitemia/Parasitemia/EEOver.fs +++ b/Parasitemia/Parasitemia/EEOver.fs @@ -718,4 +718,10 @@ let EEOverlapArea (e1: Types.Ellipse) (e2: Types.Ellipse) : (float * float[] * f | _ -> -1.0 if nintpts = 0 then Some (area, [||], [||]) - else Some (area, xint.[..nintpts-1], yint.[..nintpts-1]) \ No newline at end of file + else + let xTransform = Array.zeroCreate nintpts + let yTransform = Array.zeroCreate nintpts + for i in 0 .. (nintpts - 1) do + xTransform.[i] <- cos phi_1 * xint.[i] - sin phi_1 * yint.[i] + h1 + yTransform.[i] <- sin phi_1 * xint.[i] + cos phi_1 * yint.[i] + k1 + Some (area, xTransform, yTransform) \ No newline at end of file diff --git a/Parasitemia/Parasitemia/Ellipse.fs b/Parasitemia/Parasitemia/Ellipse.fs index 5bfcc33..62c0936 100644 --- a/Parasitemia/Parasitemia/Ellipse.fs +++ b/Parasitemia/Parasitemia/Ellipse.fs @@ -42,6 +42,8 @@ let private goldenSectionSearch (f: float -> float) (nbIter: int) (xmin: float) let x = (b + a) / 2.0 x, f x +// Ellipse.A is always equal or greater than Ellipse.B. +// Ellipse.Alpha is between 0 and Pi. 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 @@ -57,7 +59,7 @@ let ellipse (p1x: float) (p1y: float) (m1: float) (p2x: float) (p2y: float) (m2: let alpha1 = atan m1 let alpha2 = atan m2 - let r1 = sqrt (p1x**2.0 + p1y**2.0) + let r1 = sqrt (p1x ** 2.0 + p1y ** 2.0) let theta1 = atan2 p1y p1x let r2 = sqrt (p2x**2.0 + p2y**2.0) @@ -66,13 +68,13 @@ let ellipse (p1x: float) (p1y: float) (m1: float) (p2x: float) (p2y: float) (m2: let valid = 4.0 * sin (alpha1 - theta1) * (-r1 * sin (alpha1 - theta1) + r2 * sin (alpha1 - theta2)) * sin (alpha2 - theta2) * (-r1 * sin (alpha2 - theta1) + r2 * sin (alpha2 - theta2)) + - r1 * r2 * sin (alpha1 - alpha2) **2.0 * sin (theta1 - theta2) **2.0 < 0.0 + r1 * r2 * sin (alpha1 - alpha2) ** 2.0 * sin (theta1 - theta2) ** 2.0 < 0.0 if valid then let r theta = (r1 * r2 * (r1 * (cos (alpha2 + theta - theta1 - theta2) - cos (alpha2 - theta) * cos (theta1 - theta2)) * sin (alpha1 - theta1) + r2 * (-cos (alpha1 + theta - theta1 - theta2) + cos (alpha1 - theta) * cos (theta1 - theta2)) * sin (alpha2 - theta2)) * sin (theta1 - theta2)) / - (sin (alpha1 - theta1) * sin (alpha2 - theta2) * (r1 * sin (theta - theta1) - r2 * sin (theta - theta2))**2.0 - r1 * r2 * sin (alpha1 - theta) * sin (alpha2 - theta) * sin (theta1 - theta2)**2.0) + (sin (alpha1 - theta1) * sin (alpha2 - theta2) * (r1 * sin (theta - theta1) - r2 * sin (theta - theta2)) ** 2.0 - r1 * r2 * sin (alpha1 - theta) * sin (alpha2 - theta) * sin (theta1 - theta2) ** 2.0) let rabs = r >> abs @@ -118,10 +120,10 @@ let ellipse (p1x: float) (p1y: float) (m1: float) (p2x: float) (p2y: float) (m2: let rellipse theta = sqrt ( - rc**2.0 + (r1**2.0 * r2**2.0 * (r1 * (cos (alpha2 + theta - theta1 - theta2) - cos (alpha2 - theta) * cos (theta1 - theta2)) * sin (alpha1 - theta1) + r2 * (-cos (alpha1 + theta - theta1 - theta2) + cos (alpha1 - theta) * cos (theta1 - theta2)) * sin (alpha2 - theta2))**2.0 * sin (theta1 - theta2)**2.0) / - (sin (alpha1 - theta1) * sin (alpha2 - theta2) * (r1 * sin (theta - theta1) - r2 * sin (theta - theta2))**2.0 - r1 * r2 * sin (alpha1 - theta) * sin (alpha2 - theta) * sin (theta1 - theta2)**2.0)**2.0 - + rc ** 2.0 + (r1 ** 2.0 * r2 ** 2.0 * (r1 * (cos (alpha2 + theta - theta1 - theta2) - cos (alpha2 - theta) * cos (theta1 - theta2)) * sin (alpha1 - theta1) + r2 * (-cos (alpha1 + theta - theta1 - theta2) + cos (alpha1 - theta) * cos (theta1 - theta2)) * sin (alpha2 - theta2)) ** 2.0 * sin (theta1 - theta2) ** 2.0) / + (sin (alpha1 - theta1) * sin (alpha2 - theta2) * (r1 * sin (theta - theta1) - r2 * sin (theta - theta2)) ** 2.0 - r1 * r2 * sin (alpha1 - theta) * sin (alpha2 - theta) * sin (theta1 - theta2) ** 2.0) ** 2.0 - (2.0 * r1 * r2 * rc * cos (theta - psi) * (r1 * (cos (alpha2 + theta - theta1 - theta2) - cos (alpha2 - theta) * cos (theta1 - theta2)) * sin (alpha1 - theta1) + r2 * (-cos (alpha1 + theta - theta1 - theta2) + cos (alpha1 - theta) * cos (theta1 - theta2)) * sin (alpha2 - theta2)) * sin (theta1 - theta2)) / - (sin (alpha1 - theta1) * sin (alpha2 - theta2) * (r1 * sin (theta - theta1) - r2 * sin (theta - theta2))**2.0 - r1 * r2 * sin (alpha1 - theta) * sin (alpha2 - theta) * sin (theta1 - theta2)**2.0)) + (sin (alpha1 - theta1) * sin (alpha2 - theta2) * (r1 * sin (theta - theta1) - r2 * sin (theta - theta2)) ** 2.0 - r1 * r2 * sin (alpha1 - theta) * sin (alpha2 - theta) * sin (theta1 - theta2) ** 2.0)) // We search for an interval [theta_a, theta_b] and assume the function is unimodal in this interval. let r1eTheta, r1e = goldenSectionSearch rellipse accuracy_extremum_search_2 0.0 (Math.PI / 2.0) Maximum // Pi/2 and not pi because the period is Pi. @@ -187,7 +189,6 @@ let private areVectorsValid (p1x: float) (p1y: float) (p2x: float) (p2y: float) else Some (m1, m2) - let find (edges: Matrix) (xDir: Image) (yDir: Image) diff --git a/Parasitemia/Parasitemia/ImageAnalysis.fs b/Parasitemia/Parasitemia/ImageAnalysis.fs deleted file mode 100644 index feb427d..0000000 --- a/Parasitemia/Parasitemia/ImageAnalysis.fs +++ /dev/null @@ -1,138 +0,0 @@ -module ImageAnalysis - -open System -open System.Drawing - -open Emgu.CV -open Emgu.CV.Structure - -open Utils -open ImgTools -open Config -open Types - -(*type Result = { - RBCPositions : Point list - infectedRBCPositions : Point list - img: Image -}*) - -let doAnalysis (img: Image) (config: Config) : Classifier.Cell list = - - let imgFloat = img.Convert() - use scaledImg = if config.scale = 1.0 then imgFloat else imgFloat.Resize(config.scale, CvEnum.Inter.Area) - - (*use scaledImg = - if config.scale = 1.0 - then - img - else - let m = new Mat() - CvInvoke.Resize(img, m, Size(roundInt (float img.Size.Width * config.scale), roundInt (float img.Size.Height * config.scale))) - m*) - - use green = scaledImg.Item(1) - - //use green = new Matrix(scaledImg.Size) - //CvInvoke.MixChannels(scaledImg, green, [| 1; 0 |]) - - //let greenMatrix = new Matrix(green.Height, green.Width, green.DataPointer) - - //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 ] - [ 1.0f; 0.0f; -1.0f ]], Point(0, 0)) - - use xEdges = filteredGreen.Convolution(sobelKernel).Convert() - use yEdges = filteredGreen.Convolution(sobelKernel.Transpose()).Convert() - - let xEdgesData = xEdges.Data - let yEdgesData = yEdges.Data - for r in 0..xEdges.Rows-1 do - xEdgesData.[r, 0, 0] <- 0.0 - xEdgesData.[r, xEdges.Cols-1, 0] <- 0.0 - yEdgesData.[r, 0, 0] <- 0.0 - yEdgesData.[r, xEdges.Cols-1, 0] <- 0.0 - - for c in 0..xEdges.Cols-1 do - xEdgesData.[0, c, 0] <- 0.0 - xEdgesData.[xEdges.Rows-1, c, 0] <- 0.0 - yEdgesData.[0, c, 0] <- 0.0 - yEdgesData.[xEdges.Rows-1, c, 0] <- 0.0 - - use magnitudes = new Matrix(xEdges.Size) - CvInvoke.CartToPolar(xEdges, yEdges, magnitudes, new Mat()) // Compute the magnitudes (without angles). - - let min = ref 0.0 - let minLocation = ref <| Point() - let max = ref 0.0 - let maxLocation = ref <| Point() - magnitudes.MinMax(min, max, minLocation, maxLocation) - - use magnitudesByte = ((magnitudes / !max) * 255.0).Convert() // Otsu from OpenCV only support 'byte'. - use edges = new Matrix(xEdges.Size) - CvInvoke.Threshold(magnitudesByte, edges, 0.0, 1.0, CvEnum.ThresholdType.Otsu ||| CvEnum.ThresholdType.Binary) |> ignore - - logTime "Finding edges" (fun() -> - thin edges) - - logTime "Removing small connected components" (fun () -> - removeArea edges 12) - - saveMat (edges * 255.0) "edges.png" - - - let kmediansResults = KMedians.kmedians filteredGreen 1.0 - - let parasites = ParasitesMarker.find green filteredGreen kmediansResults config - - saveImg parasites.darkStain "parasites_dark_stain.png" - saveImg parasites.stain "parasites_stain.png" - saveImg parasites.infection "parasites_infection.png" - - let radiusRange = config.scale * 20.0, config.scale * 40.0 - let windowSize = roundInt (1.6 * (snd radiusRange)) - let factorNbPick = 1.5 - let ellipses = logTime "Finding ellipses" (fun () -> - Ellipse.find edges xEdges yEdges radiusRange windowSize factorNbPick) |> List.filter (fun e -> not (e.CutAVericalLine 0.0) && - not (e.CutAVericalLine (float img.Width)) && - not (e.CutAnHorizontalLine 0.0) && - not (e.CutAnHorizontalLine (float img.Height))) - - drawEllipses img ellipses (Bgr(0.0, 255.0, 255.0)) - //saveImg img "ellipses.png" - - Classifier.findCells ellipses parasites kmediansResults.fg - - // - - (*use imageHSV = scaledImage.Convert() - let H, S = match imageHSV.Split() with // Warning: H is from 0 to 179°. - | [| H; S; _|] -> H, S - | _ -> failwith "unable to split the HSV channels" - - let hueShiftValue = 175 - // Modulo operator doesn't exist on matrix thus we have to apply a function to every pixels. - let correctedH : Image = H.Convert(fun b _ _ -> - (255 - int(b) * 255 / 179 + hueShiftValue) % 256 |> byte - ) - - let correctedS : Image = S.Not() - - let filteredH = correctedH.SmoothMedian(5) - let filteredS = correctedS.SmoothMedian(5)*) - - //let greenChannel = scaledImage.Item(1) - - //let filteredImage = (gaussianFilter greenChannel config.doGSigma1) - config.doGLowFreqPercentageReduction * (gaussianFilter greenChannel config.doGSigma2) - - // let filteredImage = greenChannel.ThresholdAdaptive(Gray(255.), CvEnum.AdaptiveThresholdType.GaussianC, CvEnum.ThresholdType.Binary, 61, Gray(5.0)) - // let thresholdedImage = filteredImage.CopyBlank() - - // CvInvoke.Threshold(filteredImage, thresholdedImage, 0., 255., CvEnum.ThresholdType.Otsu ||| CvEnum.ThresholdType.BinaryInv) |> ignore - - // filteredImage <| \ No newline at end of file diff --git a/Parasitemia/Parasitemia/ImgTools.fs b/Parasitemia/Parasitemia/ImgTools.fs index 5278313..c34b9c7 100644 --- a/Parasitemia/Parasitemia/ImgTools.fs +++ b/Parasitemia/Parasitemia/ImgTools.fs @@ -10,7 +10,7 @@ open Emgu.CV.Structure open Utils // Normalize image values between 0uy and 255uy. -let normalizeAndConvert (img: Image) : Image = +let normalizeAndConvert (img: Image) : Image = let min = ref [| 0.0 |] let minLocation = ref <| [| Point() |] let max = ref [| 0.0 |] @@ -22,10 +22,10 @@ let gaussianFilter (img : Image<'TColor, 'TDepth>) (standardDeviation : float) : let size = 2 * int (ceil (4.0 * standardDeviation)) + 1 img.SmoothGaussian(size, size, standardDeviation, standardDeviation) -// Zhang and Suen algorithm. +// Zhang and Suen algorithm. // Modify 'mat' in place. -let thin (mat: Matrix) = - let neighbors = [| +let thin (mat: Matrix) = + let neighbors = [| (-1, 0) // p2 (-1, 1) // p3 ( 0, 1) // p4 @@ -41,15 +41,15 @@ let thin (mat: Matrix) = let mutable data2 = Array2D.zeroCreate h w // Return the list of neighbor values. - let neighborsValues (p1i, p1j) = - Array.map (fun (ni, nj) -> + let neighborsValues (p1i, p1j) = + Array.map (fun (ni, nj) -> let pi = p1i + ni let pj = p1j + nj if pi < 0 || pi >= h || pj < 0 || pj >= w then 0uy else data1.[pi, pj] ) neighbors // Return the number of 01 pattern in 'values' in a circular way. - let pattern01 (values: byte[]) = + let pattern01 (values: byte[]) = let mutable nb = 0 let mutable lastValue = 255uy for v in values do @@ -61,7 +61,7 @@ let thin (mat: Matrix) = then nb <- nb + 1 nb - + let mutable pixelChanged = true let mutable oddIteration = true while pixelChanged do @@ -84,16 +84,22 @@ let thin (mat: Matrix) = else data2.[i, j] <- 0uy - oddIteration <- not oddIteration + oddIteration <- not oddIteration let tmp = data1 data1 <- data2 data2 <- tmp - + + + +let pop (l: List<'a>) : 'a = + let n = l.[l.Count - 1] + l.RemoveAt(l.Count - 1) + n // Remove all 8-connected pixels with an area equal or greater than 'areaSize'. // Modify 'mat' in place. let removeArea (mat: Matrix) (areaSize: int) = - let neighbors = [| + let neighbors = [| (-1, 0) // p2 (-1, 1) // p3 ( 0, 1) // p4 @@ -102,29 +108,24 @@ let removeArea (mat: Matrix) (areaSize: int) = ( 1, -1) // p7 ( 0, -1) // p8 (-1, -1) |] // p9 - + let mat' = new Matrix(mat.Size) let w = mat'.Width let h = mat'.Height - mat.CopyTo(mat') - + mat.CopyTo(mat') + let data = mat.Data let data' = mat'.Data for i in 0..h-1 do for j in 0..w-1 do - if data'.[i, j] = 1uy + if data'.[i, j] = 1uy then let neighborhood = List<(int*int)>() let neighborsToCheck = List<(int*int)>() - neighborsToCheck.Add((i, j)) + neighborsToCheck.Add((i, j)) data'.[i, j] <- 0uy - let pop (l: List<'a>) : 'a = - let n = l.[l.Count - 1] - l.RemoveAt(l.Count - 1) - n - while neighborsToCheck.Count > 0 do let (ci, cj) = pop neighborsToCheck neighborhood.Add((ci, cj)) @@ -140,36 +141,61 @@ let removeArea (mat: Matrix) (areaSize: int) = for (ni, nj) in neighborhood do data.[ni, nj] <- 0uy +let connectedComponents (img: Image) (startPoints: List) : List = + let w = img.Width + let h = img.Height + + let pointChecked = HashSet() + let pointToCheck = List(startPoints); + + let data = img.Data + + while pointToCheck.Count > 0 do + let next = pop pointToCheck + pointChecked.Add(next) |> ignore + for ny in -1 .. 1 do + for nx in -1 .. 1 do + if ny <> 0 && nx <> 0 + then + let p = Point(next.X + nx, next.Y + ny) + if p.X >= 0 && p.X < w && p.Y >= 0 && p.Y < h && data.[p.Y, p.X, 0] > 0uy && not (pointChecked.Contains p) + then + pointToCheck.Add(p) + + List(pointChecked) + -let saveImg (img: Image<'TColor, 'TDepth>) (name: string) = - img.Save("output/" + name) - +let saveImg (img: Image<'TColor, 'TDepth>) (filepath: string) = + img.Save(filepath) -let saveMat (mat: Matrix<'TDepth>) (name: string) = + +let saveMat (mat: Matrix<'TDepth>) (filepath: string) = use img = new Image(mat.Size) mat.CopyTo(img) - saveImg img name - + saveImg img filepath + (*let drawEllipse (img: Image<'TColor, 'TDepth>) (e: Types.Ellipse) (color: 'TColor) = let e' = Ellipse(PointF(float32 e.cx, float32 e.cy), SizeF(2.0f * float32 e.a, 2.0f * float32 e.b), float32 e.alpha) img.Draw(e', color)*) -let drawLine (img: Image<'TColor, 'TDepth>) (color: 'TColor) (x0: float) (y0: float) (x1: float) (y1: float) = - let x0, y0, x1, y1 = roundInt(x0), roundInt(y0), roundInt(x1), roundInt(y1) - +let drawLine (img: Image<'TColor, 'TDepth>) (color: 'TColor) (x0: int) (y0: int) (x1: int) (y1: int) = img.Draw(LineSegment2D(Point(x0, y0), Point(x1, y1)), color, 1); -let drawEllipse (img: Image<'TColor, 'TDepth>) (e: Types.Ellipse) (color: 'TColor) = +let drawLineF (img: Image<'TColor, 'TDepth>) (color: 'TColor) (x0: float) (y0: float) (x1: float) (y1: float) = + let x0, y0, x1, y1 = roundInt(x0), roundInt(y0), roundInt(x1), roundInt(y1) + drawLine img color x0 y0 x1 y1 + +let drawEllipse (img: Image<'TColor, 'TDepth>) (e: Types.Ellipse) (color: 'TColor) = let cosAlpha = cos e.Alpha let sinAlpha = sin e.Alpha - + let mutable x0 = 0.0 let mutable y0 = 0.0 let mutable first_iteration = true - + let n = 40 let thetaIncrement = 2.0 * Math.PI / (float n) - + for theta in 0.0 .. thetaIncrement .. 2.0 * Math.PI do let cosTheta = cos theta let sinTheta = sin theta @@ -178,7 +204,7 @@ let drawEllipse (img: Image<'TColor, 'TDepth>) (e: Types.Ellipse) (color: 'TColo if not first_iteration then - drawLine img color x0 y0 x y + drawLineF img color x0 y0 x y else first_iteration <- false @@ -186,4 +212,36 @@ let drawEllipse (img: Image<'TColor, 'TDepth>) (e: Types.Ellipse) (color: 'TColo 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 + List.iter (fun e -> drawEllipse img e color) ellipses + + +let rngCell = System.Random() +let drawCell (img: Image) (drawCellContent: bool) (c: Types.Cell) = + if drawCellContent + then + let colorB = rngCell.Next(20, 70) + let colorG = rngCell.Next(20, 70) + let colorR = rngCell.Next(20, 70) + + for y in 0 .. c.elements.Height - 1 do + for x in 0 .. c.elements.Width - 1 do + if c.elements.[y, x] > 0uy + then + let dx, dy = c.center.X - c.elements.Width / 2, c.center.Y - c.elements.Height / 2 + let b = img.Data.[y + dy, x + dx, 0] |> int + let g = img.Data.[y + dy, x + dx, 1] |> int + let r = img.Data.[y + dy, x + dx, 2] |> int + img.Data.[y + dy, x + dx, 0] <- if b + colorB > 255 then 255uy else byte (b + colorB) + img.Data.[y + dy, x + dx, 1] <- if g + colorG > 255 then 255uy else byte (g + colorG) + img.Data.[y + dy, x + dx, 2] <- if r + colorR > 255 then 255uy else byte (r + colorR) + + let crossColor = match c.cellClass with + | Types.HealthyRBC -> Bgr(255.0, 0.0, 0.0) + | Types.InfectedRBC -> Bgr(0.0, 0.0, 255.0) + | Types.Peculiar -> Bgr(0.0, 0.0, 0.0) + + drawLine img crossColor (c.center.X - 3) c.center.Y (c.center.X + 3) c.center.Y + drawLine img crossColor c.center.X (c.center.Y - 3) c.center.X (c.center.Y + 3) + +let drawCells (img: Image) (drawCellContent: bool) (cells: Types.Cell list) = + List.iter (fun c -> drawCell img drawCellContent c) cells \ No newline at end of file diff --git a/Parasitemia/Parasitemia/KMedians.fs b/Parasitemia/Parasitemia/KMedians.fs index 09d6c8f..82e09cb 100644 --- a/Parasitemia/Parasitemia/KMedians.fs +++ b/Parasitemia/Parasitemia/KMedians.fs @@ -22,28 +22,28 @@ let kmedians (img: Image) (fgFactor: float) : Result = let max = ref [| 0.0 |] let maxLocation = ref <| [| Point() |] img.MinMax(min, max, minLocation, maxLocation) - + 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) - + 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 index dd0247a..3e714ac 100644 --- a/Parasitemia/Parasitemia/KdTree.fs +++ b/Parasitemia/Parasitemia/KdTree.fs @@ -39,7 +39,7 @@ type Tree<'a when 'a :> I2DCoords> = | Node of float * Tree<'a> * Tree<'a> | Leaf of 'a - static member buildTree (l: 'a list) : Tree<'a> = + static member BuildTree (l: 'a list) : Tree<'a> = let xSorted = List.toArray l let ySorted = List.toArray l Array.sortInPlaceWith cmpX xSorted @@ -64,7 +64,7 @@ type Tree<'a when 'a :> I2DCoords> = buildTreeFromSortedArray xSorted ySorted 1 - static member search (tree: Tree<'a>) (searchRegion: Region) : 'a list = + member this.Search (searchRegion: Region) : 'a list = let rec valuesFrom (tree: Tree<'a>) : 'a list = match tree with | Leaf v -> [v] @@ -94,7 +94,7 @@ type Tree<'a when 'a :> I2DCoords> = 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 + searchWithRegion this { minX = Double.MinValue; maxX = Double.MaxValue; minY = Double.MinValue; maxY = Double.MaxValue } 1 ///// Tests. TODO: to put in a unit test. @@ -118,16 +118,16 @@ let test () = Point(3.5, 1.5) Point(2.5, 0.5) ] - let tree = Tree.buildTree pts + 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. + let s1 = tree.Search { 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 } + let s2 = tree.Search { 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 } + let s3 = tree.Search { minX = 2.0; maxX = 2.0; minY = 2.0; maxY = 2.0 } Utils.dprintfn "s3: %A" s3 let test2 () = @@ -136,19 +136,19 @@ let test2 () = Point(1.0, 2.0) Point(1.0, 3.0) ] - let tree = Tree.buildTree pts + 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 } + let s1 = tree.Search { 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 } + let s2 = tree.Search { 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 } + let s3 = tree.Search { 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 } + let s4 = tree.Search { minX = 0.0; maxX = 2.0; minY = 0.0; maxY = 4.0 } Utils.dprintfn "s4: %A" s4 diff --git a/Parasitemia/Parasitemia/MainAnalysis.fs b/Parasitemia/Parasitemia/MainAnalysis.fs new file mode 100644 index 0000000..3a38c4e --- /dev/null +++ b/Parasitemia/Parasitemia/MainAnalysis.fs @@ -0,0 +1,97 @@ +module ImageAnalysis + +open System +open System.Drawing + +open Emgu.CV +open Emgu.CV.Structure + +open Utils +open ImgTools +open Config +open Types + + +let doAnalysis (img: Image) (name: string) (config: Config) : Cell list = + + let imgFloat = img.Convert() + use scaledImg = if config.scale = 1.0 then imgFloat else imgFloat.Resize(config.scale, CvEnum.Inter.Area) + + use green = scaledImg.Item(1) + + 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 ] + [ 1.0f; 0.0f; -1.0f ]], Point(0, 0)) + + use xEdges = filteredGreen.Convolution(sobelKernel).Convert() + use yEdges = filteredGreen.Convolution(sobelKernel.Transpose()).Convert() + + let xEdgesData = xEdges.Data + let yEdgesData = yEdges.Data + for r in 0..xEdges.Rows-1 do + xEdgesData.[r, 0, 0] <- 0.0 + xEdgesData.[r, xEdges.Cols-1, 0] <- 0.0 + yEdgesData.[r, 0, 0] <- 0.0 + yEdgesData.[r, xEdges.Cols-1, 0] <- 0.0 + + for c in 0..xEdges.Cols-1 do + xEdgesData.[0, c, 0] <- 0.0 + xEdgesData.[xEdges.Rows-1, c, 0] <- 0.0 + yEdgesData.[0, c, 0] <- 0.0 + yEdgesData.[xEdges.Rows-1, c, 0] <- 0.0 + + use magnitudes = new Matrix(xEdges.Size) + CvInvoke.CartToPolar(xEdges, yEdges, magnitudes, new Mat()) // Compute the magnitudes (without angles). + + let min = ref 0.0 + let minLocation = ref <| Point() + let max = ref 0.0 + let maxLocation = ref <| Point() + magnitudes.MinMax(min, max, minLocation, maxLocation) + + use magnitudesByte = ((magnitudes / !max) * 255.0).Convert() // Otsu from OpenCV only support 'byte'. + use edges = new Matrix(xEdges.Size) + CvInvoke.Threshold(magnitudesByte, edges, 0.0, 1.0, CvEnum.ThresholdType.Otsu ||| CvEnum.ThresholdType.Binary) |> ignore + + logTime "Finding edges" (fun() -> thin edges) + logTime "Removing small connected components from thinning" (fun () -> removeArea edges 12) + + let kmediansResults = KMedians.kmedians filteredGreen 1.0 + + let parasites = ParasitesMarker.find green filteredGreen kmediansResults config + + let radiusRange = config.scale * config.minRBCSize, config.scale * config.maxRBCSize + let windowSize = roundInt (1.6 * (snd radiusRange)) + let factorNbPick = 1.5 + let ellipses = logTime "Finding ellipses" (fun () -> + Ellipse.find edges xEdges yEdges radiusRange windowSize factorNbPick) + + let cells = logTime "Classifier" (fun () -> Classifier.findCells ellipses parasites kmediansResults.fg config) + + match config.debug with + | DebugOn output -> + let buildFileName postfix = System.IO.Path.Combine(output, name + postfix) + saveMat (edges * 255.0) (buildFileName " - edges.png") + saveImg parasites.darkStain (buildFileName " - parasites - dark stain.png") + saveImg parasites.stain (buildFileName " - parasites - stain.png") + saveImg parasites.infection (buildFileName " - parasites - infection.png") + + let imgEllipses = img.Copy() + drawEllipses imgEllipses ellipses (Bgr(0.0, 240.0, 240.0)) + saveImg imgEllipses (buildFileName " - ellipses.png") + + saveImg (kmediansResults.fg * 255.0) (buildFileName " - foreground.png") + + let imgCells = img.Copy() + drawCells imgCells false cells + saveImg imgCells (buildFileName " - cells.png") + + let imgCells' = img.Copy() + drawCells imgCells' true cells + saveImg imgCells' (buildFileName " - cells - full.png") + | _ -> () + + cells diff --git a/Parasitemia/Parasitemia/MatchingEllipses.fs b/Parasitemia/Parasitemia/MatchingEllipses.fs index 0964953..8a1b4f8 100644 --- a/Parasitemia/Parasitemia/MatchingEllipses.fs +++ b/Parasitemia/Parasitemia/MatchingEllipses.fs @@ -12,86 +12,90 @@ open Utils let matchingScoreThreshold1 = 0.6 let matchingScoreThreshold2 = 1.0 -type EllipseScore (matchingScore: float, e: Ellipse) = +type private EllipseScoreFlaggedKd (matchingScore: float, e: Ellipse) = let mutable matchingScore = matchingScore - member this.MatchingScore = matchingScore member this.Ellipse = e - member val Processed = false with get, set - member val Removed = false with get, set + member this.MatchingScore = matchingScore member this.AddMatchingScore(score: float) = matchingScore <- matchingScore + score - + + member val Processed = false with get, set + member val Removed = false with get, set + interface KdTree.I2DCoords with - member this.X = e.Cx - member this.Y = e.Cy + member this.X = this.Ellipse.Cx + member this.Y = this.Ellipse.Cy + type MatchingEllipses (radiusMin: float) = - let ellipses = List() + let ellipses = List() member this.Add (e: Ellipse) = - ellipses.Add(EllipseScore(0.0, e)) + ellipses.Add(EllipseScoreFlaggedKd(0.0, e)) 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 - match EEOver.EEOverlapArea e.Ellipse other.Ellipse with - | Some (commonArea, _, _) -> - 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 + if ellipses.Count = 0 then - for j in i .. nbToRemove - 1 do - ellipses.[j].Removed <- true - ellipses.RemoveRange(i, nbToRemove) + [] + else + // 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 tree.Search window do + if not other.Processed + then + let areaOther = other.Ellipse.Area + match EEOver.EEOverlapArea e.Ellipse other.Ellipse with + | Some (commonArea, _, _) -> + let matchingScore = 2.0 * commonArea / (areaE + areaOther) + if matchingScore >= matchingScoreThreshold1 + then + other.AddMatchingScore(matchingScore) + e.AddMatchingScore(matchingScore) + | _ -> () - // 5) Remove ellipses whose center is into an ellipse with a better score - for e in ellipses do - if not e.Removed + // 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(EllipseScoreFlaggedKd(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 - 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 - - List.ofSeq ellipses |> List.map (fun e -> e.Ellipse) + 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 + 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 tree.Search 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 + + List.ofSeq ellipses |> List.map (fun e -> e.Ellipse) - - diff --git a/Parasitemia/Parasitemia/Parasitemia.fsproj b/Parasitemia/Parasitemia/Parasitemia.fsproj index a14f11c..e3eff91 100644 --- a/Parasitemia/Parasitemia/Parasitemia.fsproj +++ b/Parasitemia/Parasitemia/Parasitemia.fsproj @@ -28,6 +28,7 @@ x64 bin\Debug\Parasitemia.XML false + --folder "../../../../Images/08.10.2015/test" --output "../../../Images/output" --debug pdbonly @@ -39,6 +40,7 @@ x64 bin\Release\Parasitemia.XML false + --folder "../../../../Images/08.10.2015/Without scale" --output "../../../Images/output" --debug 11 @@ -62,18 +64,18 @@ - + - + - + diff --git a/Parasitemia/Parasitemia/ParasitesMarker.fs b/Parasitemia/Parasitemia/ParasitesMarker.fs index fde0afe..e422b04 100644 --- a/Parasitemia/Parasitemia/ParasitesMarker.fs +++ b/Parasitemia/Parasitemia/ParasitesMarker.fs @@ -10,30 +10,36 @@ type Result = { stain: Image infection: Image } +// Create three binary markers : +// * 'Dark stain' corresponds to the colored pixel, it's independent of the size of the areas. +// * 'Stain' corresponds to the stain around the parasites. +// * 'Infection' corresponds to the parasite. It shouldn't contain thrombocytes. let find (green: Image) (filteredGreen: Image) (kmediansResult: KMedians.Result) (config: Config.Config) : Result = + let green = ImgTools.gaussianFilter green 1.0 + // 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 } = kmediansResult 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 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 + 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 b44c65f..ac17c82 100644 --- a/Parasitemia/Parasitemia/Program.fs +++ b/Parasitemia/Parasitemia/Program.fs @@ -1,5 +1,6 @@ module Parasitemia.Main +open System open System.IO open System.Windows open System.Windows.Media @@ -15,58 +16,108 @@ open Emgu.CV.WPF open Config -let display (window : Views.MainWindow) (img : IImage) = +let display (window : Views.MainWindow) (img : IImage) = let imgControl = window.Root.FindName("img") :?> Controls.Image imgControl.Source <- BitmapSourceConvert.ToBitmapSource(img) let log (window : Views.MainWindow) (mess : string) = let txtLog = window.Root.FindName("txtLog") :?> Controls.TextBlock txtLog.Text <- txtLog.Text + mess + "\n" - -[] -do - Utils.dprintfn "OpenCV test" - - let app = new Application() - let mainWindow = Views.MainWindow() - - Utils.log <- (fun m -> log mainWindow m) - - let config = { - scale = 1.0 - - doGSigma1 = 1.5 - doGSigma2 = 20.0 - doGLowFreqPercentageReduction = 0.9 - - darkStainLevel = 0.839 - - stainSigma = 10.0 - stainLevel = 0.9 - stainSpreadRequired = 3.0 - - infectionSigma = 2.0 - infectionLevel = 0.762 - infectionPixelsRequired = 1 } - - ///// 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_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 = Utils.logTime "Total" (fun () -> - ImageAnalysis.doAnalysis img config ) - - display mainWindow img - mainWindow.Root.Show() - app.Run() |> ignore + + +type Input = + | File of string + | Dir of string + +type RunningMode = + | CmdLine of Input * string + | Window + +type Arguments = RunningMode * bool + +let parseArgs (args: string[]) : Arguments = + + let output = Array.tryFindIndex ((=) "--output") args + + let runningMode = + match Array.tryFindIndex ((=) "--folder") args, output with + | Some i, Some i_output when i < args.Length - 2 && i_output < args.Length - 2 -> + CmdLine ((Dir args.[i+1]), args.[i_output + 1]) + | _ -> + match Array.tryFindIndex ((=) "--file") args, output with + | Some i, Some i_output when i < args.Length - 2 && i_output < args.Length - 2 -> + CmdLine ((File args.[i+1]), args.[i_output + 1]) + |_ -> + Window + + runningMode, Array.exists ((=) "--debug") args + + +[] +let main args = + match parseArgs args with + | mode, debug -> + let scale = 1. + let config = { + debug = if debug then DebugOn "." else DebugOff + + scale = scale + + // RBC size range in px at scale = 1.0. + minRBCSize = 20. + maxRBCSize = 40. + + doGSigma1 = 1.5 + doGSigma2 = 20. + doGLowFreqPercentageReduction = 0.75 + + darkStainLevel = 0.5 + + stainSigma = 10. + stainLevel = 0.9 + stainSpreadRequired = 3.0 + + infectionSigma = 2.2 + infectionLevel = 0.85 + infectionPixelsRequired = 1 + + percentageOfFgValidCell = 0.4 + + MaxDarkStainRatio = 0.1 + + minimumCellArea = 600. * scale ** 2. |> int + maxOffcenter = 0.5 } + + match mode with + | CmdLine (input, output) -> + let config = { config with debug = DebugOn output } + Utils.log <- (fun m -> Console.WriteLine m) + + Directory.CreateDirectory output |> ignore + + let files = match input with + | File file -> [ file ] + | Dir dir -> Directory.EnumerateFiles dir |> List.ofSeq + + use resultFile = new StreamWriter(new FileStream(Path.Combine(output, "results.txt"), FileMode.Create, FileAccess.Write)) + + for file in files do + try + let fileInfo = FileInfo(file) + use img = new Image(file) + let cells = Utils.logTime "Whole analyze" (fun () -> ImageAnalysis.doAnalysis img fileInfo.Name config) + let total, infected = Utils.countCells cells + fprintf resultFile "File: %s %d %d %.2f\n" file total infected (100. * (float infected) / (float total)) + with + | _ as ex -> Utils.log (sprintf "Unable to open the image '%A': %A" file ex) + 0 + + | Window -> + let app = new Application() + let mainWindow = Views.MainWindow() + + Utils.log <- (fun m -> log mainWindow m) + + //display mainWindow img + mainWindow.Root.Show() + app.Run() diff --git a/Parasitemia/Parasitemia/Types.fs b/Parasitemia/Parasitemia/Types.fs index f96573d..e93bbb0 100644 --- a/Parasitemia/Parasitemia/Types.fs +++ b/Parasitemia/Parasitemia/Types.fs @@ -1,6 +1,10 @@ module Types open System +open System.Drawing + +open Emgu.CV +open Emgu.CV.Structure type Ellipse (cx: float, cy: float, a: float, b: float, alpha: float) = member this.Cx = cx @@ -11,17 +15,41 @@ type Ellipse (cx: float, cy: float, a: float, b: float, alpha: float) = member this.Area = a * b * Math.PI // Does the ellipse contain the point (x, y)?. - member this.Contains 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 // A line is defined as : y = mx + l member this.CutALine (m: float) (l: float) : bool = - -2.0 * l ** 2.0 + a ** 2.0 + m ** 2.0 * a ** 2.0 + b ** 2.0 + m ** 2.0 * b ** 2.0 - - 4.0 * m * l * cx - 2.0 * m ** 2.0 * cx ** 2.0 + 4.0 * l * cy + 4.0 * m * cx * cy - + -2.0 * l ** 2.0 + a ** 2.0 + m ** 2.0 * a ** 2.0 + b ** 2.0 + m ** 2.0 * b ** 2.0 - + 4.0 * m * l * cx - 2.0 * m ** 2.0 * cx ** 2.0 + 4.0 * l * cy + 4.0 * m * cx * cy - 2.0 * cy ** 2.0 + (-1.0 + m ** 2.0) * (a ** 2.0 - b ** 2.0) * cos (2.0 * alpha) - 2.0 * m * (a ** 2.0 - b ** 2.0) * sin (2.0 * alpha) > 0.0 - member this.CutAVericalLine (y: float) : bool = + member this.CutAVericalLine (y: float) : bool = a ** 2.0 + b ** 2.0 - 2.0 * y ** 2.0 + 4.0 * y * cx - 2.0 * cx ** 2.0 + a ** 2.0 * cos (2.0 * alpha) - b ** 2.0 * cos (2.0 * alpha) > 0.0 - member this.CutAnHorizontalLine (x: float) : bool = - a ** 2.0 + b ** 2.0 - 2.0 * x ** 2.0 + 4.0 * x * cy - 2.0 * cy ** 2.0 - a ** 2.0 * cos (2.0 * alpha) + b ** 2.0 * cos (2.0 * alpha) > 0.0 \ No newline at end of file + member this.CutAnHorizontalLine (x: float) : bool = + a ** 2.0 + b ** 2.0 - 2.0 * x ** 2.0 + 4.0 * x * cy - 2.0 * cy ** 2.0 - a ** 2.0 * cos (2.0 * alpha) + b ** 2.0 * cos (2.0 * alpha) > 0.0 + + member this.isOutside (width: float) (height: float) = + this.Cx <= 0.0 || this.Cx >= width || + this.Cy <= 0.0 || this.Cy >= height || + this.CutAVericalLine 0.0 || this.CutAVericalLine width || + this.CutAnHorizontalLine 0.0 || this.CutAnHorizontalLine height + + +type CellClass = HealthyRBC | InfectedRBC | Peculiar + +type Cell = { + cellClass: CellClass + center: Point + elements: Matrix } + + +type Line (a: float, b: float) = + member this.A = a + member this.B = b + +type PointD (x: float, y: float) = + member this.X = x + member this.Y = y + diff --git a/Parasitemia/Parasitemia/Utils.fs b/Parasitemia/Parasitemia/Utils.fs index de16664..d6fa2f0 100644 --- a/Parasitemia/Parasitemia/Utils.fs +++ b/Parasitemia/Parasitemia/Utils.fs @@ -2,6 +2,8 @@ open System.Diagnostics +open Types + let roundInt = int << round let inline dprintfn fmt = @@ -16,4 +18,27 @@ let logTime (m: string) (f: unit -> 'a) : 'a = let res = f () sw.Stop() log <| sprintf "%A (time: %A ms)" m sw.ElapsedMilliseconds - res \ No newline at end of file + res + +let inline lineFromTwoPoints (p1: PointD) (p2: PointD) : Line = + let a = (p1.Y - p2.Y) / (p1.X - p2.X) + let b = -(p2.X * p1.Y - p1.X * p2.Y) / (p1.X - p2.X) + Line(a, b) + +let inline pointFromTwoLines (l1: Line) (l2: Line) : PointD = + let x = -(l1.B - l2.B) / (l1.A - l2.A) + let y = -(l2.A * l1.B - l1.A * l2.B) / (l1.A - l2.A) + PointD(x, y) + +let inline squaredDistanceTwoPoints (p1: PointD) (p2: PointD) = + (p1.X - p2.X) ** 2.0 + (p1.Y - p2.Y) ** 2.0 + +let distanceTwoPoints (p1: PointD) (p2: PointD) = + squaredDistanceTwoPoints p1 p2 |> sqrt + +let countCells (cells: Cell list) : int * int = + cells |> List.fold (fun (total, infected) { cellClass = cellClass } -> + match cellClass with + | HealthyRBC -> (total + 1, infected) + | InfectedRBC -> (total + 1, infected + 1) + | Peculiar -> (total, infected)) (0, 0) \ No newline at end of file -- 2.43.0