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
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 =
+ let x =
if abs x > a1
then
if x < 0.0 then -a1 else a1
else x
- let theta =
+ let theta =
if y < 0.0
then 2.0 * Math.PI - acos (x / a1)
else acos (x / a1)
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
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 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 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;
-
+ tanpts <- tanpts + 1
+ tanindex <- i
+#if DEBUG_LOG
+ printf "tanindex=%d\n" tanindex
+#endif
+
if tanpts <> 1
then
-1.0
| 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
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]
+ 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
- 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)
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)
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]
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
+ 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 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
+ then
b <- sqrt d + b
r.[1, 2] <- b
else
let mutable c = t * t * t
let mutable d = b * b - c
- if d >= 0.0
+ if d >= 0.0
then
d <- ((sqrt d) + (abs b)) ** (1.0 / 3.0)
if d <> 0.0
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
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 quad () =
- if not !quadExecuted
+ let mutable quadExecuted = false
+ let quad () =
+ 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
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 <- 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
-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
+ if a1 <= EPS || b1 <= EPS || a2 <= EPS || b2 <= EPS
then
- -1.0
+ None
else
- let phi_1 = phi_1 % Math.PI
- let phi_2 = phi_2 % Math.PI
+ 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
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<float> 5
let r = Array2D.zeroCreate<float> 3 5
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
py.[0] <- 1.0
cubicroots py r
3
-
+
elif abs cy.[2] > EPS
then
for i in 0..1 do
else
0
- let ychk = [|
- for i in 1 .. nroots do
- if abs r.[2, i] < EPS
- then
- yield r.[1, i] * b1
- |]
+#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
- let nychk = Array.length 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 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
+ 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
+ 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
- returnValue
+ None
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
-
+ 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 (area, [||], [||])
+ else
+ let xTransform = Array.zeroCreate nintpts
+ let yTransform = 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
projects / master-thesis.git / blobdiff