The main process is now complete.

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

diff --git a/Parasitemia/Parasitemia/Ellipse.fs b/Parasitemia/Parasitemia/Ellipse.fs
index 62273b5..62c0936 100644 (file)
--- a/Parasitemia/Parasitemia/Ellipse.fs
+++ b/Parasitemia/Parasitemia/Ellipse.fs
@@ -31,20 +31,22 @@ let private goldenSectionSearch (f: float -> float) (nbIter: int) (xmin: float)
         
         if fc < fd
         then
-            b <- d;
-            d <- c;
-            c <- b - gr * (b - a);
+            b <- d
+            d <- c
+            c <- b - gr * (b - a)
         else
-            a <- c;
-            c <- d;
-            d <- a + gr * (b - a);
+            a <- c
+            c <- d
+            d <- a + gr * (b - a)
     
     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;
+    let accuracy_extremum_search_1 = 4
+    let accuracy_extremum_search_2 = 3
 
     // p3 as the referencial.
     let p1x = p1x - p3x
@@ -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.
@@ -139,7 +141,7 @@ let ellipse (p1x: float) (p1y: float) (m1: float) (p2x: float) (p2y: float) (m2:
         let cx = cx + p3x
         let cy = cy + p3y
 
-        Some { cx = cx; cy = cy; a = r1e; b = r2e; alpha = alpha }
+        Some (Types.Ellipse(cx, cy, r1e, r2e, alpha))
     else
         None
 
@@ -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>) 
@@ -195,13 +196,13 @@ let find (edges: Matrix<byte>)
          (windowSize: int) 
          (factorNbPick: float) : Types.Ellipse list =
 
-    let increment = windowSize / 4;
+    let increment = windowSize / 4
 
     let r1, r2 = radiusRange
     let radiusTolerance = (r2 - r1) * 0.2
 
-    let minimumDistance = (r2 / 1.5) ** 2.0;
-    let squaredDistance x1 y1 x2 y2 = (x1 - x2) ** 2.0 + (y1 - y2) ** 2.0;
+    let minimumDistance = (r2 / 1.5) ** 2.0
+    let squaredDistance x1 y1 x2 y2 = (x1 - x2) ** 2.0 + (y1 - y2) ** 2.0
 
     let h = edges.Height
     let w = edges.Width
@@ -214,9 +215,9 @@ let find (edges: Matrix<byte>)
     let xDirData = xDir.Data
     let yDirData = yDir.Data
 
-    let rng = Random()
+    let rng = Random(42)
     
-    let ellipses = MatchingEllipses ()
+    let ellipses = MatchingEllipses(r1)
         
     for window_i in -windowSize + increment .. increment .. h - increment do                                                  
         for window_j in -windowSize + increment .. increment .. w - increment do            
@@ -257,8 +258,8 @@ let find (edges: Matrix<byte>)
                             match areVectorsValid p1xf p1yf p2xf p2yf -xDirData.[p1y, p1x, 0] -yDirData.[p1y, p1x, 0] -xDirData.[p2y, p2x, 0] -yDirData.[p2y, p2x, 0] with
                             | Some (m1, m2) ->
                                 match ellipse p1xf p1yf m1 p2xf p2yf m2 p3xf p3yf with
-                                | Some e when e.cx > 0.0 && e.cx < (float w) - 1.0 && e.cy > 0.0 && e.cy < (float h) - 1.0 && 
-                                              e.a >= r1 - radiusTolerance && e.a <= r2 + radiusTolerance && e.b >= r1 - radiusTolerance && e.b <= r2 + radiusTolerance ->
+                                | Some e when e.Cx > 0.0 && e.Cx < (float w) - 1.0 && e.Cy > 0.0 && e.Cy < (float h) - 1.0 && 
+                                              e.A >= r1 - radiusTolerance && e.A <= r2 + radiusTolerance && e.B >= r1 - radiusTolerance && e.B <= r2 + radiusTolerance ->
                                      ellipses.Add e
                                 | _ -> ()
                             | _ -> ()