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 test1 * test2 > 0.0
  59     then TANGENT_POINT
  60     else INTERSECTION_POINT
  61
  62
  63 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) =
  64     if abs x.[0] > a1
  65     then x.[0] <- if x.[0] < 0.0 then -a1 else a1
  66
  67     let mutable theta1 =
  68         if y.[0] < 0.0
  69         then 2.0 * Math.PI - acos (x.[0] / a1)
  70         else acos (x.[0] / a1)
  71
  72     if abs x.[1] > a1
  73     then x.[1] <- if x.[1] < 0.0 then -a1 else a1
  74
  75     let mutable theta2 =
  76         if y.[1] < 0.0
  77         then 2.0 * Math.PI - acos (x.[1] / a1)
  78         else acos (x.[1] / a1)
  79
  80     if theta1 > theta2
  81     then
  82         let tmp = theta1
  83         theta1 <- theta2
  84         theta2 <- tmp
  85
  86     let xmid = a1 * cos ((theta1 + theta2) / 2.0)
  87     let ymid = b1 * sin ((theta1 + theta2) / 2.0)
  88
  89     if ellipse2tr xmid ymid aa bb cc dd ee ff > 0.0
  90     then
  91         let tmp = theta1
  92         theta1 <- theta2
  93         theta2 <- tmp
  94
  95     if theta1 > theta2
  96     then
  97         theta1 <- theta1 - 2.0 * Math.PI
  98
  99     let trsign = if (theta2 - theta1) > Math.PI then 1.0 else -1.0
 100
 101     let mutable area1 = 0.5 * (a1 * b1 * (theta2 - theta1) + trsign * abs (x.[0] * y.[1] - x.[1] * y.[0]))
 102
 103     if area1 < 0.0
 104     then
 105         area1 <- area1 + a1 * b1
 106
 107     let cosphi = cos (phi_1 - phi_2)
 108     let sinphi = sin (phi_1 - phi_2)
 109
 110     let mutable x1_tr = (x.[0] - h2_tr) * cosphi + (y.[0] - k2_tr) * -sinphi
 111     let mutable y1_tr = (x.[0] - h2_tr) * sinphi + (y.[0] - k2_tr) * cosphi
 112     let mutable x2_tr = (x.[1] - h2_tr) * cosphi + (y.[1] - k2_tr) * -sinphi
 113     let mutable y2_tr = (x.[1] - h2_tr) * sinphi + (y.[1] - k2_tr) * cosphi
 114
 115     if abs x1_tr > a2
 116     then
 117         x1_tr <- if x1_tr < 0.0 then -a2 else a2
 118
 119     if y1_tr < 0.0
 120     then
 121         theta1 <- 2.0 * Math.PI - acos (x1_tr / a2)
 122     else
 123         theta1 <- acos (x1_tr / a2)
 124
 125     if abs x2_tr > a2
 126     then
 127         x2_tr <- if x2_tr < 0.0 then -a2 else a2
 128
 129     if y2_tr < 0.0
 130     then
 131         theta2 <- 2.0 * Math.PI - acos (x2_tr / a2)
 132     else
 133         theta2 <- acos (x2_tr / a2)
 134
 135     if theta1 > theta2
 136     then
 137         let tmp = theta1
 138         theta1 <- theta2
 139         theta2 <- tmp
 140
 141     let xmid = a2 * cos ((theta1 + theta2) / 2.0)
 142     let ymid = b2 * sin ((theta1 + theta2) / 2.0)
 143
 144     let cosphi = cos (phi_2 - phi_1)
 145     let sinphi = sin (phi_2 - phi_1)
 146     let xmid_rt = xmid * cosphi + ymid * -sinphi + h2_tr
 147     let ymid_rt = xmid * sinphi + ymid * cosphi + k2_tr
 148
 149     if (xmid_rt * xmid_rt) / (a1 * a1) + (ymid_rt * ymid_rt) / (b1 * b1) > 1.0
 150     then
 151         let tmp = theta1
 152         theta1 <- theta2
 153         theta2 <- tmp
 154
 155     if theta1 > theta2
 156     then
 157         theta1 <- theta1 - 2.0 * Math.PI
 158
 159     let trsign = if theta2 - theta1 > Math.PI then 1.0 else -1.0
 160
 161     let mutable area2 = 0.5 * (a2 * b2 * (theta2 - theta1) + trsign * abs (x1_tr * y2_tr - x2_tr * y1_tr))
 162     if area2 < 0.0
 163     then
 164         area2 <- area2 + a2 * b2
 165
 166     area1 + area2
 167
 168
 169 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 =
 170     let mutable tanpts = 0
 171     let mutable tanindex = 0
 172     for i in 0..2 do
 173         if istanpt xint.[i] yint.[i] a1 b2 aa bb cc dd ee ff = TANGENT_POINT
 174         then
 175             tanpts <- tanpts + 1;
 176             tanindex <- i;
 177
 178     if tanpts <> 1
 179     then
 180         -1.0
 181     else
 182         match tanindex with
 183         | 0 ->
 184             xint.[0] <- xint.[2]
 185             yint.[0] <- yint.[2]
 186         | 1 ->
 187             xint.[1] <- xint.[2]
 188             yint.[1] <- yint.[2]
 189         | _ ->
 190             ()
 191         twointpts xint yint a1 b1 phi_1 a2 b2 h2_tr k2_tr phi_2 aa bb cc dd ee ff
 192
 193
 194 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 =
 195     let a1b1 = a1 * b1
 196     let a2b2 = a2 * b2
 197     let area_1 = Math.PI * a1b1
 198     let area_2 = Math.PI * a2b2
 199
 200     let theta = Array.zeroCreate 4
 201
 202     for i in 0..3 do
 203         if abs xint.[i] > a1
 204         then
 205             xint.[i] <- if xint.[i] < 0.0 then -a1 else a1
 206         theta.[i] <- if yint.[i] < 0.0 then 2.0 * Math.PI - acos (xint.[i] / a1) else acos (xint.[i] / a1)
 207
 208     for j in 1..3 do
 209         let tmp0 = theta.[j]
 210         let tmp1 = xint.[j]
 211         let tmp2 = yint.[j]
 212
 213         let mutable k = j - 1
 214         while k >= 0 do
 215             if theta.[k] <= tmp0
 216             then
 217                 k <- -1
 218             else
 219                 theta.[k+1] <- theta.[k]
 220                 xint.[k+1] <- xint.[k]
 221                 yint.[k+1] <- yint.[k]
 222                 k <- k - 1
 223
 224         theta.[k+1] <- tmp0
 225         xint.[k+1] <- tmp1
 226         yint.[k+1] <- tmp2
 227
 228     let area1 = 0.5 * abs ((xint.[2] - xint.[0]) * (yint.[3] - yint.[1]) - (xint.[3] - xint.[1]) * (yint.[2] - yint.[0]))
 229
 230     let cosphi = cos (phi_1 - phi_2)
 231     let sinphi = sin (phi_1 - phi_2)
 232
 233     let theta_tr = Array.zeroCreate 4
 234     let xint_tr = Array.zeroCreate 4
 235     let yint_tr = Array.zeroCreate 4
 236
 237     for i in 0..3 do
 238         xint_tr.[i] <- (xint.[i] - h2_tr) * cosphi + (yint.[i] - k2_tr) * -sinphi
 239         yint_tr.[i] <- (xint.[i] - h2_tr) * sinphi + (yint.[i] - k2_tr) * cosphi
 240
 241         if abs xint_tr.[i] > a2
 242         then
 243             xint_tr.[i] <- if xint_tr.[i] < 0.0 then -a2 else a2
 244
 245         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)
 246
 247     let xmid = a1 * cos ((theta.[0] + theta.[1]) / 2.0)
 248     let ymid = b1 * sin ((theta.[0] + theta.[1]) / 2.0)
 249
 250     let mutable area2, area3, area4, area5 = 0.0, 0.0, 0.0, 0.0
 251
 252     if ellipse2tr xmid ymid aa bb cc dd ee ff < 0.0
 253     then
 254         area2 <- 0.5 * (a1b1 * (theta.[1] - theta.[0]) - abs (xint.[0] * yint.[1] - xint.[1] * yint.[0]))
 255         area4 <- 0.5 * (a2b2 * (theta_tr.[2] - theta_tr.[1]) - abs (xint_tr.[1] * yint_tr.[2] - xint_tr.[2] * yint_tr.[1]))
 256
 257         if theta_tr.[3] > theta_tr.[0]
 258         then
 259             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]))
 260         else
 261             area5 <- 0.5 * (a2b2 * (theta_tr.[0] - theta_tr.[3]) - abs (xint_tr.[3] * yint_tr.[0] - xint_tr.[0] * yint_tr.[3]))
 262     else
 263         area2 <- 0.5 * (a1b1 * (theta.[2] - theta.[1]) - abs (xint.[1] * yint.[2] - xint.[2] * yint.[1]))
 264         area3 <- 0.5 * (a1b1 * (theta.[0] - (theta.[3] - 2.0 * Math.PI)) - abs (xint.[3] * yint.[0] - xint.[0] * yint.[3]))
 265         area4 <- 0.5 * (a2b2 * (theta_tr.[1] - theta_tr.[0]) - abs (xint_tr.[0] * yint_tr.[1] - xint_tr.[1] * yint_tr.[0]))
 266         area5 <- 0.5 * (a2b2 * (theta_tr.[3] - theta_tr.[2]) - abs (xint_tr.[2] * yint_tr.[3] - xint_tr.[3] * yint_tr.[2]))
 267
 268     if area5 < 0.0
 269     then
 270         area5 <- area5 + area_2
 271
 272     if area4 < 0.0
 273     then
 274         area4 <- area4 + area_2
 275
 276     if area3 < 0.0
 277     then
 278         area3 <- area3 + area_1
 279
 280     if area2 < 0.0
 281     then
 282         area2 <- area2 + area_1
 283
 284     area1 + area2 + area3 + area4 + area5
 285
 286
 287 let private quadroots (p: float[]) (r: float[,]) =
 288     let mutable b = -p.[1] / (2.0 * p.[0])
 289     let c = p.[2] / p.[0]
 290     let mutable d = b * b - c
 291
 292     if d >= 0.0
 293     then
 294         if b > 0.0
 295         then
 296             b <- sqrt d + b
 297             r.[1, 2] <- b
 298         else
 299             b <- -sqrt d + b
 300             r.[1, 2] <- b
 301         r.[1, 1] <- c / b
 302         r.[2, 1] <- 0.0
 303         r.[2, 2] <- 0.0
 304     else
 305         d <- sqrt -d
 306         r.[2, 1] <- d
 307         r.[2, 2] <- -d
 308         r.[1, 1] <- b
 309         r.[1, 2] <- b
 310
 311 let private cubicroots (p: float[]) (r: float[,]) =
 312     if p.[0] <> 1.0 then
 313         for k in 1..3 do
 314             p.[k] <- p.[k] / p.[0]
 315         p.[0] <- 1.0
 316     let s = p.[1] / 3.0
 317     let mutable t = s * p.[1]
 318     let mutable b = 0.5 * (s * (t / 1.5 - p.[2]) + p.[3])
 319     t <- (t - p.[2]) / 3.0
 320     let mutable c = t * t * t
 321     let mutable d = b * b - c
 322
 323     if d >= 0.0
 324     then
 325         d <- ((sqrt d) + (abs b)) ** (1.0 / 3.0)
 326         if d <> 0.0
 327         then
 328             if b > 0.0
 329             then b <- -d
 330             else b <- d
 331             c <- t / b
 332         d <- sqrt(0.75) * (b - c)
 333         r.[2, 2] <- d
 334         b <- b + c
 335         c <- -0.5 * b - s
 336         r.[1, 2] <- c
 337         if b > 0.0 && s <= 0.0 || b < 0.0 && s > 0.0
 338         then
 339             r.[1, 1] <- c
 340             r.[2, 1] <- -d
 341             r.[1, 3] <- b - s
 342             r.[2, 3] <- 0.0
 343         else
 344             r.[1, 1] <- b - s
 345             r.[2, 1] <- 0.0
 346             r.[1, 3] <- c
 347             r.[2, 3] <- -d
 348     else
 349         if b = 0.0
 350         then d <- (atan 1.0) / 1.5
 351         else d <- atan ((sqrt -d) / (abs b)) / 3.0
 352
 353         if b < 0.0
 354         then b <- 2.0 * (sqrt t)
 355         else b <- -2.0 * (sqrt t)
 356
 357         c <- (cos d) * b
 358         t <- -(sqrt 0.75) * (sin d) * b - 0.5 * c
 359         d <- -t - c - s
 360         c <- c - s
 361         t <- t - s
 362
 363         if abs c > abs t
 364         then
 365             r.[1, 3] <- c
 366         else
 367             r.[1, 3] <- t
 368             t <- c
 369
 370         if abs d > abs t
 371         then
 372             r.[1, 2] <- d
 373         else
 374             r.[1, 2] <- t
 375             t <- d
 376
 377         r.[1, 1] <- t
 378         for k in 1..3 do
 379             r.[2, k] <- 0.0
 380
 381
 382 let private biquadroots (p: float[]) (r: float[,]) =
 383     if p.[0] <> 1.0
 384     then
 385         for k in 1..4 do
 386             p.[k] <- p.[k] / p.[0]
 387         p.[0] <- 1.0
 388     let e = 0.25 * p.[1]
 389     let b = ref (2.0 * e)
 390     let c = ref (!b ** 2.0)
 391     let mutable d = 0.75 * !c
 392     b := p.[3] + !b *(!c - p.[2])
 393     let mutable a = p.[2] - d
 394     c := p.[4] + e * (e * a - p.[3])
 395     a <- a - d
 396
 397     let quadExecuted = ref false
 398     let quad () =
 399         if not !quadExecuted
 400         then
 401             p.[2] <- !c / !b
 402             quadroots p r
 403             for k in 1..2 do
 404                 for j in 1..2 do
 405                     r.[j, k+2] <- r.[j, k]
 406             p.[1] <- -p.[1]
 407             p.[2] <- !b
 408             quadroots p r
 409             for k in 1..4 do
 410                 r.[1,k] <- r.[1,k] - e
 411             quadExecuted := true
 412
 413     p.[1] <- 0.5 * a
 414     p.[2] <- (p.[1] * p.[1] - !c) * 0.25
 415     p.[3] <- !b * !b / -64.0
 416     if p.[3] < 0.0
 417     then
 418         cubicroots p r
 419         let mutable k = 1
 420         while k < 4 do
 421             if r.[2, k] = 0.0 && r.[1, k] > 0.0
 422             then
 423                 d <- r.[1, k] * 4.0
 424                 a <- a + d
 425                 if a >= 0.0 && !b >= 0.0
 426                 then
 427                     p.[1] <- sqrt d
 428                 elif a <= 0.0 && !b <= 0.0
 429                 then
 430                     p.[1] <- sqrt d
 431                 else
 432                     p.[1] <- -(sqrt d)
 433                 b := 0.5 * (a + !b / p.[1])
 434                 quad ()
 435                 k <- 4
 436             k <- k + 1
 437
 438     if not !quadExecuted && p.[2] < 0.0
 439     then
 440         b := sqrt !c
 441         d <- !b + !b - a
 442         p.[1] <- 0.0
 443         if d > 0.0
 444         then
 445             p.[1] <- sqrt d
 446     elif not !quadExecuted
 447     then
 448         if p.[1] > 0.0
 449         then
 450             b := (sqrt p.[2]) * 2.0 + p.[1]
 451         else
 452             b := -(sqrt p.[2]) * 2.0 + p.[1]
 453
 454         if !b <> 0.0
 455         then
 456             p.[1] <- 0.0
 457         else
 458             for k in 1..4 do
 459                 r.[1, k] <- -e
 460                 r.[2, k] <- 0.0
 461             quadExecuted := true
 462
 463     quad ()
 464
 465
 466 open Types
 467
 468 let EEOverlapArea (e1: Ellipse) (e2: Ellipse) : float =
 469     let { cx = h1; cy = k1; a = a1; b = b1; alpha = phi_1 } = e1
 470     let { cx = h2; cy = k2; a = a2; b = b2; alpha = phi_2 } = e2
 471
 472     if a1 <= EPS || b1 <= EPS || a2 <= EPS || b2 <= EPS
 473     then
 474         -1.0
 475     else
 476         let phi_1 = phi_1 % Math.PI
 477         let phi_2 = phi_2 % Math.PI
 478         let h2_tr, k2_tr, phi_2r =
 479             let cosphi = cos phi_1
 480             let sinphi = sin phi_1
 481             (h2 - h1) * cosphi + (k2 - k1) * sinphi, (h1 - h2) * sinphi + (k2 - k1) * cosphi, (phi_2 - phi_1) % (2.0 * Math.PI)
 482
 483         let cosphi = cos phi_2r
 484         let cosphi2 = cosphi ** 2.0
 485         let sinphi = sin phi_2r
 486         let sinphi2 = sinphi ** 2.0
 487         let cosphisinphi = 2.0 * cosphi * sinphi
 488         let a22 = a2 ** 2.0
 489         let b22 = b2 ** 2.0
 490         let tmp0 = (cosphi * h2_tr + sinphi * k2_tr) / a22
 491         let tmp1 = (sinphi * h2_tr - cosphi * k2_tr) / b22
 492         let tmp2 = cosphi * h2_tr + sinphi * k2_tr
 493         let tmp3 = sinphi * h2_tr - cosphi * k2_tr
 494
 495         let aa = cosphi2 / a22 + sinphi2 / b22
 496         let bb = cosphisinphi / a22 - cosphisinphi / b22
 497         let cc = sinphi2 / a22 + cosphi2 / b22
 498         let dd = -2.0 * cosphi * tmp0 - 2.0 * sinphi * tmp1
 499         let ee = -2.0 * sinphi * tmp0 + 2.0 * cosphi * tmp1
 500         let ff = tmp2 * tmp2 / a22 + tmp3 * tmp3 / b22 - 1.0
 501
 502         let cy = [|
 503             (a1 * (a1 * aa - dd) + ff) * (a1 * (a1 * aa + dd) + ff)
 504             2.0 * b1 * (a1 * a1 * (aa * ee - bb * dd) + ee * ff)
 505             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)
 506             2.0 * b1 * (b1 * b1 * cc * ee + a1 * a1 * (bb * dd - aa * ee))
 507             a1 * a1 * a1 * a1 * aa * aa + b1 * b1 * (a1 * a1 * (bb * bb - 2.0 * aa * cc) + b1 * b1 * cc * cc)
 508             |]
 509
 510         let py = Array.zeroCreate<float> 5
 511         let r = Array2D.zeroCreate<float> 3 5
 512
 513         let nroots =
 514             if abs cy.[4] > EPS
 515             then
 516                 for i in 0 .. 3 do
 517                     py.[4-i] <- cy.[i] / cy.[4]
 518                 py.[0] <- 1.0
 519                 biquadroots py r
 520                 4
 521
 522             elif abs cy.[3] > EPS
 523             then
 524                 for i in 0..2 do
 525                     py.[3 - i] <- cy.[i] / cy.[3]
 526                 py.[0] <- 1.0
 527                 cubicroots py r
 528                 3
 529
 530             elif abs cy.[2] > EPS
 531             then
 532                 for i in 0..1 do
 533                     py.[2-i] <- cy.[i] / cy.[2]
 534                 py.[0] <- 1.0
 535                 quadroots py r
 536                 2
 537
 538             elif abs cy.[1] > EPS
 539             then
 540                 r.[1, 1] <- -cy.[0] / cy.[1]
 541                 r.[2, 1] <- 0.0
 542                 1
 543
 544             else
 545                 0
 546
 547         let ychk = [|
 548             for i in 1 .. nroots do
 549                 if abs r.[2, i] < EPS
 550                 then
 551                     yield r.[1, i] * b1
 552         |]
 553         Array.sortInPlace ychk
 554
 555         let nychk = Array.length ychk
 556         let mutable nintpts = 0
 557
 558         let xint = Array.zeroCreate 4
 559         let yint = Array.zeroCreate 4
 560
 561         let mutable returnValue = 0.0
 562
 563         let mutable i = 0
 564         while returnValue = 0.0 && i < nychk do
 565             if not (i < nychk - 1 && abs (ychk.[i] - ychk.[i+1]) < EPS / 2.0)
 566             then
 567                 let x1 = if abs ychk.[i] > b1 then 0.0 else a1 * sqrt (1.0 - (ychk.[i] * ychk.[i]) / (b1 * b1))
 568                 let x2 = -x1
 569
 570                 if abs (ellipse2tr x1 ychk.[i] aa bb cc dd ee ff) < EPS
 571                 then
 572                     nintpts <- nintpts + 1
 573                     if nintpts > 4
 574                     then
 575                         returnValue <- -1.0
 576                     else
 577                         xint.[nintpts-1] <- x1
 578                         yint.[nintpts-1] <- ychk.[i]
 579
 580                 if returnValue <> -1.0 && abs (ellipse2tr x2 ychk.[i] aa bb cc dd ee ff) < EPS && abs (x2 - x1) > EPS
 581                 then
 582                     nintpts <- nintpts + 1
 583                     if nintpts > 4
 584                     then
 585                         returnValue <- -1.0
 586                     else
 587                         xint.[nintpts-1] <- x2
 588                         yint.[nintpts-1] <- ychk.[i]
 589             i <- i + 1
 590
 591
 592         if returnValue = -1.0
 593         then
 594             returnValue
 595         else
 596             match nintpts with
 597             | 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
 598             | 2 -> match istanpt xint.[0] yint.[0] a1 b1 aa bb cc dd ee ff with
 599                    | TANGENT_POINT -> nointpts a1 b1 a2 b2 h1 k1 h2 k2 phi_1 phi_2 h2_tr k2_tr aa bb cc dd ee ff
 600                    | INTERSECTION_POINT -> twointpts xint yint a1 b1 phi_1 a2 b2 h2_tr k2_tr phi_2 aa bb cc dd ee ff
 601             | 3 -> threeintpts xint yint a1 b1 phi_1 a2 b2 h2_tr k2_tr phi_2 aa bb cc dd ee ff
 602             | 4 -> fourintpts xint yint a1 b1 phi_1 a2 b2 h2_tr k2_tr phi_2 aa bb cc dd ee ff
 603             | _ -> -1.0
 604