Remove the removing of cells bases on standard deviation.

[master-thesis.git] / Parasitemia / ParasitemiaCore / Ellipse.fs

diff --git a/Parasitemia/ParasitemiaCore/Ellipse.fs b/Parasitemia/ParasitemiaCore/Ellipse.fs
index 4e48230..01e1bec 100644 (file)
--- a/Parasitemia/ParasitemiaCore/Ellipse.fs
+++ b/Parasitemia/ParasitemiaCore/Ellipse.fs
@@ -14,6 +14,9 @@ open Config
 open MatchingEllipses
 open Const
 
+// This is a ratio of the biggest radius.
+let minimumDistanceBetweenDrawnPoints = 0.6
+
 /// <summary>
 /// Try to build an ellipse from three points and two tangents passing by the first and the second point.
 /// 'Ellipse.A' is always equal or greater than Ellipse.B.
@@ -74,30 +77,18 @@ let ellipse (p1x: float) (p1y: float) (m1: float) (p2x: float) (p2y: float) (m2:
 
             Some (Types.Ellipse(float32 cx, float32 cy, float32 majorAxis, float32 minorAxis, float32 phi'))
 
-
-let private vectorRotation (p1x: float32) (p1y: float32) (v1x: float32) (v1y: float32) (px: float32) (py: float32) : float32 =
-    let mutable rotation = 1.f
-    if p1y > py
-    then
-        if v1x > 0.f
-        then
-            rotation <- -1.f
-    elif p1y < py
+let inline private vectorRotation (px: float32) (py: float32) (vx: float32) (vy: float32) (p0x: float32) (p0y: float32) : float32 =
+    if py > p0y
     then
-        if v1x < 0.f
-        then
-            rotation <- -1.f
-    elif p1x > px
+        if vx > 0.f then -1.f else 1.f
+    elif py < p0y
     then
-        if v1y < 0.f
-        then
-            rotation <- -1.f
-    elif p1x < px
+        if vx < 0.f then -1.f else 1.f
+    elif px > p0x
     then
-        if v1y > 0.f
-        then
-            rotation <- -1.f
-    rotation
+        if vy < 0.f then -1.f else 1.f
+    else // p1x < px
+        if vy > 0.f then -1.f else 1.f
 
 let private areVectorsValid (p1x: float32) (p1y: float32) (p2x: float32) (p2y: float32) (v1x: float32) (v1y: float32) (v2x: float32) (v2y: float32) : (float32 * float32) option =
     let m1 = -v1x / v1y
@@ -105,30 +96,35 @@ let private areVectorsValid (p1x: float32) (p1y: float32) (p2x: float32) (p2y: f
 
     let b1 = -m1 * p1x + p1y
     let b2 = -m2 * p2x + p2y
-    let px = -((b1 - b2) / (m1 - m2))
-    let py = -((m2 * b1 - m1 * b2) / (m1 - m2))
+    let p0x = -(b1 - b2) / (m1 - m2)
+    let p0y = -(m2 * b1 - m1 * b2) / (m1 - m2)
 
-    let rot1 = vectorRotation p1x p1y v1x v1y px py
-    let rot2 = vectorRotation p2x p2y v2x v2y px py
+    let rot1 = vectorRotation p1x p1y v1x v1y p0x p0y
+    let rot2 = vectorRotation p2x p2y v2x v2y p0x p0y
 
     if rot1 = rot2
     then
         None
     else
-        let alpha1 = atan2 (p1y - py) (p1x - px)
-        let alpha2 = atan2 (p2y - py) (p2x - px)
+        let alpha1 =
+            let alpha1' = atan2 (p1y - p0y) (p1x - p0x)
+            if alpha1' < 0.f then 2.f * PI + alpha1' else alpha1'
 
-        let alpha1' = if alpha1 < 0.f then 2.f * PI + alpha1 else alpha1
-        let alpha2' = if alpha2 < 0.f then 2.f * PI + alpha2 else alpha2
+        let alpha2 =
+            let alpha2' = atan2 (p2y - p0y) (p2x - p0x)
+            if alpha2' < 0.f then 2.f * PI + alpha2' else alpha2'
 
-        let diff = rot1 * alpha1' + rot2 * alpha2'
+        let diff = rot1 * alpha1 + rot2 * alpha2
 
         if diff > PI || (diff < 0.f && diff > -PI)
         then
-            None
+            Some (m1, m2)
         else
-        Some (m1, m2)
+            None
 
+/// <summary>
+/// Build a set of ellipses as a 'MatchingEllipses' object by finding ellipses with the given edges and gradient.
+/// </summary>
 let find (edges: Matrix<byte>)
          (xGradient: Matrix<float32>)
          (yGradient: Matrix<float32>)
@@ -139,13 +135,15 @@ let find (edges: Matrix<byte>)
 
     // We choose a window size for which the biggest ellipse can always be fitted in.
     let windowSize = roundInt (2.f * r2)
-    let factorNbPick = config.Parameters.factorNbPick
+    let factorNbValidPick = config.Parameters.factorNbValidPick
+    let factorNbMaxPick = config.Parameters.factorNbMaxPick
+    let nbPickElementsMin = config.Parameters.nbPickElementsMin
 
     let increment = windowSize / (int incrementWindowDivisor)
 
     let radiusTolerance = (r2 - r1) * 0.2f
 
-    let squaredMinimumDistance = (float r2 / 1.5) ** 2.
+    let squaredMinimumDistance = (float config.RBCRadius.Pixel * minimumDistanceBetweenDrawnPoints) ** 2.
     let inline squaredDistance x1 y1 x2 y2 = (x1 - x2) ** 2. + (y1 - y2) ** 2.
 
     let h = edges.Height
@@ -186,10 +184,11 @@ let find (edges: Matrix<byte>)
                     if edgesData.[i, j] = 1uy
                     then currentElements.Add(Point(j, i))
 
-            if currentElements.Count >= 10
+            if currentElements.Count >= nbPickElementsMin
             then
-                let mutable nbOfPicks = (float currentElements.Count) * factorNbPick |> int
-                while nbOfPicks > 0 do
+                let mutable nbOfPicks = (float currentElements.Count) * factorNbMaxPick |> int
+                let mutable nbOfValidPicks = (float currentElements.Count) * factorNbValidPick |> int
+                while nbOfPicks > 0 && nbOfValidPicks > 0 do
                     let p1 = currentElements.[rng.Next(currentElements.Count)]
                     let p2 = currentElements.[rng.Next(currentElements.Count)]
                     let p3 = currentElements.[rng.Next(currentElements.Count)]
@@ -206,9 +205,10 @@ let find (edges: Matrix<byte>)
                             match areVectorsValid (float32 p1xf) (float32 p1yf) (float32 p2xf) (float32 p2yf) -xDirData.[p1.Y, p1.X] -yDirData.[p1.Y, p1.X] -xDirData.[p2.Y, p2.X] -yDirData.[p2.Y, p2.X] with
                             | Some (m1, m2) ->
                                 match ellipse p1xf p1yf (float m1) p2xf p2yf (float m2) p3xf p3yf with
-                                | Some e when e.Cx > 0.f && e.Cx < w_f - 1.f && e.Cy > 0.f && e.Cy < h_f - 1.f &&
+                                | Some e when e.Cx > 0.f && e.Cx <  w_f - 1.f && e.Cy > 0.f && e.Cy < h_f - 1.f &&
                                               e.A >= r1 - radiusTolerance && e.A <= r2 + radiusTolerance && e.B >= r1 - radiusTolerance && e.B <= r2 + radiusTolerance ->
-                                     ellipses.Add e
+                                    nbOfValidPicks <- nbOfValidPicks - 1
+                                    ellipses.Add e
                                 | _ -> ()
                             | _ -> ()