Parasitemia/ParasitemiaCore/EEOver.fs

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