X-Git-Url: http://git.euphorik.ch/?a=blobdiff_plain;f=Parasitemia%2FParasitemia%2FEEOver.fs;fp=Parasitemia%2FParasitemia%2FEEOver.fs;h=0000000000000000000000000000000000000000;hb=4bfa3cbdc6145e6944f02e24829ab2ef3a851ac1;hp=2bd13b2f6faeb744e0f0d62bfb83bf985965df11;hpb=48ecdfc43001c444eff6ad442986049384674af2;p=master-thesis.git diff --git a/Parasitemia/Parasitemia/EEOver.fs b/Parasitemia/Parasitemia/EEOver.fs deleted file mode 100644 index 2bd13b2..0000000 --- a/Parasitemia/Parasitemia/EEOver.fs +++ /dev/null @@ -1,724 +0,0 @@ -module EEOver - -open System - -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 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 - - 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 - -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 = - if abs x > a1 - then - if x < 0.0 then -a1 else a1 - else x - - let theta = - if y < 0.0 - then 2.0 * Math.PI - acos (x / a1) - else acos (x / a1) - - let eps_radian = 0.1 - - let x1 = a1 * cos (theta + eps_radian) - let y1 = b1 * sin (theta + eps_radian) - let x2 = a1 * cos (theta - eps_radian) - let y2 = b1 * sin (theta - eps_radian) - - let test1 = ellipse2tr x1 y1 aa bb cc dd ee ff - let test2 = ellipse2tr x2 y2 aa bb cc dd ee ff - -#if DEBUG_LOG - printf "\t\t--- debug istanpt with (x,y)=(%f, %f), A1=%f, B1=%f\n" x y a1 b1 - printf "theta=%f\n" theta - printf "eps_Radian=%f\n" eps_radian - printf "(x1, y1)=(%f, %f)\n" x1 y1 - printf "(x2, y2)=(%f, %f)\n" x2 y2 - printf "test1=%f\n" test1 - printf "test2=%f\n" test2 -#endif - - 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 mutable theta1 = - 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 - - let mutable theta2 = - if y.[1] < 0.0 - then 2.0 * Math.PI - acos (x.[1] / a1) - else acos (x.[1] / a1) - - 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 - let tmp = theta1 - theta1 <- theta2 - theta2 <- tmp - - 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 DEBUG_LOG - printf "TWO area1=%f\n" area1 -#endif - area1 <- area1 + a1 * b1 - - let cosphi = cos (phi_1 - phi_2) - let sinphi = sin (phi_1 - phi_2) - - let mutable x1_tr = (x.[0] - h2_tr) * cosphi + (y.[0] - k2_tr) * -sinphi - let mutable y1_tr = (x.[0] - h2_tr) * sinphi + (y.[0] - k2_tr) * cosphi - 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 - x1_tr <- if x1_tr < 0.0 then -a2 else a2 - - 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 - x2_tr <- if x2_tr < 0.0 then -a2 else a2 - - 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 - let tmp = theta1 - theta1 <- theta2 - theta2 <- tmp - - let xmid = a2 * cos ((theta1 + theta2) / 2.0) - let ymid = b2 * sin ((theta1 + theta2) / 2.0) - - let cosphi = cos (phi_2 - phi_1) - let sinphi = sin (phi_2 - phi_1) - 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 - let tmp = theta1 - theta1 <- theta2 - theta2 <- tmp - - 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 DEBUG_LOG - printf "TWO area2=%f\n" area2 -#endif - area2 <- area2 + a2 * b2 - - 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 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 - tanpts <- tanpts + 1 - tanindex <- i -#if DEBUG_LOG - printf "tanindex=%d\n" tanindex -#endif - - if tanpts <> 1 - then - -1.0 - else - match tanindex with - | 0 -> - xint.[0] <- xint.[2] - yint.[0] <- yint.[2] - | 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 - - -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 area_2 = Math.PI * a2b2 - - let theta = Array.zeroCreate 4 - - for i in 0 .. 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 - printf "k=%d: Theta = %f, xint=%f, yint=%f\n" k theta.[k] xint.[k] yint.[k] -#endif - - for j in 1 .. 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 - k2 <- k + 1 - k <- -1 - else - theta.[k+1] <- theta.[k] - xint.[k+1] <- xint.[k] - yint.[k+1] <- yint.[k] - k <- k - 1 - k2 <- k + 1 - - theta.[k2] <- tmp0 - xint.[k2] <- tmp1 - yint.[k2] <- tmp2 - - -#if DEBUG_LOG - printf "AFTER sorting\n" - 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) - let sinphi = sin (phi_1 - phi_2) - - let theta_tr = Array.zeroCreate 4 - let xint_tr = Array.zeroCreate 4 - let yint_tr = Array.zeroCreate 4 - - for i in 0..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 - 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) - let ymid = b1 * sin ((theta.[0] + theta.[1]) / 2.0) - - 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 - 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 - 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])) - else - area2 <- 0.5 * (a1b1 * (theta.[2] - theta.[1]) - abs (xint.[1] * yint.[2] - xint.[2] * yint.[1])) - 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 -#if DEBUG_LOG - printf "\n\t\t-------------> area5 is negativ (%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 -#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 -#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 -#endif - area2 <- area2 + area_1 - -#if DEBUG_LOG - printf "\narea1=%f, area2=%f area3=%f, area4=%f, area5=%f\n\n" area1 area2 area3 area4 area5 -#endif - - area1 + area2 + area3 + area4 + area5 - -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 - b <- sqrt d + b - r.[1, 2] <- b - else - b <- -sqrt d + b - r.[1, 2] <- b - r.[1, 1] <- c / b - r.[2, 1] <- 0.0 - r.[2, 2] <- 0.0 - else - d <- sqrt -d - r.[2, 1] <- d - r.[2, 2] <- -d - r.[1, 1] <- b - r.[1, 2] <- b - -let private cubicroots (p: float[]) (r: float[,]) = - if p.[0] <> 1.0 then - for k in 1..3 do - p.[k] <- p.[k] / p.[0] - p.[0] <- 1.0 - let s = p.[1] / 3.0 - let mutable t = s * p.[1] - let mutable b = 0.5 * (s * (t / 1.5 - p.[2]) + p.[3]) - t <- (t - p.[2]) / 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 - 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 - r.[1, 1] <- c - r.[2, 1] <- -d - r.[1, 3] <- b - s - r.[2, 3] <- 0.0 - else - r.[1, 1] <- b - s - r.[2, 1] <- 0.0 - 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 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 - d <- -t - c - s - c <- c - s - t <- t - s - - if abs c > abs t - then - r.[1, 3] <- c - else - r.[1, 3] <- t - t <- c - - 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 - r.[2, k] <- 0.0 - -let private biquadroots (p: float[]) (r: float[,]) = - if p.[0] <> 1.0 - then - for k in 1..4 do - p.[k] <- p.[k] / p.[0] - p.[0] <- 1.0 - let e = 0.25 * p.[1] - 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]) - a <- a - d - - let mutable quadExecuted = false - let inline quad () = - if not quadExecuted - then - 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 - quadroots p r - for k in 1..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 - cubicroots p r - let mutable k = 1 - while k < 4 do - 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 - p.[1] <- sqrt d - elif a <= 0.0 && b <= 0.0 - then - p.[1] <- sqrt d - else - p.[1] <- -(sqrt d) - b <- 0.5 * (a + b / p.[1]) - quad () - k <- 4 - k <- k + 1 - - if not quadExecuted && p.[2] < 0.0 - then - b <- sqrt c - d <- b + b - a - p.[1] <- 0.0 - if d > 0.0 - then - p.[1] <- sqrt d - 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 - p.[1] <- 0.0 - else - for k in 1..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 = - 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 - 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 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 - - let cosphi = cos phi_2r - let cosphi2 = cosphi ** 2.0 - let sinphi = sin phi_2r - let sinphi2 = sinphi ** 2.0 - let cosphisinphi = 2.0 * cosphi * sinphi - let a22 = a2 ** 2.0 - let b22 = b2 ** 2.0 - let tmp0 = (cosphi * h2_tr + sinphi * k2_tr) / a22 - let tmp1 = (sinphi * h2_tr - cosphi * k2_tr) / b22 - let tmp2 = cosphi * h2_tr + sinphi * k2_tr - let tmp3 = sinphi * h2_tr - cosphi * k2_tr - - let aa = cosphi2 / a22 + sinphi2 / b22 - let bb = cosphisinphi / a22 - cosphisinphi / b22 - let cc = sinphi2 / a22 + cosphi2 / b22 - let dd = -2.0 * cosphi * tmp0 - 2.0 * sinphi * tmp1 - let ee = -2.0 * sinphi * tmp0 + 2.0 * cosphi * tmp1 - let ff = tmp2 * tmp2 / a22 + tmp3 * tmp3 / b22 - 1.0 - - let cy = [| - (a1 * (a1 * aa - dd) + ff) * (a1 * (a1 * aa + dd) + ff) - 2.0 * b1 * (a1 * a1 * (aa * ee - bb * dd) + ee * ff) - a1 * a1 * ((b1 * b1 * (2.0 * aa * cc - bb * bb) + dd * dd - 2.0 * aa * ff) - 2.0 * a1 * a1 * aa * aa) + b1 * b1 * (2.0 * cc * ff + ee * ee) - 2.0 * b1 * (b1 * b1 * cc * ee + a1 * a1 * (bb * dd - aa * ee)) - a1 * a1 * a1 * a1 * aa * aa + b1 * b1 * (a1 * a1 * (bb * bb - 2.0 * aa * cc) + b1 * b1 * cc * cc) - |] - -#if DEBUG_LOG - for i in 0..4 do - printf "cy[%d]=%f\n" i cy.[i] -#endif - - let py = Array.zeroCreate 5 - let r = Array2D.zeroCreate 3 5 - - let nroots = - if abs cy.[4] > EPS - then - for i in 0 .. 3 do - py.[4-i] <- cy.[i] / cy.[4] - py.[0] <- 1.0 -#if DEBUG_LOG - for i in 0..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 - 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 - py.[2-i] <- cy.[i] / cy.[2] - py.[0] <- 1.0 - quadroots py r - 2 - - elif abs cy.[1] > EPS - then - r.[1, 1] <- -cy.[0] / cy.[1] - r.[2, 1] <- 0.0 - 1 - - else - 0 - -#if DEBUG_LOG - printf "nroots = %d\n" nroots -#endif - - 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 -#endif - Array.sortInPlace ychk - -#if DEBUG_LOG - printf "nychk=%d\n" ychk.Length - for j in 0 .. ychk.Length - 1 do - printf "\t j=%d, ychk=%f\n" j ychk.[j] -#endif - - let mutable nintpts = 0 - - let xint = Array.zeroCreate 4 - let yint = Array.zeroCreate 4 - - let mutable returnValue = 0.0 - - let mutable i = 0 - while returnValue = 0.0 && i < nychk do -#if DEBUG_LOG - 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 DEBUG_LOG - printf "check intersecting points. nintps is %d" nintpts -#endif - - let x1 = if abs ychk.[i] > b1 then 0.0 else a1 * sqrt (1.0 - (ychk.[i] * ychk.[i]) / (b1 * b1)) - 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) -#endif - - 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 - returnValue <- -1.0 - else - xint.[nintpts-1] <- x1 - yint.[nintpts-1] <- ychk.[i] -#if DEBUG_LOG - 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 - 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 - returnValue <- -1.0 - else - xint.[nintpts-1] <- x2 - yint.[nintpts-1] <- ychk.[i] - -#if DEBUG_LOG - printf "nintpts=%d, x1=%f, xint=%f, i=%d, yint=%f\n" nintpts x1 x2 i ychk.[i] -#endif - -#if DEBUG_LOG - else - 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 - else - let area = - 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 -> -#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 -> -#if DEBUG_LOG - printf "check twointpts\n" -#endif - twointpts xint yint a1 b1 phi_1 a2 b2 h2_tr k2_tr phi_2 aa bb cc dd ee ff - | 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, [||], [||]) - else - let xTransform : float32[] = Array.zeroCreate nintpts - let yTransform : float32[] = Array.zeroCreate nintpts - for i in 0 .. (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