X-Git-Url: http://git.euphorik.ch/?p=master-thesis.git;a=blobdiff_plain;f=Parasitemia%2FParasitemia%2FEEOver.fs;h=1a6bb5fbac723a7b27282e5ed69407c242154141;hp=52b59636eeadbe631cc74dcb44b706254e56deab;hb=5ac2dedf8ead5275ac216e0b41829ab39c843800;hpb=b070295cf67b2025164a34b6594e84f0d771cdc9

diff --git a/Parasitemia/Parasitemia/EEOver.fs b/Parasitemia/Parasitemia/EEOver.fs
index 52b5963..1a6bb5f 100644
--- a/Parasitemia/Parasitemia/EEOver.fs
+++ b/Parasitemia/Parasitemia/EEOver.fs
@@ -7,7 +7,7 @@ let private EPS = 1.0e-5
 let inline private ellipse2tr (x: float) (y: float) (aa: float) (bb: float) (cc: float) (dd: float) (ee: float) (ff: float) : float =
     aa * x * x + bb * x * y + cc * y * y + dd * x + ee * y + ff
 
-let private nointpts (a1: float) (b1: float) (a2: float) (b2: float) (h1: float) (k1: float) (h2: float) (k2: float) (phi_1: float) (phi_2: float) (h2_tr: float) (k2_tr: float) (aa: float) (bb: float) (cc: float) (dd: float) (ee: float) (ff: float) = 
+let private nointpts (a1: float) (b1: float) (a2: float) (b2: float) (h1: float) (k1: float) (h2: float) (k2: float) (phi_1: float) (phi_2: float) (h2_tr: float) (k2_tr: float) (aa: float) (bb: float) (cc: float) (dd: float) (ee: float) (ff: float) =
     let a1b1 = a1 * b1
     let a2b2 = a2 * b2
     let area_1 = Math.PI * a1b1
@@ -34,13 +34,13 @@ let private nointpts (a1: float) (b1: float) (a2: float) (b2: float) (h1: float)
 type private PointType = TANGENT_POINT | INTERSECTION_POINT
 
 let private istanpt (x: float) (y: float) (a1: float) (b1: float) (aa: float) (bb: float) (cc: float) (dd: float) (ee: float) (ff: float) : PointType =
-    let x = 
+    let x =
         if abs x > a1
         then
             if x < 0.0 then -a1 else a1
         else x
 
-    let theta = 
+    let theta =
         if y < 0.0
         then 2.0 * Math.PI - acos (x / a1)
         else acos (x / a1)
