The main process is now complete.
authorGreg Burri <greg.burri@gmail.com>
Mon, 14 Dec 2015 01:29:31 +0000 (02:29 +0100)
committerGreg Burri <greg.burri@gmail.com>
Mon, 14 Dec 2015 01:29:31 +0000 (02:29 +0100)
15 files changed:
Parasitemia/Parasitemia/Classifier.fs
Parasitemia/Parasitemia/Config.fs
Parasitemia/Parasitemia/EEOver.fs
Parasitemia/Parasitemia/Ellipse.fs
Parasitemia/Parasitemia/ImageAnalysis.fs [deleted file]
Parasitemia/Parasitemia/ImgTools.fs
Parasitemia/Parasitemia/KMedians.fs
Parasitemia/Parasitemia/KdTree.fs
Parasitemia/Parasitemia/MainAnalysis.fs [new file with mode: 0644]
Parasitemia/Parasitemia/MatchingEllipses.fs
Parasitemia/Parasitemia/Parasitemia.fsproj
Parasitemia/Parasitemia/ParasitesMarker.fs
Parasitemia/Parasitemia/Program.fs
Parasitemia/Parasitemia/Types.fs
Parasitemia/Parasitemia/Utils.fs

index 24ea4a7..37507c5 100644 (file)
 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<byte> }
+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<Gray, byte>) : Cell list =
+
+
+let findCells (ellipses: Ellipse list) (parasites: ParasitesMarker.Result) (fg: Image<Gray, byte>) (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<Point>()
+                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<byte>(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 })
index 75aa95f..8a172b8 100644 (file)
@@ -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
index 640d752..52b5963 100644 (file)
@@ -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
index 5bfcc33..62c0936 100644 (file)
@@ -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<byte>)
          (xDir: Image<Gray, float>) 
          (yDir: Image<Gray, float>) 
