--- /dev/null
+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 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
+ 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
+ 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 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)
+
+ for j in 1..3 do
+ let tmp0 = theta.[j]
+ let tmp1 = xint.[j]
+ let tmp2 = yint.[j]
+
+ let mutable k = j - 1
+ while k >= 0 do
+ if theta.[k] <= tmp0
+ then
+ k <- -1
+ else
+ theta.[k+1] <- theta.[k]
+ xint.[k+1] <- xint.[k]
+ yint.[k+1] <- yint.[k]
+ k <- k - 1
+
+ theta.[k+1] <- tmp0
+ xint.[k+1] <- tmp1
+ yint.[k+1] <- tmp2
+
+ 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]))
+ 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
+ area5 <- area5 + area_2
+
+ if area4 < 0.0
+ then
+ area4 <- area4 + area_2
+
+ if area3 < 0.0
+ then
+ area3 <- area3 + area_1
+
+ if area2 < 0.0
+ then
+ area2 <- area2 + area_1
+
+ 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 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 a = p.[2] - d
+ c := p.[4] + e * (e * a - p.[3])
+ a <- a - d
+
+ let quadExecuted = ref false
+ let 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 ()
+
+
+open Types
+
+let EEOverlapArea (e1: Ellipse) (e2: Ellipse) : float =
+ let { cx = h1; cy = k1; a = a1; b = b1; alpha = phi_1 } = e1
+ let { cx = h2; cy = k2; a = a2; b = b2; alpha = phi_2 } = e2
+
+ if a1 <= EPS || b1 <= EPS || a2 <= EPS || b2 <= EPS
+ then
+ -1.0
+ else
+ 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
+ (h2 - h1) * cosphi + (k2 - k1) * sinphi, (h1 - h2) * sinphi + (k2 - k1) * cosphi, (phi_2 - phi_1) % (2.0 * Math.PI)
+
+ 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)
+ |]
+
+ let py = Array.zeroCreate<float> 5
+ let r = Array2D.zeroCreate<float> 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
+ 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
+
+ let ychk = [|
+ for i in 1 .. nroots do
+ if abs r.[2, i] < EPS
+ then
+ yield r.[1, i] * b1
+ |]
+ Array.sortInPlace ychk
+
+ let nychk = Array.length ychk
+ 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 not (i < nychk - 1 && abs (ychk.[i] - ychk.[i+1]) < EPS / 2.0)
+ then
+ 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 abs (ellipse2tr x1 ychk.[i] aa bb cc dd ee ff) < EPS
+ then
+ nintpts <- nintpts + 1
+ if nintpts > 4
+ then
+ returnValue <- -1.0
+ else
+ xint.[nintpts-1] <- x1
+ yint.[nintpts-1] <- ychk.[i]
+
+ 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 nintpts > 4
+ then
+ returnValue <- -1.0
+ else
+ xint.[nintpts-1] <- x2
+ yint.[nintpts-1] <- ychk.[i]
+ i <- i + 1
+
+
+ if returnValue = -1.0
+ then
+ returnValue
+ else
+ 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 -> 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 -> 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
+
projects / master-thesis.git / blobdiff