@@ -105,7 +105,7 @@ let private twointpts (x: float[]) (y: float[]) (a1: float) (b1: float) (phi_1:
     if theta1 > theta2
     then
         theta1 <- theta1 - 2.0 * Math.PI
-    
+
     let trsign = if (theta2 - theta1) > Math.PI then 1.0 else -1.0
 
     let mutable area1 = 0.5 * (a1 * b1 * (theta2 - theta1) + trsign * abs (x.[0] * y.[1] - x.[1] * y.[0]))
@@ -182,7 +182,7 @@ let private twointpts (x: float[]) (y: float[]) (a1: float) (b1: float) (phi_1:
     area1 + area2
 
 
-let private threeintpts (xint: float[]) (yint: float[]) (a1: float) (b1: float) (phi_1: float) (a2: float) (b2: float) (h2_tr: float) (k2_tr: float) (phi_2: float) (aa: float) (bb: float) (cc: float) (dd: float) (ee: float) (ff: float) : float = 
+let private threeintpts (xint: float[]) (yint: float[]) (a1: float) (b1: float) (phi_1: float) (a2: float) (b2: float) (h2_tr: float) (k2_tr: float) (phi_2: float) (aa: float) (bb: float) (cc: float) (dd: float) (ee: float) (ff: float) : float =
     let mutable tanpts = 0
     let mutable tanindex = 0
     for i in 0..2 do
@@ -205,7 +205,7 @@ let private threeintpts (xint: float[]) (yint: float[]) (a1: float) (b1: float)
         | 1 ->
             xint.[1] <- xint.[2]
             yint.[1] <- yint.[2]
-        | _ -> 
+        | _ ->
             ()
         twointpts xint yint a1 b1 phi_1 a2 b2 h2_tr k2_tr phi_2 aa bb cc dd ee ff
 
@@ -244,8 +244,8 @@ let private fourintpts (xint: float[]) (yint: float[]) (a1: float) (b1: float) (
             else
                 theta.[k+1] <- theta.[k]
                 xint.[k+1] <- xint.[k]
-                yint.[k+1] <- yint.[k]            
-                k <- k - 1                
+                yint.[k+1] <- yint.[k]
+                k <- k - 1
                 k2 <- k + 1
 
         theta.[k2] <- tmp0
@@ -258,7 +258,7 @@ let private fourintpts (xint: float[]) (yint: float[]) (a1: float) (b1: float) (
     for k in 0..3 do
         printf "k=%d:  Theta = %f, xint=%f, yint=%f\n" k theta.[k] xint.[k] yint.[k]
 #endif
-                
+
     let area1 = 0.5 * abs ((xint.[2] - xint.[0]) * (yint.[3] - yint.[1]) - (xint.[3] - xint.[1]) * (yint.[2] - yint.[0]))
 
     let cosphi = cos (phi_1 - phi_2)
@@ -275,7 +275,7 @@ let private fourintpts (xint: float[]) (yint: float[]) (a1: float) (b1: float) (
         if abs xint_tr.[i] > a2
         then
             xint_tr.[i] <- if xint_tr.[i] < 0.0 then -a2 else a2
-        
+
         theta_tr.[i] <- if yint_tr.[i] < 0.0 then 2.0 * Math.PI - acos (xint_tr.[i] / a2) else acos (xint_tr.[i] / a2)
 
     let xmid = a1 * cos ((theta.[0] + theta.[1]) / 2.0)
@@ -298,9 +298,9 @@ let private fourintpts (xint: float[]) (yint: float[]) (a1: float) (b1: float) (
         area3 <- 0.5 * (a1b1 * (theta.[0] - (theta.[3] - 2.0 * Math.PI)) - abs (xint.[3] * yint.[0] - xint.[0] * yint.[3]))
         area4 <- 0.5 * (a2b2 * (theta_tr.[1] - theta_tr.[0]) - abs (xint_tr.[0] * yint_tr.[1] - xint_tr.[1] * yint_tr.[0]))
         area5 <- 0.5 * (a2b2 * (theta_tr.[3] - theta_tr.[2]) - abs (xint_tr.[2] * yint_tr.[3] - xint_tr.[3] * yint_tr.[2]))
-        
+
     if area5 < 0.0
-    then 
+    then
 #if DEBUG_LOG
         printf "\n\t\t-------------> area5 is negativ (%f). Add: pi*A2*B2=%f <------------\n" area5 area_2
 #endif
@@ -321,7 +321,7 @@ let private fourintpts (xint: float[]) (yint: float[]) (a1: float) (b1: float) (
         area3 <- area3 + area_1
 
     if area2 < 0.0
-    then    
+    then
 #if DEBUG_LOG
         printf "\n\t\t-------------> area2 is negativ (%f). Add: pi*A2*B2=%f <------------\n" area2 area_1
 #endif
@@ -334,7 +334,7 @@ let private fourintpts (xint: float[]) (yint: float[]) (a1: float) (b1: float) (
     area1 + area2 + area3 + area4 + area5
 
 
-let private quadroots (p: float[]) (r: float[,]) = 
+let private quadroots (p: float[]) (r: float[,]) =
     let mutable b = -p.[1] / (2.0 * p.[0])
     let c = p.[2] / p.[0]
     let mutable d = b * b - c
@@ -342,7 +342,7 @@ let private quadroots (p: float[]) (r: float[,]) =
     if d >= 0.0
     then
         if b > 0.0
-        then 
+        then
             b <- sqrt d + b
             r.[1, 2] <- b
         else
@@ -370,7 +370,7 @@ let private cubicroots (p: float[]) (r: float[,]) =
     let mutable c = t * t * t
     let mutable d = b * b - c
 
-    if d >= 0.0 
+    if d >= 0.0
     then
         d <- ((sqrt d) + (abs b)) ** (1.0 / 3.0)
         if d <> 0.0
@@ -427,7 +427,7 @@ let private cubicroots (p: float[]) (r: float[,]) =
         r.[1, 1] <- t
         for k in 1..3 do
             r.[2, k] <- 0.0
-        
+
 
 let private biquadroots (p: float[]) (r: float[,]) =
     if p.[0] <> 1.0
@@ -445,7 +445,7 @@ let private biquadroots (p: float[]) (r: float[,]) =
     a <- a - d
 
     let quadExecuted = ref false
-    let quad () = 
+    let quad () =
         if not !quadExecuted
         then
             p.[2] <- !c / !b
@@ -482,7 +482,7 @@ let private biquadroots (p: float[]) (r: float[,]) =
                     p.[1] <- -(sqrt d)
                 b := 0.5 * (a + !b / p.[1])
                 quad ()
-                k <- 4            
+                k <- 4
             k <- k + 1
 
     if not !quadExecuted && p.[2] < 0.0
@@ -517,17 +517,17 @@ let EEOverlapArea (e1: Types.Ellipse) (e2: Types.Ellipse) : (float * float[] * f
     let h1, k1, a1, b1, phi_1 = e1.Cx, e1.Cy, e1.A, e1.B, e1.Alpha
     let h2, k2, a2, b2, phi_2 = e2.Cx, e2.Cy, e2.A, e2.B, e2.Alpha
 
-    if a1 <= EPS || b1 <= EPS || a2 <= EPS || b2 <= EPS 
+    if a1 <= EPS || b1 <= EPS || a2 <= EPS || b2 <= EPS
     then
         None
     else
-        let phi_1 = phi_1 % Math.PI
-        let phi_2 = phi_2 % Math.PI
+        let phi_1 = phi_1 % Math.PI //(if phi_1 > Math.PI / 2.0 then phi_1 - Math.PI else phi_1) % Math.PI
+        let phi_2 = phi_2 % Math.PI //(if phi_2 > Math.PI / 2.0 then phi_2 - Math.PI else phi_2) % Math.PI
         let h2_tr, k2_tr, phi_2r =
             let cosphi = cos phi_1
             let sinphi = sin phi_1
             (h2 - h1) * cosphi + (k2 - k1) * sinphi, (h1 - h2) * sinphi + (k2 - k1) * cosphi, (phi_2 - phi_1) % (2.0 * Math.PI)
-        
+
 #if DEBUG_LOG
         printf "H2_TR=%f, K2_TR=%f, PHI_2R=%f\n" h2_tr k2_tr phi_2r
 #endif
@@ -562,7 +562,7 @@ let EEOverlapArea (e1: Types.Ellipse) (e2: Types.Ellipse) : (float * float[] * f
 #if DEBUG_LOG
         for i in 0..4 do
             printf "cy[%d]=%f\n" i cy.[i]
-#endif	
+#endif
 
         let py = Array.zeroCreate<float> 5
         let r = Array2D.zeroCreate<float> 3 5
@@ -587,7 +587,7 @@ let EEOverlapArea (e1: Types.Ellipse) (e2: Types.Ellipse) : (float * float[] * f
                 py.[0] <- 1.0
                 cubicroots py r
                 3
-        
+
             elif abs cy.[2] > EPS
             then
                 for i in 0..1 do
@@ -652,7 +652,7 @@ let EEOverlapArea (e1: Types.Ellipse) (e2: Types.Ellipse) : (float * float[] * f
 #if DEBUG_LOG
                 printf "\tx1=%f, y1=%f, A=%f. B=%f ---> ellipse2tr(x1)= %f\n" x1 ychk.[i] a1 b1 (ellipse2tr x1 ychk.[i] aa bb cc dd ee ff)
                 printf "\tx2=%f, y1=%f, A=%f. B=%f ---> ellipse2tr(x2) %f\n" x2 ychk.[i] a1 b1 (ellipse2tr x2 ychk.[i] aa bb cc dd ee ff)
-#endif		
+#endif
 
                 if abs (ellipse2tr x1 ychk.[i] aa bb cc dd ee ff) < EPS
                 then
@@ -661,7 +661,7 @@ let EEOverlapArea (e1: Types.Ellipse) (e2: Types.Ellipse) : (float * float[] * f
                     printf "first if x1. acc nintps=%d\n" nintpts
 #endif
                     if nintpts > 4
-                    then 
+                    then
                         returnValue <- -1.0
                     else
                         xint.[nintpts-1] <- x1
@@ -677,7 +677,7 @@ let EEOverlapArea (e1: Types.Ellipse) (e2: Types.Ellipse) : (float * float[] * f
                     printf "first if x2. nintps=%d, Dx=%f (eps2=%f) \n" nintpts (abs (x2 - x1)) EPS
 #endif
                     if nintpts > 4
-                    then 
+                    then
                         returnValue <- -1.0
                     else
                         xint.[nintpts-1] <- x2
@@ -692,8 +692,8 @@ let EEOverlapArea (e1: Types.Ellipse) (e2: Types.Ellipse) : (float * float[] * f
                 printf "i=%d, multiple roots: %f  <--------> %f. continue\n" i ychk.[i] ychk.[i-1]
 #endif
             i <- i + 1
-                    
-        
+
+
         if returnValue = -1.0
         then
             None
@@ -702,13 +702,13 @@ let EEOverlapArea (e1: Types.Ellipse) (e2: Types.Ellipse) : (float * float[] * f
                 match nintpts with
                 | 0 | 1 -> nointpts a1 b1 a2 b2 h1 k1 h2 k2 phi_1 phi_2 h2_tr k2_tr aa bb cc dd ee ff
                 | 2 -> match istanpt xint.[0] yint.[0] a1 b1 aa bb cc dd ee ff with
-                       | TANGENT_POINT -> 
+                       | TANGENT_POINT ->
 #if DEBUG_LOG
                             printf "one point is tangent\n"
 #endif
                             nointpts a1 b1 a2 b2 h1 k1 h2 k2 phi_1 phi_2 h2_tr k2_tr aa bb cc dd ee ff
 
-                       | INTERSECTION_POINT -> 
+                       | INTERSECTION_POINT ->
 #if DEBUG_LOG
                             printf "check twointpts\n"
 #endif
@@ -718,7 +718,7 @@ let EEOverlapArea (e1: Types.Ellipse) (e2: Types.Ellipse) : (float * float[] * f
                 | _ -> -1.0
             if nintpts = 0
             then Some (area, [||], [||])
-            else 
+            else
                 let xTransform = Array.zeroCreate nintpts
                 let yTransform = Array.zeroCreate nintpts
                 for i in 0 .. (nintpts - 1) do