diff --git a/Parasitemia/Parasitemia/ImageAnalysis.fs b/Parasitemia/Parasitemia/ImageAnalysis.fs
deleted file mode 100644 (file)
index feb427d..0000000
+++ /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<Bgr, byte>
-}*)
-
-let doAnalysis (img: Image<Bgr, byte>) (config: Config) : Classifier.Cell list =
-
-    let imgFloat = img.Convert<Bgr, float32>()
-    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<byte>(scaledImg.Size)
-    //CvInvoke.MixChannels(scaledImg, green, [| 1; 0 |])
-
-    //let greenMatrix = new Matrix<byte>(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<Gray, float>()
-    use yEdges = filteredGreen.Convolution(sobelKernel.Transpose()).Convert<Gray, float>()
-            
-    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<float>(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<byte>() // Otsu from OpenCV only support 'byte'.
-    use edges = new Matrix<byte>(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<Hsv, uint8>()
-    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<Gray, byte> = H.Convert(fun b _ _ ->
-        (255 - int(b) * 255 / 179 + hueShiftValue) % 256 |> byte
-    )
-   
-    let correctedS : Image<Gray, byte> = 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
index 5278313..c34b9c7 100644 (file)
@@ -10,7 +10,7 @@ open Emgu.CV.Structure
 open Utils
 
 // Normalize image values between 0uy and 255uy.
-let normalizeAndConvert (img: Image<Gray, float32>) : Image<Gray, byte> = 
+let normalizeAndConvert (img: Image<Gray, float32>) : Image<Gray, byte> =
     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<byte>) =    
-    let neighbors = [| 
+let thin (mat: Matrix<byte>) =
+    let neighbors = [|
         (-1,  0) // p2
         (-1,  1) // p3
         ( 0,  1) // p4
@@ -41,15 +41,15 @@ let thin (mat: Matrix<byte>) =
     let mutable data2 = Array2D.zeroCreate<byte> 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<byte>) =
         then
             nb <- nb + 1
         nb
-         
+
     let mutable pixelChanged = true
     let mutable oddIteration = true
     while pixelChanged do
@@ -84,16 +84,22 @@ let thin (mat: Matrix<byte>) =
                 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<byte>) (areaSize: int) =
-    let neighbors = [| 
+    let neighbors = [|
         (-1,  0) // p2
         (-1,  1) // p3
         ( 0,  1) // p4
@@ -102,29 +108,24 @@ let removeArea (mat: Matrix<byte>) (areaSize: int) =
         ( 1, -1) // p7
         ( 0, -1) // p8
         (-1, -1) |] // p9
-    
+
     let mat' = new Matrix<byte>(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<byte>) (areaSize: int) =
                     for (ni, nj) in neighborhood do
                         data.[ni, nj] <- 0uy
 
+let connectedComponents (img: Image<Gray, byte>) (startPoints: List<Point>) : List<Point> =
+    let w = img.Width
+    let h = img.Height
+
+    let pointChecked = HashSet<Point>()
+    let pointToCheck = List<Point>(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<Point>(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<Gray, 'TDeph>(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<Bgr, byte>) (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<Bgr, byte>) (drawCellContent: bool) (cells: Types.Cell list) =
+    List.iter (fun c -> drawCell img drawCellContent c) cells
\ No newline at end of file
index 09d6c8f..82e09cb 100644 (file)
@@ -22,28 +22,28 @@ let kmedians (img: Image<Gray, float32>) (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<Gray, float32>(img.Size)
     let mutable d_fg = new Image<Gray, float32>(img.Size)
     let mutable fg = new Image<Gray, byte>(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 }
index dd0247a..3e714ac 100644 (file)
@@ -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 (file)
index 0000000..3a38c4e
--- /dev/null
@@ -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<Bgr, byte>) (name: string) (config: Config) : Cell list =
+
+    let imgFloat = img.Convert<Bgr, float32>()
+    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<Gray, float>()
+    use yEdges = filteredGreen.Convolution(sobelKernel.Transpose()).Convert<Gray, float>()
+
+    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<float>(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<byte>() // Otsu from OpenCV only support 'byte'.
+    use edges = new Matrix<byte>(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
index 0964953..8a1b4f8 100644 (file)
@@ -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<EllipseScore>()
+    let ellipses = List<EllipseScoreFlaggedKd>()
         
     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<EllipseScore> 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<EllipseScoreFlaggedKd> 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)
 
-        
-        
         
 
index a14f11c..e3eff91 100644 (file)
@@ -28,6 +28,7 @@
     <PlatformTarget>x64</PlatformTarget>
     <DocumentationFile>bin\Debug\Parasitemia.XML</DocumentationFile>
     <Prefer32Bit>false</Prefer32Bit>
+    <StartArguments>--folder "../../../../Images/08.10.2015/test" --output "../../../Images/output" --debug</StartArguments>
   </PropertyGroup>
   <PropertyGroup Condition=" '$(Configuration)|$(Platform)' == 'Release|AnyCPU' ">
     <DebugType>pdbonly</DebugType>
@@ -39,6 +40,7 @@
     <PlatformTarget>x64</PlatformTarget>
     <DocumentationFile>bin\Release\Parasitemia.XML</DocumentationFile>
     <Prefer32Bit>false</Prefer32Bit>
+    <StartArguments>--folder "../../../../Images/08.10.2015/Without scale" --output "../../../Images/output" --debug</StartArguments>
   </PropertyGroup>
   <PropertyGroup>
     <MinimumVisualStudioVersion Condition="'$(MinimumVisualStudioVersion)' == ''">11</MinimumVisualStudioVersion>
     <Compile Include="NumericUpDown.xaml.fs" />
     <Resource Include="MainWindow.xaml" />
     <Compile Include="MainWindow.xaml.fs" />
-    <Compile Include="Utils.fs" />
     <Compile Include="Types.fs" />
+    <Compile Include="Utils.fs" />
     <Compile Include="ImgTools.fs" />
     <Compile Include="Config.fs" />
     <Compile Include="KMedians.fs" />
     <Compile Include="EEOver.fs" />
     <Compile Include="ParasitesMarker.fs" />
     <Compile Include="KdTree.fs" />
-    <Compile Include="Classifier.fs" />
     <Compile Include="MatchingEllipses.fs" />
+    <Compile Include="Classifier.fs" />
     <Compile Include="Ellipse.fs" />
-    <Compile Include="ImageAnalysis.fs" />
+    <Compile Include="MainAnalysis.fs" />
     <Compile Include="Program.fs" />
     <None Include="App.config" />
     <Content Include="packages.config" />
index fde0afe..e422b04 100644 (file)
@@ -10,30 +10,36 @@ type Result = {
     stain: Image<Gray, byte>
     infection: Image<Gray, byte> }
 
+// 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<Gray, float32>) (filteredGreen: Image<Gray, float32>) (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<Gray, float32>()
     let greenWithoutBg = green.Copy()
     greenWithoutBg.SetValue(Gray(0.0), fg.Not())
 
-    let findSmears (sigma: float) (level: float) : Image<Gray, byte> =       
+    let findSmears (sigma: float) (level: float) : Image<Gray, byte> =
         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 }
-    
+
 
 
 
index b44c65f..ac17c82 100644 (file)
@@ -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"
-   
-[<System.STAThreadAttribute>]
-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<Bgr, byte>("../../../../src/Tests_hough/images/1508133543-7-4-120-0001.png")
-    //use img = new Image<Bgr, byte>("../../../../src/Tests_hough/images/rbc_single.png")
-    //use img = new Image<Bgr, byte>("../../../../src/Tests_hough/images/rbc_single_oblong_4.png")
-    //use img = new Image<Bgr, byte>("../../../../src/Tests_hough/images/strange_rbc_1.png")
-    //use img = new Image<Bgr, byte>("../../../../src/Tests_hough/images/rbc_single_blurred.png")   
-    use img = new Image<Bgr, byte>("../../../../src/Tests_hough/images/lot.png")
-
-
-    ///// PARASITES /////
-    //use img = new Image<Bgr, byte>("../../../../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
+
+
+[<EntryPoint>]
+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<Bgr, byte>(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()
index f96573d..e93bbb0 100644 (file)
@@ -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<byte> }
+
+
+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
+
index de16664..d6fa2f0 100644 (file)
@@ -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