Version 1.0.10
[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 then
19 if (h2_tr * h2_tr) / (a1 * a1) + (k2_tr * k2_tr) / (b1 * b1) < 1.0 then
20 area_2
21 else
22 0.0
23
24 elif relsize < 0.0 then
25 if ff < 0.0 then
26 area_1
27 else
28 0.0
29
30 else
31 if abs (h1 - h2) < EPS && abs (k1 - k2) < EPS && abs (area_1 - area_2) < EPS then
32 area_1
33 else
34 0.0
35
36 type private PointType = TANGENT_POINT | INTERSECTION_POINT
37
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 =
39 let x =
40 if abs x > a1 then
41 if x < 0.0 then -a1 else a1
42 else
43 x
44
45 let theta =
46 if y < 0.0 then
47 2.0 * Math.PI - acos (x / a1)
48 else
49 acos (x / a1)
50
51 let eps_radian = 0.1
52
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)
57
58 let test1 = ellipse2tr x1 y1 aa bb cc dd ee ff
59 let test2 = ellipse2tr x2 y2 aa bb cc dd ee ff
60
61 #if DEBUG_LOG
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
69 #endif
70
71 if test1 * test2 > 0.0 then
72 TANGENT_POINT
73 else
74 INTERSECTION_POINT
75
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
79
80 let mutable theta1 =
81 if y.[0] < 0.0 then
82 2.0 * Math.PI - acos (x.[0] / a1)
83 else
84 acos (x.[0] / a1)
85
86 if abs x.[1] > a1 then
87 x.[1] <- if x.[1] < 0.0 then -a1 else a1
88
89 let mutable theta2 =
90 if y.[1] < 0.0 then
91 2.0 * Math.PI - acos (x.[1] / a1)
92 else
93 acos (x.[1] / a1)
94
95 if theta1 > theta2 then
96 let tmp = theta1
97 theta1 <- theta2
98 theta2 <- tmp
99
100 let xmid = a1 * cos ((theta1 + theta2) / 2.0)
101 let ymid = b1 * sin ((theta1 + theta2) / 2.0)
102
103 if ellipse2tr xmid ymid aa bb cc dd ee ff > 0.0 then
104 let tmp = theta1
105 theta1 <- theta2
106 theta2 <- tmp
107
108 if theta1 > theta2 then
109 theta1 <- theta1 - 2.0 * Math.PI
110
111 let trsign = if (theta2 - theta1) > Math.PI then 1.0 else -1.0
112
113 let mutable area1 = 0.5 * (a1 * b1 * (theta2 - theta1) + trsign * abs (x.[0] * y.[1] - x.[1] * y.[0]))
114
115 if area1 < 0.0 then
116 #if DEBUG_LOG
117 printf "TWO area1=%f\n" area1
118 #endif
119 area1 <- area1 + a1 * b1
120
121 let cosphi = cos (phi_1 - phi_2)
122 let sinphi = sin (phi_1 - phi_2)
123
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
128
129 if abs x1_tr > a2 then
130 x1_tr <- if x1_tr < 0.0 then -a2 else a2
131
132 if y1_tr < 0.0 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 then
138 x2_tr <- if x2_tr < 0.0 then -a2 else a2
139
140 if y2_tr < 0.0 then
141 theta2 <- 2.0 * Math.PI - acos (x2_tr / a2)
142 else
143 theta2 <- acos (x2_tr / a2)
144
145 if theta1 > theta2 then
146 let tmp = theta1
147 theta1 <- theta2
148 theta2 <- tmp
149
150 let xmid = a2 * cos ((theta1 + theta2) / 2.0)
151 let ymid = b2 * sin ((theta1 + theta2) / 2.0)
152
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
157
158 if (xmid_rt * xmid_rt) / (a1 * a1) + (ymid_rt * ymid_rt) / (b1 * b1) > 1.0 then
159 let tmp = theta1
160 theta1 <- theta2
161 theta2 <- tmp
162
163 if theta1 > theta2 then
164 theta1 <- theta1 - 2.0 * Math.PI
165
166 let trsign = if theta2 - theta1 > Math.PI then 1.0 else -1.0
167
168 let mutable area2 = 0.5 * (a2 * b2 * (theta2 - theta1) + trsign * abs (x1_tr * y2_tr - x2_tr * y1_tr))
169 if area2 < 0.0 then
170 #if DEBUG_LOG
171 printf "TWO area2=%f\n" area2
172 #endif
173 area2 <- area2 + a2 * b2
174
175 area1 + area2
176
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
180 for i = 0 to 2 do
181 if istanpt xint.[i] yint.[i] a1 b2 aa bb cc dd ee ff = TANGENT_POINT then
182 tanpts <- tanpts + 1
183 tanindex <- i
184 #if DEBUG_LOG
185 printf "tanindex=%d\n" tanindex
186 #endif
187
188 if tanpts <> 1 then
189 -1.0
190 else
191 match tanindex with
192 | 0 ->
193 xint.[0] <- xint.[2]
194 yint.[0] <- yint.[2]
195 | 1 ->
196 xint.[1] <- xint.[2]
197 yint.[1] <- yint.[2]
198 | _ ->
199 ()
200 twointpts xint yint a1 b1 phi_1 a2 b2 h2_tr k2_tr phi_2 aa bb cc dd ee ff
201
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 =
203 let a1b1 = a1 * b1
204 let a2b2 = a2 * b2
205 let area_1 = Math.PI * a1b1
206 let area_2 = Math.PI * a2b2
207
208 let theta = Array.zeroCreate 4
209
210 for i = 0 to 3 do
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)
214
215 #if DEBUG_LOG
216 for k = 0 to 3 do
217 printf "k=%d: Theta = %f, xint=%f, yint=%f\n" k theta.[k] xint.[k] yint.[k]
218 #endif
219
220 for j = 1 to 3 do
221 let tmp0 = theta.[j]
222 let tmp1 = xint.[j]
223 let tmp2 = yint.[j]
224
225 let mutable k = j - 1
226 let mutable k2 = 0
227 while k >= 0 do
228 if theta.[k] <= tmp0 then
229 k2 <- k + 1
230 k <- -1
231 else
232 theta.[k+1] <- theta.[k]
233 xint.[k+1] <- xint.[k]
234 yint.[k+1] <- yint.[k]
235 k <- k - 1
236 k2 <- k + 1
237
238 theta.[k2] <- tmp0
239 xint.[k2] <- tmp1
240 yint.[k2] <- tmp2
241
242
243 #if DEBUG_LOG
244 printf "AFTER sorting\n"
245 for k = 0 to 3 do
246 printf "k=%d: Theta = %f, xint=%f, yint=%f\n" k theta.[k] xint.[k] yint.[k]
247 #endif
248
249 let area1 = 0.5 * abs ((xint.[2] - xint.[0]) * (yint.[3] - yint.[1]) - (xint.[3] - xint.[1]) * (yint.[2] - yint.[0]))
250
251 let cosphi = cos (phi_1 - phi_2)
252 let sinphi = sin (phi_1 - phi_2)
253
254 let theta_tr = Array.zeroCreate 4
255 let xint_tr = Array.zeroCreate 4
256 let yint_tr = Array.zeroCreate 4
257
258 for i = 0 to 3 do
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
261
262 if abs xint_tr.[i] > a2 then
263 xint_tr.[i] <- if xint_tr.[i] < 0.0 then -a2 else a2
264
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)
266
267 let xmid = a1 * cos ((theta.[0] + theta.[1]) / 2.0)
268 let ymid = b1 * sin ((theta.[0] + theta.[1]) / 2.0)
269
270 let mutable area2, area3, area4, area5 = 0.0, 0.0, 0.0, 0.0
271
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]))
276
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]))
279 else
280 area5 <- 0.5 * (a2b2 * (theta_tr.[0] - theta_tr.[3]) - abs (xint_tr.[3] * yint_tr.[0] - xint_tr.[0] * yint_tr.[3]))
281 else
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]))
286
287 if area5 < 0.0 then
288 #if DEBUG_LOG
289 printf "\n\t\t-------------> area5 is negative (%f). Add: pi*A2*B2=%f <------------\n" area5 area_2
290 #endif
291 area5 <- area5 + area_2
292
293 if area4 < 0.0 then
294 #if DEBUG_LOG
295 printf "\n\t\t-------------> area4 is negative (%f). Add: pi*A2*B2=%f <------------\n" area4 area_2
296 #endif
297 area4 <- area4 + area_2
298
299 if area3 < 0.0 then
300 #if DEBUG_LOG
301 printf "\n\t\t-------------> area3 is negative (%f). Add: pi*A2*B2=%f <------------\n" area3 area_1
302 #endif
303 area3 <- area3 + area_1
304
305 if area2 < 0.0 then
306 #if DEBUG_LOG
307 printf "\n\t\t-------------> area2 is negative (%f). Add: pi*A2*B2=%f <------------\n" area2 area_1
308 #endif
309 area2 <- area2 + area_1
310
311 #if DEBUG_LOG
312 printf "\narea1=%f, area2=%f area3=%f, area4=%f, area5=%f\n\n" area1 area2 area3 area4 area5
313 #endif
314
315 area1 + area2 + area3 + area4 + area5
316
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
321
322 if d >= 0.0 then
323 if b > 0.0 then
324 b <- sqrt d + b
325 r.[1, 2] <- b
326 else
327 b <- -sqrt d + b
328 r.[1, 2] <- b
329 r.[1, 1] <- c / b
330 r.[2, 1] <- 0.0
331 r.[2, 2] <- 0.0
332 else
333 d <- sqrt -d
334 r.[2, 1] <- d
335 r.[2, 2] <- -d
336 r.[1, 1] <- b
337 r.[1, 2] <- b
338
339 let private cubicroots (p : float[]) (r : float[,]) =
340 if p.[0] <> 1.0 then
341 for k = 1 to 3 do
342 p.[k] <- p.[k] / p.[0]
343 p.[0] <- 1.0
344 let s = p.[1] / 3.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
350
351 if d >= 0.0 then
352 d <- (sqrt d + abs b) ** (1.0 / 3.0)
353 if d <> 0.0 then
354 if b > 0.0 then
355 b <- -d
356 else
357 b <- d
358 c <- t / b
359 d <- sqrt(0.75) * (b - c)
360 r.[2, 2] <- d
361 b <- b + c
362 c <- -0.5 * b - s
363 r.[1, 2] <- c
364 if b > 0.0 && s <= 0.0 || b < 0.0 && s > 0.0 then
365 r.[1, 1] <- c
366 r.[2, 1] <- -d
367 r.[1, 3] <- b - s
368 r.[2, 3] <- 0.0
369 else
370 r.[1, 1] <- b - s
371 r.[2, 1] <- 0.0
372 r.[1, 3] <- c
373 r.[2, 3] <- -d
374 else
375 if b = 0.0 then
376 d <- (atan 1.0) / 1.5
377 else
378 d <- atan ((sqrt -d) / (abs b)) / 3.0
379
380 if b < 0.0 then
381 b <- 2.0 * (sqrt t)
382 else
383 b <- -2.0 * (sqrt t)
384
385 c <- (cos d) * b
386 t <- -(sqrt 0.75) * (sin d) * b - 0.5 * c
387 d <- -t - c - s
388 c <- c - s
389 t <- t - s
390
391 if abs c > abs t then
392 r.[1, 3] <- c
393 else
394 r.[1, 3] <- t
395 t <- c
396
397 if abs d > abs t then
398 r.[1, 2] <- d
399 else
400 r.[1, 2] <- t
401 t <- d
402
403 r.[1, 1] <- t
404 for k = 1 to 3 do
405 r.[2, k] <- 0.0
406
407 let inline private biquadroots (p : float[]) (r : float[,]) =
408 if p.[0] <> 1.0 then
409 for k = 1 to 4 do
410 p.[k] <- p.[k] / p.[0]
411 p.[0] <- 1.0
412 let e = 0.25 * p.[1]
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])
419 a <- a - d
420
421 let mutable quadExecuted = false
422 let inline quad () =
423 if not quadExecuted then
424 p.[2] <- c / b
425 quadroots p r
426 for k = 1 to 2 do
427 for j = 1 to 2 do
428 r.[j, k+2] <- r.[j, k]
429 p.[1] <- -p.[1]
430 p.[2] <- b
431 quadroots p r
432 for k = 1 to 4 do
433 r.[1,k] <- r.[1,k] - e
434 quadExecuted <- true
435
436 p.[1] <- 0.5 * a
437 p.[2] <- (p.[1] * p.[1] - c) * 0.25
438 p.[3] <- b * b / -64.0
439 if p.[3] < 0.0 then
440 cubicroots p r
441 let mutable k = 1
442 while k < 4 do
443 if r.[2, k] = 0.0 && r.[1, k] > 0.0 then
444 d <- r.[1, k] * 4.0
445 a <- a + d
446 if a >= 0.0 && b >= 0.0 then
447 p.[1] <- sqrt d
448 elif a <= 0.0 && b <= 0.0 then
449 p.[1] <- sqrt d
450 else
451 p.[1] <- -(sqrt d)
452 b <- 0.5 * (a + b / p.[1])
453 quad ()
454 k <- 4
455 k <- k + 1
456
457 if not quadExecuted && p.[2] < 0.0 then
458 b <- sqrt c
459 d <- b + b - a
460 p.[1] <- 0.0
461 if d > 0.0 then
462 p.[1] <- sqrt d
463 elif not quadExecuted then
464 if p.[1] > 0.0 then
465 b <- (sqrt p.[2]) * 2.0 + p.[1]
466 else
467 b <- -(sqrt p.[2]) * 2.0 + p.[1]
468
469 if b <> 0.0 then
470 p.[1] <- 0.0
471 else
472 for k = 1 to 4 do
473 r.[1, k] <- -e
474 r.[2, k] <- 0.0
475 quadExecuted <- true
476
477 quad ()
478
479 /// <summary>
480 /// Return a tuple (area, x intersections, y intersections).
481 /// </summary>
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
485
486 if a1 <= EPS || b1 <= EPS || a2 <= EPS || b2 <= EPS then
487 None
488 else
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)
495
496 #if DEBUG_LOG
497 printf "H2_TR=%f, K2_TR=%f, PHI_2R=%f\n" h2_tr k2_tr phi_2r
498 #endif
499
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
505 let a22 = a2 ** 2.0
506 let b22 = b2 ** 2.0
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
511
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
518
519 let cy = [|
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)
525 |]
526
527 #if DEBUG_LOG
528 for i = 0 to 4 do
529 printf "cy[%d]=%f\n" i cy.[i]
530 #endif
531
532 let py = Array.zeroCreate<float> 5
533 let r = Array2D.zeroCreate<float> 3 5
534
535 let nroots =
536 if abs cy.[4] > EPS then
537 for i = 0 to 3 do
538 py.[4-i] <- cy.[i] / cy.[4]
539 py.[0] <- 1.0
540 #if DEBUG_LOG
541 for i to 0 to 4 do
542 printf "py[%d]=%f\n" i py.[i]
543 #endif
544 biquadroots py r
545 4
546
547 elif abs cy.[3] > EPS then
548 for i = 0 to 2 do
549 py.[3 - i] <- cy.[i] / cy.[3]
550 py.[0] <- 1.0
551 cubicroots py r
552 3
553
554 elif abs cy.[2] > EPS then
555 for i = 0 to 1 do
556 py.[2-i] <- cy.[i] / cy.[2]
557 py.[0] <- 1.0
558 quadroots py r
559 2
560
561 elif abs cy.[1] > EPS then
562 r.[1, 1] <- -cy.[0] / cy.[1]
563 r.[2, 1] <- 0.0
564 1
565
566 else
567 0
568
569 #if DEBUG_LOG
570 printf "nroots = %d\n" nroots
571 #endif
572
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
578 nychk <- nychk + 1
579 #if DEBUG_LOG
580 printf "ROOT is Real, i=%d --> %f (B1=%f)\n" i r.[1, i] b1
581 #endif
582 Array.sortInPlace ychk
583
584 #if DEBUG_LOG
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]
588 #endif
589
590 let mutable nintpts = 0
591
592 let xint = Array.zeroCreate 4
593 let yint = Array.zeroCreate 4
594
595 let mutable returnValue = 0.0
596
597 let mutable i = 0
598 while returnValue = 0.0 && i < nychk do
599 #if DEBUG_LOG
600 printf "------------->i=%d (nychk=%d)\n" i nychk
601 #endif
602
603 if not (i < nychk - 1 && abs (ychk.[i] - ychk.[i+1]) < EPS / 2.0) then
604 #if DEBUG_LOG
605 printf "check intersecting points. nintps is %d" nintpts
606 #endif
607
608 let x1 = if abs ychk.[i] > b1 then 0.0 else a1 * sqrt (1.0 - (ychk.[i] * ychk.[i]) / (b1 * b1))
609 let x2 = -x1
610
611 #if DEBUG_LOG
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)
614 #endif
615
616 if abs (ellipse2tr x1 ychk.[i] aa bb cc dd ee ff) < EPS then
617 nintpts <- nintpts + 1
618 #if DEBUG_LOG
619 printf "first if x1. acc nintps=%d\n" nintpts
620 #endif
621 if nintpts > 4 then
622 returnValue <- -1.0
623 else
624 xint.[nintpts-1] <- x1
625 yint.[nintpts-1] <- ychk.[i]
626 #if DEBUG_LOG
627 printf "nintpts=%d, xint=%f, x2=%f, i=%d, yint=%f\n" nintpts x1 x2 i ychk.[i]
628 #endif
629
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
632 #if DEBUG_LOG
633 printf "first if x2. nintps=%d, Dx=%f (eps2=%f) \n" nintpts (abs (x2 - x1)) EPS
634 #endif
635 if nintpts > 4 then
636 returnValue <- -1.0
637 else
638 xint.[nintpts-1] <- x2
639 yint.[nintpts-1] <- ychk.[i]
640
641 #if DEBUG_LOG
642 printf "nintpts=%d, x1=%f, xint=%f, i=%d, yint=%f\n" nintpts x1 x2 i ychk.[i]
643 #endif
644
645 #if DEBUG_LOG
646 else
647 printf "i=%d, multiple roots: %f <--------> %f. continue\n" i ychk.[i] ychk.[i-1]
648 #endif
649 i <- i + 1
650
651 if returnValue = -1.0 then
652 None
653 else
654 let area =
655 match nintpts with
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
658 | TANGENT_POINT ->
659 #if DEBUG_LOG
660 printf "one point is tangent\n"
661 #endif
662 nointpts a1 b1 a2 b2 h1 k1 h2 k2 phi_1 phi_2 h2_tr k2_tr aa bb cc dd ee ff
663
664 | INTERSECTION_POINT ->
665 #if DEBUG_LOG
666 printf "check twointpts\n"
667 #endif
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
671 | _ -> -1.0
672
673 if area = -1.0 then
674 None
675 elif nintpts = 0 then
676 Some (float32 area, [||], [||])
677 else
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)