9a77f48ab44d95f6fc5a0787d89c8abc06c6030f
[master-thesis.git] / 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
19 then
20 if (h2_tr * h2_tr) / (a1 * a1) + (k2_tr * k2_tr) / (b1 * b1) < 1.0
21 then area_2
22 else 0.0
23
24 elif relsize < 0.0
25 then
26 if ff < 0.0
27 then area_1
28 else 0.0
29
30 else
31 if abs (h1 - h2) < EPS && abs (k1 - k2) < EPS && abs (area_1 - area_2) < EPS
32 then area_1
33 else 0.0
34
35 type private PointType = TANGENT_POINT | INTERSECTION_POINT
36
37 let private istanpt (x: float) (y: float) (a1: float) (b1: float) (aa: float) (bb: float) (cc: float) (dd: float) (ee: float) (ff: float) : PointType =
38 let x =
39 if abs x > a1
40 then
41 if x < 0.0 then -a1 else a1
42 else x
43
44 let theta =
45 if y < 0.0
46 then 2.0 * Math.PI - acos (x / a1)
47 else acos (x / a1)
48
49 let eps_radian = 0.1
50
51 let x1 = a1 * cos (theta + eps_radian)
52 let y1 = b1 * sin (theta + eps_radian)
53 let x2 = a1 * cos (theta - eps_radian)
54 let y2 = b1 * sin (theta - eps_radian)
55
56 let test1 = ellipse2tr x1 y1 aa bb cc dd ee ff
57 let test2 = ellipse2tr x2 y2 aa bb cc dd ee ff
58
59 #if DEBUG_LOG
60 printf "\t\t--- debug istanpt with (x,y)=(%f, %f), A1=%f, B1=%f\n" x y a1 b1
61 printf "theta=%f\n" theta
62 printf "eps_Radian=%f\n" eps_radian
63 printf "(x1, y1)=(%f, %f)\n" x1 y1
64 printf "(x2, y2)=(%f, %f)\n" x2 y2
65 printf "test1=%f\n" test1
66 printf "test2=%f\n" test2
67 #endif
68
69 if test1 * test2 > 0.0
70 then TANGENT_POINT
71 else INTERSECTION_POINT
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 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 =
185 let mutable tanpts = 0
186 let mutable tanindex = 0
187 for i = 0 to 2 do
188 if istanpt xint.[i] yint.[i] a1 b2 aa bb cc dd ee ff = TANGENT_POINT
189 then
190 tanpts <- tanpts + 1
191 tanindex <- i
192 #if DEBUG_LOG
193 printf "tanindex=%d\n" tanindex
194 #endif
195
196 if tanpts <> 1
197 then
198 -1.0
199 else
200 match tanindex with
201 | 0 ->
202 xint.[0] <- xint.[2]
203 yint.[0] <- yint.[2]
204 | 1 ->
205 xint.[1] <- xint.[2]
206 yint.[1] <- yint.[2]
207 | _ ->
208 ()
209 twointpts xint yint a1 b1 phi_1 a2 b2 h2_tr k2_tr phi_2 aa bb cc dd ee ff
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 = 0 to 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 = 0 to 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 = 1 to 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 = 0 to 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 = 0 to 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 negative (%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 negative (%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 negative (%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 negative (%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 = 1 to 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 = 1 to 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 = 1 to 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 = 1 to 2 do
451 for j = 1 to 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 = 1 to 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 = 1 to 4 do
506 r.[1, k] <- -e
507 r.[2, k] <- 0.0
508 quadExecuted <- true
509
510 quad ()
511
512 /// <summary>
513 /// Return a tuple (area, x intersections, y intersections).
514 /// </summary>
515 let EEOverlapArea (e1: Types.Ellipse) (e2: Types.Ellipse) : (float32 * float32[] * float32[]) option =
516 let h1, k1, a1, b1, phi_1 = float e1.Cx, float e1.Cy, float e1.A, float e1.B, float e1.Alpha
517 let h2, k2, a2, b2, phi_2 = float e2.Cx, float e2.Cy, float e2.A, float e2.B, float e2.Alpha
518
519 if a1 <= EPS || b1 <= EPS || a2 <= EPS || b2 <= EPS
520 then
521 None
522 else
523 let phi_1 = phi_1 % Math.PI
524 let phi_2 = phi_2 % Math.PI
525 let h2_tr, k2_tr, phi_2r =
526 let cosphi = cos phi_1
527 let sinphi = sin phi_1
528 (h2 - h1) * cosphi + (k2 - k1) * sinphi, (h1 - h2) * sinphi + (k2 - k1) * cosphi, (phi_2 - phi_1) % (2.0 * Math.PI)
529
530 #if DEBUG_LOG
531 printf "H2_TR=%f, K2_TR=%f, PHI_2R=%f\n" h2_tr k2_tr phi_2r
532 #endif
533
534 let cosphi = cos phi_2r
535 let cosphi2 = cosphi ** 2.0
536 let sinphi = sin phi_2r
537 let sinphi2 = sinphi ** 2.0
538 let cosphisinphi = 2.0 * cosphi * sinphi
539 let a22 = a2 ** 2.0
540 let b22 = b2 ** 2.0
541 let tmp0 = (cosphi * h2_tr + sinphi * k2_tr) / a22
542 let tmp1 = (sinphi * h2_tr - cosphi * k2_tr) / b22
543 let tmp2 = cosphi * h2_tr + sinphi * k2_tr
544 let tmp3 = sinphi * h2_tr - cosphi * k2_tr
545
546 let aa = cosphi2 / a22 + sinphi2 / b22
547 let bb = cosphisinphi / a22 - cosphisinphi / b22
548 let cc = sinphi2 / a22 + cosphi2 / b22
549 let dd = -2.0 * cosphi * tmp0 - 2.0 * sinphi * tmp1
550 let ee = -2.0 * sinphi * tmp0 + 2.0 * cosphi * tmp1
551 let ff = tmp2 * tmp2 / a22 + tmp3 * tmp3 / b22 - 1.0
552
553 let cy = [|
554 (a1 * (a1 * aa - dd) + ff) * (a1 * (a1 * aa + dd) + ff)
555 2.0 * b1 * (a1 * a1 * (aa * ee - bb * dd) + ee * ff)
556 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)
557 2.0 * b1 * (b1 * b1 * cc * ee + a1 * a1 * (bb * dd - aa * ee))
558 a1 * a1 * a1 * a1 * aa * aa + b1 * b1 * (a1 * a1 * (bb * bb - 2.0 * aa * cc) + b1 * b1 * cc * cc)
559 |]
560
561 #if DEBUG_LOG
562 for i = 0 to 4 do
563 printf "cy[%d]=%f\n" i cy.[i]
564 #endif
565
566 let py = Array.zeroCreate<float> 5
567 let r = Array2D.zeroCreate<float> 3 5
568
569 let nroots =
570 if abs cy.[4] > EPS
571 then
572 for i = 0 to 3 do
573 py.[4-i] <- cy.[i] / cy.[4]
574 py.[0] <- 1.0
575 #if DEBUG_LOG
576 for i to 0 to 4 do
577 printf "py[%d]=%f\n" i py.[i]
578 #endif
579 biquadroots py r
580 4
581
582 elif abs cy.[3] > EPS
583 then
584 for i = 0 to 2 do
585 py.[3 - i] <- cy.[i] / cy.[3]
586 py.[0] <- 1.0
587 cubicroots py r
588 3
589
590 elif abs cy.[2] > EPS
591 then
592 for i = 0 to 1 do
593 py.[2-i] <- cy.[i] / cy.[2]
594 py.[0] <- 1.0
595 quadroots py r
596 2
597
598 elif abs cy.[1] > EPS
599 then
600 r.[1, 1] <- -cy.[0] / cy.[1]
601 r.[2, 1] <- 0.0
602 1
603
604 else
605 0
606
607 #if DEBUG_LOG
608 printf "nroots = %d\n" nroots
609 #endif
610
611 let ychk = Array.init nroots (fun _ -> Double.MaxValue)
612 let mutable nychk = 0
613 for i = 1 to nroots do
614 if abs r.[2, i] < EPS
615 then
616 ychk.[nychk] <- r.[1, i] * b1
617 nychk <- nychk + 1
618 #if DEBUG_LOG
619 printf "ROOT is Real, i=%d --> %f (B1=%f)\n" i r.[1, i] b1
620 #endif
621 Array.sortInPlace ychk
622
623 #if DEBUG_LOG
624 printf "nychk=%d\n" ychk.Length
625 for j = 0 to ychk.Length - 1 do
626 printf "\t j=%d, ychk=%f\n" j ychk.[j]
627 #endif
628
629 let mutable nintpts = 0
630
631 let xint = Array.zeroCreate 4
632 let yint = Array.zeroCreate 4
633
634 let mutable returnValue = 0.0
635
636 let mutable i = 0
637 while returnValue = 0.0 && i < nychk do
638 #if DEBUG_LOG
639 printf "------------->i=%d (nychk=%d)\n" i nychk
640 #endif
641
642 if not (i < nychk - 1 && abs (ychk.[i] - ychk.[i+1]) < EPS / 2.0)
643 then
644 #if DEBUG_LOG
645 printf "check intersecting points. nintps is %d" nintpts
646 #endif
647
648 let x1 = if abs ychk.[i] > b1 then 0.0 else a1 * sqrt (1.0 - (ychk.[i] * ychk.[i]) / (b1 * b1))
649 let x2 = -x1
650
651 #if DEBUG_LOG
652 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)
653 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)
654 #endif
655
656 if abs (ellipse2tr x1 ychk.[i] aa bb cc dd ee ff) < EPS
657 then
658 nintpts <- nintpts + 1
659 #if DEBUG_LOG
660 printf "first if x1. acc nintps=%d\n" nintpts
661 #endif
662 if nintpts > 4
663 then
664 returnValue <- -1.0
665 else
666 xint.[nintpts-1] <- x1
667 yint.[nintpts-1] <- ychk.[i]
668 #if DEBUG_LOG
669 printf "nintpts=%d, xint=%f, x2=%f, i=%d, yint=%f\n" nintpts x1 x2 i ychk.[i]
670 #endif
671
672 if returnValue <> -1.0 && abs (ellipse2tr x2 ychk.[i] aa bb cc dd ee ff) < EPS && abs (x2 - x1) > EPS
673 then
674 nintpts <- nintpts + 1
675 #if DEBUG_LOG
676 printf "first if x2. nintps=%d, Dx=%f (eps2=%f) \n" nintpts (abs (x2 - x1)) EPS
677 #endif
678 if nintpts > 4
679 then
680 returnValue <- -1.0
681 else
682 xint.[nintpts-1] <- x2
683 yint.[nintpts-1] <- ychk.[i]
684
685 #if DEBUG_LOG
686 printf "nintpts=%d, x1=%f, xint=%f, i=%d, yint=%f\n" nintpts x1 x2 i ychk.[i]
687 #endif
688
689 #if DEBUG_LOG
690 else
691 printf "i=%d, multiple roots: %f <--------> %f. continue\n" i ychk.[i] ychk.[i-1]
692 #endif
693 i <- i + 1
694
695 if returnValue = -1.0
696 then
697 None
698 else
699 let area =
700 match nintpts with
701 | 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
702 | 2 -> match istanpt xint.[0] yint.[0] a1 b1 aa bb cc dd ee ff with
703 | TANGENT_POINT ->
704 #if DEBUG_LOG
705 printf "one point is tangent\n"
706 #endif
707 nointpts a1 b1 a2 b2 h1 k1 h2 k2 phi_1 phi_2 h2_tr k2_tr aa bb cc dd ee ff
708
709 | INTERSECTION_POINT ->
710 #if DEBUG_LOG
711 printf "check twointpts\n"
712 #endif
713 twointpts xint yint a1 b1 phi_1 a2 b2 h2_tr k2_tr phi_2 aa bb cc dd ee ff
714 | 3 -> threeintpts xint yint a1 b1 phi_1 a2 b2 h2_tr k2_tr phi_2 aa bb cc dd ee ff
715 | 4 -> fourintpts xint yint a1 b1 phi_1 a2 b2 h2_tr k2_tr phi_2 aa bb cc dd ee ff
716 | _ -> -1.0
717
718 if area = -1.0
719 then
720 None
721 elif nintpts = 0
722 then
723 Some (float32 area, [||], [||])
724 else
725 let xTransform : float32[] = Array.zeroCreate nintpts
726 let yTransform : float32[] = Array.zeroCreate nintpts
727 for i in 0 .. (nintpts - 1) do
728 xTransform.[i] <- float32 <| cos phi_1 * xint.[i] - sin phi_1 * yint.[i] + h1
729 yTransform.[i] <- float32 <| sin phi_1 * xint.[i] + cos phi_1 * yint.[i] + k1
730 Some (float32 area, xTransform, yTransform)