Peculiar elements overlap RBC now.
[master-thesis.git] / Parasitemia / ParasitemiaCore / Classifier.fs
index 0012a4f..ea3f9dd 100644 (file)
@@ -10,10 +10,12 @@ open Emgu.CV.Structure
 open Types
 open Utils
 
+type CellState = RBC = 1 | Removed = 2 | Peculiar = 3
+
 type private EllipseFlaggedKd (e: Ellipse) =
     inherit Ellipse (e.Cx, e.Cy, e.A, e.B, e.Alpha)
 
-    member val Removed = false with get, set
+    member val State = CellState.RBC with get, set
 
     interface KdTree.I2DCoords with
         member this.X = this.Cx
@@ -38,32 +40,45 @@ let findCells (ellipses: Ellipse list) (parasites: ParasitesMarker.Result) (widt
 
         // Return 'true' if the point 'p' is owned by e.
         // The lines represents all intersections with other ellipses.
-        let pixelOwnedByE (p: PointF) (e: Ellipse) (others: (Ellipse * Line) list) =
+        let pixelOwnedByE (p: PointF) (e: EllipseFlaggedKd) (neighbors: (EllipseFlaggedKd * PointF * PointF) list) =
             e.Contains p.X p.Y &&
             seq {
                 let c = PointF(e.Cx, e.Cy)
 
-                for e', d1 in others do
-                    let d2 = lineFromTwoPoints c p
-                    let c' = PointF(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 not (Single.IsInfinity d2.A)
+                for e', d1 in neighbors
+                              |> List.choose (fun (otherE, p1, p2) ->
+                                                if otherE.State = CellState.Removed
+                                                then None
+                                                else Some (otherE, Utils.lineFromTwoPoints p1 p2)) do
+                    if e'.State = e.State // Peculiar vs peculiar or RBC vs RBC.
                     then
-                        let p' = Utils.pointFromTwoLines d1 d2
-                        let delta, delta' =
-                            let dx1, dx2 = (c.X - p.X), (c.X - p'.X)
-                            // To avoid rounding issue.
-                            if abs dx1 < 0.01f || abs dx2 < 0.01f then c.Y - p.Y, c.Y - p'.Y else dx1, dx2
-
-                        // Yield 'false' when the point is owned by another ellipse.
-                        if case1
+                        let d2 = lineFromTwoPoints c p
+                        let c' = PointF(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 not (Single.IsInfinity d2.A)
                         then
-                            yield sign delta <> sign delta' || Utils.squaredDistanceTwoPoints c p' > Utils.squaredDistanceTwoPoints c p
+                            let p' = Utils.pointFromTwoLines d1 d2
+                            let delta, delta' =
+                                let dx1, dx2 = (c.X - p.X), (c.X - p'.X)
+                                // To avoid rounding issue.
+                                if abs dx1 < 0.01f || abs dx2 < 0.01f then c.Y - p.Y, c.Y - p'.Y else dx1, dx2
+
+                            // Yield 'false' when the point is owned by another ellipse.
+                            if case1
+                            then
+                                yield sign delta <> sign delta' || Utils.squaredDistanceTwoPoints c p' > Utils.squaredDistanceTwoPoints c p
+                            else
+                                yield sign delta = sign delta' && Utils.squaredDistanceTwoPoints c p' < Utils.squaredDistanceTwoPoints c p
                         else
-                            yield sign delta = sign delta' && Utils.squaredDistanceTwoPoints c p' < Utils.squaredDistanceTwoPoints c p
+                            yield case1
+
+                    elif e.State = CellState.Peculiar // A peculiar always win against a RBC.
+                    then
+                        yield true
                     else
-                        yield case1
+                        yield not <| e'.Contains p.X p.Y
+
             } |> Seq.forall id
 
         let ellipses = ellipses |> List.map EllipseFlaggedKd
@@ -71,7 +86,7 @@ let findCells (ellipses: Ellipse list) (parasites: ParasitesMarker.Result) (widt
         // 1) Associate touching ellipses with each ellipses and remove ellipse with more than two intersections.
         let tree = KdTree.Tree.BuildTree ellipses
         let neighbors (e: EllipseFlaggedKd) : (EllipseFlaggedKd * PointF * PointF) list =
-            if not e.Removed
+            if e.State <> CellState.Removed
             then
                 tree.Search (searchRegion e)
                     // We only keep the ellipses touching 'e'.
@@ -80,24 +95,23 @@ let findCells (ellipses: Ellipse list) (parasites: ParasitesMarker.Result) (widt
                         then
                             match EEOver.EEOverlapArea e otherE with
                             | Some (_, px, _) when px.Length > 2 ->
-                                otherE.Removed <- true
+                                otherE.State <- CellState.Removed
                                 None
                             | Some (area, px, py) when area > 0.f && px.Length = 2 ->
                                 Some (otherE, PointF(px.[0], py.[0]), PointF(px.[1], py.[1]))
                             | _ ->
                                 None
                         else
-                            None )
+                            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
+        let ellipsesWithNeigbors = ellipses |> List.map (fun e -> e, neighbors e)
 
         // 2) Remove ellipses touching the edges.
         let widthF, heightF = float32 width, float32 height
         for e in ellipses do
-            if e.isOutside widthF heightF then e.Removed <- true
+            if e.isOutside widthF heightF then e.State <- CellState.Removed
 
         // 3) Remove ellipses with a high standard deviation (high contrast).
         // Obsolete. It was useful when the ellipses result quality wasn't good.
@@ -126,7 +140,7 @@ let findCells (ellipses: Ellipse list) (parasites: ParasitesMarker.Result) (widt
         // 4) Remove ellipses with little area.
         let minArea = config.RBCRadius.MinArea
         for e, neighbors in ellipsesWithNeigbors do
-            if not e.Removed
+            if e.State <> CellState.Removed
             then
                 let minX, minY, maxX, maxY = ellipseWindow e
 
@@ -134,13 +148,37 @@ let findCells (ellipses: Ellipse list) (parasites: ParasitesMarker.Result) (widt
                 for y in (if minY < 0 then 0 else minY) .. (if maxY >= height then height - 1 else maxY) do
                     for x in (if minX < 0 then 0 else minX) .. (if maxX >= width then width - 1 else maxX) do
                         let p = PointF(float32 x, float32 y)
-                        if pixelOwnedByE p e (neighbors |> List.choose (fun (otherE, p1, p2) -> if otherE.Removed then None else Some (otherE :> Ellipse, Utils.lineFromTwoPoints p1 p2)))
+                        if pixelOwnedByE p e neighbors
                         then
                             area <- area + 1
 
                 if area < int minArea
                 then
-                    e.Removed <- true
+                    e.State <- CellState.Removed
+
+        // 5) Define non-rbc (peculiar) cells.
+        let darkStainData = parasites.darkStain.Data
+        ellipsesWithNeigbors
+        |> List.choose (fun (e, neighbors) ->
+            if e.State = CellState.Removed
+            then
+                None
+            else
+                let mutable darkStainPixels = 0
+                let mutable nbElement = 0
+                let minX, minY, maxX, maxY = ellipseWindow e
+                for y in minY .. maxY do
+                    for x in minX .. maxX do
+                        let p = PointF(float32 x, float32 y)
+                        if pixelOwnedByE p e neighbors
+                        then
+                            nbElement <- nbElement + 1
+                            if darkStainData.[y, x, 0] > 0uy
+                            then
+                                darkStainPixels <- darkStainPixels + 1
+
+                if float darkStainPixels > config.Parameters.maxDarkStainRatio * (float nbElement) then Some e else None)
+        |> List.iter (fun e -> e.State <- CellState.Peculiar)
 
         // 5) Define pixels associated to each ellipse and create the cells.
         let diameterParasiteSquared = (2.f * config.RBCRadius.ParasiteRadius) ** 2.f |> roundInt
@@ -152,7 +190,7 @@ let findCells (ellipses: Ellipse list) (parasites: ParasitesMarker.Result) (widt
 
         ellipsesWithNeigbors
         |> List.choose (fun (e, neighbors) ->
-            if e.Removed
+            if e.State = CellState.Removed
             then
                 None
             else
@@ -161,14 +199,13 @@ let findCells (ellipses: Ellipse list) (parasites: ParasitesMarker.Result) (widt
                 let nucleusPixels = List<Point>()
                 let parasitePixels = List<Point>()
 
-                let mutable darkStainPixels = 0
                 let mutable nbElement = 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 = PointF(float32 x, float32 y)
-                        if pixelOwnedByE p e (neighbors |> List.choose (fun (otherE, p1, p2) -> if otherE.Removed then None else Some (otherE :> Ellipse, Utils.lineFromTwoPoints p1 p2)))
+                        if pixelOwnedByE p e neighbors
                         then
                             elements.[y - minY, x - minX] <- 1uy
                             nbElement <- nbElement + 1
@@ -181,20 +218,18 @@ let findCells (ellipses: Ellipse list) (parasites: ParasitesMarker.Result) (widt
                             then
                                 parasitePixels.Add(Point(x, y))
 
-                            if darkStainData.[y, x, 0] > 0uy
-                            then
-                                darkStainPixels <- darkStainPixels + 1
-
-                let mutable parasiteArea = 0
-                if nucleusPixels.Count > 0
-                then
-                    for parasitePixel in parasitePixels do
-                        if nucleusPixels.Exists(fun p -> pown (p.X - parasitePixel.X) 2 + pown (p.Y - parasitePixel.Y) 2 <= diameterParasiteSquared)
-                        then
-                            parasiteArea <- parasiteArea + 1
+                let parasiteArea =
+                    if nucleusPixels.Count > 0
+                    then
+                        seq {
+                            for parasitePixel in parasitePixels do
+                                if nucleusPixels.Exists(fun p -> pown (p.X - parasitePixel.X) 2 + pown (p.Y - parasitePixel.Y) 2 <= diameterParasiteSquared)
+                                then yield 1 } |> Seq.sum
+                    else
+                        0
 
                 let cellClass =
-                    if float darkStainPixels > config.Parameters.maxDarkStainRatio * (float nbElement)
+                    if e.State = CellState.Peculiar
                     then
                         Peculiar