* Treat some special cases when ellipses intersecting.

* Some adjustments.

diff --git a/Parasitemia/Parasitemia/Classifier.fs b/Parasitemia/Parasitemia/Classifier.fs
index 99b6edc..7f11e13 100644 (file)
--- a/Parasitemia/Parasitemia/Classifier.fs
+++ b/Parasitemia/Parasitemia/Classifier.fs
@@ -29,10 +29,10 @@ let findCells (ellipses: Ellipse list) (parasites: ParasitesMarker.Result) (img:
         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 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.RBCMaxRadius) * config.Parameters.scale
-                                          KdTree.maxX = e.Cx + (e.A + config.RBCMaxRadius) * config.Parameters.scale
-                                          KdTree.minY = e.Cy - (e.A + config.RBCMaxRadius) * config.Parameters.scale
-                                          KdTree.maxY = e.Cy + (e.A + config.RBCMaxRadius) * config.Parameters.scale }
+        let searchRegion (e: Ellipse) = { KdTree.minX = e.Cx - (e.A + config.RBCMaxRadius)
+                                          KdTree.maxX = e.Cx + (e.A + config.RBCMaxRadius)
+                                          KdTree.minY = e.Cy - (e.A + config.RBCMaxRadius)
+                                          KdTree.maxY = e.Cy + (e.A + config.RBCMaxRadius) }

 
         // The minimum window to contain a given ellipse.
         let ellipseWindow (e: Ellipse) =
 
         // The minimum window to contain a given ellipse.
         let ellipseWindow (e: Ellipse) =


