Parasitemia/Parasitemia/EEOver.fs

   1 module EEOver
   2
   3 open System
   4
   5 let private EPS = 1.0e-5
   6
   7 let inline private ellipse2tr (x: float) (y: float) (aa: float) (bb: float) (cc: float) (dd: float) (ee: float) (ff: float) : float =
   8     aa * x * x + bb * x * y + cc * y * y + dd * x + ee * y + ff
   9
  10 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) =
  11     let a1b1 = a1 * b1
  12     let a2b2 = a2 * b2
  13     let area_1 = Math.PI * a1b1
  14     let area_2 = Math.PI * a2b2
  15     let relsize = a1b1 - a2b2
  16
  17     if relsize > 0.0
  18     then
  19         if (h2_tr * h2_tr) / (a1 * a1) + (k2_tr * k2_tr) / (b1 * b1) < 1.0
  20         then area_2
  21         else 0.0
  22
  23     elif relsize < 0.0
  24     then
  25         if ff < 0.0
  26         then area_1
  27         else 0.0
  28
  29     else
  30         if abs (h1 - h2) < EPS && abs (k1 - k2) < EPS && abs (area_1 - area_2) < EPS
  31         then area_1
  32         else 0.0
  33
  34 type private PointType = TANGENT_POINT | INTERSECTION_POINT
  35
  36 let private istanpt (x: float) (y: float) (a1: float) (b1: float) (aa: float) (bb: float) (cc: float) (dd: float) (ee: float) (ff: float) : PointType =
  37     let x =
  38         if abs x > a1
  39         then
  40             if x < 0.0 then -a1 else a1
  41         else x
  42
  43     let theta =
  44         if y < 0.0
  45         then 2.0 * Math.PI - acos (x / a1)
  46         else acos (x / a1)
  47
  48     let eps_radian = 0.1
  49
  50     let x1 = a1 * cos (theta + eps_radian)
  51     let y1 = b1 * sin (theta + eps_radian)
  52     let x2 = a1 * cos (theta - eps_radian)
  53     let y2 = b1 * sin (theta - eps_radian)
  54
  55     let test1 = ellipse2tr x1 y1 aa bb cc dd ee ff
  56     let test2 = ellipse2tr x2 y2 aa bb cc dd ee ff
  57
  58 #if DEBUG_LOG
  59     printf "\t\t--- debug istanpt with (x,y)=(%f, %f), A1=%f, B1=%f\n" x y a1 b1
  60     printf "theta=%f\n" theta
  61     printf "eps_Radian=%f\n" eps_radian
  62     printf "(x1, y1)=(%f, %f)\n" x1 y1
  63     printf "(x2, y2)=(%f, %f)\n" x2 y2
  64     printf "test1=%f\n" test1
  65     printf "test2=%f\n" test2
  66 #endif
  67
  68     if test1 * test2 > 0.0
  69     then TANGENT_POINT
  70     else INTERSECTION_POINT
  71
  72
  73 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) =
  74     if abs x.[0] > a1
  75     then x.[0] <- if x.[0] < 0.0 then -a1 else a1
  76
  77     let mutable theta1 =
  78         if y.[0] < 0.0
  79         then 2.0 * Math.PI - acos (x.[0] / a1)
  80         else acos (x.[0] / a1)
  81
  82     if abs x.[1] > a1
  83     then x.[1] <- if x.[1] < 0.0 then -a1 else a1
  84
  85     let mutable theta2 =
  86         if y.[1] < 0.0
  87         then 2.0 * Math.PI - acos (x.[1] / a1)
  88         else acos (x.[1] / a1)
  89
  90     if theta1 > theta2
  91     then
  92         let tmp = theta1
  93         theta1 <- theta2
  94         theta2 <- tmp
  95
  96     let xmid = a1 * cos ((theta1 + theta2) / 2.0)
  97     let ymid = b1 * sin ((theta1 + theta2) / 2.0)
  98
  99     if ellipse2tr xmid ymid aa bb cc dd ee ff > 0.0
 100     then
 101         let tmp = theta1
 102         theta1 <- theta2
 103         theta2 <- tmp
 104
 105     if theta1 > theta2
 106     then
 107         theta1 <- theta1 - 2.0 * Math.PI
 108
 109     let trsign = if (theta2 - theta1) > Math.PI then 1.0 else -1.0
 110
 111     let mutable area1 = 0.5 * (a1 * b1 * (theta2 - theta1) + trsign * abs (x.[0] * y.[1] - x.[1] * y.[0]))
 112
 113     if area1 < 0.0
 114     then
 115 #if DEBUG_LOG
 116         printf "TWO area1=%f\n" area1
 117 #endif
 118         area1 <- area1 + a1 * b1
 119
 120     let cosphi = cos (phi_1 - phi_2)
 121     let sinphi = sin (phi_1 - phi_2)
 122
 123     let mutable x1_tr = (x.[0] - h2_tr) * cosphi + (y.[0] - k2_tr) * -sinphi
 124     let mutable y1_tr = (x.[0] - h2_tr) * sinphi + (y.[0] - k2_tr) * cosphi
 125     let mutable x2_tr = (x.[1] - h2_tr) * cosphi + (y.[1] - k2_tr) * -sinphi
 126     let mutable y2_tr = (x.[1] - h2_tr) * sinphi + (y.[1] - k2_tr) * cosphi
 127
 128     if abs x1_tr > a2
 129     then
 130         x1_tr <- if x1_tr < 0.0 then -a2 else a2
 131
 132     if y1_tr < 0.0
 133     then
 134         theta1 <- 2.0 * Math.PI - acos (x1_tr / a2)
 135     else
 136         theta1 <- acos (x1_tr / a2)
 137
 138     if abs x2_tr > a2
 139     then
 140         x2_tr <- if x2_tr < 0.0 then -a2 else a2
 141
 142     if y2_tr < 0.0
 143     then
 144         theta2 <- 2.0 * Math.PI - acos (x2_tr / a2)
 145     else
 146         theta2 <- acos (x2_tr / a2)
 147
 148     if theta1 > theta2
 149     then
 150         let tmp = theta1
 151         theta1 <- theta2
 152         theta2 <- tmp
 153
 154     let xmid = a2 * cos ((theta1 + theta2) / 2.0)
 155     let ymid = b2 * sin ((theta1 + theta2) / 2.0)
 156
 157     let cosphi = cos (phi_2 - phi_1)
 158     let sinphi = sin (phi_2 - phi_1)
 159     let xmid_rt = xmid * cosphi + ymid * -sinphi + h2_tr
 160     let ymid_rt = xmid * sinphi + ymid * cosphi + k2_tr
 161
 162     if (xmid_rt * xmid_rt) / (a1 * a1) + (ymid_rt * ymid_rt) / (b1 * b1) > 1.0
 163     then
 164         let tmp = theta1
 165         theta1 <- theta2
 166         theta2 <- tmp
 167
 168     if theta1 > theta2
 169     then
 170         theta1 <- theta1 - 2.0 * Math.PI
 171
 172     let trsign = if theta2 - theta1 > Math.PI then 1.0 else -1.0
 173
 174     let mutable area2 = 0.5 * (a2 * b2 * (theta2 - theta1) + trsign * abs (x1_tr * y2_tr - x2_tr * y1_tr))
 175     if area2 < 0.0
 176     then
 177 #if DEBUG_LOG
 178         printf "TWO area2=%f\n" area2
 179 #endif
 180         area2 <- area2 + a2 * b2
 181
 182     area1 + area2
 183
 184
 185 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 =
 186     let mutable tanpts = 0
 187     let mutable tanindex = 0
 188     for i in 0..2 do
 189         if istanpt xint.[i] yint.[i] a1 b2 aa bb cc dd ee ff = TANGENT_POINT
 190         then
 191             tanpts <- tanpts + 1
 192             tanindex <- i
 193 #if DEBUG_LOG
 194     printf "tanindex=%d\n" tanindex
 195 #endif
 196
 197     if tanpts <> 1
 198     then
 199         -1.0
 200     else
 201         match tanindex with
 202         | 0 ->
 203             xint.[0] <- xint.[2]
 204             yint.[0] <- yint.[2]
 205         | 1 ->
 206             xint.[1] <- xint.[2]
 207             yint.[1] <- yint.[2]
 208         | _ ->
 209             ()
 210         twointpts xint yint a1 b1 phi_1 a2 b2 h2_tr k2_tr phi_2 aa bb cc dd ee ff
 211
 212
 213 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 =
 214     let a1b1 = a1 * b1
 215     let a2b2 = a2 * b2
 216     let area_1 = Math.PI * a1b1
 217     let area_2 = Math.PI * a2b2
 218
 219     let theta = Array.zeroCreate 4
 220
 221     for i in 0..3 do
 222         if abs xint.[i] > a1
 223         then
 224             xint.[i] <- if xint.[i] < 0.0 then -a1 else a1
 225         theta.[i] <- if yint.[i] < 0.0 then 2.0 * Math.PI - acos (xint.[i] / a1) else acos (xint.[i] / a1)
 226
 227 #if DEBUG_LOG
 228     for k in 0..3 do
 229         printf "k=%d:  Theta = %f, xint=%f, yint=%f\n" k theta.[k] xint.[k] yint.[k]
 230 #endif
 231
 232     for j in 1..3 do
 233         let tmp0 = theta.[j]
 234         let tmp1 = xint.[j]
 235         let tmp2 = yint.[j]
 236
 237         let mutable k = j - 1
 238         let mutable k2 = 0
 239         while k >= 0 do
 240             if theta.[k] <= tmp0
 241             then
 242                 k2 <- k + 1
 243                 k <- -1
 244             else
 245                 theta.[k+1] <- theta.[k]
 246                 xint.[k+1] <- xint.[k]
 247                 yint.[k+1] <- yint.[k]
 248                 k <- k - 1
 249                 k2 <- k + 1
 250
 251         theta.[k2] <- tmp0
 252         xint.[k2] <- tmp1
 253         yint.[k2] <- tmp2
 254
 255
 256 #if DEBUG_LOG
 257     printf "AFTER sorting\n"
 258     for k in 0..3 do
 259         printf "k=%d:  Theta = %f, xint=%f, yint=%f\n" k theta.[k] xint.[k] yint.[k]
 260 #endif
 261
 262     let area1 = 0.5 * abs ((xint.[2] - xint.[0]) * (yint.[3] - yint.[1]) - (xint.[3] - xint.[1]) * (yint.[2] - yint.[0]))
 263
 264     let cosphi = cos (phi_1 - phi_2)
 265     let sinphi = sin (phi_1 - phi_2)
 266
 267     let theta_tr = Array.zeroCreate 4
 268     let xint_tr = Array.zeroCreate 4
 269     let yint_tr = Array.zeroCreate 4
 270
 271     for i in 0..3 do
 272         xint_tr.[i] <- (xint.[i] - h2_tr) * cosphi + (yint.[i] - k2_tr) * -sinphi
 273         yint_tr.[i] <- (xint.[i] - h2_tr) * sinphi + (yint.[i] - k2_tr) * cosphi
 274
 275         if abs xint_tr.[i] > a2
 276         then
 277             xint_tr.[i] <- if xint_tr.[i] < 0.0 then -a2 else a2
 278
 279         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)
 280
 281     let xmid = a1 * cos ((theta.[0] + theta.[1]) / 2.0)
 282     let ymid = b1 * sin ((theta.[0] + theta.[1]) / 2.0)
 283
 284     let mutable area2, area3, area4, area5 = 0.0, 0.0, 0.0, 0.0
 285
 286     if ellipse2tr xmid ymid aa bb cc dd ee ff < 0.0
 287     then
 288         area2 <- 0.5 * (a1b1 * (theta.[1] - theta.[0]) - abs (xint.[0] * yint.[1] - xint.[1] * yint.[0]))
 289         area3 <- 0.5 * (a1b1 * (theta.[3] - theta.[2]) - abs (xint.[2] * yint.[3] - xint.[3] * yint.[2]))
 290         area4 <- 0.5 * (a2b2 * (theta_tr.[2] - theta_tr.[1]) - abs (xint_tr.[1] * yint_tr.[2] - xint_tr.[2] * yint_tr.[1]))
 291
 292         if theta_tr.[3] > theta_tr.[0]
 293         then
 294             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]))
 295         else
 296             area5 <- 0.5 * (a2b2 * (theta_tr.[0] - theta_tr.[3]) - abs (xint_tr.[3] * yint_tr.[0] - xint_tr.[0] * yint_tr.[3]))
 297     else
 298         area2 <- 0.5 * (a1b1 * (theta.[2] - theta.[1]) - abs (xint.[1] * yint.[2] - xint.[2] * yint.[1]))
 299         area3 <- 0.5 * (a1b1 * (theta.[0] - (theta.[3] - 2.0 * Math.PI)) - abs (xint.[3] * yint.[0] - xint.[0] * yint.[3]))
 300         area4 <- 0.5 * (a2b2 * (theta_tr.[1] - theta_tr.[0]) - abs (xint_tr.[0] * yint_tr.[1] - xint_tr.[1] * yint_tr.[0]))
 301         area5 <- 0.5 * (a2b2 * (theta_tr.[3] - theta_tr.[2]) - abs (xint_tr.[2] * yint_tr.[3] - xint_tr.[3] * yint_tr.[2]))
 302
 303     if area5 < 0.0
 304     then
 305 #if DEBUG_LOG
 306         printf "\n\t\t-------------> area5 is negativ (%f). Add: pi*A2*B2=%f <------------\n" area5 area_2
 307 #endif
 308         area5 <- area5 + area_2
 309
 310     if area4 < 0.0
 311     then
 312 #if DEBUG_LOG
 313         printf "\n\t\t-------------> area4 is negativ (%f). Add: pi*A2*B2=%f <------------\n" area4 area_2
 314 #endif
 315         area4 <- area4 + area_2
 316
 317     if area3 < 0.0
 318     then
 319 #if DEBUG_LOG
 320         printf "\n\t\t-------------> area3 is negativ (%f). Add: pi*A2*B2=%f <------------\n" area3 area_1
 321 #endif
 322         area3 <- area3 + area_1
 323
 324     if area2 < 0.0
 325     then
 326 #if DEBUG_LOG
 327         printf "\n\t\t-------------> area2 is negativ (%f). Add: pi*A2*B2=%f <------------\n" area2 area_1
 328 #endif
 329         area2 <- area2 + area_1
 330
 331 #if DEBUG_LOG
 332     printf "\narea1=%f, area2=%f area3=%f, area4=%f, area5=%f\n\n" area1 area2 area3 area4 area5
 333 #endif
 334
 335     area1 + area2 + area3 + area4 + area5
 336
 337
 338 let private quadroots (p: float[]) (r: float[,]) =
 339     let mutable b = -p.[1] / (2.0 * p.[0])
 340     let c = p.[2] / p.[0]
 341     let mutable d = b * b - c
 342
 343     if d >= 0.0
 344     then
 345         if b > 0.0
 346         then
 347             b <- sqrt d + b
 348             r.[1, 2] <- b
 349         else
 350             b <- -sqrt d + b
 351             r.[1, 2] <- b
 352         r.[1, 1] <- c / b
 353         r.[2, 1] <- 0.0
 354         r.[2, 2] <- 0.0
 355     else
 356         d <- sqrt -d
 357         r.[2, 1] <- d
 358         r.[2, 2] <- -d
 359         r.[1, 1] <- b
 360         r.[1, 2] <- b
 361
 362 let private cubicroots (p: float[]) (r: float[,]) =
 363     if p.[0] <> 1.0 then
 364         for k in 1..3 do
 365             p.[k] <- p.[k] / p.[0]
 366         p.[0] <- 1.0
 367     let s = p.[1] / 3.0
 368     let mutable t = s * p.[1]
 369     let mutable b = 0.5 * (s * (t / 1.5 - p.[2]) + p.[3])
 370     t <- (t - p.[2]) / 3.0
 371     let mutable c = t * t * t
 372     let mutable d = b * b - c
 373
 374     if d >= 0.0
 375     then
 376         d <- ((sqrt d) + (abs b)) ** (1.0 / 3.0)
 377         if d <> 0.0
 378         then
 379             if b > 0.0
 380             then b <- -d
 381             else b <- d
 382             c <- t / b
 383         d <- sqrt(0.75) * (b - c)
 384         r.[2, 2] <- d
 385         b <- b + c
 386         c <- -0.5 * b - s
 387         r.[1, 2] <- c
 388         if b > 0.0 && s <= 0.0 || b < 0.0 && s > 0.0
 389         then
 390             r.[1, 1] <- c
 391             r.[2, 1] <- -d
 392             r.[1, 3] <- b - s
 393             r.[2, 3] <- 0.0
 394         else
 395             r.[1, 1] <- b - s
 396             r.[2, 1] <- 0.0
 397             r.[1, 3] <- c
 398             r.[2, 3] <- -d
 399     else
 400         if b = 0.0
 401         then d <- (atan 1.0) / 1.5
 402         else d <- atan ((sqrt -d) / (abs b)) / 3.0
 403
 404         if b < 0.0
 405         then b <- 2.0 * (sqrt t)
 406         else b <- -2.0 * (sqrt t)
 407
 408         c <- (cos d) * b
 409         t <- -(sqrt 0.75) * (sin d) * b - 0.5 * c
 410         d <- -t - c - s
 411         c <- c - s
 412         t <- t - s
 413
 414         if abs c > abs t
 415         then
 416             r.[1, 3] <- c
 417         else
 418             r.[1, 3] <- t
 419             t <- c
 420
 421         if abs d > abs t
 422         then
 423             r.[1, 2] <- d
 424         else
 425             r.[1, 2] <- t
 426             t <- d
 427
 428         r.[1, 1] <- t
 429         for k in 1..3 do
 430             r.[2, k] <- 0.0
 431
 432
 433 let private biquadroots (p: float[]) (r: float[,]) =
 434     if p.[0] <> 1.0
 435     then
 436         for k in 1..4 do
 437             p.[k] <- p.[k] / p.[0]
 438         p.[0] <- 1.0
 439     let e = 0.25 * p.[1]
 440     let mutable b = 2.0 * e
 441     let mutable c = b ** 2.0
 442     let mutable d = 0.75 * c
 443     b <- p.[3] + b *(c - p.[2])
 444     let mutable a = p.[2] - d
 445     c <- p.[4] + e * (e * a - p.[3])
 446     a <- a - d
 447
 448     let mutable quadExecuted = false
 449     let quad () =
 450         if not quadExecuted
 451         then
 452             p.[2] <- c / b
 453             quadroots p r
 454             for k in 1..2 do
 455                 for j in 1..2 do
 456                     r.[j, k+2] <- r.[j, k]
 457             p.[1] <- -p.[1]
 458             p.[2] <- b
 459             quadroots p r
 460             for k in 1..4 do
 461                 r.[1,k] <- r.[1,k] - e
 462             quadExecuted <- true
 463
 464     p.[1] <- 0.5 * a
 465     p.[2] <- (p.[1] * p.[1] - c) * 0.25
 466     p.[3] <- b * b / -64.0
 467     if p.[3] < 0.0
 468     then
 469         cubicroots p r
 470         let mutable k = 1
 471         while k < 4 do
 472             if r.[2, k] = 0.0 && r.[1, k] > 0.0
 473             then
 474                 d <- r.[1, k] * 4.0
 475                 a <- a + d
 476                 if a >= 0.0 && b >= 0.0
 477                 then
 478                     p.[1] <- sqrt d
 479                 elif a <= 0.0 && b <= 0.0
 480                 then
 481                     p.[1] <- sqrt d
 482                 else
 483                     p.[1] <- -(sqrt d)
 484                 b <- 0.5 * (a + b / p.[1])
 485                 quad ()
 486                 k <- 4
 487             k <- k + 1
 488
 489     if not quadExecuted && p.[2] < 0.0
 490     then
 491         b <- sqrt c
 492         d <- b + b - a
 493         p.[1] <- 0.0
 494         if d > 0.0
 495         then
 496             p.[1] <- sqrt d
 497     elif not quadExecuted
 498     then
 499         if p.[1] > 0.0
 500         then
 501             b <- (sqrt p.[2]) * 2.0 + p.[1]
 502         else
 503             b <- -(sqrt p.[2]) * 2.0 + p.[1]
 504
 505         if b <> 0.0
 506         then
 507             p.[1] <- 0.0
 508         else
 509             for k in 1..4 do
 510                 r.[1, k] <- -e
 511                 r.[2, k] <- 0.0
 512             quadExecuted <- true
 513
 514     quad ()
 515
 516 // Return a tuple (area, x intersections, y intersections)
 517 let EEOverlapArea (e1: Types.Ellipse) (e2: Types.Ellipse) : (float * float[] * float[]) option =
 518     let h1, k1, a1, b1, phi_1 = e1.Cx, e1.Cy, e1.A, e1.B, e1.Alpha
 519     let h2, k2, a2, b2, phi_2 = e2.Cx, e2.Cy, e2.A, e2.B, e2.Alpha
 520
 521     if a1 <= EPS || b1 <= EPS || a2 <= EPS || b2 <= EPS
 522     then
 523         None
 524     else
 525         let phi_1 = phi_1 % Math.PI //(if phi_1 > Math.PI / 2.0 then phi_1 - Math.PI else phi_1) % Math.PI
 526         let phi_2 = phi_2 % Math.PI //(if phi_2 > Math.PI / 2.0 then phi_2 - Math.PI else phi_2) % Math.PI
 527         let h2_tr, k2_tr, phi_2r =
 528             let cosphi = cos phi_1
 529             let sinphi = sin phi_1
 530             (h2 - h1) * cosphi + (k2 - k1) * sinphi, (h1 - h2) * sinphi + (k2 - k1) * cosphi, (phi_2 - phi_1) % (2.0 * Math.PI)
 531
 532 #if DEBUG_LOG
 533         printf "H2_TR=%f, K2_TR=%f, PHI_2R=%f\n" h2_tr k2_tr phi_2r
 534 #endif
 535
 536         let cosphi = cos phi_2r
 537         let cosphi2 = cosphi ** 2.0
 538         let sinphi = sin phi_2r
 539         let sinphi2 = sinphi ** 2.0
 540         let cosphisinphi = 2.0 * cosphi * sinphi
 541         let a22 = a2 ** 2.0
 542         let b22 = b2 ** 2.0
 543         let tmp0 = (cosphi * h2_tr + sinphi * k2_tr) / a22
 544         let tmp1 = (sinphi * h2_tr - cosphi * k2_tr) / b22
 545         let tmp2 = cosphi * h2_tr + sinphi * k2_tr
 546         let tmp3 = sinphi * h2_tr - cosphi * k2_tr
 547
 548         let aa = cosphi2 / a22 + sinphi2 / b22
 549         let bb = cosphisinphi / a22 - cosphisinphi / b22
 550         let cc = sinphi2 / a22 + cosphi2 / b22
 551         let dd = -2.0 * cosphi * tmp0 - 2.0 * sinphi * tmp1
 552         let ee = -2.0 * sinphi * tmp0 + 2.0 * cosphi * tmp1
 553         let ff = tmp2 * tmp2 / a22 + tmp3 * tmp3 / b22 - 1.0
 554
 555         let cy = [|
 556             (a1 * (a1 * aa - dd) + ff) * (a1 * (a1 * aa + dd) + ff)
 557             2.0 * b1 * (a1 * a1 * (aa * ee - bb * dd) + ee * ff)
 558             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)
 559             2.0 * b1 * (b1 * b1 * cc * ee + a1 * a1 * (bb * dd - aa * ee))
 560             a1 * a1 * a1 * a1 * aa * aa + b1 * b1 * (a1 * a1 * (bb * bb - 2.0 * aa * cc) + b1 * b1 * cc * cc)
 561             |]
 562
 563 #if DEBUG_LOG
 564         for i in 0..4 do
 565             printf "cy[%d]=%f\n" i cy.[i]
 566 #endif
 567
 568         let py = Array.zeroCreate<float> 5
 569         let r = Array2D.zeroCreate<float> 3 5
 570
 571         let nroots =
 572             if abs cy.[4] > EPS
 573             then
 574                 for i in 0 .. 3 do
 575                     py.[4-i] <- cy.[i] / cy.[4]
 576                 py.[0] <- 1.0
 577 #if DEBUG_LOG
 578                 for i in 0..4 do
 579                     printf "py[%d]=%f\n" i py.[i]
 580 #endif
 581                 biquadroots py r
 582                 4
 583
 584             elif abs cy.[3] > EPS
 585             then
 586                 for i in 0..2 do
 587                     py.[3 - i] <- cy.[i] / cy.[3]
 588                 py.[0] <- 1.0
 589                 cubicroots py r
 590                 3
 591
 592             elif abs cy.[2] > EPS
 593             then
 594                 for i in 0..1 do
 595                     py.[2-i] <- cy.[i] / cy.[2]
 596                 py.[0] <- 1.0
 597                 quadroots py r
 598                 2
 599
 600             elif abs cy.[1] > EPS
 601             then
 602                 r.[1, 1] <- -cy.[0] / cy.[1]
 603                 r.[2, 1] <- 0.0
 604                 1
 605
 606             else
 607                 0
 608
 609 #if DEBUG_LOG
 610         printf "nroots = %d\n" nroots
 611 #endif
 612
 613         let ychk = Array.init nroots (fun _ -> Double.MaxValue)
 614         let mutable nychk = 0
 615         for i in 1 .. nroots do
 616             if abs r.[2, i] < EPS
 617             then
 618                 ychk.[nychk] <- r.[1, i] * b1
 619                 nychk <- nychk + 1
 620 #if DEBUG_LOG
 621                 printf "ROOT is Real,  i=%d --> %f (B1=%f)\n" i r.[1, i] b1
 622 #endif
 623         Array.sortInPlace ychk
 624
 625 #if DEBUG_LOG
 626         printf "nychk=%d\n" ychk.Length
 627         for j in 0 .. ychk.Length - 1 do
 628             printf "\t j=%d, ychk=%f\n" j ychk.[j]
 629 #endif
 630
 631         let mutable nintpts = 0
 632
 633         let xint = Array.zeroCreate 4
 634         let yint = Array.zeroCreate 4
 635
 636         let mutable returnValue = 0.0
 637
 638         let mutable i = 0
 639         while returnValue = 0.0 && i < nychk do
 640 #if DEBUG_LOG
 641             printf "------------->i=%d (nychk=%d)\n" i nychk
 642 #endif
 643
 644             if not (i < nychk - 1 && abs (ychk.[i] - ychk.[i+1]) < EPS / 2.0)
 645             then
 646 #if DEBUG_LOG
 647                 printf "check intersecting points. nintps is %d"  nintpts
 648 #endif
 649
 650                 let x1 = if abs ychk.[i] > b1 then 0.0 else a1 * sqrt (1.0 - (ychk.[i] * ychk.[i]) / (b1 * b1))
 651                 let x2 = -x1
 652
 653 #if DEBUG_LOG
 654                 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)
 655                 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)
 656 #endif
 657
 658                 if abs (ellipse2tr x1 ychk.[i] aa bb cc dd ee ff) < EPS
 659                 then
 660                     nintpts <- nintpts + 1
 661 #if DEBUG_LOG
 662                     printf "first if x1. acc nintps=%d\n" nintpts
 663 #endif
 664                     if nintpts > 4
 665                     then
 666                         returnValue <- -1.0
 667                     else
 668                         xint.[nintpts-1] <- x1
 669                         yint.[nintpts-1] <- ychk.[i]
 670 #if DEBUG_LOG
 671                     printf "nintpts=%d, xint=%f, x2=%f, i=%d, yint=%f\n" nintpts x1 x2 i ychk.[i]
 672 #endif
 673
 674                 if returnValue <> -1.0 && abs (ellipse2tr x2 ychk.[i] aa bb cc dd ee ff) < EPS && abs (x2 - x1) > EPS
 675                 then
 676                     nintpts <- nintpts + 1
 677 #if DEBUG_LOG
 678                     printf "first if x2. nintps=%d, Dx=%f (eps2=%f) \n" nintpts (abs (x2 - x1)) EPS
 679 #endif
 680                     if nintpts > 4
 681                     then
 682                         returnValue <- -1.0
 683                     else
 684                         xint.[nintpts-1] <- x2
 685                         yint.[nintpts-1] <- ychk.[i]
 686
 687 #if DEBUG_LOG
 688                     printf "nintpts=%d, x1=%f, xint=%f, i=%d, yint=%f\n" nintpts x1 x2 i ychk.[i]
 689 #endif
 690
 691 #if DEBUG_LOG
 692             else
 693                 printf "i=%d, multiple roots: %f  <--------> %f. continue\n" i ychk.[i] ychk.[i-1]
 694 #endif
 695             i <- i + 1
 696
 697
 698         if returnValue = -1.0
 699         then
 700             None
 701         else
 702             let area =
 703                 match nintpts with
 704                 | 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
 705                 | 2 -> match istanpt xint.[0] yint.[0] a1 b1 aa bb cc dd ee ff with
 706                        | TANGENT_POINT ->
 707 #if DEBUG_LOG
 708                             printf "one point is tangent\n"
 709 #endif
 710                             nointpts a1 b1 a2 b2 h1 k1 h2 k2 phi_1 phi_2 h2_tr k2_tr aa bb cc dd ee ff
 711
 712                        | INTERSECTION_POINT ->
 713 #if DEBUG_LOG
 714                             printf "check twointpts\n"
 715 #endif
 716                             twointpts xint yint a1 b1 phi_1 a2 b2 h2_tr k2_tr phi_2 aa bb cc dd ee ff
 717                 | 3 -> threeintpts xint yint a1 b1 phi_1 a2 b2 h2_tr k2_tr phi_2 aa bb cc dd ee ff
 718                 | 4 -> fourintpts xint yint a1 b1 phi_1 a2 b2 h2_tr k2_tr phi_2 aa bb cc dd ee ff
 719                 | _ -> -1.0
 720             if nintpts = 0
 721             then Some (area, [||], [||])
 722             else
 723                 let xTransform = Array.zeroCreate nintpts
 724                 let yTransform = Array.zeroCreate nintpts
 725                 for i in 0 .. (nintpts - 1) do
 726                     xTransform.[i] <- cos phi_1 * xint.[i] - sin phi_1 * yint.[i] + h1
 727                     yTransform.[i] <- sin phi_1 * xint.[i] + cos phi_1 * yint.[i] + k1
 728                 Some (area, xTransform, yTransform)