-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 =
+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
let area_2 = Math.PI * a2b2
let relsize = a1b1 - a2b2
- if relsize > 0.0
- then
- if (h2_tr * h2_tr) / (a1 * a1) + (k2_tr * k2_tr) / (b1 * b1) < 1.0
- then area_2
- else 0.0
+ if relsize > 0.0 then
+ if (h2_tr * h2_tr) / (a1 * a1) + (k2_tr * k2_tr) / (b1 * b1) < 1.0 then
+ area_2
+ else
+ 0.0
- elif relsize < 0.0
- then
- if ff < 0.0
- then area_1
- else 0.0
+ elif relsize < 0.0 then
+ if ff < 0.0 then
+ area_1
+ else
+ 0.0
else
- if abs (h1 - h2) < EPS && abs (k1 - k2) < EPS && abs (area_1 - area_2) < EPS
- then area_1
- else 0.0
+ if abs (h1 - h2) < EPS && abs (k1 - k2) < EPS && abs (area_1 - area_2) < EPS then
+ area_1
+ else
+ 0.0
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 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 =
- if abs x > a1
- then
+ if abs x > a1 then
if x < 0.0 then -a1 else a1
- else x
+ else
+ x
let theta =
- if y < 0.0
- then 2.0 * Math.PI - acos (x / a1)
- else acos (x / a1)
+ if y < 0.0 then
+ 2.0 * Math.PI - acos (x / a1)
+ else
+ acos (x / a1)
let eps_radian = 0.1
printf "test2=%f\n" test2
#endif
- if test1 * test2 > 0.0
- then TANGENT_POINT
- else INTERSECTION_POINT
+ if test1 * test2 > 0.0 then
+ TANGENT_POINT
+ else
+ INTERSECTION_POINT
-let private twointpts (x: float[]) (y: 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) =
- if abs x.[0] > a1
- then x.[0] <- if x.[0] < 0.0 then -a1 else a1
+let private twointpts (x : float[]) (y : 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) =
+ if abs x.[0] > a1 then
+ x.[0] <- if x.[0] < 0.0 then -a1 else a1
let mutable theta1 =
- if y.[0] < 0.0
- then 2.0 * Math.PI - acos (x.[0] / a1)
- else acos (x.[0] / a1)
+ if y.[0] < 0.0 then
+ 2.0 * Math.PI - acos (x.[0] / a1)
+ else
+ acos (x.[0] / a1)
- if abs x.[1] > a1
- then x.[1] <- if x.[1] < 0.0 then -a1 else a1
+ if abs x.[1] > a1 then
+ x.[1] <- if x.[1] < 0.0 then -a1 else a1
let mutable theta2 =
- if y.[1] < 0.0
- then 2.0 * Math.PI - acos (x.[1] / a1)
- else acos (x.[1] / a1)
+ if y.[1] < 0.0 then
+ 2.0 * Math.PI - acos (x.[1] / a1)
+ else
+ acos (x.[1] / a1)
- if theta1 > theta2
- then
+ if theta1 > theta2 then
let tmp = theta1
theta1 <- theta2
theta2 <- tmp
let xmid = a1 * cos ((theta1 + theta2) / 2.0)
let ymid = b1 * sin ((theta1 + theta2) / 2.0)
- if ellipse2tr xmid ymid aa bb cc dd ee ff > 0.0
- then
+ if ellipse2tr xmid ymid aa bb cc dd ee ff > 0.0 then
let tmp = theta1
theta1 <- theta2
theta2 <- tmp
- if theta1 > theta2
- then
+ 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]))
- if area1 < 0.0
- then
+ if area1 < 0.0 then
#if DEBUG_LOG
printf "TWO area1=%f\n" area1
#endif
let mutable x2_tr = (x.[1] - h2_tr) * cosphi + (y.[1] - k2_tr) * -sinphi
let mutable y2_tr = (x.[1] - h2_tr) * sinphi + (y.[1] - k2_tr) * cosphi
- if abs x1_tr > a2
- then
+ if abs x1_tr > a2 then
x1_tr <- if x1_tr < 0.0 then -a2 else a2
- if y1_tr < 0.0
- then
+ if y1_tr < 0.0 then
theta1 <- 2.0 * Math.PI - acos (x1_tr / a2)
else
theta1 <- acos (x1_tr / a2)
- if abs x2_tr > a2
- then
+ if abs x2_tr > a2 then
x2_tr <- if x2_tr < 0.0 then -a2 else a2
- if y2_tr < 0.0
- then
+ if y2_tr < 0.0 then
theta2 <- 2.0 * Math.PI - acos (x2_tr / a2)
else
theta2 <- acos (x2_tr / a2)
- if theta1 > theta2
- then
+ if theta1 > theta2 then
let tmp = theta1
theta1 <- theta2
theta2 <- tmp
let xmid_rt = xmid * cosphi + ymid * -sinphi + h2_tr
let ymid_rt = xmid * sinphi + ymid * cosphi + k2_tr
- if (xmid_rt * xmid_rt) / (a1 * a1) + (ymid_rt * ymid_rt) / (b1 * b1) > 1.0
- then
+ if (xmid_rt * xmid_rt) / (a1 * a1) + (ymid_rt * ymid_rt) / (b1 * b1) > 1.0 then
let tmp = theta1
theta1 <- theta2
theta2 <- tmp
- if theta1 > theta2
- then
+ 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 area2 = 0.5 * (a2 * b2 * (theta2 - theta1) + trsign * abs (x1_tr * y2_tr - x2_tr * y1_tr))
- if area2 < 0.0
- then
+ if area2 < 0.0 then
#if DEBUG_LOG
printf "TWO area2=%f\n" area2
#endif
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
- if istanpt xint.[i] yint.[i] a1 b2 aa bb cc dd ee ff = TANGENT_POINT
- then
+ 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
tanindex <- i
#if DEBUG_LOG
printf "tanindex=%d\n" tanindex
#endif
- if tanpts <> 1
- then
+ if tanpts <> 1 then
-1.0
else
match tanindex with
()
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 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 area_1 = Math.PI * a1b1
let theta = Array.zeroCreate 4
- for i in 0 .. 3 do
- if abs xint.[i] > a1
- then
+ 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]
let mutable k = j - 1
let mutable k2 = 0
while k >= 0 do
- if theta.[k] <= tmp0
- then
+ if theta.[k] <= tmp0 then
k2 <- k + 1
k <- -1
else
#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 abs xint_tr.[i] > a2
- then
+ 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 mutable area2, area3, area4, area5 = 0.0, 0.0, 0.0, 0.0
- if ellipse2tr xmid ymid aa bb cc dd ee ff < 0.0
- then
+ 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]
- then
+ if theta_tr.[3] > theta_tr.[0] then
area5 <- 0.5 * (a2b2 * (theta_tr.[0] - (theta_tr.[3] - 2.0 * Math.PI)) - abs (xint_tr.[3] * yint_tr.[0] - xint_tr.[0] * yint_tr.[3]))
else
area5 <- 0.5 * (a2b2 * (theta_tr.[0] - theta_tr.[3]) - abs (xint_tr.[3] * yint_tr.[0] - xint_tr.[0] * yint_tr.[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
+ 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 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 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 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
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
- if d >= 0.0
- then
- if b > 0.0
- then
+ if d >= 0.0 then
+ if b > 0.0 then
b <- sqrt d + b
r.[1, 2] <- b
else
r.[1, 1] <- b
r.[1, 2] <- b
-let private cubicroots (p: float[]) (r: float[,]) =
+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
let mutable c = t * t * t
let mutable d = b * b - c
- if d >= 0.0
- then
- d <- ((sqrt d) + (abs b)) ** (1.0 / 3.0)
- if d <> 0.0
- then
- if b > 0.0
- then b <- -d
- else b <- d
+ if d >= 0.0 then
+ d <- (sqrt d + abs b) ** (1.0 / 3.0)
+ if d <> 0.0 then
+ if b > 0.0 then
+ b <- -d
+ else
+ b <- d
c <- t / b
d <- sqrt(0.75) * (b - c)
r.[2, 2] <- d
b <- b + c
c <- -0.5 * b - s
r.[1, 2] <- c
- if b > 0.0 && s <= 0.0 || b < 0.0 && s > 0.0
- then
+ if b > 0.0 && s <= 0.0 || b < 0.0 && s > 0.0 then
r.[1, 1] <- c
r.[2, 1] <- -d
r.[1, 3] <- b - s
r.[1, 3] <- c
r.[2, 3] <- -d
else
- if b = 0.0
- then d <- (atan 1.0) / 1.5
- else d <- atan ((sqrt -d) / (abs b)) / 3.0
+ if b = 0.0 then
+ d <- (atan 1.0) / 1.5
+ else
+ d <- atan ((sqrt -d) / (abs b)) / 3.0
- if b < 0.0
- then b <- 2.0 * (sqrt t)
- else b <- -2.0 * (sqrt t)
+ if b < 0.0 then
+ b <- 2.0 * (sqrt t)
+ else
+ b <- -2.0 * (sqrt t)
c <- (cos d) * b
t <- -(sqrt 0.75) * (sin d) * b - 0.5 * c
c <- c - s
t <- t - s
- if abs c > abs t
- then
+ if abs c > abs t then
r.[1, 3] <- c
else
r.[1, 3] <- t
t <- c
- if abs d > abs t
- then
+ if abs d > abs t then
r.[1, 2] <- d
else
r.[1, 2] <- t
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
+let inline private biquadroots (p : float[]) (r : float[,]) =
+ if p.[0] <> 1.0 then
+ for k = 1 to 4 do
p.[k] <- p.[k] / p.[0]
p.[0] <- 1.0
let e = 0.25 * p.[1]
let mutable quadExecuted = false
let inline quad () =
- if not quadExecuted
- then
+ if not quadExecuted 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
p.[1] <- 0.5 * a
p.[2] <- (p.[1] * p.[1] - c) * 0.25
p.[3] <- b * b / -64.0
- if p.[3] < 0.0
- then
+ if p.[3] < 0.0 then
cubicroots p r
let mutable k = 1
while k < 4 do
- if r.[2, k] = 0.0 && r.[1, k] > 0.0
- then
+ if r.[2, k] = 0.0 && r.[1, k] > 0.0 then
d <- r.[1, k] * 4.0
a <- a + d
- if a >= 0.0 && b >= 0.0
- then
+ if a >= 0.0 && b >= 0.0 then
p.[1] <- sqrt d
- elif a <= 0.0 && b <= 0.0
- then
+ elif a <= 0.0 && b <= 0.0 then
p.[1] <- sqrt d
else
p.[1] <- -(sqrt d)
k <- 4
k <- k + 1
- if not quadExecuted && p.[2] < 0.0
- then
+ if not quadExecuted && p.[2] < 0.0 then
b <- sqrt c
d <- b + b - a
p.[1] <- 0.0
- if d > 0.0
- then
+ if d > 0.0 then
p.[1] <- sqrt d
- elif not quadExecuted
- then
- if p.[1] > 0.0
- then
+ elif not quadExecuted then
+ if p.[1] > 0.0 then
b <- (sqrt p.[2]) * 2.0 + p.[1]
else
b <- -(sqrt p.[2]) * 2.0 + p.[1]
- if b <> 0.0
- then
+ if b <> 0.0 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)
-let EEOverlapArea (e1: Types.Ellipse) (e2: Types.Ellipse) : (float32 * float32[] * float32[]) option =
+/// <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
- if a1 <= EPS || b1 <= EPS || a2 <= EPS || b2 <= EPS
- then
+ if a1 <= EPS || b1 <= EPS || a2 <= EPS || b2 <= EPS 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 r = Array2D.zeroCreate<float> 3 5
let nroots =
- if abs cy.[4] > EPS
- then
- for i in 0 .. 3 do
+ if abs cy.[4] > EPS then
+ 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
4
- elif abs cy.[3] > EPS
- then
- for i in 0..2 do
+ elif abs cy.[3] > EPS then
+ for i = 0 to 2 do
py.[3 - i] <- cy.[i] / cy.[3]
py.[0] <- 1.0
cubicroots py r
3
- elif abs cy.[2] > EPS
- then
- for i in 0..1 do
+ elif abs cy.[2] > EPS then
+ for i = 0 to 1 do
py.[2-i] <- cy.[i] / cy.[2]
py.[0] <- 1.0
quadroots py r
2
- elif abs cy.[1] > EPS
- then
+ elif abs cy.[1] > EPS then
r.[1, 1] <- -cy.[0] / cy.[1]
r.[2, 1] <- 0.0
1
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
+ for i = 1 to nroots do
+ if abs r.[2, i] < EPS then
ychk.[nychk] <- r.[1, i] * b1
nychk <- nychk + 1
#if DEBUG_LOG
#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
printf "------------->i=%d (nychk=%d)\n" i nychk
#endif
- if not (i < nychk - 1 && abs (ychk.[i] - ychk.[i+1]) < EPS / 2.0)
- then
+ if not (i < nychk - 1 && abs (ychk.[i] - ychk.[i+1]) < EPS / 2.0) then
#if DEBUG_LOG
printf "check intersecting points. nintps is %d" nintpts
#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
- then
+ if abs (ellipse2tr x1 ychk.[i] aa bb cc dd ee ff) < EPS then
nintpts <- nintpts + 1
#if DEBUG_LOG
printf "first if x1. acc nintps=%d\n" nintpts
#endif
- if nintpts > 4
- then
+ if nintpts > 4 then
returnValue <- -1.0
else
xint.[nintpts-1] <- x1
printf "nintpts=%d, xint=%f, x2=%f, i=%d, yint=%f\n" nintpts x1 x2 i ychk.[i]
#endif
- if returnValue <> -1.0 && abs (ellipse2tr x2 ychk.[i] aa bb cc dd ee ff) < EPS && abs (x2 - x1) > EPS
- then
+ if returnValue <> -1.0 && abs (ellipse2tr x2 ychk.[i] aa bb cc dd ee ff) < EPS && abs (x2 - x1) > EPS then
nintpts <- nintpts + 1
#if DEBUG_LOG
printf "first if x2. nintps=%d, Dx=%f (eps2=%f) \n" nintpts (abs (x2 - x1)) EPS
#endif
- if nintpts > 4
- then
+ if nintpts > 4 then
returnValue <- -1.0
else
xint.[nintpts-1] <- x2
#endif
i <- i + 1
-
- if returnValue = -1.0
- then
+ if returnValue = -1.0 then
None
else
let area =
| 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
- for i in 0 .. (nintpts - 1) do
+ for i = 0 to (nintpts - 1) do
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
projects / master-thesis.git / blobdiff