1
// Translation from https://github.com/chraibi/EEOver.
2 module ParasitemiaCore.EEOver
6 let private EPS = 1.0e-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
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) =
14 let area_1 = Math.PI * a1b1
15 let area_2 = Math.PI * a2b2
16 let relsize = a1b1 - a2b2
19 if (h2_tr
* h2_tr
) / (a1
* a1
) + (k2_tr
* k2_tr
) / (b1
* b1
) < 1.0 then
24 elif
relsize < 0.0 then
31 if abs
(h1
- h2
) < EPS && abs
(k1
- k2
) < EPS && abs
(area_1 - area_2) < EPS then
36 type private PointType = TANGENT_POINT | INTERSECTION_POINT
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 =
41 if x < 0.0 then -a1
else a1
47 2.0 * Math.PI - acos
(x / a1)
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)
58 let test1 = ellipse2tr x1 y1 aa bb cc dd ee
ff
59 let test2 = ellipse2tr x2 y2 aa bb cc dd ee
ff
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
71 if test1 * test2 > 0.0 then
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
82 2.0 * Math.PI - acos
(x.[0] / a1)
86 if abs
x.[1] > a1 then
87 x.[1] <- if x.[1] < 0.0 then -a1 else a1
91 2.0 * Math.PI - acos
(x.[1] / a1)
95 if theta1 > theta2 then
100 let xmid = a1 * cos
((theta1 + theta2) / 2.0)
101 let ymid = b1
* sin
((theta1 + theta2) / 2.0)
103 if ellipse2tr xmid ymid aa bb cc dd ee
ff > 0.0 then
108 if theta1 > theta2 then
109 theta1 <- theta1 - 2.0 * Math.PI
111 let trsign = if (theta2 - theta1) > Math.PI then 1.0 else -1.0
113 let mutable area1 = 0.5 * (a1 * b1
* (theta2 - theta1) + trsign * abs
(x.[0] * y
.[1] - x.[1] * y
.[0]))
117 printf
"TWO area1=%f\n" area1
119 area1 <- area1 + a1 * b1
121 let cosphi = cos
(phi_1
- phi_2
)
122 let sinphi = sin
(phi_1
- phi_2
)
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
129 if abs
x1_tr > a2
then
130 x1_tr <- if x1_tr < 0.0 then -a2
else a2
133 theta1 <- 2.0 * Math.PI - acos
(x1_tr / a2)
135 theta1 <- acos
(x1_tr / a2)
137 if abs
x2_tr > a2 then
138 x2_tr <- if x2_tr < 0.0 then -a2 else a2
141 theta2 <- 2.0 * Math.PI - acos
(x2_tr / a2)
143 theta2 <- acos
(x2_tr / a2)
145 if theta1 > theta2 then
150 let xmid = a2 * cos
((theta1 + theta2) / 2.0)
151 let ymid = b2
* sin
((theta1 + theta2) / 2.0)
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
158 if (xmid_rt * xmid_rt) / (a1 * a1) + (ymid_rt * ymid_rt) / (b1
* b1
) > 1.0 then
163 if theta1 > theta2 then
164 theta1 <- theta1 - 2.0 * Math.PI
166 let trsign = if theta2 - theta1 > Math.PI then 1.0 else -1.0
168 let mutable area2 = 0.5 * (a2 * b2
* (theta2 - theta1) + trsign * abs
(x1_tr * y2_tr - x2_tr * y1_tr))
171 printf
"TWO area2=%f\n" area2
173 area2 <- area2 + a2 * b2
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
181 if istanpt
xint.[i
] yint
.[i
] a1 b2 aa bb cc dd ee
ff = TANGENT_POINT then
185 printf
"tanindex=%d\n" tanindex
200 twointpts
xint yint
a1 b1 phi_1
a2 b2 h2_tr k2_tr phi_2 aa bb cc dd ee
ff
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 =
205 let area_1 = Math.PI * a1b1
206 let area_2 = Math.PI * a2b2
208 let theta = Array.zeroCreate
4
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)
217 printf
"k=%d: Theta = %f, xint=%f, yint=%f\n" k
theta.[k
] xint.[k
] yint
.[k
]
225 let mutable k = j
- 1
228 if theta.[k] <= tmp0 then
232 theta.[k+1] <- theta.[k]
233 xint.[k+1] <- xint.[k]
234 yint
.[k+1] <- yint
.[k]
244 printf
"AFTER sorting\n"
246 printf
"k=%d: Theta = %f, xint=%f, yint=%f\n" k theta.[k] xint.[k] yint
.[k]
249 let area1 = 0.5 * abs
((xint.[2] - xint.[0]) * (yint
.[3] - yint
.[1]) - (xint.[3] - xint.[1]) * (yint
.[2] - yint
.[0]))
251 let cosphi = cos
(phi_1
- phi_2
)
252 let sinphi = sin
(phi_1
- phi_2
)
254 let theta_tr = Array.zeroCreate
4
255 let xint_tr = Array.zeroCreate
4
256 let yint_tr = Array.zeroCreate
4
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
262 if abs
xint_tr.[i
] > a2 then
263 xint_tr.[i
] <- if xint_tr.[i
] < 0.0 then -a2 else a2
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)
267 let xmid = a1 * cos
((theta.[0] + theta.[1]) / 2.0)
268 let ymid = b1
* sin
((theta.[0] + theta.[1]) / 2.0)
270 let mutable area2, area3
, area4
, area5
= 0.0, 0.0, 0.0, 0.0
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]))
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]))
280 area5
<- 0.5 * (a2b2 * (theta_tr.[0] - theta_tr.[3]) - abs
(xint_tr.[3] * yint_tr.[0] - xint_tr.[0] * yint_tr.[3]))
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]))
289 printf
"\n\t\t-------------> area5 is negative (%f). Add: pi*A2*B2=%f <------------\n" area5
area_2
291 area5
<- area5
+ area_2
295 printf
"\n\t\t-------------> area4 is negative (%f). Add: pi*A2*B2=%f <------------\n" area4
area_2
297 area4
<- area4
+ area_2
301 printf
"\n\t\t-------------> area3 is negative (%f). Add: pi*A2*B2=%f <------------\n" area3
area_1
303 area3
<- area3
+ area_1
307 printf
"\n\t\t-------------> area2 is negative (%f). Add: pi*A2*B2=%f <------------\n" area2 area_1
309 area2 <- area2 + area_1
312 printf
"\narea1=%f, area2=%f area3=%f, area4=%f, area5=%f\n\n" area1 area2 area3 area4 area5
315 area1 + area2 + area3
+ area4
+ area5
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
339 let private cubicroots
(p
: float[]) (r
: float[,]) =
342 p
.[k] <- p
.[k] / p
.[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
352 d <- (sqrt
d + abs
b) ** (1.0 / 3.0)
359 d <- sqrt
(0.75) * (b - c)
364 if b > 0.0 && s <= 0.0 || b < 0.0 && s > 0.0 then
376 d <- (atan
1.0) / 1.5
378 d <- atan
((sqrt
-d) / (abs
b)) / 3.0
386 t <- -(sqrt
0.75
) * (sin
d) * b - 0.5 * c
391 if abs
c > abs
t then
397 if abs
d > abs
t then
407 let inline private biquadroots (p
: float[]) (r
: float[,]) =
410 p
.[k] <- p
.[k] / p
.[0]
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])
421 let mutable quadExecuted = false
423 if not
quadExecuted then
428 r
.[j
, k+2] <- r
.[j
, k]
433 r
.[1,k] <- r
.[1,k] - e
437 p.[2] <- (p.[1] * p.[1] - c) * 0.25
438 p.[3] <- b * b / -64.0
443 if r
.[2, k] = 0.0 && r
.[1, k] > 0.0 then
446 if a >= 0.0 && b >= 0.0 then
448 elif
a <= 0.0 && b <= 0.0 then
452 b <- 0.5 * (a + b / p.[1])
457 if not
quadExecuted && p.[2] < 0.0 then
463 elif not
quadExecuted then
465 b <- (sqrt
p.[2]) * 2.0 + p.[1]
467 b <- -(sqrt
p.[2]) * 2.0 + p.[1]
480 /// Return a tuple (area, x intersections, y intersections).
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
486 if a1 <= EPS || b1 <= EPS || a2 <= EPS || b2 <= EPS then
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)
497 printf
"H2_TR=%f, K2_TR=%f, PHI_2R=%f\n" h2_tr k2_tr phi_2r
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
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
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
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)
529 printf
"cy[%d]=%f\n" i
cy.[i
]
532 let py = Array.zeroCreate
<float> 5
533 let r = Array2D.zeroCreate
<float> 3 5
536 if abs
cy.[4] > EPS then
538 py.[4-i
] <- cy.[i
] / cy.[4]
542 printf
"py[%d]=%f\n" i
py.[i
]
547 elif abs
cy.[3] > EPS then
549 py.[3 - i
] <- cy.[i
] / cy.[3]
554 elif abs
cy.[2] > EPS then
556 py.[2-i
] <- cy.[i
] / cy.[2]
561 elif abs
cy.[1] > EPS then
562 r.[1, 1] <- -cy.[0] / cy.[1]
570 printf
"nroots = %d\n" nroots
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
580 printf
"ROOT is Real, i=%d --> %f (B1=%f)\n" i
r.[1, i
] b1
582 Array.sortInPlace
ychk
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]
590 let mutable nintpts = 0
592 let xint = Array.zeroCreate
4
593 let yint = Array.zeroCreate
4
595 let mutable returnValue = 0.0
598 while returnValue = 0.0 && i < nychk do
600 printf
"------------->i=%d (nychk=%d)\n" i nychk
603 if not
(i < nychk - 1 && abs
(ychk.[i] - ychk.[i+1]) < EPS / 2.0) then
605 printf
"check intersecting points. nintps is %d" nintpts
608 let x1 = if abs
ychk.[i] > b1 then 0.0 else a1 * sqrt
(1.0 - (ychk.[i] * ychk.[i]) / (b1 * b1))
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)
616 if abs
(ellipse2tr x1 ychk.[i] aa bb cc dd ee ff) < EPS then
617 nintpts <- nintpts + 1
619 printf
"first if x1. acc nintps=%d\n" nintpts
624 xint.[nintpts-1] <- x1
625 yint.[nintpts-1] <- ychk.[i]
627 printf
"nintpts=%d, xint=%f, x2=%f, i=%d, yint=%f\n" nintpts x1 x2 i ychk.[i]
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
633 printf
"first if x2. nintps=%d, Dx=%f (eps2=%f) \n" nintpts (abs
(x2 - x1)) EPS
638 xint.[nintpts-1] <- x2
639 yint.[nintpts-1] <- ychk.[i]
642 printf
"nintpts=%d, x1=%f, xint=%f, i=%d, yint=%f\n" nintpts x1 x2 i ychk.[i]
647 printf
"i=%d, multiple roots: %f <--------> %f. continue\n" i ychk.[i] ychk.[i-1]
651 if returnValue = -1.0 then
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
660 printf
"one point is tangent\n"
662 nointpts
a1 b1 a2 b2 h1 k1
h2 k2 phi_1 phi_2 h2_tr k2_tr
aa bb cc dd ee ff
664 | INTERSECTION_POINT ->
666 printf
"check twointpts\n"
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
675 elif
nintpts = 0 then
676 Some (float32
area, [||], [||])
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)