Use float32 to reduce memory footprint.

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

diff --git a/Parasitemia/Parasitemia/EEOver.fs b/Parasitemia/Parasitemia/EEOver.fs
index 1a6bb5f..5a3c9e4 100644 (file)
--- a/Parasitemia/Parasitemia/EEOver.fs
+++ b/Parasitemia/Parasitemia/EEOver.fs
@@ -286,6 +286,7 @@ let private fourintpts (xint: float[]) (yint: float[]) (a1: float) (b1: float) (
     if ellipse2tr xmid ymid aa bb cc dd ee ff < 0.0
     then
         area2 <- 0.5 * (a1b1 * (theta.[1] - theta.[0]) - abs (xint.[0] * yint.[1] - xint.[1] * yint.[0]))
+        area3 <- 0.5 * (a1b1 * (theta.[3] - theta.[2]) - abs (xint.[2] * yint.[3] - xint.[3] * yint.[2]))
         area4 <- 0.5 * (a2b2 * (theta_tr.[2] - theta_tr.[1]) - abs (xint_tr.[1] * yint_tr.[2] - xint_tr.[2] * yint_tr.[1]))
 
         if theta_tr.[3] > theta_tr.[0]
@@ -436,33 +437,33 @@ let private biquadroots (p: float[]) (r: float[,]) =
             p.[k] <- p.[k] / p.[0]
         p.[0] <- 1.0
     let e = 0.25 * p.[1]
-    let b = ref (2.0 * e)
-    let c = ref (!b ** 2.0)
-    let mutable d = 0.75 * !c
-    b := p.[3] + !b *(!c - p.[2])
+    let mutable b = 2.0 * e
+    let mutable c = b ** 2.0
+    let mutable d = 0.75 * c
+    b <- p.[3] + b *(c - p.[2])
     let mutable a = p.[2] - d
-    c := p.[4] + e * (e * a - p.[3])
+    c <- p.[4] + e * (e * a - p.[3])
     a <- a - d
 
-    let quadExecuted = ref false
+    let mutable quadExecuted = false
     let quad () =
-        if not !quadExecuted
+        if not quadExecuted
         then
-            p.[2] <- !c / !b
+            p.[2] <- c / b
             quadroots p r
             for k in 1..2 do
                 for j in 1..2 do
                     r.[j, k+2] <- r.[j, k]
             p.[1] <- -p.[1]
-            p.[2] <- !b
+            p.[2] <- b
             quadroots p r
             for k in 1..4 do
                 r.[1,k] <- r.[1,k] - e
-            quadExecuted := true
+            quadExecuted <- true
 
     p.[1] <- 0.5 * a
-    p.[2] <- (p.[1] * p.[1] - !c) * 0.25
-    p.[3] <- !b * !b / -64.0
+    p.[2] <- (p.[1] * p.[1] - c) * 0.25
+    p.[3] <- b * b / -64.0
     if p.[3] < 0.0
     then
         cubicroots p r
@@ -472,50 +473,50 @@ let private biquadroots (p: float[]) (r: float[,]) =
             then
                 d <- r.[1, k] * 4.0
                 a <- a + d
-                if a >= 0.0 && !b >= 0.0
+                if a >= 0.0 && b >= 0.0
                 then
                     p.[1] <- sqrt d
-                elif a <= 0.0 && !b <= 0.0
+                elif a <= 0.0 && b <= 0.0
                 then
                     p.[1] <- sqrt d
                 else
                     p.[1] <- -(sqrt d)
-                b := 0.5 * (a + !b / p.[1])
+                b <- 0.5 * (a + b / p.[1])
                 quad ()
                 k <- 4
             k <- k + 1
 
-    if not !quadExecuted && p.[2] < 0.0
+    if not quadExecuted && p.[2] < 0.0
     then
-        b := sqrt !c
-        d <- !b + !b - a
+        b <- sqrt c
+        d <- b + b - a
         p.[1] <- 0.0
         if d > 0.0
         then
             p.[1] <- sqrt d
-    elif not !quadExecuted
+    elif not quadExecuted
     then
         if p.[1] > 0.0
         then
-            b := (sqrt p.[2]) * 2.0 + p.[1]
+            b <- (sqrt p.[2]) * 2.0 + p.[1]
         else
-            b := -(sqrt p.[2]) * 2.0 + p.[1]
+            b <- -(sqrt p.[2]) * 2.0 + p.[1]
 
-        if !b <> 0.0
+        if b <> 0.0
         then
             p.[1] <- 0.0
         else
             for k in 1..4 do
                 r.[1, k] <- -e
                 r.[2, k] <- 0.0
-            quadExecuted := true
+            quadExecuted <- true
 
     quad ()
 
 // Return a tuple (area, x intersections, y intersections)
-let EEOverlapArea (e1: Types.Ellipse) (e2: Types.Ellipse) : (float * float[] * float[]) option =
-    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
+let EEOverlapArea (e1: Types.Ellipse) (e2: Types.Ellipse) : (float32 * float32[] * float32[]) option =
+    let h1, k1, a1, b1, phi_1 = float e1.Cx, float e1.Cy, float e1.A, float e1.B, float e1.Alpha
+    let h2, k2, a2, b2, phi_2 = float e2.Cx, float e2.Cy, float e2.A, float e2.B, float e2.Alpha
 
     if a1 <= EPS || b1 <= EPS || a2 <= EPS || b2 <= EPS
     then
@@ -609,15 +610,16 @@ let EEOverlapArea (e1: Types.Ellipse) (e2: Types.Ellipse) : (float * float[] * f
         printf "nroots = %d\n" nroots
 #endif
 
-        let ychk = [|
-            for i in 1 .. nroots do
-                if abs r.[2, i] < EPS
-                then
-                    yield r.[1, i] * b1
+        let ychk = Array.init nroots (fun _ -> Double.MaxValue)
+        let mutable nychk = 0
+        for i in 1 .. nroots do
+            if abs r.[2, i] < EPS
+            then
+                ychk.[nychk] <- r.[1, i] * b1
+                nychk <- nychk + 1
 #if DEBUG_LOG
-                    printf "ROOT is Real,  i=%d --> %f (B1=%f)\n" i r.[1, i] b1
+                printf "ROOT is Real,  i=%d --> %f (B1=%f)\n" i r.[1, i] b1
 #endif
-        |]
         Array.sortInPlace ychk
 
 #if DEBUG_LOG
@@ -626,7 +628,6 @@ let EEOverlapArea (e1: Types.Ellipse) (e2: Types.Ellipse) : (float * float[] * f
             printf "\t j=%d, ychk=%f\n" j ychk.[j]
 #endif
 
-        let nychk = Array.length ychk
         let mutable nintpts = 0
 
         let xint = Array.zeroCreate 4
@@ -717,11 +718,11 @@ let EEOverlapArea (e1: Types.Ellipse) (e2: Types.Ellipse) : (float * float[] * f
                 | 4 -> fourintpts xint yint a1 b1 phi_1 a2 b2 h2_tr k2_tr phi_2 aa bb cc dd ee ff
                 | _ -> -1.0
             if nintpts = 0
-            then Some (area, [||], [||])
+            then Some (float32 area, [||], [||])
             else
-                let xTransform = Array.zeroCreate nintpts
-                let yTransform = Array.zeroCreate nintpts
+                let xTransform : float32[] = Array.zeroCreate nintpts
+                let yTransform : float32[] = Array.zeroCreate nintpts
                 for i in 0 .. (nintpts - 1) do
-                    xTransform.[i] <- cos phi_1 * xint.[i] - sin phi_1 * yint.[i] + h1
-                    yTransform.[i] <- sin phi_1 * xint.[i] + cos phi_1 * yint.[i] + k1
-                Some (area, xTransform, yTransform)
\ No newline at end of file
+                    xTransform.[i] <- float32 <| cos phi_1 * xint.[i] - sin phi_1 * yint.[i] + h1
+                    yTransform.[i] <- float32 <| sin phi_1 * xint.[i] + cos phi_1 * yint.[i] + k1
+                Some (float32 area, xTransform, yTransform)
\ No newline at end of file