Parasitemia/ParasitemiaCore/EEOver.fs

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