-module ParasitemiaCore.EEOver
+// Translation from https://github.com/chraibi/EEOver.
+module ParasitemiaCore.EEOver
open System
-let private EPS = 1.0e-5
+let private EPS = 1.0e-7
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 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
+ for i = 0 to 2 do
if istanpt xint.[i] yint.[i] a1 b2 aa bb cc dd ee ff = TANGENT_POINT
then
tanpts <- tanpts + 1
()
twointpts xint yint a1 b1 phi_1 a2 b2 h2_tr k2_tr phi_2 aa bb cc dd ee ff
-
let private fourintpts (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 a1b1 = a1 * b1
let a2b2 = a2 * b2
let theta = Array.zeroCreate 4
- for i in 0 .. 3 do
+ for i = 0 to 3 do
if abs xint.[i] > a1
then
xint.[i] <- if xint.[i] < 0.0 then -a1 else a1
theta.[i] <- if yint.[i] < 0.0 then 2.0 * Math.PI - acos (xint.[i] / a1) else acos (xint.[i] / a1)
#if DEBUG_LOG
- for k in 0..3 do
+ for k = 0 to 3 do
printf "k=%d: Theta = %f, xint=%f, yint=%f\n" k theta.[k] xint.[k] yint.[k]
#endif
- for j in 1 .. 3 do
+ for j = 1 to 3 do
let tmp0 = theta.[j]
let tmp1 = xint.[j]
let tmp2 = yint.[j]
#if DEBUG_LOG
printf "AFTER sorting\n"
- for k in 0..3 do
+ for k = 0 to 3 do
printf "k=%d: Theta = %f, xint=%f, yint=%f\n" k theta.[k] xint.[k] yint.[k]
#endif
let xint_tr = Array.zeroCreate 4
let yint_tr = Array.zeroCreate 4
- for i in 0..3 do
+ for i = 0 to 3 do
xint_tr.[i] <- (xint.[i] - h2_tr) * cosphi + (yint.[i] - k2_tr) * -sinphi
yint_tr.[i] <- (xint.[i] - h2_tr) * sinphi + (yint.[i] - k2_tr) * cosphi
if area5 < 0.0
then
#if DEBUG_LOG
- printf "\n\t\t-------------> area5 is negativ (%f). Add: pi*A2*B2=%f <------------\n" area5 area_2
+ printf "\n\t\t-------------> area5 is negative (%f). Add: pi*A2*B2=%f <------------\n" area5 area_2
#endif
area5 <- area5 + area_2
if area4 < 0.0
then
#if DEBUG_LOG
- printf "\n\t\t-------------> area4 is negativ (%f). Add: pi*A2*B2=%f <------------\n" area4 area_2
+ printf "\n\t\t-------------> area4 is negative (%f). Add: pi*A2*B2=%f <------------\n" area4 area_2
#endif
area4 <- area4 + area_2
if area3 < 0.0
then
#if DEBUG_LOG
- printf "\n\t\t-------------> area3 is negativ (%f). Add: pi*A2*B2=%f <------------\n" area3 area_1
+ printf "\n\t\t-------------> area3 is negative (%f). Add: pi*A2*B2=%f <------------\n" area3 area_1
#endif
area3 <- area3 + area_1
if area2 < 0.0
then
#if DEBUG_LOG
- printf "\n\t\t-------------> area2 is negativ (%f). Add: pi*A2*B2=%f <------------\n" area2 area_1
+ printf "\n\t\t-------------> area2 is negative (%f). Add: pi*A2*B2=%f <------------\n" area2 area_1
#endif
area2 <- area2 + area_1
let private cubicroots (p: float[]) (r: float[,]) =
if p.[0] <> 1.0 then
- for k in 1..3 do
+ for k = 1 to 3 do
p.[k] <- p.[k] / p.[0]
p.[0] <- 1.0
let s = p.[1] / 3.0
t <- d
r.[1, 1] <- t
- for k in 1..3 do
+ for k = 1 to 3 do
r.[2, k] <- 0.0
let private biquadroots (p: float[]) (r: float[,]) =
if p.[0] <> 1.0
then
- for k in 1..4 do
+ for k = 1 to 4 do
p.[k] <- p.[k] / p.[0]
p.[0] <- 1.0
let e = 0.25 * p.[1]
then
p.[2] <- c / b
quadroots p r
- for k in 1..2 do
- for j in 1..2 do
+ for k = 1 to 2 do
+ for j = 1 to 2 do
r.[j, k+2] <- r.[j, k]
p.[1] <- -p.[1]
p.[2] <- b
quadroots p r
- for k in 1..4 do
+ for k = 1 to 4 do
r.[1,k] <- r.[1,k] - e
quadExecuted <- true
then
p.[1] <- 0.0
else
- for k in 1..4 do
+ for k = 1 to 4 do
r.[1, k] <- -e
r.[2, k] <- 0.0
quadExecuted <- true
quad ()
-// Return a tuple (area, x intersections, y intersections)
+/// <summary>
+/// Return a tuple (area, x intersections, y intersections).
+/// </summary>
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
then
None
else
- 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 phi_1 = phi_1 % Math.PI
+ let phi_2 = phi_2 % Math.PI
let h2_tr, k2_tr, phi_2r =
let cosphi = cos phi_1
let sinphi = sin phi_1
|]
#if DEBUG_LOG
- for i in 0..4 do
+ for i = 0 to 4 do
printf "cy[%d]=%f\n" i cy.[i]
#endif
let nroots =
if abs cy.[4] > EPS
then
- for i in 0 .. 3 do
+ for i = 0 to 3 do
py.[4-i] <- cy.[i] / cy.[4]
py.[0] <- 1.0
#if DEBUG_LOG
- for i in 0..4 do
+ for i to 0 to 4 do
printf "py[%d]=%f\n" i py.[i]
#endif
biquadroots py r
elif abs cy.[3] > EPS
then
- for i in 0..2 do
+ for i = 0 to 2 do
py.[3 - i] <- cy.[i] / cy.[3]
py.[0] <- 1.0
cubicroots py r
elif abs cy.[2] > EPS
then
- for i in 0..1 do
+ for i = 0 to 1 do
py.[2-i] <- cy.[i] / cy.[2]
py.[0] <- 1.0
quadroots py r
let ychk = Array.init nroots (fun _ -> Double.MaxValue)
let mutable nychk = 0
- for i in 1 .. nroots do
+ for i = 1 to nroots do
if abs r.[2, i] < EPS
then
ychk.[nychk] <- r.[1, i] * b1
#if DEBUG_LOG
printf "nychk=%d\n" ychk.Length
- for j in 0 .. ychk.Length - 1 do
+ for j = 0 to ychk.Length - 1 do
printf "\t j=%d, ychk=%f\n" j ychk.[j]
#endif
let x2 = -x1
#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)
+ printf "\n\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
if abs (ellipse2tr x1 ychk.[i] aa bb cc dd ee ff) < EPS
#endif
i <- i + 1
-
if returnValue = -1.0
then
None
| 3 -> threeintpts xint yint a1 b1 phi_1 a2 b2 h2_tr k2_tr phi_2 aa bb cc dd ee ff
| 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 (float32 area, [||], [||])
+
+ if area = -1.0
+ then
+ None
+ elif nintpts = 0
+ then
+ Some (float32 area, [||], [||])
else
let xTransform : float32[] = Array.zeroCreate nintpts
let yTransform : float32[] = Array.zeroCreate nintpts