@@ -47,39 +47,57 @@ let findCells (ellipses: Ellipse list) (parasites: ParasitesMarker.Result) (img:
 
         // Return 'true' if the point 'p' is owned by e.
         // The lines represents all intersections with other ellipses.
 
         // 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) =
+        let pixelOwnedByE (p: PointD) (e: Ellipse) (others: (Ellipse * Line) list) =

             e.Contains p.X p.Y &&
             seq {
                 let c = PointD(e.Cx, e.Cy)
             e.Contains p.X p.Y &&
             seq {
                 let c = PointD(e.Cx, e.Cy)

-                for d1 in lines do
+                for e', d1 in others do

                     let d2 = Utils.lineFromTwoPoints c p
                     let d2 = Utils.lineFromTwoPoints c p

+                    let c' = PointD(e'.Cx, e'.Cy)
+                    let v = pointFromTwoLines d1 (lineFromTwoPoints c c')
+                    let case1 = sign (v.X - c.X) <> sign (v.X - c'.X) || Utils.squaredDistanceTwoPoints v c > Utils.squaredDistanceTwoPoints v c'

                     if d2.Valid
                     then
                         let p' = Utils.pointFromTwoLines d1 d2
                     if d2.Valid
                     then
                         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.
+                        // Yield 'false' when the point is owned by another ellipse.
+                        if case1
+                        then
+                            yield sign (c.X - p.X) <> sign (c.X - p'.X) || Utils.squaredDistanceTwoPoints c p' > Utils.squaredDistanceTwoPoints c p
+                        else
+                            yield sign (c.X - p.X) = sign (c.X - p'.X) && Utils.squaredDistanceTwoPoints c p' < Utils.squaredDistanceTwoPoints c p

                     else
                     else

-                        yield true
+                        yield case1

             } |> Seq.forall id
 
         let ellipses = ellipses |> List.map EllipseFlaggedKd
 
             } |> Seq.forall id
 
         let ellipses = ellipses |> List.map EllipseFlaggedKd
 

-        // 1) Associate touching ellipses with each ellipses.
+        // 1) Associate touching ellipses with each ellipses and remove ellipse with more than two intersections.

         let tree = KdTree.Tree.BuildTree ellipses
         let tree = KdTree.Tree.BuildTree ellipses

-        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 neighbors (e: EllipseFlaggedKd) : (EllipseFlaggedKd * PointD * PointD) list =
+            if not e.Removed
+            then
+                tree.Search (searchRegion e)
+                    // We only keep the ellipses touching 'e'.
+                    |> List.choose (fun otherE ->
+                        match EEOver.EEOverlapArea e otherE with
+                        | Some (_, px, _) when px.Length > 2 ->
+                            otherE.Removed <- true
+                            None
+                        | Some (area, px, py) when area > 0.0 && px.Length = 2 ->
+                            Some (otherE, PointD(px.[0], py.[0]), PointD(px.[1], py.[1]))
+                        | _ ->
+                            None )
+            else
+                []

 
         // We reverse the list to get the lower score ellipses first.
         let ellipsesWithNeigbors = ellipses |> List.map (fun e -> e, neighbors e) |> List.rev
 
         // 2) Remove ellipses with a high standard deviation (high contrast).
 
         // We reverse the list to get the lower score ellipses first.
         let ellipsesWithNeigbors = ellipses |> List.map (fun e -> e, neighbors e) |> List.rev
 
         // 2) Remove ellipses with a high standard deviation (high contrast).

-        let globalStdDiviation = MathNet.Numerics.Statistics.StreamingStatistics.StandardDeviation(seq {
+
+        // CvInvoke.Normalize(img, img, 0.0, 255.0, CvEnum.NormType.MinMax) // Not needed.
+
+        let globalStdDeviation = MathNet.Numerics.Statistics.Statistics.StandardDeviation(seq {

             for y in 0 .. h - 1 do
                 for x in 0 .. w - 1 do
                     yield img.Data.[y, x, 0] |> float })
             for y in 0 .. h - 1 do
                 for x in 0 .. w - 1 do
                     yield img.Data.[y, x, 0] |> float })


@@ -87,15 +105,14 @@ let findCells (ellipses: Ellipse list) (parasites: ParasitesMarker.Result) (img:
         for e in ellipses do
             let minX, minY, maxX, maxY = ellipseWindow e
 
         for e in ellipses do
             let minX, minY, maxX, maxY = ellipseWindow e
 

-            let stdDiviation = MathNet.Numerics.Statistics.StreamingStatistics.StandardDeviation(seq {
+            let stdDeviation = MathNet.Numerics.Statistics.Statistics.StandardDeviation (seq {

                 for y in (if minY < 0 then 0 else minY) .. (if maxY >= h then h - 1 else maxY) do
                     for x in (if minX < 0 then 0 else minX) .. (if maxX >= w then w - 1 else maxX) do
                         if e.Contains (float x) (float y)
                         then
                 for y in (if minY < 0 then 0 else minY) .. (if maxY >= h then h - 1 else maxY) do
                     for x in (if minX < 0 then 0 else minX) .. (if maxX >= w then w - 1 else maxX) do
                         if e.Contains (float x) (float y)
                         then

-                            yield img.Data.[y, x, 0] |> float })
+                            yield float img.Data.[y, x, 0] })

 
 

-            if stdDiviation > globalStdDiviation * config.Parameters.standardDeviationMaxRatio
-            then
+            if stdDeviation > globalStdDeviation * config.Parameters.standardDeviationMaxRatio then

                 e.Removed <- true
 
         // 3) Remove ellipses touching the edges.
                 e.Removed <- true
 
         // 3) Remove ellipses touching the edges.


@@ -113,7 +130,7 @@ let findCells (ellipses: Ellipse list) (parasites: ParasitesMarker.Result) (img:
                 for y in (if minY < 0 then 0 else minY) .. (if maxY >= h then h - 1 else maxY) do
                     for x in (if minX < 0 then 0 else minX) .. (if maxX >= w then w - 1 else maxX) do
                         let p = PointD(float x, float y)
                 for y in (if minY < 0 then 0 else minY) .. (if maxY >= h then h - 1 else maxY) do
                     for x in (if minX < 0 then 0 else minX) .. (if maxX >= w then w - 1 else maxX) do
                         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)))
+                        if pixelOwnedByE p e (neighbors |> List.choose (fun (otherE, p1, p2) -> if otherE.Removed then None else Some (otherE :> Ellipse, Utils.lineFromTwoPoints p1 p2)))

                         then
                             area <- area + 1
 
                         then
                             area <- area + 1
 


@@ -139,7 +156,7 @@ let findCells (ellipses: Ellipse list) (parasites: ParasitesMarker.Result) (img:
                 for y in minY .. maxY do
                     for x in minX .. maxX do
                         let p = PointD(float x, float y)
                 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)))
+                        if pixelOwnedByE p e (neighbors |> List.choose (fun (otherE, p1, p2) -> if otherE.Removed then None else Some (otherE :> Ellipse, Utils.lineFromTwoPoints p1 p2)))

                         then
                             elements.[y-minY, x-minX] <- 1uy
                             nbElement <- nbElement + 1
                         then
                             elements.[y-minY, x-minX] <- 1uy
                             nbElement <- nbElement + 1

diff --git a/Parasitemia/Parasitemia/Config.fs b/Parasitemia/Parasitemia/Config.fs
index 45b1ca3..5740e8b 100644 (file)
--- a/Parasitemia/Parasitemia/Config.fs
+++ b/Parasitemia/Parasitemia/Config.fs
@@ -30,7 +30,7 @@ type Parameters = {
     infectionLevel: float // > 1
 
     standardDeviationMaxRatio: float // The standard deviation of the pixel values of a cell can't be greater than standardDeviationMaxRatio * global standard deviation
     infectionLevel: float // > 1
 
     standardDeviationMaxRatio: float // The standard deviation of the pixel values of a cell can't be greater than standardDeviationMaxRatio * global standard deviation

-    minimumCellArea: float // Factor of RBC area.
+    minimumCellArea: float // Factor of the nominal RBC area.

 }
 
 type Config (param: Parameters) =
 }
 
 type Config (param: Parameters) =

diff --git a/Parasitemia/Parasitemia/Ellipse.fs b/Parasitemia/Parasitemia/Ellipse.fs
index 48c7f54..73771be 100644 (file)
--- a/Parasitemia/Parasitemia/Ellipse.fs
+++ b/Parasitemia/Parasitemia/Ellipse.fs
@@ -204,11 +204,11 @@ let private areVectorsValid (p1x: float) (p1y: float) (p2x: float) (p2y: float)
 
 
 let find (edges: Matrix<byte>)
 
 
 let find (edges: Matrix<byte>)

-         (xDir: Image<Gray, float>)
-         (yDir: Image<Gray, float>)
+         (xGradient: Image<Gray, float>)
+         (yGradient: Image<Gray, float>)

          (config: Config) : MatchingEllipses =
 
          (config: Config) : MatchingEllipses =
 

-    let r1, r2 = config.Parameters.scale * config.RBCMinRadius, config.Parameters.scale * config.RBCMaxRadius // FIXME: scale factor should be applied in Config!?
+    let r1, r2 = config.RBCMinRadius, config.RBCMaxRadius

     let windowSize = roundInt (config.Parameters.factorWindowSize * r2)
     let factorNbPick = config.Parameters.factorNbPick
 
     let windowSize = roundInt (config.Parameters.factorWindowSize * r2)
     let factorNbPick = config.Parameters.factorNbPick
 


@@ -227,8 +227,8 @@ let find (edges: Matrix<byte>)
     let currentElements = List<(int * int)>()
 
     let edgesData = edges.Data
     let currentElements = List<(int * int)>()
 
     let edgesData = edges.Data

-    let xDirData = xDir.Data
-    let yDirData = yDir.Data
+    let xDirData = xGradient.Data
+    let yDirData = yGradient.Data

 
     let rng = Random(42)
 
 
     let rng = Random(42)
 

diff --git a/Parasitemia/Parasitemia/ImgTools.fs b/Parasitemia/Parasitemia/ImgTools.fs
index 2f021c8..f567e68 100644 (file)
--- a/Parasitemia/Parasitemia/ImgTools.fs
+++ b/Parasitemia/Parasitemia/ImgTools.fs
@@ -21,6 +21,124 @@ let normalizeAndConvert (img: Image<Gray, float32>) : Image<Gray, byte> =
     ((img - (!min).[0]) / ((!max).[0] - (!min).[0]) * 255.0).Convert<Gray, byte>()
 
 
     ((img - (!min).[0]) / ((!max).[0] - (!min).[0]) * 255.0).Convert<Gray, byte>()
 
 

+let saveImg (img: Image<'TColor, 'TDepth>) (filepath: string) =
+    img.Save(filepath)
+
+
+let saveMat (mat: Matrix<'TDepth>) (filepath: string) =
+    use img = new Image<Gray, 'TDeph>(mat.Size)
+    mat.CopyTo(img)
+    saveImg img filepath
+
+
+let suppressMConnections (img: Matrix<byte>) =
+    let w = img.Width
+    let h = img.Height
+    for i in 1 .. h - 2 do
+        for j in 1 .. w - 2 do
+            if img.[i, j] > 0uy && img.Data.[i + 1, j] > 0uy && (img.Data.[i, j - 1] > 0uy && img.Data.[i - 1, j + 1] = 0uy || img.Data.[i, j + 1] > 0uy && img.Data.[i - 1, j - 1] = 0uy) then
+                img.[i, j] <- 0uy
+    for i in 1 .. h - 2 do
+        for j in 1 .. w - 2 do
+            if img.[i, j] > 0uy && img.Data.[i - 1, j] > 0uy && (img.Data.[i, j - 1] > 0uy && img.Data.[i + 1, j + 1] = 0uy || img.Data.[i, j + 1] > 0uy && img.Data.[i + 1, j - 1] = 0uy) then
+                img.[i, j] <- 0uy
+
+let findEdges (img: Image<Gray, float32>) : Matrix<byte> * Image<Gray, float> * Image<Gray, float> =
+    let w = img.Width
+    let h = img.Height
+
+    use sobelKernel =
+        new ConvolutionKernelF(array2D [[ 1.0f; 0.0f; -1.0f ]
+                                        [ 2.0f; 0.0f; -2.0f ]
+                                        [ 1.0f; 0.0f; -1.0f ]], Point(1, 1))
+
+    let xGradient = img.Convolution(sobelKernel).Convert<Gray, float>()
+    let yGradient = img.Convolution(sobelKernel.Transpose()).Convert<Gray, float>()
+
+    let xGradientData = xGradient.Data
+    let yGradientData = yGradient.Data
+    for r in 0 .. h - 1 do
+        xGradientData.[r, 0, 0] <- 0.0
+        xGradientData.[r, w - 1, 0] <- 0.0
+        yGradientData.[r, 0, 0] <- 0.0
+        yGradientData.[r, w - 1, 0] <- 0.0
+
+    for c in 0 .. w - 1 do
+        xGradientData.[0, c, 0] <- 0.0
+        xGradientData.[h - 1, c, 0] <- 0.0
+        yGradientData.[0, c, 0] <- 0.0
+        yGradientData.[h - 1, c, 0] <- 0.0
+
+    use magnitudes = new Matrix<float>(xGradient.Size)
+    CvInvoke.CartToPolar(xGradient, yGradient, magnitudes, new Mat()) // Compute the magnitudes (without angles).
+
+    let thresholdHigh, thresholdLow =
+        let sensibility = 0.1
+        use magnitudesByte = magnitudes.Convert<byte>()
+        let threshold = CvInvoke.Threshold(magnitudesByte, magnitudesByte, 0.0, 1.0, CvEnum.ThresholdType.Otsu ||| CvEnum.ThresholdType.Binary)
+        threshold + (sensibility * threshold), threshold - (sensibility * threshold)
+
+    // Non-maximum suppression.
+    use nms = new Matrix<byte>(xGradient.Size)
+    nms.SetValue(1.0)
+
+    for i in 0 .. h - 1 do
+        nms.Data.[i, 0] <- 0uy
+        nms.Data.[i, w - 1] <- 0uy
+
+    for j in 0 .. w - 1 do
+        nms.Data.[0, j] <- 0uy
+        nms.Data.[h - 1, j] <- 0uy
+
+    for i in 1 .. h - 2 do
+        for j in 1 .. w - 2 do
+            let vx = xGradient.Data.[i, j, 0]
+            let vy = yGradient.Data.[i, j, 0]
+            let angle =
+                let a = atan2 vy vx
+                if a < 0.0 then 2. * Math.PI + a else a
+
+            let mNeigbors (sign: int) : float =
+                if angle < Math.PI / 8. || angle >= 15.0 * Math.PI / 8. then magnitudes.Data.[i, j + sign]
+                elif angle < 3.0 * Math.PI / 8. then magnitudes.Data.[i + sign, j + sign]
+                elif angle < 5.0 * Math.PI / 8. then magnitudes.Data.[i + sign, j]
+                elif angle < 7.0 * Math.PI / 8. then magnitudes.Data.[i + sign, j - sign]
+                elif angle < 9.0 * Math.PI / 8. then magnitudes.Data.[i, j - sign]
+                elif angle < 11.0 * Math.PI / 8. then magnitudes.Data.[i - sign, j - sign]
+                elif angle < 13.0 * Math.PI / 8. then magnitudes.Data.[i - sign, j]
+                else magnitudes.Data.[i - sign, j + sign]
+
+            let m = magnitudes.Data.[i, j]
+            if m < mNeigbors 1 || m < mNeigbors -1 || m < thresholdLow then
+                nms.Data.[i, j] <- 0uy
+
+    // suppressMConnections nms // It's not usefull for the rest of the process (ellipse detection).
+
+    let edges = new Matrix<byte>(xGradient.Size)
+
+    // Histeresis thresholding.
+    let toVisit = Stack<Point>()
+    for i in 0 .. h - 1 do
+        for j in 0 .. w - 1 do
+            if nms.Data.[i, j] = 1uy && magnitudes.Data.[i, j] >= thresholdHigh then
+                nms.Data.[i, j] <- 0uy
+                toVisit.Push(Point(j, i))
+                while toVisit.Count > 0 do
+                    let p = toVisit.Pop()
+                    edges.Data.[p.Y, p.X] <- 1uy
+                    for i' in -1 .. 1  do
+                        for j' in -1 .. 1 do
+                            if i' <> 0 || j' <> 0 then
+                                let ni = p.Y + i'
+                                let nj = p.X + j'
+                                if ni >= 0 && ni < h && nj >= 0 && nj < w && nms.Data.[ni, nj] = 1uy then
+                                    nms.Data.[ni, nj] <- 0uy
+                                    toVisit.Push(Point(nj, ni))
+
+
+    edges, xGradient, yGradient
+
+

 let gaussianFilter (img : Image<'TColor, 'TDepth>) (standardDeviation : float) : Image<'TColor, 'TDepth> =
     let size = 2 * int (ceil (4.0 * standardDeviation)) + 1
     img.SmoothGaussian(size, size, standardDeviation, standardDeviation)
 let gaussianFilter (img : Image<'TColor, 'TDepth>) (standardDeviation : float) : Image<'TColor, 'TDepth> =
     let size = 2 * int (ceil (4.0 * standardDeviation)) + 1
     img.SmoothGaussian(size, size, standardDeviation, standardDeviation)


@@ -32,7 +150,6 @@ let drawPoints (img: Image<Gray, byte>) (points: Points) (intensity: byte) =
     for p in points do
         img.Data.[p.Y, p.X, 0] <- intensity
 
     for p in points do
         img.Data.[p.Y, p.X, 0] <- intensity
 

-

 type ExtremumType =
     | Maxima = 1
     | Minima = 2
 type ExtremumType =
     | Maxima = 1
     | Minima = 2


@@ -356,6 +473,7 @@ let areaOpen (img: Image<Gray, byte>) (area: int) =
 let areaClose (img: Image<Gray, byte>) (area: int) =
     areaOperation img area AreaOperation.Closing
 
 let areaClose (img: Image<Gray, byte>) (area: int) =
     areaOperation img area AreaOperation.Closing
 

+// A simpler algorithm than 'areaOpen' but slower.

 let areaOpen2 (img: Image<Gray, byte>) (area: int) =
     let w = img.Width
     let h = img.Height
 let areaOpen2 (img: Image<Gray, byte>) (area: int) =
     let w = img.Width
     let h = img.Height


@@ -537,21 +655,14 @@ let connectedComponents (img: Image<Gray, byte>) (startPoints: List<Point>) : Li
     List<Point>(pointChecked)
 
 
     List<Point>(pointChecked)
 
 

-let saveImg (img: Image<'TColor, 'TDepth>) (filepath: string) =
-    img.Save(filepath)
-
-
-let saveMat (mat: Matrix<'TDepth>) (filepath: string) =
-    use img = new Image<Gray, 'TDeph>(mat.Size)
-    mat.CopyTo(img)
-    saveImg img filepath
-

 let drawLine (img: Image<'TColor, 'TDepth>) (color: 'TColor) (x0: int) (y0: int) (x1: int) (y1: int) (thickness: int) =
     img.Draw(LineSegment2D(Point(x0, y0), Point(x1, y1)), color, thickness);
 
 let drawLine (img: Image<'TColor, 'TDepth>) (color: 'TColor) (x0: int) (y0: int) (x1: int) (y1: int) (thickness: int) =
     img.Draw(LineSegment2D(Point(x0, y0), Point(x1, y1)), color, thickness);
 

+

 let drawLineF (img: Image<'TColor, 'TDepth>) (color: 'TColor) (x0: float) (y0: float) (x1: float) (y1: float) (thickness: int) =
     img.Draw(LineSegment2DF(PointF(float32 x0, float32 y0), PointF(float32 x1, float32 y1)), color, thickness, CvEnum.LineType.AntiAlias);
 
 let drawLineF (img: Image<'TColor, 'TDepth>) (color: 'TColor) (x0: float) (y0: float) (x1: float) (y1: float) (thickness: int) =
     img.Draw(LineSegment2DF(PointF(float32 x0, float32 y0), PointF(float32 x1, float32 y1)), color, thickness, CvEnum.LineType.AntiAlias);
 

+

 let drawEllipse (img: Image<'TColor, 'TDepth>) (e: Types.Ellipse) (color: 'TColor) (alpha: float) =
 
     if alpha >= 1.0
 let drawEllipse (img: Image<'TColor, 'TDepth>) (e: Types.Ellipse) (color: 'TColor) (alpha: float) =
 
     if alpha >= 1.0

diff --git a/Parasitemia/Parasitemia/MainAnalysis.fs b/Parasitemia/Parasitemia/MainAnalysis.fs
index c987e63..6111f06 100644 (file)
--- a/Parasitemia/Parasitemia/MainAnalysis.fs
+++ b/Parasitemia/Parasitemia/MainAnalysis.fs
@@ -47,32 +47,7 @@ let doAnalysis (img: Image<Bgr, byte>) (name: string) (config: Config) : Cell li
     let filteredGreenWhitoutInfectionFloat = filteredGreenWhitoutInfection.Convert<Gray, float32>()
     let filteredGreenWhitoutStainFloat = filteredGreenWhitoutStain.Convert<Gray, float32>()
 
     let filteredGreenWhitoutInfectionFloat = filteredGreenWhitoutInfection.Convert<Gray, float32>()
     let filteredGreenWhitoutStainFloat = filteredGreenWhitoutStain.Convert<Gray, float32>()
 

-    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 = filteredGreenWhitoutStainFloat.Convolution(sobelKernel).Convert<Gray, float>()
-    use yEdges = filteredGreenWhitoutStainFloat.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 min = ref 0.0

     let minLocation = ref <| Point()
     let max = ref 0.0
     let maxLocation = ref <| Point()
     let minLocation = ref <| Point()
     let max = ref 0.0
     let maxLocation = ref <| Point()


@@ -82,11 +57,14 @@ let doAnalysis (img: Image<Bgr, byte>) (name: string) (config: Config) : Cell li
     use edges = new Matrix<byte>(xEdges.Size)
     CvInvoke.Threshold(magnitudesByte, edges, 0.0, 1.0, CvEnum.ThresholdType.Otsu ||| CvEnum.ThresholdType.Binary) |> ignore
 
     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 "Finding edges" (fun() -> thin edges)*)
+
+    let edges, xGradient, yGradient = ImgTools.findEdges filteredGreenWhitoutStainFloat
+

     logTime "Removing small connected components from thinning" (fun () -> removeArea edges 12)
 
     let allEllipses, ellipses = logTime "Finding ellipses" (fun () ->
     logTime "Removing small connected components from thinning" (fun () -> removeArea edges 12)
 
     let allEllipses, ellipses = logTime "Finding ellipses" (fun () ->

-        let matchingEllipses = Ellipse.find edges xEdges yEdges config
+        let matchingEllipses = Ellipse.find edges xGradient yGradient config

         matchingEllipses.Ellipses, matchingEllipses.PrunedEllipses )
 
     let cells = logTime "Classifier" (fun () -> Classifier.findCells ellipses parasites filteredGreenWhitoutStain config)
         matchingEllipses.Ellipses, matchingEllipses.PrunedEllipses )
 
     let cells = logTime "Classifier" (fun () -> Classifier.findCells ellipses parasites filteredGreenWhitoutStain config)


@@ -106,7 +84,7 @@ let doAnalysis (img: Image<Bgr, byte>) (name: string) (config: Config) : Cell li
         saveImg parasites.infection (buildFileName " - parasites - infection.png")
 
         let imgAllEllipses = img.Copy()
         saveImg parasites.infection (buildFileName " - parasites - infection.png")
 
         let imgAllEllipses = img.Copy()

-        drawEllipses imgAllEllipses allEllipses (Bgr(0.0, 240.0, 240.0)) 0.1
+        drawEllipses imgAllEllipses allEllipses (Bgr(0.0, 240.0, 240.0)) 0.05

         saveImg imgAllEllipses (buildFileName " - ellipses - all.png")
 
         let imgEllipses = img.Copy()
         saveImg imgAllEllipses (buildFileName " - ellipses - all.png")
 
         let imgEllipses = img.Copy()

diff --git a/Parasitemia/Parasitemia/MatchingEllipses.fs b/Parasitemia/Parasitemia/MatchingEllipses.fs
index 6e6218d..c39541e 100644 (file)
--- a/Parasitemia/Parasitemia/MatchingEllipses.fs
+++ b/Parasitemia/Parasitemia/MatchingEllipses.fs
@@ -12,8 +12,7 @@ open Utils
 // Do not take in account matching score below this when two ellipses are matched.
 let matchingScoreThreshold1 = 0.6
 
 // Do not take in account matching score below this when two ellipses are matched.
 let matchingScoreThreshold1 = 0.6
 

-// All ellipses with a score below this are removed.
-let matchingScoreThreshold2 = 600.
+let scaleOverlapTest = 0.8

 
 type private EllipseScoreFlaggedKd (matchingScore: float, e: Ellipse) =
     let mutable matchingScore = matchingScore
 
 type private EllipseScoreFlaggedKd (matchingScore: float, e: Ellipse) =
     let mutable matchingScore = matchingScore


@@ -36,6 +35,9 @@ type private EllipseScoreFlaggedKd (matchingScore: float, e: Ellipse) =
 type MatchingEllipses (radiusMin: float) =
     let ellipses = List<EllipseScoreFlaggedKd>()
 
 type MatchingEllipses (radiusMin: float) =
     let ellipses = List<EllipseScoreFlaggedKd>()
 

+    // All ellipses with a score below this are removed.
+    let matchingScoreThreshold2 = 25. * radiusMin // 600.
+

     member this.Add (e: Ellipse) =
         ellipses.Add(EllipseScoreFlaggedKd(0.0, e))
 
     member this.Add (e: Ellipse) =
         ellipses.Add(EllipseScoreFlaggedKd(0.0, e))
 


@@ -76,7 +78,7 @@ type MatchingEllipses (radiusMin: float) =
             // 3) Sort ellipses by their score.
             ellipses.Sort(fun e1 e2 -> e2.MatchingScore.CompareTo(e1.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.
+            // 4) Remove ellipses with 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 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


@@ -98,7 +100,7 @@ type MatchingEllipses (radiusMin: float) =
                     for other in tree.Search window do
                         if not other.Removed && other.MatchingScore < e.MatchingScore
                         then
                     for other in tree.Search window do
                         if not other.Removed && other.MatchingScore < e.MatchingScore
                         then

-                            if e.Ellipse.Scale(0.8).Contains other.Ellipse.Cx other.Ellipse.Cy
+                            if e.Ellipse.Scale(scaleOverlapTest).Contains other.Ellipse.Cx other.Ellipse.Cy

                             then
                                 other.Removed <- true
             ellipses.RemoveAll(fun e -> e.Removed) |> ignore
                             then
                                 other.Removed <- true
             ellipses.RemoveAll(fun e -> e.Removed) |> ignore

diff --git a/Parasitemia/Parasitemia/Program.fs b/Parasitemia/Parasitemia/Program.fs
index 528f2de..dafff56 100644 (file)
--- a/Parasitemia/Parasitemia/Program.fs
+++ b/Parasitemia/Parasitemia/Program.fs
@@ -69,7 +69,7 @@ let main args =
                 factorNbPick = 1.0
                 factorWindowSize = 2.0
 
                 factorNbPick = 1.0
                 factorWindowSize = 2.0
 

-                darkStainLevel = 0.3 // Lower -> more sensitive.
+                darkStainLevel = 0.25 // Lower -> more sensitive. 0.3

                 maxDarkStainRatio = 0.1
 
                 infectionArea = 0.012 // 1.2 %
                 maxDarkStainRatio = 0.1
 
                 infectionArea = 0.012 // 1.2 %


@@ -79,7 +79,7 @@ let main args =
                 stainLevel = 1.1 // Lower -> more sensitive.
                 maxStainRatio = 0.12 // 12 %
 
                 stainLevel = 1.1 // Lower -> more sensitive.
                 maxStainRatio = 0.12 // 12 %
 

-                standardDeviationMaxRatio = 0.55
+                standardDeviationMaxRatio = 0.58 // 0.55

                 minimumCellArea = 0.5 })
 
         match mode with
                 minimumCellArea = 0.5 })
 
         match mode with

diff --git a/Parasitemia/Parasitemia/Utils.fs b/Parasitemia/Parasitemia/Utils.fs
index d6fa2f0..58e914f 100644 (file)
--- a/Parasitemia/Parasitemia/Utils.fs
+++ b/Parasitemia/Parasitemia/Utils.fs
@@ -30,6 +30,10 @@ let inline pointFromTwoLines (l1: Line) (l2: Line) : PointD =
     let y = -(l2.A * l1.B - l1.A * l2.B) / (l1.A - l2.A)
     PointD(x, y)
 
     let y = -(l2.A * l1.B - l1.A * l2.B) / (l1.A - l2.A)
     PointD(x, y)
 

+let inline linePassThroughSegment (l: Line) (p1: PointD) (p2: PointD) : bool =
+    let p = pointFromTwoLines l (lineFromTwoPoints p1 p2)
+    sign (p.X - p1.X) <> sign (p.X - p2.X)
+

 let inline squaredDistanceTwoPoints (p1: PointD) (p2: PointD) =
     (p1.X - p2.X) ** 2.0 + (p1.Y - p2.Y) ** 2.0
 
 let inline squaredDistanceTwoPoints (p1: PointD) (p2: PointD) =
     (p1.X - p2.X) ** 2.0 + (p1.Y - p2.Y) ** 2.0