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         area4 <- 0.5 * (a2b2 * (theta_tr.[2] - theta_tr.[1]) - abs (xint_tr.[1] * yint_tr.[2] - xint_tr.[2] * yint_tr.[1]))
 290
 291         if theta_tr.[3] > theta_tr.[0]
 292         then
 293             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]))
 294         else
 295             area5 <- 0.5 * (a2b2 * (theta_tr.[0] - theta_tr.[3]) - abs (xint_tr.[3] * yint_tr.[0] - xint_tr.[0] * yint_tr.[3]))
 296     else
 297         area2 <- 0.5 * (a1b1 * (theta.[2] - theta.[1]) - abs (xint.[1] * yint.[2] - xint.[2] * yint.[1]))
 298         area3 <- 0.5 * (a1b1 * (theta.[0] - (theta.[3] - 2.0 * Math.PI)) - abs (xint.[3] * yint.[0] - xint.[0] * yint.[3]))
 299         area4 <- 0.5 * (a2b2 * (theta_tr.[1] - theta_tr.[0]) - abs (xint_tr.[0] * yint_tr.[1] - xint_tr.[1] * yint_tr.[0]))
 300         area5 <- 0.5 * (a2b2 * (theta_tr.[3] - theta_tr.[2]) - abs (xint_tr.[2] * yint_tr.[3] - xint_tr.[3] * yint_tr.[2]))
 301
 302     if area5 < 0.0
 303     then
 304 #if DEBUG_LOG
 305         printf "\n\t\t-------------> area5 is negativ (%f). Add: pi*A2*B2=%f <------------\n" area5 area_2
 306 #endif
 307         area5 <- area5 + area_2
 308
 309     if area4 < 0.0
 310     then
 311 #if DEBUG_LOG
 312         printf "\n\t\t-------------> area4 is negativ (%f). Add: pi*A2*B2=%f <------------\n" area4 area_2
 313 #endif
 314         area4 <- area4 + area_2
 315
 316     if area3 < 0.0
 317     then
 318 #if DEBUG_LOG
 319         printf "\n\t\t-------------> area3 is negativ (%f). Add: pi*A2*B2=%f <------------\n" area3 area_1
 320 #endif
 321         area3 <- area3 + area_1
 322
 323     if area2 < 0.0
 324     then
 325 #if DEBUG_LOG
 326         printf "\n\t\t-------------> area2 is negativ (%f). Add: pi*A2*B2=%f <------------\n" area2 area_1
 327 #endif
 328         area2 <- area2 + area_1
 329
 330 #if DEBUG_LOG
 331     printf "\narea1=%f, area2=%f area3=%f, area4=%f, area5=%f\n\n" area1 area2 area3 area4 area5
 332 #endif
 333
 334     area1 + area2 + area3 + area4 + area5
 335
 336
 337 let private quadroots (p: float[]) (r: float[,]) =
 338     let mutable b = -p.[1] / (2.0 * p.[0])
 339     let c = p.[2] / p.[0]
 340     let mutable d = b * b - c
 341
 342     if d >= 0.0
 343     then
 344         if b > 0.0
 345         then
 346             b <- sqrt d + b
 347             r.[1, 2] <- b
 348         else
 349             b <- -sqrt d + b
 350             r.[1, 2] <- b
 351         r.[1, 1] <- c / b
 352         r.[2, 1] <- 0.0
 353         r.[2, 2] <- 0.0
 354     else
 355         d <- sqrt -d
 356         r.[2, 1] <- d
 357         r.[2, 2] <- -d
 358         r.[1, 1] <- b
 359         r.[1, 2] <- b
 360
 361 let private cubicroots (p: float[]) (r: float[,]) =
 362     if p.[0] <> 1.0 then
 363         for k in 1..3 do
 364             p.[k] <- p.[k] / p.[0]
 365         p.[0] <- 1.0
 366     let s = p.[1] / 3.0
 367     let mutable t = s * p.[1]
 368     let mutable b = 0.5 * (s * (t / 1.5 - p.[2]) + p.[3])
 369     t <- (t - p.[2]) / 3.0
 370     let mutable c = t * t * t
 371     let mutable d = b * b - c
 372
 373     if d >= 0.0
 374     then
 375         d <- ((sqrt d) + (abs b)) ** (1.0 / 3.0)
 376         if d <> 0.0
 377         then
 378             if b > 0.0
 379             then b <- -d
 380             else b <- d
 381             c <- t / b
 382         d <- sqrt(0.75) * (b - c)
 383         r.[2, 2] <- d
 384         b <- b + c
 385         c <- -0.5 * b - s
 386         r.[1, 2] <- c
 387         if b > 0.0 && s <= 0.0 || b < 0.0 && s > 0.0
 388         then
 389             r.[1, 1] <- c
 390             r.[2, 1] <- -d
 391             r.[1, 3] <- b - s
 392             r.[2, 3] <- 0.0
 393         else
 394             r.[1, 1] <- b - s
 395             r.[2, 1] <- 0.0
 396             r.[1, 3] <- c
 397             r.[2, 3] <- -d
 398     else
 399         if b = 0.0
 400         then d <- (atan 1.0) / 1.5
 401         else d <- atan ((sqrt -d) / (abs b)) / 3.0
 402
 403         if b < 0.0
 404         then b <- 2.0 * (sqrt t)
 405         else b <- -2.0 * (sqrt t)
 406
 407         c <- (cos d) * b
 408         t <- -(sqrt 0.75) * (sin d) * b - 0.5 * c
 409         d <- -t - c - s
 410         c <- c - s
 411         t <- t - s
 412
 413         if abs c > abs t
 414         then
 415             r.[1, 3] <- c
 416         else
 417             r.[1, 3] <- t
 418             t <- c
 419
 420         if abs d > abs t
 421         then
 422             r.[1, 2] <- d
 423         else
 424             r.[1, 2] <- t
 425             t <- d
 426
 427         r.[1, 1] <- t
 428         for k in 1..3 do
 429             r.[2, k] <- 0.0
 430
 431
 432 let private biquadroots (p: float[]) (r: float[,]) =
 433     if p.[0] <> 1.0
 434     then
 435         for k in 1..4 do
 436             p.[k] <- p.[k] / p.[0]
 437         p.[0] <- 1.0
 438     let e = 0.25 * p.[1]
 439     let b = ref (2.0 * e)
 440     let c = ref (!b ** 2.0)
 441     let mutable d = 0.75 * !c
 442     b := p.[3] + !b *(!c - p.[2])
 443     let mutable a = p.[2] - d
 444     c := p.[4] + e * (e * a - p.[3])
 445     a <- a - d
 446
 447     let quadExecuted = ref false
 448     let quad () =
 449         if not !quadExecuted
 450         then
 451             p.[2] <- !c / !b
 452             quadroots p r
 453             for k in 1..2 do
 454                 for j in 1..2 do
 455                     r.[j, k+2] <- r.[j, k]
 456             p.[1] <- -p.[1]
 457             p.[2] <- !b
 458             quadroots p r
 459             for k in 1..4 do
 460                 r.[1,k] <- r.[1,k] - e
 461             quadExecuted := true
 462
 463     p.[1] <- 0.5 * a
 464     p.[2] <- (p.[1] * p.[1] - !c) * 0.25
 465     p.[3] <- !b * !b / -64.0
 466     if p.[3] < 0.0
 467     then
 468         cubicroots p r
 469         let mutable k = 1
 470         while k < 4 do
 471             if r.[2, k] = 0.0 && r.[1, k] > 0.0
 472             then
 473                 d <- r.[1, k] * 4.0
 474                 a <- a + d
 475                 if a >= 0.0 && !b >= 0.0
 476                 then
 477                     p.[1] <- sqrt d
 478                 elif a <= 0.0 && !b <= 0.0
 479                 then
 480                     p.[1] <- sqrt d
 481                 else
 482                     p.[1] <- -(sqrt d)
 483                 b := 0.5 * (a + !b / p.[1])
 484                 quad ()
 485                 k <- 4
 486             k <- k + 1
 487
 488     if not !quadExecuted && p.[2] < 0.0
 489     then
 490         b := sqrt !c
 491         d <- !b + !b - a
 492         p.[1] <- 0.0
 493         if d > 0.0
 494         then
 495             p.[1] <- sqrt d
 496     elif not !quadExecuted
 497     then
 498         if p.[1] > 0.0
 499         then
 500             b := (sqrt p.[2]) * 2.0 + p.[1]
 501         else
 502             b := -(sqrt p.[2]) * 2.0 + p.[1]
 503
 504         if !b <> 0.0
 505         then
 506             p.[1] <- 0.0
 507         else
 508             for k in 1..4 do
 509                 r.[1, k] <- -e
 510                 r.[2, k] <- 0.0
 511             quadExecuted := true
 512
 513     quad ()
 514
 515
 516 let EEOverlapArea (e1: Types.Ellipse) (e2: Types.Ellipse) : float =
 517     let h1, k1, a1, b1, phi_1 = e1.Cx, e1.Cy, e1.A, e1.B, e1.Alpha
 518     let h2, k2, a2, b2, phi_2 = e2.Cx, e2.Cy, e2.A, e2.B, e2.Alpha
 519
 520     if a1 <= EPS || b1 <= EPS || a2 <= EPS || b2 <= EPS
 521     then
 522         -1.0
 523     else
 524         let phi_1 = phi_1 % Math.PI
 525         let phi_2 = phi_2 % Math.PI
 526         let h2_tr, k2_tr, phi_2r =
 527             let cosphi = cos phi_1
 528             let sinphi = sin phi_1
 529             (h2 - h1) * cosphi + (k2 - k1) * sinphi, (h1 - h2) * sinphi + (k2 - k1) * cosphi, (phi_2 - phi_1) % (2.0 * Math.PI)
 530
 531 #if DEBUG_LOG
 532         printf "H2_TR=%f, K2_TR=%f, PHI_2R=%f\n" h2_tr k2_tr phi_2r
 533 #endif
 534
 535         let cosphi = cos phi_2r
 536         let cosphi2 = cosphi ** 2.0
 537         let sinphi = sin phi_2r
 538         let sinphi2 = sinphi ** 2.0
 539         let cosphisinphi = 2.0 * cosphi * sinphi
 540         let a22 = a2 ** 2.0
 541         let b22 = b2 ** 2.0
 542         let tmp0 = (cosphi * h2_tr + sinphi * k2_tr) / a22
 543         let tmp1 = (sinphi * h2_tr - cosphi * k2_tr) / b22
 544         let tmp2 = cosphi * h2_tr + sinphi * k2_tr
 545         let tmp3 = sinphi * h2_tr - cosphi * k2_tr
 546
 547         let aa = cosphi2 / a22 + sinphi2 / b22
 548         let bb = cosphisinphi / a22 - cosphisinphi / b22
 549         let cc = sinphi2 / a22 + cosphi2 / b22
 550         let dd = -2.0 * cosphi * tmp0 - 2.0 * sinphi * tmp1
 551         let ee = -2.0 * sinphi * tmp0 + 2.0 * cosphi * tmp1
 552         let ff = tmp2 * tmp2 / a22 + tmp3 * tmp3 / b22 - 1.0
 553
 554         let cy = [|
 555             (a1 * (a1 * aa - dd) + ff) * (a1 * (a1 * aa + dd) + ff)
 556             2.0 * b1 * (a1 * a1 * (aa * ee - bb * dd) + ee * ff)
 557             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)
 558             2.0 * b1 * (b1 * b1 * cc * ee + a1 * a1 * (bb * dd - aa * ee))
 559             a1 * a1 * a1 * a1 * aa * aa + b1 * b1 * (a1 * a1 * (bb * bb - 2.0 * aa * cc) + b1 * b1 * cc * cc)
 560             |]
 561
 562 #if DEBUG_LOG
 563         for i in 0..4 do
 564             printf "cy[%d]=%f\n" i cy.[i]
 565 #endif
 566
 567         let py = Array.zeroCreate<float> 5
 568         let r = Array2D.zeroCreate<float> 3 5
 569
 570         let nroots =
 571             if abs cy.[4] > EPS
 572             then
 573                 for i in 0 .. 3 do
 574                     py.[4-i] <- cy.[i] / cy.[4]
 575                 py.[0] <- 1.0
 576 #if DEBUG_LOG
 577                 for i in 0..4 do
 578                     printf "py[%d]=%f\n" i py.[i]
 579 #endif
 580                 biquadroots py r
 581                 4
 582
 583             elif abs cy.[3] > EPS
 584             then
 585                 for i in 0..2 do
 586                     py.[3 - i] <- cy.[i] / cy.[3]
 587                 py.[0] <- 1.0
 588                 cubicroots py r
 589                 3
 590
 591             elif abs cy.[2] > EPS
 592             then
 593                 for i in 0..1 do
 594                     py.[2-i] <- cy.[i] / cy.[2]
 595                 py.[0] <- 1.0
 596                 quadroots py r
 597                 2
 598
 599             elif abs cy.[1] > EPS
 600             then
 601                 r.[1, 1] <- -cy.[0] / cy.[1]
 602                 r.[2, 1] <- 0.0
 603                 1
 604
 605             else
 606                 0
 607
 608 #if DEBUG_LOG
 609         printf "nroots = %d\n" nroots
 610 #endif
 611
 612         let ychk = [|
 613             for i in 1 .. nroots do
 614                 if abs r.[2, i] < EPS
 615                 then
 616                     yield r.[1, i] * b1
 617 #if DEBUG_LOG
 618                     printf "ROOT is Real,  i=%d --> %f (B1=%f)\n" i r.[1, i] b1
 619 #endif
 620         |]
 621         Array.sortInPlace ychk
 622
 623 #if DEBUG_LOG
 624         printf "nychk=%d\n" ychk.Length
 625         for j in 0 .. ychk.Length - 1 do
 626             printf "\t j=%d, ychk=%f\n" j ychk.[j]
 627 #endif
 628
 629         let nychk = Array.length ychk
 630         let mutable nintpts = 0
 631
 632         let xint = Array.zeroCreate 4
 633         let yint = Array.zeroCreate 4
 634
 635         let mutable returnValue = 0.0
 636
 637         let mutable i = 0
 638         while returnValue = 0.0 && i < nychk do
 639 #if DEBUG_LOG
 640             printf "------------->i=%d (nychk=%d)\n" i nychk
 641 #endif
 642
 643             if not (i < nychk - 1 && abs (ychk.[i] - ychk.[i+1]) < EPS / 2.0)
 644             then
 645 #if DEBUG_LOG
 646                 printf "check intersecting points. nintps is %d"  nintpts
 647 #endif
 648
 649                 let x1 = if abs ychk.[i] > b1 then 0.0 else a1 * sqrt (1.0 - (ychk.[i] * ychk.[i]) / (b1 * b1))
 650                 let x2 = -x1
 651
 652 #if DEBUG_LOG
 653                 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)
 654                 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)
 655 #endif
 656
 657                 if abs (ellipse2tr x1 ychk.[i] aa bb cc dd ee ff) < EPS
 658                 then
 659                     nintpts <- nintpts + 1
 660 #if DEBUG_LOG
 661                     printf "first if x1. acc nintps=%d\n" nintpts
 662 #endif
 663                     if nintpts > 4
 664                     then
 665                         returnValue <- -1.0
 666                     else
 667                         xint.[nintpts-1] <- x1
 668                         yint.[nintpts-1] <- ychk.[i]
 669 #if DEBUG_LOG
 670                     printf "nintpts=%d, xint=%f, x2=%f, i=%d, yint=%f\n" nintpts x1 x2 i ychk.[i]
 671 #endif
 672
 673                 if returnValue <> -1.0 && abs (ellipse2tr x2 ychk.[i] aa bb cc dd ee ff) < EPS && abs (x2 - x1) > EPS
 674                 then
 675                     nintpts <- nintpts + 1
 676 #if DEBUG_LOG
 677                     printf "first if x2. nintps=%d, Dx=%f (eps2=%f) \n" nintpts (abs (x2 - x1)) EPS
 678 #endif
 679                     if nintpts > 4
 680                     then
 681                         returnValue <- -1.0
 682                     else
 683                         xint.[nintpts-1] <- x2
 684                         yint.[nintpts-1] <- ychk.[i]
 685
 686 #if DEBUG_LOG
 687                     printf "nintpts=%d, x1=%f, xint=%f, i=%d, yint=%f\n" nintpts x1 x2 i ychk.[i]
 688 #endif
 689
 690 #if DEBUG_LOG
 691             else
 692                 printf "i=%d, multiple roots: %f  <--------> %f. continue\n" i ychk.[i] ychk.[i-1]
 693 #endif
 694             i <- i + 1
 695
 696
 697         if returnValue = -1.0
 698         then
 699             returnValue
 700         else
 701             match nintpts with
 702             | 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
 703             | 2 -> match istanpt xint.[0] yint.[0] a1 b1 aa bb cc dd ee ff with
 704                    | TANGENT_POINT ->
 705 #if DEBUG_LOG
 706                         printf "one point is tangent\n"
 707 #endif
 708                         nointpts a1 b1 a2 b2 h1 k1 h2 k2 phi_1 phi_2 h2_tr k2_tr aa bb cc dd ee ff
 709
 710                    | INTERSECTION_POINT ->
 711 #if DEBUG_LOG
 712                         printf "check twointpts\n"
 713 #endif
 714                         twointpts xint yint a1 b1 phi_1 a2 b2 h2_tr k2_tr phi_2 aa bb cc dd ee ff
 715             | 3 -> threeintpts xint yint a1 b1 phi_1 a2 b2 h2_tr k2_tr phi_2 aa bb cc dd ee ff
 716             | 4 -> fourintpts xint yint a1 b1 phi_1 a2 b2 h2_tr k2_tr phi_2 aa bb cc dd ee ff
 717             | _ -> -1.0